连续在轨运行超过15 a的GRACE重力卫星为全球地表质量迁移的研究提供了丰富的资料^[1-3]。利用GRACE位系数法研究地表质量迁移的过程中主要需考虑3个问题：1)GRACE低阶项的改进，包括1阶项改正^[4]和2阶项替换^[5]；2)GRACE卫星编队特点导致其位系数中存在高频误差和相关性误差，由此发展出空间平滑滤波^[6]和去相关滤波^[7]；3)如何恢复由位系数的截断和空间滤波所导致的信号泄露和损失。

近几年来，如何对时变信号进行有效恢复逐渐成为国内外研究的热点，如尺度因子法^[8-9]、迭代恢复法^[10-11]和频域恢复法^[12]等，但不同的信号恢复方法的效果均存在明显差异。

本文采用CSR RL05数据，基于尺度因子法和迭代恢复法两种恢复方法，选取4个典型区域，从时间序列变化和空间分布两个角度对两种信号恢复方法的效果进行分析，并选用CSR Mascon产品进行对比。

1 数据处理方法 1.1 时变信号损失的成因及影响

利用位系数变化来获取地表质量变化(等效为水高)的基本原理为^[6]：

$ \begin{array}{l} \;\;\;\;\mathit{\Delta H}\left( {\theta, \lambda } \right) = \frac{{a{\rho _{{\rm{ave}}}}}}{{3{\rho _{\rm{w}}}}}\mathop \sum \limits_{l = 0}^{{l_{{\rm{max}}}}} \mathop \sum \limits_{m = 0}^l {{\bar P}_{lm}}\left( {{\rm{cos}}\theta } \right) \cdot \\ \frac{{2l + 1}}{{1 + {k_l}}}{W_{lm}}\left( {\mathit{\Delta }{{\bar C}_{lm}}{\rm{cos}}\left( {m\lambda } \right) + \mathit{\Delta }{{\bar S}_{lm}}{\rm{sin}}\left( {m\lambda } \right)} \right) \end{array} $

(1)

式中，θ、λ分别为地心余纬和经度，a为地球平均半径，ρ_ave为地球平均密度，ρ_w为水的密度，l、m分别为地球重力场位系数的阶、次，P_lm为规格化勒让德函数，k_l为负荷勒夫数，{ΔC_lm，ΔS_lm}为位系数的变化值，W_lm为空间平滑滤波。

在实际应用中，往往将位系数的最大阶数截断至60阶(60阶所对应的波长约为330 km^[9])，这导致所计算出的空间分辨率为1°×1°的格网(纬度相差1°约为111 km)中会缺失小于300 km的时变信号，从而造成截断误差。

为消除GRACE月重力场模型位系数中存在的高频误差，式(1)中还需加入空间平滑滤波W_lm。空间平滑滤波的实质是对不同阶次的位系数定权，当空间平滑为高斯滤波时，W_lm=W_l，即位系数所对应的权值大小仅与阶数l有关。从图 1可以看出，平滑半径一定时，随着阶数的增大其所对应的权值W_l逐渐减小，进而降低了来自高阶高次项对信号的贡献。阶数一定时，平滑半径越大，其所对应的权值就越低。当采用平滑半径为500 km的高斯平滑时，50阶以上的位系数所对应的权值便趋近于0，一定程度上加深了截断误差的影响。

另一方面，阶数相同的所有位系数均采用相同的权重W_l，相当于利用圆形高斯滤波器进行空间卷积，平滑半径r对应于卷积运算符下降到其峰值一半时的距离^[6]。这便会造成空间一点经平滑后信号向周围区域泄露，进而造成泄露误差。

分别在178.5°E~178.5°W之间模拟等效水柱高为30 cm的3个质量块，分别位于13.5°~16.5°N(低纬)、43.5°~46.5°N(中纬)、73.5°~76.5°N(高纬)，如图 2(a)所示；将这3个水层重力信号转换为截断至60阶的位系数后，如图 2(b)所示；再分别加入300 km和500 km的高斯平滑滤波后，分别如图 2(c)、(d)所示；沿179.5°W经线水柱高横截面值如图 2(e)所示。可以看出，经位系数截断后，信号明显削弱，且纬度越高，削弱程度越大。使用空间平滑滤波后，信号的削弱程度进一步增强，且平滑半径越大，信号损失越强，同时也出现信号向外界泄露的情况。因此采用合理、有效的反演方法来恢复GRACE数据解算过程中的信号损失，将会进一步提高GRACE计算结果的可靠性与准确性。

1.2 尺度因子法

尺度因子法最初由Velicogna等^[8]提出，Landerer等^[9]将其进一步改进，提出基于先验陆地地表模型(land system model，LSM)中的信号分布来确定尺度因子的方法。基本原理可以概括为：对于特定的先验模型，其任意格网点(i, j)在t时刻的等效水高变化为Δh_ij(t)，经过与GRACE相同处理方法后得到的等效水高变化为Δ$\tilde h$_ij(t)，那么在所选的时间范围内，存在尺度因子k_ij，使得其满足：

$ \mathop \sum \limits_i^N {\left( {\mathit{\Delta }{h_{ij}}\left( t \right) - {k_{ij}}\mathit{\Delta }{{\tilde h}_{ij}}\left( t \right)} \right)^2} = {\rm{min}} $

(2)

式中，N为时间序列的个数。将求得的k_ij应用到GRACE计算的时变信号中，即可获得改正后的时变信号。

基于全球陆地资料同化系统(global land data assimilation system，GLDAS)中NOAH、Mosaic、VIC和CLM四种陆地地表模型所获得的全球1°×1°尺度因子如图 3所示，所选时间段为2003-01~2013-12，采用500 km的高斯平滑。可以看出，基于不同地表模型所获取的尺度因子存在一定的差异。

1.3 迭代恢复法

迭代恢复法又称前进模型法(forward modeling method)，最初由Chen等^[10]提出，用于对局部区域的信号泄露进行改正；之后Chen等^[11]将其改进，应用区域扩展至全球。该方法的基本思想是通过迭代调整使得模拟模型逐渐逼近观测值来恢复损失信号。迭代过程在所有月份中重复进行，最终获取2003-01~2013-12时间段内共计124个月的质量恢复结果。

1.4 GRACE数据处理方法

采用CSR提供的GRACE RL05数据，数据时间段为2003-01~2013-12，共计124个月(部分月份数据缺失)，只计算前60阶的位系数。采用Swenon等^[4]计算的地心改正项加回一阶项，采用卫星激光测距(SLR)测得的C₂₀项对GRACE原始的C₂₀项进行替换。空间平滑滤波选择平滑半径为500 km的高斯平滑。为与CSR Mascon产品保持一致，本文选用Geruo13模型^[13]进行GIA改正。

选取CSR Mascon产品^[14]作为对比，该产品直接利用GRACE任务的卫星跟踪卫星技术解算全球等效水高变化，能有效克服GRACE位系数法中存在的高频误差和南北条带问题^[15]，解算结果更为接近真实时变信号。其空间分辨率为0.5°×0.5°，为保持一致，本文将其插值为1°×1°。

在时间序列拟合过程中，线性项仅考虑趋势项，周期项考虑周年项和半周年项。下文获取的年际变化为移除周期项并进行0.5 a滑动平均后的结果。

2 时变信号恢复效果对比分析 2.1 不同恢复方法效果的时间序列分析对比

不同信号恢复方法获得的区域水储量年际变化时间序列、长期趋势和周年振幅如图 4所示，图中NOAH、Mosaic、VIC和CLM为尺度因子法所选用的不同地表模型。可以看出，在进行质量恢复前，由GRACE数据直接计算获得的区域水储量长期趋势、周年振幅均处于较低水平，表明信号损失较为明显。进行质量恢复后，长期趋势、周期振幅均有所增强，损失信号在一定程度上有所恢复。但在不同的研究区域，不同恢复方法对信号的恢复程度具有明显差异。

在长江流域(图 4(a))，尺度因子法和迭代恢复法均较好地恢复了其水储量变化的长期趋势，迭代恢复法的结果(3.68±1.48 mm/a)更接近于Mascon的反演结果(3.97±1.00 mm/a)。尺度因子法对长江流域的质量恢复主要体现在周年振幅上，但采用不同的地表模型，得到的恢复效果也有明显差异。同样，在周年振幅上迭代恢复法的恢复效果，更接近于Mascon的结果。

恒河流域位于印度北部，人口密集，地下水大量开采致使该区域成为世界地下水亏损较为严重的地区之一^[2]。由于GLDAS水文模型缺少地下水信息，因此利用4种地表模型获得的尺度因子均不能较好地恢复恒河流域水储量亏损的趋势(图 4(b))。迭代恢复法恢复的长期趋势为-18.24±2.80 mm/a，高于Mascon的结果(-14.86±2.34 mm/a)。在周年振幅上，两种恢复方法获得的结果相当，均高于Mascon的结果。

亚马孙流域常年降水丰富，水储量变化具有明显的周期性。从图 4(c)中可以看出，信号恢复前后的亚马孙流域水储量年际变化并没有明显的差异，这表明，在GRACE数据处理中亚马孙流域的信号损失主要体现为周期信号的削弱。在长期趋势的恢复上，尺度因子法与迭代恢复法恢复结果基本相当，且均接近于Mascon的结果。而在周年振幅的恢复上尺度因子法恢复的信号强度要高于迭代恢复法的效果，且4种地表模型恢复的效果基本相当。

巴拉那流域位于南美洲南部，是南美洲第二大流域，可分为上、下两段。上段位于巴西南部，降水充沛，水储量丰富。下段位于阿根廷拉普拉塔平原，该部分靠近南美洲巴塔哥尼亚冰川区域^[12]，因此在进行质量恢复前，大量的冰川消融信号很可能会泄露进下段区域，影响对该流域整体的水储量变化分析。从图 4(d)中可以看出，在信号恢复前巴拉那流域水储量呈现出下降趋势(-3.52±1.06 mm/a)，这表明其附近的冰川消融信号可能泄露到了该区域。经过迭代恢复改正后，巴拉那流域水储量恢复上升趋势，约为2.54±2.38 mm/a，接近2.13±1.89 mm/a的Mascon结果。而由于地表模型缺乏冰川消融信息，尺度因子法对巴拉那流域的长期趋势恢复效果较差，恢复效果仍体现为对其周期项的恢复。

2.2 不同恢复方法效果的空间分布分析对比

选取§2.1中4个流域所处的亚洲、南美地区来对比研究不同恢复方法所获得的水储量变化空间分布效果，其中长江流域、恒河流域位于亚洲地区，亚马孙流域、巴拉那流域位于南美地区，以上4个流域的范围均在图 5、图 6中用虚线标示。

图 5给出亚洲地区包括多种水储量变化典型区域的相关研究，如华北地区、印度北部地下水长期匮乏的研究，长江流域水储量的研究以及青藏高原质量变化研究。经尺度因子法恢复后，图 5(b)~图 5(e)中水储量长期变化趋势信号得到一定程度的增强，信号的泄露现象也得到了有效的改善，如长江中下游地区、恒河流域等。但由于所采用的先验地表模型缺乏地下水信息，其在华北、印度北部的长期亏损信号未得到有效恢复，同时青藏高原内部明显的质量增长信号也未得到有效恢复。迭代恢复法则有效地恢复了印度北部、华北地区长期的亏损信号，青藏高原内部的增长信号也得到了有效的恢复，其信号恢复的空间分布基本与CSR Mascon相当，但在部分区域可能存在信号的过度恢复，如中南半岛地区等。

图 6中的亚马孙流域，尺度因子法恢复的效果在空间分布上与CSR Mascon相当，展现出较好的可视性，且4种地表模型获得的结果之间的差异较小，这表明当地表模型中的水文信息与真实水储量变化能较好地吻合时，尺度因子法能够很好地恢复信号。但同样由于冰川信息的缺乏，尺度因子法在南美洲南部地区的质量恢复效果较差。迭代恢复法则很好地恢复了南部地区冰川消融信号，与CSR Mascon结果吻合较好，但同样在部分地区也存在过量恢复的情况。

3 结语

尺度因子法易于实现，但其获得的尺度因子完全依赖于先验水文地表模型的选择。由于目前大多数地表模型缺乏地下水、冰川消融信息，因此该方法并不适用于恢复地下水亏损以及冰川消融地区。总体来讲，该方法主要用于恢复区域水储量变化的周期项，长期趋势恢复效果不明显。

迭代恢复法恢复过程相对繁琐，但在信号恢复过程中仅受GRACE原始观测值的约束，不需任何先验信息，能够有效地恢复原始信号中损失的长期趋势信号和周期信号。但该方法在部分地区可能存在信号过量恢复的现象。总体而言，该方法对时变信号恢复的效果较优，建议采用该方法进行时变信号恢复。

致谢: 感谢CSR提供GRACE Level 2、Level 3数据集，感谢GES DISC提供GLDAS数据集，感谢武汉大学测绘学院计算中心提供计算支持。