南方7省市陆地水储量动态变化主控因素研究

张昌伟，常 勇，罗 跃

（1. 河海大学地球科学与工程学院，江苏 南京 210098； 2. 东华理工大学水资源与环境工程学院，江西 南昌 330013）

0 前 言

陆地水储量（Terrestrial Water Storage， TWS）是储存在地表以及地下的全部水分（积雪、冰川、土壤水、地下水、河流湖泊水以及生物水）等［1］，是陆地和全球水循环的重要组成部分［2］，也是人类赖以生存和发展的基础。区域水资源的丰富程度一定程度上制约着当地农业和工业的发展，水储量的剧烈变化也会引发当地干旱或洪涝灾害，对当地经济和社会的发展具有重要影响。

传统对陆地水资源量的监测主要基于点状测量，受限于有限的监测点分布和复杂的地下状况，很难获得区域水储量的动态变化监测数据。自2002 年（Gravity Recovery and Climate Experiment）GRACE重力卫星发射以来，基于陆地区域重力场提取的陆地水储量距平（Terrestrial Water Storage Anomaly，TWSA），反映了地球陆地垂直方向上水储量的异常变化［3］，能有效避免传统方法的缺陷，为监测全球或区域陆地水储量变化提供了一种新的手段［4-6］。近年来，基于GRACE 卫星获取的TWSA 被广泛用于监测全球或区域尺度陆地水储量的动态变化，在认识不同地区水资源变化趋势以及影响因素方面展现出巨大的潜力。

气候变化和人类活动为影响TWSA 动态变化的主要因素［7-10］，一方面气候变化能改变水资源供给和排泄量，是决定区域水资源量丰富程度的关键因素；另一方面人类活动，如农业灌溉、修建水库、跨流域调水等，也会改变区域水循环过程，进一步影响TWSA的动态变化。在当前全球气候逐渐变暖和人类活动不断加强的背景下，分析区域TWSA 的动态变化及其主要控制影响因素对于正确认识未来TWSA的变化趋势和制定合理的水资源管理策略都具有重要意义。尽管目前开展了大量有关TWSA 变化方面的研究，但大部分研究通常以大型流域为基本研究单元［5，11］，一定程度上忽略了流域内部的差异。同时气候因子与TWSA 变化之间的关系也多采用简单的相关性分析，缺少定量分析气候因子对TWSA 变化的贡献，进而评估人类活动对TWSA 的影响程度。鉴于此，选取南方7 省市（云南省、四川省、贵州省、广西壮族自治区、重庆市、湖北省、湖南省），在网格尺度上分析TWSA 的动态变化趋势，在研究TWSA 与不同气候指标之间相关性的基础上，进一步采用偏最小二乘回归模型定量评价气候因子对TWSA 动态变化的贡献，分析南方7 省市不同区域TWSA 动态变化的影响因素，为该区域水资源管理和保护提供一定的依据。

1 研究区域、数据与方法

1.1 研究区概况

选取南方7省市为研究区，包括云南省、四川省、贵州省、广西壮族自治区、重庆市、湖北省、湖南省。这些省市区绝大部分地区位于亚热带季风气候，仅云南省南部局部区域属于热带季风气候，四川省西北部分地区属高山高原气候。研究区内地形起伏大，从向西向东逐渐由高原、山地逐渐转变为丘陵、平原，降雨量丰富，年降雨量从东南向西北方向逐渐降低。除四川省西北部地区外，研究区绝大部分地区多年平均降雨量在800 mm以上，局部地区最高年降雨量可达2 000 mm，区内仅局部山区存在常年冰川和积雪。

整个研究区涵盖了长江流域和珠江流域绝大部分上中游区域，水资源量丰富，基本无大面积地下水超采情况，但由于受亚热带季风气候的影响，降雨时空分布极不均匀，降雨多集中在4-9 月。这些省市也是我国岩溶集中分布区域［12］，地下水水位往往暴涨暴跌［13］，该地区也是我国干旱和洪涝频发地区［14，15］。研究区另一特点是水电站众多，其中大型的水电站包括位于长江流域的三峡水电站、珠江流域的龙滩水电站以及2021 年开始发电的白鹤滩水电站，其中龙滩水电站库容达273亿m3；三峡水电站库容393 亿m3。这些水电站的库容巨大，对于区域水资源量具有巨大的调蓄作用。研究这些省市陆地水储量的动态变化及其主要影响因素对于该地区水资源综合管理和预测未来气候下水资源的动态变化具有重要意义。

1.2 数据来源

研究采用的TWSA 数据来自美国德克萨斯大学（Center for Space Research，CSR）空间研究所最新公布的第二版的CSR_GRACE_GRACE-FO_RL06 重力场模型数据［16］，该产品有效克服了传统球谐系数因滤波而造成的泄露误差，相对于前期传统的球谐系数产品表现出更高的精度［17］。由于卫星2002 年初进行轨道校准以及2016年底电池老化而导致的数据误差，使得TWSA 不确定性小幅上升，其余时间皆稳定在30 mm 左右。该数据的原始精度为1°×1°，整个研究区共含有160个网格。选取2002 年4 月至2017 年6 月期间GRACE 数据作为水储量动态数据，由于GRACE 卫星前期调试以及运行期间传感器性能下降和功能不足问题，存在20 个月的数据缺失，采用李夫鹏重建的GRACE数据对缺失数据进行填补［18］。

为分析气候因素对TWSA 动态的影响，选取国家青藏高原科学数据中心发布的中国区域地面气候要素驱动数据集［19］，该数据集包含近地面气温、近地面气压、近地面空气比湿、近地面全风速、地面向下短波辐射、地面向下长波辐射、地面降水量共7个气候要素的月平均值，以下分别简称为气温、气压、比湿、风速、短波、长波以及降雨。为检查该数据集的精度，初步将研究区中该数据集1981 年至2010 年期间降雨、气温、比湿数据与中国气候数据网（http：∕∕www.nmic.cn∕）公开发布的中国气候标准值月值数据中的相关网格数据进行对比，发现每个网格的相关系数均达到0.99 以上，显示出该数据集在研究区较高的精度。除此之外，本文还选取由北京师范大学发布的全球陆表特征参量产品中的叶面积指数用于反应陆地植被的变化［20］。在这8个指标中，降雨为陆地水储量的主要输入量，而其他指标主要影响陆面蒸散发。这8种指标被选为陆地水储量的潜在气候驱动因子，均选取GRACE 数据的相同时段（2002.4-2017.6），同时采用简单平均法尺度粗化至1°×1°的经纬度网格月数据，与GRACE数据的空间和时间精度保持一致。

1.3 研究方法

采用Theil-Sen Median 方法分析每个网格陆地水储量的变化趋势和程度，对于时间序列（x1，x2，…，xn），斜率计算公式可表示为：

式中：Β为斜率；xi和xj分别为i、j时刻的TWSA 值；Median 为取中值函数。

采用基于局部回归分析（locally weighted regression，LOESS）算法的季节趋势分解（Seasonal-Trend-Loss，STL），对TWSA 进行分解，分别计算季节项与趋势项相对于总信号方差的贡献占比。STL是以鲁棒局部加权回归作为平滑方法的时间序列分解方法，分解周期性数据效果较好，适用于GRACE 数据［21］。STL可表示为：

式中：Xo为原始的时间序列；Xlong代表时间序列的长期趋势；Xsea代表时间序列中所存在的周期性变化规律；Xres为残差项。

互相关分析被用于研究不同气候指标与TWSA之间的相关性程度，前期的大量研究表明TWSA 的动态变化一般滞后于气候指标的变化，且滞后时间对于不同的气候指标或在不同的地区均存在差异。考虑到气候指标的周期性，结合袁瑞强等在全球尺度上分析TWS 与气候要素的时间滞后皆为2 或3 个月［22］，选取每个网格滞后时间3个月之内的相关系数的最大值表征该网格TWSA与气候指标之间的最大相关性。每个网格相关系数的计算方法如下：

式中：k为滞后月数；n为数据个数；Gt和Qt为TWSA 和气候指标的时间序列数据；σG和σQ分别为TWSA 和气候指标的方差。选择5%的显著性水平对相关系数进行显著性检验。

气候变化与人类活动是影响TWSA 变化的主要因素，由于人类活动在区域尺度上对TWSA 的影响很难进行统计和分析，通过定量评估气候指标对TWSA 动态变化的预测程度，可以在一定程度上间接了解人类活动对TWSA动态的影响。考虑到不同气候指标之间存在很强的相关性，选取偏最小二乘回归模型定量评估气候指标对TWSA动态的解释程度。偏最小二乘回归模型（Partial Least Squares Regression， PLSR）是多元线型回归模型的扩展，可以同时实现多元线型回归和主成分分析过程，特别适合当因变量比自变量有更多的变量，以及自变量存在多重共线性的情况。PLSR 算法主要包括以下步骤：首先，获取自变量的气候要素矩阵和因变量的TWSA 数据矩阵，进行标准化处理，得到自变量矩阵X和因变量矩阵Y。结合典型关联分析以及主成分分析思想，求出第一对主成分，并分别建立Y和X对第一对主成分的回归模型，随后用该方程的残差代替X和Y重复以上步骤，最终可建立以下偏最小二乘回归方程：

式中：ωk为主成分的轴向量；βk为Y对因变量主成分回归建模的回归系数；k即为主成分个数。

对于每个网格，TWSA 作为PLSR 模型因变量，考虑最优滞后时间的8 个气候指标作为模型的因变量，采用决定系数R2来评估模型的模拟效果和留一法交叉验证来最终确定最优模型结构。除采用PLSR 模型评估气候指标对TWSA 的预测能力外，PLSR 模型中每个气候指标的回归系数描述了PLSR 中的每个分量依赖于原始变量的权重，可用于评估不同气候指标对TWSA动态变化的贡献程度。

除采用互相关分析和PLS模型在月尺度上分析每个单元格不同气候指标与陆地水储量变化的关系外，研究还将每个气候指标和陆地水储量做13点滑动平均，一定程度上消除各指标的季节变化，分析它们在年尺度上的关系。

2 研究结果

2.1 TWSA变化的趋势以及主导模式

2002 年至2017 年间，TWSA 在南方七省市中，75%的网格多年平均水储量距平呈现明显上升趋势，在重庆、贵州和湖南省交界处上升幅度最大，平均增幅可达每年0.9 cm。在四川以及云南西部，TWSA则呈现出缓慢下降的趋势，最大下降幅度为每年0.3 cm左右。

基于STL 法，将TWSA 分解为趋势项、季节项以及残差项，通过比较每个分量与总信号的方差，可以评估这3 个分量的相对贡献。图2展示了研究区TWSA 趋势项与季节项方差的相对贡献的空间分布，展现出明显的区域规律。云南、四川中南部、广西南部以及湖南湖北二省的东部季节项的方差贡献占比皆大于50%，云南中西部占比达到最大为78.9%；趋势项的方差占比则普遍小于35%，云南中西部最小占比为11.9%，反映出TWSA在这些区域以周期性季节变化为主。四川北部、贵州、重庆以及湖南湖北二省的西部季节项的方差贡献占比皆小于50%，重庆、贵州、湖北和湖南四省市交界处占比最小达到26.1%；趋势项的方差占比则普遍在30%～50%之间，四省市交界处这一区域最大贡献可达48.5%。这表明TWSA 在这些区域受长期趋势和季节变化的共同影响。

2.2 气候指标与TWSA相关性分析

当考虑研究区气候指标与TWSA 最大相关时，月尺度下80%以上的网格TWSA 与气候指标之间均存在2～3个月的滞后时间，仅气压与TWSA 的滞后时间无明显规律。年尺度下，80%以上的网格TWSA 与降雨、湿度和短波之间仍存在1～3 个月的滞后，与其他指标基本无滞后。图3 显示了不同时间尺度下研究区内不同气候指标与TWSA 之间最大的相关系数，图中仅显示相关系数显著性水平大于5%的网格。

在月尺度上，气温、长波、降雨、叶面指数和比湿5种指标在研究区呈现相似的相关系数空间分布特征，除在重庆市、贵州省、湖南省和湖北省交界处以及四川西北部15个网格存在相关系数低于0.5 外，其他区域5 种指标与TWSA 均较高的相关性，最大相关系数可达0.95。在这5 种指标中，比湿、气温、叶面积指数与TWSA 的相关系数要高于其他2 个指标，降雨并未显示比其他指标更强相关性，仅在四省市交界处的相关性高于其他气候指标。短波、风速和气压与TWSA 的相关关系明显弱于其他指标。

当所有气候指标和TWSA 采取13 点滑动平均一定程度上消除季节变化后，从年尺度分析TWSA 与气候因素的互相关关系，不同气候指标与TWSA互相关程度出现了明显变化。比湿、长波以及叶面积指数总体上仍呈现明显正相关，但与月尺度的相关系数相比，对应网格会降低0.1～0.5，其中长波与叶面指数与TWSA 相关系数的下降最为明显，最大相关系数仅为0.71。在这四种指标中，降雨在绝大部分网格均呈现出与TWSA 最高的相关性，这也表明在年时间尺度上，TWSA在研究区80%以上的区域与降雨的关系最密切。对于其他指标，气温与TWSA 的相关关系最弱，45%网格没有通过显著性检验。90%网格的短波呈现负相关关系，气压的空间分布与之类似，负相关网格占据65%左右，但在贵州、湖南以及研究区域北部仍有正相关关系分布。风速整体上未呈现明显的规律，正、负相关都有所体现，在研究区中部和东部地区与TWSA的相关程度较高，最大相关系数值可达0.84。

从TWSA 和8 种不同指标的相关性分析可以发现，在月时间尺度上，除研究区中部区域外，TWSA与多种指标均呈很强的相关关系，而在年时间尺度上，仅降雨与比湿与TWSA呈现较强的相关性，表明不同时间尺度下控制TWSA 的动态变化因素存在一定的差异。

2.3 偏最小二乘回归模型拟合结果

图4 为PLSR 模型评价结果，其中图4（b）和图4（d）仅显示了每个网格贡献前三的气候指标，圆圈大小表示贡献度的大小。从图4中可以看出，月尺度下PLSR模型对研究区西部区域（整个云南省、四川省中南部）和东南部区域（广西南部和湖南东部）的模拟结果较高，决定系数值R2均大于0.60，最大决定系数达0.92。与月尺度下互相关分析结果类似，PLSR模型对研究区北部和中部区域的拟合结果较差，最低R2值仅为0.17，表明在这些区域仅采用气候指标不能很好的解释TWSA的月动态变化。进一步分析每个气候指标对陆地水储量变化的贡献，在65%网格中，比湿均为贡献度最大的指标，最大贡献度可达48.3%，同时气温在70%的网格尤其是研究区东部均在贡献前三的指标中，而降雨仅在40%网格中进入贡献前三。除此之外，在研究区南部和四川省，叶面积指数和长波对TWSA的变化起主要贡献，其余气候指标仅在局部网格中进入贡献前三。

相对于月尺度的模拟结果，年尺度下PLSR 模型对云南省内TWSA 预测的决定系数略微降低，但对月尺度下预测结果较差的区域拟合结果明显提高。除研究区三峡水库附近的网格外，绝大部分网格的决定系数R2均大于0.6，研究区中部部分网格R2值可达0.91，表明在研究区大部分区域TWSA 的60%以上年尺度变化可以由气候指标的波动来解释。对比月尺度下不同指标对TWSA的贡献，网格中降雨出现的频率明显增加，且降雨对TWSA 变化的贡献通常排名第一或第二，研究区内75%以上的地区，降雨仍是影响TWSA 变化的重要气候因素。除降雨外，比湿和气温也仍为半数网格中贡献前三的气候指标。在研究区中部及东部地区的网格中，风速为贡献排名第一的气候指标，其他气候指标仅在局部网格为贡献前三的指标。

3 讨 论

互相关分析结果指出除研究区中部区域外，TWSA 月动态变化与研究区内多个气候指标（除风速、短波辐射和气压外）均具有密切的相关性，结合TWSA变化主要模式分析结果，这均与这些气候指标和TWSA在亚热带气候影响下均呈现出明显的周期性季节变化有关。图5（b）为云南省一个网格（中心坐标为（25.5，100.5），具体位置见图1中TWSA 和不同气候指标在月尺度上的变化，这些指标均呈现出以年为周期的规律性变化，除风速外，在冬季出现最低值，在夏季达到最大值。然而，当采用13点滑动平均一定程度上降低不同指标的季节变化后，仅降雨和比湿与TWSA 呈现明显的相关性，这也进一步验证了TWSA月动态变化与大部分气候指标的强相关性主要与周期性的季节波动有关。当进一步采用PLSR 模型分析每个指标对TWSA月动态变化的贡献，发现相对于降雨和温度，比湿在大部分网格中均为贡献度最高的指标，表明比湿的月动态变化与TWSA的月变化最为接近。尽管研究区降雨在季风气候下也呈现一定的季节变化，但其变化主要体现在丰水期，而枯水期TWSA主要受蒸散发的影响，降雨量的变化远小于TWSA，这可能是导致降雨对TWSA 月动态变化贡献度较低的原因。在年尺度下，70%以上网格中降雨为解释TWSA 动态变化的第一贡献因子，降雨仍为大部分区域TWSA长期变化的主控因素。

图1 研究区域图Fig.1 Map of the study area

图2 TWSA变化趋势、季节项以及趋势项方差贡献占比分布Fig.2 The TWSA changing trend and variance contribution of seasonal term and trend term to TWSA variance

图3 月尺度和年尺度下TWSA与气候指标之间的最大互相关系数分布图Fig.3 Maximum correlation coefficients between TWSA and climate indicators in monthly and annual scales

图4 不同尺度下每个网格PLSR模型的决定系数值和不同尺度下每个网格对TWSA变化贡献前三的气候指标分布Fig.4 Determination coefficient of each grid simulated by PLSR in different time scales， the top three climate factors of each grid that contribute to the TWSA variation in different time scales

图5 2002-2017年不同网格TWSA与8种气候要素变化图Fig.5 Variations of TWSA and eight climate factors in two different grids from 2002 to 2017

PLSR 模型的结果显示在月尺度上除四川省西北部部分网格外，研究区中部（贵州省、重庆市、湖北省和湖南省西部）模拟结果较差，R2值低于0.45，表明这些区域的TWSA 的变化可能受到人类活动的强烈影响。研究区中分布有长江与珠江两大水系，地表水资源丰富且地表地形坡度大，该区域修建了大量的水电站，研究区中部分布有三峡和龙滩两个大型水库前期的研究也表明三峡、龙滩等特大型水库蓄水期间会导致该地区重力异常增加［23］，蓄水量变化相当于该局部区域的80%以上的TWS变化［24］，这也反映了水库蓄水会对TWSA 的变化造成显著影响。因此，该区域水库可能为研究区中部影响TWSA 月尺度变化的主要人为因素。相对而言，研究区其他区域PLSR 模型均展示了较好的模拟结果，大部分网格R2大于0.60，最大值可达0.92，表明这些区域TWSA 的变化与气候波动基本一致，受人类活动影响较小。

为进一步分析TWSA 动态过程的影响因素，分别选取了贵州省和云南省某网格的TWSA、气候指标的变化进一步分析。这两个网格PLSR 模型的R2分别为0.21 和0.88，代表了是否主要受气候波动控制下的TWSA 的变化。从图5 中可以看出，当TWSA 主要受气候波动影响时，TWSA 与大部分气候指标的变化相似，表现出明显的周期性且波动幅度较大。但是，当TWSA受到人类活动时，TWSA的周期性变化明显减弱，除部分丰水年或枯水年出现明显的峰值或谷值外，其他年份的波动都非常小，这可能主要与水库消峰填谷的调控方式有关，一定程度上明显削减了丰枯季陆地水储量的季节波动幅度，相对于未受人为活动影响的地区表现出较弱的波动。这也说明该地区水库的修建确实对于当地防洪抗旱有一定的作用。

当从年尺度分析气候指标对TWSA 动态变化的预测能力时，结果显示研究区中部绝大部分网格模拟结果均出现大幅提升，人为活动对年尺度TWSA动态变化的影响明显减弱，仅三峡水库集水区域的网格R2值仍相对较低，表明研究区大部分地区年尺度TWSA 变化仍主要受气候变化的影响。一般情况下，水库主要在丰水期拦截洪水，枯水期向下游放水，因此水库蓄水库容一般在一个水文年内主要为周期性波动，而在不同水文年内蓄水库容的变化也取决于上游来水量的多少，即与气候变化也存在一定的联系。尽管水库的连续性蓄水会导致陆地水储量的持续增加，但这种蓄水情况一般仅发生在水库运行初期，很难长时间影响陆地水储量的年尺度变化。模拟结果也显示三峡水库集水区域内的网格年尺度TWSA变化仍一定程度上受到人为活动的影响，这可能与研究时段三峡水库多次频繁的蓄水活动有关。据百度百科资料记载，三峡水库从2003 年6 月首次蓄水开始，2006 年9 月开始实行了第二次蓄水，随后多年间多次启动175 m 试验性蓄水，研究时段内这种多次人为蓄水活动可能会对TWSA年尺度上的变化造成影响。总体上，PLSR模型在不同时间尺度上的模拟结果表明，研究区中部地区除三峡水库附近区域外，研究时间内水库仅对TWSA 月尺度的变化具有明显的影响，而对于TWSA年尺度的变化影响相对较弱，气候条件的变化仍是TWSA年尺度变化的主控因素。

分析过程中，PLSR 模型被选取用于分析气候指标对TWSA的预测能力，作为一种线型回归模型，当气候指标与TWSA之间存在较强的非线性关系时，该模型很难准确评估气候指标对TWSA 贡献。我国南方地区广泛分布大量的岩溶地貌，岩溶地区由于特殊的水文地质结构［25］，相对于非岩溶区水文过程更为复杂，气候指标与TWSA之间可能存在强烈的非线性关系，采用PLSR 线型模型可能会一定程度上低估气候指标对TWSA 动态的贡献，高估人为因素对TWSA 的影响。除此之外，选取的CSR Mascon 提供的原始分辨率为1°×1°，但由于GRACE 的信号衰减，Mascon的分辨率仍然受到GRACE 的频带限制［26］，其实际的分辨率可能远小于该精度，因此最后得出的结论仅用于认识区域上的空间变化规律认识，而不适用于解释具体一个或几个网格TWSA的影响因素。

4 结 论

基于8种不同的气候指标和GRACE 重力卫星数据，从网格尺度（1°×1°）上调查了我国南方7 省市区域上TWSA 动态变化的影响因素，相对于前期研究主要以大型流域为研究单元，以网格为基本研究单元有助于认识研究区内局部区域TWSA动态变化的影响因素。采用趋势分析、主导模式分析、互相关分析和PLSR 模型详细分析了不同气候指标与TWSA 之间的相关性以及气候指标对TWSA变化的解释程度，主要得出以下结论。

（1）2002 年至2017 年期间，研究区内除云南省和四川省东部区域外，75%以上地区TWSA 均呈上升趋势，在重庆、贵州和湖南省交界处上升幅度最大，平均增幅可达每年0.9 cm。云南、四川中南部、广西南部以及湖南湖北二省的东部区域中，TWSA随时间变化的主导模式为季节周期性变化；重庆、贵州、湖北和湖南四省市交界处，由长期趋势引导的年际变化占据了TWSA变化的主导地位。

（2）除研究区中北部区域外，TWSA 月动态变化与大多数气候指标呈良好的相关性，这主要由于研究区大部分地区受季风气候影响，TWSA与气候指标均呈明显的季节变化有关，与比湿的动态变化相似度最高。而年尺度下，TWSA 的变化仅与降雨和比湿相关性较高，降雨为70%以上地区解释TWSA 动态变化贡献度最高的指标。

（3）在月时间尺度上，研究区中北部地区TWSA的动态变化明显受到水库消峰填谷调控的影响，整体上TWSA 的季节波动远低于其他地区，15%网格的气候指标仅能解释陆地水储量17%～28%的变化，而其他区域TWSA 动态变化仍主要受气候条件的控制。进一步分析TWSA 年时间尺度变化时，研究区大部分地区TWSA 的变化仍主要受气候波动的影响，水库的影响相对较小，降雨为控制TWSA 年际变化的主控因素，陆地水储量60%以上的变化可由气候波动解释，仅三峡水库附近的TWSA年尺度的变化仍明显受到人类活动的强烈影响。