论文快递｜刘家宏等：黄河流域地下水储量时空演变特征及影响因素分析

地下水是重要的淡水资源，对干旱和半干旱区工农业生产、生活和生态系统的维持具有重要的支撑［1］作用。 过去几十年，受气候变化和人类活动影响，地下水灌溉区域地下水位下降明显［2］。 气候变化对地下水系统具有直接和间接影响，直接影响表现为改变地下水补给量或消耗量，间接影响为改变地下水利用方式。 这些影响同时受土地利用等因素的显著调控［3］。鉴于缺乏足够的地下水位监测井和对复杂水文地质信息的全面了解，量化长时间尺度地下水变化及分析气候和人类活动对其造成的影响具有挑战性。

重力恢复及气候实验（GRACE）卫星于2002 年4月成功发射以来，有效监测地球重力场变化。 其继任卫星GRACE Follow－On 于2018 年发射，为地表质量变化持续监测提供支撑。 依据GRACE／GRACE Follow－On监测数据通过对地表质量变化的反演进而得到陆地水储量变化，是全球和区域尺度陆地水储量变化监测的重要手段［4－6］。 陆地水储量包含地表水、土壤水、地下水、冰川积雪等，基于水量平衡原理去除其他分量后，可得到地下水储量的变化情况，该模式在地下水研究中得到广泛应用［7］。 Feng 等［8］根据2003—2010 年GRACE 卫星监测数据和地面测量数据，基于水量平衡法估算了华北地区地下水储量变化，估算结果与同一时期监测井站估计结果一致。 Zhong 等［9］利用水量平衡法估算了2005—2011 年辽河西部地下水储量的变化，评估了该地区地下水枯竭状况。 Xu 等［10］利用GRACE 和全球陆面数据同化系统（Global Land Data Assimilation System，GLDAS）数据评估了华北平原地下水长期变化趋势，并分析了降水和南水北调工程的影响。

黄河流域位于干旱和半干旱气候区，水资源供需矛盾十分突出。 流域地表水资源开发利用率已达80%，超过流域水资源承载能力［11］。 地下水资源在流域水资源安全和生态保护方面具有至关重要的作用。气候变化背景下黄河流域气温显著升高，降水量减小。作为生态系统敏感区域，加之退耕还林还草工程的实施，黄河流域植被覆盖面积显著增加［12］。 许多学者对黄河流域地下水变化进行了相关研究。 Lin 等［13］依据水量平衡原理，利用时间稳定性分析识别了2003—2016 年黄河流域地下水波动的热点和代表性监测区域。 王希栋等［14］采用GRACE／GRACE Follow－On 和GLDAS 结合的方式评估了2003—2021 年黄河流域地下水储量变化，定性分析了降水和植被对流域上游、中游和下游水储量变化的贡献。 然而，现有研究对地下水时空变化特征讨论较多，对地下水储量空间变化主要驱动因素及驱动过程的认识存在不足。

考虑到GRACE 和GRACE Follow－On 卫星间有近1 a 时间缺测，常用的线性和样条插值在捕捉时空变化特征方面存在不足，本研究基于时间序列分解与随机森林算法，构建像元尺度流域地下水储量变化耦合模型，分析近20 a 黄河流域地下水储量时空演变规律，量化影响因素的贡献，揭示气候和人类活动对地下水变化的驱动过程，以期为黄河流域地下水资源利用和相关政策制定提供参考。

黄河流域覆盖青藏高原、黄土高原、内蒙古高原和华北平原，地势复杂，贯穿东、中、西三大地势阶梯，流域面积79.5 万km2。 流域多年平均气温为5.8 ℃、降水量为466.5 mm、潜在蒸散发量为653.6 mm［15］。 流域地下水分为松散岩类平原盆地孔隙水、松散岩类黄土高原孔隙水、碳酸盐岩裂隙岩溶水、碎屑岩孔隙裂隙水、岩浆岩与变质岩类裂隙水、冻结层水，对流域水资源循环演化有重要调节作用。 流域地下水资源划分为4 个二级区、11 个三级区、31 个四级区、109 个五级区［16］。 本文以黄河流域花园口水文站以上区域为研究对象（后文统称黄河流域），分析地下水储量变化，根据水沙特性和地形地质，将黄河流域分为7 个子流域，分别为龙羊峡以上、龙羊峡—兰州、兰州—头道拐、头道拐—龙门、龙门—三门峡、三门峡—花园口、内流区。

1.2  数据来源

1）GRACE 和GRACE Follow－On 监测数据。 采用得克萨斯州大学空间研究中心（University of Texas Center for Space Research，CSR）、美国喷气推进实验室（Jet Propulsion Laboratory，JPL）发布的Mascon RL06 V2 等效水高数据产品，时间序列为2003 年1 月到2022 年12 月，空间分辨率分别为0.25°×0.25°和0.5°×0.5°。 不同数据产品对流域水储量变化研究各有优势，相关研究表明，融合多个数据产品可有效降低水储量变化监测的不确定性［17］。 本研究以CSR 和JPL 数据的算术平均值作为地下水储量分析的基础，考虑到CSR 和JPL 数据在空间分辨率上的差异，为了更好地表征空间变异性，对JPL 产品重采样至分辨率为0.25°后，再与CSR 数据进行算术平均。 对于传感器故障导致的GRACE 数据序列缺失，本研究采用三次样条法逐像元进行插值。

2） GLDAS 数据。 由美国航空航天局和国家海洋与大气管理局合作开发，采用陆地表面模型和数据同化技术，结合卫星和地面观测数据，生成最佳陆面通量和状态量［18］。 以GLDAS Noah－V2.1 模型输出的2 m 深度土壤含水量为基础，评估黄河流域的土壤水储量变化。该数据的空间分辨率为0.25°，时间分辨率为月。

3）气象数据。 采用国家青藏高原科学数据中心发布的中国1 km 分辨率逐月降水量、逐月平均气温和潜在蒸散发量表征流域气候特征。 上述数据集的气温和降水来源于CRU 发布的全球0.5°分辨率气候数据集和WorldClim 发布的全球高分辨率气候数据集，通过Delta 空间降尺度方案在中国地区进行降尺度，并采用中国496 个独立气象观测点的数据对降尺度结果进行验证。 潜在蒸散发量以最低气温、最高气温和月均气温为输入，基于Hargreaves 公式计算得到。

4）MODIS 数据。 采用MOD13C1 归一化植被指数（NDVI）和MOD15A2H 叶面积指数（LAI）表征地表植被覆盖和植被生长情况。 数据的时间分辨率为8 d，空间分辨率为500 m。 以最大值法合成月尺度和年际尺度NDVI 和LAI，并使用S－G 滤波算法对NDVI 和LAI进行时间滤波重建，以减少噪声并插补缺失数据。

5）GLEAM 数据。 GLEAM（Global Land Evaporation Amsterdam Model）是一个评估陆地蒸散发的模型，利用Priestley－Taylor 公式和微波观测数据来分析潜在和实际蒸散发。 该模型结合了地表净辐射、近地表气温、植被光学特征和土壤水分，基于Gash 模型计算降水拦截损失，从而评估水体及冰雪覆盖区的蒸发量［19］。 本研究采用GLEAM 3.7 a数据集，该数据以遥感和再分析数据为驱动，具有月尺度时间分辨率和0.25°×0.25°的空间分辨率。

陆地水储量包含土壤水、地下水、雪水当量、冰川积雪和冠层含水等。 冰川积雪及冠层含水对黄河流域年际水储量变化影响较小，本研究不予考虑。 以GRACE 卫星陆地水储量异常与GLDAS 模拟土壤水储量异常的差异作为区域地下水储量异常。 由于GLDAS 模拟输出土壤水为当前土壤水的状态值，因此为了与水储量异常保持一致，以研究时段土壤水储量相对于2004—2009 年多年平均土壤水储量的差值作为土壤水储量的异常。

GRACE 和GRACE Follow－On 数据序列间存在11个月（2017 年7 月—2018 年5 月）的数据缺失，缺失数据序列较长，传统的插值方法在人类活动影响显著的干旱半干旱区域插值精度不足［20］。 随机森林算法是常用的高效机器学习算法，在地球科学、遥感领域得到广泛应用，尤其在全球和区域尺度应用效果较好［21］。同时，随机森林算法的树结构模型通过量化变量对树节点的贡献评估变量重要性，为分析影响因素对地下水储量变化的贡献提供了新的视角。

为了更好模拟区域地下水储量时空变化特征，采用时间序列分解（STL）和随机森林算法结合的方法构建气象水文和人类活动与地下水储量变化模型。 首先利用STL 方法［22］将地下水储量变化分解为趋势项、季节项和残差项，然后在不同时间尺度上构建随机森林模型，最后将不同时间尺度的拟合项相加得到最终地下水储量变化模拟值。 随机森林模型是一种集成模型，由多个决策树组成，每个决策树都基于机器学习对数据进行分类或回归，该模型是在0.25°×0.25°网格尺度下学习，输入变量为降水、气温、潜在蒸散发量、NDVI、LAI、陆面蒸散发量，因变量（输出结果）为地下水储量异常，模型表达式如下：

式中：GWSA（i，j）为（i，j）网格的地下水储量异常，P（i，j）、Tem（i，j）、PET（i，j）、NDVI（i，j），LAI（i，j）、ET（i，j）分别为（i，j）网格的降水、气温、潜在蒸散发量、归一化植被指数、叶面积指数、陆地蒸散发量。

SHAP（SHapley Additive exPlanations）模型是一种解释机器学习模型输出的方法，该模型借鉴博弈论中的Shapley 值概念，可量化每个特征对模型预测结果的贡献，关于该模型的详细介绍见参考文献［23］。 为了提高模型的可解释性，采用SHAP 方法评估不同时间尺度随机森林模型各自变量对因变量的贡献，并将不同时间尺度贡献之和作为变量最终贡献。 考虑到各因素对地下水储量变化影响存在明显的交互作用［24］，模型构建时对变量进行了归一化处理，并假定贡献率大于0.4 的因素为地下水储量变化的主导因素。

图1 为2003—2022 年黄河流域中上游地区地下水储量异常值年际变化情况，近20 a 地下水储量持续显著减小，变率为－3.3 mm／a。 各子流域地下水储量均呈显著减小趋势，变率为－0.9 ～－5.7 mm／a，其中三门峡—花园口地下水储量减小幅度最大，龙羊峡—兰州地下水储量减小幅度较小（见图2、表1）。 2008 年之后，兰州以下相邻子流域地下水储量差异明显增大，地下水储量空间变异性显著增强。 三门峡—花园口地下水储量变化在2004 年、2012 年、2018 年和2022 年均出现极大值，其他时期地下水储量以持续减小为主，2020 年地下水储量最小。 内流区和龙门—三门峡，龙羊峡以上和龙羊峡—兰州，地下水储量时间变化特征基本一致。

表１　各子流域地下水储量年际变率

 

图１　黄河流域中上游地区地下水储量异常值年际变化情况

 

图２　各子流域地下水储量异常值年际变化情况

 

黄河流域地下水储量年际变化空间格局如图3所示，从东到西地下水储量年际变化幅度表现为先减小后增大的格局，流域东部地下水减小速度最快。

图３　黄河流域地下水储量年际变化空间格局

 

具体而言，黄土高原中心区域地下水减小幅度较大，尤其在汾河盆地和关中平原，地下水储量减小速度为－15～－30 mm／a。 流域南部边缘（甘肃境内）和宁夏部分区域地下水储量呈增加趋势，增加速度为0 ～6 mm／a。 对年际变率进行显著性分析，结果表明趋势拟合决定系数（R2）大于0.5 的网格占比约为72%，主要分布在兰州以下区域、黄河源区中部及扎陵湖和鄂陵湖区域。

3.2  地下水储量影响因素分析

图4 为黄河流域降水年际时空变化情况（图中①为头道拐—龙门、②为内流区、③为龙门—三门峡、④为龙羊峡—兰州、⑤为兰州—头道拐、⑥为三门峡—花园口、⑦为龙羊峡以上，下同），2003—2022 年黄河流域年降水量无显著趋势，其中2015—2022 年波动较大。 相较于2015 年，2019 年流域平均年降水量增加了122 mm，到2022 年减小至接近2015 年值。

图４　黄河流域降水年际时空变化

 

相较于降水的变化，地下水储量仅在2017 年略微增大，其他时期均持续减小。 这一现象说明地下水对降水的响应具有明显的滞后性，并且降水增加对地下水储量增加的补给效应不显著。 黄河流域降水具有明显的空间变异性，其中龙门—三门峡、三门峡—花园口2 个子流域降水年际波动较为明显，2004 年、2011 年、2016 年和2021 年为降水显著增加年份，除2004 年外，地下水储量均在后一年表现为增加，说明地下水储量变化对降水变化的响应具有时间尺度效应。

气温通过蒸发影响地下水补给，潜在蒸散发表征区域最大蒸散发能力。 当气温升高时，潜在蒸散发量的增加意味着更多的水分可能从地表和土壤中蒸发，减少了能够渗透到地下的水量，从而影响地下水的自然补给。 黄河流域气温和潜在蒸散发量年际变化如图5 和图6 所示。

图５　气温年际时空变化

 

图６　潜在蒸散发量年际时空变化

 

研究时段气温具有明显的升高趋势，在2006 年、2016 年和2022 年流域平均和子流域气温达到时段极值。 潜在蒸散发量无显著变化趋势，波动明显，与气温一致性较好。 气温升高促进土壤和植被蒸发，因而相应时期的地下水储量表现为减小趋势，地下水储量对气温升高的响应速度大于对降水增加的。

黄河流域陆地蒸散发量、 归一化植被指数（NDVI）、叶面积指数（LAI）年际时空变化如图7～图9所示。 自1999 年退耕还林工程实施后，黄河流域植被绿度显著提高［25］。 本研究也表明2003—2022 年黄河流域平均NDVI 和LAI 均呈显著增大趋势，植被覆盖状况和植被生长状况持续向好。 黄河流域陆地蒸散发量也表现出明显的增加趋势。 空间区域上，龙门—花园口区域即黄土高原退耕还林核心区NDVI 和LAI 显著增大。 相应的龙门到花园口区域陆地蒸散发量也表现出明显的增加趋势。 陆地蒸散发量变化特征与前述气温和降水变化基本一致，2004 年、2008 年、2012 年和2022 年蒸散发量显著减小，对应年份降水量减小、气温下降，植被对地下水储量的影响同时受气候调节。

图７　陆地蒸散发量年际时空变化

 

图８　ＮＤＶＩ 年际时空变化

 

图９　ＬＡＩ 年际时空变化

 

3.3  影响因素贡献量化分析

图10 为各子流域降水、气温、NDVI、LAI、蒸散发量、潜在蒸散发量的SHAP 贡献值。 整体而言，各因素对黄河流域地下水储量变化的贡献从大到小依次为NDVI（0.80）、气温（0.55）、蒸散发量（0.51）、潜在蒸散发量（0.51）、LAI（0.50）和降水（0.43），表明黄河流域地下水储量减小主要由NDVI 增大、气温升高和蒸散发量增加等共同驱动。 黄河流域不同区域地下水储量变化主要驱动因素表现出显著的空间差异性。 龙羊峡以上地下水储量变化主要驱动因素为陆地蒸散发量，说明该子流域地下水储量变化主要依赖大气条件和植被的相互作用。 龙羊峡—兰州地下水储量变化主要受植被状况和覆盖（NDVI 和LAI）、气温、降水、蒸散发量和潜在蒸散发量共同影响，这些因素组合反映了植被状况、气候变化对水资源的综合影响。 兰州—头道拐地下水储量变化主导因素为NDVI，植被覆盖变化是该区域地下水变化的关键。 头道拐—三门峡地下水储量变化主要由NDVI、LAI 和气温等因素共同驱动，表明该子流域植被覆盖、植被生长状况和气候条件对地下水储量变化影响显著。 内流区地下水储量变化主要由NDVI、潜在蒸散发量、气温和陆地蒸散发量驱动，这些因素反映了复杂地表－大气相互作用、植被状况对地下水的影响。 三门峡—花园口NDVI 是地下水储量变化主要驱动因素，揭示了植被绿度提高对该区域地下水资源具有不可忽视的影响。

图１０　子流域影响因素ＳＨＡＰ 贡献值

 

3.4  讨论与建议

地下水储量变化受地下水补给量和消耗量共同影响，研究发现近20 a 黄河流域气候呈暖干化趋势，降水无明显增加，气温上升显著。 加之陆地表面蒸发量显著增加，地下水消耗量大于补给量，导致地下水储量显著减小。 虽然2015 年后黄河流域丰水年频次增加，但是地表植被覆盖面积增加、气温上升导致水汽压差增大、蒸散发量增大，加之地下水对降水响应的滞后性，地下水储量不能得到及时恢复，而存储的土壤水在下一个干旱年被消耗。 因此，黄河流域地下水储量减小是气候变暖和植被变绿共同作用的结果。 在流域不同区域植被对地下水储量变化的贡献均较为显著，其他因素贡献不同，与植被的作用方式也存在差异，应对黄河流域地下水储量减小应该采取与地域特点相一致的策略。

分析了近20 a 黄河流域地下水储量时空演变规律，量化了各因素的贡献，揭示了气候和人类活动对地下水变化的驱动过程，主要结论如下。

1）地下水储量时空变化特征：时间上，2003—2022 年黄河流域地下水储量显著减小，变率为－3.3 mm／a，其中头道拐—龙门和三门峡—花园口变率为－5.7 mm／a 和－5.5 mm／a。 空间上，地下水储量变幅表现为自东向西先减小后增大的格局，在汾河盆地和关中平原减小幅度较大（－15～－30 mm／a）。

2）因素变化特征：黄河流域气候呈暖干化趋势。2003—2022 年降水变化趋势不显著，降水量明显波动。 2015 年后丰水年频次增加，气温呈现出明显的上升趋势，且2015 年之后平均气温显著高于之前。 相比之下，潜在蒸散发量未表现出一致的变化特征。 地表植被覆盖和植被生长状况持续向好，地表蒸散发量显著增大。

3）变化驱动机制：黄河流域地下水储量的减小是植被增加、气温升高和蒸散发量变化共同驱动的结果。地下水储量空间变化主导因素差异显著，龙羊峡以上由陆地蒸散发过程驱动，龙羊峡—兰州由植被覆盖、气温、降水和蒸散发共同驱动，兰州—头道拐由植被覆盖变化驱动，头道拐—三门峡由植被覆盖变化和气温驱动，内流区由植被和气候变化驱动，三门峡—花园口由NDVI 驱动。

本研究未考虑地下水抽取及原煤开采引起的地下水排泄等人类活动对地下水储量变化的影响。 同时，黄河流域植被变绿与生态恢复工程密切相关，需构建合适的人类活动量化指标，综合考虑人类生产生活用水、提取地下水及煤炭开采和生态恢复工程，分离人类活动的正负效应。 另外，本研究在子流域划分时以地表水为依据，对水文地质条件考虑不足，后续研究应以黄河流域地下水资源分区为标准，考虑地质水文参数和土壤入渗率，从降水－土壤水－地下水三水转换关系及机理层面揭示地下水储量变化的机制。