GRACE(gravity recovery and climate experiment)卫星观测数据在大地测量学等多个地球科学取得了广泛应用^[1-4]，利用GRACE数据反演高精度、高分辨率的地球重力场模型是这些应用的基础。目前，国际全球重力场模型中心ICGEM(International Centre for Global Earth Models)已经发布的数十套GRACE时变重力场模型主要由动力学法解算^[5-11]，部分采用短弧边值法^[12-14]和加速度法^[15-16]解算，我国学者利用这些方法也解算了多种静态与时变重力场模型^[17-21]。这3种解算方法都以牛顿运动定律和万有引力定律为基础，在理论上是等价的。本质上都是根据卫星位置、星间距离和速度观测量(或其加速度和星间相对加速度导出量)相对于卫星参考轨道(或其星间距离、速度和加速度导出量)的线性摄动量建立线性观测方程进行求解。动力学法以轨道起始点的状态向量(3个坐标和3个速度分量)为初值，利用先验力模型积分计算出整个弧段的参考轨道；积分外推导致线性化误差随轨道弧段的增长而迅速增大，且观测方程系数矩阵的性质也快速变差；需要通过迭代来减小线性化误差的影响。短弧边值法以轨道弧段两个端点的位置向量为初值，利用先验力模型内插解算出整个弧段参考轨道相对于几何轨道的改正量。这样求得参考轨道与几何轨道非常接近，因此解算重力场模型时忽略计算保守力时参考轨道改正量的影响。然而，由于改正量的数目随弧段增长而迅速增加，且改正量解算方程的稳定性也迅速变差，因此短弧边值法的弧段不能太长，通常只有2~3 h。加速度法将相邻3个历元的卫星位置和星间距离观测值进行二阶差分求得卫星加速度和星间相对加速度，并表示成3个历元对应弧段参考轨道加权平均加速度的线性摄动量，因此加速度法线性摄动量的观测方程形式相对简单，然而差分计算必然放大观测误差，通常需要进行滤波处理^[22]。根据文献[23]以几何轨道为初值进行线性化的思路，文献[24]对保守力模型中的几何轨道直接引入误差改正数，并以几何轨道为初值进行线性化，在解算重力场模型的同时估计轨道改正数。因此在解算地球重力场时不需要再计算参考轨道，且理论上更加严密和自洽，并用该线性化方法改进了短弧边值法和加速度法，解算了重力场模型^{[13-14, 17]}。

为了满足相关学科的科学研究和全球气候变化分析与灾害评估的需求，时变地球重力场需要达到100 km空间分辨率、1天至数天的时间分辨率且精度提高1个数量级，因此，国际上正在研究新一代重力卫星观测计划^[25-27]，我国不仅在研制低低跟踪模式的重力卫星^[28]，且也在研究下一代重力卫星计划^[29]。下一代重力卫星计划采用两对低低跟踪重力卫星进行观测，并采用激光干涉以纳米级精度测定重力卫星的星间距离，这需要改进重力反演解算方法，必须将计算误差控制在纳米级精度。卫星重力反演方法的发展主要在以下3个方面：① 采用长弧重力反演。因长弧段积分不仅可平滑随机误差的影响，而且使微小的作用力产生可观的位置变化，理论上可以解算更高精度的重力场模型，但需要构建新的轨道计算方法，使轨道计算误差能够满足纳米级观测精度的要求，并控制轨道积分误差随弧段长度的积累；② 改进参数化方式。其中，地球重力场采用点质量模型表示，非保守力的偏差参数采用2次多项式或样条模型表示，并将非保守力变换到惯性系时对姿态引入参数进行改正；③ 构建约束解模型，主要利用海洋区域和陆地水流域等不同区域质量变化的频谱特性构建约束模型。此外，大气和海洋混叠改正模型将是影响下一代重力卫星解算精度重要因素，需要利用各种类型数据构建更精确的改正模型。由于动力学法是目前卫星重力反演的主要方法，因此研究动力学法的特点并合理改进，不仅对提高GRACE数据的重力反演精度有参考价值，而且对下一代重力卫星的数据处理也有重要意义。本文主要基于理论模型，分析讨论动力学法反演地球重力场的特点，在此基础上提出一些改进设想。

1 理论基础

地球引力位V在地固坐标系中可表示成如下球函数展开式

式中，r、θ、λ分别为球坐标的向径、余纬和经度；GM为万有引力常数与地球质量之积；R为地球平均半径；P_nm为缔合勒让德函数；u={… C_nm, S_nm …}称为重力场模型系数。由于地球引力是低轨卫星的主要作用力，如果对日月引力、大气和海洋潮汐等保守力精确建模，对大气阻力、太阳光压等非保守力精确测定，就可根据低轨卫星轨道观测值解算地球重力场模型系数。

利用GRACE卫星高精度的轨道、星间距离或速度，以及非保守力观测量解算地球重力场的理论基础是如式(2) 所示的牛顿运动方程

式中，r为惯性坐标系中的卫星位置向量，是时间t的函数；∂V(r, u, t)/∂r和f(r, dr/dt, p, t)分别为卫星单位质量所受的引力和其他摄动力；a(r, u, p, t)为两者之和；u为待估的地球重力场模型系数；p为其他待估参数。在式(2) 的力模型中，所有保守力都与卫星位置有关；但当大气阻力等非保守力用加速度计实测时，式(2) 右边的力模型与卫星速度无关。动力学法和短弧边值法都是以式(2) 积分得到的卫星速度和位置为基础建立GRACE卫星的位置和星间速度观测方程；二者的区别主要是对6个积分常数和参考轨道的处理。

2 基于轨道积分公式线性化的重力解算模型

动力学法以轨道起始点的位置r₀(t₀)和速度为初值，对式(2) 进行积分得到如下卫星位置和速度的严密表达式

式中，t₀和t为轨道的起始时刻和观测时刻；上标“·”表示对时间的一阶导数，参数x=(u^T, p^T)^T。若将(3) 式中的起始点初值、力模型参数用其近似值和改正数表示，则式(3) 可线性化为

式中，I为3×3阶单位阵；r₀、和x₀分别为初值r₀(t₀)、和参数x的近似值，相应的改正数为δr₀、和δx; A(τ)=∂a/∂r(τ)、B(τ)=及C(τ)=∂a/∂x。在t时刻的参考轨道r₀(t)和由近似值按式(6) 计算

引入辅助量

则其对时间t的一阶导数可表示为

在初始时刻t=t₀时，式(7) 和式(8) 积分的上下限相同，显然有

将式(7) 和式(8) 分别代入式(4) 和式(5) 得

式中

根据泰勒级数展开理论，式(10) 中的展开系数可以表示为Y(t)=∂r(t)/∂x、∂x、R(t)=∂r(t)/∂r₀、、以及。因此，式(7)、(8) 和(11) 中的偏导数∂r(τ)/∂x、、∂r(τ)/∂r₀、、和分别可用Y(τ)、、R(τ)、、S(τ)和表示。

由于在初始时刻式(9) 为0，这表明摄动量δx对式(10) 中初始位置和速度的摄动影响均为0。将式(10) 的第一式的r(t)用轨道观测值r_g(t)和改正数v_r(t)表示，则其观测方程为

由式(11) 可得，和S(t₀)=0，则式(12) 在初始时刻可简化为：r_g(t₀)+v_r(t₀)=r₀+δr₀。因此，如果用轨道观测值r_g(t₀)作为初值计算参考轨道，则初值的改正量δr₀必然等于观测值的改正数v_r(t₀)。两颗GRACE卫星A和B之间的相对速度观测量及其改正数与两颗卫星的速度向量和的关系如下

式中，e_AB(t)为A、B卫星视线方向的单位向量。将式(10) 的第2式代入式(13) 就可求得星间速度的观测方程。

不难验证，式(7) 和式(8) 就是动力学法关于参数x的变分方程

满足初值条件式(9) 时的解。而且，式(11) 也是动力学法关于初始位置和速度的变分方程

满足初值为单位阵时的解。由此可见，动力学法本质上是以参考轨道为初值的线性摄动解，摄动参数δx对起始点位置和速度的摄动量影响为0必然导致变分方程式(14) 的初值为0。文献[30]所讨论问题的本质上并非是卫星相对于参考轨道线性摄动问题，因此认为变分方程初值不应该为0。

3 算法特点与改进设想

由动力学法观测方程式(12) 和式(13) 可以看出，GRACE数据解算地球重力场的核心是卫星位置和速度的线性摄动公式(10)，以及参考轨道积分公式(6) 和系数计算公式(7)、式(8) 和式(11)。下面从线性化误差、参数化方式和长弧段轨道积分三个方面讨论算法的特点和改进设想。

3.1 线性化误差

由于轨道初值、力模型参数都存在误差，因此式(6) 积分计算的参考轨道r₀(t)和与实际轨道的偏差将随积分时间而快速增大，式(10) 只用线性项表示该偏差量，其位置和速度向量的线性化误差的二次项分别为

式中，O²(·)表示二次项的影响，记号∂(·)/∂x^T=[∂(·)/∂x]^T。R(t)、S(t)、和对r₀^T和的偏导数都是3×3×3解立体矩阵，对x的偏导数是m×3×3阶立体矩阵，其中，m为参数x的数目；Y(t)和对x^T的偏导数是m×3×m阶立体矩阵，对r₀^T和的偏导数都是3×3×m阶立体矩阵。根据式(7)、式(8) 和式(11)，不难导出式(16) 和式(17) 中18个偏导数的表达式。因偏导数数目太多，不给出所有这些偏导数表达式，只利用∂Y(t)/∂x^T说明线性化误差的特性。考虑到GRACE数据的非保守力是实测的，力模型中的加速度与卫星速度项无关，即式(7) 中B(τ′)=0，因此式(7) 对x的偏导数为

式(8) 应用了矩阵对向量的导数规则以及立体矩阵的运算规则，并顾及了矩阵C(τ′)仅是位置r(τ′)的函数。偏导数∂A(τ′)/∂r(τ′)和∂C(τ′)/∂r(τ′)约为A(τ′)和C(τ′)的1/r(τ′)，说明二次舍弃项约为线性项的1/r(τ′)。然而，式(7) 中Y(t)经过时间t的二次积分得到，且积分号中也包含Y(τ′)项，因此Y(t)接近按t的4次方速度递增。式(18) 包含了Y(τ′)乘积项的时间二次积分，这意味着偏导数∂Y(t)/∂x接近按时间10次方的速度递增。同理可发现接近按时间9次方的速度递增，其他各偏导数的递增速度也都约为时间的9与10次方。由此可见，尽管动力学法线性化误差的数值约为线性项的1/r，但随轨道积分弧段长度快速递增。

根据式(16) 和式(17)，欲减小线性化误差，可从两个方面着手，即控制参数改正量δr₀、、δx和各项偏导数的数值。通过迭代计算减少参数改正量是目前动力学法减少线性化误差的主要手段，但各偏导数的数值并不能通过迭代来控制。由于各偏导数的数值与线性化方式密切相关，文献[23]以几何轨道观测值为初值的线性化方法，使式(4) 和式(5) 中的轨道改正值δr(τ′)不会随时间而快速递增，可有效减小各偏导数数值随时间的积累速度；然而其轨道改正值的处理方式涉及四重积分并使常数项的误差模型更加复杂化，这可能是该方法至今不能实际应用的主要原因。由于以轨道观测值为初值进行线性化时，轨道改正值必定与平差后观测误差改正数相同，因此文献[31]将δr(τ′)表示成误差改正数v_r(τ′)，用GRACE数据成功解算了180阶次的重力场模型。然而，按文献[31]的处理方式，其误差改正数不能用未知参数改正值直接表示，导致长弧段重力反演的求逆运算量太大。因此，以轨道观测值为初值进行线性化，还需要进一步改进算法，或建立合理的迭代方法，控制各偏导数的随时间积累速度，有效减小长弧段重力解算时的线性化误差。

3.2 参数化方式

动力学法采用轨道起始点的位置和速度向量参数，重力位参数和非保守力加速度的尺度与偏差参数进行参数化。这种参数化方式导致式(10) 中参考轨道与观测轨道之差随时间快速增加，式(7)、式(8) 和式(11) 中的系数矩阵R(t)、S(t)、Y(t)按时间4次方的速度递增，、、按时间3次方的速度递增；且轨道弧段越长，不同频段重力位参数的平滑影响相差也越大；这意味着动力学法线性化观测方程关于重力位参数的系数矩阵性质也会随弧段增长而变坏。根据式(11) 不难判断，位置和速度向量的系数矩阵R(t), S(t)的数值也随弧段增长相差甚大，因为即使忽略式(11)R(t), S(t)计算公式右边积分结果的差异(因积分号内各项类同)，1天弧段R(t)与S(t)相差也要约86 400倍，这对解算线性方程显然是不利的。

改进参数化方式，有望改善观测方程的性质。改善位置和速度系数阵性质的最简单办法是改变其速度单位，若以1天为弧段长度，只要取km/s为速度单位即可。动力学法迭代收敛后，位置、速度和力模型参数改正量都趋向于0，观测方程式(12) 的常数项主要反映了观测误差的影响。由于迭代收敛后，参考轨道的平差改正量完全可以忽略，即式(4) 和式(5) 中的δr(τ′)和的影响可以忽略，因此式(7) 和式(8) 可简化为

显然，式(19) 中的Y(t)和分别按时间的二次和一次的速度递增，意味着其性质将优于式(7) 和式(8) 的计算值。因此，动力学法迭代收敛后，利用式(19) 计算观测方程系数后再解算一次，有望能够进一步提高动力学法的解算精度。

由于短弧边值法可直观理解为轨道内插问题，动力学法解算初值参数属于轨道外推问题，因此短弧边值法观测方程性质应该优于动力学法。考虑到迭代收敛后弧段中各点的轨道已经固定，边值公式已经不受弧段长度的限制，因此用短弧边值公式作最后一次解算，也有望进一步提高解算结果的精度。

3.3 轨道积分方式

由于卫星的力模型复杂，给出完整的轨道解析解非常困难，因此轨道数值积分对于动力学法是不可避免的。尽管采用目前常用数值积分器积分1天的轨道精度能满足目前GRACE重力反演的精度要求。然而下一代重力卫星的激光干涉观测精度将达到纳米量级，比目前的K波段微波测距精度提高约1000倍，因此需要更高的轨道积分精度。长弧重力反演计算可提高观测方程的信噪比，有利于提高重力场参数的解算精度，但轨道积分精度必须要满足长弧段解算的要求。因此，卫星轨道的数值积分精度是下一代重力卫星数据处理和长弧段重力反演的一个重要瓶颈问题。

式(3) 的第1式是二重积分，利用分步积分可转化成如下一重积分

利用式(20) 进行数值积分，其计算误差将明显小于式(3)。当采用多步法积分时，式(20) 积分所需的力是利用积分区间附近若干离散点的力，通过内插或外推求得。在动力学法定轨时，只已知t时刻之前的轨道，只能采用外推计算；但在卫星重力反演计算时，卫星轨道是已知的，可内插计算。基于高阶牛顿插值求得的12阶外推积分系数和内插积分系数如下。

外推系数：19 494 601/11 404 800, -99 642 413/22 809 600, 40 413 623/2 851 200, -4 955 916 683/159 667 200, 278 428 507/5 702 400, -4 496 090 419/79 833 600, 955 625 177/19 958 400, -2 374 517 119/79 833 600, 1 050 348 479/79 833 600, -627 827 071/159 667 200, 84 671/118 800, -4671/78 848

内插系数：-4139/79 833 600, 19 567/22 809 600, -107 539/13 305 600, 17 274 001/159 667 200, 6 386 783/7 983 360, 326 441/3 193 344, -53/399 168, -58 213/11 404 800, 226 637/79 833 600, -4121/4 561 920, 317/1 900 800, -317/22 809 600

其中，外推系数的数值和正负变化量要明显大于内插系数，因此，内插系数的数值积分精度将明显优于外推系数。低轨重力卫星的主要作用力是地球引力，相对于地球中心引力，C₂₀项引力约为10^-3量级，其他各项引力最大只有10^-6量级，日月引力等保守力摄动和非保守力摄动也只有10^-6量级。若将式(20) 的力模型分离成中心引力与C₂₀项引力项和剩余部分，则式(20) 可改写为

式中

如果只考虑地球中心引力和C₂₀项引力作用，可导出卫星运行的解析公式^[32]，因此理论上式(22) 可用解析公式表示，只需要对式(21) 进行数值积分计算。顾及约只有10^-6量级，用式(21) 计算的参考轨道精度将有望大幅提高。然而，实际计算时还有不少细节问题需要处理，如式(22) 中的r(τ), 是卫星的实际位置和速度，解析公式的理论轨道是由轨道初值和中心引力与C₂₀项引力计算得到，因此式(22) 积分结果与解析公式计算结果之间必定存在差值，根据文献[33]，该差值与轨道之比必须要控制在1%以内。如何改正该差值使其影响可以忽略，以及式(21) 中用几何轨道计算力模型引起的误差是否可以忽略等问题，都需要作进一步的理论研究和数值模拟分析。

4 结语

本文直接利用动力学轨道的积分公式进行线性化，导出了线性化观测方程系数矩阵的积分计算公式，阐明了传统动力学方法本质上是相对于参考轨道的线性摄动问题，因此对力模型参数的偏导数初值必定为0。利用导出的积分公式分析了线性化误差随轨道弧长快速增大，线性观测方程性质也随轨道弧长增长也变差等特点。建议通过迭代计算和进一步改进以几何轨道为初值的线性化模型，减少长弧段重力反演的线性化误差；通过改变初始速度单位改善初始状态参数系数矩阵的性质，采用迭代收敛后轨道直接计算设计矩阵，改善线性观测方程的性质。考虑到下一代重力卫星计划采用高精度星间距离观测数据和长弧段重力反演计算都需要高精度数值积分计算卫星轨道，建议保守力模型采用内插法计算，并将地球中心引力和C₂₀项引力采用解析公式计算，其他保守力和非保守力采用数值积分公式计算，以提高卫星轨道的计算精度。

本文并未涉及卫星重力反演的星座和载荷指标的优化、非保守力改正模型的精化、大气和海洋改正混叠影响等，这些问题都是下一代卫星重力计划需要研究的重要内容。