基于GRACE 卫星的黄淮海地区地下水储量变化与归因分析

王 坤，李玉龙，李琛曦，赖燕玲，杨传国,3

(1.河海大学 水文水资源学院, 江苏 南京 210098； 2.河北省地质资源环境监测与保护重点实验室, 河北 石家庄050021； 3.河海大学 水灾害防御全国重点实验室, 江苏 南京 210098)

地下水是水资源的重要组成部分，由于具有分布广、水质好、可持续开发利用等许多优点，已成为最重要的和最佳的供水选择[1]，特别是在水资源相对匮乏的北方地区。黄淮海地区是中国重要的粮食生产基地[2]，是人均水资源量较少的地区之一，也是人类活动对水资源影响最严重的区域[3]。大规模开发地下水资源导致地下水位下降严重，产生了水质恶化、地面沉降、海水入侵等一系列水安全问题[4]，已成为制约地区发展的主要因素。全面监测地下水并掌握其变化规律是当务之急。传统人工监测地下水的方法耗费人力物力，且由于监测站点布设的困难，难以实现地下水的大面积整体监测。随着重力卫星探测技术的发展，通过测量地球重力场的变化反演全球陆地水储量变化，为实现地下水的区域整体监测提供了新途径[5]。国内外诸多学者利用GRACE （Gravity Recovery and Climate ）重力卫星在不同区域开展了地下水研究，证实了在大尺度范围内反演地下水资源量变化的可靠性[6]。Famiglietti 等[7]和Scanlon 等[8]利用GRACE 数据反演加利福尼亚州中央山谷地下水储量变化，与研究区实测地下水位验证结果大体一致。国内学者验证了GRACE 与GLDAS （Global Land Data Assimilation System）数据反演黄河流域[9]、海河流域[10]、三江源[11]等流域地下水储量变化具有较好适用性。

黄淮海地区地下水储量变化引起越来越多的研究者关注，对其影响的主导因素开展了相关研究。Guo 等[12]研究发现人类活动在海河流域地下水储量变化中占主导作用，农业灌溉是地下水储量消耗的主要因素，而降水和蒸发量相对稳定；林雪敏等[13]研究表明农业用水量对京津冀地区地下水储量变化的影响程度明显大于降水量与蒸散量；秦欢欢等[14]利用模型设定不同情景模拟华北平原水资源利用状况，结果表明农业节水措施和南水北调工程可以有效恢复地下水储量。南水北调通水后黄淮海地区供水结构发生了改变，相关研究较少，亟需采用多源数据深入分析区域地下水储量变化的特征与成因。因此，本文利用最新的GRACE 数据和GLDAS 陆面同化数据，定量计算地下水储量变化，评估黄淮海地区地下水储量变化的时空特征，结合降水、蒸散发、地下水供水、跨流域调水等自然和人为因素，探讨影响地下水储量变化的原因，从而为科学评估与合理利用地下水提供依据。

1 研究区概况

本文选择黄淮海地区，包括海河流域、淮河流域及黄河流域花园口以下地区为研究区域（图1）。地理位置位于北纬31°～43°，东经112°～122°，北起燕山南麓，南抵桐柏山、大别山，西起太行山，东临渤海、黄海，总面积为67.8 万km2。黄淮海地区属温带季风气候，季节变化明显，冬季干燥寒冷，夏季高温多雨[15]；年平均气温为13.1 ℃；年平均降水量为500～800 mm；人均水资源量为462 m3，为全国平均水平的五分之一，水资源承载能力与经济社会发展矛盾突出，是南水北调东线和中线工程的受水区[16]。黄淮海地区土层深厚，土质肥沃，农业规模大，主要种植小麦和玉米。地下水开采量的持续增加，导致一些地区的地下水位下降严重。

图1 研究区位置Fig.1 The location of the study area

2 数据与方法

2.1 数据来源

（1）GRACE 重力卫星数据。GRACE 计划是美国国家航空航天局和德国航空中心的合作项目，通过卫星监测地球重力场，估算陆地水储量的变化。本文选用美国德克萨斯大学研究中心（CSR）最新发布的RL06 Mascon 版本（https://www2.csr.utexas.edu/grace/）提供的重力卫星陆地水储量变化数据。相比于传统的球谐系数，Mascon 产品的最大特点是无需做任何后处理，可以直接提供格网产品。该数据以2004—2009 年的陆地水储量等效水深为基准，空间分辨率为0.25°×0.25°，时间分辨率为1 个月，时间序列从2003 年2 月至2021 年12 月，个别缺测时段采用奇异谱分析法[17]插补缺失值。

（2）GLDAS 陆面同化数据。全球陆面数据同化系统（GLDAS）包含4 个陆面模型（Noah、CLM、Mosaic、VIC），本文采用的是GLDAS-2.1 Noah 模型陆面同化数据，可解决内陆水域无效数据的问题。从中选取总冠层水量（Canoplnt）、雪水当量（SWE）、土壤含水量（SMS）和蒸散发（Evap）数据，其中SMS 为4 层土壤含水量数据之和（0～10 cm、10～40 cm、40～100 cm、100～200 cm）。时间长度及空间分辨率与GRACE 数据一致。

（2）降水和供水量数据。数据来源于水利部和流域部委水资源公报，时间序列为2003—2021 年。选取各流域降水、地表水供水量、地下水供水量、跨流域调水和其他水源等数据，其他水源包括污水处理回水、集雨工程供水、海水淡化供水等。

2.2 研究方法

2.2.1 地下水储量变化估算 陆地水储量包括地表水、地下水、土壤水、生物水、冰雪、植被冠层水等，其中地表水与生物水相对于其他成分，其多年变化量可以忽略不计[18-19]。因此，陆地水储量变化可表示为：

由此可得到地下水储量变化量的计算：

式中：ΔGWS为地下水储量变化（cm）；ΔTWS为陆地水储量变化（cm）；ΔSMS为土壤含水量变化（cm）；ΔSWE为雪水当量变化（cm）；ΔCanoplnt为冠层水量变化（cm）。陆地水储量变化由GRACE 重力卫星反演得到，土壤含水量、雪水当量和植被冠层水的变化由GLDAS 数据计算。

2.2.2 Theil-Sen Median 和Mann-Kendall 法 Theil-Sen Median 方法又称Sen 斜率估计，是一种稳健的非参数统计的趋势计算方法[20]，用于估计地下水储量的变化速率 β，计算公式为：

式中：x为 各时段的ΔGWS；i和j为时段月份数；Median 表示取中值。若β大于0，则表示地下水储量为增长趋势，反之则为下降趋势。

Mann-Kendall 检验是一种非参数的时间序列趋势性检验方法，广泛应用于水文气象数据分析，定义检验统计量S：

使用检验统计量Z进行趋势检验，Z值计算方法为：

式中：Svar为统计量S的方差。取显著性水平α=0.05，当|Z|≥1.96 时，表示该序列通过95%置信度检验，地下水储量变化趋势显著，反之变化趋势不显著。

2.2.3 随机森林模型 随机森林（Random Forest）是一种基于统计学习理论的组合分类智能算法。基本单元是决策树，使用自抽样（Bootstrap）法从原始数据集中随机重复采样，得到多个决策树模型，并将所有结果的平均值作为最终回归结果。根据基尼指数判断多个自变量对因变量的所做贡献，可以有效进行特征重要性排序，随机森林算法具有训练速度快、精度高、不易过拟合等优点，在水文地质学领域应用越来越广泛[21]。本文采用python 语言中sklearn 库建立随机森林模型，以地下水储量变化作为因变量，将影响地下水储量变化的因素作为自变量，计算影响因素的贡献率。

基尼指数代表了模型的不纯度，基尼指数越小，则不纯度越低，特征越好。其计算公式为：

式中：K为数据集分类总数；pk为某个子集属于第k类的概率。

特征A在内部某个节点m的重要性Vim,m为分枝前后基尼指数的变化量：

式中：Gini,m为节点m处的基尼指数，Gini,l、Gini,r为分枝后出现新节点处的基尼指数；l、r为节点m处分枝后的新节点。

则特征A在随机森林模型中基尼指数评分为：

式中：Vim,i为第i棵决策树各节点处分支前后基尼指数变化量，n为节点数。

3 结果与分析

3.1 时间变化特征

根据GRACE 和GLDAS 数据，采用式（2）计算研究区各网格地下水储量的逐月变化。整个黄淮海地区相对于基准期，地下水储量在2003—2021 年总体波动范围为−24.94～5.99 cm（见图2），整体呈减少趋势，近20 年内共减少23.18 cm，平均每年减少1.22 cm。其中在2003—2010 年平均减少0.45 cm/a，在2011—2019 年间减少速度加快为1.68 cm/a，其原因可能与经济快速发展及人口增长有关，在2020 年以后以1.05 cm/a 的速度开始增加。

图2 黄淮海地区地下水储量时间变化趋势Fig.2 Temporal fluctuations in groundwater storage within the Huang-Huai-Hai region

逐一分析比较黄淮海地区内各二级子流域地下水储量变化（见图3），具体变化速度见表1。从表1 和图3 可见研究期内各子流域的地下水储量均处于持续减少状态，海河流域平均下降速率为1.75 cm/a，淮河流域平均下降速率为0.57 cm/a。在2003—2010 年间，地下水储量下降最快的是海河南系，达1.94 cm/a，而淮河流域整体处于缓慢上升阶段。2011—2019 年研究区地下水储量亏损程度扩大，亏损地区增多，除淮河下游变化不明显，各子流域呈现不同程度的下降趋势，花园口以下地区以2.92 cm/a 的速率下降，海河南系下降速率为2.86 cm/a。2020 年后，部分子流域地下水储量开始上升，花园口以下地区上升速率为7.22 cm/a。

表1 黄淮海地区子流域地下水储量变化速度Tab.1 The rates of change in groundwater storage within the Huang-Huai-Hai region

图3 黄淮海地区各子流域地下水储量变化趋势Fig.3 Temporal fluctuations in groundwater storage within the sub-basins of the Huang-Huai-Hai region

各子流域地下水储量变化在年内变化表现出差异性（图4）。滦河及冀东沿海、海河北系地区地下水储量逐月变化较为稳定，山东半岛沿海诸岛、淮河下游地区地下水储量在4—8 月份处于高位。海河南系、徒骇马颊河、花园口以下、淮河上游等主要农作物耕种地区的地下水储量下降明显，且在6—7 月份明显偏少。黄淮海地区水储量持续减少主要是因农业灌溉进行的大规模开采地下水[22]，用于农作物灌溉的地下水供水量占地下水开采总量的80%。黄淮海地区的农作物以冬小麦和夏玉米为主，5—8 月为作物需水关键期，需水量占整个生长发育过程的79%[23]，地下水下降最快，地下水储量处于最低状态。冬小麦开始大量播种时间为10 月中旬，越冬期不需要灌溉，降水基本能满足需要，此时地下水储量得到恢复。提高农业用水效率和减少地下水开采量是缓解区域地下水资源短缺的重要措施。

图4 黄淮海地区地下水储量逐月变化Fig.4 The monthly average of groundwater storage in the Huang-Huai-Hai region

3.2 空间变化趋势

从空间分布看，研究区地下水储量变化存在明显的空间差异性，整体上由东至西呈现下降幅度增大趋势，变化范围为−21.9～2.27 cm，负值为相对于基准期处于减少状态。图5 表明地下水储量减幅最大的是海河南系，主要集中在邯郸、安阳、鹤壁，淮河流域上游较为盈余；最低值出现在鹤壁市，最高值出现在淮安市。位于海河南系的白洋淀流域是华北最大的淡水湖，已成为雄安新区发展的重要生态水体，该地区地下水储量最低值为−10.75 cm，最高值为−4.07 cm，亏损较为严重。

图5 黄淮海地区多年平均地下水储量变化分布Fig.5 The spatial pattern of groundwater storage within the Huang-Huai-Hai region

基于Sen 斜率估计和Mann-Kendall 法分析黄淮海地区2003—2021 年地下水储量变化趋势，计算出每个格网的β斜率和显著性Z值，研究区88.87%区域通过置信水平0.05 的检验，其中94.2%的区域呈现出显著的负向变化，即地下水储量减少。为了进一步比较各个地区地下水储量的变化情况，将研究区域划分为极显著增加、显著增加、微显著增加、不显著增加、不显著减少、微显著减少、显著减少、极显著减少等8 个等级（图6）。从图6 可以看出：黄淮海地区地下水储量减少速率整体表现出自南向北逐步增大的特点，海河南系大部分区域处于极显著减少状态，水储量增加的地区位于淮河流域上游及中游部分区域。

图6 Sen+M-K 计算地下水储量变化趋势分布Fig.6 The trend of variations in groundwater storage calculated using Sen+M-K methods

3.3 地下水储量变化成因

3.3.1 降水、蒸散发变化影响 降水入渗补给是影响地下水储量变化最直接的因子。淮河流域降水丰沛，年均降水量为873.3 mm，充足的降水可促进地下水补给，有利于恢复地下水储量。海河流域年均降水量为532.7 mm，降雨径流经水库截留等损失后流域下渗较少，对地下水储量变化影响的并不明显（图7）。同时，由于海河流域农业部门耗水最多，降水量变化会影响灌溉量，枯水年灌溉量增多，地下水开采量会随之增多。蒸散发过程会消耗土壤水分，减少地下水的补给。黄淮海地区多年平均蒸散量表现为增加趋势，淮河流域的蒸散量由2003—2015 年的722.4 mm 增加至2016—2021 年的782.7 mm，海河流域的蒸散量由533.9 mm 增加至597.3 mm，花园口以下的蒸散量由629.7 mm 增加至678.0 mm。

图7 降水、蒸散发与地下水储量变化关系Fig.7 The interrelationship among precipitation, evapotranspiration, and groundwater storage

降水在黄淮海地区时间空间上的分布差异明显。淮河流域降水年际变化更大，直接影响地下水储量变化速度，地下水储量变化与降水变化趋势具有较明显的一致性。而海河流域降水对地下水储量变化的影响不显著，整个海河流域年降水量与地下水储量之间的相关系数仅为0.02，滦河及冀东沿海、海河北系甚至出现负相关，这表明影响海河流域地下水储量变化的主要因素不是降水量，而是人类开采的地下水量，这与之前学者研究结果一致[24]。

3.3.2 供水量变化影响 流域供水来源包括当地地表水、地下水、跨流域调水等，其中淮河流域外调水为引长江水与黄河水，海河流域外调水在2015 年之前为引黄河水，南水北调引长江水主要从2015 年开始。淮河流域地下水供水比例逐年缓慢减少，从2003 年占总供水量的32.3%降至2021 年的22%，主要原因是跨流域调水量呈增加趋势；地表水多年平均供水比例为55.8%，是其最主要的供水来源，地下水次之，平均占比26.9%（图8）。黄河花园口以下地区地下水供水比例呈减小趋势，但总供水量呈增加趋势。

图8 黄淮海地区年供水量变化Fig.8 The annual fluctuations in water supply within the Huang-Huai-Hai region

海河流域地下水是主要的供水来源，2015 年以前地下水供水占比64.2%，2015 年南水北调工程通水后开始显著减少，2019—2021 年已降到39.1%，但仍是海河流域最主要的供水来源。南水北调工程通水后，海河流域供水结构发生了改变，地下水供水量明显减少，但由于地下水多年持续超采严重，目前仍维持在近150 亿m³/a，地下水储量减少速率有所减缓，但下降趋势未得到根本遏制[10]。

3.3.3 成因分析 以降水量、蒸散量、地下水供水量和跨流域调水量为自变量建立随机森林模型，将计算所得影响因素重要性评分进行归一化处理，定量分析比较不同因子对黄淮海地区地下水储量变化的相对影响程度（图9）。结果表明，淮河流域地下水储量变化的影响因子相对贡献率为降水＞蒸散发＞地下水供水＞跨流域调水，而海河流域影响因子的相对贡献率为地下水供水＞跨流域调水＞降水＞蒸散发，黄河流域花园口以下的相对贡献率为地下水供水＞蒸散发＞降水。各个流域地下水储量变化的主导因素不同，降水是淮河流域的主导因素，贡献率28.3%；地下水供水是海河流域和黄河流域花园口以下地区的主导因素，贡献率分别达52.6%和58.0%。

图9 黄淮海地区地下水储量变化影响因素的相对贡献度Fig.9 The relative contributions of factors influencing variations in groundwater storage within the Huang-Huai-Hai region

降水和蒸发等自然因素对淮河流域地下水储量变化的影响程度显著大于海河流域，而地下水供水是影响海河流域和黄河流域花园口以下地区地下水储量下降的主导因素。研究表明人为因素对黄淮海地区地下水储量变化影响逐渐增大[25-26]，2003—2021 年间，地下水消耗远远超出了地下水系统的自然补给能力，导致整个海河流域和淮河流域北部地下水储量大幅下降。海河流域跨流域调水的相对贡献率大于淮河流域，随着地下水超采治理及南水北调供水工程持续推进，地下水储量下降速率近年来有所减缓，北京、天津等受水区城市地下水位上升明显[27-28]，但农业区地下水位仍在下降。

4 结 语

（1）2003—2021 年间，黄淮海地区地下水储量呈波动减少趋势，平均变化速率为−1.22 cm/a，其中海河流域为−1.75 cm/a，淮河流域为−0.57 cm/a；海河南系亏损程度最严重，淮河上游和下游变化不显著。2010 年代的下降速率高于2000 年代，随着南水北调工程调水量的增加，地下水储量下降速率有所减缓，沿线部分城市地下水储量有所回升。

（2）地下水储量变化存在明显的空间差异性，海河流域亏损程度大于淮河流域。黄淮海地区94.2%的区域地下水储量呈下降状态，下降速率整体表现出自南向北逐步增加的特征；海河南系为极显著减少趋势，淮河流域上游及中游部分地区呈增加趋势。

（3）黄淮海地区地下水储量变化的主导因素不同，降水对淮河流域地下水储量变化具有主导作用，而影响海河流域和黄河流域花园口以下地区最主要的因素是地下水供水量，多年持续超采造成地下水储量大幅减少。南水北调工程通水后海河流域地下水开采量减少，改善了区域供水结构，地下水储量下降得到缓解，但还需进一步提高农业用水效率，加强地下水资源管理。