[發明專利]利用星間距離插值建立全球重力場模型的方法

無效

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

1.一種利用星間距離插值建立全球重力場模型的方法，包含下列步驟：

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

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

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

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

步驟二：構建星間距離插值觀測方程

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

r(t)=r(t0)+Σj=1nαjΣξ=0j(-1)j+ξjξr(tξ),---(1)]]>

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

在(1)式兩邊同時對時間t求二階導數，可得單星軌道加速度的展開公式

r··(t)=Σj=1nαj′′Σξ=0j(-1)j+ξjξr(tξ).---(2)]]>

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

r··12(t)=Σj=1nαj′′Σξ=0j(-1)j+ξjξr12(tξ),---(3)]]>

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

將(3)式中的投影到星星連線方向(LOS)可得

e12(t)·r··12(t)=Σj=1nαj′′Σξ=0j(-1)j+ξjξe12(t)·r12(tξ).---(4)]]>

其中，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A衛星指向GRACE-B衛星的單位矢量；e₁₂(t)·r₁₂(t_ξ)可改寫為

e12(t)·r12(tξ)=e12(t)·[r12||(tξ)+r12⊥(tξ)],---(5)]]>

其中，r12||(tξ)=(r12·e12)e12]]>表示r₁₂的星星連線方向分量；r12⊥(tξ)=r12-(r12·e12)e12]]>表示r₁₂的垂直于星星連線方向分量；

通過將GRACE衛星K波段測量儀高精度的星間距離ρ₁₂e₁₂替換(r₁₂·e₁₂)e₁₂，(5)式可改寫為

e₁₂(t)·r₁₂(t_ξ)＝e₁₂(t)·[ρ₁₂(t_ξ)e₁₂(t_ξ)+{r₁₂(t_ξ)-[r₁₂(t_ξ)·e₁₂(t_ξ)]e₁₂(t_ξ)}]，(6)將(6)式代入(4)式可得

e12(t)·r··12(t)=Σj=1nαj′′Σξ=0j(-1)j+ξjξe12(t)·rρ12(tξ),---(7)]]>

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

在(7)式中，的具體形式表示如下

r··12=g120+g12T+a12+f12,---(8)]]>

其中，表示作用于雙星的相對地球擾動引力，▽表示梯度算子；a₁₂＝a₂-a₁表示除地球引力之外的作用于雙星的相對保守力；f₁₂＝f₂-f₁表示作用于雙星的相對非保守力；表示作用于雙星的相對地球中心引力

g120=-GM(r2|r2|3-r1|r1|3),---(9)]]>

其中，GM表示地球質量M和萬有引力常數G之積，表示雙星各自的地心半徑，x₁₍₂₎，y₁₍₂₎，z₁₍₂₎分別表示雙星各自位置矢量r₁₍₂₎的三個分量；

將(8)式和(9)式代入(7)式，星間距離插值觀測方程表示如下

e12(t)·▿T12(t)=e12(t)·{Σj=1nαj′′Σξ=0j(-1)j+ξjξrρ12(tξ)]]>

+GM[r2(t)|r2(t)|3-r1(t)|r1(t)|3]-a12(t)-f12(t)},---(10)]]>

其中，▽T₁₂＝▽(T₂-T₁)表示相對擾動位梯度，T(r，θ，λ)表示地球擾動位

T(r,θ,λ)=GMReΣl=2L(Rer)l+1Σm=0l(C‾lmcosmλ+S‾lmsinmλ)P‾lm(cosθ),---(11)]]>

其中，r，θ和λ分別表示衛星的地心半徑、余緯度和經度，R_e表示地球的平均半徑，L表示球函數展開的最大階數；表示規格化的Legendre函數，l表示階數，m表示次數；和表示待求的規格化引力位系數；

將(10)式按照泰勒公式展開，3點、5點、7點和9點星間距離插值公式表示如下

e12(tj)·▿T12(tj)=e12(tj)·{1(Δt)2[rρ12(tj-1)-2rρ12(tj)+rρ12(tj+1)]]]>

+GM[r2(tj)|r2(tj)|3-r1(tj)|r1(tj)|3]-a12(tj)-f12(tj)},---(12)]]>

e12(tj)·▿T12(tj)=e12(tj)·{1(Δt)2[-112rρ12(tj-2)+43rρ12(tj-1)]]>

-52rρ12(tj)+43rρ12(tj+1)-112rρ12(tj+2)]]]>

+GM[r2(tj)|r2(tj)|3-r1(tj)|r1(tj)|3]-a12(tj)-f12(tj)},---(13)]]>

e12(tj)·▿T12(tj)=e12(tj)·{1(Δt)2[190rρ12(tj-3)-320rρ12(tj-2)+32rρ12(tj-1)]]>

-4918rρ12(tj)+32rρ12(tj+1)-320rρ12(tj+2)+190rρ12(tj+3)]]]>

+GM[r2(tj)|r2(tj)|3-r1(tj)|r1(tj)|3]-a12(tj)-f12(tj)},---(14)]]>

e12(tj)·▿T12(tj)=e12(tj)·{1(Δt)2[-1560rρ12(tj-4)+8315rρ12(tj-3)-15rρ12(tj-2)+85rρ12(tj-1)]]>

-20572rρ12(tj)+85rρ12(tj+1)-15rρ12(tj+2)+8315rρ12(tj+3)-1560rρ12(tj+4)]]]>

+GM[r2(tj)|r2(tj)|3-r1(tj)|r1(tj)|3]-a12(tj)-f12(tj)},---(15)]]>

步驟三：優選不同點數星間距離插值公式

利用“步驟一”中獲取的GRACE星載K波段測量儀的星間距離、星載GPS接收機的衛星軌道位置和衛星軌道速度、以及星載加速度計的衛星非保守力觀測數據，基于“步驟二”中3點、5點、7點和9點星間距離插值公式(12)～(15)分別反演地球重力場，可得到9點星間距離插值公式(15)有利于120階GRACE衛星重力反演精度的有效提高；

步驟四：建立全球重力場模型

利用“步驟一”中獲取的GRACE星載K波段測量儀的星間距離、星載GPS接收機的衛星軌道位置和衛星軌道速度、以及星載加速度計的衛星非保守力觀測數據，基于“步驟二”中建立的9點星間距離插值公式(15)，計算獲得地球引力位系數和最終通過引力位系數的集合建立120階全球重力場模型WHIGG-GEGM01S；在120階處，大地水準面累積誤差為10.98cm，重力異常累積誤差為1.741×10^-6m/s²；

步驟五：驗證全球重力場模型的正確性和可靠性

全球重力場模型的大地水準面標準誤差為0.726m，更接近于國際已公布的全球重力場模型EIGEN-GRACE02S的大地水準面標準誤差0.735m。