[發明專利]基于雙星能量插值原理的衛星重力反演方法

無效

• 申請號: 201110204206.2
• 申請日: 2011-07-20
• 公開（公告）號: CN102393535A
• 公開（公告）日: 2012-03-28
• 發明（設計）人: 鄭偉;許厚澤;熊熊;鐘敏;劉成恕
• 申請（專利權）人: 中國科學院測量與地球物理研究所
• 主分類號: G01V7/00
• 分類號: G01V7/00
• 代理公司: 暫無信息
• 代理人: 暫無信息

1.一種基于雙星能量插值原理的衛星重力反演方法，包含下列步驟：

步驟一：對衛星觀測數據進行預處理，具體包括：

1.1)采集星載K波段測量儀得到的星間距離數據ρ₁₂：基于t檢驗準則即羅曼諾夫斯基準則，剔除星間距離數據中存在的粗大誤差；基于9階Lagrange多項式，插值獲得間斷的星間距離數據；

1.2)采集星載雙頻GPS接收機得到的衛星軌道數據，包括軌道位置r和軌道速度去除衛星軌道存在的重疊期，進行軌道數據的拼接；截掉由于定軌弱約束造成的衛星軌道數據的開始和結束時段處精度較低的數據；基于3σ準則即萊以特準則，剔除軌道數據中存在的粗大誤差；

1.3)采集星載加速度計得到的衛星非保守力數據f：基于t檢驗準則即羅曼諾夫斯基準則，剔除衛星非保守力數據中存在的粗大誤差；基于9階Lagrange多項式，插值獲得間斷的衛星非保守力數據；

步驟二：基于雙星能量插值原理反演全球重力場

2.1)在地心慣性系中，建立衛星運動方程如下，

r··(r,t′)=Fe(r,t′)+FT(r,t′)+f(r,t′),---(1)]]>

其中，表示衛星總加速度，r表示衛星的絕對位置矢量，t′表示觀測時間；F_e(r，t′)表示地球的引力；F_T(r，t′)表示三體攝動力；f(r，t′)表示作用于衛星的非保守力；在式(1)兩邊同時乘以軌道速度

r··r··(r,t′)=r··[Fe(r,t′)+FT(r,t′)]+r··f(r,t′),---(2)]]>

其中，F_e(r，t′)和F_T(r，t′)可表示為

Fe(T)(r,t′)=∂Ve(T)∂r,---(3)]]>

其中，V_e＝V₀+T_e表示地球引力位，表示地球中心引力位，GM表示地球質量M和萬有引力常數G之積，表示衛星的地心半徑，x，y，z表示位置矢量r的分量，T_e表示擾動位；V_T表示三體攝動能；

V_e(T)對時間的一階導數表示如下：

dVe(T)dt′=∂Ve(T)∂r·drdt′+∂Ve(T)∂t′·dt′dt′,---(4)]]>

將式(3)代入式(4)中可得

Fe(T)(r,t′)·r·=dVe(T)dt′-∂Ve(T)∂t′,---(5)]]>

將式(5)代入式(2)，并兩邊同時積分；在地心慣性系中，單星擾動位觀測方程如下：

T_e＝E_k-E_f+V_ω-V_T-V₀-E₀，????????(6)

其中，表示衛星的動能；表示衛星的耗散能；Vω=∫∂(Ve+VT)∂t′dt′≈-ωe(xy·-yx·)]]>表示衛星的位旋轉能，ω_e表示地球自轉角速度，表示衛星軌道速度的3個分量；E₀表示衛星系統的能量積分常數；

據式(6)，在地心慣性系中，雙星擾動位差的觀測方程建立如下：

T_e12＝E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂，????????(7)

其中，T_e12表示雙星擾動位差Te12(r1,θ1,λ1,r2,θ2,λ2)=GMReΣl=2LΣm=-ll{[(Rer2)l+1Y‾lm(θ2,λ2)-(Rer1)l+1Y‾lm(θ1,λ1)]C‾lm}]]>其中，Y‾l,m(θ,λ)=P‾l,|m|(cosθ)Qm(λ),]]>Qm(λ)=cosmλm≥0sin|m|λm<0;]]>r₁₍₂₎，θ₁₍₂₎，λ₁₍₂₎表示雙星各自的地心半徑、地心余緯度和地心經度，R_e表示地球的平均半徑；表示規格化的Legendre函數，其中l表示階數，m表示次數；表示待求的規格化地球引力位系數；表示雙星動能差，Ef12=∫(r·2·f2-r·1·f1)dt′]]>表示雙星耗散能差，Vω12=-ωe(x12y·2-y2x·12-y12x·1+x1y·12)]]>表示雙星位旋轉能差，V_T12表示雙星三體攝動能差，表示雙星中心引力位差，E₀₁₂表示雙星能量積分常數差，可通過初始位置和速度計算得到；

2.2)將GRACE衛星的K波段測量儀的精確星間距離數據ρ₁₂引入到式(7)中的項E_k12，建立雙星能量插值衛星觀測方程；包括：

在地心慣性坐標系中，基于Newton插值模型，單星軌道位置r的泰勒展開表示式，如下

r(t)=r(t0)+Σi=1nλiΣτ=0i(-1)i+τiτr(tτ),---(8)]]>

其中，表示二項式系數，t表示插值點的時間，t₀表示插值點的初始時刻，Δt表示采樣間隔，n表示插值點的個數；

在式(8)兩邊同時對t求一階導數，可得單星軌道速度的展開公式：

r·(t)=Σi=1nλi′Στ=0i(-1)i+τiτr(tτ).---(9)]]>

基于式(9)，雙星軌道速度差的展開公式表示如下

r·12(t)=Σi=1nλi′Στ=0i(-1)i+τiτr12(tτ),---(10)]]>

其中，r₁₂＝r₂-r₁和分別表示雙星相對軌道位置矢量和相對軌道速度矢量，r₁和r₂分別表示雙星絕對軌道位置矢量，和分別表示雙星絕對軌道速度矢量；

基于式(10)，雙星相對動能的展開公式表示如下：

12[r·2(t)+r·1(t)]·r·12(t)=Σi=1nλi′Στ=0i(-1)i+τiτ12[r·2(t)+r·1(t)]·r12(tτ),---(11)]]>

其中，12[r·2(t)+r·1(t)]·r12(tτ)]]>可被改寫為：

12[r·2(t)+r·1(t)]·r12(tτ)=12[r·2(t)+r·1(t)]·[r12||(tτ)+r12⊥(tτ)],---(12)]]>

其中，表示r₁₂的星星連線方向分量，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A衛星指向GRACE-B衛星的單位矢量；表示r₁₂的垂直于星星連線方向分量；

使用GRACE衛星K波段測量儀的高精度星間距離ρ₁₂e₁₂來替換(r₁₂·e₁₂)e₁₂；式(11)可改寫為：

12[r·2(t)+r·1(t)]·r·ρ12(t)=Σi=1nλi′Στ=0i(-1)i+τiτ12[r·2(t)+r·1(t)]·rρ12(tτ),---(13)]]>

其中，r_ρ12(t_τ)＝ρ₁₂(t_τ)e₁₂(t_τ)+{r₁₂(t_τ)-[r₁₂(t_τ)e₁₂(t_τ)]e₁₂(t_τ)}；

因此，2點、4點、6點和8點星間距離插值公式表示如下

r·ρ12(ti)=-12Δt[rρ12(ti-1)-rρ12(ti+1)],---(14)]]>

r·ρ12(ti)=112Δt[rρ12(ti-2)-8rρ12(ti-1)+8rρ12(ti+1)-rρ12(ti+2)],---(15)]]>

r·ρ12(ti)=-160Δt[rρ12(ti-3)-9rρ12(ti-2)+45rρ12(ti-1),---(16)]]>

-45rρ12(ti+1)+9rρ12(ti+2)-rρ12(ti+3)]]]>

r·ρ12(ti)=1Δt[1280rρ12(ti-4)-4105rρ12(ti-3)+15rρ12(ti-2)-45rρ12(ti-1)---(17)]]>

+45rρ12(ti+1)-15rρ12(ti+2)+4105rρ12(ti+3)-1280rρ12(ti+4)]]]>

聯合式(7)和式(13)，雙星能量插值觀測方程表示如下

T_e12(t)＝E_ρ12(t)-E_f12(t)+V_ω12(t)-V_T12(t)-V₀₁₂(t)-E₀₁₂(t)，????(18)

其中，Eρ12(t)=Σi=1nλi′Στ=0i(-1)i+τiτ12[r·2(t)+r·1(t)]·rρ12(tτ)]]>表示雙星動能差，V_T12(t)＝V_E12(t)+V_S12(t)+V_M12(t)，V_E12(t)表示雙星地球固體潮汐能差，V_S12(t)表示雙星太陽引力位差，V_M12(t)表示雙星月球引力位差，E_f12表示雙星耗散能差，V_ω12表示雙星位旋轉能差，V₀₁₂表示雙星中心引力位差，E₀₁₂表示雙星能量積分常數差；

2.3)基于2點、4點、6點和8點雙星能量插值觀測方程(18)分別反演120階GRACE地球重力場。