南极冰盖是地球表面最大的大陆冰川，蕴涵有全球70%的淡水资源，这里的许多冰川厚度超过4 000 m，平均厚度约为2 450 m。当今全球环境变化越来越严重，南极的冰雪消融导致全球海平面升高，从而造成海岸、河口、海湾自然生态环境失衡，给海岸带生态环境系统带来灾难，因此研究南极地区冰盖质量变化成为了南极地区科学研究的一个重大课题。

目前，研究南极冰盖质量变化的方法主要有质量平衡法、高程测量法、重力测量法^[1]。然而这3种方法都有着各自的局限性。质量平衡法无法大面积精确确定冰雪的积累量和融化量，它适用于小范围内的质量变化研究；卫星测高法研究质量变化会有较大的坡度误差，在边缘区域和坡度大的区域精度不高^[2]，且冰雪密度变化引起的体积变化并不一定和冰雪质量变化有关^[3]；重力测量法确定的质量变化精度较高，但空间分辨率较低，且受到冰川均衡调整模型误差影响^[4]。

2002年，GRACE卫星的发射极大地改善了利用时变重力场模型来研究地球物理机制。GRACE卫星由两颗相同的卫星组成，两者相距约200 km，在同一近地轨道 (轨道高度约300~500 km) 上飞行，轨道倾角89.5°，通过卫星上的微波测距系统精密测量两星之间的距离随时间的变化。GRACE通过反演高精度的月重力场模型，直接测定南极冰盖总质量变化。近年来，许多学者利用GRACE在南极冰盖质量变化方面做了很多研究。Velicogna等^[5]采用CSR发布的GRACE RL01数据估算出2002-04-2005-08南极的冰盖质量消融速率为 (152±80) km³/a; Chen等^[6]基于去相关滤波和高斯滤波采用CSR发布的GRACE RL04数据估算出2002-04-2006-09南极冰盖质量大致处于平衡状态；罗志才等^[7]利用CSR发布的GRACE RL04数据，采用去相关滤波P3M6结合300 km Fan滤波计算出南极冰盖呈现加速消融趋势，冰盖消融速率为78.3 Gt/a。本文采用CSR最新公布的GRACE RL05数据，计算2002-2015年近13年的南极冰盖质量变化情况，并比较分析了几种滤波算法和GIA模型对反演南极冰盖质量变化的影响。

1 数据与计算方法 1.1 数据准备

目前，有多家数据中心能够提供GRACE大地水准面球谐系数，包括CSR、GFZ (Geo Forscheungs Zentrum Posdam)、JPL (Jet Propulsion Labratory)、CNES/GRGS (Centre National d’ Etudes Spatiales/Groupe de Recherches) 等。本文采用的数据为CSR公布的2002-07-2015-04期间的RL05版本GRACE月重力数据，数据形式为月系数解。

这几家数据中心所发布的GRACE二级数据中，重力场高阶次的球谐系数均存在较大的误差，在空间域表现为相当严重的条带噪声，这主要是由混频误差和相关误差造成的。由于混频误差和相关误差的存在，即使使用现有的潮汐模型、海洋模型和大气模型会减小这些误差，但仍不能完全消除，因此在采用GRACE数据处理过程中，需要采用合适的滤波算法来消除或削弱这种误差的影响。

1.2 数据计算方法 1.2.1 球谐系数计算等效水量

依据Wahr等 (1998) 的理论，GRACE时变重力场反演地球表面任意一点的面密度变化可用以下公式表示^[8]：

$ \begin{array}{l} {\Delta _\sigma }\left( {\theta, \varphi } \right) = \frac{{a{\rho _{{\rm{ave}}}}}}{3}\mathop \sum \limits_{l = 0}^ \propto \mathop \sum \limits_{m = 0}^l {{\tilde P}_{lm}}\left( {{\rm{cos}}\theta } \right) \bullet \\ \frac{{2l + 1}}{{1 + {k_l}}} \times \left[{\Delta {C_{lm}}{\rm{cos}}\left( {m\varphi } \right) + \Delta {S_{lm}}\sin \left( {m\varphi } \right)} \right] \end{array} $

(1)

式中，θ和φ分别表示地心余纬和余经；a为地球平均半径；ρ_ave为地球平均密度；$ {\tilde P_{lm}} $为正则化的连带勒让德级数；k_l为负荷勒夫数；ΔC_lm和ΔS_lm为时变重力场模型的球谐系数。

质量变化通常以等效水量Δ_h(θ, φ) 表示：Δ_h(θ, φ)=Δ_σ(θ, φ)/ρ_w, ρ_w为水的密度。

1.2.2 滤波方法

由于GRACE重力场模型系数误差会随着阶数的增大而迅速增大，会导致结果中有非常明显的条带噪声，为了消除或减弱这种条带噪声，需要采用合适的滤波方法进行处理。

Wahr等 (1998) 提出用高斯滤波法进行平滑。高斯滤波平均核函数球谐系数W_lm^C和W_lm^S与对应的平均核函数球谐系数和的对应关系为：

$\left\{ {\begin{array}{*{20}{c}} {W_{lm}^C}\\ {W_{lm}^S} \end{array}} \right\} = {W_l}\left\{ {\begin{array}{*{20}{c}} {{\vartheta _{lm}}^C}\\ {{\vartheta _{lm}}^S} \end{array}} \right\}$

(2)

W_l由以下递推关系得到：

$ \left\{ \begin{array}{l} {W_0} = 1\\ {W_1} = \frac{{1 + {e^{-2b}}}}{{1-{e^{-2b}}}} - \frac{1}{b}\\ {W_{l + 1}} = - \frac{{2l + 1}}{b}{W_l} + {W_{l - 1}} \end{array} \right. $

(3)

式中，$b = \frac{{{\rm{ln}}2}}{{1 - {\rm{cos}}\left( {r/a} \right)}}$，其中r为高斯滤波半径；W_l表示与时变重力场球谐系数阶数相关的高斯滤波函数。

还有学者提出了扇形滤波^{[8, 9]}，它是一种非各向同性的高斯滤波，对重力场模型球谐系数的阶数和次数同时进行高斯平滑处理，对应的等效水量表示为：

$ \begin{array}{l} {\Delta _\sigma }\left( {\theta, \varphi } \right) = \frac{{a{\rho _{{\rm{ave}}}}}}{{3{\rho _w}}}\mathop \sum \limits_{l = 0}^ \propto \mathop \sum \limits_{m = 0}^l {W_l}{{\tilde P}_{lm}}\left( {{\rm{cos}}\theta } \right) \bullet \\ \frac{{2l + 1}}{{1 + {k_l}}} \times {W_m}\left[{\Delta {C_{lm}}\cos \left( {m\varphi } \right) + \Delta {S_{lm}}\sin \left( {m\varphi } \right)} \right] \end{array} $

(4)

式中，W_l和W_m分别表示与时变重力场模型球谐系数阶数和次数相关的高斯滤波函数。

Swenson提出了去条带滤波^[10]，由于球谐系数相同次偶数阶与奇数阶之间是分别相关的，因此去条带滤波的基本思路是从m次开始分别采用n次多项式来消除球谐系数间的相关性。

2 计算结果分析 2.1 不同滤波算法比较

如前所述，由于混频误差和相关误差的存在，GRACE二级数据中，重力场模型高阶次的球谐系数均存在较大的误差，在空间域表现为非常明显的条带噪声，为了消除或削弱这种条带噪声，需要采用合适的滤波方法进行处理。本文分析比较了滤波半径为300 km的高斯平滑滤波、扇形平滑滤波以及去条带扇形滤波方法，计算结果的效果图如图 1所示。

从图 1可以看出，未经过滤波处理时，等效水量在全球分布存在明显的条带现象，无法从中看出主要质量变化区域，高斯滤波能很大程度地消除高频误差，提高结果精度，而扇形滤波是在高斯滤波基础上进一步滤波，因此其效果更优，去条带扇形滤波是在考虑去除条带现象的基础上进行扇形滤波，因此该滤波精度更高。从图 1的效果图可以看出，这3种滤波算法都能够有效地削弱条带噪声，高斯滤波能一定程度压制去条带现象，扇形滤波比高斯滤波效果更好，而去条带扇形滤波算法消除了去条带现象，获得的结果效果更佳。

2.2 不同GIA模型差异

GRACE数据反演得到的是各种地球物理过程引起的总质量变化，扣除大气、海洋、陆地冰、冰后回弹等影响，剩下的才是真正的冰盖质量变化。已有研究表明，GIA模型是影响GRACE反演南极冰盖冰雪质量变化的关键因素，其误差量与南极冰盖冰雪质量变化量级相当。本文收集了目前最新的3种GIA模型 (Geruo13、IJ05_R2和W12a)，这3大模型反演时都将GRACE观测量作为模型约束。首先由GIA模型计算得到对应的球谐系数变化，然后采用和GRACE相同的数据处理方法得到相应模型的等效水量时空变化分布，从而可以比较分析不同GIA模型对于南极冰盖质量变化的贡献，结果见图 2。

从图 2发现，不同GIA模型之间反演得到的南极大陆质量变化存在着很大差异，尤其在坎普冰架、龙尼冰架和东南极区域差异明显。在罗斯冰架区域，Geruo13模型和W12a模型反演结果分别显示出该地区有着50 mm/a和30 mm/a的质量积累，而IJ05_R2模型反演结果显示出该地区的质量积累率为8 mm/a；在龙尼冰架地区，Geruo13模型和IJ05_R2模型结果显示该地区质量处于大致均衡状态，而W12a模型的结果则显示龙尼冰架区域处于近20 mm/a速率的质量削减状态；在整个东南极地区，Geruo13模型和IJ05_R2模型均显示轻微的质量减少，而W12a模型则显示出非常明显的质量减少。在这些区域采用不同GIA模型将导致反演的冰雪质量变化的误差较大。

2.3 不同时段南极冰盖质量变化分析

本文采用的由CSR公布的RL05版本GRACE月重力数据的时间跨度为2002-07-2015-04。为了研究不同时间段南极冰盖质量变化，根据恩德比地和毛德皇后地的冰盖质量变化时间序列，分成了2002-07-2007-12、2008-01-2011-12、2012-01-2015-04共3个时间段，然后分别计算这3个时间段和整个时间跨度南极冰盖质量变化，利用球谐系数得到南极冰盖1°×1°的等效水量时间序列，对该结果进行分析，得到了南极冰盖质量变化的趋势项，结果如图 3所示。

从图 3可以明显看出，整体上南极冰盖质量变化呈现出非常明显的地区不平衡性，西南极 (区域A) 冰盖的质量消融最为严重，年变化速率达到了80 mm/a，南极半岛 (区域B) 也存在比较明显的质量消融现象。而冰盖质量积累主要出现在恩德比地 (区域D)、毛德皇后地 (区域E) 和坎普冰流 (区域C)。已有研究主要关注在西南极、南极半岛和恩德比地区，本文研究发现，毛德皇后地和坎普冰流也出现了明显的质量变化。

比较不同时间段南极冰盖质量变化情况发现，虽然整体上不同时间段南极地区冰盖质量变化趋势是一致的，但是各区域的质量变化的量级却是不同的，尤其在恩德比地和毛德皇后地表现得非常明显。在上述3个时间段内，恩德比地在2008-01-2011-12期间的质量积累最为明显，冰盖质量积累速率达到了近80 mm/a，而在2002-07-2007-12和2012-01-2015-04，该地区的质量积累速率约为30 mm/a。而毛德皇后地在2002-07-2007-12冰盖质量呈均衡状态，在2008-01-2011-12和2012-01-2015-04冰盖质量有着40 mm/a的积累速率。

为了进一步研究其详细特征，计算了这5个区域的等效水量变化的时间序列，结果如图 4所示。从图 4中折线的斜率变化特征可以看出，不同时间段南极地区冰盖质量变化有较大差异。

事实上，Geruo13、IJ05_R2和W12a这3种GIA模型在建立过程中均采用了GRACE数据对其结果进行约束^[11], 但由于不同时间段GRACE计算结果存在差异，因此在采用GIA模型对GRACE数据进行改正时，需要考虑不同时间段的影响。

为了获取冰盖冰雪质量变化，需要从冰雪质量变化中扣除冰后回弹的影响。§2.2对3个不同模型进行了比较，发现IJ05_R2在冰架的信号最小，而实际上冰架总体处于平衡状态，因此选用IJ05_R2作为冰后回弹模型，计算得到了冰盖冰雪质量变化，结果如图 5所示。可以看出，西南极A区域在冰后回弹改正前后变化不明显，而在其他区域的冰后回弹影响在图中均有不同程度的显示。

3 结束语

1) 高斯平滑滤波、扇形平滑滤波、去相关滤波以及组合滤波都能有效地消除条带噪声，但不同滤波算法效果有所差别，组合滤波算法效果更好。

2) 南极冰盖质量变化呈现非常明显的地区不平衡性，除了已有研究发现的西南极、南极半岛、恩德比地这3大区域存在很大差异外，还发现坎普冰流和毛德皇后地存在很大差异，且不同时间段质量变化也存在很大差异。

3) 利用最新的3个GIA模型分析了其对南极冰盖质量变化的贡献，发现不同模型之间差异较大，主要体现在坎普冰流、龙尼冰架和东南极等区域。由于在建立这3个模型时均利用了GRACE数据，而不同时间段GRACE计算结果存在很大差别，因此采用GIA模型来对GRACE数据进行改正时，需要考虑不同时间段的影响，而反演GIA模型时，需要考虑更长时间的GRACE观测结果进行约束。