地球表面的质量迁移分为流体质量迁移和固体质量迁移两类.由于地球不是严格的刚体，质量季节性迁移引起的重新分布将会导致地表的季节性负荷响应、产生形变.其中流体质量季节性迁移包括了大气、海洋、陆地水（地表水、土壤水、冰雪等）各部分的综合影响，主要受气候因素驱动^{[1, 2]}.

近年来，地面点的季节性垂直形变多是通过GPS长期连续的观测获得.GPS连续观测的季节性垂直形变与大气、海洋、陆地水等地球物理因素驱动有关^{[3, 4]}，目前用陆地水等地球物理模型综合模拟地表形变受到模型精度不确定性的限制；另一方面，GPS站点的分布有限、分散，获取的单点地面垂直形变结果不能满足空间连续的形变监测需求，甚至受环境或经济因素干扰在很多区域缺乏GPS连续观测站，如图 1中恒河-澜沧江流域.

图 1

Fig. 1

图 1 背景为中国及邻区的GRACE振幅 白色圆圈表示中国国内的GPS站，白色方块表示邻区的GPS站，蓝色箭头的长度和方向分别表示站上GRACE周年垂直形变的振幅和初相位，红色箭头的长度和方向分别表示站上GPS周年垂直形变的振幅和初相位. Fig. 1 Climate-driven annual vertical deformation amplitude of the solid Earth calculated from GRACE using De-stripes in China and around The circles and squares filled with White color represent GPS sites in China and around China respectively. Amplitudes and initial phases of annual vertical deformation from GRACE as blue arrows, GPS as red arrows.

随着卫星重力探测技术发展，特别是2002年3月17日GRACE（Gravity Recoveryand Climate Experiment）卫星发射后，重力卫星观测质量重新分布的中、长波信号达到了前所未有的高精度，其动态和空间连续结果对地表负荷响应的形变监测无疑是一种新技术手段---它的监测结果包含了各种地表负荷响应，与各种地球物理模型综合模拟的地形变相比，避免了模型的精度不确定性带来的误差.2004年，Davis等^[5]利用GRACE2002年4月至2004年4月最初两年的数据基于“高斯平滑”滤波方法首次解算了南美洲亚马逊河流域的地表周年垂直形变，并与流域及附近的12个IGS台站观测的地表形变值进行了比较，发现在所有站上两类形变的振幅之比约为1；文中存在的问题是：（1）参与比较的IGS站太少；（2）采用的GRACE数据时间段较短，而且GRACE早期观测资料由于星载设备调试和轨道调整等原因质量较差^[6].2007年，vanDam等^[7]利用GRACE前34个月数据解算欧洲的地表周年垂直形变，滤波方法采用半径为5×10⁵m的“高斯平滑”，但是GRACE数据时间段依然较短.此外，2006年King等^[8]、2009年Tregoning等^[9]也在澳大利亚等区域开展过利用GRACE卫星监测地表周年形变的研究.

本文采用前72个月的GRACE卫星二级产品，基于两种不同的滤波方法，探求中国及邻区的地表周年垂直形变，并与区内42个GPS台站的观测值进行了比较，对于了解GRACE重力卫星结果监测地表形变的数据处理方法，探求重力卫星观测技术为地表形变监测服务有一定的科学意义.

GRACE月平均重力场模型的数据处理和存档主要由3个研究中心共同协作完成并发布：CSR（美国德克萨斯大学空间研究中心，Center for Space Research）、JPL（美国宇航局喷气推进实验室，Jet Propulsion Laboratory）和GFZ（德国地学中心，GeoForschungsZentrum Potsdam）.本文的月平均重力场模型取自CSR，为归一化的60阶球谐系数，时间跨度为2002年4月至2008年3月.缺失2002年6月、7月以及2003年6月的重力场模型分别用各自前后2个月的重力场模型的平均值代替.由于在反演月平均重力场模型的过程中扣除了固体潮、海潮、极潮以及非潮汐海洋和大气质量负载的变化，从而该球谐系数包含了陆地水信号、少数构造活动信号以及未被模型完全模拟的残留信号.

作为对比，我们采用与GRACE数据对应时间段的总共42个GPS台站的坐标变化时间序列（台站位置见图 1）.这些台站包括中国地壳运动观测网络的27个基准站、武汉九峰WHJF站以及SOPAC（Scripps定轨与固定台阵中心，ScrippsOrbitand Permanent Array Center）提供的区内16个IGS（国际GNSS服务，International GNSSS ervice）站，其中曼谷CUSV站的数据不在要求的时间段内，新竹TCMS站在时间段内的资料记录仅有466天，均不考虑.这些GPS站在解算过程中扣除了固体潮、海潮、极潮，但是没有扣除大气负载的影响，扣除的海潮模型中也没有包括非潮汐海洋质量负载.

为了研究气候驱动的弹性地表垂直形变并与GPS结果进行比较，我们需要加回在反演GRACE月平均重力场过程中扣除的非潮汐海洋和大气质量负载的变化.为了统一参考框架，GRACE时变重力场还必须加上一阶重力位系数，见文献[5]、[9].由于GRACE恢复的重力场模型没有一阶分量，我们采用Swenson等^[10]2008年估算的一阶重力位系数.此外，GRACE卫星的观测方式造成重力位系数C20不稳定，文中采用激光测距数据分析的C20项^[11]来替换.

为了得到GRACE月平均重力场相对某一静态（背景）重力场的时间变化，我们采用CSR2004年发布的全球重力场模型GGM02S（GRACE Gravity Model02，S表示只使用了GRACE卫星的数据）的前60阶次作为背景重力场，它在潮汐和误差参数上与GRACE月平均重力场最接近.

根据Blewitt^{[12, 13]}、Wahr^[14]的结果，推导由GRACE时变重力资料解算气候驱动的地表垂直形变：

(1)

式中，ΔU为地表垂直形变，ϕ、λ分别为地面网格点的纬度和经度，R_E为地球的平均半径，N_max为月平均重力场模型的最大阶数，本文取值为60，k_n与h_n为负荷勒夫数，本文采用SNREI地球模型表面负荷勒夫数^[15]，P_n，m是完全规格化的缔合勒让德多项式，ΔC_n，m、ΔS_n，m为n阶m次月平均重力场相对背景重力场的球谐系数之差.

为降低高阶球谐系数中的误差影响，需要引入滤波算子.本文分别测试并比较了两类滤波算子的效果，一类是Chambers^[16]提出的“去相关”滤波方法，另一类是Chambers“去相关”方法结合半径3×10⁵m“高斯平滑”的组合滤波，其中“高斯平滑”函数采用Jekeli^[17]提出的算子.分析过程中，为了提取气候驱动的流体质量重新分布产生的形变信号，假定地球表面固体质量重新分布对地表垂直形变的影响是线性的，通过扣除线性趋势项排除固体质量负载的影响.

GPS台站的坐标变化时间序列，在扣除线性趋势项排除固体质量负载的影响以外，还需要经过以下处理：（1）原始序列中去掉大于2.6倍标准偏差的异常值，再用线性插值构成完整的分析时间序列，便于下文拟合.（2）考虑到GPS每天1值的分辨率远大于GRACE每月1值的分辨率，用小波分解分析时间序列后，保留各种长周期信号；本文采用Coiflets小波，见参考文献[4].（3）带通滤波提取1月~1年周期的GPS季节性信号.

为了得到地表垂直形变的周年平均振幅和初相位，无论是3.2节中GPS观测的还是3.1节中GRACE解算的形变时间序列，都需要进行非线性最小二乘拟合，（2）式为拟合式，拟合结果见表 1.

(2)

表 1(Table 1)

表 1 GRACE和各GPS站垂直形变时间序列的拟合结果 Table 1 Fitted results of amplitude and phase derived from GRACE and GPS

表 1 GRACE和各GPS站垂直形变时间序列的拟合结果 Table 1 Fitted results of amplitude and phase derived from GRACE and GPS

式中t_i为给定的时间点，表示为相对某一参考时间历元的间隔，以“天”为单位，本文的参考时间历元选择2002年4月1日，Y（t_i）为输入的原始序列，a表示常数项，b表示线性趋势项，ω=2π/T（本文取T等于365天），A_l和φ_l分别表示周年（l=1）和半周年（l=2）的振幅和初相位.

在表 1中的GPS台站上，GPS周年振幅（以下简称“GPS振幅”）分别结合GRACE（D）和GRACE（DG）周年振幅（以下简称“GRACE振幅”），用“（GRACE振幅-GPS振幅）/GPS振幅”表示GRACE振幅偏离GPS振幅的程度（以下简称“振幅偏离程度”）.

计算发现，振幅偏离程度与GRACE滤波方法（GRACE（D）、GRACE（DG））有关.在GUAN、HLAR、JIXN、LHAS、URUM、XIAA、XIAM、YANC、WHJF、MALD、SUWN、ULAB、IRKT、SELE、POL215个台站上，分别用两类滤波计算后的振幅偏离程度差距在5%（含）以内；在LUZH、SUIY、TAIN、WUHN、XNIN站上差距在5%~10%（含）以内，在这20个站上认为两类GRACE滤波方法的效果相当；在BJFS、BJSH、DXIN、DLHA、HRBN、WUSH、XIAG、YONG、TWTF、TNML、BAN2、NTUS、HYDE、PIMO、DAEJ15个站上，用“去相关”滤波计算后的振幅偏离程度优于用组合滤波计算后的振幅偏离程度；在CHUN、KMIN、QION、SHAO、TASH、ZHNZ、IISC7个站上，用组合滤波计算后的振幅偏离程度优于用“去相关”滤波计算后的振幅偏离程度.另一方面，GRACE用“去相关”滤波计算后，振幅偏离程度大于30%的台站有15个；小于20%（含）的台站有17个，这其中小于10%（含）的台站达到12个，占70.58%.GRACE用组合滤波计算后，振幅偏离程度大于30%的台站有16个；小于20%（含）的台站有16个，而这其中小于10%（含）的台站仅有8个，只占50%.因此，GRACE用“去相关”滤波计算后的振幅总体上优于用组合滤波计算后的振幅.另外，从表 1中的周年初相位可见：GRACE用“去相关”滤波计算后的相位与GPS相位偏差小于半个月的台站有19个，偏差在半个月与2个月之间的台站有19个；GRACE用组合滤波计算后的相位与GPS相位偏差小于半个月的台站仅有12个，偏差在半个月与2个月之间的台站却有27个；因此，GRACE用“去相关”滤波计算后的相位总体上优于用组合滤波计算后的相位.从振幅和相位上综合考虑，GRACE解算中国及邻区的地表周年垂直形变时，采用“去相关”滤波方法达到了预期的结果-与GPS观测的结果总体上一致，效果优于组合滤波方法.下文的GRACE数据处理中一致使用“去相关”滤波方法.

图 1中背景表示了中国及邻区的GRACE振幅，反映了气候驱动下地表流体质量引起的地表周年垂直形变分布.GRACE振幅大小因地而异，在大陆上最小的地方位于TASH站东南，约1×10^-3m；最大的地方位于恒河-澜沧江流域，可达10×10^-3m.其中的蓝色箭头和红色箭头分别表示各台站上GRACE和GPS的振幅和初相位，可见在大多数台站上振幅与初相位都符合较好.TASH、GUAN、KMIN、HYDE站的相位相差较大，ULAB、IRKT、SELE、POL2、URUM站的GRACE振幅比GPS振幅偏小，主要因素可能与GPS解算策略、GPS观测资料的连续性或局部大气、水文过程等地球物理因素有关，具体原因有待进一步研究.

图 2列出了区内8个台站的地表季节性垂直形变时间序列.图中蓝色线（夹蓝色方形符号）表示GRACE的地表垂直形变时间序列；浅灰色线表示GPS的地表垂直形变分析时间序列，其中URUM站的虚线椭圆表示其中的分析时间序列由于缺少GPS原始观测资料而进行了线性插值；红色线表示经过小波提取和带通滤波后GPS的地表垂直形变季节性信号.

图 2

Fig. 2

图 2 区内8个GPS台站上地表垂直形变的时间序列 蓝色线（夹蓝色方形符号）表示GRACE的时间序列；浅灰色线表示GPS的分析时间序列，URUM站的虚线椭圆表示其中的分析时间序列由于缺少原始资料而进行了线性插值；红色线表示经过小波提取和带通滤波后的GPS季节性信号. Fig. 2 Time series of vertical deformation of the solid Earth derived from GRACECblue) and GPS (red) at eight GPS sites in China and around. GPS (light gray) means GPS time series in processing and linear interpolation have been done in dashed ellipse.

与地面GPS观测的地表形变具有全波段信号相比，GRACE的监测结果只使用了前60阶中长波信号，比60更高阶的信号对地表垂直形变的“贡献”有多大呢？以2002年4月至2008年3月CPC（气候预测中心，Climate Prediction Center）的月平均陆地水储量载荷^[18]为例，模拟地表流体质量载荷前60阶信号对地表垂直形变的“贡献”.将月平均陆地水储量载荷用球谐函数展开至360阶，记为资料1；将其用球谐函数展开至360阶后取前60阶，记为资料2.分别模拟资料1和资料2质量载荷引起BJFS站（图 3a）和WUHN站（图 3b）的地表垂直形变.

图 3

Fig. 3

图 3 BJFS站（a)和WUHN站（b) CPC模型确定的陆地水质量载荷引起的地表垂直形变 Fig. 3 Vertical deformation of the solid Earth derived from degree 360 and 60 CPC model in (a) BJFS and (b) WUHN

陆地水质量载荷的前60阶信号导致的地表垂直形变占360阶信号引起形变的主要部分.图 3b可见，在WUHN站，信号中第61~360阶对地表垂直形变的“贡献”很微弱；在BJFS站“贡献”只占较小一部分.

GPS监测获取的地表形变是单点的地面垂直形变；大气、海洋、陆地水等地球物理因素综合模拟的地表形变又受到模型精度不确定性的制约；卫星重力具有能够反演地表质量时空再分布的独特优势，利用其动态和空间连续结果监测地表形变是一种新的技术手段，它的监测结果包含了各种地表负荷响应，可以避免用陆地水等地球物理模型综合模拟地表形变时的精度不确定性问题.本文研究了利用GRACE卫星时变重力资料解算气候驱动的地表垂直形变的数据处理方法，并计算了中国及邻区的地表周年垂直形变.通过与区内42个GPS台站上观测到的季节性信号进行比较，发现GPS台站观测到的周年垂直形变基本上能够被GRACE卫星监测到.

由于GRACE卫星时变重力资料的高阶系数主要是噪声，直接影响解算的地表形变.为了降低高阶系数误差对解算地表垂直形变的影响，我们分别比较了两类滤波方法.综合考虑，用GRACE解算区内的地表周年垂直形变时，采用“去相关”滤波方法处理后的结果优于采用组合滤波处理后的结果.用“去相关”滤波计算后，GRACE振幅、相位和GPS振幅、相位总体上一致；少数站上GRACE和GPS的振幅或相位相差较大，可能与GPS解算策略、GPS观测资料的连续性或局部大气、水文过程等地球物理因素有关.区内GRACE振幅的最大值出现在恒河-澜沧江流域，可达10×10^-3m，在大陆上最小值位于TASH站东南，约1×10^-3m.

通过地表流体质量载荷模拟，发现陆地水质量载荷球谐函数展开后的第61~360阶信号对地表垂直形变的“贡献”比较小，因此本文只取GRACE观测的前60阶中长波信号，就能够监测地表季节性垂直形变的主要部分.本文对于探讨重力卫星结果监测地表形变的数据处理方法有一定的科学意义.

中国地壳运动观测网和武汉九峰站提供GPS数据；中国科学院测量与地球物理研究所李军助理研究员在资料收集中提供的有益帮助；Marcel.liFarran提供WIN4GMT Version2.1.4，在此一并表示感谢.