利用多元线性回归模型重构中国九大流域陆地水储量变化

杨鑫春 游 为 万祥禹 宋梦芝

1 西南交通大学地球科学与环境工程学院，成都市犀安路999号，611756 2 西南科技大学环境与资源学院，四川省绵阳市青龙大道中段59号，621010

GRACE卫星于2017-10退役，其提供的数据只更新至2017-06。但其继任者GRACE-FO卫星直到2018-05才发射升空，并于同年6月开始发布全球陆地水储量变化(terrestrial water storage changes, TWSC)信息。两代GRACE卫星数据之间存在11个月的空缺期，打破了GRACE/GRACE-FO监测TWSC的连续性，也会限制其进一步应用[1]。因此，填补11个月的数据空缺或重构连续的TWSC尤为重要。

部分学者采用机器学习方法来重构或预测连续的TWSC[1-2]。如Sun等[1]采用深度神经网络(deep neural network, DNN)、多变量季节性自回归整合移动平均模型(seasonal autoregressive integrated moving average with external variables, SARIMAX)和多元线性回归(multiple linear regression, MLR)3种机器学习方法，重构多种全球TWSC数据产品。结果表明，DNN性能略优于SARIMAX，明显优于MLR。Li等[2]联合多种统计分解、时间序列分解和机器学习方法重构全球TWSC信号，与采用整体信号重构不同，其首先分别重构TWSC时变模态上的季节项、年际变化项和残余项信号，然后将GRACE/GRACE-FO原有时变模态上的趋势项信号加回到重构的非趋势项信号，以生成重构的全频段时变模态信号，最后将GRACE/GRACE-FO原有空间模态乘以重构时变模型得到总的TWSC信号。同时，Li等[2]认为在联合使用统计分解方法时，相比于其他复杂的机器学习方法，虽然MLR方法在训练阶段表现并非最佳，但在测试阶段表现较好，因此MLR模型是一种可靠和稳健的重构方法。基于以上研究，本文将采用MLR方法，以降水、气温和水文模型模拟的陆地水储量数据作为预测参数，重构中国九大流域(边界数据由资源环境科学与数据中心(http:∥www.resdc.cn)提供)连续的TWSC；同时，对比分析整体信号重构、去趋势项信号重构和去趋势项去季节项信号重构等3种策略对重构数据质量的影响，以期提供最适合九大流域的重构策略。

1 数据处理

1.1 GRACE/GRACE-FO数据

GRACE/GRACE-FO观测数据产品具有2种不同的形式，即传统的球谐系数产品和最新的Mascon产品[1]。相较于传统的球谐系数产品，Mascon产品主要具有以下优势[3]：1)减少从陆地到海洋的泄漏误差，增加数据的信噪比；2)在处理过程中应用地球物理数据约束，几乎无经验后处理的平滑滤波要求，更方便非大地测量用户的使用。本实验采用CSR最新发布的RL06 v02 Mascon产品[4]，下文简称CSR-M。CSR-M空间分辨率为0.25°×0.25°，所选择的时间跨度为2002-04～2021-12，即共计183个月的GRACE数据、43个月的GRACE-FO数据以及11个月的数据空缺。需要指出的是，本文重点关注GRACE/GRACE-FO系统间11个月空缺数据的重构问题，因此对于因GRACE/GRACE-FO卫星仪器问题、校准活动、电池管理等原因造成的系统1～2个月的数据空缺采用简单线性插值方法进行补全[5]。

1.2 预测参数

降水和气温是陆地水储量变化的重要驱动因子[6]，故常作为重构TWSC的主要预测参数[1-2]。此外，Sun等[1]指出虽然水文模型模拟的TWSC通常不包括地下水和表面水，但其与GRACE TWSC存在强相关性，增加该变量作为预测参数所重构的TWSC性能优于只考虑降水和气温作为驱动参数的重构数据。基于此，本文采用ECMWF发布的月度气候再分析数据集ERA5-Land中降水、气温和TWSC(包括土壤水和雪水)数据作为预测参数[7]。考虑到3个预测参数与GRACE/GRACE-FO TWSC存在相位差[6]，选择时间窗口前3～0个月内的所有变量作为预测参数，即每个预测参数在2002-01～2021-12范围内将具有4个预测变量。同时，为了与CSR-M数据保持一致性，3个驱动参数的空间分辨率由0.1°×0.1°重新采样为0.25°×0.25°，并删除2004～2009年均值。

2 重构方法与策略

2.1 MLR模型

MLR模型是一种用于估计2个或多个自变量与1个因变量之间线性关系的统计方法。本文使用3个预测参数(降水、气温和水文模型模拟的TWSC)作为自变量，GRACE/GRACE-FO TWSC 作为因变量，其对应的MLR模型计算公式为：

yi=β0+β1xi1+β2xi2+β3xi3+ε

(1)

式中，yi为某个流域的GRACE/GRACE-FO TWSC 时间序列；xi1、xi2、xi3为对应流域内3个预测参数的时间序列；β0、β1、β2、β3为待估参数；ε为模型误差。本文使用183个月的GRACE数据，采用最小二乘方法估算MLR模型的待估参数，并采用43个月的GRACE-FO数据测试所确定MLR模型的可靠性。根据所确定的MLR模型，通过预测参数来重构2002-04～2021-12连续的TWSC。

2.2 重构策略

本文使用3种不同的重构策略。重构策略1为整体信号重构，是对GRACE/GRACE-FO TWSC整体信号使用式(1)直接重构生成连续的TWSC；重构策略2为去趋势项信号重构，是在使用MLR模型重构GRACE/GRACE-FO TWSC去长趋势项信号基础上加回GRACE/GRACE-FO TWSC原有的长趋势项信号以生成连续的TWSC；重构策略3为去趋势项去季节项信号重构，是在采用MLR模型重构GRACE/GRACE-FO TWSC去长趋势项去季节项信号后加回GRACE/GRACE-FO TWSC原有的长趋势项和季节项信号以生成连续的TWSC。

为获得信号的趋势项和季节项，采用最小二乘分解方法对GRACE/GRACE-FO TWSC 和预测参数的时间序列进行分解[3,8]：

Stotal=Strend+Sseasonal+Sresidual

(2)

式中，Stotal为总的原有信号；Strend为信号的长期线性趋势项；Sseasonal为信号的季节项；Sresidual为信号的余项，主要包含年际变化项、次季节项和噪声。需要说明的是，虽然11个月数据空缺会导致GRACE/GRACE-FO TWSC趋势项值存在一定偏差，但这种偏差对重构数据质量的影响可以忽略[2,8]。此外，GRACE/GRACE-FO TWSC季节项信号在重构时间范围内可认为并未发生实质性变化[2]，因此GRACE/GRACE-FO TWSC季节项信号计算也不受数据空缺的影响。

2.3 精度评定指标

本文采用皮尔逊相关系数(Pearson correlation coefficient，CC)、归一化均方根误差(normalized root mean square error，NRMSE)和纳什效率系数(Nash-Sutcliffe efficiency coefficient，NSE)等3种常用的精度指标来评定重构TWSC的性能，各指标详细计算公式见文献[1]。CC、NRMSE和NSE的取值范围分别为[-1,1]、[0,+∞)和(-∞,1]，CC和NSE数值越大、NRMSE数值越小，重构数据的性能越好。此外，采用重构的TWSC与GRACE/GRACE-FO TWSC位于测试阶段(2018-06～2021-12)的RMSE值作为重构数据的不确定性值[2]。

3 结果与讨论

图1为中国九大流域的CSR-M TWSC以及基于3种不同策略重构的TWSC，其相应的性能指标和不确定性值如图2和表1所示。基于策略1重构的TWSC在松花江辽河、黄河、淮河、东南诸河和珠江五个流域表现出良好的性能(CC/NRMSE/NSE值分别优于0.80/0.10/0.70)，但在海滦河、长江、西南诸河和内陆河四个流域表现出较差的性能。这种性能差异主要是因为本文只采用与气候变化有关的降水、气温和水文模型模拟的TWSC数据作为驱动参数，未考虑人类活动和冰川融化等其他影响因子。松花江辽河、黄河、淮河、东南诸河和珠江流域TWSC主要受气候变化影响，而海滦河、长江、西南诸河和内陆河流域TWSC除考虑气候变化影响外，还应考虑人类活动或冰川融化的影响[5]。此外，重构的TWSC在海滦河和西南诸河两个流域存在异常的不确定性，其对应数值分别高达13.55 cm和8.67 cm。因此，基于整体信号重构策略重构的TWSC不宜作为九大流域的最后重构结果。

图1 基于不同策略重构的中国九大流域TWSCFig.2 The reconstructed TWSC of nine major river basins in China based on different strategies

相较于策略1，策略2可显著提高九大流域重构数据的性能，特别是在海滦河、长江、西南诸河以及内陆河流域。例如在海滦河流域，重构数据的性能指标CC/NRMSE/NSE分别提高0.25/0.09/0.48，不确定性值降低9.71 cm。这种性能的提高得益于人类活动和冰川融化等因素主要影响GRACE/GRACE-FO TWSC趋势项信号的变化[2]，而去趋势项信号重构策略完全保留该部分信号。

相较于策略2，策略3进一步提升重构数据的质量，但提升幅度并不明显。重构数据质量的进一步提升是因为策略3完全包含原有的GRACE/GRACE-FO TWSC 趋势项和季节项信号，进而能最大限度地保留人类活动和冰川融化对TWSC的影响。此外还可以发现，基于策略3重构的TWSC在海滦河、淮河和东南诸河流域存在较大的不确定性，其值分别为3.82 cm、4.01 cm和3.66 cm，这主要受GRACE/GRACE-FO 空间分辨率(从低纬度约 200 000 km2到两极附近约 90 000 km2范围[5])限制，较小流域的TWSC可能会存在明显的信号泄露[9]。

图2 中国九大流域重构的TWSC性能指标对比Fig.2 Comparison of performance indexes of reconstructed TWSC of nine major river basins in China

表1 不同策略重构的中国九大流域TWSC的不确定性

根据以上分析，将基于策略3重构的TWSC作为九大流域最后的重构结果，表2为重构数据的性能指标。从表2可以看出，重构数据的性能在不同流域存在一定差异，这主要与各自流域的GRACE/GRACE-FO 数据信噪比以及预测参数与GRACE/GRACE-FO TWSC的相关性有关。例如在信噪比以及相关性高的西南诸河流域(图1(h)和表3)，重构数据的性能指标CC/NRMSE/NSE高达0.97/0.05/0.94，而在信噪比和相关性低的内陆河流域(图1(i)和表3)，重构数据的性能指标CC/NRMSE/NSE值仅有0.77/0.12/0.60。因此在对比不同流域的重构数据性能指标时，应综合考虑各流域的GRACE/GRACE-FO 数据信噪比以及预测参数与GRACE/GRACE-FO TWSC的相关性。

表2 基于策略3重构的中国九大流域TWSC性能

表3 预测参数与GRACE TWSC在中国九大流域的最大相关系数(时间窗口为前3～0个月)

4 结 语

本文利用MLR模型，采用3种不同的重构策略，重构中国九大流域的TWSC，得到以下结论：

1)基于整体信号重构策略的重构结果仅在TWSC主要受气候变化影响的流域(如松花江辽河、黄河、淮河、东南诸河和珠江流域)表现出良好的性能，而在人类活动或冰川融化(如海滦河、长江、西南诸河和内陆河流域)频繁的流域表现较差。

2)去趋势项信号重构策略的重构结果在九大流域均优于整体信号重构策略，在人类活动和冰川融化剧烈的流域优势更为明显。

3)去趋势项去季节项信号重构策略的重构结果在九大流域均稍优于去趋势项信号重构策略。

4)不同流域重构结果的性能还与该流域的GRACE/GRACE-FO数据信噪比以及预测参数(降水、气温、模型模拟的TWSC)与GRACE/GRACE-FO TWSC的相关性有关，信噪比和相关性高的流域(如西南诸河流域)其重构结果的性能优于信噪比和相关性低的流域(如内陆河流域)。