[发明专利]基于星间速度插值原理的卫星重力反演方法

权利要求书

1.一种基于星间速度插值原理的卫星重力反演方法，其特征如下：

步骤一：获取GRACE型卫星观测数据，具体包括：

1.1)通过星载K波段测量仪采集得到星间速度数据：基于t检验准则即罗曼诺夫斯基准则，剔除星间速度数据中存在的粗大误差；基于9阶Lagrange多项式，插值获得间断的星间速度数据；

1.2)通过星载双频GPS接收机采集得到卫星轨道位置r和卫星轨道速度数据：去除卫星轨道存在的重叠期，进行卫星轨道数据的拼接；截掉由于定轨弱约束造成的卫星轨道数据的开始和结束时段处精度较低的数据；基于3σ准则即莱以特准则，剔除卫星轨道数据中存在的粗大误差；

1.3)通过星载加速度计采集得到卫星非保守力f数据：基于t检验准则即罗曼诺夫斯基准则，剔除卫星非保守力数据中存在的粗大误差；基于9阶Lagrange多项式，插值获得间断的卫星非保守力数据；

步骤二：将GRACE型卫星K波段测量仪的高精度星间速度引入双星相对轨道速度矢量的星星连线分量，建立星间速度插值观测方程；

由GRACE型卫星双频GPS接收机的卫星轨道位置r计算由GRACE-A卫星指向GRACE-B卫星的单位矢量e₁₂＝r₁₂/|r₁₂|，其中r₁₂＝r₂-r₁表示双星相对轨道位置矢量，r₁和r₂分别表示双星绝对轨道位置矢量；

由GRACE型卫星K波段测量仪的高精度星间速度和卫星轨道速度联合计算双星相对轨道速度矢量

r·ρ12(tα)=ρ·12(tα)e12(tα)+[r·12(tα)-(r·12(tα)·e12(tα))e12(tα)],---(1)]]>

其中，表示双星相对轨道速度矢量，和分别表示双星绝对轨道速度矢量，t_α表示插值点的时间；

在地心惯性坐标系中，基于Newton插值模型，利用泰勒展开获得双星轨道加速度差的星星连线方向投影：

e12(t)·r··12(t)=Σi=1nβi′Σα=0i(-1)i+αiαe12(t)·r·ρ12(tα),---(2)]]>

其中，表示二项式系数，t表示卫星观测点的时间，t₀表示卫星观测点的初始时刻，Δt表示采样间隔，n表示插值点的个数；计算

r··12=F120+F12T+F12C+f12,---(3)]]>表示作用于双星的相对地球扰动引力，▽表示梯度算子；表示除地球引力外的作用于双星的太阳引力，月球引力，地球固体、海洋、大气和极潮摄动力，广义相对论效应摄动力相对保守力；f₁₂＝f₂-f₁表示作用于双星的大气阻力、太阳光压、地球辐射压力、卫星轨道高度及姿态控制力、经验摄动力相对非保守力；表示作用于双星的相对地球中心引力

F120=-GM(r2|r2|3-r1|r1|3),---(4)]]>

其中，GM表示地球质量M和万有引力常数G之积，表示双星各自的地心半径，x₁₍₂₎，y₁₍₂₎，z₁₍₂₎分别表示双星各自位置矢量r₁₍₂₎的三个分量；

将(3)式和(4)式结果代入(2)式，计算插值方程

e12(t)·▿T12(t)=e12(t)·{Σi=1nβi′Σα=0i(-1)i+αiα[ρ·12(tα)e12(tα)]]>

+(r·12(tα)-(r·12(tα)·e12(tα))e12(tα))],---(5)]]>

+GM(r2(t)|r2(t)|3-r1(t)|r1(t)|3)-F12C(t)-f12(t)}]]>

步骤三：选择6点星间速度插值公式，通过前后各3个卫星观测值内插获得中间位置点的卫星观测值；基于(5)式得到

e12(ti)·▿T12(ti)=e12(ti)·{-160Δt[r·ρ12(ti-3)-9r·ρ12(ti-2)+45r·ρ12(ti-1)]]>

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

+GM(r2(ti)|r2(ti)|3-r1(ti)|r1(ti)|3)-F12C(ti)-f12(ti)}]]>

其中，t_i表示第i个观测点的时间，由(1)式计算获得；▽T₁₂＝▽(T₂-T₁)表示相对扰动位梯度，T(r，θ，λ)表示地球扰动位

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

其中，r，θ和λ分别表示卫星的地心半径、余纬度和经度，R_e表示地球的平均半径，L表示球函数展开的最大阶数；表示规格化的Legendre函数，l表示阶数，m表示次数；和表示待求的规格化地球引力位系数；

步骤四：利用最小二乘原理解算6点星间速度插值矩阵，获得地球引力位系数和进而反演120阶地球重力场；

在地心惯性系中，(6)式的矩阵形式如下

y_k×1＝A_k×p·x_p×1，                 (8)

其中，y_k×1表示(6)式右边的所有项，由GRACE型卫星K波段测量仪的星间速度、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据计算获得，k表示观测点的数量；x_p×1表示最终求解的地球引力位系数向量其中引力位系数按照阶数l排列形成(次数m固定)，p＝L²+2L-3表示待求引力位系数的个数；A_k×p表示k行p列的设计矩阵，在(6)式左边除地球引力位系数外的所有部分，由GRACE型卫星GPS接收机的卫星轨道位置和卫星轨道速度观测数据计算获得；

选择最大阶数L＝120，在(8)式两边同乘最终可求解获得地球引力位系数

xp×1=(Ak×pTAk×p)-1·Ak×pTyk×1.---(9)]]>