CN103076639A - 基于残余星间速度原理反演地球重力场的方法 - Google Patents

Abstract

本发明涉及一种基于残余星间速度原理反演地球重力场的方法；通过将星载激光干涉测距仪的高精度残余星间速度观测量引入GPS接收机的残余轨道速度差分矢量的视线分量建立新型残余星间速度观测方程，进而精确和快速反演地球重力场的方法；该方法地球重力场计算精度高，卫星重力反演速度快，计算机性能要求低，敏感于重力场中高频信号，易于卫星重力反演误差分析；残余星间速度法是建立高精度和高空间分辨率全球重力场模型的关键技术。

Description

基于残余星间速度原理反演地球重力场的方法

一、技术领域

本发明涉及卫星重力学、空间大地测量学、地球动力学、宇航学等交叉技术领域，特别是涉及一种通过将星载激光干涉测距仪的高精度残余星间速度观测量（测量精度10^-9m/s）引入GPS接收机的残余轨道速度差分矢量的视线分量建立新型残余星间速度观测方程，进而精确和快速建立高阶次地球重力场模型的方法。

二、背景技术

地球重力场及其时变反映地球表层及内部物质的空间分布、运动和变化，同时决定着大地水准面的起伏和变化。因此，确定地球重力场的精细结构及其时变不仅是大地测量学、地球物理学、地球动力学、海洋学、空间科学等的需求，同时也将为寻求资源、保护环境和预测灾害提供重要的信息资源。基于GRACE（Gravity Recovery and Climate Experiment）卫星精密探测地球中长波静态和长波时变重力场的突出贡献，美国宇航局（NASA）提出了用于高精度探测地球中短波静态和中长波时变重力场的GRACE Follow-On卫星重力测量计划。GRACE Follow-On双星采用近圆（轨道离心率0.001）、近极（轨道倾角89°）和低轨道（轨道高度250km）设计，利用激光干涉测距仪精确测量星间距离（10^-8m）和星间速度（10^-9m/s），基于高轨GPS卫星确定轨道位置和轨道速度，通过星载加速度计感测作用于卫星的非保守力（10^-13m/s²），以及依靠非保守力补偿系统（DFCS）平衡非保守力（大气阻力、太阳光压、地球辐射压、轨道高度和姿态控制力等）。基于较低的卫星轨道高度和较高的关键载荷测量精度，利用GRACE Follow-On计划反演地球重力场精度较当前GRACE计划至少提高10倍。

Wolff于1969年首次提出了利用卫星跟踪卫星低低技术（SST-LL）测量地球重力场的新思想。自此以后，国内外众多学者积极投身于地球重力场反演的技术研究和数值计算之中。地球重力场模型是指地球引力位系数的集合当前，国际通用的卫星重力反演方法主要包括：动力学法、能量守恒法、卫星加速度法、解析/半解析法等。美国宇航局喷气推进实验室（NASA-JPL）、美国德克萨斯州立大学空间研究中心（CSR）、以及德国波兹坦地学研究中心（GFZ）等研究机构采用动力学法建立了全球重力场模型GGM01S/02S，EIGEN-GRACE01S/02S，EIGEN-CG01C/03C，EIGEN-GL04C/04S1，EIGEN-5C等。但动力学法的主要缺点是随着卫星轨道弧长的增加，动力模型的误差将快速累积，而且计算过程较复杂，需要高性能并行计算机支持。

为了有效克服动力学法的缺点以及进一步提高地球重力场的反演精度，本发明首次通过将星载激光干涉测距仪的高精度残余星间速度观测量（测量精度10^-9m/s）引入GPS接收机的残余轨道速度差分矢量的视线方向建立了新型残余星间速度法，进而精确和快速反演了GRACE Follow-On地球重力场。

三、发明内容

本发明的目的是：通过将星载激光干涉测距仪的高精度残余星间速度观测量引入GPS接收机的残余轨道速度差分矢量的视线分量，建立新型残余星间速度观测方程，进而精确和快速建立高阶次地球重力场模型的方法。

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

基于残余星间速度原理反演地球重力场的方法包含下列步骤：

步骤一：卫星观测数据采集

（1）通过星载激光干涉测距仪获取星间速度

通过星载GPS接收机获取双星轨道位置（r₁,r₂）和双星轨道速度通过星载加速度计获取作用于双星的非保守力（f₁,f₂）。

（2）利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置

和双星参考轨道速度

（3）参考星间速度

通过参考轨道速度计算获得

其中，

表示相对参考轨道速度矢量，

表示第一颗参考卫星指向第二颗参考卫星的参考单位矢量，

表示相对参考轨道位置矢量，

和

表示双星的参考轨道位置矢量。

（4）参考非保守力

通过DTM2000阻力温度模型计算获得。

（5）通过国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于双星的保守力（F₁,F₂）和参考保守力

步骤二：残余星间速度观测方程建立

在地心惯性系中，基于牛顿插值原理，单星轨道速度

的泰勒展开表示如下

$\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)$

其中， $\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)$ 表示二项式系数，

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

单星参考轨道速度

的泰勒展开表示如下

基于公式（1）-（2），单星残余轨道速度

的泰勒展开表示如下

$\mathrm{\δ}\stackrel{\·}{r}\left(t\right)=\mathrm{\δ}\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)\mathrm{\δ}\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right)---\left(3\right)$

其中，

基于公式（3）的一阶时间导数，单星残余轨道加速度

的泰勒展开表示如下

$\mathrm{\δ}\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)\mathrm{\δ}\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right)---\left(4\right)$

其中，

基于公式（4），双星残余轨道加速度差分

的泰勒展开表示如下

$\mathrm{\δ}{\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)\mathrm{\δ}{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)---\left(5\right)$

其中，

和

分别表示轨道速度差分矢量和轨道加速度差分矢量，

和

表示双星的轨道速度矢量，

和

表示双星的参考轨道速度矢量，

和

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

和

表示双星的参考轨道加速度矢量。

双星残余轨道加速度差分

的视线分量表示如下

$e_{12}\left(t\right)\·\mathrm{\δ}{\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)\·\mathrm{\δ}{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)---\left(6\right)$

其中，e₁₂=r₁₂/|r₁₂|表示第一颗卫星指向第二颗卫星的单位矢量，r₁₂=r₂-r₁表示双星的轨道位置差分矢量，r₁和r₂分别表示双星的轨道位置矢量。

引入激光干涉测距系统的高精度残余星间速度

进一步提高地球重力场精度，其中

和分别表示星间速度和参考星间速度；可被改写为

$\mathrm{\δ}{\stackrel{\·}{r}}_{12}={\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\left|\right|}+{\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\⊥}---\left(7\right)$

其中， ${\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\left|\right|}=(\mathrm{\δ}{\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 和 ${\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\⊥}=\mathrm{\δ}{\stackrel{\·}{r}}_{12}-(\mathrm{\δ}{\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 分别表示

的视线分量和垂向分量。

基于误差传播原理，为了有效降低

将公式（7）中的

替换为

因此，公式（6）可改写为

$e_{12}\left(t\right)\·\mathrm{\δ}{\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)\·\mathrm{\δ}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{\mathrm{\α}}\right)---\left(8\right)$

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

在公式（8）中，取插值点数n=2,4,6,8，可得到2点、4点、6点和8点残余星间速度公式，分别表示如下：

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

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

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

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

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

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

在公式（8）中，

的具体形式表示如下

$\mathrm{\δ}{\stackrel{\·\·}{r}}_{12}=\mathrm{\δ}{g}_{12}+\mathrm{\δ}{T}_{12}+\mathrm{\δ}{F}_{12}+\mathrm{\δ}{f}_{12}---\left(13\right)$

其中，δT₁₂表示作用于双星的残余地球扰动引力差；

表示除地球引力之外的残余保守力差，F₁和F₂表示作用于双星的保守力，

和

表示参考保守力；表示残余非保守力差，f₁和f₂表示作用于双星的非保守力，

和

表示参考非保守力；

表示残余地心引力差，g₁和g₂表示双星的地心引力，

和

表示参考地心引力

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

分别表示双星的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎表示轨道位置矢量r₁₍₂₎的3个分量。

通过将公式（13）和（14）代入（8），残余星间速度观测方程表示如下

$e_{12}\left(t\right)\·\mathrm{\δ}{T}_{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)\·\{\mathrm{\δ}{\stackrel{\·}{\mathrm{\ρ}}}_{12}\left(t_{\mathrm{\α}}\right)e_{12}\left(t_{\mathrm{\α}}\right)$

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

$-e_{12}\left(t\right)\·[\mathrm{\δ}{g}_{12}+\mathrm{\δ}{F}_{12}\left(t\right)+\mathrm{\δ}{f}_{12}\left(t\right)]$

其中，

表示残余地球扰动位的一阶梯度，V₁和V₂表示地球扰动位，

和表示参考地球扰动位

$V(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(16\right)$

其中，r,θ，λ分别表示地心半径、地心余纬度和地心经度，R_e表示地球平均半径；表示正规化的缔合Legendre函数，l表示阶数，m表示次数；

和

表示待估的地球引力位系数。

步骤三：地球重力场反演

基于残余星间速度观测方程（15），利用激光干涉测距仪的星间速度、GPS接收机的轨道位置和轨道速度、加速度计的非保守力、以及双星的参考非保守力、保守力和参考保守力，求解地球引力位系数

和

最终通过地球引力位系数的集合建立全球重力场模型。

本发明是基于新型残余星间速度法有利于精确和快速反演地球重力场的特点而设计的，优点是：1）地球重力场计算精度高；2）卫星重力反演速度快；3）计算机性能要求低；4）敏感于重力场中高频信号；5）易于卫星重力反演误差分析。

四、附图说明

图1表示卫星实测轨道（实线）和参考轨道（虚线）。

图2表示基于不同的相关系数和相同的采样间隔5s模拟的观测值色噪声，(a)星间距离，(b)星间速度。

图3表示基于相关系数0.85和采样间隔5s模拟的观测值色噪声，(a)星间距离，(b)星间速度。

图4表示基于不同的相关系数和相同的采样间隔5s模拟的视线方向色噪声，(a)轨道位置，(b)轨道速度。

图5表示基于相关系数0.95和采样间隔5s模拟的视线方向色噪声，(a)轨道位置，(b)轨道速度。

图6表示基于相关系数0.90和采样间隔5s模拟的视线方向非保守力色噪声。

图7表示基于不同的星间速度插值点数反演累计大地水准面精度。

图8表示基于6点残余星间速度法以及利用不同的相关系数和相同的采样间隔5s反演地球引力位系数精度，（a）星间速度，（b）轨道位置和轨道速度，（c）非保守力。

图9表示基于不同关键载荷精度指标反演地球引力位系数精度。

图10表示基于残余星间速度法反演GRACE Follow-On累计大地水准面精度。

五、具体实施方式

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

实施例一：

基于残余星间速度原理反演地球重力场的方法包含下列步骤：

步骤一：卫星观测数据采集

（1）通过星载激光干涉测距仪获取星间速度

通过星载GPS接收机获取双星轨道位置（r₁,r₂）和双星轨道速度

通过星载加速度计获取作用于双星的非保守力（f₁,f₂）。

（2）如图1所示，利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置

和双星参考轨道速度

（3）参考星间速度

通过参考轨道速度

计算获得

其中，表示相对参考轨道速度矢量，

表示第一颗参考卫星指向第二颗参考卫星的参考单位矢量，

表示相对参考轨道位置矢量，和表示双星的参考轨道位置矢量。

（4）参考非保守力通过DTM2000阻力温度模型计算获得。

（5）通过国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于双星的保守力（F₁,F₂）和参考保守力

其中步骤（4）、（5）的计算方法已在【郑伟,许厚泽,钟敏,员美娟,周旭华,彭碧波.卫星跟踪卫星测量模式中星载加速度计高低灵敏轴分辨率指标优化设计论证.地球物理学报,2009,52(11):27122720.】和【Tapley B,Ries J,Bettadpur S,Chambers D,Cheng M,Condi F,Gunter B,Kang Z,Nagel P，Pastor R,Pekker T,Poole S,Wang F.GGM02-An improved Earth gravity field model fromGRACE.Journal of Geodesy,2005,79(8):467478.】中公开。

步骤二：残余星间速度观测方程建立

在地心惯性系中，基于牛顿插值原理，单星轨道速度

的泰勒展开表示如下

$\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(17\right)$

其中， $\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)$ 表示二项式系数，

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

单星参考轨道速度

的泰勒展开表示如下

基于公式（17）-（18），单星残余轨道速度

的泰勒展开表示如下

$\mathrm{\δ}\stackrel{\·}{r}\left(t\right)=\mathrm{\δ}\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)\mathrm{\δ}\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right)---\left(19\right)$

其中，

基于公式（19）的一阶时间导数，单星残余轨道加速度

的泰勒展开表示如下

$\mathrm{\δ}\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)\mathrm{\δ}\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right)---\left(20\right)$

其中，

基于公式（20），双星残余轨道加速度差分

的泰勒展开表示如下

$\mathrm{\δ}{\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)\mathrm{\δ}{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)---\left(21\right)$

其中，

和

分别表示轨道速度差分矢量和轨道加速度差分矢量，

和

表示双星的轨道速度矢量，

和

表示双星的参考轨道速度矢量，

和

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

和

表示双星的参考轨道加速度矢量。

双星残余轨道加速度差分

的视线分量表示如下

$e_{12}\left(t\right)\·\mathrm{\δ}{\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)\·\mathrm{\δ}{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)---\left(22\right)$

其中，e₁₂=r₁₂/|r₁₂|表示第一颗卫星指向第二颗卫星的单位矢量，r₁₂=r₂-r₁表示双星的轨道位置差分矢量，r₁和r₂分别表示双星的轨道位置矢量。

引入激光干涉测距系统的高精度残余星间速度进一步提高地球重力场精度，其中

和

分别表示星间速度和参考星间速度；

可被改写为

$\mathrm{\δ}{\stackrel{\·}{r}}_{12}={\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\left|\right|}+{\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\⊥}---\left(23\right)$

其中， ${\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\left|\right|}=(\mathrm{\δ}{\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 和 ${\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\⊥}=\mathrm{\δ}{\stackrel{\·}{r}}_{12}-(\mathrm{\δ}{\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 分别表示

的视线分量和垂向分量。

基于误差传播原理，为了有效降低

将公式（7）中的

替换为

因此，公式（22）可改写为

$e_{12}\left(t\right)\·\mathrm{\δ}{\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)\·\mathrm{\δ}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{\mathrm{\α}}\right)---\left(24\right)$

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

在公式（24）中，取插值点数n=2,4,6,8，得到2点、4点、6点和8点残余星间速度公式，分别表示如下：

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

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

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

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

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

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

在公式（24）中，的具体形式表示如下

$\mathrm{\δ}{\stackrel{\·\·}{r}}_{12}=\mathrm{\δ}{g}_{12}+\mathrm{\δ}{T}_{12}+\mathrm{\δ}{F}_{12}+\mathrm{\δ}{f}_{12}---\left(29\right)$

其中，δT₁₂表示作用于双星的残余地球扰动引力差；

表示除地球引力之外的残余保守力差，F₁和F₂表示作用于双星的保守力，

和

表示参考保守力；

表示残余非保守力差，f₁和f₂表示作用于双星的非保守力，

和

表示参考非保守力；表示残余地心引力差，g₁和g₂表示双星的地心引力，

和

表示参考地心引力

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

分别表示双星的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎表示轨道位置矢量r₁₍₂₎的3个分量。

通过将公式（29）和（30）代入（24），残余星间速度观测方程表示如下

$e_{12}\left(t\right)\·\mathrm{\δ}{T}_{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)\·\{\mathrm{\δ}{\stackrel{\·}{\mathrm{\ρ}}}_{12}\left(t_{\mathrm{\α}}\right)e_{12}\left(t_{\mathrm{\α}}\right)$

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

$-e_{12}\left(t\right)\·[\mathrm{\δ}{g}_{12}+\mathrm{\δ}{F}_{12}\left(t\right)+\mathrm{\δ}{f}_{12}\left(t\right)]$

其中，

表示残余地球扰动位的一阶梯度，V₁和V₂表示地球扰动位，

和

表示参考地球扰动位

$V(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(32\right)$

其中，r,θ，λ分别表示地心半径、地心余纬度和地心经度，R_e表示地球平均半径；

表示正规化的缔合Legendre函数，l表示阶数，m表示次数；

和

表示待估的地球引力位系数。

步骤三：地球重力场反演

基于残余星间速度观测方程（31），利用激光干涉测距仪的星间速度、GPS接收机的轨道位置和轨道速度、星载加速度计的非保守力、以及双星的参考非保守力、保守力和参考保守力等数据，求解地球引力位系数和

最终通过地球引力位系数的集合建立全球重力场模型。

本发明是基于新型残余星间速度法有利于精确和快速反演地球重力场的特点而设计的，优点是：1）地球重力场计算精度高；2）卫星重力反演速度快；3）计算机性能要求低；4）敏感于重力场中高频信号；5）易于卫星重力反演误差分析。

实施例二：

由于卫星观测值不是相互独立，而具有明显的相关性。因此，在反演下一代GRACE Follow-On地球重力场时，在卫星观测值中引入了具有相关性的色噪声。

基于残余星间速度原理反演地球重力场的方法包含下列步骤：

步骤一：卫星观测数据采集

（1）通过星载激光干涉测距仪获取星间速度

通过星载GPS接收机获取双星轨道位置（r₁,r₂）和双星轨道速度

通过星载加速度计获取作用于双星的非保守力（f₁,f₂）。

（2）如图1所示，利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置

和双星参考轨道速度

（3）参考星间速度

通过参考轨道速度计算获得

其中，

表示相对参考轨道速度矢量，

表示第一颗参考卫星指向第二颗参考卫星的参考单位矢量，

表示相对参考轨道位置矢量，

和

表示双星的参考轨道位置矢量。

（4）参考非保守力

通过DTM2000阻力温度模型计算获得

（5）通过国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于双星的保守力（F₁,F₂）和参考保守力

其中步骤（4）、（5）的计算方法已在【郑伟,许厚泽,钟敏,员美娟,周旭华,彭碧波.卫星跟踪卫星测量模式中星载加速度计高低灵敏轴分辨率指标优化设计论证.地球物理学报,2009,52(11):2712-2720.】和【Tapley B,Ries J,Bettadpur S,Chambers D,Cheng M,Condi F,Gunter B,Kang Z,Nagel P，Pastor R,Pekker T,Poole S,Wang F.GGM02-An improved Earth gravity field model fromGRACE.Journal of Geodesy,2005,79(8):467478.】中公开。

步骤二：残余星间速度观测方程建立

在地心惯性系（ECI）中，基于牛顿插值原理，单星实测轨道速度

的泰勒展开表示如下

$\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(33\right)$

其中， $\left(\begin{array}{c}\mathrm{\β}\\ i\end{array}\right)$ 表示二项式系数，

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

单星参考轨道速度的泰勒展开表示如下

基于公式（33）-（34），单星残余轨道速度

的泰勒展开表示如下

$\mathrm{\δ}\stackrel{\·}{r}\left(t\right)=\mathrm{\δ}\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)\mathrm{\δ}\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right)---\left(35\right)$

其中，

基于公式（35）的一阶时间导数，单星残余轨道加速度

的泰勒展开表示如下

$\mathrm{\δ}\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)\mathrm{\δ}\stackrel{\·}{r}\left(t_{\mathrm{\α}}\right)---\left(36\right)$

其中，

基于公式（36），双星残余轨道加速度差分

的泰勒展开表示如下

$\mathrm{\δ}{\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)\mathrm{\δ}{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)---\left(37\right)$

其中，

和

分别表示轨道速度差分矢量和轨道加速度差分矢量，

和

表示双星的轨道速度矢量，和

表示双星的参考轨道速度矢量，

和

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

和

表示双星的参考轨道加速度矢量。

双星残余轨道加速度差分

的视线分量表示如下

$e_{12}\left(t\right)\·\mathrm{\δ}{\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)\·\mathrm{\δ}{\stackrel{\·}{r}}_{12}\left(t_{\mathrm{\α}}\right)---\left(38\right)$

其中，e₁₂=r₁₂/|r₁₂|表示第一颗卫星指向第二颗卫星的单位矢量，r₁₂=r₂-r₁表示双星的轨道位置差分矢量，r₁和r₂分别表示双星的轨道位置矢量。

由于GPS低精度的轨道测量，假如

被直接使用于公式（38），地球重力场精度将无法实质性提高。因此，激光干涉测距系统的高精度残余星间速度

的有效引入是进一步提高地球重力场精度的有效手段，其中

和

分别表示实测星间速度和参考星间速度。

可被改写为

$\mathrm{\δ}{\stackrel{\·}{r}}_{12}={\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\left|\right|}+{\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\⊥}---\left(39\right)$

其中， ${\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\left|\right|}=(\mathrm{\δ}{\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 和 ${\mathrm{\δ}\stackrel{\·}{r}}_{12}^{\⊥}=\mathrm{\δ}{\stackrel{\·}{r}}_{12}-(\mathrm{\δ}{\stackrel{\·}{r}}_{12}\·e_{12})e_{12}$ 分别表示

的视线分量和垂向分量。

基于误差传播原理，为了有效降低

将公式（39）中的

替换为

因此，公式（38）可改写为

$e_{12}\left(t\right)\·\mathrm{\δ}{\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)\·\mathrm{\δ}{\stackrel{\·}{r}}_{\mathrm{\ρ}12}\left(t_{\mathrm{\α}}\right)---\left(40\right)$

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

基于公式（40），取插值点数n=2,4,6,8，得到2点、4点、6点和8点残余星间速度公式，分别表示如下

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

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

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

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

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

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

在公式（40）中，

的具体形式表示如下

$\mathrm{\δ}{\stackrel{\·\·}{r}}_{12}=\mathrm{\δ}{g}_{12}+\mathrm{\δ}{T}_{12}+\mathrm{\δ}{F}_{12}+\mathrm{\δ}{f}_{12}---\left(45\right)$

其中，δT₁₂表示作用于双星的残余地球扰动引力差；

表示除地球引力之外的残余保守力差，F₁和F₂表示作用于双星的保守力，

和表示参考保守力；

表示残余非保守力差，f₁和f₂表示作用于双星的非保守力，

和

表示参考非保守力；表示残余地心引力差，g₁和g₂表示双星的地心引力，

和

表示参考地心引力

其中，GM表示地球质量M和万有引力常数G的乘积，分别表示双星的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎表示轨道位置矢量r₁₍₂₎的3个分量。

通过将公式（45）代入（40），残余星间速度观测方程表示如下：

$e_{12}\left(t\right)\·\mathrm{\δ}{T}_{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)\·\{\mathrm{\δ}{\stackrel{\·}{\mathrm{\ρ}}}_{12}\left(t_{\mathrm{\α}}\right)e_{12}\left(t_{\mathrm{\α}}\right)$

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

$-e_{12}\left(t\right)\·[\mathrm{\δ}{g}_{12}+\mathrm{\δ}{F}_{12}\left(t\right)+\mathrm{\δ}{f}_{12}\left(t\right)]$

其中，

表示残余地球扰动位的一阶梯度，V₁和V₂表示地球扰动位，

和表示参考地球扰动位

$V(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(48\right)$

其中，r,θ，λ分别表示地心半径、地心余纬度和地心经度，R_e表示地球平均半径；

表示正规化的缔合Legendre函数，l表示阶数，m表示次数；和

表示待估的地球引力位系数。

本发明的核心技术是将残余轨道速度差分矢量

分解为沿视线方向分量

和垂直方向分量

由于沿视线方向单位矢量e₁₂的点乘作用，垂向分量误差被有效降低。为了有效提高地球重力场反演精度，本发明用高精度的残余星间速度观测量

替换原视线分量。

步骤三：卫星观测值色噪声模拟

基于残余星间速度观测方程（47），利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟了GRACE Follow-On-A/B双星的实测轨道和参考轨道。GRACE Follow-On-A/B双星的模拟参数为先验地球重力场模型EGM2008、观测时间30天和采样间隔5s。

卫星观测值不是相互独立，而具有明显的相关性。因此，在卫星观测值中引入了具有相关性的色噪声。基于Gauss-Markov模型，卫星观测值色噪声公式表示如下

$\left\{\begin{array}{c}{\mathrm{\ϵ}}_0={\mathrm{\δ}}_0\\ {\mathrm{\ϵ}}_1=\mathrm{\μ}{\mathrm{\ϵ}}_0+\sqrt{1-{\mathrm{\μ}}^2}{\mathrm{\δ}}_1\\ {\mathrm{\ϵ}}_2=\mathrm{\μ}{\mathrm{\ϵ}}_1+\sqrt{1-{\mathrm{\μ}}^2}{\mathrm{\δ}}_2\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_i=\mathrm{\μ}{\mathrm{\ϵ}}_{i-1}+\sqrt{1-{\mathrm{\μ}}^2}{\mathrm{\δ}}_i\end{array}\right.---\left(49\right)$

其中，μ表示相关系数；δ_i(i=1,2,…,k)表示正态分布的随机白噪声（μ=0），k表示卫星观测值的数量；ε_i(i＝1,2,…,k)表示色噪声（0<μ<1）。

3.1）星间距离和星间速度的色噪声

如图2(a)所示，基于不同的相关系数0≤μ≤0.99和相同的采样间隔5s，按照公式（49）数值模拟了激光干涉测距系统星间距离的色噪声10^-8m。

基于6点牛顿插值公式（43），星间速度色噪声表示如下

${\mathrm{\&epsiv;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}12}\left(t_i\right)=-\frac{1}{60\mathrm{\&Delta;t}}\&CenterDot;[{\mathrm{\&epsiv;}}_{\mathrm{\&rho;}12}\left(t_{i-3}\right)-9{\mathrm{\&epsiv;}}_{\mathrm{\&rho;}12}\left(t_{i-2}\right)+45{\mathrm{\&epsiv;}}_{\mathrm{\&rho;}12}\left(t_{i-1}\right)$ $\left(50\right)$

$-45{\mathrm{\&epsiv;}}_{\mathrm{\&rho;}12}\left(t_{i+1}\right)+9{\mathrm{\&epsiv;}}_{\mathrm{\&rho;}12}\left(t_{i+2}\right)-{\mathrm{\&epsiv;}}_{\mathrm{\&rho;}12}\left(t_{i+3}\right)]$

其中，ε_ρ12和

分别表示星间距离和星间速度的色噪声。

图2(b)表示基于不同的相关系数0≤μ≤0.99和相同的采样间隔5s计算的星间速度色噪声

研究结果表明：随着相关系数逐渐增加0~0.99，星间速度色噪声逐渐减小2.338×10^-9~0.263×10^-9m/s。因此，相关系数的适当增加有利于提高星间速度的精度。由于星间速度的测量精度约为10^-9m/s，因此本发明基于相关系数0.80~0.90和采样间隔5s计算了星间速度的色噪声。图3表示基于美国喷气推进实验室公布的GRACE-Level-1B实测数据以及利用相关系数0.85计算的星间距离和星间速度色噪声。

3.2）轨道位置和轨道速度的色噪声

图4表示基于不同的相关系数0≤μ≤0.99和相同的采样间隔5s，分别利用公式（49）和（50）数值模拟的视线方向轨道位置和轨道速度的色噪声。结果表明：随着相关系数的逐渐增加，轨道位置和轨道速度的色噪声逐渐减小。基于GRACE Follow-On卫星的轨道位置精度10^-5m和轨道速度精度10^-7m/s，利用相关系数0.90~0.99数值模拟了轨道位置和轨道速度的色噪声。图5表示基于相关系数0.95和采样间隔5s数值模拟的视线方向轨道位置和轨道速度的色噪声。另外，由于重力卫星轨道三轴分量几乎为等精度测量，因此卫星轨道垂向和径向分量的色噪声特性类似于视线分量。

3.3）非保守力的色噪声

基于美国喷气推进实验室公布的GRACE-Level-1B实测数据，星载加速度计非保守力的相关系数约为0.85~0.95。图6表示基于相关系数0.90和采样间隔5s，利用公式（49）模拟的加速度计视线方向的非保守力色噪声10^-13m/s²。

步骤四：地球重力场反演

基于残余星间速度观测方程（47），利用激光干涉测距仪的星间速度、GPS接收机的轨道位置和轨道速度、星载加速度计的非保守力、以及双星的参考非保守力、保守力和参考保守力等经过色噪声处理的数据，求解地球引力位系数

和

最终通过地球引力位系数的集合建立全球重力场模型。

如图7所示，圆圈线、十字线、实线和虚线分别表示基于观测时间30天、采样间隔5s和相关系数（星间速度0.85、轨道位置和轨道速度0.95、以及非保守力0.90），利用2点、4点、6点和8点残余星间速度公式，累计大地水准面精度。当相关系数和采样间隔一定，插值点数的适当增加有利于地球重力场精度的提高。研究结果表明：第一，在120阶内，基于2点残余星间速度公式的累计大地水准面精度远低于基于4点、6点和8点残余星间速度公式的精度；第二，基于6点残余星间速度公式的累计大地水准面精度远高于基于2点、4点和8点残余星间速度公式的精度。总而言之，6点残余星间速度公式是建立高精度和髙阶次全球重力场模型的优选方法之一。

图8(a)表示基于卫星观测时间30天和采样间隔5s，通过6点残余星间速度公式，利用相关系数（轨道位置和轨道速度0.95、非保守力0.90、星间速度0.80、0.85和0.90）反演地球引力位系数的精度；图8(b)表示利用相关系数（星间速度0.85、非保守力0.90、轨道位置和轨道速度0.90、0.95和0.99）反演地球引力位系数的精度；图8(c)表示利用相关系数（星间速度0.85、轨道位置和轨道速度0.95、非保守力0.85、0.90和0.95）反演地球引力位系数的精度。研究结果表明：相关系数对地球重力场精度的影响在不同频段具有不同的特性。第一，在地球重力场长波段，随着相关系数减小，地球重力场精度逐渐提高。基于高相关性的卫星观测值，地球重力场长波段信号强度将被降低，因此在一定程度上将损失地球重力场低频信号精度。第二，在地球重力场中长波段，随着相关系数的增加，由于卫星观测值的误差逐步减小，因此地球重力场精度逐渐提高。总而言之，相关系数的优化设计是进一步提高地球重力场精度的关键技术。

如图9所示，实线、虚线、十字线和圆圈线分别表示基于关键载荷误差（激光干涉测距系统的星间速度1×10^-9m/s、GPS接收机的轨道位置3×10^-5m和轨道速度1×10^-7m/s、加速度计的非保守力3×10^-13m/s²），通过6点残余星间速度公式，利用卫星观测时间30天、采样间隔5s和相关系数（星间速度0.85、轨道位置和轨道速度0.95、非保守力0.90）反演地球引力位系数的精度。以基于星间速度误差1×10^-9m/s反演地球重力场精度为标准，本发明分别基于轨道位置误差10^-6m~10^-4m、轨道速度误差10^-8m/s~10^-6m/s和非保守力误差10^-14m/s²~10^-12m/s²反演了地球重力场。研究结果表明：星间速度误差1×10^-9m/s对地球重力场精度的影响与轨道位置误差(2~3)×10^-5m、轨道速度误差(0.5~1)×10^-7m/s和非保守力误差(2~3)×10^-13m/s²相匹配。因此，残余星间速度法是获得重力卫星系统的关键载荷匹配精度指标的有效方法。

如图10所示，十字线表示德国波兹坦地学研究中心（GFZ）公布的120阶GRACE-only地球重力场模型EIGEN-GRACE02S的实测精度，在120阶处累计大地水准面精度为1.893×10^-1m。虚线和实线表示基于6点残余星间速度公式，通过卫星观测时间30天、采样间隔5s和相关系数（星间速度0.85、轨道位置和轨道速度0.95、非保守力0.90），基于关键载荷匹配精度指标（激光干涉测距系统的星间速度1×10^-9m/s、GPS接收机的轨道位置3×10^-5m和轨道速度1×10^-7m/s、加速度计的非保守力3×10^-13m/s²），分别利用星间速度和残余星间速度观测值，反演GRACE Follow-On地球重力场精度；在120阶处，GRACE Follow-On累计大地水准面精度分别为1.693×10^-4m和1.475×10^-4m。研究结果表明：（1）基于GRACE Follow-On卫星（虚线）反演地球重力场精度较当前GRACE卫星（十字线）至少提高一个数量级；（2）重力双星的残余星间速度观测量（实线）较星间速度（虚线）对地球重力场反演精度更敏感；（3）残余星间速度法是有效抑制地球重力场的高频噪声，进而精确和快速建立髙阶次地球重力场模型的有效方法。

Claims (1)

1.一种基于残余星间速度原理反演地球重力场的方法，其特征如下：

步骤一：卫星观测数据采集

（1）通过星载激光干涉测距仪获取星间速度

通过星载GPS接收机获取双星轨道位置（r₁,r₂）和双星轨道速度

通过星载加速度计获取作用于双星的非保守力（f₁，f₂）；

（2）利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置

和双星参考轨道速度

（3）参考星间速度

通过参考轨道速度

计算获得

其中，

表示相对参考轨道速度矢量，

表示第一颗参考卫星指向第二颗参考卫星的参考单位矢量，

表示相对参考轨道位置矢量，

和

表示双星的参考轨道位置矢量；

（4）参考非保守力

通过DTM2000阻力温度模型计算获得；

（5）通过国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于双星的保守力（F₁,F₂）和参考保守力

步骤二：残余星间速度观测方程建立

在地心惯性系中，基于牛顿插值原理，单星轨道速度

的泰勒展开表示如下

$\stackrel{\&CenterDot;}{r}\left(t\right)=\stackrel{\&CenterDot;}{r}\left(t_0\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)\underset{\mathrm{\&alpha;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&alpha;}}\left(\begin{array}{c}i\\ \mathrm{\&alpha;}\end{array}\right)\stackrel{\&CenterDot;}{r}\left(t_{\mathrm{\&alpha;}}\right)---\left(1\right)$

其中， $\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)$ 表示二项式系数，t表示计算点的时刻，t₀表示插值点的初始时刻，Δt表示采样间隔，n表示插值点的数量；

单星参考轨道速度

的泰勒展开表示如下

基于公式（1）-（2），单星残余轨道速度

的泰勒展开表示如下

$\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}\left(t\right)=\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}\left(t_0\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)\underset{\mathrm{\&alpha;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&alpha;}}\left(\begin{array}{c}i\\ \mathrm{\&alpha;}\end{array}\right)\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}\left(t_{\mathrm{\&alpha;}}\right)---\left(3\right)$

其中，

基于公式（3）的一阶时间导数，单星残余轨道加速度

的泰勒展开表示如下

$\mathrm{\&delta;}\stackrel{\&CenterDot;\&CenterDot;}{r}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&alpha;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&alpha;}}\left(\begin{array}{c}i\\ \mathrm{\&alpha;}\end{array}\right)\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}\left(t_{\mathrm{\&alpha;}}\right)---\left(4\right)$

其中，

基于公式（4），双星残余轨道加速度差分

的泰勒展开表示如下

$\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&alpha;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&alpha;}}\left(\begin{array}{c}i\\ \mathrm{\&alpha;}\end{array}\right)\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\left(t_{\mathrm{\&alpha;}}\right)---\left(5\right)$

其中，

和

分别表示轨道速度差分矢量和轨道加速度差分矢量，

和表示双星的轨道速度矢量，和

表示双星的参考轨道速度矢量，

和

表示双星的轨道加速度矢量，和

表示双星的参考轨道加速度矢量；

双星残余轨道加速度差分

的视线分量表示如下

$e_{12}\left(t\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&alpha;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&alpha;}}\left(\begin{array}{c}i\\ \mathrm{\&alpha;}\end{array}\right)e_{12}\left(t\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\left(t_{\mathrm{\&alpha;}}\right)---\left(6\right)$

其中，e₁₂=r₁₂/|r₁₂|表示第一颗卫星指向第二颗卫星的单位矢量，r₁₂=r₂-r₁表示双星的轨道位置差分矢量，r₁和r₂分别表示双星的轨道位置矢量；

引入激光干涉测距系统的高精度残余星间速度

进一步提高地球重力场精度，其中

和

分别表示星间速度和参考星间速度；

可被改写为

$\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}={\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}}_{12}^{\left|\right|}+{\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}}_{12}^{\&perp;}---\left(7\right)$

其中， ${\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}}_{12}^{\left|\right|}=(\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\&CenterDot;e_{12})e_{12}$ 和 ${\mathrm{\&delta;}\stackrel{\&CenterDot;}{r}}_{12}^{\&perp;}=\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}-(\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\&CenterDot;e_{12})e_{12}$ 分别表示

的视线分量和垂向分量；

基于误差传播原理，为了有效降低

将公式（7）中的替换为

因此，公式（6）可改写为

$e_{12}\left(t\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&alpha;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&alpha;}}\left(\begin{array}{c}i\\ \mathrm{\&alpha;}\end{array}\right)e_{12}\left(t\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&alpha;}}\right)---\left(8\right)$

其中， $\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&alpha;}}\right)=\mathrm{\&delta;}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_{12}\left(t_{\mathrm{\&alpha;}}\right)e_{12}\left(t_{\mathrm{\&alpha;}}\right)+\{\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\left(t_{\mathrm{\&alpha;}}\right)-[\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\left(t_{\mathrm{\&alpha;}}\right)\&CenterDot;e_{12}\left(t_{\mathrm{\&alpha;}}\right)\left]e_{12}\left(t_{\mathrm{\&alpha;}}\right)\right\};$

在公式（8）中，取插值点数n=2,4,6,8，得到2点、4点、6点和8点残余星间速度公式

$e_{12}\left(t_i\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}\left(t_i\right)=-\frac{e_{12}\left(t_i\right)}{2\mathrm{\&Delta;t}}\&CenterDot;[\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-1}\right)-\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+1}\right)]---\left(9\right)$

$e_{12}\left(t_i\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}\left(t_i\right)=\frac{e_{12}\left(t_i\right)}{12\mathrm{\&Delta;t}}\&CenterDot;[\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-2}\right)-8\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-1}\right)+8\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+1}\right)-\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+2}\right)]---\left(10\right)$

$e_{12}\left(t_i\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}\left(t_i\right)=-\frac{e_{12}\left(t_i\right)}{60\mathrm{\&Delta;t}}\&CenterDot;[\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-3}\right)-9\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-2}\right)+45\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-1}\right)---\left(11\right)$

$-45\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+1}\right)+9\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+2}\right)-\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+3}\right)]$

$e_{12}\left(t_i\right)\&CenterDot;\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}\left(t_i\right)=\frac{e_{12}\left(t_i\right)}{\mathrm{\&Delta;t}}\&CenterDot;[\frac{1}{280}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-4}\right)-\frac{4}{105}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-3}\right)+\frac{1}{5}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-2}\right)-\frac{4}{5}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i-1}\right)---\left(12\right)$

$+\frac{4}{5}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+1}\right)-\frac{1}{5}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+2}\right)+\frac{4}{105}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+3}\right)-\frac{1}{280}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_{i+4}\right)]$

在公式（8）中，

的具体形式表示如下

$\mathrm{\&delta;}{\stackrel{\&CenterDot;\&CenterDot;}{r}}_{12}=\mathrm{\&delta;}{g}_{12}+\mathrm{\&delta;}{T}_{12}+\mathrm{\&delta;}{F}_{12}+\mathrm{\&delta;}{f}_{12}---\left(13\right)$

其中，δT₁₂表示作用于双星的残余地球扰动引力差；表示除地球引力之外的残余保守力差，F₁和F₂表示作用于双星的保守力，

和

表示参考保守力；

表示残余非保守力差，f₁和f₂表示作用于双星的非保守力，

和

表示参考非保守力；

表示残余地心引力差，g₁和g₂表示双星的地心引力，

和表示参考地心引力

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

分别表示双星的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎表示轨道位置矢量r₁₍₂₎的3个分量；

通过将公式（13）和（14）代入（8），残余星间速度观测方程表示如下

$e_{12}\left(t\right)\&CenterDot;\mathrm{\&delta;}{T}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&beta;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&alpha;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&alpha;}}\left(\begin{array}{c}i\\ \mathrm{\&alpha;}\end{array}\right)e_{12}\left(t\right)\&CenterDot;\{\mathrm{\&delta;}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_{12}\left(t_{\mathrm{\&alpha;}}\right)e_{12}\left(t_{\mathrm{\&alpha;}}\right)$

$+[\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\left(t_{\mathrm{\&alpha;}}\right)-(\mathrm{\&delta;}{\stackrel{\&CenterDot;}{r}}_{12}\left(t_{\mathrm{\&alpha;}}\right)\&CenterDot;e_{12}\left(t_{\mathrm{\&alpha;}}\right))e_{12}\left(t_{\mathrm{\&alpha;}}\right)]\}---\left(15\right)$

$-e_{12}\left(t\right)\&CenterDot;[\mathrm{\&delta;}{g}_{12}+\mathrm{\&delta;}{F}_{12}\left(t\right)+\mathrm{\&delta;}{f}_{12}\left(t\right)]$

其中，

表示残余地球扰动位的一阶梯度，V₁和V₂表示地球扰动位，和

表示参考地球扰动位

$V(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}{\left(\frac{R_e}{r}\right)}^{l+1}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;}){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\cos \mathrm{\&theta;}\right)---\left(16\right)$

其中，r,θ，λ分别表示地心半径、地心余纬度和地心经度，R_e表示地球平均半径；

表示正规化的缔合Legendre函数，l表示阶数，m表示次数；

和

表示待估的地球引力位系数；

步骤三：地球重力场反演

基于残余星间速度观测方程（15），利用激光干涉测距仪的星间速度、GPS接收机的轨道位置和轨道速度、加速度计的非保守力、以及双星的参考非保守力、保守力和参考保守力，求解地球引力位系数

和

最终通过地球引力位系数的集合建立全球重力场模型。

Images