基于GRACE 卫星数据的中国陆地水储量变化和影响因素分析

2023-10-08石振君朱秀芳唐谊娟

干旱区地理(汉文版) 2023年9期

石振君，朱秀芳，唐谊娟

（1.北京师范大学遥感科学国家重点实验室，北京 100875；2.宁夏回族自治区遥感调查院，宁夏银川 750021）

陆地水储量通常包括地下水、土壤水、地表水（湖泊、河流和水库）、冰川与积雪等部分［1］。陆地水储量的变化通常被认为是一个关键的水文指标，它反映了整个地区水文输入（如降雨、降雪、冰川融化和径流）与输出（如蒸散、人类用水和径流）之间的净效应，影响下渗与基流过程。掌握陆地水储量变化的时空特征，了解其变化的原因有助于水资源的可持续性综合管理。

陆地水储量变化的增减会改变水文循环，威胁生态系统与人口的可持续用水［1］。美国国家航空航天局（NASA）和德国航空航天中心（DLR）于2002年3 月联合开发了重力恢复和气候实验（GRACE）卫星，利用GRACE 卫星观察到的地球重力场变化，可反演地球系统中的陆地水储量变化，这大大方便了对陆地水储量变化的研究［2］。迄今为止，不少学者针对陆地水储量变化的时空特征和影响因素开展了研究。相关研究指出中国陆地地区北方陆地水储量变化以减少为主，南方以增加为主［3-4］，云南省、天山山区以及西北地区陆地水储量变化整体上呈现下降趋势，且天山山区以及西北地区的陆地水储量变化呈现的明显季节变化主要由降水造成［5-7］。Yang等［8］、Wei等［9］分别对塔里木河和柴达木盆地陆地水储量进行了影响因素分析，均指出降水变化是所研究区域陆地水储量变化的主要驱动因素。李晓英等［10］从趋势性、相关性及时空变化特征3 个方面分析了长江流域陆地水储量变化以及与植被覆盖变化的关系，指出长江流域归一化植被指数与陆地水储量变化之间有强相关性。Meng等［11］从地表水平衡的角度出发，研究发现青藏高原陆地水储量的变化主要归因于降水和蒸散的变化。Zhong等［12］在海河流域的研究表明每年由气候驱动的陆地水储量虽呈现一定的波动上升趋势，但长期的人为取水导致陆地水储量持续下降。

以上研究大多在中国的局部区域开展，还没有看到对整个中国地区的陆地水储量变化的影响因素分析。为此，本文使用GRACE 卫星反演得到的陆地水储量变化数据，以中国陆地为研究区域，探究近十几年来陆地水储量变化，综合分析气象因素、下垫面和社会因素对陆地水储量变化的影响，为我国水资源管理提供适用的评估工具和科学决策支撑信息。

1 数据与方法

1.1 研究区概况

中国位于亚洲的东部、太平洋的西岸，地势呈西高东低的阶梯状分布。中国经纬度跨度大，受到气候、地势等因素的共同影响，中国的水资源分布呈现出东多西少、南多北少的格局。中国的东部和南部地区主要是亚热带季风和温带季风气候，降水呈现出明显的季节性特征，夏季和秋季降水多，春季和冬季降水少；西部和北部地区主要为温带大陆性气候和高原山地气候，全年降水均比较少。

1.2 数据来源

本文所用到的数据包括陆地水储量变化数据、气象数据、地表数据和社会经济数据。

陆地水储量变化数据［13］来自国家青藏高原科学数据中心，为2003—2019年的中国区域陆地水储量变化月数据，空间分辨率为0.25°。该数据集连接了CSR RL06 Mascon 和JPL RL06 Mascon 数据在GRACE 和GRACE-FO 之间的11 个月的间断期，但并未对GRACE 和GRACE-FO 各自数据之间存在的个别月份的缺失数据进行弥补。在使用数据前，使用线性插值方法［14］对个别月份的缺失数据进行了插补，得到了完整的陆地水储量变化时序数据。

气象数据包括降水量、温度和标准化降水蒸散发指数（SPEI）。降水量和温度数据均为来自国家气象数据中心的2003—2019年的月数据，空间分辨率为0.25°。SPEI 数据为来自全球标准化降水蒸散指数数据库的2003—2019年的月数据，空间分辨率为0.5°。

地表数据包括土地利用现状图、归一化植被指数（NDVI）和数字高程模型（DEM）数据。土地利用现状图为来自于资源环境科学与数据中心的1 km分辨率的数据，时间分别为2005、2010 年和2015年。NDVI 数据为NASA 的陆地过程分布式数据档案中心提供的2003—2019 年的MODIS 的陆地标准产品，空间分辨率为1 km，时间分辨率为8 d。NDVI数据包括年度NDVI 数据集和月度NDVI 数据集。月度NDVI数据集通过最大值合成法生成，年度NDVI 数据集是在月度数据基础上采用平均值合成法生成。DEM 为基于最新的SRTM V4.1 数据经重采样生成的分辨率为1 km 的数据。社会经济数据包括国内生产总值（GDP）数据和人口数据，2 类数据均来自与资源环境科学与数据中心，分辨率均为1 km，所用数据的时间分别为2005、2010 年和2015年。在使用前，对上述数据进行了重投影和重采样处理，重采样后的分辨率为0.25°。

1.3 研究方法

（1）曼-肯德尔法（Mann-Kendall，M-K）趋势检验

M-K 趋势检验是一种非参数方法［15］，其目的是评估变量随着时间的变化是否有单调上升或下降的趋势［16-17］。M-K 检验不需要样本遵循特定的分布，也不会受少数异常值的干扰，即使数据有缺失值或者存在值低于一个或者多个检测的限制，也可以使用M-K检验。

（2）经验正交函数（EOF）分析方法

EOF 分析方法是提取数据主要特征量的一种方法［18-19］。EOF分解得到的特征向量对应的是空间样本，也称作空间特征向量或者空间模态，反映地理要素的空间分布特点；主成分对应的是时间变化，也称时间系数，反映相应空间模态随时间的变化而发生的变化。

（3）皮尔逊（Pearson）相关系数

Pearson相关系数是一种线性相关系数，也是最常用的一种相关系数［20］。其值介于-1到1之间，绝对值越大表明相关性越强。值大于0 时，表明2 个变量正相关；值小于0 时，表明2 个变量负相关；当等于0时，表明2个变量不是线性相关的，但是可能存在其他方式的相关。

（4）随机森林

随机森林是基于Bagging 框架通过集成学习的思想将多棵树进行集成的一种算法，其基本单元是决策树，本质是集成学习方法。随机森林模型不易发生过拟合问题，具有较强的泛化能力。它可以用于评估各个特征对于结果的重要性［21］。

（5）地理探测分析

地理探测器是探测因子的空间分异性，并揭示其驱动力的一组方法［22］。地理探测器的核心思想是如果某个自变量对某个因变量有重要的影响，那么该自变量和因变量的空间分布应该具有相似性。地理探测器既可以用于数值型数据，也可以用于定性数据。地理探测器通过计算和比较各单因子的q值以及两因子叠加后的q值，判断两因子是否存在交互作用，以及交互作用的强弱、方向、线性还是非线性等。地理探测器包括分异及因子探测、交互作用探测、风险区探测和生态探测4个子探测器。

1.4 技术路线

本文的技术路线（图1）主要包括3个内容：基于M-K 趋势检验的陆地水储量变化的趋势分析，基于EOF的陆地水储量变化的时空特征分析，基于Pearson 相关、随机森林和地理探测器的陆地水储量变化影响因素分析。其中，陆地水储量变化的影响因素分析的主要步骤如下：

图1 技术路线Fig.1 Technical route

（1）影响因子选择：包括气象因素（气温、降水、SPEI）、下垫面因素（不透水层占比、水体占比、NDVI、高程、坡度）和社会经济因素（GDP、人口）3个方面10个影响因子。

（2）Pearson相关分析：在10个影响因子中，NDVI、SPEI、降水量和温度会随着时间的变化而发生较大变化，且存在时间序列的数据，其他6个因子在短时间内的变化量很小，且没有时间序列的数据。因此，基于Pearson 的相关分析分成了2 部分，一部分是时间维的分析，另一部分是空间维的分析。时间维的分析是考虑NDVI、SPEI、降水量、温度与陆地水储量变化存在滞后效应，逐栅格计算NDVI、SPEI、降水量、温度与陆地水储量变化在0～3个月滞后时间内的相关系数，并找出最大滞相关系数对应的滞后月数。空间维的分析是首先逐栅格计算时间序列数据（包括NDVI、SPEI、降水量、温度和陆地水储量变化）的多年均值，然后以所有格网为样本，分别计算陆地水储量变化和10个影响因子的Pearson 相关系数，以相关系数大小来反映影响因子和陆地水储量变化的相关程度。

（3）基于随机森林进行影响因子重要性排序：对随机森林中某个树节点的样本方差与分裂后2个叶子节点样本方差进行差值运算，用该差值来评价某个特征的重要性程度，差值越高表示该特征重要性越大。计算得到模型所有特征重要性后，对所有特征作归一化处理，对于某一模型，所有特征重要性值累计为1，特征重要性具体的值本身并无意义，因此分析时关注的是特征重要性之间的相对差异。

（4）基于地理探测器探究单个影响因子及因子间的交互作用对陆地水储量变化空间差异的解释力：在使用地理探测器计算之前，先对自变量进行离散化处理，按照干旱等级划分标准对平均SPEI进行分级，对其他9个影响因子（年均NDVI、年均降水量、年均温度、人口、GDP、不透水层占比、水体占比，高程和坡度）用自然间断点分级法重分类为6 类。然后以离散化后的因子为x，以地表水储量变化趋势值为y，进行因子探测、交互探测和生态探测分析。因子探测用来确定10 个影响因子多大程度上解释了陆地水储量变化的空间分异；交互探测用于评估2个自变量共同作用时是否会增加或减弱对因变量的解释力；生态探测用于比较2 个影响因子之间对陆地水储量变化的空间分布的影响是否有显著的差异。

2 结果与分析

2.1 陆地水储量变化趋势

由2003—2019 年中国陆地水储量变化趋势图（图2）可以看出，大部分地区陆地水储量变化趋势明显，其中陆地水储量变化显著增加的区域主要分布在松花江、嫩江和松嫩平原附近以及柴达木盆地-长江-东南沿海条带上，陆地水储量变化显著减少的区域主要分布在中国西南以及新疆-黄土高原-华北平原条带上。

图2 中国陆地水储量变化趋势Fig.2 Change trend of terrestrial water storage in China

2.2 陆地水储量变化时空特征

对2003—2019年的204幅（每月1幅影像）陆地水储量变化影像进行EOF分析。EOF的前3个模态的解释方差分别为47.98%、19.14%、16.89%，累计方差贡献约为84.01%。由图3a～c可以看出，第1模态表现出中国的陆地水储量变化随纬度带交替反相变化的特征，从高纬到低纬呈现出水储量变化“高-低-高-低”的变化特征；第2模态表现出中国南北方陆地水储量变化反相变化的特征；第3 模态显示出四川盆地、横断山脉、云贵高原、柴达木盆地陆地水储量变化比较高，而新疆部分区域和中国东南部水储量变化明显偏低。结合时间系数看（图3d～f），第1模态的时间系数在研究时间段的前9 a均为负值，而后6 a均为正值，数值呈现逐渐增大的特征，说明随着时间的变化，中国陆地的水储量变化由第1 模态反应的由北到南的“高-低-高-低”特征的反向逐渐演变到由北到南的“高-低-高-低”的空间格局特征。第2模态和第3模态反应出的特征具有震荡变化的特点，总体来说每年的上半年的时间系数为负值，下半年的时间系数为正值，说明上半年的各月的陆地水储量与第3 模态反应的特征相反，而下半年与第3模态反应的特征一致。

图3 中国陆地水储量变化的EOF分析Fig.3 EOF analysis of terrestrial water storage change in China

2.3 影响因素分析

2.3.1 基于地理探测器的影响因子重要性分析因子探测、交互探测和生态探测的结果如表1～3 所示。在单个因子中，降水量对陆地水储量变化解释力最强，q值为0.276；水体占比对陆地水储量变化解释力最弱，q值仅为0.001。在因子交互作用中，对陆地水储量变化解释力最强的是降水量和温度的交互作用，q值达到0.461；其次是降水量和DEM的交互作用，q值为0.377；水体占比和坡度的交互作用对陆地水储量变化的解释力最弱，q值为0.037。综合比较因子探测和交互探测的结果，可得到因子两两交互作用相较于单个因子对于陆地水储量变化的解释力均为增加。生态探测结果表明，降水量与DEM、GDP、NDVI、人口，SPEI 与GDP、人口，温度与DEM、GDP、NDVI、人口、SPEI 对陆地水储量变化空间分布的影响有显著的差异，其余2个因子之间对陆地水储量变化空间分布的影响无显著差异。

表1 因子探测结果Tab.1 Factor detection results

表2 因子交互探测结果Tab.2 Factor interactive detection results

表3 生态探测结果Tab.3 Ecological detection results

2.3.2 基于随机森林的影响因子重要性分析由基于随机森林的影响因子重要性分析结果（表4）可知，10个影响因子对陆地水储量变化的重要性排序为降水量＞温度＞SPEI＞人口＞GDP＞DEM＞NDVI＞坡度＞不透水层占比＞水体占比，降水量因子对陆地水储量变化的影响最大，特征重要性值为0.335，水体占比因子对陆地水储量变化的影响最小，特征重要性值仅为0.032。在因子类型层面分析，气象因子对于陆地水储量变化的影响最大，其次是社会经济因子，下垫面因子对于陆地水储量变化的影响最小。

表4 基于随机森林的影响因子重要性评价结果Tab.4 Evaluation results of importance of influencing factors based on random forests

2.3.3 基于Pearson相关分析的影响因子重要性分析由空间维上陆地水储量变化和各影响因子相关性分析的结果（表5）可知，从全区来看，DEM、人口和水体占比因子与陆地水储量变化不存在线性相关关系；不透水层、坡度、GDP与陆地水储量变化呈负相关，且不透水层与陆地水储量变化的负相关性最强，相关系数为-0.16；NDVI、降水量、SPEI和温度与陆地水储量变化呈正相关关系，且降水量与陆地水储量变化的正相关性最强，相关系数为0.339。

表5 Pearson相关分析结果Tab.5 Pearson correlation analysis results

陆地水储量变化对气候因素和NDVI的响应滞后时间分布如图4所示。大部分区域陆地水储量变化与NDVI在当月以及滞后1个月时相关性最强，陆地水储量变化与NDVI 存在2～3 个月滞后的区域主要位于中国中部地区（图4a）。中国大部分地区的陆地水储量变化与降水量存在1 个月的滞后（图4b）。具体来说，中国西部靠北地区陆地水储量变化与当月降水量相关性最高，中国东北部、东南部以及西部部分地区陆地水储量变化在滞后1个月时与降水量相关性最高，中国中部地区陆地水储量变化在滞后2 个月时与降水量相关性最高，华北平原中部、河北以及河北北部地区陆地水储量变化在在滞后2个月时与降水量相关性最高。中国大部分地区的陆地水储量变化与温度存在0或3个月的滞后（图4c）。具体来说，中国东北部、东南部以及西北部地区陆地水储量变化与当月的温度相关性最高，青藏高原中部和南部陆地水储量变化在滞后2个月时与温度相关性最高，云南省、华北平原以及华北平原西部地区陆地水储量变化在滞后3个月时与温度相关性最高。中国大部分地区的陆地水储量变化与SPEI 存在1 或3 个月的滞后（图4d）。具体来说，中国东南部、东北部以及中部靠北地区陆地水储量变化在滞后1个月时与SPEI相关性最高，云南省西北部以及青藏高原东南部陆地水储量变化在滞后3个月时与SPEI相关性最高。

图4 陆地水储量变化与影响因子之间的滞相关分析结果Fig.4 Lag correlation analysis results between terrestrial water storage change and influencing factors

由表6可知，在中国陆地水储量变化与NDVI存在0个月滞后的栅格数目最多，占比达到42.0%；与降水量存在1 个月滞后的栅格数目最多，占比达到42.1%；与温度存在0 个月滞后的栅格数目最多，占比达到43.7%；与SPEI 存在1 个月滞后的栅格数目最多，占比达到30.8%。

表6 不同滞后时长下的栅格数目占比Tab.6 Proportion of grid number under different lag length /%

2.3.4 3种重要性分析结果比较对比地理探测器、随机森林和Pearson 相关分析3 种方法对各个影响因子对陆地水储量变化影响程度大小的排序（表1、表4 和表5）可知，地理探测器和随机森林得到的结果总体相差不大，但和Pearson相关分析结果相比相差较大，原因在于Pearson相关分析只能反映2个变量间线性相关程度的强弱，而随机森林和地理探测器并不局限于线性关系。结合3种方法得到的结果可知，气象因素对陆地水储量变化的影响最大，下垫面因素对陆地水储量变化的影响最小，在气象因素中降水量对陆地水储量变化最为重要。

3 讨论

本文从气象、社会经济和下垫面3个方面入手，定量分析了陆地水储量变化的影响因素，研究发现降水量是导致陆地水储量变化的主要原因，这与Wei等［9］、Meng等［11］研究结论一致。但需注意的是，降水量增加并不能完全代表陆地水储量增加，例如新疆部分地区降水量增加明显，但蒸散发能力明显高于降水量［23］，且新疆的农业绝大部分为灌溉农业，农业灌溉在各种人类用水中占比最高，因此陆地水储量整体仍然呈下降趋势；类似的，华北地区虽降水量变化不明显，因工业、农业大量开采地下水，导致陆地水储量不断减少［24］。此外，经分析得到黄土高原地区陆地水储量变化呈显著减少趋势，但相关研究表明黄土高原近年来植被增加显著，主要原因是退耕还林还草工程实施以来黄土高原深层土壤储水量严重亏损［25］。

本文结论是以中国陆地地区作为整体进行影响因素分析所得到的，对于各个地区，因气候不同、经济社会发展程度不同、地理因素不同，各个影响因子以及因子的交互作用对各地的影响不同，后续研究将进一步分析提取中国陆地水储量变化热点区，以热点区为单位探究同一因素对不同地区的作用效果，深层次研究各热点区引起陆地水储量变化的主导因素。在选择影响因素时，因可获取数据有限，未考虑生态输水产生的影响［26］以及工业用水、农业用水等［27］的影响，如何科学补充影响因子，使之更具完整性，仍需进一步考虑。