GRACE(Gravity Recovery and Climate Experiment)卫星发射于2002年3月，通过精密监测同一轨道上的两颗卫星间的距离变化来观测地球重力场，得到比CHAMP(Challenging Minisatellite Payload)重力卫星更高精度的静态地球重力场^[1]；此外,GRACE卫星还能够提供短周期的时变地球重力场。因此，GRACE卫星可用于地球活动的动态研究，包括地球内物质运动、冰川融化、大洋环流等^[2]。

GRACE卫星的Level 2数据产品是将GRACE卫星Level 1的数据解算后得到的重力场解，并截断高阶项，经过多重误差改正，去除或者削弱了大气、海洋、潮汐的影响，精度较高，也是使用最广泛的数据产品^[3]。目前，Level 2数据主要由JPL、CSR、GFZ^[3]三家研究机构处理并分发，分别提供了60~120阶不等的月度全球重力场解。另外还有一些研究机构也根据Level 1数据处理并发布了月度全球重力场模型，如Bonn大学的ITG^[4]、Delft大学的DMT-1^[5]、Bern大学的AIUB^[6]、Graz科技大学的ITSG-Grace2014 ^[7]、CNES/GRGS^[8]以及同济大学的Tongji-GRACE01^[9]等。除月重力场解，也有其他周期的数据产品，CSR提供了30阶的周重力场解，CNES/GRGS也提供50阶的10 d重力场解^[10]。

在使用GRACE Level 2 数据时,还要考虑其中包含的系统误差。首先，GRACE得到的重力场球谐系数中的二阶系数C₂₀并不准确，所以需要用卫星激光测距(satellite laser ranging,SLR)观测得到的C₂₀替换^[11]；其次，根据Swenson等的研究^[12]，GRACE的重力场球谐系数的一阶系数定义为0，在恢复高纬度质量变化和大尺度海洋质量交换时会有显著误差，所以应该再加上一个地心模型来修正一阶系数；最后还需要根据实际研究的问题考虑物理模型改正，比如冰后回弹效应^[13]等。

除以上通过物理模型消除的系统误差外，由于GRACE卫星数据测量的是南北向轨道上两颗卫星的距离变化，所以最终得到的重力场解中存在明显的南北向条带；重力场误差也随着重力场球谐系数阶数的增加而增大；同时，在很多研究中需要关注局部地区的重力变化，因此需要对GRACE得到的全球重力数据进行局部化处理。在应用GRACE Level 2数据时还有大量的数据再处理工作需要完成，针对如何消除条带误差、高频误差和局部化处理等问题，研究者们提出了许多方法。

一般从GRACE Level 1数据处理得到Level 2数据时的各项模型误差改正过程称为预处理过程 ^{[14, 15, 16]}，而针对GRACE Level 2的数据处理则称为后处理过程^{[17, 18, 19]}。本文中，针对GRACE数据的时序分析和局部化应用分析过程均包含在后处理过程中。

基于不同的思路和目的，后处理过程大体上可以分为两种处理策略：一种是对球谐系数直接进行频谱滤波，包含从最通用的各向同性高斯滤波^[20]到更加复杂的与阶次^{[18, 21]}甚至与时间^{[22, 23]}相关的滤波方式；另一种充分顾及具体误差^{[24, 25]}，需要独特的先验信息，尤其是误差协方差矩阵，在最小二乘反演过程中计算出更准确的月重力场解。本文将分类论述这些不同的后处理算法研究进展，并展望未来后处理算法的方向。

不需要先验信息的算法都不依赖月份与信号和噪声的特定结构，可以直接应用于不同月份的数据中，具有普适性，无需先验信息或误差模型，计算量也相对较小，计算速度快。例如应用最广泛的高斯滤波就可以不加先验信息直接应用于任何一个月度重力场解上。

由于GRACE重力卫星月重力场球谐系数的高阶项误差过大，基本无法恢复出有效信息，所以通常只截取重力场前60~120阶的球谐系数作为Level 2数据。即使如此，保留下的高阶系数仍有显著误差。高斯滤波是解决这一问题最简单方便的方法，通过引入滤波因子降低高阶系数的权重可以降低重力场整体误差^[20]。

通过对阶数进行高斯滤波能有效削弱高阶项误差。但由于高阶系数中高阶高次系数的误差大于高阶低次系数的误差，GRACE噪声并不是各向同性的，所以仅针对阶数进行滤波的各向同性高斯滤波是不够准确的。随后研究者提出了针对阶次分别进行高斯滤波的非各向同性滤波算法^{[21, 26]}。

Han^[23]针对不同方向选取了不同的高斯滤波半径，在纬度方向上有更高的分辨率，建立了一个同时针对阶数和次数进行高斯滤波的滤波器；Zhang^[28]的扇形滤波法与Han的方法思路类似，也同时对阶次进行高斯滤波，但使用的高斯滤波核函数不同。相比高斯滤波，这两种非各向同性滤波方法都能更有效地削弱南北向的条带误差，且同时有效保留时变重力信号，处理之后的重力场与其他地球物理模型更加一致。

对于南北向条带误差，基于高斯滤波的方式效果有限，尤其在平滑半径较小时。由于条带误差的存在反映了球谐系数中存在系统性的相关误差，所以直接消除模型系数之间的相关误差可以更好地去除条带^{[18, 27, 28]}。这类算法在实际中通常和高斯滤波结合使用。

Swenson和Wahr^[18]首先提出了一种去相关滤波的方法。该方法首先平滑任意给定次数m的球谐系数，当固定重力场模型系数的次数m时，对于偶数阶系数之间和奇数阶系数之间存在的相关误差，分别使用一个滑动窗口多项式最小二乘拟合法来消除。根据经验对这个滑动窗口滤波器设置不同的窗口宽度和次数参数来进行滤波。随后，Chambers^[28]和Chen^[29]基于类似的方法保留部分低阶次球谐系数并用多项式拟合高阶次球谐系数以去除系统误差。Duan^[30]比较了以上三种去相关算法，并与他本人提出的改进后的去相关算法对比；由于改进的算法中窗函数与阶数次数都相关，从全球重力场数据来看，较好地平衡了削弱条带误差，保留了地球物理信号。

去相关算法是GRACE后处理过程中应用非常广泛的一种算法，在大地水准面变化^[31]、海水质量变化和冰川均衡调整^[32]等研究中均被用于减少误差，起到了良好作用；在南极冰盖质量变化^{[33, 34]}和河流流域水储量变化^{[25, 37]}等方面也得到了广泛应用。

此外，去相关算法处理过的GRACE数据中能探测到特大地震的同震和震后重力变化，比如2004年苏门答腊地震^{[29, 38]}，2010年智利地震^[39]等。Heki和Matsuo在研究2010年智利地震^[40]和2011年日本东北大地震^[41]时都结合扇形滤波和去相关方法处理GRACE数据，通过对震中附近某一点的重力时序分析可以明显观察到同震重力变化。

詹金刚^[42]认为这种算法在距离赤道很远的地方内效果显著，但是赤道两侧区域效果提升并不十分理想；他在Wahr的基础上改进了去相关方法^[42]，针对滑动窗的特点使用反向边界延拓技术；改进后的方法在赤道两侧区域也有不错的去条带误差效果。

GRACE产品及其地球物理模型都是以空间地图的时序为表现形式，在使用这些数据时图上会被某种特定的噪声图案占据。所以需要识别出这些主要图案，并且移除不必要的小尺度噪声信号。主成分分析正是这样一种处理方式。

主成分分析法作为一种数据压缩方法可以用于数据重建，能找出数据集中相对较小，同时尽可能减少冗余的独立模数，并且能客观地分析出数据变化规律。因为此方法不需要先假设数据呈周期性变化，也被称为经验正交函数法(emperical orthogonal functions，EOF)。经验正交函数可由一组格网值的时间序列确定，将其信号分离为空间模态和时间系数。EOF法既可以在空间上平滑，也可以在时间序列上拟合。

Schrama和Wouters^{[43, 44]}研究了如何应用EOF法处理GRACE月数据，并且给出了比较合适的参数配置。处理之后去除条带效应的效果较为明显。

相对于空间上的滤波，EOF法主要用于处理时间序列，比较适合发掘GRACE时变重力场的季度^[45]、年度变化规律^{[46, 47]}。GRACE数据也能客观反映这种季节性变化的分布规律，并提取出地球物理信号^{[48, 49]}。

由于超导重力仪也是连续观测的^[50]，有研究者将超导重力仪的时序观测资料和GRACE时变重力场数据结合起来分析^{[51, 52]}，EOF法是比较两者特征和效果的关键。不过，关于联合这两种数据的可行性，Crossley^[53]和Camp^[52]得到了截然相反的结论。

小波分解法，也叫小波多尺度分析法，是一种位场分离的方法，可以将重力场数据分解到不同的尺度空间中，不同的尺度大小反映不同对象引起的重力变化，这种思路正好适用于研究局部地区的重力变化。

Panet等^[54]在利用GRACE月度重力数据分析2004年12月和2005年3月的苏门达腊地区地震的质量迁移过程时，利用小波分解法处理GRACE数据得到了水准面的变化，结合地壳和上地幔的密度变化以及垂直形变的研究分析了震后信号的意义。之后他又使用了另一种连续小波分解法^[55]，有效地分离出苏门答腊地震的震后重力变化信号。Mikhailo等^{[56, 57]}利用小波分解法处理GRACE数据，分析比较了近10 a的三个特大地震的同震以及震后重力变化时空分布。

Fengler等^[58]人利用三次多项式生成的小波分析GRACE数据，研究了地球上不同地区的时变重力位，可发现时序中存在明显的季度变化，并推测这种变化是由大尺度地下水质量变化导致。陈国雄等^[59]利用GRACE卫星数据分析研究了汶川及其周边的重力场变化情况。由于地下水、构造形变以及深部物质流动在空间尺度上具有差异，表现为不同波长成分。因此，利用小波多尺度分析方法能够得到不同层次的扰动或者地壳地幔密度变化，能反映不同深度处的质量迁移。

很多时候研究区域并不是整个地球而是局部地区，而GRACE重力数据是包含全球的重力场信息,由全球球谐系数构成的信号在整个球面上正交，将这组球谐系数直接用在局部地区不够理想，因为只有球面上一小部分的信号被考虑进来，在局部不满足正交性，从而给信号的提取带来了困难。利用去相关算法等处理方法不仅会削弱误差信号，也会削弱地震水文等变化引起的信号，所以如何有效提高目标区域的信噪比成为一个难题。

在数学上，空间局部化和频谱局部化无法兼顾，不过在这两者中间可以找到最优。如果在球面上设置一组正交函数,将信号按照需求最大化也集中在局部区域和局部频域,这样投影在这组正交函数上的信号在局部地区能量最大,泄漏到外部的能量最小。我们称这种基函数为Slepian基函数^[60]，这种处理方法为Slepian方法，可以同时在空间上和频域上得到最优。

Simons等^[60]将Slepian方法引入到地球物理领域中，随后他应用这种方法估算地球表面位场^[61]，并发现在研究冰川质量变化^[62]和冰盖质量变化^[63]时，该方法能凸显一些较细微的重力变化，因此较为适合GRACE数据，所以也被Han^{[64, 65, 66]}、Wang^{[67, 68]}、Dai^[69]和Camobitti^[70]等应用于地震的同震重力变化研究中。利用GRACE观测的同震重力变化图形和正演模型很类似，并且能一定程度上反演震源参数，证明了该方法在监测地震活动中的有效性。不过GRACE重力数据对于地壳中发生的逆冲地震的约束作用明显强于更深的比如上地幔处的逆中地震^[64]。

基于先验信息的后处理方法需要对特定问题做特定处理，比如设置陆地边界，建立误差估计模型等，并且往往需要确定噪声信号和目标信号。这类方法通常需要建立误差协方差矩阵，但GRACE产品并不直接提供这类数据，需要自行建模，这样的GRACE重力场可能需要地面模型进行约束。除了使用前述滤波方法减少南北向的条带误差，也可以通过添加先验信息并正规化的方式稳定这个病态问题^{[24, 71, 72, 73]}。

由于在削弱误差的同时不可避免地会削弱信号本身，所以正规化矩阵的选择必须注意，要在误差消除和信号衰减中取得平衡。Save等^[74]确定正规化矩阵的方法与其他同类研究不同，没有使用信号特征作为先验信息，而是使用GRACE数据解的误差特征。根据地面覆盖情况将月数据分为两类，根据不同的误差特征分别设计不同的正规化矩阵。处理后的GRACE重力场的误差有了明显减少，南北向条带误差也被削弱，同时适配后的残差分析表明通过这种处理方法损失的地球物理信号小于GRACE噪声信号的数量级，也就是说正规化产生的信号损失是可以忽略的。

部分基于先验信息的算法和无需先验信息的算法有类似性质，只不过前者的平滑核函数考虑了误差协方差矩阵，比如Kusche等^[75]讨论的另一种去相关算法以及非均衡的滤波器；Davis等^[76]设计的一种基于统计学的滤波器，在空间域该滤波器和高斯平滑有类似的性质；Klees等^[77]进一步优化基于经验正交函数法的滤波器等。实验表明上述后处理算法均能更有效地去除南北向条带误差，并且减少信号损失，比不考虑先验信息的滤波器表现更好，更能适应每一种误差情况。

总的来说，在考虑每组数据具体的误差情况时，设计的后处理方法能够尽可能多地消除噪声信号的同时保留信号。

经过十多年的研究，GRACE数据的后处理方法出现了很多不同的发展方向，以上归纳的算法在数据处理思路上有很大区别，在处理对象上也有不同。高斯滤波和非各向同性滤波方法，使用时不必考虑不同月份的数据差异；位系数去相关法、小波分解法和Slepian方法都是针对单独月份的重力场数据处理计算的滤波方法，使用时不考虑不同月份的数据差异；而基于先验信息的处理方法有时则需要考虑误差与时间的相关性，因而需要同时对整个时序进行处理^[76]。主成分分析法作为一种时空处理手段，也可以整体对时序的重力场球谐系数进行处理^[44]。

不管是否考虑信号和噪声的特定结构，都可以得到精度相当的结果。不需要先验信息的方法相对考虑先验信息的方法来说计算量小，容易评估效果，有较强的普适性。考虑先验信息的方法很复杂，误差协方差矩阵不仅需要自行建模且计算量也非常大，计算过程中也需要很高超的技巧，如果误差模型与时间相关在时序变长的情况下计算还将进一步复杂化。

滤波方法通常会在压制噪声信号的同时对所需要的目标信号产生影响，不利于提取地球物理变化目标。在某些情况下，滤波导致的信号损失会非常严重，比如Velicogna和Wahr^[78]估计利用最优平均核函数滤波法在计算格陵兰的质量平衡时可能会有50%的信号损失，Swenson和Wahr^[79]也发现在使用去相关方法研究里海的质量变化时，变化的振幅有60%的衰减。这种信号衰减的问题如果通过引入误差协方差矩阵，振幅衰减会减弱。

所以，当外部信息足够准确可靠的时候，根据先验信息处理的算法更好。对于无需经验信息的滤波方法，也应该灵活应用发挥其快速方便的优点。未来更有效的算法应该结合两类算法的优点，一方面利用普适的滤波器进行基本处理；另一方面针对具体研究区域增加先验信息的约束，比如Slepian方法本身可以在有效提高信噪比的同时保留目标信号，如果再考虑先验信息也许可以进一步提高滤波效果。另外，对于目标信号容易被淹没在背景噪声中的情况，比如研究同震或震后重力变化时，则应尽量减少使用会削弱地震信号的滤波算法，而应考虑先验信息以保留更多的有效信号。