高精度的地球重力场模型对经济、军事以及地球科学研究等具有十分重要的作用。为此，国外发射了多颗具有重力场探测效能的地球观测卫星，如CHAMP、GRACE、GOCE以及海洋测高卫星等。其中，CHAMP、GRACE和GOCE统称为重力卫星，GRACE卫星相对于CHAMP和GOCE卫星有更加广泛的应用前景，如海洋、南极冰盖、陆地水等地球科学方面的应用(许厚泽 等，2012)。GRACE卫星已于2017年10月坠落，其后续卫星GRACE Follow-on已于2018年5月发射，近期Earth system mass transport mission(e.motion)和Next Generation satellite Gravimetry Mission(NGGM)卫星计划正在设计研究阶段(Panet 等，2013；Roland和Thomas，2015)。该类卫星的基本原理是通过高低跟踪和低低跟踪距离测量来探测地球重力场所引起的卫星星间距离的变化，并以此来反演地球重力场，其中高低跟踪测量主要通过星载GPS接收机来完成，低低跟踪则主要依靠微波测距完成。对于GRACE卫星而言，卫星低低跟踪数据的精度对重力场模型的精度起着至关重要的作用，因此在新的卫星计划中，除了微波测距外，还将搭载高精度的激光测距来进行更高精度的测距实验。由此，便产生这样的问题：激光测距对重力场模型恢复精度的提升会带来多大帮助，特别是当其他载荷的精度没有提升时。

事实上，关于以上问题，当前已有不少研究和结论。Loomis等人(2012)通过模拟仿真讨论了块质量异常体所引起的重力场变化的探测精度问题，其结论表明：激光测距对块质量异常体的探测精度的提升帮助有限，只有背景场的时变误差得到削弱后，才会明显改善该文所提出的质量异常体的探测精度。显然该文主要讨论的对象为时变重力场。Elsaka等人(2014)讨论了7种不同轨道构型下的重力场模型探测精度，其中一种构型(该文中以GRACE Follow-on表示)用激光来进行测距，其余载荷和GRACE一样，另外的区别是将轨道降低了40 km，此种情况下的重力场恢复精度和GRACE相比精度提高了约10倍，由于该结论受到了轨道降低的影响，因此难以估计激光单独的贡献占多大量级。此外多篇国外文献讨论了采用激光后的重力卫星计划，如Panet等人(2013)、Gruber等人(2014)和Sheard等人(2012)，但少有文献单独评估激光的贡献。

围绕下一代的低低跟踪重力卫星，国内多家单位进行了相应的指标论证研究。庞振兴等人(2012)利用谱分析方法估计了GRACE Follow-On的空间分辨率，该文指出后者比前者能够以更高的空间分辨率和至少高一个量级的精度来恢复地球的时变重力场和静态重力场，后者所采用的技术除了激光测距之外，还采用了无拖曳技术，并且提高了定轨、定速及加速度计的精度，因此此文中重力场精度的提升是各类因素的一个综合效应，而不是激光测距单独的贡献；郑伟等人(2010)利用解析法估计了GRACE Follow-On的地球重力场精度，但该文将定轨、定速误差与星间测距误差及加速度计误差作了等权处理，这与实际的重力场解算不符，很显然定轨定速数据的权重应低于星间测距和加速度计数据的权重；赵倩等人(2014)利用仿真分析讨论了观测误差对新一代重力卫星任务的影响，指出定轨精度主要影响重力场模型低阶项的反演精度，星间距变率主要影响中高阶。

综合以上的文献，发现已有文献很少单独评估激光测距精度本身的贡献，很多结论是各类新技术的综合效应，同时不同文献讨论的对象不同，有的关注的是时变重力场，而有的结论则主要针对的是静态重力场。时变重力场强调实时性，一般需要以较短的时间为单位来解算重力场，如以月为单位，这就对各类时变背景场模型的精度提出了更高的要求。例如：当前所有的GRACE时变重力场模型均存在南北条带误差问题，而静态重力场则不存在，这表明静态重力场由于能够使用更多的数据，因此能够降低对各类背景场模型的精度要求。另外一点也需要注意的是定轨定速数据的权重问题，本文认为二者与星间距变率数据和加速度计数据的权重应该有所区别，对于不同的解算方法也会有差别。针对以上问题，本文限制了相应的前提，即：(1)主要讨论地球静态重力场模型的精度，而不考虑重力场解算过程中的时变误差；(2)主要评估采用激光测距后，相对已有的GRACE卫星而言，其余载荷和轨道均不变的情况，重力场模型的精度提升有多大。在以上约束条件下，本文主要通过高精度的数值反演方法讨论本文开头所提出的问题。

GRACE卫星反演重力场的方法主要有Kaula线性摄动法、动力学法、基线法、能量守恒法、加速度法和短弧长法等。从理论上讲动力学方法来自于对卫星轨道方程组的线性化，理论上相对是严密的，精度相对其他方法更高，因此本研究采用动力学方法进行模拟计算。

动力学积分法(Tapley，1989)的基本原理是基于牛顿运动方程，将卫星轨道数据线性化表达成状态转移矩阵(即任意时刻的状态向量(卫星或速度)对初始时刻状态向量的偏导数矩阵)和参数敏感矩阵(即任意时刻的状态向量对未知参数的偏导数，如引力场位系数、加速度校正参数和力模型未知参数)的形式，来建立卫星跟踪数据与重力场位系数间的关系式，然后基于最小二乘准则解算引力场位系数。动力学定轨的微分方程可用以下形式表示：

$\left\{ \begin{gathered} \frac{{{\rm{d}}{{r}}\left(t \right)}}{{{\rm{d}}t}} = {\dot{ r}}(t) \\ \frac{{{\rm{d}}{\dot{ r}}\left(t \right)}}{{{\rm{d}}t}} = {{f}}\left(t \right) \\ {{r}}\left({{t_0}} \right) = {{{r}}_0} \\ {\dot{ r}}\left({{t_0}} \right) = {{{\dot{ r}}}_0} \\ \end{gathered} \right.$

(1)

式中， ${{r}}\left(t \right)$ 、 ${\dot{ r}}\left(t \right)$ 分别为t时刻的位置和速度向量， ${{f}}\left(t \right)$ 为卫星所受到的全部摄动力向量， ${{{r}}_0}$ 、 ${{\dot{ r}}_{{0}}}$ 是初始时刻的位置和速度。若令 ${{X}}\left(t \right) = \left(\begin{gathered} {{r}}\left(t \right) \\ {\dot{ r}}\left(t \right) \\ \end{gathered} \right), {{F}}\left(t \right) =$ $ \left(\begin{gathered} {\dot{ r}}\left(t \right) \\ {{f}}\left(t \right) \\ \end{gathered} \right),{{X}}\left({{t_0}} \right) = \left(\begin{gathered} {{{r}}_0}\left(t \right) \\ {{{\dot{ r}}}_0}\left(t \right) \\ \end{gathered} \right)$ 则式(1)可写成：

${{X}}\left(t \right) = {{X}}\left({{t_0}} \right) + \int\limits_t {{{F}}\left(t \right)} {\rm{d}}t$

(2)

$\left\{ \begin{gathered} \frac{{{\rm{d}}{{X}}\left(t \right)}}{{{\rm{d}}t}} = {{F}}\left(t \right) \\ {{X}}\left({{t_0}} \right) = {{{X}}_0} \\ \end{gathered} \right.$

(3)

对式(2)进行积分，可以得到任意时刻的状态向量 ${{X}}\left(t \right)$ ，但是由于初始状态向量是未知数，而力模型 ${{F}}\left(t \right)$ 难以精确确定，包含未知参数，如加速度观测值的尺度和偏差参数，所以很难求解式(2)。因此我们通过将其在参考轨道进行Taylor展开，将其表达成参数敏感矩阵和状态转移矩阵的形式：

${{X}}\left(t \right) = {{{X}}_ {\rm{ref}}}\left(t \right) + {{S}}\left(t \right)\delta {{p}} + {{Q}}\left(t \right)\delta {{X}}\left({{t_0}} \right)$

(4)

其矩阵形式为

$\delta {{X}}\left(t \right) = \left({{{Q}}\left(t \right) {{S}}\left(t \right) } \right)\left(\begin{gathered} \delta {{X}}\left({{t_0}} \right) \\ \delta {{p}} \\ \end{gathered} \right)$

(5)

式中， ${{{X}}_ {\rm{ref}}}\left(t \right)$ 为利用先验力模型和相同的初始状态向量数值积分出参考轨道的状态向量， $\delta {{X}}\left(t \right)$ 为实际轨道与参考轨道之差。 ${{Q}}\left(t \right) = \left({\begin{array}{*{20}{c}} {{{\varPhi }}\left(t \right)} \\ {{\dot{ \varPhi }}\left(t \right)} \end{array}} \right)$ 为状态转移矩阵， ${{{\varPhi }}\left(t \right)}$ 与 ${{\dot{ \varPhi }}\left(t \right)} $ 分别表示对位置和速度的状态转移矩阵， $\delta {{X}}\left({{t_0}} \right)$ 为初始时刻的状态改正量。 ${{S}}\left(t \right) = \left({\begin{array}{*{20}{c}} {{{{B}}}\left(t \right)} \\ {{\dot{ { B}}}\left(t \right)} \end{array}} \right)$ 为参数敏感矩阵，B(t)与 ${{\dot{ {B}}}\left(t \right)} $ 分别表示对位置与速度的参数敏感矩阵， $\delta {{P}}$ 为力模型改正参数，一般包括位系数参数、加速度计观测数据的尺度和偏差参数和大气阻力参数等。参数敏感矩阵和状态转移矩阵很难得到解析解，通常采用数值积分求解其微分方程。通常的数值微分方法主要有龙格库塔、Adam预测校正和高斯-杰克森方法等。

GRACE卫星主要测量两颗卫星之间的速度在两颗卫星的水平线上的分速度，即星间变率。星间变率测量的精度高于GPS测量位置差分得到的速度。相比于卫星精密轨道数据，星间变率的优势是对地球重力场中高阶次较为敏感，并且能够减少SST测量系统振荡器的不稳定性。但星间变率与重力场位系数等参数之间是非线性的，本文在于锦海等人(2017)推导的轨道扰动方程组的形式上进行星间变率公式线性化的推导，现直接给出卫星轨道扰动方程推导得到的式(6)和式(7)，分别用A和B来表示GRACE的两颗卫星，方程如下：

${{{r}}^{\rm{A}}}(t) = {{r}}_ {\rm{ref}}^{\rm{A}}(t) + {{{\varPhi }}^{\rm{A}}}(t)\delta {{{X}}^{\rm{A}}}\left({{t_0}} \right) + {{{{B}}}^{\rm{A}}}\left(t \right)\delta {{p}}$

(6)

${{{r}}^{\rm{B}}}(t) = {{r}}_ {\rm{ref}}^{\rm{B}}(t) + {{{\varPhi }}^{\rm{B}}}(t)\delta {{{X}}^{\rm{B}}}\left({{t_0}} \right) + {{{{B}}}^{\rm{B}}}\left(t \right)\delta {{p}}$

(7)

其导数可以写为

${{\dot{ r}}^{\rm{A}}}(t) = {\dot{ r}}_ {\rm{ref}}^{\rm{A}}(t) + {{\dot{ \varPhi }}^{\rm{A}}}(t)\delta {{{X}}^{\rm{A}}}\left({{t_0}} \right) + {{\dot{ { B}}}^{\rm{A}}}\left(t \right)\delta {{p}}$

(8)

${{\dot{ r}}^{\rm{B}}}(t) = {\dot{ r}}_ {\rm{ref}}^{\rm{B}}(t) + {{\dot{ \varPhi }}^{\rm{B}}}(t)\delta {{{X}}^{\rm{B}}}\left({{t_0}} \right) + {{\dot{ {B}}}^{\rm{B}}}\left(t \right)\delta {{p}}$

(9)

式中， ${{r}}_ {\rm{ref}}^{\rm{A}}(t)$ 和 ${{r}}_ {\rm{ref}}^{\rm{B}}(t)$ 是数值计算得到的参考轨道，而 ${{{r}}^{\rm{A}}}(t)$ 和 ${{{r}}^{\rm{B}}}(t)$ 是GPS观测值。分别将式(6)减去式(7)，将式(8)减去式(9)可得星间距和星间距变率方程：

$ \begin{aligned}{{\rho }}(t) = & {{{\rho }}_0}(t) + {{{\varPhi }}^{\rm{A}}}\left(t \right)\delta {{{X}}^{\rm{A}}}\left({{t_0}} \right) - {{{\varPhi }}^{\rm{B}}}\left(t \right)\delta {{{X}}^{\rm{B}}}\left({{t_0}} \right) + \\ & \left({{{{B}}^{\rm{A}}}\left(t \right) - {{{B}}^{\rm{B}}}\left(t \right)} \right)\delta {{p}}\end{aligned}$

(10)

$\begin{aligned}{\dot{ \rho }}(t) = & {{\dot{ \rho }}_0}(t) + {{\dot{ \varPhi }}^{\rm{A}}}\left(t \right)\delta {{{X}}^{\rm{A}}}\left({{t_0}} \right) - {{\dot{ \varPhi }}^{\rm{B}}}\left(t \right)\delta {{{X}}^{\rm{B}}}\left({{t_0}} \right) + \\ & \left({{{{\dot{ {B}}}}^{\rm{A}}}\left(t \right) - {{{\dot{ {B}}}}^{\rm{B}}}\left(t \right)} \right)\delta {{p}}\end{aligned}$

(11)

式中， ${{\rho }}(t) = {{{r}}^{\rm{A}}}(t) - {{{r}}^{\rm{B}}}(t), {{\rho }}{(t)_0} = {{r}}_ {\rm{ref}}^{\rm{A}}(t) - {{r}}_ {\rm{ref}}^{\rm{B}}(t)$ ，在式(10)和式(11)中 ${{\rho }}_0$ 、 $\dot{{\rho }}_0$ 分别是 ${{\rho }}$ 、 $\dot{{\rho }}$ 的主项，而 $\delta {{{X}}^{\rm{A}}}\left( {{t_0}} \right)$ 、 $\delta {{{X}}^{\rm{B}}}\left( {{t_0}} \right)$ 在实际使用时其实质是GRACE卫星轨道的测量误差，因此 $\delta {{{X}}^{\rm{A}}}\left({{t_0}} \right), \delta {{{X}}^{\rm{B}}}\left({{t_0}} \right)$ 也是小量。由GRACE卫星的星间距和星间距变率的实际测量值可知，其量级大约为 $ {d _0} \approx d \approx 250 $ km， $ {\dot d _0} \approx \dot d \approx 1.0 $ m/s，这里 $ {d _0} = \left| {{{{\rho }}_0}} \right|$ ， $d = \left| {{\rho }} \right|$ 。对于GRACE卫星而言，其观测值为 $D(t) = {\dot{ \rho }} \cdot \displaystyle\frac{{{\rho }}}{d }$ 。为了得到观测方程，需要对 ${\dot{ \rho }} \cdot \displaystyle\frac{{{\rho }}}{d }$ 进行线性化展开。为使推导方便，引入记号：

$\begin{aligned}{{\varepsilon }}(t) = & {{{\varPhi }}^{\rm{A}}}\left(t \right)\delta {{{X}}^{\rm{A}}}\left({{t_0}} \right) - {{{\varPhi }}^{\rm{B}}}\left(t \right)\delta {{{X}}^{\rm{B}}}\left({{t_0}} \right) + \\ & \left({{{{B}}^{\rm{A}}}\left(t \right) - {{{ B}}^{\rm{B}}}\left(t \right)} \right)\delta {{p}}\end{aligned}$

(12)

$\begin{aligned}{\dot{{ \varepsilon }}}(t) = & {\dot{{ \varPhi }}^{\rm{A}}}\left(t \right)\delta {{{X}}^{\rm{A}}}\left({{t_0}} \right) - {{\dot{ \varPhi }}^{\rm{B}}}\left(t \right)\delta {{{X}}^{\rm{B}}}\left({{t_0}} \right) + \\ & \left({{{{\dot{ {B}}}}^{\rm{A}}}\left(t \right) - {{{\dot{ {B}}}}^{\rm{B}}}\left(t \right)} \right)\delta {{p}}\end{aligned}$

(13)

而 $ {{\varepsilon }}{{(t)}}$ 相对于 $ {{{\rho }}_0}$ 是小量，因此

${{\rho }}(t) = {{{\rho }}_0}(t) + {{\varepsilon }}(t)$

(14)

${\dot{ \rho }}(t) = {{\dot{ \rho }}_0}(t) + {\dot{ \varepsilon }}(t)$

(15)

由星间距矢量和星间距标量的关系及式(12)可以得

${d ^2} = {{{\rho }}^2}(t) = {[{{{\rho }}_0}(t) + {{\varepsilon }}(t)]^2} = d _0^2(t) + 2{{{\rho }}_0}(t) \cdot {{\varepsilon }}(t) + {\varepsilon ^2}$

(16)

根据式(11)可知 ${\varepsilon ^2}$ 是小量，所以略去 ${\varepsilon ^2}$ 量级后式(14)可以写为

$d (t) = {d _0}(t) + \frac{{{{{\rho }}_0}(t)}}{{{d _0}(t)}}{{\varepsilon }}(t)$

(17)

由式(17)变形可得

$\frac{1}{{d (t)}} = \frac{1}{{{d _0}(t)}}\left[ {1 - \frac{{{{{\rho }}_0}(t) \cdot {{\varepsilon }}(t)}}{{d _0^2(t)}}} \right]$

(18)

代回观测值 $D(t) = {\dot{ \rho }} \cdot \displaystyle\frac{{{\rho }}}{d }$ ，便有

$D(t) = \frac{{[{{{\dot{ \rho }}}_0}(t) + {\dot{ \varepsilon }}(t)] \cdot [{{{\rho }}_0}(t) + {{\varepsilon }}(t)]}}{{{d _0}(t)}}\left[ {1 - \frac{{{{{\rho }}_0}(t) \cdot {{\varepsilon }}(t)}}{{d _0^2(t)}}} \right]$

(19)

在式(18)中再进一步地舍去 ${\varepsilon ^2}$ 量级小量后，便得线性化结果如下：

$\begin{aligned} D(t) = & {D_0}(t) + \frac{{{{{\dot{ \rho }}}_0}(t) \cdot {{\varepsilon }}(t) + {{{\rho }}_0}(t) \cdot {\dot{ \varepsilon }}(t)}}{{{d _0}(t)}} - \\ & \frac{{[{{{\dot{ \rho }}}_0}(t) \cdot {{{\rho }}_0}(t)][{{{\rho }}_0}(t) \cdot {{\varepsilon }}(t)]}}{{d _0^3(t)}} {\text{ = }} \\ &{D_0}(t) + \left({\frac{{{{{\dot{ \rho }}}_0}(t)}}{{{d _0}(t)}} - \frac{{[{{{\dot{ \rho }}}_0}(t) \cdot {{{\rho }}_0}(t)]{{{\rho }}_0}(t)}}{{d _0^3(t)}}} \right) \cdot \\ & {{\varepsilon }}(t) + \frac{{{{{\dot{ \rho }}}_0}(t)}}{{{d _0}(t)}} \cdot {\dot{ \varepsilon }}(t) \end{aligned}$

(20)

这里 ${D_0}(t) = {{\dot{ \rho }}_0} \cdot \displaystyle\frac{{{{{\rho }}_0}}}{{{d _0}}}$ 为已知值(利用参考轨道计算)。将式(12)、式(13)代入式(19)便是关于待求参数的观测方程，因此式(20)仅是待求参数的线性方程。在观测式(20)中带有下标“0”的项都是从参考轨道计算得到已知量。

由于两颗卫星在同一轨道上一前一后分别运行，若单独使用星间变率数据，则势必会造成局部法方程的欠定问题。因此必须要将卫星精密轨道数据与星间变率数据一起参与计算才可解算正确模型。一方面有卫星精密轨道数据提供卫星位置信息来降低局部法方程欠定的问题。另一方面轨道观测值对卫星重力场的低阶项较为敏感，能够提供一定精度的重力场信息。由于轨道的相对误差比较大，因此会将误差传递到法方程中。因此必须加相应的权比降低其对方程的影响，综上所述，必须融合卫星精密轨道数据及星间变率观测值来反演GRACE卫星的重力场，而实质是通过一定的比例关系将精密轨道方程与星间变率方程相加后进行求解，因此可以直接建立起两个法方程之间的关系。本文直接给出星间变率和精密轨道最终的法方程可以写为式(21)和式(22)：

${{{N}}_ {\rm{kbrr}}}{{x}} = {{{L}}_ {\rm{kbrr}}}$

(21)

${{{N}}_ {\rm{orbit}}}{{x}} = {{{L}}_ {\rm{orbit}}}$

(22)

式中， ${{x}}$ 代表两颗卫星的初始状态参数、重力场位系数和加速度计标定参数等，下标“kbrr”代表与星间变率相关的量，下标“orbit”代表与卫星精密轨道相关的量。 ${\delta _{{\text{orbit}}}}$ 是轨道观测精度， ${\delta _ {\rm{kbrr}}}$ 是星间变率观测精度。用 $\alpha =\displaystyle \frac{{\delta _ {\rm{orbit}}^2}}{{\delta _ {\rm{kbrr}}^2}}$ 来表示两者系数之比，将式(21)乘以α与式(22)相加可以得到融合星间变率和精密轨道数据的法方程，即：

$\left({\alpha {N_ {\rm{kbrr}}} + {N_ {\rm{orbit}}}} \right)x = \alpha {L_ {\rm{kbrr}}} + {L_ {\rm{orbit}}}$

(23)

由于星间变率测量精度要高于GPS测量轨道的精度，如果选取不恰当的权，则高精度的重力场的计算精度一定会被污染或者由于法方程条件数过高而使方程欠定，而权比的确定实质上是由卫星的观测精度决定。本文不加推导的利用王长青(2015)提供的轨道与星间变率的权比 $\alpha {\text{ = }}{10^{10}}$ 进行计算。

通过数值实验来讨论论文开头提出的问题。(1)首先来验证有关程序的正确性；(2)讨论不同情况下反演精度。

(1)先用数值积分方法模拟出GRACE 1个月的轨道数据，模拟GRACE的轨道参数与真实轨道参数相同，然后并进一步得到星间距变率数据，再利用上节所述的方法进行反演计算。为了减小计算量，模拟轨道使用EGM08模型前60阶完整阶次来当作真实场，数值积分方法我们使用Adam预测校正方法，所使用的软件为Intel FORTRAN 64位软件。在解算时不加入任何误差，得到了一组重力场参数，将真实的重力场与反演计算得到的重力场相减，即式(24)画出了KUALA曲线误差图(Kaula，1966)，误差量级越小代表计算的精度越高，结果如图1所示。

${\varsigma _n} = \sqrt {\frac{1}{{2n + 1}}\sum\limits_{m = 0}^n {\left[ {{{\left({{C_{nm}} - C_{nm}^{\rm{ref}}} \right)}^2} + {{\left({{S_{nm}} - S_{nm}^{\rm{ref}}} \right)}^2}} \right]} } $

(24)

式中，n代表重力场计算阶数， ${\varsigma _n}$ 代表第n阶的KUALA误差值， ${C_{nm}}$ 、 ${S_{nm}}$ 、 $C_{nm}^{\rm{ref}}$ 和 $S_{nm}^{\rm{ref}}$ 分别表示计算得到的位系数和参考场的位系数。

图1中，“EGM08”代表真实重力场EGM08模型的KUALA曲线，“EGM08-EGM96”代表真实重力场EGM08模型减去参考模型EGM96以后得到的曲线，即不同重力场模型之间的误差曲线，“Error-free”代表利用动力学方法反演得到的重力场与模型EGM08相减后的误差曲线，该曲线表明了本实验计算程序的正确性，能够满足本研究的计算精度要求。

(2)解算流程与计算方法与(1)一致。区别在于首先将重力场的阶数增加到120阶。然后对观测数据加入几组误差，误差如表1所示，即GRACE卫星的设计精度指标，以及GRACE Follow On 的测量精度指标，主要是提高了星间变率的测量精度，结果如图2所示。

图2中“EGM08”代表重力场模型EGM08的KUALA曲线，“GRACE”代表GRACE卫星的设计精度指标，“Follow-1、Follow-2、Follow-3、Follow-4”分别代表星间变率计算精度提高到5.0×10^–7 m/s、1.0×10^–7 m/s、5.0×10^–8 m/s、1.0×10^–8 m/s。该计算结果表明当星间距变率提高到1.0×10^–8 m/s时，相对于GRACE现有的精度指标，重力场恢复精度能够提高1个量级。

本文将图2中的误差曲线图转换为图3的累计大地水准面误差图，各曲线的表示与图2一致。通过分析图3可知：保持GRACE现有的轨道和其他载荷精度，当星间变率精度依次从1.0×10^–6 m/s提高到5.0×10^–7 m/s、1.0×10^–7 m/s、5.0×10^–8 m/s、1.0×10^–8 m/s时，累积大地水准面误差则在120阶时从85.14 cm依次降低为33.09 cm、7.33 cm、3.70 cm、3.59 cm。这表明当采用高精度的激光测距后，采用低低卫卫跟踪模式，地球静态重力场模型的探测精度有望比GRACE提高1个量级。上述结论与冉将军等人(2012，2015)的结论一致。但星间变率精度为5.0×10^–8 m/s和1.0×10^–8 m/s大地水准面累积误差和KUALA曲线图几乎无差别。

本文主要讨论了基于轨道扰动方程对星间变率进行了线性化的推导，并通过一定的权值将线性化的星间变率方程组与精密轨道方程组融合，通过实例验证了公式及程序的正确性。通过仿真计算验证了激光测距对静态重力场恢复精度提升的贡献。结果表明：当激光测距的精度达到50 nm/s时，轨道和其余指标与GRACE卫星相同时，静态重力场的恢复精度可提高约1个量级。但当激光测距精度提高至10 nm/s时，计算结果与精度为50 nm/s的计算结果无明显差别，这表明过高的星间变率测量精度相对于其他指标而言精度有冗余。因此，建议在其他测量技术精度未有提高的前提下，只需将激光测距的精度提高至50 nm/s即可。

值得说明的是：以上工作所讨论的对象为地球静态重力场，而非时变重力场。静态重力场反演可利用多年的数据，如1年甚至5年，因此部分时变误差可被大大削弱和平滑。而对于时变重力场，由于其解算以月甚至更短时间为单位，对时变环境误差的扣除提出了更高的要求(例如：固体潮、海潮、大气潮等)，因此仅仅提高星间测距，对于时变重力场的解算不如静态重力场的帮助明显。时变重力场的解算除了提高测距精度外，还需发射多颗不同轨道倾角的卫星来提供观测数据。而无论是对于静态重力场，还是时变重力场的解算，通过激光测距来提高星间距离测量精度，是提高下一代重力卫星反演精度的关键因素之一。因此，及时开展基于激光测距技术的下一代重力卫星关键技术研究具有重要意义。

本文未能对时变重力场的精度提升问题进行系统的分析，下一步可对重力卫星测量的其他精度指标以及背景场模型(如海潮、大气潮等)误差对时变重力场恢复精度的影响进行研究，从而为大幅提升时变重力场的精度提供可行的卫星设计方案。