摘要：引文格式：刘鹏, 段虎荣, 张成浩, 等. 融合GNSS和GRACE/GFO数据反演西南地区陆地水储量变化[J]. 测绘通报, 2024(10): 39-45.doi: 10.13474/j.cnki.11-2246.2024.1007.
本文内容来源于《测绘通报》2024年第10期，审图号：GS京（2024）1967号
融合GNSS和GRACE/GFO数据反演西南地区陆地水储量变化
刘鹏1, 段虎荣1,2, 张成浩1, 王金驰1, 梁文康1
1. 西安科技大学测绘科学与技术学院, 陕西 西安 710054;
2. 中国科学院精密测量科学与技术创新 研究院大地测量与地球动力学国家重点实验室, 湖北 武汉 430077
基金项目：大地测量与地球动力学国家重点实验室开放基金(SKLGED2022-5-2);陕西省自然科学基础研究计划(2021JQ562)
关键词：GNSS, GRACE/GFO, 融合反演, 陆地水储量变化, 中国西南地区
引文格式：刘鹏, 段虎荣, 张成浩, 等. 融合GNSS和GRACE/GFO数据反演西南地区陆地水储量变化[J]. 测绘通报, 2024(10): 39-45.doi: 10.13474/j.cnki.11-2246.2024.1007.
摘要
摘要：为解决单一使用GNSS或GRACE/GFO等技术反演区域陆地储水量变化存在的不足,本文借鉴分解-重构法的思想,提出了一种融合GNSS和GRACE/ GFO反演区域陆地储水量变化的方法。利用来自中国陆态网西南地区的90个GNSS测站时间序列数据,融合GNSS和GRACE/GFO反演了中国西南地区2012—2022年的陆地水储量变化,并通过空间分布分析对比其他融合方法及降水数据,验证了本文融合方法的性能。结果表明,融合反演结果在空间分布上综合了GNSS和GRACE/GFO的优点;在时间尺度上,GNSS、GRACE/GFO的反演结果均滞后于降水2个月,而融合反演结果仅滞后于降水1个月,且相应的相关系数也高于GNSS、GRACE/GFO反演结果。因此,本文方法为融合多源数据以获取更可靠的区域陆地水储量变化提供了一种参考。
正文
陆地蓄水(terrestrial water storage，TWS)是指陆地上存储的所有形式的水的总和，是全球和区域水循环的重要组成部分，对于深入研究气候变化和水循环信息具有重要意义。虽然水文模型可以较好地监测陆地水储量时空变化，但现有水文模型缺乏对所有水成分(如深层地下水等)的全面建模，这导致通过水文模型无法获取完整的水文信息[1]目前，重力恢复与气候试验(The gravity recovery and climate experiment，GRACE) 及其后续试验(GRACE follow-on，GFO)已被广泛应用于反演陆地水储量变化[2-4]。然而，GRACE/GFO的空间分辨率为330 km，时间分辨率为1月，这使其难以监测高频或局部区域的质量变化。为了消除高阶噪声和条带误差，通常会对GRACE/GFO原始数据进行滤波处理[3, 5]，这也会导致原始信息丢失。此外，GRACE与GFO衔接期间有着11个月的数据空缺，对如此长时间缺失数据的插值难免会引入人为误差。以上因素不可避免地阻碍了GRACE/GFO在实际区域水文监测中的应用。随着全球导航卫星系统(global navigation satellite system，GNSS)的快速发展，其定位精度可达到毫米级，能有效监测区域地表质量变化导致的位移变化[6]。目前，GNSS被广泛应用于监测区域陆地水储量变化，其结果与利用GRACE/GFO反演的陆地水储量变化有着较高的一致性[7-8]。然而，GNSS测站的密度直接影响其推算陆地水储量变化的可靠性和分辨率。此外，GNSS垂直位移时间序列中包含复杂的非水文信息(如孔隙弹性效应、板块构造运动等)，对反演结果有很大的影响，准确分离GNSS垂直位移中的非水文信息仍是一大难题[9]因此，GNSS和GRACE/GFO在反演区域陆地储水量变化方面具有互补优势[10]。文献[11]利用GNSS和GRACE/GFO数据融合反演了长江流域的陆地水储量变化，结果表明融合反演结果与长江流域降水数据的相关性优于仅用GNSS或GRACE/GFO反演的结果。与上述不同，本文借鉴分解-重构法(decomposition and reconstruction method，DRM)的思想，将仅GNSS、GRACE/GFO反演的陆地水储量变化分别分解为“背景场”和“剩余场”，再将GNSS的“剩余场”与GRACE/GFO的“背景场”相结合，以获取更可靠的反演结果。同时反演2012—2022年西南地区陆地水储量变化，再结合基于广义三角帽方法(the generalized three-cornered hat，GTCH)的融合结果及降水数据，验证该融合方法的性能。1 数据源1.1 GNSS时间序列数据
使用国家地震科学数据中心提供的GNSS站点垂直位移时间序列，获取了西南地区2012年1月—2022年12月共90个GNSS站点的单天解时间坐标序列，如图 1所示。该数据扣除了海洋潮、固体潮及极潮的影响。由于水文负荷主要体现在垂直位移上[12]，因此本文仅使用GNSS垂直位移进行反演。图 1 西南地区GNSS站点分布
图选项
为获得由水文负荷引起的季节性位移形变，利用德国地学研究中心(Geo Forschungs Zentrum，GFZ)提供的产品，扣除了非潮汐大气负荷(on-tidal atmospheric loading，NTAL)和非潮汐海洋负荷(non-tidal oceanic loading，NTOL)引起的垂直位移[13]。利用最小二乘法对GNSS时间序列进行拟合，以去除线性趋势和阶跃[14]。由于实际的GNSS数据难免有缺失，因此本文使用克里金-卡尔曼滤波插值方法对缺失数据进行补全[15]。图 2展示了部分GNSS测站(GZGY、QHBM、SCJU和YNYL)数据补全前后的垂直位移时间序列，其中蓝色线条为补全前；红色线条为补全后。可以看到，处理过程能有效地抑制高频噪声，并填补缺失数据。图 2 处理前后的GNSS垂直位移时间序列
图选项
1.2 GRACE/GFO数据
使用空间研究中心(Center for Space Research，CSR)发布的GRACE和GFO Level-2 RL06 GSM数据，时间跨度为2012年1月—2022年12月(下文中均用GRACE表示)。采用球谐系数对GSM月重力场模型进行正则化，使用GRACE提供的卫星激光测距数据替代C20和C30系数[16]。使用文献[5]提供的球谐系数代替一阶系数，利用ICE6G-D模型修正冰川效应。采用P4M6去条纹方法[3]和平均半径为300 km的高斯滤波校正南北条纹和残留的高阶噪声。此外，对于少量缺失的数据，采用线性插值方法补全。将GRACE在西南地区约36个伪观测站计算的陆地水储量变化转换为地表垂直位移，再与GNSS垂直位移数据融合反演得到西南地区陆地水储量变化。1.3 水文模型与降水数据
使用的水文模型和降水数据均来自全球陆地数据同化系统(Global Land Data Assimilation System，GLDAS)。水文模型采用空间分辨率为0.25°×0.25°，时间分辨率为1月的GLDAS NOAH 2.1数据，将模型中的冰雪、树冠水、四层土壤水(0~2 m)求和后再进行去均值处理即可得到地表水储量变化[17]根据弹性负荷理论[18]，GNSS垂直位移与陆地水储量变化的关系可表示为11为以等效水高(equivalent water height，EWH)表示的利用GNSS反演的陆地水储量变化值；e1。其中，ψ为质量变化源到待求点之间的角距离；Me、R分别为地球的质量和半径；hl为径向位移负荷勒夫数；P式中，L为拉普拉斯先验矩阵；β为正则化参数，采用交叉验证方法确定。2.2 分解-重构法反演水储量变化
将GRACE监测的陆地水储量变化转换为垂直位移，再与GNSS垂直位移相融合，反演区域陆地水储量变化。类似于式(1)，基于GRACE计算的地面点处垂直位移与陆地水储量变化的方程可表示为
式中，d22212(4)
式中，G1(1b)、G1(bl)分别表示采用基于1~b阶、b~l阶勒让德多项式计算的GNSS格林函数矩阵；同理，G2(1b)、G2(bl)为GRACE格林函数矩阵。则融合反演的x最终解可表示为
式中，x1(bl)为GNSS反演水储量变化的“剩余场”，表示为(G1(bl)TG1(bl)+βL)-1G1(bl)TG2(1b)+βL)-1G2(1b)T2.3 广义三角帽方法
广义三角帽方法可在缺少数据先验信息的情况下，评估多组数据集的不确定性，该方法常用于评估不同GRACE时变重力场模型的不确定度[20]。本文采用GTCH方法对GNSS和GRACE数据反演的TWS变化的不确定性进行评估，再根据评估的方差进行融合。对于不同数据，其加权融合的结果可以表示为式中，N为不同数据的个数；X为不同数据反演得到的TWS变化；w对应的权重，其计算方法为3 结果与讨论3.1 GNSS水文负荷位移的周年变化
西南地区多年平均垂直形变呈现出明显的空间差异。如图 3(a)所示，云南地区的GNSS垂直位移振幅最为显著，为6~10 mm，而在川东北、重庆、贵州部分地区，GNSS的垂直位移振幅较小，原因可能是西南地区处于东北信风带，导致东部降水减小，西部降水增多。
图 3 GNSS站点处水文负荷引起的位移振幅、相位
图选项
为揭示季节性储水量最大的时间，本文计算了GNSS时间序列的周年相位(如图 3(b)所示)。可以看出，西南地区GNSS站点处的周年相位在7—9月达到峰值，但南北差异较大，整体上由北向南逐渐增大，在云南南部达到最大值，原因可能是云南南部属热带季风性气候，降水更强。此外，云南南部主要为山地地形，导致降水长期在山谷地区聚集，因此达到峰值的时间相比其他地区更晚。
3.2 反演模型的不确定性
勒让德函数阶数的取值直接影响反演结果的可靠性，为得到最优的融合反演结果，本文计算了不同阶数下2012—2022年GNSS反演的陆地水储量变化(如图 4所示)。考虑GLDAS中缺乏地下水、未建模的地表水等，难以监测到高频的水波动信号，因此将解算的水储量变化作为“背景场”的参考。
图 4 不同阶数下模拟GNSS反演水储量周年振幅试验
图选项
由图 4可知，当阶数为25~30时，由于低阶的勒让德函数不能很好地描述地球表面的复杂形状和位移场变化，因此边界处失真现象明显。随着阶数的增大，边界效应减弱，当阶数为40~45时，信号逐渐趋于稳定，反映出更多的局部细节信息。当阶数为35时，GNSS反演的陆地水储量变化与模拟值(如图 4(a)所示)最为接近，因此，本文选择勒让德函数为35阶时的陆地水储量变化作为“背景场”。下文中，将GNSS、GRACE反演陆地水储量变化结果的“背景场”和“剩余场”记为“GNSS-background” “GNSS-residual”和“GRACE-background” “GRACE-residual”。
3.3 水储量变化的空间分布
确定勒让德函数的阶数取值后，可得出“GNSS-residual”和“GRACE-background”，如图 5所示。可以发现“GRACE-background”与“GNSS-background”(即图 4(d))在空间分布上大体一致，进一步表明本文选择勒让德函数阶数为35阶时的陆地水储量变化作为“背景场”是可靠的。图 5(a)是“GNSS-residual”，可以看出，空间变化明显的区域与图 3(a)中GNSS站点位移振幅较大的地区一致，因为其主要反映陆地水储量变化空间分布的局部细节信息。
图 5 GNSS反演陆地水储量周年振幅的“剩余场”及GRACE的“背景场”
图选项
为更好地评估融合反演方法的性能，本文还分别计算了单一GNSS、单一GRACE及基于GTCH融合方法反演的陆地水储量变化空间分布，如图 6所示。由图 6(a)—(b)可以看出，在空间分布上，GNSS、GRACE的反演结果具有较好的一致性，但GNSS结果表现出更强的振幅，最大约达300 mm。原因可能是GNSS位移中包含了其他非水文信息，对反演结果有着一定的影响，而GRACE由于时空分辨率低，难以监测到高频或局部区域的质量变化，因此表现出较低的振幅变化。图 6(c)显示，基于DRM的融合反演方法结果更接近GRACE而非GNSS的估算，表明融合反演中GRACE的贡献大于GNSS，但融合反演结果比GRACE的局部细节变化更强，原因是GNSS和GRACE观测在时空分辨率、空间覆盖和灵敏度上存在差异。对于GTCH融合方法的反演结果(如图 6(d)所示)，可以看出，GRACE对融合反演的贡献几乎占到了全部，融合反演结果与GRACE反演结果几乎一样，表明基于GTCH的融合方法无法体现出GNSS的优势。由此看出，基于DRM的融合反演结果包含了GNSS和GRACE的综合贡献，一方面，GRACE抑制了GNSS中其他非水文信号造成的影响；另一方面，GNSS在一定程度上弥补了GRACE对于区域质量变化不敏感的缺点。
图 6 西南地区2012—2022年陆地水储量周年振幅
图选项
为了进一步研究不同反演方法在监测局部地区陆地水储量变化方面的性能，选择了4个GNSS站点处(即图 1中的GZGY、SCLH、SCNC、YNLC)2012—2022年的陆地水储量变化时间序列，并计算了GNSS、GRACE反演结果与融合反演结果的相关系数，结果如图 7、表 1所示。由图 7可知，这3种反演结果在变化趋势上基本一致。由表 1可以看出，在GZGY和SCNC站点，GRACE与融合反演结果的相关系数高于GNSS，表明GRACE对融合反演的贡献大于GNSS，原因可能是这两个站点周围GNSS台站稀疏，GNSS监测陆地水储量变化的性能弱于GRACE；与之相反的是，由于SCLH和YNLC站点处GNSS台站密集，GNSS监测陆地水储量变化的性能强于GRACE，因此GNSS对融合反演的贡献大于GRACE。值得注意的是，在YNLC站点处，GNSS、GRACE与融合反演结果的相关系数均高于0.9，原因可能是YNLC站点处的陆地水储量变化幅度较大，GNSS和GRACE均能很好地监测到陆地水储量的变化。总体而言，融合反演结果在空间分布上综合了GNSS和GRACE的优点，一定程度上克服了单一技术反演的局限性。
图 7 西南地区2012—2022年在GZGY、SCLH、SCNC、YNLC站点处GNSS、GRACE及融合反演陆地水储量变化时间序列(蓝色线条为GNSS，黑色线条为GRACE，红色线条为融合反演)
图选项
3.4 时变水储量变化特征
本文研究区位于低纬度地区，陆地水储量的主要供应来源是降水[21]。为了研究西南地区陆地水储量变化的时间特征，分别计算了GNSS、GRACE及融合反演结果在区域内的平均EWH，并结合区域内平均降水量进行了综合对比(如图 8所示)。可以看出，3种反演结果表现出很强的一致性，均在每年的8—10月达到高峰，1—3月出现低谷，呈明显的季节性变化。对比降水数据可知，西南地区的降水主要集中在夏季，但这3种反演结果均滞后于降水，因为降水需要通过一系列复杂的水动力过程(如蒸散发、径流等)才会转换为部分陆地水储量。图 8 西南地区2012—2022年平均陆地水储量变化及降水时间序列
图选项
为更好地分析陆地水储量时间变化与降水量之间的关系，计算了GNSS、GRACE及融合反演的陆地水储量变化与不同滞后时间降水量之间的相关系数(见表 2)。可以看出，GNSS和GRACE的陆地水储量变化均滞后于降水2个月，而融合反演的陆地水储量变化仅滞后于降水1个月，且在无滞后时，融合反演结果与降水的相关系数也高于单独反演。总体而言，这3种反演结果与降水之间都表现出很强的相关性(最高相关系数均在0.8以上)，融合反演结果与降水之间的相关性最强。这表明相较于单独使用GNSS或GRACE反演，融合反演方法能更迅速且准确地监测到由降水引起的陆地水储量变化。
4 结论
本文提供了一种融合GNSS和GRACE反演区域陆地水储量变化的方法(即DRM)，并利用来自陆态网的90个GNSS站点时间序列数据和来自CSR GSM解的GRACE数据，反演了2012年1月—2022年12月西南地区陆地水储量变化。通过与其他融合方法(即GTCH)及降水数据的对比验证了DRM融合反演方法的性能。结论如下。
(1) 在陆地水储量变化空间分布上，GNSS、GRACE及基于DRM的融合反演结果具有较高的一致性。基于DRM的融合反演结果综合了GNSS和GRACE各自的优点，提高了反演结果的可靠性。
(2) GNSS、GRACE和基于DRM的融合反演结果具有一致的季节性变化。GNSS、GRACE反演结果与降水的最佳相关性出现在滞后降水2个月时，相关系数分别为0.84和0.82；而基于DRM的融合反演结果仅滞后于降水1个月，相关系数为0.89，高于仅GNSS和仅GRACE的反演结果。这表明基于DRM的融合反演方法对陆地水储量变化的响应比GNSS和GRACE更为敏捷和准确。
本文初步验证了DRM融合反演方法的性能，该方法将GNSS和GRACE反演陆地水储量变化的优势相结合，一定程度上弥补了各自的不足之处，为融合多源数据反演陆地水储量变化提供了一定的参考。
作者简介：刘鹏(1997—),男,硕士生,研究方向为卫星大地测量。E-mail:pliu126@126.com
通信作者：段虎荣,E-mail:duanhurong@126.com
资讯
来源：测绘学报