CN103091723B - 降低重力卫星质心调整误差对地球重力场精度影响的方法 - Google Patents

Abstract

本发明公开了一种降低重力卫星质心调整误差对地球重力场精度影响的方法；基于新型扰动能量原理，通过建立的重力卫星质心调整观测方程，有效降低重力卫星质心调整误差对地球引力位系数精度影响，进而高精度和高空间分辨率反演地球重力场；该方法卫星重力反演精度高，质心调整误差对重力场精度影响小，易于卫星重力系统需求分析，卫星观测方程物理含义明确，计算机性能要求低。因此，扰动能量法是建立高精度和高空间分辨率地球重力场模型的有效方法。

Description

降低重力卫星质心调整误差对地球重力场精度影响的方法

一、技术领域

本发明涉及卫星大地测量学、卫星重力学、航空航天等交叉技术领域，特别是涉及一种基于新型扰动能量原理，通过建立重力卫星质心调整观测方程，有效降低重力卫星质心调整误差对地球引力位系数精度影响的方法，进而高精度和高空间分辨率反演地球重力场。

二、背景技术

自1957年10月4日成功发射第一颗人造卫星Sputnik-1以来，国内外的许多学者在利用卫星技术精密探测地球重力场方面取得了辉煌成就。21世纪是人类利用SST（Satellite-to-Satellite Tracking）和SGG（Satellite Gravity Gradiometry）技术提升对“数字地球”认知能力的新纪元。地球重力场及其时变特性反映地球表层及内部物质的空间分布、运动和变化，同时决定着大地水准面的起伏和变化。因此，确定地球重力场的精细结构及其时变不仅是大地测量学、地震学、海洋学、空间科学、国防建设等的需求，同时也将为寻求资源、保护环境和预测灾害提供重要的信息资源。

GRACE由美国宇航局（NASA）和德国航天局（DLR）共同研制开发，已于2002年3月17日发射升空。如图1所示，GRACE初始轨道高度为500km，采用近圆和极地轨道设计，轨道倾角89°，轨道离心率e<0.004，在10年的飞行使命中轨道高度由500km降到300km。GRACE-A/B双星采用卫星跟踪卫星高低/低低相结合的飞行模式（SST-HL/LL），除利用高轨GPS卫星对低轨双星精密跟踪定位，同时两颗低轨卫星在同一轨道平面内前后相互跟踪（星间距离220±50km）编队飞行，并利用共轨双星轨道摄动之差高精度测量地球重力场。GRACE卫星利用冷气微推进器和磁力矩器辅助双频GPS接收机精密定轨，利用K波段高频链路高精度测量星间距离和星间速度，利用高精度SuperSTAR静电悬浮加速度计测量作用于卫星的非保守力。

在地心惯性系中研究卫星绕地球的运动规律，通常将卫星视为质点。因此，在卫星飞行中作用于卫星的非保守力可以等效为作用于卫星的质点处。在卫星重力测量中，为了将地球引力从卫星受到的合外力中有效分离，作用于卫星的非保守力的精确扣除是能否反演高精度和高空间分辨率地球重力场的关键技术，因此GRACE星载加速度计检验质量的质心要求精确定位于卫星体的质心处。由于卫星在实际飞行中，卫星体的质心和星载加速度计检验质量的质心实时存在偏移，因此质心偏差的研究是加速度计能否将作用于卫星体的非保守力精确扣除的关键技术。GRACE卫星体和星载加速度计检验质量的质心偏差源主要来自于两个方面：第一，地面安装误差源：由于在地面安装时卫星体质心和加速度计检验质量质心存在偏移，导致了GRACE卫星加速度计的静电力和作用于卫星的非保守力存在固有偏差。第二，在轨飞行误差源：由于空间环境（温度、压力等）的复杂性导致在轨飞行的卫星发生形变以及对卫星进行实时轨道和姿态控制引起喷气燃料消耗（每2～3min喷气1次，每次喷气时间200～300ms），将会导致GRACE卫星体和星载加速度计检验质量的质心存在实时偏差。由于GRACE卫星体和星载加速度计检验质量的质心偏差和卫星姿态测量具有耦合效应，因此在地球重力场反演时会同时将卫星姿态测量误差引入卫星观测方程。GRACE星体和加速度计检验质量的质心偏差以及卫星姿态测量误差的引入必将会在加速度计的三轴测量中附加扰动误差，从而影响地球重力场反演的精度。因此，GRACE星体和星载加速度计检验质量质心偏差的系统研究是提高地球重力场反演精度的重要保证。

国际卫星重力测量计划的成功实施对我国既存在机遇又不乏挑战。我国应尽快汲取国外长期积累的成功经验，积极推动我国卫星重力测量计划的实施，加快自主研制重力卫星的步伐，通过卫星重力测量计划的实现带动相关领域（地学、航天、电子、通信、材料等）的发展。目前国内大地测量学界紧跟国际卫星重力测量的动态，积极投身于重力卫星需求分析的研究当中。如果GRACE星体和星载加速度计检验质量的质心调整精度设计合理，在保证利用加速度计精确扣除非保守力的前提下，可以适当降低星载加速度计和质心调整系统研制的难度以及避免不必要的人力、物力和财力的浪费。基于此目的，本发明开展了GRACE星体和星载加速度计检验质量的不同质心偏差影响GRACE地球重力场精度的研究论证，不仅为我国卫星重力测量计划中质心调整精度的优化选取提供了理论基础和技术保证，同时为将来研究工作的深入开展搭建了桥梁和纽带。

三、发明内容

本发明的目的是：基于新型扰动能量原理，通过建立重力卫星质心调整观测方程，有效降低重力卫星质心调整误差对地球引力位系数精度影响，进而高精度和高空间分辨率反演地球重力场。

本发明采用能量守恒原理反演地球重力场，并在此基础上引入参考数据，通过测量数据和参考数据的残差建立新型扰动能量观测方程，进一步提高了地球重力场的反演精度。

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

一种降低重力卫星质心调整误差对地球重力场精度影响的方法，包含下列步骤：

第一步，通过真实重力卫星采集测量数据和利用理想重力卫星获取参考数据，

所述测量数据包括，

通过K波段测距仪获取星间速度数据通过星载GPS接收机获取双星轨道位置数据（r₁,r₂）和双星轨道速度数据通过星载加速度计获取作用于双星的非保守力数据（f₀₁,f₀₂），通过星敏感器获取角速度数据Ω=(Ω_x,Ω_y,Ω_z)；

所述参考数据包括，

通过星载GPS接收机获取的双星轨道位置数据（r₁,r₂）和双星轨道速度数据的初始点，利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置数据和双星参考轨道速度数据

参考星间速度数据通过参考轨道速度数据计算获得其中，表示相对参考轨道速度矢量，表示第一颗卫星指向第二颗卫星的参考单位矢量，表示相对参考轨道位置矢量；

参考非保守力数据通过DTM2000阻力温度模型计算获得；

分别将测量数据和参考数据代入国际公布模型DE-405、IERS96和CSR4.0，联合计算获取作用于双星的三体摄动能差数据V_T12和参考三体摄动能差数据

第二步，计算各个重力卫星由于单个卫星体和星载加速度计检验质量的质心偏移而在加速度计三轴附加的非保守力Δf

$\mathrm{\&Delta;f}=\stackrel{\&CenterDot;\&CenterDot;}{L}+\mathrm{\&omega;}\&times;(\mathrm{\&omega;}\&times;L)+2\mathrm{\&omega;}\&times;\stackrel{\&CenterDot;}{L}+\stackrel{\&CenterDot;}{\mathrm{\&omega;}}\&times;L-K\&CenterDot;L---\left(1\right)$

其中，L表示卫星体质心和星载加速度计检验质量质心偏移在星体坐标系O_S-X_SY_SZ_S中的位置矢量，和分别表示L对时间的一阶导数和二阶导数；ω=R^TΩR表示在星体坐标系中重力卫星绕星体坐标系（X_S,Y_S,Z_S）轴旋转的角速度矩阵，Ω表示在地心惯性系中重力卫星绕地心坐标系（X_I,Y_I,Z_I）轴旋转的角速度矩阵

$\mathrm{\&omega;}=\left(\begin{array}{ccc}0& -{\mathrm{\&omega;}}_z& {\mathrm{\&omega;}}_y\\ {\mathrm{\&omega;}}_z& 0& -{\mathrm{\&omega;}}_x\\ -{\mathrm{\&omega;}}_y& {\mathrm{\&omega;}}_x& 0\end{array}\right)$

R表示由星体坐标系到地心惯性系的转换矩阵

$R=\left(\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right)$

$R_{11}=\frac{x_I}{\left|r\right|},$ $R_{12}=\frac{y_I}{\left|r\right|},$ $R_{13}=\frac{z_I}{\left|r\right|}$

$R_{21}=\frac{(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)z_I-(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)y_I}{\left|r\right|\left|n\right|}$

$R_{22}=\frac{(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)x_I-(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)z_I}{\left|r\right|\left|n\right|}$

$R_{23}=\frac{(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)y_I-(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)x_I}{\left|r\right|\left|n\right|}$

$R_{31}=\frac{y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I}{\left|n\right|},$ $R_{32}=\frac{z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I}{\left|n\right|},$ $R_{33}=\frac{x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I}{\left|n\right|}$

$\left|r\right|=\sqrt{x_I^2+y_I^2+z_I^2}$

$\left|n\right|=\sqrt{{(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)}^2-{(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)}^2+{(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)}^2}$

其中，x_I,y_I,z_I分别表示重力卫星在地心惯性坐标系中位置矢量的3个分量，分别表示重力卫星在地心惯性坐标系中速度矢量的3个分量；K表示在星体坐标系中地球引力位的梯度

$K=\left(\begin{array}{ccc}\frac{{\&PartialD;}^2{V}_S}{\&PartialD;x_S^2}& \frac{{\&PartialD;}^2{V}_s}{{\&PartialD;x}_S\&PartialD;y_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_S{\&PartialD;z}_S}\\ \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S{\&PartialD;x}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S^2}& \frac{{\&PartialD;}^2{V}_s}{{\&PartialD;y}_S{\&PartialD;z}_S}\\ \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;x}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;y}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S^2}\end{array}\right)$

其中，V_S=GMr表示中心引力位，GM表示地球质量M和万有引力常数G之积，x_S,y_S,z_S分别表示重力卫星在星体坐标系中位置矢量的3个分量

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_S^2}=\frac{\mathrm{GM}({6x}_S-r)}{r^4},$ $\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S^2}=\frac{\mathrm{GM}({6y}_S-r)}{r^4},$ $\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S^2}=\frac{\mathrm{GM}({6z}_S-r)}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_s{\&PartialD;y}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S{\&PartialD;x}_S}=\frac{3\mathrm{GM}{x}_S{y}_S}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_s{\&PartialD;z}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;x}_S}=\frac{3\mathrm{GM}{x}_S{z}_S}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_s{\&PartialD;z}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;y}_S}=\frac{3\mathrm{GM}{y}_S{z}_S}{r^4}.$

第三步，反演地球重力场

计算作用于卫星体的非保守力矢量和为和将位于星体系的和转化到地心系计算作用于双星的非保守力和将位于星体系的和转化到地心系计算作用于双星的参考非保守力 ${f_1}^o={\mathrm{Rf}}_{01}^o$ 和 ${f_2}^o={\mathrm{Rf}}_{02}^o;$

在地心惯性系中，双星扰动位差观测方程建立如下

T_e12=E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂         （2）

在地心惯性系中，双星参考扰动位差观测方程建立如下

$T_{e12}^o=E_{k12}^o-E_{f12}^o+V_{\mathrm{\&omega;}12}^o-V_{T12}^o-V_{012}^o-E_{012}^o---\left(3\right)$

其中，带上角标“o”的参数表示参考值；

由公式（2）－公式（3）得到双星相对扰动位差观测方程

δT_e12=δE_k12-δE_f12+δV_ω12-δV_T12-δV₀₁₂-δE₀₁₂       （4）其中，表示双星相对扰动位差，T_e12(r,θ,λ)表示地球扰动位

$T_{e12}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L_0}{\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}}_{\mathrm{lm}}(\mathrm{\&theta;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\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.;$ 表示双星各自的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎分别表示双星各自位置矢量r₁₍₂₎的三个分量，θ₁和θ₂分别表示双星各自的地心余纬度，λ₁和λ₂分别表示双星各自的地心经度；R_e表示地球的平均半径，L₀表示地球扰动位按球函数展开的最大阶数；表示规格化的Legendre函数，l表示阶数，m表示次数；表示待求的规格化地球引力位系数；利用参考数据同样计算参考地球扰动位

双星相对扰动位差观测方程（4）右边第一项表示双星的扰动动能差， $E_{k12}=\frac{1}{2}({\stackrel{\&CenterDot;}{r}}_2+{\stackrel{\&CenterDot;}{r}}_1)\&CenterDot;\{{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_{12}{e}_{12}+[{\stackrel{\&CenterDot;}{r}}_{12}-({\stackrel{\&CenterDot;}{r}}_{12}\&CenterDot;e_{12})e_{12}\left]\right\}$ 表示双星的动能差，表示双星的相对速度矢量，e₁₂=r₁₂/|r₁₂|表示由第一颗卫星指向第二颗卫星的单位方向矢量，r₁₂=r₂-r₁表示双星的相对位置矢量；利用参考数据同样计算双星的参考动能差

第二项表示双星的扰动耗散能差，表示双星的耗散能差；利用参考数据同样计算双星的参考耗散能差

第三项表示双星的扰动旋转能差，表示双星的旋转能差，ω_e表示地球的自转角速度，x₁₂=x₂-x₁和y₁₂=y₂-y₁分别表示双星的相对轨道位置，和分别表示双星的相对轨道速度；利用参考数据同样计算双星的参考旋转能差

第四项表示双星的扰动三体摄动能差，V_T12表示双星的三体摄动能差，表示双星的参考三体摄动能差；

第五项表示双星的扰动中心引力位差，V₀₁₂=GM/r₂-GM/r₁表示双星的中心引力位差；利用参考轨道数据同样计算双星的参考中心引力位差 $V_{012}^o=\mathrm{GM}/{r_2}^o-\mathrm{GM}/{r_1}^o;$

最后一项表示双星的扰动能量积分常数差，E₀₁₂表示双星的能量积分常数差，通过卫星的初始位置和初始速度计算得到；表示双星的参考能量积分常数差，通过卫星的参考初始位置和参考初始速度计算得到；

基于最小二乘法解算双星相对扰动位差观测方程（4），进而获得地球引力位系数

本发明是基于扰动能量法有利于快速反演高精度和高空间分辨率地球重力场的特点而设计的，优点是：

1）卫星重力反演精度高；

2）质心调整误差对重力场精度影响小；

3）易于卫星重力系统需求分析；

4）卫星观测方程物理含义明确；

5）计算机性能要求低。

四、附图说明

图1表示GRACE重力卫星。

图2表示GRACE星体和加速度计检验质量的质心偏移。

图3表示基于扰动能量法利用不同质心调整精度反演地球引力位系数精度。

图4表示基于未引入参考数据和引入参考数据反演地球重力场精度的对比。

五、具体实施方式

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

降低重力卫星质心调整误差对地球重力场精度影响的方法包含下列步骤：

步骤一：重力卫星测量数据采集

（1）通过K波段测距仪获取星间速度数据通过星载GPS接收机获取双星轨道位置数据（r₁,r₂）和双星轨道速度数据通过星载加速度计获取作用于双星的非保守力数据（f₁,f₂），通过星敏感器获取角速度数据Ω=(Ω_x,Ω_y,Ω_z)。

（2）通过星载GPS接收机获取的双星轨道位置数据（r₁,r₂）和双星轨道速度数据的初始点，利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置数据和双星参考轨道速度数据

（3）参考星间速度数据通过参考轨道速度数据计算获得其中，表示相对参考轨道速度矢量，表示第一颗卫星指向第二颗卫星的参考单位矢量，表示相对参考轨道位置矢量。

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

（5）分别将测量数据和参考数据代入国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于双星的三体摄动能差数据V_T12和参考三体摄动能差数据 $V_{T12}^o.$

步骤二：重力卫星质心调整观测方程建立

如图2所示，在地心惯性系O_I-X_IY_IZ_I中，GRACE卫星体质心和星载加速度计检验质量质心的位置矢量关系如下

r_A=r_S+RL             （5）其中，r_A和r_S分别表示星载加速度计质心和卫星体质心在地心惯性系中的位置矢量，L表示卫星体质心和星载加速度计检验质量质心偏移在星体坐标系O_S-X_SY_SZ_S中的位置矢量，和分别表示L对时间的一阶导数和二阶导数，R表示由星体坐标系到地心惯性系的转换矩阵

$R_{11}=\frac{x_I}{\left|r\right|},$ $R_{12}=\frac{y_I}{\left|r\right|},$ $R_{13}=\frac{z_I}{\left|r\right|}$

$R_{21}=\frac{(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)z_I-(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)y_I}{\left|r\right|\left|n\right|}$

$R_{22}=\frac{(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)x_I-(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)z_I}{\left|r\right|\left|n\right|}$

$R_{23}=\frac{(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)y_I-(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)x_I}{\left|r\right|\left|n\right|}$

$R_{31}=\frac{y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I}{\left|n\right|},$ $R_{32}=\frac{z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I}{\left|n\right|},$ $R_{33}=\frac{x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I}{\left|n\right|}$

$\left|r\right|=\sqrt{x_I^2+y_I^2+z_I^2}$

$\left|n\right|=\sqrt{{(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)}^2-{(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)}^2+{(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)}^2}$

其中，x_I,y_I,z_I分别表示重力卫星在地心惯性坐标系中位置矢量的3个分量，分别表示重力卫星在地心惯性坐标系中速度矢量的3个分量；

在公式（5）两边同时对时间t求导，可得速度运动方程

${\stackrel{\&CenterDot;}{r}}_A={\stackrel{\&CenterDot;}{r}}_S+\mathrm{\&Omega;RL}+R\stackrel{\&CenterDot;}{L}---\left(6\right)$

其中，Ω=(Ω_x,Ω_y,Ω_z)表示在地心惯性系中GRACE绕（X_I,Y_I,Z_I）轴旋转的角速度矩阵。

在公式（6）两边同时对时间t求导，可得加速度运动方程

${\stackrel{\&CenterDot;\&CenterDot;}{r}}_A={\stackrel{\&CenterDot;\&CenterDot;}{r}}_S+R\stackrel{\&CenterDot;\&CenterDot;}{L}+{\mathrm{\&Omega;}}^2\mathrm{RL}+2\mathrm{\&Omega;R}\stackrel{\&CenterDot;}{L}+\stackrel{\&CenterDot;}{\mathrm{\&Omega;}}\mathrm{RL}---\left(7\right)$

在公式（7）两边同时左乘R^T，可得星体坐标系中的加速度运动方程

${\stackrel{\&CenterDot;\&CenterDot;}{r}}_A^0={\stackrel{\&CenterDot;\&CenterDot;}{r}}_S^0+\stackrel{\&CenterDot;\&CenterDot;}{L}+{\mathrm{\&omega;}}^{2}L+2\mathrm{\&omega;}\stackrel{\&CenterDot;}{L}+\stackrel{\&CenterDot;}{\mathrm{\&omega;}}L---\left(6\right)$

其中，和分别表示加速度计检验质量质心和卫星体质心在星体坐标系中的加速度矢量，ω=R^TΩR表示在星体坐标系中GRACE绕（X_S,Y_S,Z_S）轴旋转的角速度矩阵

$\mathrm{\&omega;}=\left(\begin{array}{ccc}0& -{\mathrm{\&omega;}}_z& {\mathrm{\&omega;}}_y\\ {\mathrm{\&omega;}}_z& 0& -{\mathrm{\&omega;}}_x\\ -{\mathrm{\&omega;}}_y& {\mathrm{\&omega;}}_x& 0\end{array}\right)$

表示在星体坐标系中GRACE绕（X_S,Y_S,Z_S）轴旋转的角加速度矩阵

$\stackrel{\&CenterDot;}{\mathrm{\&omega;}}=\left(\begin{array}{ccc}0& -{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_z& {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_y\\ {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_z& 0& -{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_x\\ -{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_y& {\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_x& 0\end{array}\right)$

ω²=R^TΩ²R表示在星体坐标系中作用于GRACE的离心角加速度矩阵

${\mathrm{\&omega;}}^2=\left(\begin{array}{ccc}-{\mathrm{\&omega;}}_y^2-{\mathrm{\&omega;}}_z^2& {\mathrm{\&omega;}}_x{\mathrm{\&omega;}}_y& {\mathrm{\&omega;}}_x{\mathrm{\&omega;}}_z\\ {\mathrm{\&omega;}}_y{\mathrm{\&omega;}}_x& -{\mathrm{\&omega;}}_x^2-{\mathrm{\&omega;}}_z^2& {\mathrm{\&omega;}}_y{\mathrm{\&omega;}}_z\\ {\mathrm{\&omega;}}_z{\mathrm{\&omega;}}_x& {\mathrm{\&omega;}}_z{\mathrm{\&omega;}}_y& -{\mathrm{\&omega;}}_x^2-{\mathrm{\&omega;}}_y^2\end{array}\right)$

在星体坐标系中，将公式（8）由矩阵形式改写为矢量形式

${\stackrel{\&CenterDot;\&CenterDot;}{r}}_A^0-{\stackrel{\&CenterDot;\&CenterDot;}{r}}_S^0=\stackrel{\&CenterDot;\&CenterDot;}{L}+\mathrm{\&omega;}\&times;(\mathrm{\&omega;}\&times;L)+2\mathrm{\&omega;}\&times;\stackrel{\&CenterDot;}{L}+\stackrel{\&CenterDot;}{\mathrm{\&omega;}}\&times;L---\left(9\right)$

其中，ω×(ω×L)和分别表示作用于GRACE卫星的惯性离心力和科里奥利力。

在星体坐标系中，GRACE卫星的动力学方程为

${\stackrel{\&CenterDot;\&CenterDot;}{r}}_S^0=\&dtri;V_S+f_S---\left(10\right)$

其中，和fS分别表示作用于GRACE卫星体质心的引力位梯度和非保守力的真值。

在星体坐标系中，GRACE星载加速度计检验质量的动力学方程为

${\stackrel{\&CenterDot;\&CenterDot;}{r}}_A^0=\&dtri;V_A+f_A---\left(11\right)$

其中，和f_A分别表示作用于GRACE星载加速度计检验质量质心的引力位梯度和静电力（等效为非保守力测量值）。

由公式（11）－公式（10）得

${\stackrel{\&CenterDot;\&CenterDot;}{r}}_A^0-{\stackrel{\&CenterDot;\&CenterDot;}{r}}_S^0=(\&dtri;V_A-\&dtri;V_S)+(f_A-f_S)---\left(12\right)$

其中，Δf=f_A-f_S表示由于GRACE卫星体和星载加速度计检验质量的质心偏移引起的加速度计非保守力测量值和真值的偏差；在卫星质心O_S处将星载加速度计检验质量质心的引力位梯度按泰勒展开（取零阶和一阶项）

$\&dtri;V_A\&ap;\&dtri;V_S+{\&dtri;}^2{V}_S\&CenterDot;L---\left(13\right)$ 其中，表示在星体坐标系中引力位的梯度

$K=\left(\begin{array}{ccc}\frac{{\&PartialD;}^2{V}_S}{\&PartialD;x_S^2}& \frac{{\&PartialD;}^2{V}_s}{{\&PartialD;x}_S\&PartialD;y_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_S{\&PartialD;z}_S}\\ \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S{\&PartialD;x}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S^2}& \frac{{\&PartialD;}^2{V}_s}{{\&PartialD;y}_S{\&PartialD;z}_S}\\ \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;x}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;y}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S^2}\end{array}\right)---\left(14\right)$

其中，V_S=GMr表示中心引力位，GM表示地球质量M和万有引力常数G之积，x_S,y_S,z_S分别表示GRACE在星体坐标系中位置矢量的3个分量

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_S^2}=\frac{\mathrm{GM}({6x}_S-r)}{r^4},$ $\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S^2}=\frac{\mathrm{GM}({6y}_S-r)}{r^4},$ $\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S^2}=\frac{\mathrm{GM}({6z}_S-r)}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_s{\&PartialD;y}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S{\&PartialD;x}_S}=\frac{3\mathrm{GM}{x}_S{y}_S}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_s{\&PartialD;z}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;x}_S}=\frac{3\mathrm{GM}{x}_S{z}_S}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_s{\&PartialD;z}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;y}_S}=\frac{3\mathrm{GM}{y}_S{z}_S}{r^4}.$

由公式（9）和公式（12）联合可得，由于GRACE卫星体和星载加速度计检验质量的质心偏移而在加速度计三轴附加的非保守力

$\mathrm{\&Delta;f}=f_A-f_S=\stackrel{\&CenterDot;\&CenterDot;}{L}+\mathrm{\&omega;}\&times;(\mathrm{\&omega;}\&times;L)+2\mathrm{\&omega;}\&times;\stackrel{\&CenterDot;}{L}+\stackrel{\&CenterDot;}{\mathrm{\&omega;}}\&times;L-K\&CenterDot;L---\left(15\right).$

步骤三：地球重力场反演

本发明基于能量守恒原理反演地球重力场，同时引入参考数据，通过测量数据和参考数据的残差建立新型扰动能量观测方程，进一步提高了地球重力场的反演精度。

在地心惯性系中，双星扰动位差观测方程建立如下

T_e12=E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂          （16）

在地心惯性系中，双星参考扰动位差观测方程建立如下

$T_{e12}^o=E_{k12}^o-E_{f12}^o+V_{\mathrm{\&omega;}12}^o-V_{T12}^o-V_{012}^o-E_{012}^o---\left(17\right)$

其中，带上角标“o”的参数表示参考值。

由公式（16）－公式（17）可得

δT_e12=δE_k12-δE_f12+δV_ω12-δV_T12-δV₀₁₂-δE₀₁₂       （18）其中，表示双星相对扰动位差，T_e12(r,θ,λ)表示地球扰动位

$T_{e12}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L_0}{\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}}_{\mathrm{lm}}(\mathrm{\&theta;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\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.;$ 表示双星各自的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎分别表示双星各自位置矢量r₁₍₂₎的三个分量，θ₁和θ₂分别表示双星各自的地心余纬度，λ₁和λ₂分别表示双星各自的地心经度；R_e表示地球的平均半径，L₀表示地球扰动位按球函数展开的最大阶数；表示规格化的Legendre函数，l表示阶数，m表示次数；表示待求的规格化地球引力位系数；利用参考数据同样计算参考地球扰动位

双星相对扰动位差观测方程（18）右边第一项表示双星的扰动动能差， $E_{k12}=\frac{1}{2}({\stackrel{\&CenterDot;}{r}}_2+{\stackrel{\&CenterDot;}{r}}_1)\&CenterDot;\{{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_{12}{e}_{12}+[{\stackrel{\&CenterDot;}{r}}_{12}-({\stackrel{\&CenterDot;}{r}}_{12}\&CenterDot;e_{12})e_{12}\left]\right\}$ 表示双星的动能差，表示双星的相对速度矢量，e₁₂=r₁₂/|r₁₂|表示由第一颗卫星指向第二颗卫星的单位方向矢量，r₁₂=r₂-r₁表示双星的相对位置矢量；利用参考数据同样计算双星的参考动能差第二项表示双星的扰动耗散能差，表示双星的耗散能差；利用参考数据同样计算双星的参考耗散能差第三项表示双星的扰动旋转能差，表示双星的旋转能差，ω_e表示地球的自转角速度，x₁₂=x₂-x₁和y₁₂=y₂-y₁分别表示双星的相对轨道位置，和分别表示双星的相对轨道速度；利用参考数据同样计算双星的参考旋转能差第四项表示双星的扰动三体摄动能差，V_T12表示双星的三体摄动能差，表示双星的参考三体摄动能差；第五项表示双星的扰动中心引力位差，V₀₁₂=GM/r₂-GM/r₁表示双星的中心引力位差；利用参考轨道数据同样计算双星的参考中心引力位差 $V_{012}^o=\mathrm{GM}/{r_2}^o-\mathrm{GM}/{r_1}^o;$ 最后一项 $\mathrm{\&delta;}{E}_{012}=E_{012}-E_{012}^o$ 表示双星的扰动能量积分常数差，E₀₁₂表示双星的能量积分常数差，通过卫星的初始位置和初始速度计算得到；表示双星的参考能量积分常数差，通过卫星的参考初始位置和参考初始速度计算得到。

GRACE星体和加速度计的质心调整精度影响地球重力场精度的具体原理和计算步骤如下：计算作用于卫星体的非保守力矢量和为和将位于星体系的和转化到地心系计算作用于双星的非保守力和将位于星体系的和转化到地心系计算作用于双星的参考非保守力和并代入公式（18）的双星耗散能差E_f12，然后基于最小二乘法解算公式（18），进而获得地球引力位系数在重力卫星飞行中，由于需要利用质心调节装置实时补偿星体和加速度计检验质量的质心偏移，因此会在加速度计三轴测量中引入新的误差源—质心偏移误差，进而影响地球重力场反演的精度。

如图3所示，实线表示德国波兹坦地学研究中心（GFZ）公布的120阶EIGEN-GRACE02S地球重力场模型的精度；圆圈线、三角线、虚线和十字线分别表示基于新型扰动能量法，当质心调整精度设计为0m、5×10^-5m、1×10^-4m和5×10^-4m时，反演地球引力位系数精度。结果表明：第一，在120阶处，当质心调整精度设计为0m时，反演地球引力位系数精度为4.911×10^-10；当质心调整精度分别设计为5×10^-5m、1×10^-4m和5×10^-4m时，反演精度依次降低至6.551×10^-10、1.171×10^-9和4.760×10^-9。第二，以EIGEN-GRACE02S模型的地球引力位系数精度为标准，当质心调整精度设计为(5～10)×10^-5m时，其与K波段测距仪、GPS接收机、SuperSTAR加速度计等GRACE关键载荷的精度指标相匹配，对地球重力场反演精度的影响较小。

图4表示基于未引入参考数据（公式(16)）和引入参考数据（公式(18)）分别反演GRACE地球重力场精度对比。虚线表示基于未引入参考数据反演120阶GRACE地球重力场精度；在120阶处，累计大地水准面精度为1.762×10^-1m。实线表示基于引入参考数据反演120阶GRACE地球重力场精度；在120阶处，累计大地水准面精度为6.113×10^-2m。研究结果表明：基于引入参考数据反演120阶GRACE累计大地水准面精度较基于未引入参考数据的反演精度平均提高2～3倍。因此，引入参考数据可以大幅提高地球重力场反演精度，有利于建立下一代高精度和高空间分辨率地球重力场模型，优化选取卫星重力测量计划中质心调整精度，以及降低星载加速度计和质心调整系统研制的难度。

Claims (1)

1.一种降低重力卫星质心调整误差对地球重力场精度影响的方法，包含下列步骤：

第一步，采集重力卫星测量数据和获取参考数据，

所述测量数据包括，

通过K波段测距仪获取星间速度数据通过星载GPS接收机获取双星轨道位置数据（r₁,r₂）和双星轨道速度数据通过星载加速度计获取作用于双星的非保守力数据（f₀₁,f₀₂），通过星敏感器获取角速度数据Ω=(Ω_x,Ω_y,Ω_z)；

所述参考数据包括，

通过星载GPS接收机获取的双星轨道位置数据（r₁,r₂）和双星轨道速度数据的初始点，利用9阶Runge-Kutta线性单步法和12阶Adams-Cowell线性多步法数值模拟公式获取双星参考轨道位置数据和双星参考轨道速度数据

参考星间速度数据通过参考轨道速度数据计算获得其中，表示相对参考轨道速度矢量，表示第一颗卫星指向第二颗卫星的参考单位矢量，表示相对参考轨道位置矢量；

参考非保守力数据通过DTM2000阻力温度模型计算获得；

通过国际公布模型DE-405、IERS96和CSR4.0联合计算获取作用于双星的三体摄动能差数据V_T12和参考三体摄动能差数据

第二步，计算各个重力卫星由于单个卫星体和星载加速度计检验质量的质心偏移而在加速度计三轴附加的非保守力Δf

$\mathrm{\&Delta;f}=\stackrel{\&CenterDot;\&CenterDot;}{L}+\mathrm{\&omega;}\&times;(\mathrm{\&omega;}\&times;L)+2\mathrm{\&omega;}\&times;\stackrel{\&CenterDot;}{L}+\stackrel{\&CenterDot;}{\mathrm{\&omega;}}\&times;L-K\&CenterDot;L---\left(1\right)$ 其中，L表示卫星体质心和星载加速度计检验质量质心偏移在星体坐标系O_S-X_SY_SZ_S中的位置矢量，和分别表示L对时间的一阶导数和二阶导数；ω=R^TΩR表示在星体坐标系中重力卫星绕星体坐标系（X_S,Y_S,Z_S）轴旋转的角速度矩阵，Ω表示在地心惯性系中重力卫星绕地心坐标系（X_I,Y_I,Z_I）轴旋转的角速度矩阵

$\mathrm{\&omega;}=\left(\begin{array}{ccc}0& -{\mathrm{\&omega;}}_z& {\mathrm{\&omega;}}_y\\ {\mathrm{\&omega;}}_z& 0& -{\mathrm{\&omega;}}_x\\ -{\mathrm{\&omega;}}_y& {\mathrm{\&omega;}}_x& 0\end{array}\right)$

R表示由星体坐标系到地心惯性系的转换矩阵

$R=\left(\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right)$

$R_{11}=\frac{x_I}{\left|r\right|},$ $R_{12}=\frac{y_I}{\left|r\right|},$ $R_{13}=\frac{z_I}{\left|r\right|}$

$R_{21}=\frac{(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)z_I-(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)y_I}{\left|r\right|\left|n\right|}$

$R_{22}=\frac{(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)x_I-(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)z_I}{\left|r\right|\left|n\right|}$

$R_{23}=\frac{(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)y_I-(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)x_I}{\left|r\right|\left|n\right|}$

$R_{31}=\frac{y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I}{\left|n\right|},$ $R_{32}=\frac{z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I}{\left|n\right|},$ $R_{33}=\frac{x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I}{\left|n\right|}$

$\left|r\right|=\sqrt{x_I^2+y_I^2+z_I^2}$

$\left|n\right|=\sqrt{{(y_I{\stackrel{\&CenterDot;}{z}}_I-z_I{\stackrel{\&CenterDot;}{y}}_I)}^2-{(z_I{\stackrel{\&CenterDot;}{x}}_I-x_I{\stackrel{\&CenterDot;}{z}}_I)}^2+{(x_I{\stackrel{\&CenterDot;}{y}}_I-y_I{\stackrel{\&CenterDot;}{x}}_I)}^2}$

其中，x_I,y_I,z_I分别表示重力卫星在地心惯性坐标系中位置矢量的3个分量，分别表示重力卫星在地心惯性坐标系中速度矢量的3个分量；K表示在星体坐标系中地球引力位的梯度

$K=\left(\begin{array}{ccc}\frac{{\&PartialD;}^2{V}_S}{\&PartialD;x_S^2}& \frac{{\&PartialD;}^2{V}_s}{{\&PartialD;x}_S\&PartialD;y_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_S{\&PartialD;z}_S}\\ \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S{\&PartialD;x}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S^2}& \frac{{\&PartialD;}^2{V}_s}{{\&PartialD;y}_S{\&PartialD;z}_S}\\ \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;x}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;y}_S}& \frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S^2}\end{array}\right)$

其中，V_S=GMr表示中心引力位，GM表示地球质量M和万有引力常数G之积，x_S,y_S,z_S分别表示重力卫星在星体坐标系中位置矢量的3个分量

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_S^2}=\frac{\mathrm{GM}({6x}_S-r)}{r^4},$ $\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S^2}=\frac{\mathrm{GM}({6y}_S-r)}{r^4},$ $\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S^2}=\frac{\mathrm{GM}({6z}_S-r)}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_s{\&PartialD;y}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_S{\&PartialD;x}_S}=\frac{3\mathrm{GM}{x}_S{y}_S}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;x}_s{\&PartialD;z}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;x}_S}=\frac{3\mathrm{GM}{x}_S{z}_S}{r^4},$

$\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;y}_s{\&PartialD;z}_S}=\frac{{\&PartialD;}^2{V}_S}{{\&PartialD;z}_S{\&PartialD;y}_S}=\frac{3\mathrm{GM}{y}_S{z}_S}{r^4}.$

第三步，反演地球重力场

计算作用于卫星体的非保守力矢量和为和将位于星体系的和转化到地心系计算作用于双星的非保守力和将位于星体系的和转化到地心系计算作用于双星的参考非保守力 ${f_1}^o=R{f}_{01}^o$ 和 ${f_2}^o=R{f}_{02}^o;$

在地心惯性系中，双星扰动位差观测方程建立如下

T_e12=E_k12-E_f12+V_ω12-V_T12-V₀₁₂-E₀₁₂          （2）

在地心惯性系中，双星参考扰动位差观测方程建立如下

$T_{e12}^o=E_{k12}^o-E_{f12}^o+V_{\mathrm{\&omega;}12}^o-V_{T12}^o-V_{012}^o-E_{012}^o---\left(3\right)$ 其中，带上角标“o”的参数表示参考值；

由公式（2）－公式（3）得到双星相对扰动位差观测方程

δT_e12=δE_k12-δE_f12+δV_ω12-δV_T12-δV₀₁₂-δE₀₁₂        （4）其中，表示双星相对扰动位差，T_e12(r,θ,λ)表示地球扰动位

$T_{e12}(r,\mathrm{\&theta;},\mathrm{\&lambda;})=\frac{\mathrm{GM}}{R_e}\underset{l=2}{\overset{L_0}{\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}}_{\mathrm{lm}}(\mathrm{\&theta;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\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.;$ 表示双星各自的地心半径，x₁₍₂₎,y₁₍₂₎,z₁₍₂₎分别表示双星各自位置矢量r₁₍₂₎的三个分量，θ₁和θ₂分别表示双星各自的地心余纬度，λ₁和λ₂分别表示双星各自的地心经度；R_e表示地球的平均半径，L₀表示地球扰动位按球函数展开的最大阶数；表示规格化的Legendre函数，l表示阶数，m表示次数；表示待求的规格化地球引力位系数；利用参考数据同样计算参考地球扰动位

双星相对扰动位差观测方程（4）右边第一项表示双星的扰动动能差， $E_{k12}=\frac{1}{2}({\stackrel{\&CenterDot;}{r}}_2+{\stackrel{\&CenterDot;}{r}}_1)\&CenterDot;\{{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_{12}{e}_{12}+[{\stackrel{\&CenterDot;}{r}}_{12}-({\stackrel{\&CenterDot;}{r}}_{12}\&CenterDot;e_{12})e_{12}\left]\right\}$ 表示双星的动能差，表示双星的相对速度矢量，e₁₂=r₁₂/|r₁₂|表示由第一颗卫星指向第二颗卫星的单位方向矢量，r₁₂=r₂-r₁表示双星的相对位置矢量；利用参考数据同样计算双星的参考动能差

第二项表示双星的扰动耗散能差，表示双星的耗散能差；利用参考数据同样计算双星的参考耗散能差

第三项表示双星的扰动旋转能差，表示双星的旋转能差，ω_e表示地球的自转角速度，x₁₂=x₂-x₁和y₁₂=y₂-y₁分别表示双星的相对轨道位置，和分别表示双星的相对轨道速度；利用参考数据同样计算双星的参考旋转能差

第四项表示双星的扰动三体摄动能差，V_T12表示双星的三体摄动能差，表示双星的参考三体摄动能差；

第五项表示双星的扰动中心引力位差，V₀₁₂=GM/r₂-GM/r₁表示双星的中心引力位差；利用参考轨道数据同样计算双星的参考中心引力位差 $V_{012}^o=\mathrm{GM}/{r_2}^o-\mathrm{GM}/{r_1}^o;$

最后一项表示双星的扰动能量积分常数差，E₀₁₂表示双星的能量积分常数差，通过卫星的初始位置和初始速度计算得到；表示双星的参考能量积分常数差，通过卫星的参考初始位置和参考初始速度计算得到；

基于最小二乘法解算双星相对扰动位差观测方程（4），进而获得地球引力位系数