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

摘要

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

说明书

一、技术领域

本发明涉及卫星大地测量学、地球物理学、空间科学等交叉技术领域，特 别是涉及一种通过将星载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 from  WuHan 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}{\Sigma }}\left(\begin{array}{c}\alpha \\ j\end{array}\right)\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r\left(t_{\xi }\right),---\left(1\right)$

其中，表示二项式系数，t表示插值点的时间，t₀表示插值点的 初始时刻，Δt表示采样间隔，n表示插值点的个数。

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

$\stackrel{··}{r}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r\left(t_{\xi }\right).---\left(2\right)$

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

${\stackrel{··}{r}}_{12}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r_{12}\left(t_{\xi }\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}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)e_{12}\left(t\right)·r_{12}\left(t_{\xi }\right).---\left(4\right)$

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

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

其中，$r_{12}^{\left|\right|}\left(t_{\xi }\right)=(r_{12}·e_{12})e_{12}$表示r₁₂的星星连线方向分量；$r_{12}^{\perp }\left(t_{\xi }\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}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)e_{12}\left(t\right)·r_{\rho 12}\left(t_{\xi }\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}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r_{\rho 12}\left(t_{\xi }\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,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\Sigma }}({\bar{C}}_{\mathrm{lm}}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\cos \theta \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{\Delta t}\right)}^2}\right[r_{\rho 12}\left(t_{j-1}\right)-{2r}_{\rho 12}\left(t_j\right)+r_{\rho 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{\Delta t}\right)}^2}\right[-\frac{1}{12}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{4}{3}{r}_{\rho 12}\left(t_{j-1}\right)$

$-\frac{5}{2}{r}_{\rho 12}\left(t_j\right)+\frac{4}{3}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{1}{12}{r}_{\rho 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{\Delta t}\right)}^2}\right[\frac{1}{90}{r}_{\rho 12}\left(t_{j-3}\right)-\frac{3}{20}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{3}{2}{r}_{\rho 12}\left(t_{j-1}\right)$

$-\frac{49}{18}{r}_{\rho 12}\left(t_j\right)+\frac{3}{2}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{3}{20}{r}_{\rho 12}\left(t_{j+2}\right)+\frac{1}{90}{r}_{\rho 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{\Delta t}\right)}^2}\right[-\frac{1}{560}{r}_{\rho 12}\left(t_{j-4}\right)+\frac{8}{315}{r}_{\rho 12}\left(t_{j-3}\right)-\frac{1}{5}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{8}{5}{r}_{\rho 12}\left(t_{j-1}\right)$

$-\frac{205}{72}{r}_{\rho 12}\left(t_j\right)+\frac{8}{5}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{1}{5}{r}_{\rho 12}\left(t_{j+2}\right)+\frac{8}{315}{r}_{\rho 12}\left(t_{j+3}\right)-\frac{1}{560}{r}_{\rho 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}{\Sigma }}\left(\begin{array}{c}\alpha \\ j\end{array}\right)\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r\left(t_{\xi }\right),---\left(1\right)$

其中，表示二项式系数，t表示插值点的时间，t₀表示插值点的 初始时刻，Δt表示采样间隔，n表示插值点的个数。

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

$\stackrel{··}{r}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r\left(t_{\xi }\right).---\left(2\right)$

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

${\stackrel{··}{r}}_{12}\left(t\right)=\underset{j=1}{\overset{n}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r_{12}\left(t_{\xi }\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}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)e_{12}\left(t\right)·r_{12}\left(t_{\xi }\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_{\xi }\right)=e_{12}\left(t\right)·[r_{12}^{\left|\right|}\left(t_{\xi }\right)+r_{12}^{\perp }\left(t_{\xi }\right)],---\left(5\right)$

其中，$r_{12}^{\left|\right|}\left(t_{\xi }\right)=(r_{12}·e_{12})e_{12}$表示r₁₂的星星连线方向分量；$r_{12}^{\perp }\left(t_{\xi }\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}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)e_{12}\left(t\right)·r_{\rho 12}\left(t_{\xi }\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}{\Sigma }}{\left(\begin{array}{c}\alpha \\ j\end{array}\right)}^{\prime \prime }\underset{\xi =0}{\overset{j}{\Sigma }}{(-1)}^{j+\xi }\left(\begin{array}{c}j\\ \xi \end{array}\right)r_{\rho 12}\left(t_{\xi }\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,\theta ,\lambda )=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\Sigma }}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\Sigma }}({\bar{C}}_{\mathrm{lm}}\cos \mathrm{m\lambda }+{\bar{S}}_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\cos \theta \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{\Delta t}\right)}^2}\right[r_{\rho 12}\left(t_{j-1}\right)-{2r}_{\rho 12}\left(t_j\right)+r_{\rho 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{\Delta t}\right)}^2}\right[-\frac{1}{12}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{4}{3}{r}_{\rho 12}\left(t_{j-1}\right)$

$-\frac{5}{2}{r}_{\rho 12}\left(t_j\right)+\frac{4}{3}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{1}{12}{r}_{\rho 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{\Delta t}\right)}^2}\right[\frac{1}{90}{r}_{\rho 12}\left(t_{j-3}\right)-\frac{3}{20}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{3}{2}{r}_{\rho 12}\left(t_{j-1}\right)$

$-\frac{49}{18}{r}_{\rho 12}\left(t_j\right)+\frac{3}{2}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{3}{20}{r}_{\rho 12}\left(t_{j+2}\right)+\frac{1}{90}{r}_{\rho 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{\Delta t}\right)}^2}\right[-\frac{1}{560}{r}_{\rho 12}\left(t_{j-4}\right)+\frac{8}{315}{r}_{\rho 12}\left(t_{j-3}\right)-\frac{1}{5}{r}_{\rho 12}\left(t_{j-2}\right)+\frac{8}{5}{r}_{\rho 12}\left(t_{j-1}\right)$

$-\frac{205}{72}{r}_{\rho 12}\left(t_j\right)+\frac{8}{5}{r}_{\rho 12}\left(t_{j+1}\right)-\frac{1}{5}{r}_{\rho 12}\left(t_{j+2}\right)+\frac{8}{315}{r}_{\rho 12}\left(t_{j+3}\right)-\frac{1}{560}{r}_{\rho 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是 正确的和可靠的。新型星间距离插值法是反演高精度和高空间分辨率地球重力 场的优化方法之一。