反演地表质量变化的附有方差约束的径向点质量方法

尹恒， 游为， 范东明， 方伟浩， 万祥禹， 宋梦芝

西南交通大学地球科学与环境工程学院， 成都 611756

0 引言

2002年以来，由美德联合实施的GRACE重力卫星以前所未有的精度和时空分辨率观测时变地球重力场信号，提供了研究地球表层质量变化的直接观测手段(Wahr et al., 2004; Tapley et al., 2004b).近20年来，基于GRACE时变重力场进行的全球或局部陆地水储量变化研究(Rodell et al., 2007；冯伟, 2013；Zhong et al., 2018)、极地冰盖和高山冰川质量变化研究(Luthcke et al., 2006b, 2013；Baur, 2013；朱传东等, 2015；高春春等, 2016)、海平面变化研究(Chambers, 2006；Uebbing et al., 2019)等取得了丰硕成果，促进了大地测量学、地球物理学、海洋学和冰川学等领域的发展，因此利用重力卫星观测数据研究全球质量变化具有重大意义.

目前利用时变地球重力场表达地表物质迁移的方法主要分为球谐位系数法(Wahr et al., 1998)、Mascon法(Rowlands et al., 2005；Luthcke et al., 2006a)、点质量法(Baur and Sneeuw, 2011)和径向基函数(杨帆, 2017)等方法.由于GRACE卫星本身的仪器载荷误差、近极近圆飞行轨道和高频背景摄动力模型误差等问题，使得高阶次球谐位系数被噪声严重污染，导致直接使用球谐位系数反演的地表质量变化呈现严重的南北条带误差(Tapley et al., 2004a).许多学者研究了众多滤波平滑技术，如常用的多项式拟合去相关算法和高斯滤波方法(Wahr et al., 1998；Swenson and Wahr, 2006).虽然滤波方法的实现使得利用球谐位系数法计算地表质量变化被广泛应用，但滤波在过滤噪声的同时不可避免地会过滤掉有用的信号甚至改变信号的空间分布，相比于对球谐位系数施加平滑滤波，对地球物理信号施加时间和空间约束能得到更高精度和空间分辨率的地表质量变化反演结果(Luthcke et al., 2006a)，因此有学者提出了另一种计算地表质量迁移的Mascon方法.

Mascon方法最早被用于月球重力场模型的确定(Muller and Sjogren, 1968)，其基本思路是将研究区域划分为若干质量块并假设每个质量块内部质量变化相等，直接以每个质量块的质量变化作为未知参数，建立卫星轨道和星间距离变率等原始观测值与未知参数的函数关系.Rowlands等(2005)、Luthcke等(2006a)和Rodell等(2007)将地表质量块与GRACE卫星K波段星间距离变率建立联系，利用时空约束方程求解地表质量变化，通过与水文模型和卫星测高对比，验证了该方法的有效性.Watkins等(2015)利用星间距离变率作为观测值，引入GLDAS (Global Land Data Assimilation System)、卫星测高和ECCO2(Estimating the Circulation and Climate of the Ocean Phase 2)海洋模型等外部物理信号建立时空约束方程反演等效水高变化，Save等(2016)利用球谐系数法反演结果作为先验信息构建约束矩阵反演全球质量变化.

有学者发现也可直接建立地表质量块与重力卫星所受摄动加速度(或重力位)的函数关系来求解地表质量变化，该类方法被称为点质量方法.目前点质量方法主要有扰动位点质量方法(Han et al., 2005)、径向点质量方法(Baur and Sneeuw, 2011)以及三维加速度点质量方法(苏勇等, 2017)，其观测值分别为扰动位变化、径向加速度变化以及三维加速度变化.Han等(2005)使用扰动位点质量方法反演了亚马逊流域和奥里诺科河流域的质量变化，表明其空间分辨率优于球谐位系数法结果.Baur和Sneeuw(2011)、Baur(2013)利用径向点质量方法计算了格陵兰岛冰盖质量变化，验证了该方法的有效性.郭飞霄等(2018)、苏勇等(2019)和魏伟等(2020)分别利用附有空间约束的径向点质量方法和三维点质量方法反演了亚马逊流域和华北平原陆地水储量变化，对比分析了不同空间约束的效果，表明采用的各类空间约束均能提高反演结果精度，其中线性形式的空间约束效果最好.Chen等(2016,2020)利用径向点质量方法通过截断奇异值分解和Tikhonov正则化方法计算青藏高原冰川冻土消融情况，提高了计算结果的信噪比.Forsberg等(2017)利用径向点质量方法计算了格陵兰岛和南极大陆冰盖融化情况对海平面上升的贡献.Ran等(2018)以地表质量变化先验误差信息为约束条件，利用基于特征值分解的径向点质量方法计算了格陵兰岛冰盖消融情况.大量研究结果已验证了点质量方法的可行性和有效性，但由于约束矩阵的缺陷，其反演陆地水储量变化结果仍受条带误差和海陆信号泄漏的影响.

本文仿照CSR Mascon正则化矩阵的构造思路，几乎不加入任何外部数据的质量变化信息，直接以GRACE时变重力场模型的计算结果作为先验质量变化信息，对CSR Mascon正则化矩阵的构造思路进行若干改进，推导出了附有方差约束的径向点质量方法.具体改进的地方包括：(1)以单位矩阵为正则化矩阵的径向点质量方法计算的质量变化的RMS作为先验质量变化信息；(2)考虑各种质量变化信号的正则化参数和正则化矩阵不唯一，构造了迭代正则化的处理策略.

1 数据与方法

1.1 数据处理与验证数据

本文选用CSR发布的2003-01—2014-11共133个月的Level-2 Release 06 GSM时变重力场模型数据反演地表质量变化, 并将球谐位系数截断至60阶.由于GRACE时变重力场解算的一阶项为0(C10,C11,S10)，采用Sun等(2016)解算的一阶项替换原有一阶项(Swenson et al., 2008； Sun et al., 2016).由于GRACE解算的C20项精度较低，采用SLR (Satellite Laser Ranging,卫星激光测距)解算的C20替代GRACE球谐位系数中的C20项(Cheng et al., 2013； Loomis et al., 2020).

采用最新的GLDAS_NOAH_M.2.1水文模型(简称GLDAS2.1)用于评估本文反演的陆地水储量变化精度.该模型是一组地表水和能量场的数据集，由GSFC(Goddard Space Flight Center,戈达德航天中心)和NCEP(National Centers for Environmental Prediction,美国国家环境预报中心)联合研制(Rodell et al., 2004)，本文采用其提供的空间分辨率为0.25°×0.25°，时间分辨率为1个月的模式数据.

为更好地评价反演结果，本文采用CSR和JPL RL06 Mascon(以下简称CSRM和JPLM)模型进行全球、陆地、冰冻圈信号综合对比.由于上述Mascon模型均扣除了2004年1月至2009年12月的均值，为与MASCON模型保持一致，本文所用球谐位系数也扣除相同时段的球谐位系数均值(Watkins et al., 2015；Save et al., 2016).

1.2 空间格网

径向点质量方法通过将地表点质量变化与空间点重力扰动建立函数关系求解区域或全球质量变化，因此需要事先将地球表面划分为不同质量块，并确保所有质量块均匀地覆盖在地球表面.为确保划分地表质量块面积的均匀性，本文采用基于二十面体划分的等面积格网，格网划分方法可参考Majewski等(2002).为了有效减少海陆泄漏误差，本文按上述方法将地球表面划分为40962个等面积格网，每个格网面积约为12450 km2，格网点之间距离约为120 km，即格网点空间分辨率约为1°.图1为本文用于径向点质量方法的等面积海陆格网空间分布图，其中紫色表示陆地格网，蓝色表示海洋格网(图中白色部分为地图投影问题，真实格网已完全覆盖地球表面).

1.3 径向点质量模型方法基本原理

地球表层质量变化引起的重力扰动径向分量可由球谐位系数变化量 ΔClm,ΔSlm表示为(Baur and Sneeuw, 2011)：

×[ΔClmcos(mλ)+ΔSl msin(mλ)],(1)

图2 径向点质量方法几何示意图Fig.2 Radial point-mass modeling geometry

如图2所示，在地心坐标系中有质量变化为 δmj的地面点Pj(λj,φj,a)，根据牛顿第二运动定律，其引起的外部空间点Si(λi,φi,r)的径向重力扰动δgij可表示为：

(2)

cosψij=sinφisinφj+cosφicosφjcos(λi-λj),(3)

式中，lij是空间点Si到地面点Pj的距离，ψij为空间点Si到地面点Pj的球面角距.根据图2中的几何关系，特定区域内所有地表点质量变化引起的某空间点径向重力扰动可表示为：

(4)

式(4)左端为时变重力场球谐位系数计算的径向重力扰动伪观测值，式(4)右端δmj为待求的单个地表质量变化未知数.当利用等效水高变化来表示地表质量变化时，可设δmj=ρw×sj×Δhj，其中ρw=1025 kg·m-3为水的密度，sj为单个格网面积，Δhj为等效水高变化，将单位换算为cm后，式(4)表示为：

i=1,…,s.

(5)

联合式(1)和式(5)即可得利用径向点质量方法计算地球表面等效水高变化的函数模型，其中i=1,…,s为空间点的个数，j=1,…,p为研究区域地面格网点的个数，本文令空间点个数s和地面格网点个数p均为40962.进一步，将式(5)线性化可得线性方程：

y=Ax+e,(6)

式中y(s×1)表示观测值向量，A(s×j)表示设计矩阵，x(j×1)为待求未知数向量，e(s×1)为观测值拟合残差.

1.4 参数估计方法

由于利用卫星高度处的重力扰动观测值计算地表质量变化本质上是一个重力向下延拓问题，会使设计矩阵A呈现病态性，因此直接用最小二乘法则求解式(6)不能获得稳定解(Kusche and Klees, 2002).为了得到稳定的估值，本文利用Tikhonov正则化原理，引入正则化参数和正则化矩阵作为约束条件，构造以下估计准则(王振杰, 2003)：

min{‖y-Ax‖2+α‖Mx‖2},(7)

将式(7)对x求偏导数并联合式(6)可得待估参数的解为：

(8)

式中α是正则化参数，用于平衡观测值残差和约束条件；M为正则化矩阵，为方程提供约束条件，‖·‖是欧式2-范数.

与最小二乘解相比, 式(8)右端求逆部分加入了αMTM一项，由于它的作用抑制了法方程的病态性, 使得(ATA+αMTM)的求逆变得稳定, 因此式(8)可以得到稳定的估值.式(8)的结果与正则化参数α有关,过小的正则化参数会导致解仍然不稳定，过大的参数会使反演结果过度平滑.常见选取正则化参数的方法有L-曲线法(Hansen, 1992)、U-曲线法(Krawczyk-Stańdo and Rudnicki,2007)、GCV(Generalized Cross-Validation，广义交叉验证)方法(Xu, 2009)和RMSE最小法(Xu et al., 2006)等.本文使用Hansen(1992)针对线性病态问题提出的L曲线法估计正则化参数，以‖y-Ax‖为横坐标，‖x‖为纵坐标获得L曲线图，L曲线上曲率最大的点对应参数即为L曲线法确定的正则化参数，但需要注意的是此方法确定的正则化参数仅是近似最优的.

作为约束条件的正则化矩阵M有多种构造方法，最初的方法是直接将M设为单位矩阵(Baur and Sneeuw, 2011; 郭飞霄等, 2017)，后来有学者利用空间平滑条件设计正则化矩阵，以此对地表质量变化点进行空间平滑约束(郭飞霄等, 2018; 魏伟等, 2020).直接利用单位矩阵或空间平滑条件的约束矩阵均能有效去除部分条带误差，但其结果仍受部分条带误差影响(苏勇等, 2019).

1.5 无空间约束径向点质量方法反演地表质量变化

本节利用CSR RL06 GSM数据通过正则化矩阵为单位矩阵的无空间约束径向点质量方法(以下简称PM)计算2003-01—2014-11共133个月的全球等效水高，图3为PM反演各月份等效水高的L曲线图，不同月份对应的L曲线有所差异，取所有L曲线均值作为平均正则化曲线.

图3 2003-01—2014-11不同月份L曲线(绿)和平均L曲线(红)Fig.3 L curve (green) and average L curve (red) in different months from January 2003 to November 2014

图3中的红色曲线为所有月份平均L曲线，当正则化参数为1×10-19时，L曲线曲率最大，因此将其作为所有月份的正则化参数.为了比较不同正则化参数的效果，同时选取了大于和小于最佳正则化参数反演结果作为对比，绘制了如图4的2010-01等效水高分布图，为了更直观地对比反演效果，同时绘制了滤波半径为300 km的高斯滤波反演结果分布图.

对比图4b和4d可以看出，当正则化参数为1×10-19(最优正则化参数)时，PM反演结果仍受条带误差和泄漏误差影响，其条带误差与滤波半径为300 km的高斯滤波反演结果相当.当增大正则化参数至3×10-19时(图4c)，PM反演结果信号强度随之减小，说明过大的正则化参数会抑制反演结果的信号强度，但与此同时条带误差也在减小.当减小正则化参数为5×10-20时(图4a)，信号强度和条带误差都显著增大.从以上分析可以看出PM仍受条带误差和泄漏误差影响，通过增大正则化参数可减少条带误差，但也会严重削弱信号.

2 附有方差约束的径向点质量方法

针对无空间约束的PM仍受条带误差和海陆信号泄漏影响的缺陷，本文利用无空间约束的PM反演结果设计正则化矩阵，构建附有方差约束的径向点质量方法，旨在通过先验信息的约束，消除条带误差和海陆泄漏误差.

2.1 正则化矩阵

VCPM需要借助质量块先验方差信息设计正则化矩阵，Watkins等(2015)利用陆地水文模型、卫星测高和海底压力模型等多源数据的方差信息构造约束矩阵，Save等(2016)利用GRACE球谐位系数法反演结果估计质量块的方差信息构建正则化矩阵.式(9)为本文所用的正则化矩阵结构，由于Save等(2016)指出质量块之间的协方差信息会给反演结果带来横向误差，因此仅取质量块方差信息设计正则化矩阵.

(9)

图4 不同正则化参数的点质量方法(a,b,c)和高斯滤波半径为300 km的球谐位系数法(d)反演的2010-01等效水高分布图Fig.4 Equivalent water height estimated by (a,b,c) point-mass method with different regularization parameters or (d) spherical harmonic potential coefficient method with Gaussian filter radius of 300 km in January 2010

式中σi,i=1,2,…,n为第i个质量块的方差信息，通过计算第i个质量块多个月的RMS值得到.Save等(2016)设计了两种正则化矩阵，一种是直接以2004—2009年等效水高RMS设计的0409F，另一种是在0409F基础上再结合其他单个月的等效水高信息设计的附有时间相关性的0409V.结果表明：由于0409V包含了2004—2009年以外的陆地水储量变化信息，因此以0409V反演的2004—2009年时间段以外的等效水高解效果优于0409F反演结果，而在2004—2009年时间段以内，两个正则化矩阵反演结果相同.本文直接以2003—2014年全时段PM反演结果构建正则化矩阵，包含了时间段内的所有地表水储量变化信号特征，因此不再构建时间相关的正则化矩阵.

本文的正则化矩阵分两步设计，第一步：设计初始正则化矩阵，获得中间解，这一步主要是为了消除海陆泄漏误差，同时获得海洋区域质量块真实方差信息；第二步：利用第一步获得的真实海洋方差信息和初始陆地质量变化信息构建迭代正则化矩阵，通过迭代的方法逐步提取质量变化信号.对于第一步的初始正则化矩阵，本文与Save等(2016)的方法完全相同，第二步中Save等(2016)采用扣除了陆地振幅信号和高纬冰川趋势信号的剩余信号构建正则化矩阵，利用该矩阵一次性反演全球等效水高信号，本文为充分提取观测值中的物理信号设计了三个不同的正则化矩阵迭代地反演全球等效水高.

2.2 初始正则化矩阵设计

Save等(2016)利用Save等(2012)正则化的时变重力场模型经200 km高斯滤波反演的等效水高结果构建初始正则化矩阵，本文使用信噪比更高的PM计算结果构建初始正则化矩阵.图5a为利用PM反演的133个月等效水高计算的均方根值(RMS)分布图，从中可以直观地看到许多代表质量变化的物理信号：如陆地冰川冰盖融化信号，各大流域季节性水储量变化信号，印度西北部、华北平原、美国加利福利亚地区地下水消减信号等；除了陆地信号之外，海洋区域也能看出部分质量变化较强的信号，其中苏门答腊、东京和智利等海底地震信号最为明显，里海、黑海、红海、澳大利亚卡奔塔利亚湾等封闭半封闭海域质量变化也较为明显.同时能看到在长期趋势较大的冰川冰盖融化区域，海陆泄漏误差极为严重，如南极、格陵兰岛、阿拉斯加、巴塔哥尼亚冰原等(谷延超, 2018).

为了消除海陆泄漏误差，同时保证封闭、半封闭海域信号和地震海域信号的恢复，对图5a的RMS信号做如下特殊处理：(1)保持陆地区域RMS信号不变；(2)保留质量变化明显的海底地震区域RMS(苏门答腊地震、东京地震和智利地震)以及封闭、半封闭海域RMS(里海、红海、黑海、哈德孙湾、福克斯湾、澳大利亚卡奔塔利亚湾、泰国湾、北冰洋部分信号明显海域)；(3)将其余海洋区域RMS均设置为3 cm(因为大多数海洋区域RMS小于3 cm)(Save et al., 2016; 谷延超, 2018).

图5b为经以上特殊处理后的RMS信号分布图，将此RMS格网点元素作为式(9)M-1的对角线元素，即获得初始正则化矩阵.将以上初始化矩阵代入式(8)，用L曲线法寻找最佳正则化参数，即可获得利用初始正则化矩阵求得的地表质量变化，称为中间解.图6a和6b分别为2010-01的无空间约束点质量解和利用初始正则化矩阵求得的中间解，通过对比可以看出：中间解在初始正则化矩阵的约束作用下基本消除海陆泄漏误差，条带误差也大幅减弱；此外，中间解等效水高信号强度和空间分辨率也大大提高.信号增强现象在质量变化较大的区域最为明显，如：亚马逊流域、赞比西河流域以及南北极冰川冰盖融化区域；但也能看出在质量变化较小的区域，信号几乎没有得到增强，如：青藏高原、华北平原和美国加利福利亚区域等.出现上述信号恢复情况不一致的原因是：全球质量变化信号复杂，仅使用单一的正则化参数不能同时兼顾质量变化强的信号和弱的信号的特征.

图5 利用PM(Point-Mass，点质量)方法反演的2003-01—2014-11共133个月等效水高计算的(a)RMS分布图和(b)初始正则化矩阵RMS分布图Fig.5 (a) RMS distribution and (b) initial regularization matrix RMS distribution of equivalent water height calculated from January 2003 to November 2014 by PM

2.3 迭代正则化矩阵设计

图5b的初始正则化矩阵RMS分布图反映的是全球各陆地区域质量变化大小的分布，从图中可以直观地看出，全球质量变化大小分布极度不均匀，要消除质量变化较强区域的条带误差，需要较大的正则化参数进行约束，但对于信号较小的区域，过大的正则化参数会削弱其真实信号(Save et al., 2016)，因此单一的正则化参数很难平衡全球质量变化差异很大的情况.为兼顾全球质量变化信号较强区域的条带误差消除，以及质量变化较小区域的信号提取，可构建迭代的附有方差约束的径向点质量方法，逐步提取时变重力场模型信号(Luthcke et al., 2013； Arendt et al., 2013).本文第一次迭代的正则化矩阵为利用133个月的中间解直接计算的RMS，其包含所有的质量变化信号；第二次用于迭代的正则化矩阵为扣除全球陆地区域等效水高季节性信号和高纬度冰盖区域及海底地震区域长期趋势信号后剩余信号计算的RMS，其主要包括低纬度地区趋势信号较弱的质量变化信号；第三次迭代的是扣除所有陆地区域季节性信号和所有陆地及地震海域长期趋势信号后剩余信号计算的RMS，其主要包括海洋信号和黑海、红海等内陆海域信号计算的RMS.图7a、7b和7c分别为三次迭代的正则化矩阵RMS分布图.

图6 分别利用单位矩阵和初始正则化矩阵的点质量方法解算的2010-01(a)无约束解和(b)中间解等效水高分布图Fig.6 Equivalent water height distribution of January 2010 (a) unconstrained solution and (b) intermediate solution calculated by point-mass method using unit matrix and initial regularization matrix respectively

图7 用于VCPM方法的不同迭代次数的正则化矩阵RMS分布图Fig.7 Regularized matrix RMS distribution for different iterations of VCPM

其中第一次迭代所用的观测值是式(1)直接计算出的重力扰动；第二次迭代的观测值是利用式(6)计算出的第一次迭代剩余的重力扰动残差，其包含的是第一次迭代未提取的重力扰动信号；第三次迭代的观测值是前两次迭代剩余的重力扰动残差，包含的是前两次迭代剩余的微小重力扰动信号.

迭代正则化矩阵按照上述方法设计的原因是：第一次迭代的正则化矩阵包含了所有质量变化特征，因此需要较大的正则化参数来抑制质量变化信号较大区域(如格陵兰岛)的条带误差，因此第一次迭代主要提取的是质量变化较大的信号，而质量变化较小的信号会留在重力扰动观测值内，留待下一次迭代提取；第二次迭代的正则化矩阵由于扣除了高纬冰川融化趋势和陆地振幅等信号，因此仅需要较小的正则化参数就能抑制条带误差，本次迭代就可以很容易的提取出重力扰动残差中剩余的质量变化较小的信号(如华北平原地下水信号)；第三次迭代的正则化矩阵扣除了所有的陆地振幅和趋势信号，因此可以使用更小的正则化参数提取前两次未完全提取的海洋信号和内陆海洋信号，因此三次迭代提取的信号累加即可得到观测值中所有的质量变化信号.

Save等(2016)利用扣除高纬冰川趋势信号和陆地振幅信号的剩余信号RMS构造正则化矩阵，以此正则化矩阵一次性提取观测值中蕴含的质量变化信号，因为这样可以较好解决单一正则化参数不能平衡正则化矩阵中质量变化大小极不均匀的问题，但是这样做也会导致先验方差信息不真实的问题；而本文利用迭代的方法解决单一正则化参数的适应性问题，避免了先验方差信息不真实的情况.

2.4 附有方差约束的径向点质量方法反演结果

根据上述正则化矩阵的设计思路，附有方差约束的径向点质量方法的计算流程图如图8所示，利用VCPM方法计算了2003-01—2014-11的全球等效水高变化，以2010-01为例，分析此方法单个月的反演效果.图9为不同迭代次数的等效水高变化分布图，其中图9a，9b，9c分别是第一次到第三次迭代的反演结果，图9d为三次迭代结果之和，也即VCPM计算的2010-01等效水高变化结果.为了反映不同迭代次数提取的重力扰动物理信号，图10给出了2010-01重力扰动初始观测值和不同迭代次数后的重力扰动残差.可以看到:第一次迭代提取了重力扰动观测值中包含的大多数质量变化信号，其中质量变化较大的信号(如南极、亚马逊、奥卡万戈三角洲等信号)已经基本提取完毕，但是从图10b第一次迭代的重力扰动残差中可以看出，质量变化较小的信号(如华北平原、青藏高原、堪察加半岛和黑海)仍有部分未被提取；因此第二次迭代提取的主要是质量变化较小的信号，如：黑海、华北平原、青藏高原、非洲西部热带雨林等；第三次迭代主要提取的是前两次迭代剩余的微小海洋信号和红海、黑海等内陆海信号，此时图10d第三次迭代后的重力扰动残差中已基本不可见质量变化的物理信号，剩余的信号主要是条带误差，说明VCPM经过三次迭代能够充分提取重力扰动观测值中蕴含的质量变化信号.

图8 附有方差约束的径向点质量方法的计算流程图Fig.8 Flow chart of variance-constrained radial point-mass method

为验证三次迭代相对于一次迭代的优势以及三次迭代提取的信号是否为真实信号，如图11所示，分别利用一次迭代和三次迭代VCPM的2003-01—2014-11的全球质量变化计算了相对于CSRM和JPLM两种Mascon数据的RMSE在20个流域的按面积加权平均值，另外也计算了CSRM和JPLM之间的RMSE在各流域的加权平均值，并将其作为参考值.从图中可以看出，在除塔里木河流域之外的绝大多数流域，三次迭代的VCPM与两种Mascon数据的RMSE平均值均要小于一次迭代的VCPM与两种Mascon数据的RMSE平均值，这说明三次迭代的VCPM相对于一次迭代的VCPM要更加接近CSRM和JPLM这两种Mascon数据，因此也能说明后两次迭代提取了更多代表质量变化的真实信号，三次迭代能够有效解决单一的正则化参数无法适应全球质量变化不均匀的问题.其次还可以看出三次迭代VCPM与CSRM的RMSE在各流域的加权平均值均要小于CSRM与JPLM的相应RMSE加权平均值，这说明VCPM与CSRM在这些流域之间的差异要小于CSRM与JPLM之间的差异，也能说明三次迭代的VCPM能够达到与这两种Mascon同等的精度.

图9 利用VCPM方法反演的2010-01不同迭代次数和三次迭代之和的等效水高分布图Fig.9 Equivalent water height distribution of different iterations and the sum of three iterations estimated by VCPM in January 2010

图10 2010-01重力扰动初始值和不同迭代次数的重力扰动残差Fig.10 Initial observation of gravity disturbation and residual of gravity disturbation with different iterations in January 2010

图11 不同方法反演的2003-01—2014-11等效水高之间的RMSE在20个流域的按面积加权平均值柱状图Fig.11 Area-weighted mean values of RMSE between equivalent water heights estimated by different methods in 20 river basins from January 2003 to November 2014

图12 VCPM与CSRM分别反演的2008年1、3、5月等效水高分布图Fig.12 Equivalent water height distribution estimated by VCPM and CSRM methods in January, March, May 2008 respectively

为直观展示VCPM反演质量变化的去条带和泄漏误差效果，将VCPM反演的等效水高与CSRM进行对比.图12为 VCPM与CSRM分别反演的2008年1、3、5月等效水高分布图，可以看出，两种方法均已扣除条带误差及海陆泄漏误差的影响，提高了信号精度和空间分辨率，并且反演的等效水高分布基本一致，信号强度在大多数区域基本相同.以上对比分析表明本文设计的VCPM方法相对于传统滤波方法和PM方法更有优势，在抑制重力场模型噪声和获取高精度高分辨率等效水高变化效果上与CSRM模型方法相当，后文将继续在多个尺度上比较各反演结果优缺点.

3 不同方法反演结果比较与分析

为了具体评价VCPM反演地表质量变化的效果，引入CSRM、JPLM和GLDAS2.1水文模型进行对比分析.同时为了比较VCPM方法相对于传统滤波方法的优势，引入300 km高斯和去相关组合滤波算法(简称DS300)(Chen et al., 2007).由于下载的CSRM和JPLM数据中包含GAD产品，因此需要将其扣除；同时CSRM和JPLM均已经扣除冰川均衡调整GIA模型ICE-6G_D(VM5a)，因此VCPM和DS300也扣除该模型(Peltier et al., 2018).

3.1 全球尺度信号对比

为比较不同方法反演等效水高在多个月上的一致性，分别计算了2003-01—2014-11的VCPM与CSRM、JPLM、DS300数据和CSRM与JPLM、PM数据的均方根误差(RMSE)，如图13所示.同时表1给出了相应在全球、陆地和海洋区域的RMSE按面积加权平均值.从图13a中可以看出，VCPM与CSRM一致性较好，海洋区域相差较小，除苏门答腊、东京地震海域外的其余海域RMSE均很小；陆地区域中欧亚大陆、北美洲、澳大利亚和非洲北部等大部分地区RMSE很小，一致性较高；RMSE较大的区域主要分布在南极、格陵兰岛、阿拉斯加和巴塔哥尼亚等趋势信号较大的区域，除此之外，恒河流域、印度河流域和亚马逊流域北部RMSE也较大.CSRM与JPLM的RMSE分布图和VCPM与CSRM的RMSE分布图基本一致，RMSE较大区域与上述列举的位置基本相同.另外从图13c和13e中可以看出，DS300和PM在格陵兰岛等冰川冰盖融化区域存在巨大泄漏误差.

从表1的RMSE加权平均值来看，VCPM与CSRM的RMSE加权平均值在全球和陆地区域均最小(2.12 cm和4.16 cm)，优于CSRM与JPLM在全球和陆地区域的RMSE加权平均值(2.22 cm和4.62 cm).这可能是因为VCPM和CSRM均只采用GRACE观测数据导出的正则化矩阵反演等效水高，没有利用任何外部数据，而JPLM采用了GLDAS、卫星测高等外部物理模型设计正则化矩阵，另外VCPM与CSRM均为1°分辨率而JPLM为3°分辨率，因此在细节信号上VCPM与CSRM更为一致.在海洋区域RMSE加权平均值最小的是CSRM与JPLM，说明这两种数据在海洋区域反演结果最接近；在海洋区域RMSE最大的是CSRM与PM，这是因为PM在最优正则化参数下条带误差仍然较大且存在泄漏误差.VCPM与DS300在陆地区域的RMSE加权平均值最大，这是因为DS300由于平滑滤波导致了大量的信号泄漏.从以上对比分析可以看出，VCPM在利用先验方差信息设计的正则化矩阵约束下能有效去除条带误差和抑制泄漏误差.

表1 不同数据之间均方根误差按面积加权的均值(cm)Table 1 Area-weighted mean value of root mean square error between different data (cm)

本文也利用不同数据计算了全球长时间序列等效水高的周年振幅和长期趋势，如图14和图15所示.从周年振幅上来看，VCPM与CSRM和JPLM在空间分布和信号强度上表现出极好的一致性，尤其在亚马逊流域、赞比西河流域、印度北部区域、湄公河流域、格陵兰岛、阿拉斯加等区域信号相似.DS300由于滤波平滑作用，其周年振幅被严重削弱，PM较DS300稍好，但其高纬度区域振幅信号也被严重削弱.从长期趋势来看，DS300信号泄漏最为严重，很难分辨出信号较小的趋势信号.PM除冰川冰盖融化区域外，其余区域信号泄漏均较小，能清晰分辨出各种长期趋势变化的细节信息，这反映出PM方法优于DS300方法.VCPM、CSRM和JPLM三种数据在长期趋势上基本一致，其中VCPM和CSRM一致性略优于JPLM，这是因为地球质量变化信号中长期趋势信号更加复杂，而JPLM空间分辨率为3°，VCPM和CSRM空间分辨率均为1°，能够反映更多细节信号.从整体来看三者均能够有效提取出地球上主要的质量变化信号，如南极、格陵兰岛、阿拉斯加、巴塔哥尼亚及亚洲高山等地区的冰雪消融信号；青藏高原降雨增加信号、印度西北部地下水枯竭、中国华北平原地下水消耗、美国加利福利亚地区地下水消耗和德克萨斯州持续干旱信号、里海海平面下降、非洲奥卡万戈三角洲降雨增加信号等水储量变化信号；2004年苏门答腊地震、2010年智利地震和2011年东京地震信号(Chen et al., 2014; Scanlon et al., 2016; Peltier et al., 2018; Rodell et al., 2018;谷延超, 2018).从以上分析可以看出，在反演等效水高变化的周年振幅和长期趋势方面，VCPM方法明显优于DS300和PM方法，能够达到官方发布的Mascon数据的效果.

图13 不同方法反演的2003-01—2014-11全球等效水高之间的均方根误差Fig.13 Root mean square errors between global equivalent water heights estimated from January 2003 to November 2014 by different methods

图14 不同方法反演的2003-01—2014-11全球等效水高周年振幅Fig.14 Annual amplitude of global equivalent water height estimated by different methods from January 2003 to November 2014

图15 不同方法反演的2003-01—2014-11全球等效水高长期趋势Fig.15 Trends of global equivalent water height estimated by different methods from January 2003 to November 2014

3.2 陆地流域信号对比

为评价VCPM模型反演区域陆地水储量变化的可靠性，选取了亚马逊流域(Amazon)、奥里诺科河流域(Orinoco)、赞比西河流域(Zambeze)、育空河流域(Yukon)、印度河流域(Indus)、华北平原(NCP)等陆地流域和里海(Caspian)、红海(Redsea)等内陆海域作为研究对象进行对比分析.图16为不同数据在以上不同区域反演的2003-01—2014-11等效水高时间序列图，表2给出了将所有数据扣除季节项和趋势项后得到的CSRM与其余数据的时间序列相关系数，表3和表4分别给出了不同数据在各流域反演结果的长期趋势和周年振幅.

表2 CSRM与其余不同数据均扣除季节项和趋势项后的时间序列相关系数Table 2 Correlation coefficients of time series between CSRM and other different data after deducting the seasonal items and trend items

表3 不同方法反演的2003-01—2014-11等效水高在不同区域的长期趋势(cm·a-1)Table 3 Trends of equivalent water height in different regions estimated by different methods from January 2003 to November 2014 (cm·a-1)

结合上述图表可以看到：VCPM、CSRM和JPLM三种数据在各流域一致性都比较好，扣除季节项和趋势项后的时间序列相关系数在除红海之外的各区域仍能达到0.8以上，红海区域相关系数较差的原因可能是该区域较为狭小，各方法反演结果的空间分布差异较大；DS300和PM在亚马逊河、奥里诺科河、赞比西河流域区域与以上三种数据一致性较好，在华北平原、里海、红海区域一致性较差.这是因为亚马逊河、奥里诺科河、赞比西河流域均为季节性信号较强而且面积较大的区域，PM和DS300信号虽然会有泄漏误差，但大多数信号都泄漏在区域内；而红海区域虽然季节性信号占主要地位，但该区域面积较小且形状狭长，因此PM和DS300信号极易泄漏到区域外；华北平原和里海为趋势信号较强而季节性信号较弱区域，由于信号泄漏对趋势信号影响较大，因此PM和DS300在该区域与VCPM、CSRM、JPLM一致性较差.在育空区域，VCPM与JPLM一致性更好，而CSRM与PM、DS300一致性较好，表现出了较大的负增长趋势信号，这可能是阿拉斯加冰雪融化信号泄漏至该区域引起的.

图16 不同方法反演的2003-01—2014-11等效水高在不同区域的时间序列Fig.16 Time series of equivalent water height in different regions estimated by different methods from January 2003 to November 2014

表4 不同方法反演的2003-01—2014-11等效水高在不同区域的周年振幅(cm)Table 4 Annual amplitude of equivalent water height in different regions estimated by different methods from January 2003 to November 2014 (cm)

GLDAS2.1模型反映的振幅信号和长期趋势信号在多个区域与其他模型有较大差异，这是因为该模型主要反映的是浅层土壤水、地表积雪及植被冠层水等质量变化，而忽略了深层土壤水、地下水和地表湖泊等影响(詹金刚等, 2011).从华北平原、印度河流域的时间序列可以看出，GLDAS2.1时间序列趋势和其他模型相反，表明这两个区域存在地下水亏损的现象(Chen et al., 2014;魏伟等, 2020).

3.3 冰冻圈信号对比

由于全球气候变暖导致南极和格陵兰岛存在巨大的冰雪消融趋势信号，PM方法和传统滤波方法在这两个区域均存在巨大的泄漏误差，因此本文利用各数据在南极和格陵兰岛进行对比分析以验证VCPM消除泄漏误差的能力.由于CSRM和JPLM均选用了ICE6G_D(VM05)模型进行GIA改正，因此本文其他模型也进行了同样的GIA改正，其中DS300扣除的GIA模型经过了与之相同的截断与滤波处理(Peltier et al., 2018).经计算，ICE6G_D(VM05)在南极的趋势改正值约为98.8 Gt/a，在格陵兰岛改正值约为2.5 Gt/a.

图17为不同方法反演的2003-01—2014-11南极区域质量变化长期趋势，DS300由于滤波平滑作用存在很大泄漏误差，其趋势信号也被极大削弱；PM泄漏误差也较为严重，反映出利用单位矩阵作为正则化矩阵在处理泄漏误差方面的缺陷.VCPM由于先验方差的约束作用极大程度地消除了南极区域的泄漏误差，同时提高了南极区域质量变化的空间分辨率和信号精度，从南极半岛、毛德皇后地、恩比德地、威尔克斯地等区域反映的趋势信号来看，VCPM反演结果与CSRM和JPLM一致性极好.图18为不同数据计算的南极区域质量变化时间序列，从整体来看VCPM、CSRM、JPLM三种数据时间序列基本一致，VCPM与CSRM的相关系数为0.997，与JPLM的相关系数为0.991，均有很好的相关性.由于DS300和PM两种数据在南极区域泄漏误差很大，因此对应时间序列与以上三种数据差异较大.通过对时间序列进行最小二乘拟合，得到了如表5所示的南极质量变化速度与加速度，其中VCPM模型反演的南极冰盖融化速率为-139.8 Gt/a，与CSRM(-130.4 Gt/a)和JPLM(-123.1 Gt/a)的融化速率差值均在20 Gt/a以内.三种方法均探测到南极冰川呈现加速融化的趋势，VCPM反演的加速度最大(-9.8 Gt/a2)，CSRM反演的加速度最小(-9.1 Gt/a2)，其差值也在8%以内，可能是因为VCPM利用迭代方法充分提取了观测值中所包含的信号.

表5 不同方法反演的2003-01—2014-11南极质量变化速度与加速度Table 5 Velocity and acceleration of mass variations in Antarctic estimated by different methods from January 2003 to November 2014

图19为不同方法反演的格陵兰岛质量变化长期趋势和加速度空间分布图，与南极区域类似，DS300和PM在格陵兰岛区域也存在很大泄漏误差，而VCPM、CSRM和JPLM均有效地去除了泄漏误差且趋势空间分布较为一致，三种数据均反映出格陵兰岛外围海岸线附近存在较大的质量亏损趋势，且以东南和西北区域最为严重，是因为这两个区域均存在流速较快的外流冰川；格陵兰岛北部、西北和东南部分质量损失趋势比之较小，是因为这三个区域外流冰川较少，冰雪融化在其质量损失中占主要作用(Rignot and Mouginot, 2012；Velicogna et al., 2014).三种数据在格陵兰岛中部地区均表现出了质量积累趋势，这是因为该地区地处高纬且海拔较高引起了降雪积累，但其空间分布和量级有所差异(Khan et al., 2015).图20为不同方法计算的格陵兰岛区域质量变化时间序列，表6给出了相应的不同方法反演的质量变化速度与加速度.VCPM、CSRM和JPLM三种数据的时间序列一致性较好，JPLM质量变化速率最大为-308.3 Gt/a，CSRM和VCPM分别为-297.8 Gt/a和-288.6 Gt/a，三者差距在7%以内，其中VCPM与Velicogna等(2020)结果(-261 Gt/a)更为接近；在加速度方面，三模型反演结果无明显差别.从以上分析可以看出利用质量块质量变化的先验方差对点质量方法进行合理约束可有效改正冰川区域的泄漏误差，提升反演地表质量变化的效果.

图17 不同方法反演的2003-01—2014-11南极等效水高长期趋势Fig.17 Trends of equivalent water height in Antarctic estimated by different methods from January 2003 to November 2014

图18 不同方法反演的2003-01—2014-11南极质量变化时间序列Fig.18 Time series of mass variations in Antarctic estimated by different methods from January 2003 to November 2014

图19 不同方法反演的2003-01—2014-11格陵兰岛等效水高长期趋势Fig.19 Trends of equivalent water height in Greenland estimated by different methods from January 2003 to November 2014

图20 不同方法反演的2003-01—2014-11格陵兰岛质量变化时间序列Fig.20 Time series of mass variations in Greenland estimated by different methods from January 2003 to November 2014

表6 不同方法反演的2003-01—2014-11格陵兰岛质量变化速度与加速度Table 6 Velocity and acceleration of mass variations in Greenland estimated by different methods from January 2003 to November 2014

4 结论

本文利用附有方差约束的径向点质量方法反演了2003-01—2014-11全球等效水高变化，结合PM、DS300、CSRM、JPLM和GLDAS2.1模型在多个尺度进行了对比分析，得出如下结论：

(1)PM在最佳正则化参数下的反演结果仍会受到条带误差影响，同时存在很大泄漏误差.本文设计的附有方差约束的径向点质量方法通过不同正则化矩阵迭代提取重力扰动观测值中所蕴含的质量变化信号，能够有效抑制海陆泄漏误差，消除条带误差的影响，同时充分提取观测信号中不同变化尺度的信号.

(2)在全球尺度的对比中，VCPM与CSRM的RMSE在全球和陆地尺度上最小，海洋区域仅略低于CSRM与JPLM的RMSE；VCPM、CSRM和JPLM反演的周年振幅和长期趋势信号在空间分布和信号强度上较为一致，均能反映出高精度高分辨率的地球各区域质量变化物理信号，表明VCPM在全球尺度上精度与Mascon数据相当.

(3)在各陆地流域和两个内陆海的时间序列对比中，VCPM与CSRM、JPLM时间序列一致性均较高，在扣除季节项和趋势项后，除红海之外的其他区域相关系数均能达到0.86以上，在以季节性信号为主的区域，相关系数可达0.98以上，红海区域的相关系数较低可能是该区域较小且较为狭长，导致各方法反演的质量变化空间分布不一致引起的；同时发现泄漏误差对于长期趋势信号影响更为明显，PM和DS300在长期趋势较大区域信号泄漏更大.

(4)在南极和格陵兰岛两大冰盖区域，VCPM能够有效抑制海陆泄漏误差，在南极区域VCPM、CSRM和JPLM反演的质量减小趋势分别为-139.8 Gt/a、-130.4 Gt/a和-123.4 Gt/a，在格陵兰岛区域分别为-288.6 Gt/a、-308.3 Gt/a和-297.8 Gt/a，VCPM在反演冰盖质量变化方面相对于DS300和PM方法有极大提升，表明本文设计的方差约束可以有效提高点质量方法反演地表质量变化的效果.

致谢感谢CSR提供的GRACE时变重力场数据和Mascon数据，JPL提供的Mascon数据，NASA(National Aeronautics and Space Administration，美国国家航空航天局)提供的GLDAS2.1数据.感谢审稿专家提供的宝贵意见.