利用卫星轨道数据反演重力场模型的常用方法主要有能量守恒法、天体力学法、短弧积分法和加速度法。能量守恒法对卫星速度精度要求较高，且建立的标量观测方程损失了卫星状态矢量的方向信息，天体力学法计算复杂且耗时大，短弧积分法不利于恢复低阶次重力场模型，而加速度法基于牛顿第二运动定律直接建立加速度与位系数的线性关系，从而估算位系数，其原理简单、直观，相比其他方法具有一定的优势。但卫星加速度值不是直接观测量，而是由卫星轨道经过两次数值微分计算得到的间接观测量，因此加速度法的关键是如何利用卫星轨道数据微分得到所需精度的加速度数据。沈云中^[1]、Reubelt^[2-4]、徐天河^[5]、钟波^[6]均采用牛顿数值微分公式计算卫星加速度值，并分别应用于CHAMP卫星模拟和实测数据以及地球重力场和稳态海洋环流探测 (gravity field and steady-state ocean circulation explorer，GOCE) 模拟数据解算中。牛顿微分公式形式上简单，但在实际应用中特别是编程上不易实现且公式繁琐，因此本文引入新的实用算法进行微分计算。孙海燕等^[7]进行全球定位系统 (global positioning system, GPS) 卫星摄动力数值分析由精密星历计算卫星加速度时，采用连续17个历元的星历作拟合，并计算中间节点处的加速度；游为^[8]针对力模型的积分计算问题，采用基于移动窗口的9次多项式内插法进行内插并积分。鉴于以上思路, 本文引入移动中心窗口多项式数值微分算法来计算卫星瞬时加速度，通过模拟数据分析多项式阶数、观测点数、采样间隔对微分结果的影响，并采用GOCE卫星2009年11月1日至2010年1月10日的实测几何法轨道 (precise kinematic orbit，PKI) 数据^[9]恢复至100阶GOCE重力场模型GOCE-PKI-AA，验证了该算法的有效性，最后通过阶误差指标评定计算模型的精度。

假定卫星几何位置中X方向 (Y、Z方向以此类推) 的时间序列可由N阶多项式插值公式表达为：

式中，r(t₀+τ_i) 表示观测历元的位置分量；c_j为位置系数；多项式阶数N=2n；观测点数M≥N+1；Δt为采样间隔；t₀表示计算时刻；τ_i表示计算时刻与实际观测时刻的差值。

为了保证观测历元关于中心时刻对称分布，设定插值多项式阶数为偶数，观测点数为奇数。这样也可以降低多项式插值的截断误差，提高微分计算精度。式 (2) 第一式为式 (1) 的矩阵形式，式中待定系数c =[c₀ c₁ c₂ … c_N+1]^T; 设计算时刻位于中心处，此时系数阵A为范德蒙特方阵，其矩阵元素表达式为式 (2) 第二式。

利用经典最小二乘平差法，可得：

式中，W=(A^TPA)^-1A^TP。其中，P为观测量的权阵，本文设为单位矩阵。当给定多项式阶数、观测点数、采样间隔时，W阵元素恒定不变。可见待定系数实质上是观测量的线性组合。对式 (1) 进行两次微分，可得：

仅考虑中心时刻，此时τ_i=0，则w=(0 0 2 0 … 0)^T。将式 (3) 代入式 (4) 得：

式中，e=(w^T·(A^TPA)^-1A^TP)^T=P^TA(A^TPA)^-1·w。由式 (5) 可知，计算时刻的加速度值实际上只与多项式插值待定系数中的c₂元素有关。可见加速度值仅是常数W阵中第3行元素值与观测量的线性组合。当给定多项式阶数、观测点数和采样间隔时，微分算子e向量元素均为常数，从而观测弧度只需要计算一次，这样可以大大提高计算速率。上述过程称为移动中心窗口多项式数值微分算法，如图 1所示。该算法的输入量为卫星轨道数据，即星历位置，输出量为中心时刻的瞬时加速度值。其实质是滤波过程，由图 1可知，在某计算时刻左右两边各取相同观测点数，加上该中心时刻的观测点构成一个完整的窗口函数；再结合多项式插值对窗口内的观测点进行最小二乘拟合，通过线性组合滤波核函数与观测值来计算中心点；随后移动窗口，接着计算下一个观测点，以此类推。其中，滤波核函数为常数阵W中第3行元素值的两倍，最终加速度值可以表示为：

图 1(Figure 1)

图 1 移动中心窗口多项式数值微分算法示意图[10] Figure 1 Application Scheme of Numerical Differential Algorithm with Polynomial Kernel of a Moving Window Function

式中，常数阵E满足下式：

由于GOCE卫星运行时加速的真值未知，利用轨道数据采用移动中心窗口多项式数值微分算法获得的卫星加速度计算值，难以判定其精度。因此笔者首先采用截断至180阶次的EGM2008模型以1 s的采样间隔进行数值积分模拟一天GOCE卫星轨道，共计86 400个历元，同时基于EGM2008模型模拟了该轨道的卫星加速度，将其作为真值分析算法的精度与稳定性，模拟卫星轨道的星下点轨迹如图 2所示。为了使添加到模拟轨道的噪声更接近GOCE卫星实测数据的情况，本文将实测的几何轨道与简化动力学轨道求差后的结果加入到模拟轨道中。

图 2(Figure 2)

图 2 模拟一天GOCE卫星轨道星下点轨迹 Figure 2 Subsatellite Track of Simulated GOCE Satellite Data with One Day

如前所述，移动中心窗口多项式数值微分算法与多项式阶数、观测点数和采样间隔有关。无噪声和有噪声情形时不同因素对微分加速度精度的影响如图 3所示。结果表明：①一定范围内，随着阶数和采样间隔的增大，加速度精度逐渐提高，但阶数过高，信号过度拟合且算法稳定性差，采样间隔过大，信号过度平滑导致精度降低；②随着观测点数增大，精度急剧下降，这主要是微分算法对误差的高频放大效应所引起的，无论是无噪声还是有噪声，均使得10 s采样间隔8阶9点多项式微分后加速度精度最高，无噪声的情况达到10^-9 m/s²，有噪声时精度为10^-5 m/s²，故笔者在下文采纳8阶9点10 s采样间隔的移动窗口多项式数值微分公式来计算加速度值。

图 3(Figure 3)

图 3 不同因素计算加速度残差图 Figure 3 Residual Accelerations Diagram with Different Factors

本文基于上述算法，采用欧空局官方发布的2009年11月1日至2010年1月10日共71天的GOCE卫星实测几何法轨道数据恢复至100阶重力场模型GOCE-PKI-AA。以高精度的EIGEN-5C作为参考模型，计算其与求解位系数的差异，并对差值取以10为底的对数，如图 4所示。由图可知，GOCE-PKI-AA模型在低阶扇谐位系数精度较高，但低次位系数精度偏低，主要是由于GOCE卫星轨道倾角为96.7°，两极附近各有约6.7°的空白区域，因此影响带谐位系数的解算精度；GOCE-PKI-AA模型在75~100阶计算精度略低，主要原因是两次微分导致高频误差放大，同时模型最大阶次选择100阶，则在反演时接近100阶次的系数会受到高频误差向低频映射的影响^[11]。

图 4(Figure 4)

图 4 GOCE-PKI-AA与EIGEN-5C模型位系数差值 Figure 4 Coefficient Diffrences Between GOCE-PKI-AA and EIGEN-5C

本文采用不同重力场模型位系数的阶误差评价反演的重力场模型精度，考虑到极空白问题导致低次位系数精度较低，因此利用阶方差评定精度时不考虑5次以下位系数，其计算公式^[12]为：

这里选择EIGEN-5C为参考模型，图 5刻画了不同模型的阶误差，从图中看出GOCE-PKI-AA模型在100阶以前的计算精度优于EIGEN-CHAMP03S模型，主要是由于GOCE卫星轨道高度低于CHAMP卫星轨道，可以感应出较强的重力场信号，同时GOCE卫星定轨精度高于CHAMP；考虑到欧空局发布的第一代纯GOCE卫星重力场模型GO-CONS-GCF-2-TIM-R1^[13]与本文GOCE-PKI-AA模型采用的数据段相同，对比分析可知低阶部分两种模型具有良好的一致性，但在高阶部分，GO-CONS-GCF-2-TIM-R1模型计算精度明显优于GOCE-PKI-AA模型，主要是因为GO-CONS-GCF-2-TIM-R1模型联合了高低卫星跟踪卫星数据和重力梯度观测数据，后者对中高频信号更为敏感。

图 5(Figure 5)

图 5 不同模型阶误差图 Figure 5 Degree Errors of Different Gravity Field Models

对于模拟数据，采用8阶9点10 s采样间隔的移动中心窗口多项式微分算法得到的加速度精度最佳。采用实测的几何法轨道数据验证了移动中心窗口多项式数值微分算法的有效性，并且计算模型与参考模型在80阶以前具有良好的一致性，但带谐系数精度较低，主要是由于GOCE卫星没有全球覆盖，存在极空白问题。利用GOCE卫星几何法轨道数据反演的重力场模型GOCE-PKI-AA在100阶之前较EIGEN-CHAMP03S模型更接近于EIGEN-5C模型，15阶以前较GO-CONS-GCF-2-TIM-R1模型更接近于EIGEN-5C模型。