基于载荷误差分析原理的卫星重力反演方法

摘要

本发明涉及一种精密测量地球重力场的方法，特别是一种基于载荷误差分析原理精确建立K波段测距仪的星间距离误差、GPS接收机的卫星轨道位置误差和轨道速度误差、以及星载加速度计的非保守力误差联合影响累计大地水准面精度的误差模型，进而精确和快速反演地球重力场的方法；该方法地球重力场反演精度高，保证计算精度的前提下有效提高反演速度，卫星重力反演过程简单，计算机性能要求低，卫星观测方程物理含义明确；基于载荷误差分析原理的卫星重力反演方法是计算高精度和高空间分辨率地球重力场的有效方法。

说明书

一、技术领域

本发明涉及卫星重力学、空间大地测量学、航空航天等交叉技术领域，特 别是涉及一种基于载荷误差分析原理精确和快速反演地球重力场的方法。

二、背景技术

21世纪是人类利用SST-HL/LL（Satellite-to-Satellite Tracking in the  High-Low/Low-Low Mode）和SGG（Satellite Gravity Gradiometry）提升对数字 地球认知能力的新纪元。地球重力场及其时变反映地球表层及内部物质的空间 分布、运动和变化，同时决定着大地水准面的起伏和变化。因此，确定地球重 力场的精细结构及其时变不仅是大地测量学、地球物理学、地震学、海洋学、 空间科学、国防建设等的需求，同时也将为寻求资源、保护环境和预测灾害提 供重要的信息资源。

GRACE（Gravity Recovery and Climate Experiment）双星采用近圆和近极地 轨道设计，由美国宇航局（NASA）和德国航天局（DLR）共同研制开发。GRACE 利用K波段测距仪高精度测量星间距离，利用高轨GPS（Global Positioning  System）卫星对低轨双星精密跟踪定位，利用高精度SuperSTAR加速度计测量 作用于双星的非保守力。GRACE系统既包含两组SST-HL，同时以差分原理测 定两个低轨卫星之间的相互运动，因此得到的静态和动态全球重力场精度比 CHAMP（Challenging Minisatellite Payload）至少高一个数量级，同时为将来 GOCE（Gravity Field and Steady-State Ocean Circulation Explorer）卫星重力梯度 测量奠定了坚实基础。

早在20世纪60年代，Baker首次提出了利用SST恢复地球重力场的重要思 想。自此以后，国际大地测量学界的许多学者都积极投身于地球重力场恢复的 方法与算法的理论研究和数值计算之中。在众多方法中，按照卫星观测方程的 建立和求解的不同可分为解析法和数值法。解析法是指通过分析地球重力场和 卫星观测数据的关系建立卫星观测方程模型，进而估计地球重力场的精度。解 析法的优点是卫星观测方程物理含义明确，易于误差分析且可快速求解高阶地 球重力场；缺点是由于在建立卫星观测方程模型时作了不同程度的近似，因此 求解精度较低。数值法是指通过分析地球引力位系数和卫星观测数据的关系建 立卫星观测方程，并通过最小二乘法拟合出地球引力位系数。数值法的优点是 地球重力场求解精度较高；缺点是求解速度较慢且对计算机要求较高。不同于 以前的技术，本发明基于载荷误差分析法建立了K波段测距仪的星间距离、GPS 接收机的轨道位置和轨道速度、以及加速度计的非保守力误差联合影响累积大 地水准面的误差模型，基于关键载荷精度指标的匹配关系论证了误差模型的可 靠性，基于美国宇航局喷气推进实验室（NASA-JPL）公布的2009年的 GRACE-Level-1B实测误差数据，有效和快速地反演了120阶GRACE地球重力 场精度。

三、发明内容

本发明的目的是：基于载荷误差分析法较大程度优化地球重力场反演速度， 而且进一步提高地球重力场反演精度。

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

基于载荷误差分析原理的卫星重力反演方法包含下列步骤：

步骤一：卫星关键载荷数据采集

1.1）通过星载K波段测距仪获取星间距离误差数据δρ₁₂；

1.2）通过星载GPS接收机获取轨道位置误差数据δr和轨道速度误差数据

1.3）通过星载加速度计获取非保守力误差数据δf；

步骤二：关键载荷误差模型建立

2.1）K波段测距仪的星间距离误差模型

基于能量守恒定律，卫星观测方程可表示为

$\frac{1}{2}{\stackrel{·}{r}}^2=V+C---\left(1\right)$

其中，表示卫星的瞬时速度，表示卫星的平均速度，GM 表示地球质量M和万有引力常数G之积，r表示由卫星质心到地心之间的距离， 表示由地球扰动位引起的速度变化；V=V₀+T表示地球引力位，V₀表示中心 引力位，T表示扰动位；C表示能量积分常数；公式（1）可变形为

$\frac{1}{2}{({\stackrel{·}{r}}_0+\Delta \stackrel{·}{r})}^2=V_0+T+C---\left(2\right)$

由于忽略二阶小量且公式（2）可变形为

$T={\stackrel{·}{r}}_0\Delta \stackrel{·}{r}---\left(3\right)$

扰动位方差和速度变化方差的关系为

${\sigma }^2\left(\mathrm{\delta T}\right)={\stackrel{·}{r}}_0^2{\sigma }^2\left(\delta \stackrel{·}{r}\right)---\left(4\right)$

表示K波段测距仪的星间速度，表示星间速度的变化量； 星间速度的方差表示为

${\sigma }^2\left(\delta {\stackrel{·}{\rho }}_{12}\right)\approx 2[{\sigma }^2\left(\delta \stackrel{·}{r}\right)-\mathrm{cov}(\Delta {\stackrel{·}{r}}_1,\Delta {\stackrel{·}{r}}_2)]---\left(5\right)$

其中，表示协方差函数，$\mathrm{cov}(\Delta {\stackrel{·}{r}}_1,{\Delta \stackrel{·}{r}}_2)=\underset{l=2}{\overset{L}{\Sigma }}{\sigma }_l^2\left(\delta \stackrel{·}{r}\right)P_l\left(\cos \theta \right),$P_l(cosθ) 表示Legendre函数，l表示阶数，θ表示地心角；公式（5）可变形为

${\sigma }_l^2\left(\delta {\stackrel{·}{\rho }}_{12}\right)\approx 2{\sigma }_l^2\left(\delta \stackrel{·}{r}\right)[1-P_l\left(\cos \theta \right)]---\left(6\right)$

由于因此，公式（6）可表示为

${\sigma }_l^2\left(\delta {\rho }_{12}\right)\approx 2{\left(\mathrm{\Delta t}\right)}^2{\sigma }_l^2\left(\delta \stackrel{·}{r}\right)[1-P_l\left(\cos \theta \right)]---\left(7\right)$

其中，δρ₁₂表示K波段测距仪的星间距离误差，Δt表示采样间隔；

地球扰动位T(r,φ，λ)表示为

$T(r,\phi ,\lambda )=\frac{\mathrm{GM}}{r}\underset{l=2}{\overset{L}{\Sigma }}\underset{m=0}{\overset{l}{\Sigma }}\left[{\left(\frac{R_e}{r}\right)}^l(C_{\mathrm{lm}}\cos \mathrm{m\lambda }+S_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\sin \phi \right)\right]---\left(8\right)$

其中，φ表示地心纬度，λ表示地心经度，R_e表示地球的平均半径，L表示地球 扰动位按球函数展开的最大阶数；表示规格化的Legendre函数，m表 示次数；C_lm,S_lm表示待求的规格化地球引力位系数；

地球扰动位的方差表示为

${\sigma }_l^2\left(\mathrm{\delta T}\right)=\underset{m=0}{\overset{l}{\Sigma }}{[\frac{1}{4\pi }\int \int \mathrm{\delta T}(r,\phi ,\lambda ){\bar{Y}}_{\mathrm{lm}}(\phi ,\lambda )\cos \mathrm{\phi d\phi d\lambda }]}^2---\left(9\right)$

其中，${\bar{Y}}_{\mathrm{lm}}(\phi ,\lambda )={\bar{P}}_{l\left|m\right|}\left(\sin \phi \right)Q_m\left(\lambda \right),$$Q_m\left(\lambda \right)=\left(\begin{array}{cc}\cos \mathrm{m\lambda }& m\ge 0\\ \sin \left|m\right|\lambda & m<0\end{array}\right);$

基于球谐函数的正交性，公式（9）可化简为

${\sigma }_l^2\left(\mathrm{\delta T}\right)={\left(\frac{\mathrm{GM}}{R_e}\right)}^2{\left(\frac{R_e}{r}\right)}^{2l+2}\underset{m=0}{\overset{l}{\Sigma }}(\delta C_{\mathrm{lm}}^2+\delta S_{\mathrm{lm}}^2)---\left(10\right)$

其中，δC_lm,δS_lm表示地球引力位系数精度；

大地水准面高的方差为

${\sigma }_l^2\left(\delta N_{{\rho }_{12}}\right)=R_e^2\underset{m=0}{\overset{l}{\Sigma }}(\delta C_{\mathrm{lm}}^2+\delta S_{\mathrm{lm}}^2)---\left(11\right)$

联合公式（10）和公式（11），可得和的关系式

${\sigma }_l^2\left(\delta N_{{\rho }_{12}}\right)=R_e^2{\left(\frac{R_e}{\mathrm{GM}}\right)}^2{\left(\frac{r}{R_e}\right)}^{2l+2}{\sigma }_l^2\left(\mathrm{\delta T}\right)---\left(12\right)$

联合公式（4）、（7）和（12），可得累积大地水准面误差和星间距离误差之间的 关系式

$\delta N_{{\rho }_{12}}=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\delta {\rho }_{12}\right)\right\}}---\left(13\right)$

2.2）GPS接收机的轨道位置误差模型

卫星向心加速度和瞬时速度的关系式表示为

$\stackrel{··}{r}=\frac{{\stackrel{·}{r}}^2}{r}---\left(\text{14}\right)$

其中，表示在星星连线方向的投影；公式（14）可变形 为

${\stackrel{··}{r}}_{{\rho }_{12}}=\frac{\sin (\theta /2)}{r}{\stackrel{·}{r}}^2---\left(15\right)$

在公式（15）两边同时微分可得

$d{\stackrel{··}{r}}_{{\rho }_{12}}=\frac{2\stackrel{·}{r}\sin (\theta /2)}{r}d\stackrel{·}{r}---\left(16\right)$

由于且忽略二阶小量在公式（16）两边同乘时间t可得

$d{\stackrel{·}{r}}_{{\rho }_{12}}=\sqrt{\frac{4\mathrm{GM}{\sin }^2(\theta /2)}{r^3}}\mathrm{dr}---\left(17\right)$

基于公式（17）且星间距离误差δρ₁₂和轨道位置误差δr的关系表示 为

$\delta {\rho }_{12}=\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}---\left(18\right)$

将公式（18）代入公式（13）可得累积大地水准面误差和轨道位置误差之间的 关系式

$\delta N_r=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}\right)\right\}}---\left(19\right)$

2.3）GPS接收机的轨道速度误差模型

卫星加速度在星星连线方向投影和卫星加速度之间的关系为

${\stackrel{··}{r}}_{{\rho }_{12}}=\stackrel{··}{r}\sin (\theta /2)---\left(\text{20}\right)$

其中，表示K波段测距仪的星间加速度；在公式（20）两边同 时微分并乘时间t可得

$d{\stackrel{·}{\rho }}_{12}=2\sin (\theta /2)d\stackrel{·}{r}---\left(21\right)$

基于公式（21）且星间距离误差δρ₁₂和轨道速度误差之间的关系 为

$\delta {\rho }_{12}=2\Delta t\sin (\theta /2)\delta \stackrel{·}{r}---\left(22\right)$

将公式（22）代入公式（13），可得累积大地水准面误差和轨道速度误差之间的 关系式

$\delta N_{\stackrel{·}{r}}=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(2\Delta t\sin (\theta /2)\delta \stackrel{·}{r}\right)\right\}}---\left(\text{23}\right)$

2.4）加速度计的非保守力误差模型

星间速度误差和非保守力误差δf的关系表示为

$\delta {\stackrel{·}{\rho }}_{12}=\sqrt{{\int \left(\mathrm{\delta f}\right)}^2\mathrm{dt}}---\left(24\right)$

由于公式（24）表示如下

$\delta {\rho }_{12}=\mathrm{\Delta t}\sqrt{\int {\left(\mathrm{\delta f}\right)}^2\mathrm{dt}}---\left(25\right)$

将公式（25）代入公式（13）可得累积大地水准面误差和非保守力误差之间的 关系式

$\delta N_f=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\mathrm{\Delta t}\sqrt{\int {\left(\mathrm{\delta f}\right)}^2\mathrm{dt}}\right)\right\}}---\left(26\right)$

2.5）关键载荷联合误差模型

联合公式（13）、（19）、（23）和（26），可得星间距离、轨道位置、轨道速 度和非保守力误差联合影响累积大地水准面的误差模型

$\delta N_c=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\mathrm{\delta \eta }\right)\right\}}---\left(27\right)$

其中，

$\mathrm{\delta \eta }=\sqrt{{\sigma }_l^2\left(\delta {\rho }_{12}\right)+{\sigma }_l^2\left(\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}\right)+{\sigma }_l^2\left(2\Delta t\sin (\theta /2)\delta \stackrel{·}{r}\right)+{\sigma }_l^2\left(\mathrm{\Delta t}\sqrt{\int {\left(\mathrm{\delta f}\right)}^2\mathrm{dr}}\right)},$表示K波段测距仪的星间距离方差，${\sigma }_l^2\left(\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}\right)$表示 GPS接收机的轨道位置方差，表示GPS接收机的轨道速度方 差，表示加速度计的非保守力方差；

步骤三：地球重力场反演

基于载荷误差分析法，利用星间距离误差数据δρ₁₂、轨道位置误差数据δr和 轨道速度误差数据以及非保守力误差数据δf反演累积大地水准面误差的过 程如下：

第一，首先以0.5°×0.5°为网格分辨率，在地球表面的经度0°~360°和纬度 -90°~90°范围内绘制网格；其次，在地球表面的轨迹点位置依次加入δη；最后， 将分布于地球表面的δη平均归算于划分的网格点δη(φ，λ)；

第二，将δη(φ，λ)按球谐函数展开为

$\mathrm{\delta \eta }(\phi ,\lambda )=\underset{l=0}{\overset{L}{\Sigma }}\underset{m=0}{\overset{l}{\Sigma }}\left[(C_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}\cos \mathrm{m\lambda }+S_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\sin \phi \right)\right]---\left(28\right)$

其中，表示δη(φ，λ)按球函数展开的系数

$(C_{{\mathrm{\delta \eta }}_{\mathrm{lm}}},S_{{\mathrm{\delta \eta }}_{\mathrm{lm}}})=\frac{1}{4\pi }\int \int \left[\mathrm{\delta \eta }(\phi ,\lambda ){\bar{Y}}_{\mathrm{lm}}(\phi ,\lambda )\cos \mathrm{\phi d\phi d\lambda }\right]---\left(29\right)$

δη在各阶处的方差表示为

${\sigma }_l^2\left(\mathrm{\delta \eta }\right)=\underset{m=0}{\overset{l}{\Sigma }}(C_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}^2+S_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}^2)---\left(30\right)$

基于公式（29）计算将公式（30）代入公式（27），可有效和快速反 演地球重力场精度。

本发明是基于载荷误差分析法有利于精确和快速反演地球重力场的特点而 设计的，优点为：

1）地球重力场反演精度高；

2）保证计算精度的前提下有效提高反演速度；

3）卫星重力反演过程简单；

4）计算机性能要求低；

5）卫星观测方程物理含义明确。

四、附图说明

图1表示GRACE双星在轨飞行示意图。

图2表示GRACE-A卫星的全球轨道分布图。

图3表示GRACE卫星关键载荷误差在全球地表的分布(单位:μm/s)。

图4表示GRACE关键载荷匹配精度指标论证。

图5表示基于载荷误差分析法反演GRACE累积大地水准面误差。

五、具体实施方式

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

基于载荷误差分析原理的卫星重力反演方法：

步骤一：卫星关键载荷数据采集

1.1）通过星载K波段测距仪获取星间距离误差数据δρ₁₂；

1.2）通过星载GPS接收机获取轨道位置误差数据δr和轨道速度误差数据

1.3）通过星载加速度计获取非保守力误差数据δf。

步骤二：关键载荷误差模型建立

2.1）K波段测距仪的星间距离误差模型

基于能量守恒定律，卫星观测方程可表示为

$\frac{1}{2}{\stackrel{·}{r}}^2=V+C---\left(31\right)$

其中，表示卫星的瞬时速度，表示卫星的平均速度，GM 表示地球质量M和万有引力常数G之积，r表示由卫星质心到地心之间的距离， 即r取平均值，r＝R_e(地球平均半径)+H(平均轨道高度)，表示由扰动位引起 的速度变化；V=V₀+T表示地球引力位，V₀表示中心引力位，T表示扰动位；C 表示能量积分常数。公式（31）可变形为

$\frac{1}{2}{({\stackrel{·}{r}}_0+\Delta \stackrel{·}{r})}^2=V_0+T+C---\left(32\right)$

由于忽略二阶小量（近似程度约10^-10）且公式（32）可变形 为

$T={\stackrel{·}{r}}_0\Delta \stackrel{·}{r}---\left(33\right)$

扰动位方差和速度变化方差的关系为

${\sigma }^2\left(\mathrm{\delta T}\right)={\stackrel{·}{r}}_0^2{\sigma }^2\left(\delta \stackrel{·}{r}\right)---\left(34\right)$

如图1所示，O_I-X_I Y_I Z_I表示地心惯性系；θ表示地心角，对于GRACE双星， θ=2°；表示K波段测距仪的星间速度，表示星间速度的变化 量。星间速度的方差表示为

${\sigma }^2\left(\delta {\stackrel{·}{\rho }}_{12}\right)\approx 2[{\sigma }^2\left(\delta \stackrel{·}{r}\right)-\mathrm{cov}(\Delta {\stackrel{·}{r}}_1,\Delta {\stackrel{·}{r}}_2)]---\left(35\right)$

其中，表示协方差函数，$\mathrm{cov}(\Delta {\stackrel{·}{r}}_1,{\stackrel{·}{r}}_2)=\underset{l=2}{\overset{L}{\Sigma }}{\sigma }_l^2\left(\delta \stackrel{·}{r}\right)P_l\left(\cos \theta \right),$P_l(cosθ) 表示Legendre函数，l表示阶数。公式（35）可变形为

${\sigma }_l^2\left(\delta {\stackrel{·}{\rho }}_{12}\right)\approx 2{\sigma }_l^2\left(\delta \stackrel{·}{r}\right)[1-P_l\left(\cos \theta \right)]---\left(36\right)$

由于因此，公式（36）可表示为

${\sigma }_l^2\left(\delta {\rho }_{12}\right)\approx 2{\left(\mathrm{\Delta t}\right)}^2{\sigma }_l^2\left(\delta \stackrel{·}{r}\right)[1-P_l\left(\cos \theta \right)]---\left(37\right)$

其中，δρ₁₂表示K波段测距仪的星间距离误差，Δt表示采样间隔。

地球扰动位T(r,φ，λ)表示为

$T(r,\phi ,\lambda )=\frac{\mathrm{GM}}{r}\underset{l=2}{\overset{L}{\Sigma }}\underset{m=0}{\overset{l}{\Sigma }}\left[{\left(\frac{R_e}{r}\right)}^l(C_{\mathrm{lm}}\cos \mathrm{m\lambda }+S_{\mathrm{lm}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\sin \phi \right)\right]---\left(38\right)$

其中，φ表示地心纬度，λ表示地心经度，R_e表示地球的平均半径，L表示地球 扰动位按球函数展开的最大阶数；表示规格化的Legendre函数，m表 示次数；C_lm,S_lm表示待求的规格化地球引力位系数。

地球扰动位的方差表示为

${\sigma }_l^2\left(\mathrm{\delta T}\right)=\underset{m=0}{\overset{l}{\Sigma }}{[\frac{1}{4\pi }\int \int \mathrm{\delta T}(r,\phi ,\lambda ){\bar{Y}}_{\mathrm{lm}}(\phi ,\lambda )\cos \mathrm{\phi d\phi d\lambda }]}^2---\left(39\right)$

其中，${\bar{Y}}_{\mathrm{lm}}(\phi ,\lambda )={\bar{P}}_{l\left|m\right|}\left(\sin \phi \right)Q_m\left(\lambda \right),$$Q_m\left(\lambda \right)=\left(\begin{array}{cc}\cos \mathrm{m\lambda }& m\ge 0\\ \sin \left|m\right|\lambda & m<0\end{array}\right).$

基于球谐函数的正交性，公式（39）可化简为

${\sigma }_l^2\left(\mathrm{\delta T}\right)={\left(\frac{\mathrm{GM}}{R_e}\right)}^2{\left(\frac{R_e}{r}\right)}^{2l+2}\underset{m=0}{\overset{l}{\Sigma }}(\delta C_{\mathrm{lm}}^2+\delta S_{\mathrm{lm}}^2)---\left(40\right)$

其中，δC_lm,δS_lm表示地球引力位系数精度。

大地水准面高的方差为

${\sigma }_l^2\left(\delta N_{{\rho }_{12}}\right)=R_e^2\underset{m=0}{\overset{l}{\Sigma }}(\delta C_{\mathrm{lm}}^2+\delta S_{\mathrm{lm}}^2)---\left(41\right)$

联合公式（40）和（41），可得和的关系式

${\sigma }_l^2\left(\delta N_{{\rho }_{12}}\right)=R_e^2{\left(\frac{R_e}{\mathrm{GM}}\right)}^2{\left(\frac{r}{R_e}\right)}^{2l+2}{\sigma }_l^2\left(\mathrm{\delta T}\right)---\left(42\right)$

联合公式（34）、（37）和（42），可得累积大地水准面误差和星间距离误差之间 的关系式

$\delta N_{{\rho }_{12}}=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\delta {\rho }_{12}\right)\right\}}---\left(43\right)$

2.2）GPS接收机的轨道位置误差模型

如图1所示，卫星向心加速度和瞬时速度的关系式表示为

$\stackrel{··}{r}=\frac{{\stackrel{·}{r}}^2}{r}---\left(\text{44}\right)$

其中，表示在星星连线方向的投影。公式（44）可变形 为

${\stackrel{··}{r}}_{{\rho }_{12}}=\frac{\sin (\theta /2)}{r}{\stackrel{·}{r}}^2---\left(45\right)$

在公式（45）两边同时微分可得

$d{\stackrel{··}{r}}_{{\rho }_{12}}=\frac{2\stackrel{·}{r}\sin (\theta /2)}{r}d\stackrel{·}{r}---\left(46\right)$

由于且忽略二阶小量（近似程度约10^-10），在公式（46）两边 同乘时间t可得

$d{\stackrel{·}{r}}_{{\rho }_{12}}=\sqrt{\frac{4\mathrm{GM}{\sin }^2(\theta /2)}{r^3}}\mathrm{dr}---\left(47\right)$

基于公式（47）且星间距离误差δρ₁₂和轨道位置误差δr的关系表示 为

${\delta }_{{\rho }_{12}}=\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}---\left(48\right)$

将公式（48）代入公式（43）可得累积大地水准面误差和轨道位置误差之间的 关系式

$\delta N_r=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}\right)\right\}}---\left(49\right)$

2.3）GPS接收机的轨道速度误差模型

如图1所示，卫星加速度在星星连线方向投影和卫星加速度之间的关 系为

${\stackrel{··}{r}}_{{\rho }_{12}}=\stackrel{··}{r}\sin (\theta /2)---\left(50\right)$

其中，表示K波段测距仪的星间加速度。在公式（50）两边同 时微分并乘时间t可得

$d{\stackrel{·}{\rho }}_{12}=2\sin (\theta /2)d\stackrel{·}{r}---\left(51\right)$

基于公式（51）且星间距离误差δρ₁₂和轨道速度误差之间的关系 为

$\delta {\rho }_{12}=2\Delta t\sin (\theta /2)\delta \stackrel{·}{r}.---\left(52\right)$

将公式（52）代入（43），可得累积大地水准面误差和轨道速度误差之间的关系 式

$\delta N_{\stackrel{·}{r}}=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(2\Delta t\sin (\theta /2)\delta \stackrel{·}{r}\right)\right\}}.---\left(\text{53}\right)$

2.4）加速度计的非保守力误差模型

由于双星受到的主要非保守力和星间速度近似同向，且非保守力通常表现 为累积误差特性，据平方误差积分准则，星间速度误差和非保守力误差δf的 关系表示为

$\delta {\stackrel{·}{\rho }}_{12}=\sqrt{{\int \left(\mathrm{\delta f}\right)}^2\mathrm{dt}}---\left(51\right)$

由于公式（54）表示为

$\delta {\rho }_{12}=\mathrm{\Delta t}\sqrt{{\int \left(\mathrm{\delta f}\right)}^2\mathrm{dt}}---\left(55\right)$

将公式（55）代入公式（43）可得累积大地水准面误差和非保守力误差之间的 关系式

$\delta N_f=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\mathrm{\Delta t}\sqrt{\int {\left(\mathrm{\delta f}\right)}^2\mathrm{dt}}\right)\right\}}---\left(56\right)$

2.5）关键载荷联合误差模型

联合公式（43）、（49）、（53）和（56），可得星间距离、轨道位置、轨道速 度和非保守力误差联合影响累积大地水准面的误差模型

$\delta N_c=R_e\sqrt{\underset{l=2}{\overset{L}{\Sigma }}\left\{\frac{1}{2{\left(\mathrm{\Delta t}\right)}^2[1-P_l\left(\cos \theta \right)]}\frac{R_e}{\mathrm{GM}}{\left(\frac{r}{R_e}\right)}^{2l+1}{\sigma }_l^2\left(\mathrm{\delta \eta }\right)\right\}}---\left(57\right)$

其中，

$\mathrm{\delta \eta }=\sqrt{{\sigma }_l^2\left(\delta {\rho }_{12}\right)+{\sigma }_l^2\left(\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}\right)+{\sigma }_l^2\left(2\Delta t\sin (\theta /2)\delta \stackrel{·}{r}\right)+{\sigma }_l^2\left(\mathrm{\Delta t}\sqrt{\int {\left(\mathrm{\delta f}\right)}^2\mathrm{dr}}\right)},$表示K波段测距仪的星间距离方差，${\sigma }_l^2\left(\sqrt{\frac{4\mathrm{GM}{\left(\mathrm{\Delta t}\right)}^2{\sin }^2(\theta /2)}{r^3}}\mathrm{\delta r}\right)$表示 GPS接收机的轨道位置方差，表示GPS接收机的轨道速度方 差，表示加速度计的非保守力方差。该模型不仅加入K波段 测量系统的星间速度、GPS接收机的轨道位置和加速度计的非保守力误差，还 考虑GPS接收机的轨道速度误差对地球重力场精度的影响，基于此联合误差模 型可更高精度和高空间分辨率地反演了GRACE地球重力场。

步骤三：地球重力场反演

基于载荷误差分析法，利用K波段测距仪的星间距离误差数据δρ₁₂、GPS 接收机的轨道位置误差数据δr和轨道速度误差数据以及加速度计的非保守 力误差数据δf反演累积大地水准面误差的过程如下

第一，首先以0.5°×0.5°为网格分辨率，在地球表面的经度（0°~360°）和纬 度（-90°~90°）范围内绘制网格；其次，按照GRACE卫星轨道（如图2所示） 在地球表面的轨迹点位置依次加入δη；最后，如图3所示，将分布于地球表面 的δη平均归算于划分的网格点δη(φ，λ)处，其中横坐标和纵坐标分别表示经度 和纬度，颜色表示平均归算于网格点处的误差值δη(φ，λ)的大小（μm/s）。

第二，将δη(φ，λ)按球谐函数展开为

$\mathrm{\delta \eta }(\phi ,\lambda )=\underset{l=0}{\overset{L}{\Sigma }}\underset{m=0}{\overset{l}{\Sigma }}\left[(C_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}\cos \mathrm{m\lambda }+S_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}\sin \mathrm{m\lambda }){\bar{P}}_{\mathrm{lm}}\left(\sin \phi \right)\right]---\left(58\right)$

其中，表示δη(φ，λ)按球函数展开的系数

$(C_{{\mathrm{\delta \eta }}_{\mathrm{lm}}},S_{{\mathrm{\delta \eta }}_{\mathrm{lm}}})=\frac{1}{4\pi }\int \int \left[\mathrm{\delta \eta }(\phi ,\lambda ){\bar{Y}}_{\mathrm{lm}}(\phi ,\lambda )\cos \mathrm{\phi d\phi d\lambda }\right]---\left(59\right)$

δη在各阶处的方差表示为

${\sigma }_l^2\left(\mathrm{\delta \eta }\right)=\underset{m=0}{\overset{l}{\Sigma }}(C_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}^2+S_{{\mathrm{\delta \eta }}_{\mathrm{lm}}}^2)---\left(60\right)$

基于公式（59）计算将公式（60）代入公式（57），可有效和快 速反演全球重力场精度。如图4所示，实线、圆圈线、十字线和虚线分别表示 单独引入K波段测距仪的星间距离误差1×10^-5m、GPS接收机的轨道位置误差 3×10^-2m和轨道速度误差3×10^-5m/s、以及加速度计的非保守力误差3×10^-10m/s²反演累积大地水准面误差。基于GRACE关键载荷精度指标的匹配关系，据图中 4条曲线在各阶处的符合性可验证本发明基于载荷误差分析法建立的误差模型 是可靠的。

如图5所示，虚线表示德国波兹坦地学研究中心（GFZ）公布的120阶 EIGEN-GRACE02S全球重力场模型的实测精度，在120阶处反演累计大地水准 面精度为18.938cm；实线表示基于关键载荷联合误差模型反演累计大地水准面 的精度，在120阶处累计大地水准面精度为18.825cm。通过两条曲线在各阶处 的符合性可知，载荷误差分析法是反演高精度和高空间分辨率全球重力场的有 效方法之一。

以上具体实施方式仅为本发明的一种实施示例，其描述较为具体和详细， 但不能因此而理解为对本发明专利范围的限制。其具体实施步骤顺序和模型参 数可根据实际需要进行相应的调整。应当指出的是，对于本领域的普通技术人 员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都 属于本发明的保护范围。