CN102305949A - 利用星间距离插值建立全球重力场模型的方法 - Google Patents

Abstract

本发明公开了一种利用星间距离插值建立全球重力场模型的方法，通过将星载K波段测量仪的精确星间距离引入双星相对轨道位置矢量的星星连线分量，构建星间距离插值卫星观测方程，进而建立高精度和高空间分辨率的全球重力场模型。该方法卫星重力反演精度高，易于感测中高频重力场信号，利于重力卫星系统误差分析，观测方程物理含义明确，计算机性能要求低。

Description

利用星间距离插值建立全球重力场模型的方法

一、技术领域

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

二、背景技术

如图1所示，重力卫星在地球重力场作用下绕地球作近圆极轨运动，若精密定轨必须知道精确的地球重力场参数；反之，精确测定卫星轨道摄动，利用摄动跟踪观测数据又可以提高地球重力场参数的精度，两者相辅相成。在大地测量领域，地球重力场对研究地球形状和精确求定地面控制点的三维坐标起着重要作用；在固体地球物理学中，基于地球重力场可以研究地球的内部构造和板块运动；在海洋学中，为了研究海面地形，揭示洋流和环流的活动规律也需应用地球重力场数据。因此，本世纪地球重力场反演精度的进一步提高不仅是大地测量、地球物理、地震预报、海洋勘探、空间技术、航空航天等相关学科发展的迫切需求，同时也将为全人类寻求资源、保护环境和预测灾害提供重要的地球空间信息。

在利用重力卫星观测数据反演地球重力场的众多方法中，按引力位系数解算方式的差异可分为空域法和时域法。

1、空域法：不直接处理空间位置相对不规则的卫星轨道采样点的观测值，而将这些观测值归算到以卫星平均轨道高度为半径的球面上利用快速傅立叶变换(FFT)技术进行网格化处理，将问题转化为某类型边值问题的解，常见的准解析法、最小二乘配置法等属于空域法的范畴。优点是因网格点数固定从而方程维数一定，且可以利用快速傅立叶变换技术进行快速批量处理，因此极大地降低了计算量；缺点是在进行网格化处理中作了近似处理，且不能对色噪声进行处理。

2、时域法：将卫星观测数据按时间序列处理，卫星星历值直接表示成引力位系数的函数，由最小二乘等方法直接反求引力位系数。优点是直接对卫星观测数据进行处理，不需作任何近似，求解精度较高且能有效处理色噪声；缺点是随着卫星观测数据的增多，观测方程数量剧增，极大地增加了计算量。过去由于地球重力场反演方法的历史局限性和当时计算机技术发展的限制，为了减少计算量，因此空域法较为盛行。然而，由于空域法做了许多人为性的假设，存在许多潜在的弊端且随着近年来计算机技术的飞速发展及各种快速算法的广泛应用，计算量的大小不再是制约地球重力场反演精度的重要因素，时域法的优点正逐渐体现于地球重力场反演之中。时域法主要包括四种类型：(1)Kaula线性摄动法；(2)基于加速度观测值的数值微分法；(3)基于数值积分的轨道动力学法；(4)基于能量守恒定律的能量法。国内外研究表明，Kaula线性摄动法和基于加速度观测值的数值微分法只适合于求解低阶地球重力场且计算精度较低，现在最为盛行的是轨道动力学法和能量守恒法。轨道动力学法的优点是求解精度较高；缺点是观测数据运算量较大、求解过程复杂程度较高且反演较高阶重力场时需要高性能的并行计算机支持；能量守恒法的优点是观测方程物理含义明确且易于地球重力场的敏感度分析，在保证求解精度的前提下计算量大大降低，通常采用PC计算机可完成高阶地球重力场的快速求解；缺点是对卫星速度的测量精度要求较高。

不同于国内外已有的卫星重力反演法，本发明首次通过在相对轨道位置的星星连线分量中引入星载K波段测量仪(如图2所示)的高精度星间距离，利用新型和精确的星间距离插值卫星重力反演法构建了全球重力场模型。基于美国宇航局喷气推进实验室(NASA-JPL)公布的GRACE Level-1B实测数据建立了新型全球重力场模型WHIGG-GEGM01S(GRACE Earth’s Gravity Model fromWuHan Institute of Geodesy and Geophysics)，进而验证了新型星间距离插值卫星重力反演法的正确性和有效性。由于我国自主研制和正在建设的首期地球重力卫星系统预计于国家“十二五”规划末期发射升空，因此星间距离插值法以其独特的优越性将成为我国高精度和高空间分辨率地球重力场反演的优选方法之一。

三、发明内容

本发明的目的是：由于当前GPS全球定位系统的卫星轨道位置精度相对较低，通过将星载K波段测量仪的精确星间距离引入双星相对轨道位置矢量的星星连线分量构建新型星间距离插值卫星观测方程，进而快速建立高精度和高空间分辨率的全球重力场模型。

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

1、一种利用星间距离插值建立全球重力场模型的方法，包含下列步骤：

步骤一：对GRACE卫星观测数据进行预处理，具体包括

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

1.2)采集星载双频GPS接收机得到的卫星轨道数据，包括轨道位置r和轨道速度

：为了保证卫星轨道数据的精度和连续性，去除卫星轨道存在的重叠期，进行卫星轨道数据的拼接；截掉由于定轨弱约束造成的卫星轨道数据的开始和结束时段处精度较低的数据；基于3σ准则即莱以特准则，剔除卫星轨道数据中存在的粗大误差。

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

步骤二：构建星间距离插值观测方程

在地心惯性坐标系中，基于Newton插值模型，单星轨道位置r的泰勒展开表示如下

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

其中，

表示二项式系数，

t表示插值点的时间，t₀表示插值点的初始时刻，Δt表示采样间隔，n表示插值点的个数。

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

的展开公式

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

基于(2)式，双星轨道加速度差分的展开公式表示如下

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

其中，r₁₂＝r₂-r₁和

分别表示双星相对轨道位置矢量和相对轨道加速度矢量，r₁和r₂分别表示双星绝对轨道位置矢量，

和

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

将(3)式中的

投影到星星连线方向可得

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

其中，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A卫星指向GRACE-B卫星的单位矢量。e₁₂(t)·r₁₂(t_ξ)可改写为

$e_{12}\left(t\right)\·r_{12}\left(t_{\mathrm{\ξ}}\right)=e_{12}\left(t\right)\·[r_{12}^{\left|\right|}\left(t_{\mathrm{\ξ}}\right)+r_{12}^{\⊥}\left(t_{\mathrm{\ξ}}\right)],---\left(5\right)$

其中， $r_{12}^{\left|\right|}\left(t_{\mathrm{\ξ}}\right)=(r_{12}\·e_{12})e_{12}$ 表示r₁₂的星星连线方向分量； $r_{12}^{\⊥}\left(t_{\mathrm{\ξ}}\right)=r_{12}-(r_{12}\·e_{12})e_{12}$ 表示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)式可得

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

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

在(7)式中，

的具体形式表示如下

${\stackrel{\·\·}{r}}_{12}=g_{12}^0+g_{12}^T+a_{12}+f_{12},---\left(8\right)$

其中，

表示作用于双星的相对地球扰动引力，▽表示梯度算子；a₁₂＝a₂-a₁表示除地球引力之外的作用于双星的相对保守力；f₁₂＝f₂-f₁表示作用于双星的相对非保守力；

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

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

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

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

将(8)式和(9)式代入(7)式，星间距离插值观测方程表示如下

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

$+\mathrm{GM}[\frac{r_2\left(t\right)}{{\left|r_2{\left(t\right)}_2\right|}^3}-\frac{r_1\left(t\right)}{{\left|r_1{\left(t\right)}_1\right|}^3}]-a_{12}\left(t\right)-f_{12}\left(t\right)\},---\left(10\right)$

其中，▽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(11\right)$

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

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

和

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

将(10)式按照泰勒公式展开，3点、5点、7点和9点星间距离插值公式表示如下

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

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\},---\left(12\right)$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[-\frac{1}{12}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{4}{3}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{5}{2}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{4}{3}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{1}{12}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)],---\left(13\right)$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[\frac{1}{90}{r}_{\mathrm{\ρ}12}\left(t_{j-3}\right)-\frac{3}{20}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{3}{2}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{49}{18}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{3}{2}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{3}{20}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)+\frac{1}{90}{r}_{\mathrm{\ρ}12}\left(t_{j+3}\right)],---\left(14\right)$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[-\frac{1}{560}{r}_{\mathrm{\ρ}12}\left(t_{j-4}\right)+\frac{8}{315}{r}_{\mathrm{\ρ}12}\left(t_{j-3}\right)-\frac{1}{5}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{8}{5}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{205}{72}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{8}{5}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{1}{5}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)+\frac{8}{315}{r}_{\mathrm{\ρ}12}\left(t_{j+3}\right)-\frac{1}{560}{r}_{\mathrm{\ρ}12}\left(t_{j+4}\right)]$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}.---\left(15\right)$

步骤三：优选不同点数星间距离插值公式

利用“步骤一”中获取的GRACE星载K波段测量仪的星间距离、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，基于“步骤二”中3点、5点、7点和9点星间距离插值公式(12)～(15)，分别计算获得地球引力位系数

和

进而反演地球重力场；结果表明9点星间距离插值公式(15)有利于120阶GRACE卫星重力反演精度的有效提高。

步骤四：建立全球重力场模型

利用“步骤一”中获取的GRACE星载K波段测量仪的星间距离、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，基于“步骤二”中建立的9点星间距离插值公式(15)，计算获得地球引力位系数

和

最终通过引力位系数的集合建立120阶全球重力场模型；在120阶处，大地水准面累积误差为10.98cm，重力异常累积误差为1.741×10^-6m/s²。

步骤五：验证新型全球重力场模型的正确性和可靠性

新型全球重力场模型WHIGG-GEGM01S的大地水准面标准误差RMS(Root-Mean-Square)＝0.726m，更接近于国际已公布的全球重力场模型EIGEN-GRACE02S的大地水准面标准误差RMS＝0.735m。

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

四、附图说明

图1表示GRACE-A/B重力双星在轨相互跟踪飞行示意图。

图2表示GRACE星载K波段测量仪的高精度星间测距原理。

图3表示基于3点、5点、7点和9点新型星间距离插值公式分别反演地球重力场精度，其中横坐标表示地球引力位按球函数展开的阶数，纵坐标表示大地水准面累积误差(单位：m)。

图4表示基于新建立的全球重力场模型WHIGG-GEGM01S绘制的全球大地水准面高分布图，其中横坐标表示经度(0°～360°)，纵坐标表示纬度(-90°～90°)，颜色条表示全球大地水准面的高度起伏(单位：m)。

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

五、具体实施方式

以下结合附图，对本发明作进一步的详细说明。

利用星间距离插值建立全球重力场模型的方法包含下列步骤，其中，全球重力场模型是指地球引力位按球谐函数展开系数的集合

步骤一：对GRACE(Gravity Recovery And Climate Experiment)(即重力场恢复与气候实验)卫星观测数据进行预处理

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波段测量仪(如图2所示)的高精度星间距离引入双星相对轨道位置矢量的星星连线分量，构建新型和精确的星间距离插值卫星观测方程。

在地心惯性坐标系中，基于Newton插值模型，单星轨道位置r的泰勒展开表示如下

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

其中，

表示二项式系数，

t表示插值点的时间，t₀表示插值点的初始时刻，Δt表示采样间隔，n表示插值点的个数。

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

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

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

的展开公式表示如下

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

其中，r₁₂＝r₂-r₁和

分别表示双星相对轨道位置矢量和相对轨道加速度矢量，r₁和r₂分别表示双星绝对轨道位置矢量，

和

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

将(3)式中的

投影到星星连线方向(LOS)可得

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

其中，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A卫星指向GRACE-B卫星的单位矢量。

由于当前GPS定轨精度较低，因此在卫星观测方程中直接使用r₁₂将无法实质性提高地球重力场反演的精度。因此，星载K波段测量仪高精度星间距离观测值的使用是进一步提高地球重力场反演精度的有效途径。e₁₂(t)·r₁₂(t_ξ)可改写为

$e_{12}\left(t\right)\·r_{12}\left(t_{\mathrm{\ξ}}\right)=e_{12}\left(t\right)\·[r_{12}^{\left|\right|}\left(t_{\mathrm{\ξ}}\right)+r_{12}^{\⊥}\left(t_{\mathrm{\ξ}}\right)],---\left(5\right)$

其中， $r_{12}^{\left|\right|}\left(t_{\mathrm{\ξ}}\right)=(r_{12}\·e_{12})e_{12}$ 表示r₁₂的星星连线方向分量； $r_{12}^{\⊥}\left(t_{\mathrm{\ξ}}\right)=r_{12}-(r_{12}\·e_{12})e_{12}$ 表示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)式可得

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

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

在(7)式中，

的具体形式表示如下

${\stackrel{\·\·}{r}}_{12}=g_{12}^0+g_{12}^T+a_{12}+f_{12},---\left(8\right)$

其中，

表示作用于双星的相对地球扰动引力，▽表示梯度算子；a₁₂＝a₂-a₁表示除地球引力之外的作用于双星的相对保守力，包括太阳引力、月球引力、地球固体、海洋、大气和极潮摄动力、以及广义相对论效应摄动力，可通过国际已公布的太阳引力模型、月球引力模型、地球固体、海洋、大气和极潮摄动力模型、以及广义相对论效应摄动力模型获得；f₁₂＝f₂-f₁表示作用于双星的相对非保守力，包括大气阻力、太阳光压、地球辐射压力、卫星轨道高度及姿态控制力、以及经验摄动力；

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

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

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

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

将(8)式和(9)式代入(7)式，星间距离插值方程表示如下

$e_{12}\left(t\right)\·\▿T_{12}\left(t\right)=e_{12}\left(t\right)\·\{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\α}\\ j\end{array}\right)}^{\′\′}\underset{\mathrm{\ξ}=0}{\overset{j}{\mathrm{\Σ}}}{(-1)}^{j+\mathrm{\ξ}}\left(\begin{array}{c}j\\ \mathrm{\ξ}\end{array}\right)r_{\mathrm{\ρ}12}\left(t_{\mathrm{\ξ}}\right),---\left(10\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}]-a_{12}\left(t\right)-f_{12}\left(t\right)\}$

其中，▽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(11\right)$

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

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

和

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

将(10)式按照泰勒公式展开，3点、5点、7点和9点星间距离插值公式表示如下

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

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\},---\left(12\right)$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[-\frac{1}{12}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{4}{3}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{5}{2}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{4}{3}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{1}{12}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)],---\left(13\right)$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[\frac{1}{90}{r}_{\mathrm{\ρ}12}\left(t_{j-3}\right)-\frac{3}{20}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{3}{2}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{49}{18}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{3}{2}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{3}{20}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)+\frac{1}{90}{r}_{\mathrm{\ρ}12}\left(t_{j+3}\right)],---\left(14\right)$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[-\frac{1}{560}{r}_{\mathrm{\ρ}12}\left(t_{j-4}\right)+\frac{8}{315}{r}_{\mathrm{\ρ}12}\left(t_{j-3}\right)-\frac{1}{5}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{8}{5}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{205}{72}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{8}{5}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{1}{5}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)+\frac{8}{315}{r}_{\mathrm{\ρ}12}\left(t_{j+3}\right)-\frac{1}{560}{r}_{\mathrm{\ρ}12}\left(t_{j+4}\right)],---\left(15\right)$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\}$

步骤三：优选不同点数星间距离插值公式

如图3所示，利用“步骤一”中获取的GRACE星载K波段测量仪的星间距离、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，基于“步骤二”中3点、5点、7点和9点星间距离插值公式(12)～(15)，分别计算获得地球引力位系数

和

进而反演地球重力场；结果表明：基于较优的信噪比，9点星间距离插值公式(15)有利于120阶GRACE卫星重力反演精度的有效提高。第一，在120阶内，基于3点星间距离插值公式反演地球重力场的精度远低于分别基于5点、7点和9点星间距离插值公式的反演精度。原因分析如下：首先，由于(12)式的左边是点域值，而右边是平均值，因此(12)式的左右两边几乎不相等；其次，由于3点星间距离插值公式的插值点数较少，因此无法提供足够的地球重力场反演插值信息。第二，在120阶内，基于9点星间距离插值公式反演地球重力场的精度高于分别基于3点、5点和7点星间距离插值公式的反演精度。原因分析如下：由于随着插值点数的增多，卫星观测值的信息量逐渐增加，因此基于9点星间距离插值公式反演地球重力场的精度高于分别基于3点、5点和7点星间距离插值公式的反演精度。综上所述，9点星间距离插值公式是有效提高120阶地球重力场反演精度的较优选择。

步骤四：建立新型全球重力场模型WHIGG-GEGM01S

如图4所示，利用“步骤一”中获取的GRACE星载K波段测量仪的星间距离、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，基于“步骤二”中建立的9点星间距离插值公式(15)，计算获得地球引力位系数

和

最终通过引力位系数的集合建立120阶新型全球重力场模型WHIGG-GEGM01S。图5中实线(大地水准面累积误差)和虚线(重力异常累积误差)分别表示利用预处理的GRACE星载K波段测量仪的星间距离、星载GPS接收机的卫星轨道位置和卫星轨道速度、以及星载加速度计的卫星非保守力观测数据，基于本发明构建的新型和精确9点星间距离插值法，反演GRACE地球重力场的精度；在120阶处大地水准面累积误差为10.98cm，重力异常累积误差为1.741×10^-6m/s²。

步骤五：验证新型全球重力场模型WHIGG-GEGM01S的正确性和可靠性

结合美国、欧洲和澳大利亚的GPS/Levelling数据以及国际已公布的GRACE全球重力场模型，验证建立的新型全球重力场模型WHIGG-GEGM01S的正确性和可靠性。具体计算过程如下：(1)基于本发明新建立的全球重力场模型WHIGG-GEGM01S计算获得全球大地水准面高N₁(经度、纬度和正高)；(2)基于国际已公布的GPS/Levelling数据计算获得全球大地水准面高N₂；(3)通过ΔN＝N₂-N₁评价全球重力场模型WHIGG-GEGM01S的正确性和可靠性。研究结果表明：相对于其它全球重力场模型EIGEN-GRACE01S(标准误差RMS＝0.851m)、EIGEN-CG03C(RMS＝0.481m)、EIGEN-GL04S1(RMS＝0.393m)和EIGEN-5C(RMS＝0.346m)，本发明新建立的全球重力场模型WHIGG-GEGM01S的大地水准面标准误差RMS(Root-Mean-Square)＝0.726m，更接近于已公布全球重力场模型EIGEN-GRACE02S(RMS＝0.735m)的大地水准面标准误差。综上所述，本发明新建立的全球重力场模型WHIGG-GEGM01S是正确的和可靠的。新型星间距离插值法是反演高精度和高空间分辨率地球重力场的优化方法之一。

Claims (1)

1.一种利用星间距离插值建立全球重力场模型的方法，包含下列步骤：

步骤一：对GRACE卫星观测数据进行预处理，具体包括

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

1.2)采集星载双频GPS接收机得到的卫星轨道数据，包括轨道位置(r)和轨道速度

为了保证卫星轨道数据的精度和连续性，去除卫星轨道存在的重叠期，进行卫星轨道数据的拼接；截掉由于定轨弱约束造成的卫星轨道数据的开始和结束时段处精度较低的数据；基于3σ准则即莱以特准则，剔除卫星轨道数据中存在的粗大误差；

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

步骤二：构建星间距离插值观测方程

在地心惯性坐标系中，基于Newton插值模型，单星轨道位置r的泰勒展开表示如下

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

其中，

表示二项式系数，

t表示插值点的时间，t₀表示插值点的初始时刻，Δt表示采样间隔，n表示插值点的个数；

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

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

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

的展开公式表示如下

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

其中，r₁₂＝r₂-r₁和

分别表示双星相对轨道位置矢量和相对轨道加速度矢量，r₁和r₁分别表示双星绝对轨道位置矢量，

和

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

将(3)式中的投影到星星连线方向(LOS)可得

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

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

$e_{12}\left(t\right)\·r_{12}\left(t_{\mathrm{\ξ}}\right)=e_{12}\left(t\right)\·[r_{12}^{\left|\right|}\left(t_{\mathrm{\ξ}}\right)+r_{12}^{\⊥}\left(t_{\mathrm{\ξ}}\right)],---\left(5\right)$

其中， $r_{12}^{\left|\right|}\left(t_{\mathrm{\ξ}}\right)=(r_{12}\·e_{12})e_{12}$ 表示r₁₂的星星连线方向分量； $r_{12}^{\⊥}\left(t_{\mathrm{\ξ}}\right)=r_{12}-(r_{12}\·e_{12})e_{12}$ 表示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)式可得

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

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

在(7)式中，

的具体形式表示如下

${\stackrel{\·\·}{r}}_{12}=g_{12}^0+g_{12}^T+a_{12}+f_{12},---\left(8\right)$

其中，

表示作用于双星的相对地球扰动引力，▽表示梯度算子；a₁₂＝a₂-a₁表示除地球引力之外的作用于双星的相对保守力；f₁₂＝f₂-f₁表示作用于双星的相对非保守力；

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

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

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

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

将(8)式和(9)式代入(7)式，星间距离插值观测方程表示如下

$e_{12}\left(t\right)\·\▿T_{12}\left(t\right)=e_{12}\left(t\right)\·\{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(\begin{array}{c}\mathrm{\α}\\ j\end{array}\right)}^{\′\′}\underset{\mathrm{\ξ}=0}{\overset{j}{\mathrm{\Σ}}}{(-1)}^{j+\mathrm{\ξ}}\left(\begin{array}{c}j\\ \mathrm{\ξ}\end{array}\right)r_{\mathrm{\ρ}12}\left(t_{\mathrm{\ξ}}\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}]-a_{12}\left(t\right)-f_{12}\left(t\right)\},---\left(10\right)$

其中，▽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(11\right)$

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

和

表示待求的规格化引力位系数；

将(10)式按照泰勒公式展开，3点、5点、7点和9点星间距离插值公式表示如下

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

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\},---\left(12\right)$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[-\frac{1}{12}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{4}{3}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{5}{2}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{4}{3}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{1}{12}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)]$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\},---\left(13\right)$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[\frac{1}{90}{r}_{\mathrm{\ρ}12}\left(t_{j-3}\right)-\frac{3}{20}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{3}{2}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{49}{18}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{3}{2}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{3}{20}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)+\frac{1}{90}{r}_{\mathrm{\ρ}12}\left(t_{j+3}\right)]$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\},---\left(14\right)$

$e_{12}\left(t_j\right)\·\▿T_{12}\left(t_j\right)=e_{12}\left(t_j\right)\·\left\{\frac{1}{{\left(\mathrm{\Δt}\right)}^2}\right[-\frac{1}{560}{r}_{\mathrm{\ρ}12}\left(t_{j-4}\right)+\frac{8}{315}{r}_{\mathrm{\ρ}12}\left(t_{j-3}\right)-\frac{1}{5}{r}_{\mathrm{\ρ}12}\left(t_{j-2}\right)+\frac{8}{5}{r}_{\mathrm{\ρ}12}\left(t_{j-1}\right)$

$-\frac{205}{72}{r}_{\mathrm{\ρ}12}\left(t_j\right)+\frac{8}{5}{r}_{\mathrm{\ρ}12}\left(t_{j+1}\right)-\frac{1}{5}{r}_{\mathrm{\ρ}12}\left(t_{j+2}\right)+\frac{8}{315}{r}_{\mathrm{\ρ}12}\left(t_{j+3}\right)-\frac{1}{560}{r}_{\mathrm{\ρ}12}\left(t_{j+4}\right)]$

$+\mathrm{GM}[\frac{r_2\left(t_j\right)}{{\left|r_2\left(t_j\right)\right|}^3}-\frac{r_1\left(t_j\right)}{{\left|r_1\left(t_j\right)\right|}^3}]-a_{12}\left(t_j\right)-f_{12}\left(t_j\right)\},---\left(15\right)$

步骤三：优选不同点数星间距离插值公式

利用“步骤一”中获取的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。

Images