GRACE重力卫星监测煤矿开采区地下水变化研究

李 舒，师鹏飞，谷晓伟，付新峰，倪深海，张 楠

（1. 河南省黄河流域生态环境保护与修复重点实验室，河南 郑州 450003；2. 黄河水利委员会黄河水利科学研究院，河南 郑州 450003；3. 河海大学 水文水资源学院，江苏 南京 210098；4. 南京水利科学研究院，江苏 南京 210029）

1 研究背景

地下水是支撑陕西省社会经济发展和维系生态系统健康的重要水源，其储量变化对陕西省农业、城镇生活、工业及生态环境的用水量响应较为敏感。而陕西省作为煤炭资源主产区，煤矿开采这一人类活动对地下水水储量的扰动作用也十分明显。采用传统的地下水水位、给水度、有效孔隙度等水文地质参数推测地下水水储量的变化需要开展大量的地质钻孔工作，既费时又费力［1-2］。区域内的地下水水位观测井布设的密度较小，而且有些煤矿的观测井并未纳入政府统一管理，含水层受采矿扰动后，通过地下水模型模拟区域地下水水储量变化亦存在难度［3-5］。

近年来，卫星遥感技术的发展，为监测地下水水储量变化提供了新的途径。GRACE（Gravity Recovery and Climate Experiment）重力卫星于2002年3月发射升空，该卫星可建立月尺度下的全球重力场模型，进而反演全球水储量变化的分布情况［6］。全球水储量变化主要由陆地水储量变化引起，陆地水储量主要是由地表水、土壤水、地下水和冰川融雪等4 部分组成。近年来，国内外学者利用GRACE重力卫星成功反演了全球和华北平原、松花江流域等大、中尺度流域的陆地水储量变化情况［7-9］，约旦、青藏高原、加利福尼亚中央山谷等区域内地下水储量变化情况［10-12］，以及冰川融雪的水储量变化情况［13］。这些研究主要采用以下方法反演地下水和冰川融雪的水储量变化情况，并评估GRACE监测数据在研究区的适用性［14］。首先，利用气候水文模型模拟地表水、土壤水储量变化，其次，从GRACE反演出的陆地水储量变化量中扣除地表水和土壤水储量变化量，最后，将扣除后的陆地水储量变化量与模型模拟结果或是实测数据进行比较，评估反演结果的合理性和准确性。

国内外学者还利用GRACE 数据反演了不同气候条件下区域的陆地水储量和地下水水储量。例如，Naveed等［15］采用GRACE反演出了月尺度的地下水水储量变化，可为湿润气候条件下印度河上游的水行政主管部门提供决策所需的基础数据。Ouyang等［16］通过比较GRACE反演出的地下水水储量变化和农田扩张面积的关系，发现GRACE数据可支撑半湿润区域下的冻融作物水资源管理。Pascal等［17］发现GRACE 无法监测墨西哥中部农业灌溉区的退水情况。Ahmed 等［18］利用GRACE 数据反演了地处干旱区的努比亚沙岩含水层水储量的变化量，发现该含水层水储量以极快的速率下降。曹艳萍等［19］反演了黑河中上游地下水储量，发现与实测地下水数据的相关性较好。Soumendra 等［20］利用印度2005—2013年15000多口地下水观测井数据，验证了基于质量密集模型和球谐系数模型解译GRACE重力场反演地下水水储量的适用性，研究结果表明采用质量密集模型解译法在区域尺度反演效果更好。另外，还有一些学者在德克萨斯州［21］、中亚［22］、乍得湖流域［23］和底格里斯-幼发拉底河流域［24］等干旱区开展了一系列的应用研究，结果表明GRACE可有效监测干旱区地下水水储量变化量。除此之外，Liu等［25］利用GRACE 结合统计模型预测了黄河下游枯水期的流量，减少了预报误差，增长了预见期。Liu等［26］将GRACE与干旱指数结合分析了海河流域干旱情况。

当前，大规模、大范围煤矿开采活动导致流域地表-地下水循环发生了重要变化［5］，仅采用传统的地下水水位观测井数据尚无法全面刻画区域尺度下地下水的变化情况。GRACE在流域和区域地下水储量变化监测中的成功运用，特别是在干旱区的成功运用，有力指导了流域水资源的科学配置，为研究大规模大范围煤矿开采活动下的地下水储量变化提供了新的途径和先进手段。为此，本文选择典型的煤矿开采区—窟野河流域，采用GRACE重力卫星、结合GLDAS数据，监测窟野河流域地下水水储量变化，以揭示煤矿开采影响下的窟野河流域地下水水储量变化，并拓展GRACE的应用范围。

2 研究区概述

窟野河（图1）是黄河中游的一级支流，河长242 km，流域面积8706 km2。该流域位于黄土高原和毛乌素沙漠之间，地理范围为北纬38°22′—39°50′，东经109°28′—110°45′，主要由风沙区、基岩出露区和黄土覆盖区组成，其中，流域上中游以前两种地貌为主，流域下游以后一种地貌为主［27］。受大陆季风气候影响，流域多年平均温度为7.9℃，降水量为415 mm，降水多集中于6—9月，占年降水总量的75%～81%。流域多年平均径流量为5.44亿m3，多集中于7—10月，占年径流量的57%～74%［28］。

图1 窟野河流域气象站、水文站及矿区分布

通过对15个煤矿开采区的水文地质资料分析发现，矿井涌水的主要来源为萨拉乌苏组含水层（第四系潜水）和侏罗系烧变岩含水层。萨拉乌苏组含水层主要分布在流域中游，而在其他地方主要呈分散状、薄层分布，分布面积为336.58 km2，补给源主要来自降雨，入渗系数在0.03 ～0.25之间。侏罗系烧变岩含水层主要分布在窟野河西边的支流两侧，补给源主要来自萨拉乌苏组含水层，构造有利部位可形成强富水区，其渗透系数一般可达100 m/d。经测算，研究区内煤矿的矿井涌水量一般为500 ～1000 m3/d［29］。2009年研究区内的煤矿开采量达到2.03亿t，研究发现，流域内地下水发生了较大变化，径流量为1.3亿m3，是1954—2019年的最小径流量。2009年之后，研究区开始大规模应用“保水开采”技术。因此，本研究以2009年为典型年，分析煤矿开采影响下的地下水储量变化。

3 数据与方法

3.1 数据来源GRACE 数据来源于CSR（Center for Space Research，The University of Texas at Austin）、GFZ（GeoForschungsZentrum）和JPL（Jet Propulsion Laboratory）三个机构。数据为level-3的球谐系数产品（RL-05）［30-31］，空间尺度为1°×1°，时间尺度为月，时间范围为2002年4月到2017年1月。该数据已经过高斯平滑和去条纹滤波技术处理［32］，由于电源问题产生的空白数据，通过线性差值方法处理［33］。GLDAS数据来源于NASA官网中的NOAH［34］、VIC［35］、CLM［36］和Mosaic［37］水文模型的模拟结果，为匹配GRACE数据的时空分辨率，采用分辨率为1°×1°的月数据。

降雨数据来源于“中华人民共和国水文年鉴”，包括2004—2009年20个雨量站、1个水文站实测的日降雨资料，蒸发数据来源于气象部门，包括2004—2009年3个气象站日实测蒸发量。

3.2 方法地下水水储量变化通常受气候变化和人类活动双重影响，为阐明煤矿开采活动对地下水水储量的影响，采用以下步骤进行处理。

（1）初步分析监测要素合理性。分析GRACE 反演的陆地水储量和站点监测（降雨、蒸散发）要素的相关性，初步验证GRACE监测结果的合理性。

（2）计算地下水储量变化。将GRACE反演的陆地水储量减去GLDAS模拟的地表径流量、土壤水储量、叶冠层储水量、冰雪水储量，得到地下水水储量变化。具体计算方法如下：

受GRACE产品空间分辨率影响，选取覆盖研究区的2个网格数据（网格分布在东经109°—111°，北纬38°～40°之间，中心点分别为东经110°、北纬38°30'和东经110°、北纬39°30'），并以此两个网格2004年—2009年为基准期的陆地水储量变化的平均值作为研究区的陆地水储量变化量。尽管2个网格已超出研究区范围，但是研究期内超出研究区范围的地方并没有大规模煤矿开采活动，因此，对研究结论影响较小。根据水量平衡原理，从GRACE反演的陆地水储量中分离计算出地下水水储量变化量，陆地水储量的数学表达式如下：

式中：TWS为陆地水储量；SMS为土壤水储量；SWS为地表水体储水量；SWE为冰雪水量；CWS为叶冠层储水量；GWS为地下水水储量。

GRACE重力卫星反演出的陆地水储量是变化量，根据水量平衡原理，式（1）可改写为：

式中Δ 表示某一时期的变化量，此处为月变化量。

流域内SWS变化主要包括河流、湖泊、水库中水储量的变化，许多研究忽略了这项的变化量，会导致反演出的地下水水储量的变化量与实际情况不一致［38-39］。但对于窟野河流域，由于煤矿开采导致地表水体储水量变化相比于地下水和土壤水储量变化要小一个量级，而且GLDAS中的地表水体储水量未考虑湖泊与水库［40-41］，因此，可以忽略地表水体储水量变化项。地下水水储量变化ΔGWS为：

式中：ΔTWS项为来源于三个机构（CSR、GFZ和JPL）发布的GRACE数据的平均值，ΔSMS、ΔSWE、ΔCWS来源于GLDAS 4个水文模型（NOAH、VIC、CLM和Mosaic）模拟结果的平均值，对4种模型结果取平均值是当前研究中较为常用的方法，取平均值后的结果具有较强的可靠性［42］。

（3）验证结果合理性。采用相关系数法，将步骤（2）计算的地下水储量变化ΔGWS与李舒等［43］采用地表-地下水耦合模型模拟的煤矿开采影响下的地下水水储量结果进行相关分析，验证GRACE监测煤矿开采区地下水储量的合理性。此外，由于地下水与土壤水间通常存在很强的水力联系，通过分析土壤水储量变化ΔSMS与地下水储量变化ΔGWS间的相关性，验证地下水储量计算结果的合理性。具体计算方法如下：

式中：ΔGWSi为地下水储量月变化值；为2009年地下水储量变化月变化值的平均值。当计算ΔSMS与地下水储量变化ΔGWS间的相关性时，Yi和分别为GLDAS 模拟的ΔSMS逐月变化值、ΔSMS2009年月变化值的平均值。当计算地下水储量变化ΔGWS与采用地表-地下水耦合模型模拟的地下水水储量结果的相关性时，Yi和分别表示采用地表-地下水耦合模型模拟的地下水水储量逐月变化值、模拟期月变化的平均值。

4 结果分析

4.1 GRACE反演的陆地水储量合理性分析将CSR、GFZ和JPL三个机构发布的数据进行比较，从图2中可以看出，三个机构提供的2009年反演结果差距较大，其中，JPL 1—7月的陆地水储量比另外两个机构数据的平均值大了16.7 mm，而GFZ机构8—10月的陆地水储量比另外两个机构数据的平均值大了7.4 mm。可知，不同机构发布的GRACE 产品数据存在较大差异性，换言之，GRACE 反演的陆地水储量数据存在较大不确定性。研究表明［45］采用三个机构陆地水储量的平均值作为GRACE最终反演结果，是减小误差的有效方法。对三个机构的数据进行平均处理后，陆地水储量在1—3月，3—9月，9—12月均呈现减少-增加的趋势。与2004—2009年的基准期相比，2009年陆地水储量减少量为15.5 mm/月。

图2 2009年GRACE反演的陆地水储量结果

由于GRACE反演出的陆地水储量受到气候变化和人类活动的双重影响，所以先定性分析2009年研究区降雨量和蒸发量相较于基准期（2004—2009年）的变化趋势。从图3中可以看出，2009年研究区蒸发量为2280 mm，大于基准期年平均蒸发量2143 mm；降雨量为374 mm，略小于基准期的年平均降雨量388 mm；相对基准期，2009年降水量小、蒸发量大，GRACE反演的陆地水储量相较于基准期也较小，说明从年尺度的变化趋势上看，反演结果基本合理。

图3 窟野河2004—2009年降雨量和蒸发量

2009年月降雨量和蒸发量的变化值与GRACE 反演结果的对比，如图4 所示。仅考虑气候变化，当蒸发量大于（小于）降雨量时，陆地水储量应呈现减少（增加）的趋势。但是，4月、7月、8月和11月反演的陆地水储量的变化量与月降雨量和蒸发量与上述规律正好相反，说明GRACE反演结果与月降雨量和蒸发量的趋势匹配度一般。工业用水、生活用水、水保措施以及煤矿开采等人类活动对陆地水储量也有较大的影响，进一步说明GRACE反演出的陆地水储量的变化趋势是合理的。

图4 研究区2009年GRACE反演结果与月降雨量和蒸发量比较（紫色部分说明月降雨量与蒸发量变化结果重叠）

4.2 GRACE与GLDAS模拟结果对比分析

4.2.1 GRACE反演的陆地水储量与GLDAS模拟结果对比分析 前文定性分析了GRACE反演的陆地水储量的合理性，本节结合GLDAS水文模型定量分析GRACE反演结果的合理性。GLDAS给出了4种水文模型计算出的土壤水量，冰雪水量和叶冠层储水量等结果，通过与2004—2009年基准期相比，2009年研究区GLDAS的模拟结果见图5。

图5 GLDAS模拟结果与GRACE反演结果

由图5可知，4种模型模拟结果的变化趋势基本一致，变化值存在差异。例如，Mosaic模型模拟出的土壤水量变化量在1—4月与其他三种模型差距较大，VIC模拟出的土壤水量变化量在5—6月与其他三种模型差距较大，4种模型在9—12月均不一致；VIC模拟出的冰雪水量变化量与其他三种模型模拟结果均不一致；Noah模拟出的叶冠层储水变化量在11月与其他三种模型差距较大，CLM模拟出的叶冠层储水变化量在8月与其他三种模型差距较大。

为了减少模型模拟结果的不确定性，当前研究中多采用对4种模型模拟结果取平均值的方法，以平均值作为GLDAS模拟结果［42］。从图5（a）中可以看出土壤水量变化量的峰值出现在9月，为46.9 mm，谷值出现在6月，为-17 mm，2009年土壤水量平均增加了14 mm/月。图5（b）中可以看出冰雪水量变化量的峰值出现在11月，为0.97 mm，2009年冰雪水量平均减少了0.13 mm/月。图5（c）中可以看出叶冠层储水变化量的峰值出现在8月，为0.07 mm，谷值出现在1月，为-0.03 mm，2009年叶冠层储水变化量基本无变化。

将GRACE 反演的陆地水储量变化量与GLDAS 模拟结果进行比较（图5（d）），可知两者的变化趋势基本一致，但是两者变化值存在较大差异，相关系数为0.45，这是因为GLDAS中的计算结果未包含地下水储量的变化情况。结果表明，研究区地下水水储量的变化是造成GRACE与GLDAS变化值不一致的主要原因。

4.2.2 GRACE反演的地下水水储量与GLDAS模拟结果对比分析 根据式（3），利用GRACE联合GLDAS反演出的研究区地下水水储量变化结果，如图6所示。将地下水水储量与GLDAS模拟结果进行比较，发现地下水水储量与GLDAS模拟值在趋势上呈反对称分布的特点，相关系数达到0.84，并且通过95%的置信度检验。而研究区土壤水量在GLDAS 模拟结果中的比重达到95%左右。由此可知，研究区土壤水变化量与地下水水储量变化量存在着较强的相关性。

图6 GRACE反演地下水水储量与GLADS模拟结果的比较

GRACE反演出的研究区2009年地下水水储量的减少量为29.4 mm，年内分布的峰值出现在6月，为-8.6 mm，谷值出现在9月，为-63.2 mm。变化趋势呈现“N”型，即2—6月呈现增加的趋势，随后呈现减少的趋势，到9月底后又呈现增加的趋势。地下水水储量变化幅度很明显，这与研究区40%的耗水来源于地下水的现状也是相符的，特别是夏季耗水量激增也造成了6—9月地下水水储量减少。

4.3 GRACE监测煤矿开采区地下水水储量的结果为了进一步分析GRACE反演煤矿开采区地下水水储量变化的适用性，将反演结果与李舒等［43-44］利用地表-地下水耦合模型（SWAT-VISUAL MODFLOW）计算出的2009年研究区煤矿开采对地下水影响的结果进行比较。地表-地下水耦合模型基于水量平衡原理，设置了有、无煤矿开采两种情景，计算煤矿开采对地下水水量的影响，比较结果见图7。

图7 GRACE反演的地下水水储量变化量与耦合模型模拟结果比较

对GRACE反演的地下水水储量与耦合模型模拟结果进行相关性分析，结果表明两组数据相关系数为0.64，并且通过了95%的置信度检验，说明煤矿开采是导致研究区地下水变化的一个重要因素。此外，王童等［46］研究发现窟野河流域植被覆盖度2004年至2009年，由45%快速增加至58%，而研究区2009年较基准期的蒸发量也是增加的，植物的蒸腾作用成为地下水排泄的一个重要途径，说明研究区地下水水储量2009年呈现减少趋势的原因是大规模的煤矿开采和植树造林造成的。

但是，两组数据月尺度下的等效水深除5月、7月和11月较接近，其他月份差距较大。年尺度下，GRACE 反演出的地下水水储量变化量为-29.4 mm/月，大于耦合模型模拟结果的-23.2 mm/月。从结果来看，GRACE 反演的煤矿开采区的地下水水储量的月均值变化量较为合理，但是GRACE 反演出的逐月值结果较差。

5 讨论

5.1 煤矿开采对地下水水储量的影响机理分析通过对研究区内的榆林市2004—2009年的耗水量进行分析，发现其中约40%的耗水量来源于地下水，将耗水量按榆林市面积换算为等效水深，2009年地下水耗水量相较于基准期减少了0.04 mm/月［47］，耗水量减少但是地下水水储量却呈现减少的趋势，说明直接取用地下水对地下水水储量的影响较小。而研究区土壤水变化量与地下水水储量变化量存在着较强的相关性，直接说明研究区内的煤矿开采造成的导水裂隙带导通了土壤水与地下水交换的通道，加强了土壤水与地下水的交换频率和交换量，这与国内相关研究成果也是一致的［5］。

煤矿开采过程中，导水裂隙带由下至上逐渐发育，最终发育至地表，形成地裂缝，最终完全导通煤层以上的含水层，在此过程中为了采矿安全需进行矿井水的疏干外排，进而加速了地下水与地表水的交换频率和交换量。而耦合模型模拟出的煤矿开采情景下的地下水水量也呈下降趋势，第一阶段（3—10月）缓慢下降，下降率为3.4 mm/月；第二阶段（10—11月）迅速下降，下降率为15 mm/月；第三阶段（11—12月）稳定下降，下降率为0.2 mm/月，反映了导水裂隙带发育过程中，地下水水储量的变化特征，即采矿初期矿井疏干水缓慢增加，采矿中期矿井疏干水迅速增大，采矿末期矿井疏干水逐渐稳定。而GRACE 反演的地下水水储量与耦合模型模拟结果的相关系数为0.64，间接说明GRACE反演结果也能反映出导水裂隙带加强了地表水-地下水之间的交换频率和交换量。

5.2 GRACE监测煤矿开采区地下水水储量的适用性分析GRACE反演的地下水水储量变化幅度很明显，这与研究区40%的耗水来源于地下水的现状是相符的，特别是夏季耗水量激增也造成了6—9月地下水水储量减少。因此，说明GRACE反演的煤矿开采区的地下水水储量的月均值变化量还是较为合理的，但是GRACE反演出的逐月值与耦合模型模拟结果相比较差，主要原因有以下3个：（1）研究区除了煤矿开采外还有大面积的植树造林活动，而耦合模型并未考虑该影响；（2）GRACE 反演陆地水储量时选择的2 个网格中还包含了秃尾河流域，而耦合模型模拟范围仅包含窟野河流域，但是2010年后秃尾河流域内的煤矿开采活动才大规模出现，对结果影响较小；（3）耦合模型模拟结果是2009年的煤矿开采对地下水的影响量，与GRACE反演结果的基准期不一致，但是煤炭产量在2004—2009年间增长速度较缓慢，因此对结果影响较小。

6 结论与展望

通过比较GRACE重力卫星反演的陆地水总储量和地下水水储量的变化量与GLDAS中4种水文模型模拟结果以及降雨、蒸发等气象因素的相关性，分析了反演结果的合理性。将GRACE重力卫星反演的地下水水储量的变化量与地表-地下水耦合模型模拟结果进行了比较，分析了GRACE重力卫星在监测煤矿开采区地下水水储量的适用性。研究结果表明：

（1）GRACE 重力卫星可反演干旱半干旱区大规模煤矿开采区地下水水储量的月均值变化量，但是无法准确的反演逐月的地下水水储量的变化量。

（2）煤矿开采区的土壤水与地下水存在着较强的相关性，煤矿开采造成的导水裂隙带导通了土壤水与地下水交换的通道，加强了土壤水与地下水的交换频率和交换量。

（3）2018年5月GRACE-Follow on 计划实施，新一代的卫星较上一代增加了激光干涉测距仪，该仪器可大幅增加数据观测精度，未来随着数据精度的提高必将使其成为监测流域煤矿开采对地下水影响的一个新手段。

GRACE在大流域反演地下水水储量已有很多成果，但是在中、小流域反演地下水水储量的研究还较少，主要是因为其空间分辨率较低。因此，如何提高中、小流域尺度下监测地下水及其变化的精度还有待进一步研究。