CN103064128B - 基于星间距离误差模型的地球重力场恢复方法 - Google Patents

Abstract

本发明涉及一种地球重力场精密测量方法，特别是一种基于星间距离误差模型原理的地球重力场恢复方法；该方法基于星间距离误差影响累计大地水准面精度的关系建立星间距离误差模型，进而使用该星间距离误差模型来精确和快速恢复当前GRACE和下一代GRACE-II地球重力场。该方法对地球重力场恢复精度高，较大程度提高解算速度，易于开展高阶重力场误差分析，卫星观测方程物理含义明确，计算机性能要求低。因此，星间距离误差模型法是恢复高精度和高空间分辨率地球重力场的有效方法。

Description

基于星间距离误差模型的地球重力场恢复方法

技术领域

本发明涉及卫星重力学、大地测量学、空间科学等交叉技术领域，特别是涉及一种基于星载激光干涉测距仪的星间距离误差影响累计大地水准面精度的关系来建立激光干涉测距仪星间距离误差模型，进而使用这种星间距离误差模型来精确和快速恢复当前GRACE和下一代GRACE-II地球重力场的方法。

背景技术

地球重力场及其时变反映地球表层及内部物质的空间分布、运动和变化，同时决定着大地水准面的起伏和变化。因此，确定地球重力场的精细结构及其时变不仅是大地测量学、海洋学、空间科学等的需求，同时也将为寻求资源、保护环境和预测灾害提供重要的信息资源。重力恢复和气候实验卫星（GRACE）的成功发射以及下一代GRACE-II卫星的即将发射昭示着人类将迎来一个前所未有的卫星重力探测新时代。基于GRACE双星高精度感测中长波地球重力场的优秀表现，美国宇航局（NASA）提出了又一项专用于中短波地球重力场精密探测的GRACE-II未来卫星计划。如图1所示，GRACE-II双星预期采用近圆、近极地和低轨道设计，利用激光干涉测距仪高精度感测星间距离（测量精度10^-8m）。因此，下一代GRACE-II得到的静态和时变地球重力场精度比目前GRACE至少高一个数量级。

在卫星重力恢复的众多方法中，按引力位系数解算方法的差异可分为空域法和时域法。空域法是指不直接处理空间位置相对不规则的卫星轨道采样点的观测值，而将这些观测值归算到以卫星平均轨道高度为半径的球面上利用快速傅立叶变换（FFT）进行网格化处理，将问题转化为某类型边值问题的解，如准解析法、最小二乘配置法等属于空域法的范畴。优点是因网格点数固定从而方程维数一定，且可以利用FFT方法进行快速批量处理，因此极大地降低了计算量；缺点是在进行网格化处理中作了不同程度的近似计算，且不能对色噪声进行处理。时域法是指将卫星观测数据按时间序列处理，卫星星历值直接表示成地球引力位系数的函数，由最小二乘等方法直接反求引力位系数。优点是直接对卫星观测数据进行处理，不需作任何近似，求解精度较高且能有效处理色噪声；缺点是随着卫星观测数据的增多，观测方程数量剧增，极大地增加了计算量。过去由于地球重力场恢复方法的历史局限性和当时计算机技术发展的限制，为了减少计算量，因此空域法较为盛行，Colombo（1989）、Sanso（1995）、Reguzzoni（2003）、Sharifi（2006）等在此方面开展了广泛研究。然而，由于空域法做了许多人为性的假设，存在许多潜在的弊端且随着近年来计算机技术的飞速发展及各种快速算法的广泛应用，计算量的大小不再是制约地球重力场恢复精度的重要因素，时域法的优点正逐渐体现于卫星重力反演之中，Hanetal.（2002）、Reigber（2002）、SchwintzerandReigber（2002）等学者直接利用时域法反演了高精度的地球重力场。时域法主要包括：Kaula线性摄动法、数值微分法、动力学法、能量守恒法、卫星加速度法等。国内外研究表明，Kaula线性摄动法和数值微分法只适合于求解低阶地球重力场且计算精度较低，因此目前基本上已无人问津，现在最为盛行的是动力学法和能量守恒法。动力学法的优点是求解精度较高；缺点是观测数据运算量较大、求解过程复杂程度较高且反演较高阶重力场（L>100阶）时需要高性能的并行计算机支持；能量守恒法的优点是观测方程物理含义明确且易于地球重力场的敏感度分析，在保证求解精度的前提下计算量大大降低，通常采用PC计算机可完成高阶地球重力场的快速求解；缺点是对卫星速度的测量精度要求较高。

为了有效综合已有卫星重力恢复方法的优点，本发明提出了基于星间距离误差模型精确和快速恢复当前GRACE和下一代GRACE-II地球重力场的新技术，并精确和快速地恢复了120阶GRACE和360阶GRACE-II全球重力场。

发明内容

本发明的目的是：基于星间距离误差模型法较大程度加快计算速度，而且进一步提高当前GRACE和下一代GRACE-II地球重力场恢复的精度。

为达到上述目的，本发明提供了一种星间距离误差模型的地球重力场恢复方法，包括如下步骤：

步骤1：采集重力恢复和气候实验卫星的关键载荷数据：通过星载测距仪获取星间距离误差数据δρ₁₂，通过星载GPS接收机获取轨道位置数据r；

步骤2：通过星间距离误差数据δρ₁₂与累计大地水准面精度的关系，建立星间距离误差模型；

步骤3：基于所述星间距离误差模型，通过所采集的卫星关键载荷数据，对地球重力场进行恢复；其中，所述步骤3包括：

步骤3.1：利用9阶Runge-Kutta线性单步法结合12阶Adams-Cowell线性多步法数值积分公式模拟卫星的星历；

步骤3.2：确定参考球面网格分辨率，在地球表面的经度λ和纬度φ范围内按所确定的参考球面网格分辨率绘制网格，按照卫星轨道在地球表面的轨迹点位置依次加入星间距离误差δρ₁₂(φ，λ)；

步骤3.3：基于所述星间距离误差模型和星间距离误差数据δρ₁₂恢复地球重力场。

本发明是基于激光干涉测距仪星间距离误差模型法有利于精确和快速恢复下一代GRACE-II地球重力场的特点而设计的，优点是：

1）地球重力场恢复精度高；

2）较大程度提高解算速度；

3）易于开展高阶重力场误差分析；

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

5）计算机性能要求低。

附图说明

图1为下一代GRACE-II卫星重力计划的测量原理图。

图2表示T(r,φ，λ),和的功率谱。

图3表示GRACE-II卫星轨道。

图4表示基于星间距离误差模型法恢复GRACE和GRACE-II累计大地水准面精度对比。

具体实施方式

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

基于星间距离误差模型的地球重力场恢复方法包括如下具体步骤：

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

（1）通过星载测距仪获取星间距离误差数据δρ₁₂；对于GRACE卫星采用星载K波段测距仪获取星间距离误差数据δρ₁₂，对于GRACE-II卫星采用星载激光干涉测距仪获取星间距离误差数据δρ₁₂。

（2）通过星载GPS接收机获取轨道位置数据r。

步骤2：星间距离误差模型建立

地球扰动位T(r,φ，λ)按球谐函数展开表示如下

$T(r,\mathrm{\φ},\mathrm{\λ})=\frac{\mathrm{GM}}{r}\underset{l=2}{\overset{L}{\mathrm{\Σ}}}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}\left[{\left(\frac{R_e}{r}\right)}^l({\stackrel{\‾}{C}}_{\mathrm{lm}}\cos \mathrm{m\λ}+{\stackrel{\‾}{S}}_{\mathrm{lm}}\sin \mathrm{m\λ}){\stackrel{\‾}{P}}_{\mathrm{lm}}\left(\sin \mathrm{\φ}\right)\right]---\left(1\right)$

其中，r,φ和λ分别表示卫星轨道的地心半径、地心纬度和地心经度，R_e表示地球的平均半径，GM表示地球质量M和万有引力常数G的乘积，L表示地球引力位按球谐函数展开的最大阶数，表示l阶和m次的缔合勒让德函数，表示待估的正规化地球引力位系数。

T(r,φ，λ)的功率谱表示如下

$P_l^2\left[T(r,\mathrm{\φ},\mathrm{\λ})\right]=\underset{m=0}{\overset{l}{\mathrm{\Σ}}}{[\frac{1}{4\mathrm{\π}}\∫\∫T(r,\mathrm{\φ},\mathrm{\λ}){\stackrel{\‾}{Y}}_{\mathrm{lm}}(\mathrm{\φ},\mathrm{\λ})\cos \mathrm{\φd\φd\λ}]}^2---\left(2\right)$

其中， ${\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}(\mathrm{\&phi;},\mathrm{\&lambda;})={\stackrel{\&OverBar;}{P}}_{l\left|m\right|}\left(\sin \mathrm{\&phi;}\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..$

基于球谐函数的正交归一性，公式（2）可被简化为

$P_l^2\left[T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\right]={\left(\frac{\mathrm{GM}}{R_e}\right)}^2{\left(\frac{R_e}{r}\right)}^{2l+2}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}^2+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}^2)---\left(3\right)$

大地水准面高的功率谱表示如下

$P_l^2\left[N\right]=R_e^2\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}^2+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}^2)---\left(4\right)$

联合公式（3）和（4），和的关系式表示如下

$P_l^2\left[N\right]=R_e^2{\left(\frac{R_e}{\mathrm{GM}}\right)}^2{\left(\frac{r}{R_e}\right)}^{2l+2}{P}_l^2\left[T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\right]---\left(5\right)$

在球坐标系中，T(r,φ，λ)对φ和λ的偏微分表示如下

$\left\{\begin{array}{c}\frac{\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&PartialD;\mathrm{\&phi;}}=\frac{\mathrm{GM}}{r}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}\left[{\left(\frac{R_e}{r}\right)}^l({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;})({\stackrel{\&OverBar;}{P}}_{l,m+1}\left(\sin \mathrm{\&phi;}\right)-\mathrm{mtg\&phi;}{\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\sin \mathrm{\&phi;}\right))\right]\\ \frac{\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&PartialD;\mathrm{\&lambda;}}=\frac{\mathrm{GM}}{r}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}\left[{\left(\frac{R_e}{r}\right)}^l(-m{\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}\sin \mathrm{m\&lambda;}+m{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}\cos \mathrm{m\&lambda;}){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\sin \mathrm{\&phi;}\right)\right]\end{array}\right.---\left(6\right)$

如图2所示，三角线、圆圈线和十字线分别表示 $P_l^2[\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})/\&PartialD;\mathrm{\&phi;}]$ 和 $P_l^2[\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})/\&PartialD;\mathrm{\&lambda;}].$ $P_l^2[\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})/\&PartialD;\mathrm{\&lambda;}]$ 和 $P_l^2\left[T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\right]$ 的关系式表示如下

$P_l^2[\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})/\&PartialD;\mathrm{\&lambda;}]={\left(\frac{\mathrm{GM}}{R_e}\right)}^2{\left(\frac{R_e}{r}\right)}^{2l+2}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}{m}^2({\stackrel{\&OverBar;}{C}}_{\mathrm{lm}}^2+{\stackrel{\&OverBar;}{S}}_{\mathrm{lm}}^2)=\frac{l^2}{2}{P}_l^2\left[T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\right]---\left(7\right)$

基于球对称性， $P_l^2[\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})/\&PartialD;\mathrm{\&phi;}]$ 和 $P_l^2[\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})/\&PartialD;\mathrm{\&lambda;}]$ 相等

$P_l^2\left[\frac{\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&PartialD;\mathrm{\&phi;}}\right]=P_l^2\left[\frac{\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&PartialD;\mathrm{\&lambda;}}\right]=\frac{l^2}{2}{P}_l^2\left[T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\right]---\left(8\right)$

基于能量守恒法，单星观测方程表示如下

$\frac{1}{2}{\stackrel{\&CenterDot;}{r}}^2=V_0+T+C---\left(9\right)$

其中，表示速度，V₀表示中心引力位，C表示能量常数。

双星观测方程表示如下

$\frac{1}{2}({\stackrel{\&CenterDot;}{r}}_2+{\stackrel{\&CenterDot;}{r}}_1){\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_{12}=T_2-T_1---\left(10\right)$

其中，和表示卫星的绝对速度，表示星间速度，T₁和T₂表示双星的地球扰动位。

在公式（10）两边同时乘以采样间隔Δt可得

$\frac{1}{2}({\stackrel{\&CenterDot;}{r}}_2+{\stackrel{\&CenterDot;}{r}}_1){\mathrm{\&rho;}}_{12}=(T_2-T_1)\mathrm{\&Delta;t}---\left(11\right)$

其中，表示卫星的平均速度；ρ₁₂=r₁₂·e₁₂表示星间距离，r₁₂=r₂-r₁表示双星的相对位置，e₁₂=r₁₂/|r₁₂|表示由第一颗卫星指向第二颗卫星的单位矢量；表示地球扰动位差分，

星间距离ρ₁₂的功率谱表示如下

$P_l^2\left[{\mathrm{\&rho;}}_{12}\right]=\frac{r{\left(\mathrm{\&Delta;t}\right)}^2}{\mathrm{GM}}{P}_l^2\left[\frac{\&PartialD;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&PartialD;\mathrm{\&phi;}}\right]{\left(\mathrm{\&Delta;\&phi;}\right)}^2---\left(12\right)$

联合公式（5）、（8）和（12），累计大地水准面精度和星间距离误差δρ₁₂之间的关系式表示如下

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

步骤3：地球重力场恢复

基于星间距离误差模型法，通过星间距离误差数据δρ₁₂，恢复GRACE和GRACE-II累计大地水准面精度的过程如下：

第一步，利用9阶Runge-Kutta线性单步法结合12阶Adams-Cowell线性多步法数值积分公式模拟了GRACE-II双星的星历。模拟轨道如图3所示，模拟过程耗时约2小时。

第二步，以0.5°×0.5°为网格分辨率，在地球表面的经度λ（0°～360°）和纬度φ（-90°~90°）范围内绘制网格，按照卫星轨道在地球表面的轨迹点位置依次加入星间距离误差δρ₁₂(φ，λ)。本发明利用不同的参考球面网格分辨率0.1°×0.1°～10°×10°分别恢复了地球重力场精度。结果表明：随着网格分辨率的增加，虽然网格化误差逐渐减小，但是计算耗时却大幅度提高。权衡利弊，我们选择网格分辨率0.5°×0.5°，在保证地球重力场恢复精度的前提下，可有效提高计算速度。

第三步，将星间距离误差δρ₁₂(φ，λ)按球谐函数展开为

${\mathrm{\&delta;\&rho;}}_{12}(\mathrm{\&phi;},\mathrm{\&lambda;})=\underset{l=0}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}\left[(C_{{\mathrm{\&delta;\&rho;}}_{\mathrm{lm}}}\cos \mathrm{m\&lambda;}+S_{{\mathrm{\&delta;\&rho;}}_{\mathrm{lm}}}\sin \mathrm{m\&lambda;}){\stackrel{\&OverBar;}{P}}_{\mathrm{lm}}\left(\sin \mathrm{\&phi;}\right)\right]---\left(14\right)$

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

$(C_{{\mathrm{\&delta;\&rho;}}_{\mathrm{lm}}},S_{{\mathrm{\&delta;\&rho;}}_{\mathrm{lm}}})=\frac{1}{4\mathrm{\&pi;}}\&Integral;\&Integral;\left[{\mathrm{\&delta;\&rho;}}_{12}(\mathrm{\&phi;},\mathrm{\&lambda;}){\stackrel{\&OverBar;}{Y}}_{\mathrm{lm}}(\mathrm{\&phi;},\mathrm{\&lambda;})\cos \mathrm{\&phi;d\&phi;d\&lambda;}\right]---\left(15\right)$

δρ₁₂的方差表示如下

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

将公式（16）代入（13），基于δρ₁₂可精确和快速地确定全球重力场的精度。

图4表示基于星间距离误差模型法分别利用GRACE以及相同星间距离50km和不同轨道高度的GRACE-II卫星恢复累计大地水准面精度对比；其中对于GRACE卫星采用实测数据进行处理，并将处理结果与国际公布的实测结果进行对比依此来验证所建立的星间距离误差模型的有效性和准确性；对于GRACE-II卫星则采用数值模拟数据来估计本发明应用于GRACE-II卫星的实际效果。

对GRACE卫星实测数据处理的结果与国际公布的实测结果对比如图4所示，十字线表示德国波兹坦地学研究中心（GFZ）公布的120阶EIGEN-GRACE02S全球重力场模型的实测精度，在120阶处，累计大地水准面精度为1.839×10^-1m；虚细线表示本发明基于2009年美国宇航局喷气推进实验室（NASA-JPL）公布的星载K波段测距仪的星间距离实测误差数据（测量精度10^-5m）恢复累计大地水准面的精度，在120阶处，累计大地水准面精度为1.826×10^-1m。通过两条曲线在各阶处的符合性可知，本发明建立的基于星间距离误差模型的地球重力场恢复方法是可靠的。

对于检验之后的星间距离误差模型，采用数值模拟数据对于GRACE-II卫星的应用效果进行估算的结果如图4所示，生成时间长度30天和采样间隔10秒的星间距离正态分布随机白噪声来代替需要GRACE-II卫星实测的星间距离误差数据δρ₁₂，实粗线、虚粗线和实细线分别表示基于激光干涉测距仪星间距离误差模型法（星间距离测量精度10^-8m），利用卫星轨道高度250km、350km和450km恢复GRACE-II累计大地水准面的精度。在360阶处，当卫星轨道高度选择为250km，累计大地水准面的误差为5.263×10^-2m；当卫星轨道高度选择为350km，累计大地水准面的误差提高了189倍；当卫星轨道高度选择为450km，累计大地水准面的误差提高了35622倍。结果表明：随着卫星轨道高度逐渐增加（250－450km），全球重力场的精度迅速降低。因此，GRACE-II（~250km）全球重力场的精度较GRACE（~450km）至少高一个数量级的最主要原因是较大程度地降低了GRACE-II卫星的轨道高度，从而全球重力场信号随卫星轨道高度增加的衰减效应得到了有效抑制。

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

Claims (6)

1.一种基于星间距离误差模型的地球重力场恢复方法，其特征在于包括如下步骤：

步骤1：采集卫星关键载荷数据，通过星载测距仪获取星间距离误差数据δρ₁₂，通过星载GPS接收机获取轨道位置数据r；

步骤2：通过星间距离误差数据δρ₁₂与累计大地水准面精度的关系，建立星间距离误差模型；

步骤3：基于所述星间距离误差模型，通过所采集的卫星关键载荷数据，对地球重力场进行恢复；其中，所述步骤3包括：

步骤3.1：利用9阶Runge-Kutta线性单步法结合12阶Adams-Cowell线性多步法数值积分公式模拟卫星的星历；

步骤3.2：确定参考球面网格分辨率，在地球表面的经度λ和纬度φ范围内按所确定的参考球面网格分辨率绘制网格，按照卫星轨道在地球表面的轨迹点位置依次加入星间距离误差δρ₁₂(φ,λ)；

步骤3.3：基于所述星间距离误差模型和星间距离误差数据δρ₁₂恢复地球重力场；

所述步骤2为：

将地球扰动位T(r,φ,λ)按球谐函数展开表示为

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

其中，r,φ和λ分别表示卫星轨道的地心半径、地心纬度和地心经度，R_e表示地球的平均半径，GM表示地球质量M和万有引力常数G的乘积，L表示地球引力位按球谐函数展开的最大阶数，表示l阶和m次的缔合勒让德函数，表示待估的正规化地球引力位系数；

T(r,φ,λ)的功率谱表示为

$P_l^2\&lsqb;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\&rsqb;=\underset{m=0}{\overset{l}{\&Sigma;}}{\&lsqb;\frac{1}{4\mathrm{\&pi;}}\&Integral;\&Integral;T(r,\mathrm{\&phi;},\mathrm{\&lambda;}){\stackrel{\&OverBar;}{Y}}_{lm}(\mathrm{\&phi;},\mathrm{\&lambda;})cos\mathrm{\&phi;}d\mathrm{\&phi;}d\mathrm{\&lambda;}\&rsqb;}^2---\left(2\right)$

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

基于球谐函数的正交归一性，公式（2）被简化为

${P_l}^2\&lsqb;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\&rsqb;={\left(\frac{GM}{R_e}\right)}^2{\left(\frac{R_e}{r}\right)}^{2l+2}\underset{m=0}{\overset{l}{\&Sigma;}}({\stackrel{\&OverBar;}{C}}_{lm}^2+{\stackrel{\&OverBar;}{S}}_{lm}^2)---\left(3\right)$

大地水准面高的功率谱表示为

${P_l}^2\&lsqb;N\&rsqb;=R_e^2\underset{m=0}{\overset{l}{\&Sigma;}}({\stackrel{\&OverBar;}{C}}_{lm}^2+{\stackrel{\&OverBar;}{S}}_{lm}^2)---\left(4\right)$

联合公式（3）和公式（4），P_l ²[N]和P_l ²[T(r,φ,λ)]的关系式表示为

${P_l}^2\&lsqb;N\&rsqb;=R_e^2{\left(\frac{R_e}{GM}\right)}^2{\left(\frac{r}{R_e}\right)}^{2l+2}{P_l}^2\&lsqb;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\&rsqb;---\left(5\right)$

在球坐标系中，T(r,φ,λ)对φ和λ的偏微分表示为

$\left\{\begin{array}{c}\frac{\&part;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&part;\mathrm{\&phi;}}=\frac{GM}{r}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}\&lsqb;{\left(\frac{R_e}{r}\right)}^l({\stackrel{\&OverBar;}{C}}_{lm}\cos m\mathrm{\&lambda;}+{\stackrel{\&OverBar;}{S}}_{lm}\sin m\mathrm{\&lambda;})\left({\stackrel{\&OverBar;}{P}}_{l,m+1}\right(sim\mathrm{\&phi;})-mtg\mathrm{\&phi;}{\stackrel{\&OverBar;}{P}}_{lm}\left(\sin \mathrm{\&phi;}\right))\&rsqb;\\ \frac{\&part;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&part;\mathrm{\&lambda;}}=\frac{GM}{r}\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=0}{\overset{l}{\mathrm{\&Sigma;}}}\&lsqb;{\left(\frac{R_e}{r}\right)}^l(-m{\stackrel{\&OverBar;}{C}}_{lm}\sin m\mathrm{\&lambda;}+m{\stackrel{\&OverBar;}{S}}_{lm}\cos m\mathrm{\&lambda;}){\stackrel{\&OverBar;}{P}}_{lm}\left(sim\mathrm{\&phi;}\right)\&rsqb;\end{array}\right.---\left(6\right)$

和P_l ²[T(r，φ，λ)]的关系式表示为

${P_l}^2\&lsqb;\&part;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})/\&part;\mathrm{\&lambda;}\&rsqb;={\left(\frac{GM}{R_e}\right)}^2{\left(\frac{R_e}{r}\right)}^{2l+2}\underset{m=0}{\overset{l}{\&Sigma;}}({\stackrel{\&OverBar;}{C}}_{lm}^2+{\stackrel{\&OverBar;}{S}}_{lm}^2)=\frac{l_2}{2}{P_l}^2\&lsqb;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\&rsqb;---\left(7\right)$

基于球对称性，和相等

${P_l}^2\&lsqb;\frac{\&part;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&part;\mathrm{\&phi;}}\&rsqb;=P^2\&lsqb;\frac{\&part;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&part;\mathrm{\&lambda;}}\&rsqb;=\frac{l^2}{2}{P_l}^2\&lsqb;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})\&rsqb;---\left(8\right)$

基于能量守恒法，单星观测方程表示为

$\frac{1}{2}{\stackrel{\&CenterDot;}{r}}^2=V_0+T+C---\left(9\right)$

其中，表示速度，V₀表示中心引力位，C表示能量常数；

双星观测方程表示为

$\frac{1}{2}({\stackrel{\&CenterDot;}{r}}_2+{\stackrel{\&CenterDot;}{r}}_1){\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_{12}=T_2-T_1---\left(10\right)$

其中，和表示卫星的绝对速度，表示星间速度，T₁和T₂表示双星的地球扰动位；

在公式（10）两边同时乘以采样间隔Δt可得

$\frac{1}{2}({\stackrel{\&CenterDot;}{r}}_2+\stackrel{\&CenterDot;}{r}){\mathrm{\&rho;}}_{12}=(T_2-T_1)\mathrm{\&Delta;}t---\left(11\right)$

其中，表示卫星的平均速度；ρ₁₂=r₁₂·e₁₂表示星间距离，r₁₂=r₂-r₁表示双星的相对位置，e₁₂=r₁₂/|r₁₂|表示由第一颗卫星指向第二颗卫星的单位矢量；表示地球扰动位差分，

星间距离ρ₁₂的功率谱表示为

${P_l}^2\&lsqb;{\mathrm{\&rho;}}_{12}\&rsqb;=\frac{r{\left(\mathrm{\&Delta;}t\right)}^2}{GM}{P_l}^2\&lsqb;\frac{\&part;T(r,\mathrm{\&phi;},\mathrm{\&lambda;})}{\&part;\mathrm{\&phi;}}\&rsqb;{\left(\mathrm{\&Delta;}\mathrm{\&phi;}\right)}^2---\left(12\right)$

联合公式（5）、（8）和（12），累计大地水准面精度和星间距离误差δρ₁₂之间的关系式表示为

${\mathrm{\&delta;N}}_{{\mathrm{\&rho;}}_{12}}=\sqrt{\frac{R_e^3{r}^2{\left(\mathrm{\&Delta;}t\right)}^2}{{\mathrm{GM\&rho;}}_{12}^2}\underset{l=2}{\overset{L}{\&Sigma;}}\&lsqb;\frac{2}{l^2}{\left(\frac{r}{R_e}\right)}^{2l+1}{\mathrm{\&sigma;}}_l^2\left({\mathrm{\&delta;\&rho;}}_{12}\right)\&rsqb;}---\left(13\right);$

所述步骤3.3为：

基于星间距离误差模型，通过星间距离误差数据δρ₁₂，恢复地球重力场的过程如下：

将星间距离误差δρ₁₂(φ,λ)按球谐函数展开为

${\mathrm{\&delta;\&rho;}}_{12}(\mathrm{\&phi;},\mathrm{\&lambda;})=\underset{l=0}{\overset{L}{\&Sigma;}}\underset{m=0}{\overset{l}{\&Sigma;}}\&lsqb;(C_{{\mathrm{\&delta;\&rho;}}_{lm}}\cos m\mathrm{\&lambda;}+S_{{\mathrm{\&delta;\&rho;}}_{lm}}\sin m\mathrm{\&lambda;}){\stackrel{\&OverBar;}{P}}_{lm}\left(sin\mathrm{\&phi;}\right)\&rsqb;---\left(14\right)$

其中，表示δρ₁₂(φ,λ)按球函数展开的系数

$(C_{{\mathrm{\&delta;\&rho;}}_{lm}},S_{{\mathrm{\&delta;\&rho;}}_{lm}})=\frac{1}{4\mathrm{\&pi;}}\&Integral;\&Integral;\&lsqb;{\mathrm{\&delta;\&rho;}}_{12}(\mathrm{\&phi;},\mathrm{\&lambda;}){\stackrel{\&OverBar;}{Y}}_{lm}(\mathrm{\&phi;},\mathrm{\&lambda;})cos\mathrm{\&phi;}d\mathrm{\&phi;}d\mathrm{\&lambda;}\&rsqb;---\left(15\right)$

δρ₁₂的方差表示为

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

将公式（16）代入公式（13），即可基于星间距离误差δρ₁₂恢复地球重力场。

2.如权利要求1所述的基于星间距离误差模型的地球重力场恢复方法，其特征在于：所述步骤3.2中所确定参考球面网格分辨率为0.1°×0.1°～10°×10°。

3.如权利要求2所述的基于星间距离误差模型的地球重力场恢复方法，其特征在于：所述步骤3.2中所确定参考球面网格分辨率为0.5°×0.5°。

4.如权利要求1-3中任一项所述的基于星间距离误差模型的地球重力场恢复方法，其特征在于：所述卫星为GRACE卫星或者下一代GRACE-II卫星。

5.如权利要求4所述的基于星间距离误差模型的地球重力场恢复方法，其特征在于：当所述卫星为GRACE卫星时，采用星载K波段测距仪获取星间距离误差数据δρ₁₂；当所述卫星为下一代GRACE-II卫星时，采用星载激光干涉测距仪获取星间距离误差数据δρ₁₂。

6.如权利要求5所述的基于星间距离误差模型的地球重力场恢复方法，其特征在于：当所述卫星为下一代GRACE-II卫星时，所述步骤2还包括：

使用所述星间距离误差模型并利用GRACE卫星所采集的实测数据进行地球重力场恢复；

将地球重力场恢复结果与EIGEN-GRACE02S地球重力场模型的实测精度进行对比；

根据对比结果对所述星间距离误差模型进行检验。