地球重力场及其时变特性反应了地球层及内部的密度分布和物质运动状态，获得高精度重力场在地球科学领域具有深远的意义.重力卫星以其快速、高效和几近全球覆盖的特性，在确定高精度重力场方面正发挥着重要作用.基于卫卫跟踪技术的GRACE 重力卫星可提供高精度和高空间分辨率的静态及时变地球重力场，其信息可有效揭示地球系统质量重新分布，具有巨大的科学价值和应用前景^[1].为实现在GRACE 任务结束后对地球物质变迁的连续监测，国际上相继展开对GRACE-followon计划的论证，预期以更低的轨道、更高精度的星间测距系统来获取更高空间分辨率的地球重力场，其结果将极大的推进人类对地球物质变迁、重大自然灾害监测、测绘及地球动力学机制等方面的认知^[2].

卫卫跟踪技术是当前地球重力场测量最有价值和最具应用前景的方法之一，我国亦准备实施自主知识产权的重力卫星测量计划.精密星间测距系统是低低卫卫跟踪重力卫星的关键有效载荷.GRACE卫星携带的K 波段测距系统是一微米量级的测距系统，GRACE-followon计划又提出了一种更高精度的激光干涉测距系统.精密星间测距系统通过高精度星间距离观测感应地球重力场分布，由于更高阶次的星间距离变率和距离加速度对地球重力场的高频信息更加敏感，地面数据处理中心需对测距系统直接测得的高采样率下的星间距离观测值进行滤波处理，得到后处理的星间距离、距离变率和距离加速度数据，用于后续地球重力场模型解算^[3].目前国内外许多科研机构和学者投身于精密星间测距系统观测数据的处理技术研究当中，对测距系统观测模式、误差源模型、测距系统精度与地球重力场恢复精度的匹配关系等进行了研究和论证^[4-10].对于星间测距系统滤波技术的研究，MacArthur等采用了滑动平均算法实现滤波处理^[11]，Thomas等给出了基于最小二乘拟合原理和窗函数法的滤波器设计^[12]，GRACE 地面数据处理中心采用了矩形窗函数自卷积技术构造滤波器^[13]，用来实现星间瞬时距离观测数据的去噪和降采样处理.本文拟从精密星间测距系统观测值与地球重力场频谱关系的角度，建立星间距离观测值关于重力位系数的敏感矩阵，详细分析其频谱特性，得出测距观测值可感应的地球重力场信号的有效频带，最后结合美国当前GRACE 卫星和下一代GRACEFollow-On 卫星计划、以及我国首期卫星跟踪卫星重力测量计划，给出低通滤波器的通带截止频率、通带增益波纹和频率采样率技术指标设计方案.该滤波器设计方案是基于星间测距观测值敏感矩阵与地球重力场频谱关系的建立与分析的基础上的，能最大限度保留地球重力场信息，避免了单一从仪器观测高频噪声抑制的角度设计降采样滤波器.重力卫星精密星间测距系统滤波器指标设计及优化论证属于系统预研工作的范畴，可直接用于我国地面重力数据处理系统建设.

地球重力场引力位V的球谐函数展开式为

(1)

式中，

(2)

(3)

其中，φ、λ 为空间点的地心纬度与经度；P_nm为Legendre缔合函数；C_nm、S_nm为地球引力位系数，r_e为地球平均半径.

卫星S₁、S₂ 的加速度矢量可示为^[14]:

(4)

其中，E^T 为从地固系到惯性系的转换矩阵，∇为地球引力位的梯度，V_i为卫星S_i处的地球引力位，a_{nc_i}为摄动加速度，a_{ε_i}为仪器观测噪声.

设卫星S₁ 与S₂ 之间的视线单位矢量为e₁₂，则星间距离加速度(相对视线加速度)可近似示为^[15]:

(5)

将公式(1)、(4)代入(5)可得基于地球引力位系数和卫星轨道坐标的距离加速度$\ddot{\rho }$ 的解析达式.

敏感矩阵用来定量评价模型参数对模型结果产生的影响^[16]，距离加速度对n阶m级重力位系数的敏感矩阵T_nm可征距离加速度对地球重力场的敏感程度，表达式为^[12]

(6)

式中，

(7)

(8)

由式(6)可知敏感矩阵T_nm为距离加速度$\ddot{\rho }$ 按n阶m级重力位系数C_bnm球谐展开的系数项，即对于给定的卫星轨道，距离加速度$\ddot{\rho }$ 可示为

(9)

依据频谱分析的线性叠加原理，$\ddot{\rho }$ 的频谱可由敏感矩阵T_nm的频谱线性叠加得到，我们可以通过敏感矩阵T_nm来研究距离加速度观测值$\ddot{\rho }$ 的频谱特性.由信号频谱的时域微分定理知，星间距离、距离变率相应的敏感矩阵的频谱可由上述T_nm的频谱分别乘以因子1/(2πif)²、1/2πif得到^[17].由于精密星间测距系统直接测得的物理量为瞬时星间距离，卫星地面数据处理中心通过滤波和降采样将该观测值转换为后处理星间距离、距离变率和距离加速度，滤波技术的参数指标设计应基于重力卫星星间距离观测值的频谱特征.

本文采取GRACE 卫星开普勒轨道根数模拟卫星运行16个周期内的理论轨道数据，并依据式(6)计算了敏感矩阵中T_2，0、T_10，0、T_10，5 等12个不同阶次的敏感因子的时间序列，对其分别作离散傅里叶变换并乘以因子1/(2πif)² 得到星间距离敏感矩阵的频谱.敏感矩阵中重力位二阶导数的计算采用了非奇异的算法^[18]，有效地避免地球两极附近的奇异性.采用的卫星轨道参数指标^[19]如表 1，重力卫星星间距离敏感矩阵频谱如图 1，其振幅谱相关特征参数如表 2.

表 1(Table 1)

表 1 卫星轨道参数 Table 1 Parameters for satellite orbit

表 1 卫星轨道参数 Table 1 Parameters for satellite orbit

表 2(Table 2)

表 2 振幅谱相关特征参数 Table 2 Characteristic parameters of amplitude spectrum

表 2 振幅谱相关特征参数 Table 2 Characteristic parameters of amplitude spectrum

图 1

Fig. 1

图 1 重力卫星星间距离观测值敏感矩阵频谱 Fig. 1 Spectrum of sensitivity factors for inter-satellite range measurements of gravity satellite

由图 1和表 2可知，对于n阶带谐系数，敏感矩阵T的频谱峰值个数为n/2(若n为奇数，频谱峰值个数为(n+1)/2)；主频可近似用n/rev(rev≈5430s)估算，有效频带范围近似达为2/rev~n/rev；频谱峰值随着频率的增加而增大，主频为有效频带上限值.对于n阶扇谐系数、田谐系数，敏感矩阵T的主频、峰值等频谱特征参数没有统一的解析式定量达，在此仅作定性分析.对于n阶球谐系数，敏感矩阵T_n，0 的频谱以高频信息为主，主频峰值约比T_n，m(m≠0)大1~2个数量级，随着m的增加T_n，m的高频信息逐渐减弱，主频峰值大致呈减小的趋势，T_n，n的主频峰值为最小，以低频信息为主(如图 2).上述敏感矩阵频谱成分随次数m的变化趋势明:对于n阶球谐系数，重力卫星星间距离观测值的主频可近似用敏感矩阵T_n，0 的主频估算；地球重力场信号有效频带上限由T_n，0 的频谱决定，下限由T_n，n的频谱决定.

由上述讨论可知，敏感矩阵的主频由球谐函数的空间周期和卫星轨道时间周期共同决定.n阶带谐函数二阶导数在0~π之间有n个零点，幅值从中心点(θ=2/π)到两极(θ=0，π)逐渐增大，在乘以加权因子$\sqrt{\cos \theta }$后幅值可近似用正弦函数sin(nθ)拟合^[20].当卫星以圆极轨道运行时，其轨道方程在纬度方向上为θ=2π_t/rev，则当卫星在仅由n阶带谐函数所确定的重力场模型中运行时，其敏感矩阵的时序曲线拟合函数为sin(2πnt/rev)，亦即其主频为n/rev.另外，由n阶勒让德函数的θ 倍角函数公式不难得到上述n阶带谐函数对应的敏感矩阵频峰个数的结论(P_n(cosθ)的二阶导数由n/2 个(若n为奇数，个数为(n+1)/2)加权正弦函数叠加而成^[21]).

对于n阶扇谐系数、田谐系数，敏感矩阵T的频谱特征同样可由勒让德函数的性质得到.对于同一阶数n，缔合勒让德函数P_nm(cosθ)随着m的增大，零点个数逐渐减小，幅值逐渐向中心(θ=0)集中；当m=n时，零点消失，幅值集中在中心附近较狭小的区域^[22].球谐函数相当于将上述仅与纬度相关的缔合勒让德函数乘上了一个与经度相关的权重cos(mλ)或sin(mλ)，其幅值在单位球面上的分布可看作是围绕赤道分布的正负相间的条带，该条带随着m的增大而逐渐向赤道集中^[23]，相应地，敏感矩阵T频谱成分逐渐由高频向低频过渡.为验证上述结论，笔者模拟计算了圆极轨道运行的卫星在仅由n阶m次球谐函数确定的重力场模型中的敏感矩阵，发现对于同一阶数n，不同级数m对应的敏感矩阵的频带宽度大致相同，但随着m的增大，其主频峰值的频率点逐渐向零频靠近，说明低频信息逐渐显著，而高频信息逐渐削弱.上述结论与 1 和表 2 中揭示的现象相吻合.

图 2

Fig. 2

图 2 重力卫星星间距离观测值敏感矩阵(10阶球谐系数)的频谱峰值 Fig. 2 Peak value of spectrum for sensitivity factors of intersatellite range measurements of gravity satellite (10 degrees)

为最大限度保留星间测距观测值中的地球重力场信息，滤波器通带截止频率应大于观测信号有效频带的上限频率值.由第二节基于GRACE 卫星开普勒轨道根数的敏感矩阵频谱分析结果知:当地球重力场用n阶球谐函数逼近时，星间距离观测值的有效信号频带上限主要由T_n，0 的主频值决定.实际观测中还应考虑卫星轨道摄动加速度和仪器观测噪声的影响.

由GRACE 卫星摄动加速度频谱分析结果知^[24-25]:卫星正常轨道加速度频谱的主频幅值量级为10^-2m·s^-2；地球潮汐、日月引力以及其它行星的第三体引力引起的保守力摄动加速度频谱的主频幅值量级为10^-7~10^-14m·s^-2；大气阻力、太阳辐射压和地球辐射压引起的非保守力摄动加速度频谱的主频幅值量级为10^-7~10^-8m·s^-2.摄动加速度频谱中含有与正常轨道加速度观测值同频的信号，但由于其主频幅值比后者小5~12 个量级，在确定星间测距观测信号有效频带范围时可直接用基于正常轨道的加速度敏感矩阵进行分析.另外，星间测距系统观测误差相关研究结果明^[26]:由振荡器噪声、系统噪声和多路径效应引起的观测噪声主要以高频噪声为主，对星间观测值的数据质量及重力场模型恢复的精度有很大的影响，在滤波器的阻带截止频率的指标设计中应予以考虑，但仪器观测噪声中与正常轨道加速度同频的噪声和上述重力场模型恢复需要的有效信号混合在一起，无法用滤波的方法彻底分离.基于上述讨论，本文中基于卫星正常轨道的星间观测值敏感矩阵频谱得到的地球重力场有效频带中同时包含了卫星轨道摄动引起的同频信号，以及不可避免的混合了星载仪器观测噪声中与有效信号同频的噪声.

基于美国当前GRACE 卫星运行参数^[19]和下一代GRACEFollow-On卫星计划设计参数^[27]，本文计算了100阶地球重力场模型对应的星间测距信号有效频带分别为:3.523×10^-4Hz~1.760×10^-2Hz和3.767×10^-4Hz~1.880×10^-2Hz, 低通滤波器通带截止频率应至少大于该频带上限.我国首期卫星跟踪卫星重力测量计划滤波器的指标设计中除了考虑上述有效信号频带之外，还要综合考虑滤波器长度、滤波后数据采样率等因素来确定合理的通带截止频率.

滤波器通带增益为以一定容限的允许误差来征的通带幅度响应指标.GRACE 观测任务中重力场模型大地水准面起伏容限ε=0.01cm^[28]，则n阶重力场位系数对应的大地水准面起伏相对误差容限为由第二节知，n阶带谐系数敏感矩阵可近似将上述大地水准面起伏相对误差转换为星间测距观测值的相对振幅容限.基于Kaula准则^[29]，可得通带增益容限A_n^err的近似值如下:

(10)

图 3 给出了100 阶地球重力场模型(Kaula准则和GGM02S模型^[30])滤波器的通带增益，可以看出，基于Kaula准则的星间测距观测值相对振幅容限与卫星实测模型有很好的拟合度(J2项除外).对于重力场低阶系数J2项的精密确定目前已有大量的研究结果^[31-32]，本文采取静态地球重力场模型EGM96^[33]给出的结果0.484165×10^-3，估算J2项对应的频率点0.00037Hz处增益容限为:

图 3

Fig. 3

图 3 100阶地球重力场模型低通滤波器通带增益 Fig. 3 Pass-band gain of low-pass filter for100 degree gravity field model

f单位为H_{z, αp}单位为dB，该结果与基于Kaula准则的通带增益容限A_n^err的估算结果较为接近.为使滤波器在信号通带内的幅度响应满足重力场模型解算对位系数阶方差的要求，在第n阶球谐系数对应的频率点处的通带增益应尽量小于A_n^err，滤波器最大通带增益应小于max(A_n^err).

n阶重力场模型低通滤波器频率响应容限如图 4所示(中给出了n=100的低通滤波器通带允许逼近误差的极限).在设计满足上述技术指标的滤波器时，我们必须确定离散时间线性系统的系统函数，该系统的频率响应应落在预先给定的容限内，即通带增益应该在1+A_n^err和1-A_n^err之间变化(图 4中灰色区域).我国首期卫星跟踪卫星重力测量计划的大地水准面精度设计指标比当前的GRACE 计划高1~2个数量级，相应地，通带增益容限要求也更加苛刻.在实际滤波器设计中应在地球重力场模型所要求的频率响应容限内综合考虑滤波器阶数、通带波纹特性等因素尽量逼近理想滤波器.

图 4

Fig. 4

图 4 100阶地球重力场模型低通滤波器频率响应容限 Fig. 4 Tolerance of low-pass filterS frequency responsefor 100 degree gravity field models

根据尼奎斯特采样定理，当可检测的最高阶频谱信号一定，为保证此频率信号的正常恢复，要求卫星重力观测保持相应的数据采样率.对于最高阶为nmax 的重力场模型，其空间分辨率为λ = πr_e/n_max.实际采集数据时观测采样率应高于标称采样率，若要求卫星重力数据有1.0°×1.0°分辨率，则一般应在1.0°×1.0°格网内均匀进行4点重力测量，即经度和纬度方向的空间采样率应达到相应阶次球谐函数相邻零点值区域内每0.25 相位至少有一个采样点，也即满足n_max/rev·s≤0.25，式中s为频率采样率.相应地，滤波器频率采样率至少大于4nmax/rev才能满足恢复n阶地球重力场模型的需求.

我国首期卫星跟踪卫星重力测量计划重力场空间分辨率设计指标为66km, 可恢复300 阶重力场模型，基于上述指标设计方案，星间观测值滤波器频率采样率至少应大于0.2236 Hz.恢复更高阶次的重力场模型，总是要求相应的更密集更均匀的全球空间覆盖的观测数据，但更高的数据采样率对系统软硬件提出更苛刻的要求.在滤波器设计中，应综合考虑用户对重力场模型的空间分辨率的需求和系统软硬件条件，确定合适的滤波器频率采样率.

利用重力卫星可获得全球一致的高精度、高分辨率的地球重力场及其时变参数，丰富人们对地球重力场的认知程度和拓展地球重力场的应用领域^[34-35]，科学的数据处理技术是获得高质量地球重力场的前提.本文基于重力卫星精密星间测距测量模式，对星间测距系统观测值敏感矩阵的频谱特性作了定性分析，在讨论了星间测距观测值中能够反应地球重力场信息的有效信号频带的基础上，给出如下低通滤波器的技术指标设计方案:(1)当地球重力场用n阶球谐函数逼近时，低通滤波器的通带截止频率应至少大于n阶带谐系数敏感矩阵因子T_n，0的主频，该主频由球谐函数的空间周期和卫星轨道时间周期共同决定；(2)低通滤波器在有效频带各频率点处的通带增益容限可用基于Kaula准则的星间测距观测值相对振幅容限估算，滤波器通带内的频率响应应落在预先给定的容限内；(3)滤波器频率采样率至少大于4n_max/rev才能满足恢复n阶地球重力场模型的需求.

由于滤波器阻带技术指标与星上载荷的仪器噪声水平、卫星运行空间环境背景噪声等因素有关，本文只讨论了基于地球重力场信号频谱特性的通带截止频率、通带增益和频率采样率的设计指标，随着我国自主产权的重力卫星工作的全面开展，有望得到更精细的数据处理分析系统技术指标论证结果.