CN102393535B - 基于双星能量插值原理的卫星重力反演方法 - Google Patents

Abstract

本发明涉及卫星大地测量学、地球物理学、空间科学等交叉技术领域，提出一种基于双星能量插值原理的卫星重力反演方法。通过将星载K波段测量仪的高精度星间距离插值引入双星相对轨道动能项，进而建立双星能量插值观测方程，旨在精确和快速反演全球重力场。该方法卫星重力反演精度较高，观测方程物理含义明确，利于重力卫星系统误差分析，易于感测中高频重力场信号，计算机性能要求较低。由于可实质性提高120阶GRACE地球重力场的感测精度，因此本发明提出的双星能量插值法是建立高精度和高空间分辨率的新一代全球重力场模型的有效方法。

Description

基于双星能量插值原理的卫星重力反演方法

一、技术领域

本发明涉及卫星大地测量学、地球物理学、空间科学等交叉技术领域，特别是涉及一种为克服目前全球卫星导航系统(GNSS)的动态定轨精度相对较低的缺点，通过将星载K波段测量仪的高精度星间距离观测值引入双星相对轨道动能项，基于插值原理建立新型双星能量插值观测方程，进而精确和快速反演全球重力场的方法。

二、背景技术

卫星重力测量技术的实现是继美国GPS星座成功构建之后在大地测量领域的又一项创新和突破，它既不同于传统的车载、船载和机载测量，也不同于卫星测高和轨道摄动分析，而是通过卫星跟踪卫星技术(SST)和卫星重力梯度技术(SGG)反演高精度和高空间分辨率的地球重力场。如图1所示，重力卫星CHAMP、GRACE和GOCE的成功发射昭示着人类将迎来一个前所未有的卫星重力探测时代。确定地球重力场的精细结构及其时变不仅是大地测量、地球物理、海洋勘探应用等的需求，同时也为全人类寻求资源、保护环境和预测灾害提供重要的信息资源。第一，重力变化(重力异常)决定于地质构造和矿产资源(如煤炭、石油、天然气等)分布的不均匀性。由于我国远景资源储备中有78％的石油和93％的天然气亟待开发，因此精密重力测量有利于我国将来的资源勘探。第二，由于地震、火山、海啸、龙卷风、泥石流等自然灾害发生前后，地球重力场将发生明显的变化和异常，因此精密重力测量能帮助我们预测自然灾害的发生，进而降低人员伤亡和经济损失。

目前基于卫星重力观测数据确定地球重力场的众多卫星重力反演方法包括：(1)Kaula线性摄动法和轨道数值微分法。国内外研究表明，Kaula线性摄动法和基于加速度观测值的数值微分法只适合于求解低阶地球重力场且计算精度较低，现在较为盛行的是轨道动力学法和能量守恒法。(2)轨道动力学法。优点是求解精度较高；缺点是观测数据运算量较大、求解过程复杂程度较高且反演较高阶重力场时需要高性能的并行计算机支持。(3)能量守恒法。优点是观测方程物理含义明确且易于地球重力场的敏感度分析，在保证求解精度的前提下计算量大大降低，通常采用PC计算机可完成高阶地球重力场的快速求解；缺点是对卫星速度的测量精度要求较高。

本发明不同于前人已构建的地球重力场反演方法，首次通过在相对轨道动能中引入GRACE卫星(如图2所示)K波段测量仪的高精度星间距离，基于插值原理精确和快速建立了新型双星能量插值观测方程；并基于美国宇航局喷气推进实验室(NASA-JPL)公布的1年的GRACE Level-1B卫星观测数据，首次建立了全球重力场模型WHIGG-GEGM03S(GRACE Earth’s Gravity Model fromWuhan Institute of Geodesy and Geophysics)，进而检验了新型双星能量插值卫星重力反演法的可靠性和有效性。由于我国自主研制和正在建设的首期地球重力卫星系统预计于国家“十二五”规划末期发射升空，因此双星能量插值法以其独特的优越性(如反演精度高、计算速度快等)将成为我国高精度和高阶次全球重力场模型建立的优选方法之一。

三、发明内容

本发明的目的是：通过将星载K波段测量仪的高精度星间距离插值引入双星相对轨道动能项，建立新型双星能量插值观测方程，旨在精确和快速反演全球重力场。

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

一种基于双星能量插值原理的卫星重力反演方法，包含下列步骤：

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

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

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

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

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

步骤二：基于双星能量插值原理反演全球重力场

2.1)在地心惯性系中，建立卫星运动方程如下，

$\stackrel{\·\·}{r}(r,t^{\′})=F_e(r,t^{\′})+F_T(r,t^{\′})+f(r,t^{\′}),---\left(1\right)$

其中，

表示卫星总加速度，r表示卫星的绝对位置矢量，t′表示观测时间；F_e(r，t′)表示地球的引力；F_T(r，t′)表示三体摄动力；f(r，t′)表示作用于卫星的非保守力。在式(1)两边同时乘以轨道速度

$\stackrel{\·}{r}\·\stackrel{\·\·}{r}(r,t^{\′})=\stackrel{\·}{r}\·[F_e(r,t^{\′})+F_T(r,t^{\′})]+\stackrel{\·}{r}\·f(r,t^{\′}),---\left(2\right)$

其中，F_e(r，t′)和F_T(r，t′)可表示为

$F_{e\left(T\right)}(r,t^{\′})=\frac{{\∂V}_{e\left(T\right)}}{\∂r},---\left(3\right)$

其中，V_e＝V₀+T_e表示地球引力位，表示地球中心引力位，GM表示地球质量M和万有引力常数G之积，

表示卫星的地心半径，x，y，z表示位置矢量r的分量，T_e表示扰动位；V_T表示三体摄动能。

V_e(T)对时间的一阶导数表示如下：

$\frac{{\mathrm{dV}}_{e\left(T\right)}}{{\mathrm{dt}}^{\′}}=\frac{{\∂V}_{e\left(T\right)}}{\∂r}\·\frac{\mathrm{dr}}{{\mathrm{dt}}^{\′}}+\frac{{\∂V}_{e\left(T\right)}}{{\∂t}^{\′}}\·\frac{{\mathrm{dt}}^{\′}}{{\mathrm{dt}}^{\′}},---\left(4\right)$

将式(3)代入式(4)中可得

$F_{e\left(T\right)}(r,t^{\′})\·\stackrel{\·}{r}=\frac{{\mathrm{dV}}_{e\left(T\right)}}{{\mathrm{dt}}^{\′}}-\frac{{\∂V}_{e\left(T\right)}}{{\∂t}^{\′}},---\left(5\right)$

将式(5)代入式(2)，并两边同时积分；在地心惯性系中，单星扰动位观测方程如下：

T_e＝E_k-E_f+V_ω-V_T-V₀-E₀，        (6)

其中，

表示卫星的动能；

表示卫星的耗散能； $V_{\mathrm{\ω}}=\∫\frac{\∂(V_e+V_T)}{{\∂t}^{\′}}{\mathrm{dt}}^{\′}\≈-{\mathrm{\ω}}_e(x\stackrel{\·}{y}-y\stackrel{\·}{x})$ 表示卫星的位旋转能，ω_e表示地球自转角速度，

表示卫星轨道速度的3个分量；E₀表示卫星系统的能量积分常数。

据式(6)，在地心惯性系中，双星扰动位差的观测方程建立如下：

T_e12＝E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂，        (7)

其中，T_e12表示双星扰动位差 $T_{e12}(r_1,{\mathrm{\&theta;}}_1,{\mathrm{\&lambda;}}_1,r_2,{\mathrm{\&theta;}}_2,{\mathrm{\&lambda;}}_2)=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=-l}{\overset{l}{\mathrm{\&Sigma;}}}\left\{\right[{\left(\frac{R_e}{r_2}\right)}^{l+1}{\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}({\mathrm{\&theta;}}_2,{\mathrm{\&lambda;}}_2)-{\left(\frac{R_e}{r_1}\right)}^{l+1}{\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}({\mathrm{\&theta;}}_1,{\mathrm{\&lambda;}}_1)\left]{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\right\}$ 其中， ${\stackrel{\&OverBar;}{Y}}_{l,m}(\mathrm{\&theta;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{l,\left|m\right|}\left(\cos \mathrm{\&theta;}\right)Q_m\left(\mathrm{\&lambda;}\right),$ $Q_m\left(\mathrm{\&lambda;}\right)=\left\{\begin{array}{cc}\cos \mathrm{m\&lambda;}& m\&GreaterEqual;0\\ \sin \left|m\right|\mathrm{\&lambda;}& m<0\end{array}\right.;$ r₁₍₂₎，θ₁₍₂₎，λ₁₍₂₎表示双星各自的地心半径、地心余纬度和地心经度，R_e表示地球的平均半径；

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

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

表示双星动能差， $E_{f12}=\&Integral;({\stackrel{\&CenterDot;}{r}}_2\&CenterDot;f_2-{\stackrel{\&CenterDot;}{r}}_1\&CenterDot;f_1){\mathrm{dt}}^{\&prime;}$ 表示双星耗散能差， $V_{\mathrm{\&omega;}12}=-{\mathrm{\&omega;}}_e(x_{12}{\stackrel{\&CenterDot;}{y}}_2-y_2{\stackrel{\&CenterDot;}{x}}_{12}-y_{12}{\stackrel{\&CenterDot;}{x}}_1+x_1{\stackrel{\&CenterDot;}{y}}_{12})$ 表示双星位旋转能差，V_T12表示双星三体摄动能差，

表示双星中心引力位差，E₀₁₂表示双星能量积分常数差，可通过初始位置和速度计算得到。

2.2)将GRACE卫星的K波段测量仪的精确星间距离数据ρ₁₂引入到式(7)中的项E_k12，建立双星能量插值卫星观测方程；包括：

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

$r\left(t\right)=r\left(t_0\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)r\left(t_{\mathrm{\&tau;}}\right),---\left(8\right)$

其中，

表示二项式系数，

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

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

的展开公式：

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

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

的展开公式表示如下

${\stackrel{\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)r_{12}\left(t_{\mathrm{\&tau;}}\right),---\left(10\right)$

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

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

和

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

基于式(10)，双星相对动能的展开公式表示如下：

$\frac{1}{2}{[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;{\stackrel{\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]{\&CenterDot;r}_{12}\left(t_{\mathrm{\&tau;}}\right),---\left(11\right)$

其中， $\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{12}\left(t_{\mathrm{\&tau;}}\right)$ 可被改写为：

$\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{12}\left(t_{\mathrm{\&tau;}}\right)=\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;[r_{12}^{\left|\right|}\left(t_{\mathrm{\&tau;}}\right)+r_{12}^{\&perp;}\left(t_{\mathrm{\&tau;}}\right)],---\left(12\right)$

其中，

表示r₁₂的星星连线方向分量，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A卫星指向GRACE-B卫星的单位矢量；

表示r₁₂的垂直于星星连线方向分量。

使用GRACE卫星K波段测量仪的高精度星间距离ρ₁₂e₁₂来替换(r₁₂·e₁₂)e₁₂；式(11)可改写为：

$\frac{1}{2}{[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]{\&CenterDot;r}_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&tau;}}\right),---\left(13\right)$

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

因此，2点、4点、6点和8点星间距离插值公式表示如下

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

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

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

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

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

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

联合式(7)和式(13)，双星能量插值观测方程表示如下

T_e12(t)＝E_ρ12(t)-E_f12(t)+V_ω12(t)-V_T12(t)-V₀₁₂(t)-E₀₁₂(t)，    (18)

其中， $E_{\mathrm{\&rho;}12}{\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]{\&CenterDot;r}_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&tau;}}\right)$ 表示双星动能差；V_T12(t)＝V_E12(t)+V_S12(t)+V_M12(t)，V_E12(t)表示双星地球固体潮汐能差，V_S12(t)表示双星太阳引力位差，V_M12(t)表示双星月球引力位差，E_f12表示双星耗散能差，V_ω12表示双星位旋转能差，V₀₁₂表示双星中心引力位差，E₀₁₂表示双星能量积分常数差。

2.3)基于2点、4点、6点和8点双星能量插值观测方程(18)分别反演120阶GRACE地球重力场。

其中，在120阶内，基于2点双星能量插值观测方程反演地球重力场的精度比分别基于4点、6点和8点双星能量插值观测方程的反演精度约低2倍；在120阶内，基于6点双星能量插值观测方程反演地球重力场的精度分别高于基于2点、4点和8点双星能量插值观测方程的反演精度；基于美国宇航局喷气推进实验室NASA-JPL公布的2009年的GRACE Level-1B卫星观测数据，构建120阶高精度全球重力场模型WHIGG-GEGM03S；检验全球重力场模型WHIGG-GEGM03S的正确性和可靠性，其中：(1)基于建立的全球重力场模型WHIGG-GEGM03S计算得到全球大地水准面高N₁(经度、纬度和正高)；(2)基于国际已公布的GPS/Levelling数据计算获得全球大地水准面高N₂；(3)通过ΔN＝N₂-N₁评价全球重力场模型WHIGG-GEGM03S的正确性和可靠性。

本发明是基于新型双星能量插值卫星重力反演法有利于快速建立高精度全球重力场模型的特点而设计的。优点是：1)卫星重力反演精度较高；2)观测方程物理含义明确；3)利于重力卫星系统误差分析；4)易于感测中高频重力场信号；5)计算机性能要求较低。

四、附图说明

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

图2美国宇航局(NASA)和德国航天局(DLR)合作发射的GRACE重力双星。

图3双星能量插值观测方程中的动能、耗散能、位旋转能、地球固体潮汐能、太阳引力位、月球引力位和地球中心引力位。

图4分别基于2点、4点、6点和8点星间距离插值公式反演地球重力场的精度，其中横坐标表示地球引力位按球函数展开的阶数，纵坐标表示地球引力位系数的精度。

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

图6本发明新建立的和国际已公布的全球重力场模型精度对比，其中横坐标表示地球引力位按球函数展开的阶数，纵坐标表示地球引力位系数的精度。

五、具体实施方式

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

基于新型双星能量插值原理的卫星重力反演方法包含下列步骤：

步骤一：卫星观测数据的预处理

1.1 采集星载K波段测量仪得到的星间距离数据ρ₁₂

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

(2)基于9阶Lagrange多项式，插值获得间断的星间距离数据。

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

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

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

(3)基于3σ准则(莱以特准则)，有效剔除轨道数据中存在的粗大误差。1.3采集星载加速度计得到的卫星非保守力数据f

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

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

步骤二：基于新型双星能量插值原理精确反演全球重力场

2.1)由于当前全球卫星导航系统(GNSS)的定轨精度相对较低，因此通过将GRACE卫星K波段测量仪的高精度星间距离插值引入双星相对轨道动能项，建立新型和精确的双星能量插值观测方程。

在地心惯性系(ECI)中，卫星运动方程建立如下，

$\stackrel{\&CenterDot;\&CenterDot;}{r}(r,t^{\&prime;})=F_e(r,t^{\&prime;})+F_T(r,t^{\&prime;})+f(r,t^{\&prime;}),---\left(1\right)$

其中，表示卫星总加速度，r表示卫星的绝对位置矢量，t′表示观测时间；F_e(r，t′)表示地球的引力；F_T(r，t′)表示三体摄动力，包括太阳引力、月球引力、地球固体潮汐力、以及广义相对论效应；f(r，t′)表示作用于卫星的非保守力，包括大气阻力、太阳光压、地球辐射压、以及轨道控制力。在式(1)两边同时乘以轨道速度

$\stackrel{\&CenterDot;}{r}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;}{r}(r,t^{\&prime;})=\stackrel{\&CenterDot;}{r}\&CenterDot;[F_e(r,t^{\&prime;})+F_T(r,t^{\&prime;})]+\stackrel{\&CenterDot;}{r}\&CenterDot;f(r,t^{\&prime;}),---\left(2\right)$

其中，F_e(r，t′)和F_T(r，t′)可表示为

$F_{e\left(T\right)}(r,t^{\&prime;})=\frac{{\&PartialD;V}_{e\left(T\right)}}{\&PartialD;r},---\left(3\right)$

其中，V_e＝V₀+T_e表示地球引力位，

表示地球中心引力位，GM表示地球质量M和万有引力常数G之积，表示卫星的地心半径，x，y，z表示位置矢量r的分量，T_e表示扰动位；V_T表示三体摄动能。

V_e(T)对时间的一阶导数表示如下：

$\frac{{\mathrm{dV}}_{e\left(T\right)}}{{\mathrm{dt}}^{\&prime;}}=\frac{{\&PartialD;V}_{e\left(T\right)}}{\&PartialD;r}\&CenterDot;\frac{\mathrm{dr}}{{\mathrm{dt}}^{\&prime;}}+\frac{{\&PartialD;V}_{e\left(T\right)}}{{\&PartialD;t}^{\&prime;}}\&CenterDot;\frac{{\mathrm{dt}}^{\&prime;}}{{\mathrm{dt}}^{\&prime;}},---\left(4\right)$

将式(3)代入式(4)中可得

$F_{e\left(T\right)}(r,t^{\&prime;})\&CenterDot;\stackrel{\&CenterDot;}{r}=\frac{{\mathrm{dV}}_{e\left(T\right)}}{{\mathrm{dt}}^{\&prime;}}-\frac{{\&PartialD;V}_{e\left(T\right)}}{{\&PartialD;t}^{\&prime;}},---\left(5\right)$

将式(5)代入式(2)，并两边同时积分。在地心惯性系中，单星扰动位观测方程如下：

T_e＝E_k-E_f+V_ω-V_T-V₀-E₀，        (6)

其中，

表示卫星的动能；

表示卫星的耗散能； $V_{\mathrm{\&omega;}}=\&Integral;\frac{\&PartialD;(V_e+V_T)}{{\&PartialD;t}^{\&prime;}}{\mathrm{dt}}^{\&prime;}\&ap;-{\mathrm{\&omega;}}_e(x\stackrel{\&CenterDot;}{y}-y\stackrel{\&CenterDot;}{x})$ 表示卫星的位旋转能，ω_e表示地球自转角速度，

表示卫星轨道速度矢量的3个分量；E₀表示卫星系统的能量积分常数。

据式(6)，在地心惯性系中，双星扰动位差的观测方程建立如下：

T_e12＝E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂，        (7)

其中，T_e12表示双星扰动位差 $T_{e12}(r_1,{\mathrm{\&theta;}}_1,{\mathrm{\&lambda;}}_1,r_2,{\mathrm{\&theta;}}_2,{\mathrm{\&lambda;}}_2)=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=-l}{\overset{l}{\mathrm{\&Sigma;}}}\left\{\right[{\left(\frac{R_e}{r_2}\right)}^{l+1}{\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}({\mathrm{\&theta;}}_2,{\mathrm{\&lambda;}}_2)-{\left(\frac{R_e}{r_1}\right)}^{l+1}{\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}({\mathrm{\&theta;}}_1,{\mathrm{\&lambda;}}_1)\left]{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\right\}$ 其中， ${\stackrel{\&OverBar;}{Y}}_{l,m}(\mathrm{\&theta;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{l,\left|m\right|}\left(\cos \mathrm{\&theta;}\right)Q_m\left(\mathrm{\&lambda;}\right),$ $Q_m\left(\mathrm{\&lambda;}\right)=\left\{\begin{array}{cc}\cos \mathrm{m\&lambda;}& m\&GreaterEqual;0\\ \sin \left|m\right|\mathrm{\&lambda;}& m<0\end{array}\right..$ r₁₍₂₎，θ₁₍₂₎，λ₁₍₂₎表示双星各自的地心半径、地心余纬度和地心经度，R_e表示地球的平均半径；

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

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

表示双星动能差， $E_{f12}=\&Integral;({\stackrel{\&CenterDot;}{r}}_2\&CenterDot;f_2-{\stackrel{\&CenterDot;}{r}}_1\&CenterDot;f_1){\mathrm{dt}}^{\&prime;}$ 表示双星耗散能差， $V_{\mathrm{\&omega;}12}=-{\mathrm{\&omega;}}_e(x_{12}{\stackrel{\&CenterDot;}{y}}_2-y_2{\stackrel{\&CenterDot;}{x}}_{12}-y_{12}{\stackrel{\&CenterDot;}{x}}_1+x_1{\stackrel{\&CenterDot;}{y}}_{12})$ 表示双星位旋转能差，V_T12表示双星三体摄动能差，表示双星中心引力位差，E₀₁₂表示双星能量积分常数差，可通过初始位置和速度计算得到。

2.2)因为目前国际全球卫星导航系统(GNSS)定轨精度相对较低，所以在卫星观测方程(7)中直接使用卫星动能项E_k12将无法实质性提高地球重力场的感测精度。因此，在卫星观测方程中引入星载K波段测量仪的高精度星间距离是有效提高卫星重力反演精度的关键。不同于已有卫星重力反演方法，本发明将星载K波段测量仪的精确星间距离ρ₁₂引入卫星动能项E_k12，进而建立了新型双星能量插值卫星观测方程。

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

$r\left(t\right)=r\left(t_0\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)r\left(t_{\mathrm{\&tau;}}\right),---\left(8\right)$

其中，

表示二项式系数，

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

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

的展开公式

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

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

的展开公式表示如下

${\stackrel{\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)r_{12}\left(t_{\mathrm{\&tau;}}\right),---\left(10\right)$

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

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

和

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

基于式(10)，双星相对动能的展开公式表示如下

$\frac{1}{2}{[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;{\stackrel{\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]{\&CenterDot;r}_{12}\left(t_{\mathrm{\&tau;}}\right),---\left(11\right)$

其中， $\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{12}\left(t_{\mathrm{\&tau;}}\right)$ 可被改写为

$\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{12}\left(t_{\mathrm{\&tau;}}\right)=\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;[r_{12}^{\left|\right|}\left(t_{\mathrm{\&tau;}}\right)+r_{12}^{\&perp;}\left(t_{\mathrm{\&tau;}}\right)],---\left(12\right)$

其中，

表示r₁₂的星星连线方向分量，e₁₂＝r₁₂/|r₁₂|表示由GRACE-A卫星指向GRACE-B卫星的单位矢量；

表示r₁₂的垂直于星星连线方向分量。

本发明引入了星载GRACE卫星K波段测量仪的高精度星间距离ρ₁₂e₁₂来替换(r₁₂·e₁₂)e₁₂。因此，式(11)可改写为

$\frac{1}{2}{[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]{\&CenterDot;r}_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&tau;}}\right),---\left(13\right)$

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

基于式(13)，2点、4点、6点和8点星间距离插值公式表示如下

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

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

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

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

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

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

联合式(7)和式(13)，新型双星能量插值观测方程表示如下

T_e12(t)＝E_ρ12(t)-E_f12(t)+V_ω12(t)-V_T12(t)-V₀₁₂(t)-E₀₁₂(t)，    (18)

其中，如图3所示， $E_{\mathrm{\&rho;}12}{\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]{\&CenterDot;r}_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&tau;}}\right)$ 表示双星动能差；V_T12(t)＝V_E12(t)+V_S12(t)+V_M12(t)，V_E12(t)表示双星地球固体潮汐能差，V_S12(t)表示双星太阳引力位差，V_M12(t)表示双星月球引力位差；E_f12表示双星耗散能差，V_ω12表示双星位旋转能差，V₀₁₂表示双星中心引力位差，E₀₁₂表示双星能量积分常数差。

2.3)如图4所示，基于2点、4点、6点和8点双星能量插值观测方程(18)分别反演了120阶GRACE地球重力场，结果表明：基于较优的信噪比，6点双星能量插值观测方程可实质性提高卫星重力反演精度。第一，在120阶内，基于2点双星能量插值观测方程反演地球重力场的精度较分别基于4点、6点和8点双星能量插值观测方程的反演精度约低2倍。原因分析如下：首先，由于式(14)的左边

是点域值，而右边

是平均值，因此式(14)的左右两边几乎不相等；其次，由于2点双星能量插值观测方程的插值点数较少，因此无法提供足够的地球重力场反演插值信息。第二，在120阶内，基于6点双星能量插值观测方程反演地球重力场的精度高于分别基于2点、4点和8点双星能量插值观测方程的反演精度。原因分析如下：由于随着插值点数的增多，卫星观测值的信息量逐渐增加，因此基于6点双星能量插值观测方程反演地球重力场的精度高于分别基于2点和4点双星能量插值观测方程的反演精度。但是，随着插值点数的增多，卫星观测值的误差量也同时增加，因此基于8点双星能量插值观测方程反演地球重力场的精度低于基于6点双星能量插值观测方程的反演精度。综上所述，6点双星能量插值观测方程是有效提高120阶地球重力场反演精度的较优选择。

2.4)基于美国宇航局喷气推进实验室(NASA-JPL)公布的2009年的GRACELevel-1B卫星观测数据，构建了120阶新型和高精度全球重力场模型WHIGG-GEGM03S。图5中实线和虚线表示利用NASA-JPL公布的2009年的GRACE Level-1B卫星观测数据和基于新型双星能量插值法，分别采用低精度的轨道位置观测值(见式(7))和星载K波段测量仪的高精度星间距离观测值(见式(18))反演GRACE地球重力场的精度，在120阶处累计大地水准面精度分别为5.445×10¹m和1.496×10^-1m。通过2条曲线的对比可知，新型双星能量插值法是反演高精度和高空间分辨率地球重力场的优选方法之一。

2.5)结合中国、美国、德国和澳大利亚的GPS/Levelling数据以及国际已公布的GRACE全球重力场模型，检验了新建立的全球重力场模型WHIGG-GEGM03S的正确性和可靠性。具体计算过程如下：(1)基于本发明新建立的全球重力场模型WHIGG-GEGM03S计算获得全球大地水准面高N₁(经度、纬度和正高)；(2)基于国际已公布的GPS/Levelling数据计算获得全球大地水准面高N₂；(3)通过ΔN＝N₂-N₁评价全球重力场模型WHIGG-GEGM03S的正确性和可靠性。研究结果表明：相对于其它全球重力场模型EIGEN-CHAMP03S(RMS＝0.863m)、EIGEN-CG03C(RMS＝0.382m)、EIGEN-GL04S1(RMS＝0.635m)和EIGEN-5C(RMS＝0.361m)，本发明新建立的全球重力场模型WHIGG-GEGM03S(RMS＝0.769m)的大地水准面误差更接近于全球重力场模型EIGEN-GRACE02S(RMS＝0.752m)的误差。另外，如图6所示，本发明对比了全球重力场模型EIGEN-CHAMP03S、EIGEN-GRACE02S、EIGEN-CG03C、EIGEN-GL04S1、EIGEN-5C和WHIGG-GEGM03S的精度。综上所述，本发明新建立的全球重力场模型WHIGG-GEGM03S是正确的和可靠的。

Claims (4)

1.一种基于双星能量插值原理的卫星重力反演方法，包含下列步骤：

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

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

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

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

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

步骤二：基于双星能量插值原理反演全球重力场

2.1)在地心惯性系中，建立卫星运动方程如下，

$\stackrel{\&CenterDot;\&CenterDot;}{r}(r,t^{\&prime;})=F_e(r,t^{\&prime;})+F_T(r,t^{\&prime;})+f(r,t^{\&prime;}),$ (1)其中，

表示卫星总加速度，r表示卫星的绝对位置矢量，t′表示观测时间；F_e(r，t′)表示地球的引力；F_T(r，t′)表示三体摄动力；f(r，t′)表示作用于卫星的非保守力；在式(1)两边同时乘以轨道速度

，

$\stackrel{\&CenterDot;}{r}\&CenterDot;\stackrel{\&CenterDot;\&CenterDot;}{r}(r,t^{\&prime;})=\stackrel{\&CenterDot;}{r}\&CenterDot;[F_e(r,t^{\&prime;})+F_T(r,t^{\&prime;})]+\stackrel{\&CenterDot;}{r}\&CenterDot;f(r,t^{\&prime;}),$ (2)其中，F_e(r，t′)和F_T(r，t′)可表示为

$F_{e\left(T\right)}(r,t^{\&prime;})=\frac{\&PartialD;V_{e\left(T\right)}}{\&PartialD;r},$ (3)其中，V_e＝V₀+T_e表示地球引力位，

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

表示卫星的地心半径，x，y，z表示位置矢量r的分量，T_e表示扰动位；V_T表示三体摄动能；

V_e(T)对时间的一阶导数表示如下：

$\frac{d{V}_{e\left(T\right)}}{d{t}^{\&prime;}}=\frac{\&PartialD;V_{e\left(T\right)}}{\&PartialD;r}\&CenterDot;\frac{\mathrm{dr}}{d{t}^{\&prime;}}+\frac{\&PartialD;V_{e\left(T\right)}}{\&PartialD;t^{\&prime;}}\&CenterDot;\frac{d{t}^{\&prime;}}{d{t}^{\&prime;}},$ (4)将式(3)代入式(4)中可得

$F_{e\left(T\right)}(r,t^{\&prime;})\&CenterDot;\stackrel{\&CenterDot;}{r}=\frac{d{V}_{e\left(T\right)}}{d{t}^{\&prime;}}-\frac{\&PartialD;V_{e\left(T\right)}}{\&PartialD;t^{\&prime;}},$ (5)将式(5)代入式(2)，并两边同时积分；在地心惯性系中，单星扰动位观测方程如下：

T_e＝E_k-E_f+V_ω-V_T-V₀-E₀，(6)其中，

表示卫星的动能；表示卫星的耗散能；

表示卫星的位旋转能，ω_e表示地球自转角速度，

表示卫星轨道速度的3个分量；E₀表示卫星系统的能量积分常数；

据式(6)，在地心惯性系中，双星扰动位差的观测方程建立如下：

T_e12＝E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂，(7)其中，T_e12表示双星扰动位差

$T_{e12}(r_1,{\mathrm{\&theta;}}_1,{\mathrm{\&lambda;}}_1,r_2,{\mathrm{\&theta;}}_2,{\mathrm{\&lambda;}}_2)=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=-l}{\overset{l}{\mathrm{\&Sigma;}}}\left\{\right[{\left(\frac{R_e}{r_2}\right)}^{l+1}{\stackrel{-}{Y}}_{\mathrm{lm}}({\mathrm{\&theta;}}_2,{\mathrm{\&lambda;}}_2)-{\left(\frac{R_e}{r_1}\right)}^{l+1}{\stackrel{-}{Y}}_{\mathrm{lm}}({\mathrm{\&theta;}}_1,{\mathrm{\&lambda;}}_1)\left]{\stackrel{-}{C}}_{\mathrm{lm}}\right\}$ 其中， ${\stackrel{-}{Y}}_{l,m}(\mathrm{\&theta;},\mathrm{\&lambda;})={\stackrel{-}{P}}_{l,\left|m\right|}\left(\cos \mathrm{\&theta;}\right)Q_m\left(\mathrm{\&lambda;}\right),$ $Q_m\left(\mathrm{\&lambda;}\right)=\{\begin{array}{cc}\cos \mathrm{m\&lambda;}& m\&GreaterEqual;0\\ \sin \left|m\right|\mathrm{\&lambda;}& m<0\end{array};$ r₁(2)，θ₁(2)，λ₁(2)表示双星各自的地心半径、地心余纬度和地心经度，R_e表示地球的平均半径；表示规格化的Legendre函数，其中l表示阶数，m表示次数；

表示待求的规格化地球引力位系数；表示双星动能差，

表示双星耗散能差，

表示双星位旋转能差，V_T12表示双星三体摄动能差，

表示双星中心引力位差，E_O12表示双星能量积分常数差，可通过初始位置和速度计算得到；

2.2)将GRACE卫星的K波段测量仪的精确星间距离数据ρ₁₂引入到式(7)中的项E_k12，建立双星能量插值卫星观测方程；包括：

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

$r\left(t\right)=r\left(t_0\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)r\left(t_{\mathrm{\&tau;}}\right),$ (8)

其中，

表示二项式系数，

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

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

的展开公式：

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

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

的展开公式表示如下

${\stackrel{\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)r_{12}\left(t_{\mathrm{\&tau;}}\right),$ (10)其中，r₁₂＝r₂-r₁和分别表示双星相对轨道位置矢量和相对轨道速度矢量，r₁和r₂分别表示双星绝对轨道位置矢量，

和

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

基于式(10)，双星相对动能的展开公式表示如下：

$\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;{\stackrel{\&CenterDot;}{r}}_{12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{12}\left(t_{\mathrm{\&tau;}}\right),---\left(11\right)$

其中， $\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{12}\left(t_{\mathrm{\&tau;}}\right)$ 可被改写为：

$\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{12}\left(t_{\mathrm{\&tau;}}\right)=\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;[r_{12}^{\left|\right|}\left(t_{\mathrm{\&tau;}}\right)+r_{12}^{\&perp;}\left(t_{\mathrm{\&tau;}}\right)],$ (12)其中，

表示r₁₂的星星连线方向分量，

表示由GRACE-A卫星指向GRACE-B卫星的单位矢量；

表示r₁₂的垂直于星星连线方向分量；

使用GRACE卫星K波段测量仪的高精度星间距离ρ₁₂e₁₂来替换式(11)可改写为：

$\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{.}{r}}_1\left(t\right)]\&CenterDot;{\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&tau;}}\right),$ (13)其中，r_ρ12(t_τ)＝ρ₁₂(t_τ)e₁₂(t_τ)+{r₁₂(t_τ)-[r₁₂(t_τ)e₁₂(t_τ)]e₁₂(t_τ)}；

因此，2点、4点、6点和8点星间距离插值公式表示如下

${\stackrel{\&CenterDot;}{r}}_{12}\left(t_i\right)=-\frac{1}{2\mathrm{\&Delta;t}}[r_{\mathrm{\&rho;}12}\left(t_{i-1}\right)-r_{\mathrm{\&rho;}12}\left(t_{i+1}\right)],$ (14) ${\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_i\right)=\frac{1}{12\mathrm{\&Delta;t}}[r_{\mathrm{\&rho;}12}\left(t_{i-2}\right)-8r_{\mathrm{\&rho;}12}\left(t_{i-1}\right)+8r_{\mathrm{\&rho;}12}\left(t_{i+1}\right)-r_{\mathrm{\&rho;}12}\left(t_{i+2}\right)],$ (15) ${\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_i\right)=-\frac{1}{60\mathrm{\&Delta;t}}[r_{\mathrm{\&rho;}12}\left(t_{i-3}\right)-{9r}_{\mathrm{\&rho;}12}\left(t_{i-2}\right)+45r_{\mathrm{\&rho;}12}\left(t_{i-1}\right),$ (16) ${\stackrel{\&CenterDot;}{r}}_{\mathrm{\&rho;}12}\left(t_i\right)=\frac{1}{\mathrm{\&Delta;t}}[\frac{1}{280}{r}_{\mathrm{\&rho;}12}\left(t_{i-4}\right)-\frac{4}{105}{r}_{\mathrm{\&rho;}12}\left(t_{i-3}\right)+\frac{1}{5}{r}_{\mathrm{\&rho;}12}\left(t_{i-2}\right)-\frac{4}{5}{r}_{\mathrm{\&rho;}12}\left(t_{i-1}\right)$ ，(17)

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

联合式(7)和式(13)，双星能量插值观测方程表示如下

T_e12(t)＝E_ρ12(t)-E_f12(t)+V_ω12(t)-V_T12(t)-V₀₁₂(t)-E₀₁₂(t)，(18)其中， $E_{\mathrm{\&rho;}12}\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(\begin{array}{c}\mathrm{\&lambda;}\\ i\end{array}\right)}^{\&prime;}\underset{\mathrm{\&tau;}=0}{\overset{i}{\mathrm{\&Sigma;}}}{(-1)}^{i+\mathrm{\&tau;}}\left(\begin{array}{c}i\\ \mathrm{\&tau;}\end{array}\right)\frac{1}{2}[{\stackrel{\&CenterDot;}{r}}_2\left(t\right)+{\stackrel{\&CenterDot;}{r}}_1\left(t\right)]\&CenterDot;r_{\mathrm{\&rho;}12}\left(t_{\mathrm{\&tau;}}\right)$ 表示双星动能差，

$-45r_{\mathrm{\&rho;}12}\left(t_{i+1}\right)+9r_{\mathrm{\&rho;}12}\left(t_{i+2}\right)-r_{\mathrm{\&rho;}12}\left(t_{i+3}\right)]$ V_T12(t)＝V_E12(t)+V_S12(t)+V_M12(t)，V_E12(t)表示双星地球固体潮汐能差，V_S12(t)表示双星太阳引力位差，V_M12(t)表示双星月球引力位差，E_f12(t)表示双星耗散能差，V_ω12(t)表示双星位旋转能差，V₀₁₂(t)表示双星中心引力位差，E₀₁₂(t)表示双星能量积分常数差；

2.3)基于2点、4点、6点和8点双星能量插值观测方程(18)分别反演120阶GRACE地球重力场。

2.如权利要求1所述的方法，其特征在于，在120阶内，基于6点双星能量插值观测方程反演地球重力场的精度分别高于基于2点、4点和8点双星能量插值观测方程的反演精度。

3.如权利要求1所述的方法，其特征在于，基于美国宇航局喷气推进实验室NASA-JPL公布的2009年的GRACE Level-1B卫星观测数据，构建120阶高精度全球重力场模型WHIGG-GEGM03S。

4.如权利要求3所述的方法，其特征在于，检验全球重力场模型WHIGG-GEGM03S的正确性和可靠性，其中：(1)基于建立的全球重力场模型WHIGG-GEGM03S计算得到全球大地水准面高N₁；(2)基于国际已公布的GPS/Levelling数据计算获得全球大地水准面高N₂；(3)通过ΔN＝N₂-N₁评价全球重力场模型WHIGG-GEGM03S的正确性和可靠性。

Images