CN102313905A - 基于星间速度插值原理的卫星重力反演方法 - Google Patents

Abstract

本发明涉及一种地球重力场精密测量的方法，特别是一种基于星间速度插值原理的卫星重力反演方法；通过将高精度的星间速度引入双星相对轨道速度矢量的星星连线分量，建立星间速度插值观测方程，进而精确和快速反演地球重力场；该方法卫星重力反演精度高，卫星观测方程物理含义明确，易于重力卫星系统敏感度分析，利于探测中短波重力场信号，对计算机性能要求低；星间速度插值法是解算高精度和高空间分辨率地球重力场的有效方法。

Description

基于星间速度插值原理的卫星重力反演方法

一、技术领域

本发明涉及卫星大地测量学、地球物理学、空间科学等交叉技术领域，特别是涉及一种通过将GRACE型卫星K波段测量仪的高精度星间速度引入双星相对轨道速度矢量的星星连线分量建立新型星间速度插值观测方程，进而精确和快速反演地球重力场的方法。

二、背景技术

如图1所示，地球重力场及其随时间的变化反映地球表层及内部物质的空间分布、运动和变化，同时决定着大地水准面的起伏和变化。因此，确定地球重力场的精细结构及其时变不仅是大地测量、地球物理、海洋勘探、地震预报、空间科学、惯性导航、航天技术、工业应用等的需求，同时也将为全人类寻求资源、保护环境和预测灾害提供重要的信息资源。第一，重力变化决定于地质构造和矿产资源(如煤炭、石油、天然气等)分布的不均匀性。由于我国远景资源储备中有78％的石油和93％的天然气亟待开发，因此精密重力测量有利于我国将来的资源勘探。第二，由于地震(如2005年3月印尼苏门答腊8.5级地震、2008年5月中国汶川8.0级地震、2010年2月智利8.8级地震、2011年3月日本福岛9.0级地震等)、火山、海啸、龙卷风、泥石流等自然灾害发生前后，地球重力场将发生明显的变化和异常，因此精密重力测量能帮助我们预测自然灾害的发生，进而降低经济损失和保障人员安全。

早在20世纪60年代，Baker首次提出了利用卫星跟踪卫星测量技术(SST)确定地球重力场的重要思想(如图2所示)。自此以后，国际大地测量学界的许多学者都积极投身于地球重力场反演的方法与算法的理论研究和数值计算之中。在利用卫星重力测量数据反演地球重力场的众多方法中，按引力位系数解算方式的差异可分为空域法和时域法。时域法主要包括5种类型：(1)轨道摄动法；(2)动力学法；(3)能量守恒法；(4)半解析法；(5)卫星加速度法。但直到2000年以后，专用于地球重力场精密探测的重力卫星CHAMP(ChallengingMinisatellite Payload)、GRACE(Gravity Recovery and Climate Experiment)和GOCE(Gravity Field and Steady-State Ocean Circulaion Explorer)的成功发射及实际观测数据的有效处理应用，上述方法与算法才真正进入实际操作阶段(如图3所示)。

卫星加速度法包括轨道加速度法和星间加速度法，是指以卫星轨道加速度或星间加速度为观测量建立卫星观测方程，基于最小二乘法、预处理共轭梯度法(PCCG)等解算大型线性超定方程组，进而恢复地球引力位系数，最终目的是反演高精度和高空间分辨率的地球重力场。卫星轨道加速度法通常用于解算CHAMP单星地球重力场。Reubelt et al.(2003)基于牛顿插值的轨道加速度法确定了CHAMP地球重力场；Ditmar et al.(2006)提出了基于3点插值加权平均加速度法计算地球重力场的思想，并利用1年的CHAMP卫星数据建立了DEOSCHAMP-01C 70地球重力场模型。星间加速度法适于确定GRACE双星地球重力场，优点是观测方程形式简单、物理含义明确、易于重力场的敏感度分析、在保证求解精度的前提下计算量大大降低，通常采用PC计算机可完成高阶重力场的快速求解；缺点是采用的数值微分算法在一定程度上损失了低频地球重力场的精度。Rummel(1979)基于经典的星间加速度法解算了地球重力场；Liu(2008)基于3点星间距离联合法(3RC)解算了GRACE地球重力场。在以往研究建立的卫星观测方程中，由于引入了星载K波段测量仪的高精度星间距离、星间速度和星间加速度观测量，因此地球重力场反演的精度得到了较大程度的改善。由于星载K波段测量仪星间速度是地球重力场精度提高的关键因素，同时也是主要误差源，因此开展星载K波段测量仪星间速度对地球重力场精度影响因素的分析研究是进一步提高地球重力场精度的关键。

不同于前人提出的卫星重力反演法，本发明首次通过在相对轨道速度的星星连线分量中引入星载K波段测量仪的高精度星间速度，建立了新型星间速度插值观测方程，并通过与德国地学研究中心(GFZ)公布的EIGEN-GRACE02S地球重力场模型精度的符合性验证了其正确性和可靠性。由于我国自主研制和正在建设的首期卫星跟踪卫星测量模式的地球重力卫星系统预计于国家“十二五”规划末期发射升空，因此，星间速度插值法(IRRI)以其独特的优越性将成为我国将来高精度和高空间分辨率地球重力场反演的优选方法之一。

三、发明内容

本发明的目的是：通过将星载K波段测量仪的高精度星间速度引入双星相对轨道速度矢量的星星连线分量，建立新型星间速度插值观测方程，进而精确和快速反演地球重力场。

为达到上述目的，本发明采用了如下技术方案：

基于星间速度插值原理的精确和快速卫星重力反演方法包含下列步骤：

(1)利用GRACE型卫星观测数据并进行预处理，具体包括：

1.1)采集星载K波段测量仪得到的星间速度

数据：基于t检验准则即罗曼诺夫斯基准则，剔除星间速度数据中存在的粗大误差；基于9阶Lagrange多项式，插值获得间断的星间速度数据；

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

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

(2)通过将GRACE型卫星K波段测量仪的高精度星间速度引入双星相对轨道速度矢量的星星连线分量，建立新型星间速度插值观测方程。

(3)利用预处理后的GRACE型卫星观测数据，对比论证2点、4点、6点和8点星间速度插值公式对地球重力场反演精度的影响。

(4)利用最小二乘原理解算6点星间速度插值公式矩阵，获得地球引力位系数，进而反演120阶地球重力场。

本发明是基于新型星间速度插值法有利于快速反演高精度和高空间分辨率地球重力场的特点而设计的，优点是：1)卫星重力反演精度高；2)卫星观测方程物理含义明确；3)易于重力卫星系统敏感度分析；4)利于感测中短波重力场信号；5)对计算机性能要求低。

四、附图说明

图1表示由于地球内部、表面和空间质量(陆地、海洋、大气等)的迁移引起地球重力场的变化。

图2表示卫星跟踪卫星高低/低低观测模式的测量原理；不仅以高轨GPS卫星对低轨重力双星实时精密跟踪定轨(SST-HL)，同时以差分原理测定低轨双星之间的相互运动(SST-LL)。

图3表示国际成功发射的三颗重力卫星CHAMP(德国航天局，2000-07-15)、GRACE(美国宇航局和德国航天局，2002-03-17)和GOCE(欧洲空间局，2009-03-17)。

图4表示论证2点、4点、6点和8点星间速度插值公式对地球重力场精度的影响，其中横坐标表示地球引力位按球函数展开的阶数，纵坐标表示大地水准面累积误差(单位：m)。

图5表示基于新型星间速度插值法反演地球重力场的精度，其中横坐标表示地球引力位按球函数展开的阶数，纵坐标表示大地水准面累积误差(单位：m)。

五、具体实施方式

以下结合附图，对本发明的具体实施方式作进一步的说明。

基于星间速度插值原理的卫星重力反演方法包含下列步骤：

步骤一：获取GRACE型卫星观测数据

1.1采集星载K波段测量仪得到的星间速度

数据

(1)基于t检验准则(罗曼诺夫斯基准则)，有效剔除星间速度数据中存在的粗大误差；

(2)基于9阶Lagrange多项式，插值获得间断的星间速度数据。1.2采集星载双频GPS接收机得到的卫星轨道位置r和轨道速度

数据

(1)为了保证卫星轨道数据的精度和连续性，有效去除卫星轨道存在的重叠期，进而完成卫星轨道数据的拼接；

(2)有效截掉由于定轨弱约束造成的卫星轨道数据的开始和结束时段处精度较低的数据；

(3)基于3σ准则(莱以特准则)，有效剔除卫星轨道数据中存在的粗大误差。

1.3采集星载加速度计得到的卫星非保守力f数据

(1)基于t检验准则(罗曼诺夫斯基准则)，有效剔除卫星非保守力数据中存在的粗大误差；

(2)基于9阶Lagrange多项式，插值获得间断的卫星非保守力数据。

步骤二：由于目前GPS卫星轨道速度测量精度相对较低，因此通过将GRACE型卫星K波段测量仪的高精度星间速度引入双星相对轨道速度矢量的星星连线分量，建立新型星间速度插值观测方程。

在地心惯性坐标系中，基于Newton插值模型，单星轨道速度

的泰勒展开表示如下

$\stackrel{\·}{r}\left(t\right)=\stackrel{\·}{r}\left(t_0\right)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right)\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right),---\left(1\right)$

其中，

表示二项式系数，

t表示卫星观测点的时间，t₀表示卫星观测点的初始时刻，Δt表示采样间隔，n表示插值点的个数，t_α表示插值点的时间。

在(1)式两边同时对时间t求一阶导数，可得单星轨道加速度

的展开公式

$\stackrel{\·\·}{r}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)}^{\′}\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right)\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right).---\left(2\right)$

基于(2)式，双星轨道加速度差

的展开公式表示如下

${\stackrel{\·\·}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)}^{\′}\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right){\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right),---\left(3\right)$

其中，

表示双星相对轨道加速度矢量，

和

分别表示双星绝对轨道加速度矢量；

表示双星相对轨道速度矢量，

和

分别表示双星绝对轨道速度矢量。

将(3)式中的

投影到星星连线方向可得

$e_{12}\left(t\right)\·{\stackrel{\·\·}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)}^{\′}\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right)e_{12}\left(t\right)\·{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right),---\left(4\right)$

其中，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A卫星指向GRACE-B卫星的单位矢量，r₁₂＝r₂-r₁表示双星相对轨道位置矢量，r₁和r₂分别表示双星绝对轨道位置矢量。

如(4)式所示，假如直接利用

反演地球重力场，其结果只相当于同一轨道上两个单颗CHAMP卫星相互跟踪，对地球重力场反演精度的提高没有本质贡献。因此，引入GRACE型卫星K波段测量仪的高精度星间速度

是大幅度提高地球重力场反演精度的关键。基于以上原因，可将(4)式等号右边的

改写为

$e_{12}\·{\stackrel{\·}{r}}_{12}=e_{12}\·({\stackrel{\·}{r}}_{12}^{\left|\right|}+{\stackrel{\·}{r}}_{12}^{\⊥}),---\left(5\right)$

其中， ${\stackrel{\·}{r}}_{12}^{\left|\right|}=({\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 表示

的星星连线方向分量， ${\stackrel{\·}{r}}_{12}^{\⊥}={\stackrel{\·}{r}}_{12}-({\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 表示

的垂直于星星连线方向分量。

据误差原理可知，(5)式中的

远大于为了降低

本发明引入了GRACE型卫星K波段测量仪的高精度星间速度

通过将

替换(5)式中的

再代入(4)式可得

$e_{12}\left(t\right)\·{\stackrel{\·\·}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)}^{\′}\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right)e_{12}\left(t\right)\·{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{\mathrm{\α}}\right),---\left(6\right)$

其中， ${\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{\mathrm{\α}}\right)={\stackrel{\·}{\mathrm{\ρ}}}_{12}\left(t_{\mathrm{\α}}\right)e_{12}\left(t_{\mathrm{\α}}\right)+[{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)-({\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)\·e_{12}\left(t_{\mathrm{\α}}\right))e_{12}\left(t_{\mathrm{\α}}\right)].$

在(6)式中，

的具体形式表示如下

${\stackrel{\·\·}{r}}_{12}=F_{12}^0+F_{12}^T+F_{12}^C+f_{12},---\left(7\right)$

其中，

表示作用于双星的相对地球扰动引力，▽表示梯度算子；

表示除地球引力外的作用于双星的所有相对保守力(太阳引力，月球引力，地球固体、海洋、大气和极潮摄动力，广义相对论效应摄动力等)；f₁₂＝f₂-f₁表示作用于双星的所有相对非保守力(大气阻力、太阳光压、地球辐射压力、卫星轨道高度及姿态控制力、经验摄动力等)；

表示作用于双星的相对地球中心引力

$F_{12}^0=-\mathrm{GM}(\frac{r_2}{{\left|r_2\right|}^3}-\frac{r_1}{{\left|r_1\right|}^3}),---\left(8\right)$

其中，GM表示地球质量M和万有引力常数G之积，

表示双星各自的地心半径，x₁₍₂₎，y₁₍₂₎，z₁₍₂₎分别表示双星各自位置矢量r₁₍₂₎的三个分量。

将(7)式和(8)式代入(6)式，星间速度插值观测方程表示如下

$e_{12}\left(t\right)\·\▿T_{12}\left(t\right)=e_{12}\left(t\right)\·\left\{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)}^{\′}\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right)\right[{\stackrel{\·}{\mathrm{\ρ}}}_{12}\left(t_{\mathrm{\α}}\right)e_{12}\left(t_{\mathrm{\α}}\right)$

$+({\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)-({\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)\·e_{12}\left(t_{\mathrm{\α}}\right))e_{12}\left(t_{\mathrm{\α}}\right))],---\left(9\right)$

$+\mathrm{GM}(\frac{r_2\left(t\right)}{{\left|r_2\left(t\right)\right|}^3}-\frac{r_1\left(t\right)}{{\left|r_1\left(t\right)\right|}^3})-F_{12}^C\left(t\right)-f_{12}\left(t\right)\}$

其中，t表示卫星观测点时间，t_α表示插值点时间；▽T₁₂＝▽(T₂-T₁)表示相对扰动位梯度，T(r，θ，λ)表示地球扰动位

$T(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right),---\left(10\right)$

其中，r，θ和λ分别表示卫星的地心半径、余纬度和经度，R_e表示地球的平均半径，L表示球函数展开的最大阶数；

表示规格化的Legendre函数，l表示阶数，m表示次数；

和

表示待求的规格化地球引力位系数。

将(9)式按泰勒公式展开，建立2点、4点、6点和8点星间速度插值公式，计算原理为选取观测值时间序列中的2i+1个连续观测点，将第i+1个观测点前后各i个卫星观测值代入星间速度插值观测方程中等号的右边项，最终插值获得中间位置第i+1个点的观测值，公式表示如下

$e_{12}\left(t_i\right)\·\▿T_{12}\left(t_i\right)=e_{12}\left(t_i\right)\·\{-\frac{1}{2\mathrm{\Δt}}\·[{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-1}\right)-{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+1}\right)]$

$+\mathrm{GM}(\frac{r_2\left(t_i\right)}{{\left|r_2\left(t_i\right)\right|}^3}-\frac{r_1\left(t_i\right)}{{\left|r_1\left(t_i\right)\right|}^3})-F_{12}^C\left(t_i\right)-f_{12}\left(t_i\right)\},---\left(11\right)$

$e_{12}\left(t_i\right)\·\▿T_{12}\left(t_i\right)=e_{12}\left(t_i\right)\·\{\frac{1}{12\mathrm{\Δt}}\·[{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-2}\right)-8{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-1}\right)+8{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+1}\right)-{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+2}\right)]$

$+\mathrm{GM}(\frac{r_2\left(t_i\right)}{{\left|r_2\left(t_i\right)\right|}^3}-\frac{r_1\left(t_i\right)}{{\left|r_1\left(t_i\right)\right|}^3})-F_{12}^C\left(t_i\right)-f_{12}\left(t_i\right)\},---\left(12\right)$

$e_{12}\left(t_i\right)\·\▿T_{12}\left(t_i\right)=e_{12}\left(t_i\right)\·\{-\frac{1}{60\mathrm{\Δt}}\·[{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-3}\right)-9{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-2}\right)+45{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-1}\right)$

$-45{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+1}\right)+9{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+2}\right)-{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+3}\right)],---\left(13\right)$

$+\mathrm{GM}(\frac{r_2\left(t_i\right)}{{\left|r_2\left(t_i\right)\right|}^3}-\frac{r_1\left(t_i\right)}{{\left|r_1\left(t_i\right)\right|}^3})-F_{12}^C\left(t_i\right)-f_{12}\left(t_i\right)\}$

$e_{12}\left(t_i\right)\·\▿T_{12}\left(t_i\right)=e_{12}\left(t_i\right)\·\{\frac{1}{\mathrm{\Δt}}\·[\frac{1}{280}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-4}\right)-\frac{4}{105}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-3}\right)+\frac{1}{5}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-2}\right)-\frac{4}{5}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-1}\right)$

$+\frac{4}{5}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+1}\right)-\frac{1}{5}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+2}\right)+\frac{4}{105}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+3}\right)-\frac{1}{280}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+4}\right)],---\left(14\right)$

$+\mathrm{GM}(\frac{r_2\left(t_i\right)}{{\left|r_2\left(t_i\right)\right|}^3}-\frac{r_1\left(t_i\right)}{{\left|r_1\left(t_i\right)\right|}^3})-F_{12}^C\left(t_i\right)-f_{12}\left(t_i\right)\}$

其中，t_i表示卫星的观测点时间，i＝1，2，3......。

步骤三：将“步骤一”中获取的GRACE型卫星K波段测量仪的星间速度、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，分别代入“步骤二”中建立的2点、4点、6点和8点星间速度插值公式(11)～(14)，对比论证了不同点数星间速度插值公式对地球重力场精度的影响。

如图4所示，基于较优的信噪比，6点星间速度插值公式有利于卫星重力反演精度的提高；6点星间速度插值原理为通过前后各3个卫星观测值内插获得中间位置点的卫星观测值。第一，在120阶内，基于2点星间速度插值公式反演地球重力场的精度远低于分别基于4点、6点和8点星间速度插值公式的反演精度。原因分析如下：首先，由于(11)式的左边是点域值，而右边是平均值，因此(11)式的左右两边几乎不相等；其次，由于2点星间速度插值公式的插值点数较少，因此无法提供足够的地球重力场反演插值信息。第二，在120阶内，基于6点星间速度插值公式(13)反演地球重力场的精度高于分别基于2点、4点和8点星间速度插值公式的反演精度。原因分析如下：由于随着插值点数的增多，卫星观测值的信息量逐渐增加，因此基于6点星间速度插值公式反演地球重力场的精度高于分别基于2点和4点星间速度插值公式的反演精度。但是，随着插值点数的增多，卫星观测值的误差量也同时增加，因此基于8点星间速度插值公式反演地球重力场的精度低于基于6点星间速度插值公式的反演精度。综上所述，6点星间速度插值公式是提高120阶地球重力场反演精度的较优选择。

步骤四：将“步骤一”中获取的GRACE型卫星K波段测量仪的星间速度、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，代入“步骤二”中建立的6点星间速度插值公式(13)，通过解算最小二乘超定方程组(方程的个数大于未知数的个数)获得地球引力位系数

和

进而反演120阶地球重力场。

在地心惯性系中，6点星间速度插值公式(13)的矩阵形式如下

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

其中，y_k×1表示(13)式右边的所有项，由GRACE型卫星K波段测量仪的星间速度、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据计算获得，k表示观测点的数量；x_p×1表示最终求解的地球引力位系数向量

其中引力位系数按照阶数l排列形成(次数m固定)，p＝L²+2L-3表示待求引力位系数的个数；A_k×p表示k行p列的设计矩阵，在(13)式左边除地球引力位系数

外的所有部分，由GRACE型卫星GPS接收机的卫星轨道位置和卫星轨道速度观测数据计算获得。

在(15)式两边同乘

可得最终求解的地球引力位系数

$x_{p\×1}={\left(A_{k\×p}^T{A}_{k\×p}\right)}^{-1}\·A_{k\×p}^T{y}_{k\×1}.---\left(16\right)$

图5中星号线表示德国GFZ公布的120阶EIGEN-GRACE02S全球重力场模型的精度，在120阶处大地水准面累积误差为18.938cm；虚线表示基于新型和精确星间速度插值法，利用预处理后的GRACE型卫星K波段测量仪的星间速度、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，反演GRACE地球重力场精度，在120阶处大地水准面累积误差为11.401cm。通过2条曲线的符合性可知，星间速度插值法是解算高精度和高空间分辨率地球重力场的有效方法。

Claims (1)

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波段测量仪的高精度星间速度

和卫星轨道速度

联合计算双星相对轨道速度矢量

${\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{\mathrm{\α}}\right)={\stackrel{\·}{\mathrm{\ρ}}}_{12}\left(t_{\mathrm{\α}}\right)e_{12}\left(t_{\mathrm{\α}}\right)+[{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)-({\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)\·e_{12}\left(t_{\mathrm{\α}}\right))e_{12}\left(t_{\mathrm{\α}}\right)],---\left(1\right)$

其中，

表示双星相对轨道速度矢量，

和

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

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

$e_{12}\left(t\right)\·{\stackrel{\·\·}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)}^{\′}\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right)e_{12}\left(t\right)\·{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{\mathrm{\α}}\right),---\left(2\right)$

其中，

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

${\stackrel{\·\·}{r}}_{12}=F_{12}^0+F_{12}^T+F_{12}^C+f_{12},---\left(3\right)$ 表示作用于双星的相对地球扰动引力，▽表示梯度算子；

表示除地球引力外的作用于双星的太阳引力，月球引力，地球固体、海洋、大气和极潮摄动力，广义相对论效应摄动力相对保守力；f₁₂＝f₂-f₁表示作用于双星的大气阻力、太阳光压、地球辐射压力、卫星轨道高度及姿态控制力、经验摄动力相对非保守力；

表示作用于双星的相对地球中心引力

$F_{12}^0=-\mathrm{GM}(\frac{r_2}{{\left|r_2\right|}^3}-\frac{r_1}{{\left|r_1\right|}^3}),---\left(4\right)$

其中，GM表示地球质量M和万有引力常数G之积，

表示双星各自的地心半径，x₁₍₂₎，y₁₍₂₎，z₁₍₂₎分别表示双星各自位置矢量r₁₍₂₎的三个分量；

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

$e_{12}\left(t\right)\·\▿T_{12}\left(t\right)=e_{12}\left(t\right)\·\left\{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)}^{\′}\underset{\mathrm{\α}=0}{\overset{i}{\mathrm{\Σ}}}{(-1)}^{i+\mathrm{\α}}\left(\begin{array}{c}i\\ \mathrm{\α}\end{array}\right)\right[{\stackrel{\·}{\mathrm{\ρ}}}_{12}\left(t_{\mathrm{\α}}\right)e_{12}\left(t_{\mathrm{\α}}\right)$

$+({\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)-({\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)\·e_{12}\left(t_{\mathrm{\α}}\right))e_{12}\left(t_{\mathrm{\α}}\right))],---\left(5\right)$

$+\mathrm{GM}(\frac{r_2\left(t\right)}{{\left|r_2\left(t\right)\right|}^3}-\frac{r_1\left(t\right)}{{\left|r_1\left(t\right)\right|}^3})-F_{12}^C\left(t\right)-f_{12}\left(t\right)\}$

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

$e_{12}\left(t_i\right)\·\▿T_{12}\left(t_i\right)=e_{12}\left(t_i\right)\·\{-\frac{1}{60\mathrm{\Δt}}[{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-3}\right)-9{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-2}\right)+45{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i-1}\right)$

$-45{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+1}\right)+9{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+2}\right)-{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{i+3}\right)],---\left(6\right)$

$+\mathrm{GM}(\frac{r_2\left(t_i\right)}{{\left|r_2\left(t_i\right)\right|}^3}-\frac{r_1\left(t_i\right)}{{\left|r_1\left(t_i\right)\right|}^3})-F_{12}^C\left(t_i\right)-f_{12}\left(t_i\right)\}$

其中，t_i表示第i个观测点的时间，

由(1)式计算获得；▽T₁₂＝▽(T₂-T₁)表示相对扰动位梯度，T(r，θ，λ)表示地球扰动位

$T(r,\mathrm{\θ},\mathrm{\λ})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\θ}\right),---\left(7\right)$

其中，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)式两边同乘

最终可求解获得地球引力位系数

$x_{p\×1}={\left(A_{k\×p}^T{A}_{k\×p}\right)}^{-1}\·A_{k\×p}^T{y}_{k\×1}.---\left(9\right)$

Images