地球重力场及其时变效应反映了地球质量的空间分布及其迁移变化。全球性环境问题如海平面上升、冰川融化以及干旱等都与地球表层质量迁移紧密相关，研究地球表层质量迁移对监测全球环境和气候变化具有重要意义。随着GRACE卫星的成功发射，卫星重力测量技术为监测全球地表质量迁移提供了一种有效途径，目前已广泛应用于监测陆地水储量变化、冰川融化、海水质量变化等领域研究^[1-4]。

当前，利用GRACE卫星测量数据反演地表质量变化的方法按照反演原理不同，主要可分为两大类：一类是文献[5]建立的球谐系数法，采用GRACE卫星Level-2时变重力场模型，通过计算时变重力场模型相对于稳态重力场模型的球谐系数改变量反演得到质量变化，该方法计算简单应用也最广泛，但是也存在着反演结果时空分辨率低、南北条带噪声显著和信号泄露等缺陷^[6-8]；另一类方法为Mascon方法，采用GRACE卫星Level-1B数据，将研究区域划分为规则网格，视每个网格区域内具有相同质量变化，分离不同网格区域质量变化对重力场变化的影响，建立网格区域质量变化和星间距离变率之间的直接函数关系。研究表明Mascon方法反演结果具有更高的时空分辨率^[9-12]，但是与球谐系数法相比，Mascon方法计算相对复杂。

为减少反演数值计算复杂程度，国外学者基于Mascon方法原理，提出了径向点质量模型方法。该方法采用Level-2时变重力场模型数据，将空间点扰动重力向下延拓，建立与地表质量变化直接关系，进而反演得到质量变化。文献[13]建立了地球表层点质量变化引起空间点扰动重力变化的数学模型；文献[14]采用径向点质量模型方法反演了格陵兰区域冰盖质量变化，反演结果空间分辨率有所提高，其中病态问题采用Tikhonov正则化求解，正则化参数通过L曲线法确定；文献[15]建立了三维点质量模型方法，并通过模拟试验验证了方法的有效性；文献[16]采用径向点质量模型方法反演了冰岛区域冰川质量变化；文献[17]采用径向点质量模型方法反演了格陵兰区域冰盖质量变化，病态问题采用截断奇异值分解方法处理，并根据反演结果统计分析了反演参数的最优选取。

径向点质量模型方法本质上是重力数据的向下延拓，属于病态问题，解算时法方程组系数矩阵奇异，采用最小二乘法无法求得稳定解，通常采用正则化方法处理^[14]，正则化参数一般采用L曲线或广义交叉核实法(generalized cross-validation，GCV)确定^[18-20]。但是，L曲线法在解算中可能存在过度平滑问题，而GCV方法亦存在最小值点不明显的问题^[21]。针对传统正则化方法存在的不足，本文利用相邻网格区域质量变化不存在阶跃变化的特性，对相邻网格区域质量变化之间引入空间约束方程条件，构建虚拟观测方程组，由于引入空间约束条件后存在两类观测方程，而两类观测方程初始权比不确定，故引入赫尔默特方差分量估计方法自适应修正两类观测方程权比，使反演结果更加稳定和符合实际。

1 径向点质量模型方法基本原理 1.1 数学模型

地球重力场由正常重力场和扰动重力场两部分组成。地球表层质量变化会引起空间点的扰动重力径向分量值的改变。而扰动重力径向分量值的改变可由时变重力场模型球谐系数的改变量来表示

式中，G为引力常数；M为地球质量；a为地球平均半径；l和m分别为重力场模型的阶数和次数；λ、φ分别为地心经度和地心余纬；r=a+h为空间点到地心的向径，h为空间点到地面的高度；l_max是重力场球谐模型截断阶数；{ΔC_lm, ΔS_lm}={C_lm－C_lm，S_lm－S_lm}是时变重力场模型{C_lm, S_lm}相对于参考重力场模型{C_lm, S_lm}的球谐系数改变量；k_l是l阶负荷勒夫数。

如图 1所示，地球表面一点P质量变化δm_j引起的空间一点S的扰动重力径向分量改变值δg_ρ^ij为

式中，d_ij是空间点S和地球表面点P的距离。根据图 1几何关系，可得地球表层质量变化引起的空间点S的扰动重力径向分量改变值为^[12]

式中，n是研究区域质量异常的点的个数；r_i为空间点S到地心的向径。进一步，可将待求质量变化划为等效水柱高形式。设δm_j=ρ×s_j×Δh_j，ρ为水密度，λ_j1、λ_j2为研究区域经度范围，φ_j1、φ_j2为研究区域纬度范围，s_j为研究区域面积，具体形式为

于是，可将式(3)所求区域质量变化换算成等效水柱高度Δh_j形式

联立式(1)和式(5)即得径向点质量模型方法数学模型，其中i=1、2、…、k，为空间离散点的个数。

1.2 参数估计

联立式(1)和式(5)，整理后得到直接观测方程组

式中，e₁为误差向量；y₁=[δg_ρⁱ]_k×1为观测向量可按式(1)计算求得；x=[Δh_j]_n×1是待求量；矩阵A₁(k×n)为线性化后系数矩阵，可按式(5)右端计算得到，权阵P₁简化为单位阵I^{[14, 16]}。由于点质量模型方法解算时法方程组系数矩阵A₁^TP₁A₁奇异^[14]，最小二乘解不可靠。通常采用Tikhonov正则化方法得到稳定解

式中，α≥0是正则化参数。从式(7)可以看出，待估参数解与正则化参数密切相关，求解的关键在于如何选择合适的正则化参数。

2 附加空间约束的径向点质量模型方法 2.1 附加空间约束

病态问题可通过增加约束信息得到适定的解，这类约束信息可以是待估参数的线性等式约束，也可是待估参数的二次范数约束。Tikhonov正则化方法增加待估参数间二次范数约束，通过正则化因子使观测残差和参数范数间达到平衡，从而求得稳定正则化解^[20]。鉴于L曲线法在解算中可能存在过度平滑问题，而GCV方法亦存在最小值点不明显的问题。为避免此类情况发生，针对径向点质量模型方法，从实际角度考虑，空间上临近网格区域受相似环境影响，相邻网格内质量变化在数值上不应存在阶跃性变化。因此，可对不同网格间质量变化构建空间约束条件，具体形式如下^{[10, 22]}

即认为任意两个网格区域质量变化Δh_n和Δh_k(n≠k)相等，而相等的“可能性”满足一个与空间距离有关的指数函数^[23]，如式(9)所示，式中d_nk是相应的两个网格中心的空间距离；D是相关距离，表征空间约束程度，与研究区域网格划分以及时变重力场模型数据分辨率有关，一般取值为网格宽度^[10]。从式(9)可知，两网格间距离越小，其质量变化在数值上接近的可能性越大；反之差异越大。

将空间约束条件式改写成观测值为0的虚拟观测方程式(10)，虚拟观测方程对应的权值为W_nk。经整理后，得到虚拟观测方程组式(11)，其中y₂=[0 … 0]^T，e₂为误差向量。系数矩阵A₂具体形式如式(12)

虚拟观测方程组权阵为对角阵P₂，对角线元素由式(9)计算，具体形式如下

2.2 赫尔默特方差分量估计

联立式(6)和式(11)即可解算得到反演结果。显然，增加空间约束条件后，反演解算实质上是联合求解两类观测方程组，但其初始权比往往不够合理，直接采用广义最小二乘准则得到的解不准确。在大地测量数据处理中，在初始权比不准确情况下，通常采用方差分量估计确定不同类型观测值的权^[24]。为求得更准确和稳定的解，本文引入验后赫尔默特方差分量估计自适应确定两类观测方程组的权比，将病态问题转化为两类观测值的测量平差问题^[25]。

赫尔默特方法分量估计是根据各次平差得到改正数向量估计方差分量，修正平差所用的权值，然后再以修正后的权值重新平差，直至趋近到权比合理为止。联合式(6)和式(11)，即整体观测方程组式(14)，设、分别是两类观测的方差分量估值

对式(14)进行赫尔默特方差分量估计，经过第J次迭代有

式中

通常两类观测方程权方差初值不等，通过若干次迭代修正后可得到

于是得到两类观测方程的权比为

再根据广义最小二乘准则，可得待估参数解

对比式(7)和式(19)可知，赫尔默特方差分量估计得到的两类观测方程权比与Tikhonov正则化因子的作用类似，都是调节最小二乘估计和稳定泛函约束间的平衡以提高待估参数解的可靠性^[16]。

3 区域陆地水储量变化反演结果与分析 3.1 试验数据及处理过程

由于南美大陆河流众多，水量丰富，拥有亚马孙河、巴拉那河、奥里诺科河等河流，特别是亚马孙河是世界上流域面积最广、流量最大的河流。故选择南美大陆区域(30°S—10°N，35°W—80°W)为研究区域，采用GFZ(Geo Forschungd Zentrum)发布Realease 05版本2006年1月—2008年12月共36个月的GRACE卫星时变重力场模型数据反演该区域陆地水储量变化，以验证附加空间约束的径向点质量模型方法的有效性。考虑到计算区域的边缘效应，实际选取计算区域比研究区域扩大了5°，计算时地面区域划分为边长为2°的正方形网格，空中扰动重力计算划分为边长为1°的正方形网格。GRACE时变重力场模型误差随阶数增大而增大，因此反演时球谐系数截断至60阶^[26]。由于时变重力场模型C₂₀项精度较差，反演时采用卫星激光测距(SLR)观测的C₂₀项替换。试验中，首先对36个月的时变重力场模型求平均值作为稳态背景场，然后计算每个月时变重力场模型相对于稳态背景场的球谐系数改变量，再由每个月时变重力场模型球谐系数改变量根据第2节所述方法对研究区域陆地水储量变化进行反演计算，构造空间约束条件时只对直接相邻的不同网格之间建立空间约束方程。

3.2 反演结果验证与分析

为了方便描述，附加空间约束的径向点质量模型方法记为PM-VCE；无空间约束条件的径向点质量模型方法记为PM-L，采用Tikhonov正则化方法，其中正则化参数采用L曲线方法求解；球谐系数法记为SH，采用400 km半径的扇形滤波；GLDAS(global land data assimilation system)水文模型结果记为GLDAS。图 2所示为4种方法所得南美大陆区域陆地水储量变化反演结果，第1列至第3列分别为2006年2月、2007年2月和2008年2月反演结果。

从图 2结果可以看出，4种方法所得到的2006年2月、2007年2月以及2008年2月的南美大陆区域陆地水储量变化反演结果在整体上具有相似性。亚马孙流域的巴西中西部、玻利维亚北部以及秘鲁东部区域，陆地水储量变化呈明显增大变化；而在奥里诺科河流域的巴西北部、哥伦比亚东部和委内瑞拉南部地区陆地水储量变化呈显著减小变化。PM-VCE方法反演结果与PM-L、SH方法反演结果在空间分布和数值上具有很好的一致性，水储量变化最大和最小的区域一致。与GLDAS水文模型结果相比，PM-VCE所得结果与GLDAS模型结果在变化趋势和空间分布上也具有较好的一致性。由于GLDAS水文模型是以雷达和微波辐射计资料为主，主要反映了浅层地表土壤湿度及地表积雪随时间变化，忽略了深层土壤水、地下水和地表湖泊等影响，故GLDAS模型结果与GRACE时变重力场反演结果在数值上有所差异^[27]。图 2南美大陆区域陆地水储量变化反演结果验证了PM-VCE方法的正确性。

表 1对2006年2月、2007年2月以及2008年2月的南美大陆区域陆地水储量变化PM-VCE方法反演结果进行了统计。由表 1可知，增加空间约束条件后，采用赫尔默特方差分量估计得到两类观测方程权比与Tikhonov正则化因子作用类似，调节了最小二乘和稳定泛函约束间的平衡，使得法方程组条件数明显减小、病态性减弱，进而得到待估参数的稳定解。

对2006年1月至2008年12月间南美大陆区域陆地水储量变化反演结果时间序列进行了统计分析。图 3所示为不同方法所得2006—2008年南美大陆区域陆地水储量变化反演结果时间序列。从图 3可以看出，陆地水储量变化呈现明显年周期变化，PM-VCE、PM-L和SH方法所得反演结果时间序列变化特征相符一致，并与GLDAS水文模型结果在变化整体趋势上也基本一致，其中SH方法反演结果振幅最大。

表 2为4种方法所得反演结果时间序列年周期特征，可看出，GLDAS模型结果时间序列周年振幅最小，这是因为GLDAS模型是浅层土壤水含量模型结果；PM-VCE方法和PM-L方法结果时间序列周年振幅一致，SH方法结果时间序列周年振幅最大，其中PM-VCE方法结果时间序列周年振幅最接近GLDAS模型结果；并且，SH、PM-L以及PM-VCE所得时间序列周年相位结果基本一致，GLDAS模型结果相位与其他3种方法结果相差大约30°，图 3时间序列结果也反映了该相位差异。

表 3为4种方法所得反演结果时间序列信号相关性统计。从结果可见，PM-VCE方法与PM-L、SH方法时间序列结果相关系数相一致，均为0.99。而PM-VCE、PM-L以及SH方法结果与GLDAS模型结果相关性基本一致，相关系数分别为0.77、0.78和0.79。

4 结论

针对径向点质量模型方法解算中的病态问题，本文从实际角度考虑构建空间约束条件，对不同网格区域质量变化建立空间约束的虚拟观测方程组，联合原始观测方程进行反演求解，并采用赫尔默特方差分量估计自适应修正两类不同观测方程权比，使得反演结果更符合实际和稳定。根据南美大陆区域陆地水储量变化反演试验结果可知：

(1) 引入空间约束条件后，增加了待估参数的约束信息，降低了法方程组病态程度，法方程组条件数明显减小。

(2) 附加空间约束的径向点质量模型方法与球谐系数法、GLDAS模型所得南美大陆区域陆地水储量变化反演结果在空间分布、数值变化上整体相一致，验证了本文方法的正确性，表明该方法可有效应用于区域地表质量变化反演，为利用卫星重力监测地表质量变化提供了一种新途径。

(3) 针对径向点质量模型方法，考虑到Tikhonov正则化中正则化参数确定方法的不稳定性，本文通过引入空间约束条件构建虚拟观测方程组，并采用赫尔默特方差分量估计方法自适应修正虚拟观测方程组和原始观测方程组权比，避免了L曲线等方法的缺陷。

当前GRACE卫星即将结束使命任务，未来随着GRACE Follow-on计划的实施，星间测距精度将得到大幅度提高，时变重力场模型精度和空间分辨率将得到有效提升。考虑相关距离、网格划分以及平差模型等参数组合对附加空间约束的径向点质量模型方法的影响，后续将进一步研究并改进解算方法，以期提高区域地表质量变化反演结果精度和空间分辨率。

致谢: 感谢中国科学院测量与地球物理研究所冯伟副研究员、杨帆博士，德国洪堡学者陈秋杰博士对本文提出的宝贵意见。