强烈人工干预下的地下水流动系统变化，已经成为影响局部乃至区域生态水文平衡的主要驱动力^[1].干旱-半干旱地区的煤炭资源开发往往会改变区域地下水循环模式，破坏地下水与湖泊、河流等地表水体之间的水力联系^[2-5]，而湖泊对维持生态平衡、保护区域栖息物种等具有重要意义^[6-11].国内外学者们对于煤碳资源开采对湖泊的影响进行了大量研究.赵宁等和杨立彬等^[12-13]研究发现鄂尔多斯盆地的风沙区依赖地下水补给的湖泊多达500个，由于建坝、打井抽取地下水和煤炭资源开采等人为因素的影响，引发区域地下水位大幅下降、部分湖泊湿地变干，如中国第一大沙漠湖泊-红碱淖的水域面积减少了26.66 km²；王文科等^[14]研究发现荒漠化地带的陕北矿区内河流及湖泊40 % ~70 %的流量来自地下水的补给，而在大规模煤矿山开发过程中，地下水普遍下降10~12 m，导致区域河水流量及湖泊水面大幅减少；马雄德等^[15]研究发现，在1990 -2011年期间，沙漠滩地的陕西榆神府矿区内水体总体呈现减少的趋势，境内水体面积从162.63 km²降低到103.84 km²；Pöschke等^[16]对斯坦克林湖的研究中发现只有一部分集水区的地下水补给湖泊，而其余的则向更深的含水层泄漏；Rudnick等^[17]研究发现湖泊的地下水排放可以在湖泊的水平衡中发挥主要作用；Hajati等^[18]对汉普顿湖的研究中发现在低于平均降水量的情况下，该湖经历了一次倒流，在此期间该湖主要转变为补给湖，然后于2012年初返回一个过流湖.

综上所述，学者们对于湖泊与地下水之间的演化关系进行了一定的研究，但定量^{[16, 19-20]}评价地下水对湖泊的补给作用是很困难的，现有研究主要集中于定性描述煤矿开采引起的地下水位变化对湖泊水量的影响，有鉴于此，本研究以内蒙古呼伦贝尔草原伊敏露天煤矿为研究对象，在对区域地下水流系统调查的基础上，利用遥感技术、地下水流系统理论及水均衡原理，探讨煤矿开采活动对地下水湖泊系统的演变影响，构建地下水位与湖泊面积之间的关系模型，定量揭示采矿引起的地下水位变化对湖泊面积(水量)的影响.

1 研究区概况 1.1 地理地貌

研究区位于大兴安岭西坡，伊敏盆地南部，伊敏河中下游地区，隶属于内蒙古自治区呼伦贝尔市鄂温克自治旗，面积约100 km².伊敏煤田地貌呈一半封闭型的盆地，东西两侧为丘陵，南部为台地，盆地内为冲积平原，海拔在651~673 m，南高北低，东西两侧稍高，中心低洼^[10](图 1).

1.2 气象水文

研究区属于典型的中温带大陆性季风气候，四季分明，冬冷夏凉.年均气温-1.9℃，境内降水自东南向西北递减，多年年均降水量为354.73 mm，年均蒸发量为1318.3 mm.

区内最主要的河流为伊敏河，苇子坑河等其余河流均为其支流，伊敏河由南向北从盆地东部穿过，全长2000 km，流域面积可达9000 km²，河床宽约60 m，水深0.50~2.50 m.目前伊敏露天煤矿区内地下水补给主要依赖于大气降水、伊敏河以及盆地内零星分布的湖泊.天然状况下，伊敏河水与地下水的关系为丰水季节伊敏河水补给地下水，枯水季节地下水补给河水.在煤矿开采的影响下，伊敏河受影响河段与地下水的关系是伊敏河水补给地下水.由于伊敏矿区疏排地下水，改变了区域地下水循环特征.目前地下水排泄方式主要是向露天矿疏干区排泄.地表水则以径流的方式，汇集于盆地低洼处，或渗透于地层中或以蒸发的形式排泄^[21].

2 数据来源与处理方法

遥感影像采用Landsat MSS、TM和OLI数据(表 1)；水文地质资料来源于黑龙江省伊敏煤田精查报告、伊敏煤矿水文地质调查报告、伊敏煤矿矿区地下水位动态监测报告及课题组2016 -2018年的实地调查；气象数据来源于当地气象局.

水体指数法是实现水体快速提取的简单便捷方法之一^[22].伊敏盆地的湖泊数量较多，大小不一，故采用水体指数法并结合目视解译的方式进行提取.本文利用在ENVI5.1+IDL平台开发出的归一化水体指数(Normal difference water index，NDWI)从上述4期Landsat系列遥感数据中提取伊敏盆地水体信息并绘制湖泊水体图，然后结合实地调查和收集的数据，分析煤矿开采35年来地下水湖泊系统的演化规律.

3 区域地下水湖泊关系模型 3.1 湖泊水均衡方程

湖水的补给包括降雨、可能的地表径流及地下水补给，湖泊的排泄主要是蒸发、可能的地表径流排泄以及工农业用水.因此，一般湖泊的水均衡方程为：

$ \frac{{\partial V}}{{\partial t}} = A\left( t \right)P\left( t \right) + f\left( t \right) - A\left( t \right)E\left( t \right) - Q\left( t \right) - {W_{\rm{q}}}\left( t \right) $

(1)

式中，V为湖泊的湖水体积，m³；∂V/∂t为单位时间内的湖泊水体体积变化率，m³/a；A为t年份的湖泊水面面积，m²；P为t年份的年降水量，m；f (t)为产汇流函数；E为t年份的年蒸发量，m；Q为t年份的湖泊人工抽取水量及地表径流排泄量，m³/a；W_q为单位时间内地下水与湖泊的交换量，m³/a.

3.2 草原圆台型湖泊面积与地下水位关系方程

前人研究表明直接测量地下水与湖泊的交换量不仅耗时、耗力还容易出错^[23]，本次研究利用水均衡原理、达西定律及几何学知识建立地下水位湖泊面积之间的关系方程.由于研究区湖泊四周地势平坦，形态大多近似圆台体，故以圆台体湖泊类型为例，湖泊水体体积示意见图 2.

湖泊水体体积可表示为：

$ V = \frac{1}{3}\left( {{\rm{ \mathsf{ π} }} \cdot r \cdot R_1^2 \cdot {H_{\rm{A}}} - {\rm{ \mathsf{ π} }} \cdot R_2^2 \cdot {H_{\rm{B}}}} \right) = \frac{1}{3}\left( {A \cdot {H_{\rm{A}}} - B \cdot {H_{\rm{B}}}} \right) $

(2)

则H_A为：

$ {H_{\rm{A}}} = \frac{{{H_{\rm{B}}}}}{{\sqrt {B/A} }} $

(3)

式中，V为湖泊水体体积，m³；R₁为湖泊水面半径，m；A为湖泊水面积，m²；R₂为湖底面半径，m；H_A为圆锥体高度，m；B为湖底面积，m²；H_B为湖底圆锥体的高，m. H_B和B为定值，则知道了A就可以求出H_A，反之亦然.

由于圆台体湖泊类型一般为图 2所示，故以湖泊水面和湖底所构成的圆锥体顶部为(0，0)，垂直于湖泊水面的方向为纵坐标，平行于湖泊水面的方向为横坐标建立直角坐标系(图 3)，则单位时间内地下水与湖泊的交换量为：

$ {W_{\rm{q}}}\left( t \right) = \frac{{{H_{\rm{A}}}\left( t \right) - H\left( t \right)}}{M}K \cdot A\left( t \right) $

(4)

式中，H为t年份湖泊周围地下水位平均高程，m；M为湖底淤泥层平均厚度，m；H_A为t年份的湖泊水位高程(即等于圆锥体高度)，m；K为湖泊淤泥层垂向渗透系数，m/d.

把公式(4)代入公式(1)得到如下方程：

$ \frac{{\partial V}}{{\partial t}} = A\left( t \right)P\left( t \right) + f\left( t \right) - A\left( t \right)E\left( t \right) - Q\left( t \right) - \frac{{{H_{\rm{A}}}\left( t \right) - H\left( t \right)}}{M}K \cdot A\left( t \right) $

(5)

把公式(2)和公式(3)代入公式(5)得如下方程：

$ \begin{array}{l} \frac{{A\left( {t + \Delta t} \right)}}{{\Delta t}} = \frac{{2\sqrt {A\left( t \right)B} \left[ {P\left( t \right) - E\left( t \right)} \right]}}{{{H_{\rm{B}}}}} + 2\frac{{\sqrt B }}{{{H_{\rm{B}}}\sqrt {A\left( t \right)} }}\left[ {f\left( t \right) - Q\left( t \right)} \right]\\ - \frac{{2K\left[ {A\left( t \right){H_{\rm{B}}} - \sqrt {A\left( t \right)B} H\left( t \right)} \right]}}{{{H_{\rm{B}}} \cdot M}} \end{array} $

(6)

公式(6)提供了一种利用地下水位变化求解湖泊面积的方法.

研究区的湖泊均属于封闭型湖泊，主要由降水和地下水补给，无地表径流出入，排泄方式以蒸发和渗漏为主.根据实地调查资料可知，湖泊周围为牧区，无工农业用水，湖泊均为咸水湖不适合牲畜饮用，因此公式(6)可简化为：

$ \frac{{A\left( {t + \Delta t} \right) - A\left( t \right)}}{{\Delta t}} = \frac{{2\sqrt {A\left( t \right)B} \left[ {P\left( t \right) - E\left( t \right)} \right]}}{{{H_{\rm{B}}}}} - K\frac{{2\left[ {A\left( t \right){H_{\rm{B}}} - \sqrt {A\left( t \right)B} H\left( t \right)} \right]}}{{{H_{\rm{B}}} \cdot M}} $

(7)

根据湖泊水位与地下水位高程的关系，将研究区的地下水湖泊系统分为地下水位高于湖泊水位，地下水补给湖泊；地下水位低于湖泊水位，高于湖泊底部淤泥质底板高程，湖泊补给地下水且补给量与地下水位有关；地下水位低于湖泊底部淤泥质隔水层底板高程，湖泊补给地下水但补给量与地下水位无关3种类型，对应的关系方程为：

$ \frac{{A\left( {t + \Delta t} \right) - A\left( t \right)}}{{\Delta t}} = \\\left\{ {\begin{array}{*{20}{c}} {\frac{{2\sqrt {A\left( t \right)B} \left( {P\left( t \right) - E\left( t \right)} \right)}}{{{H_{\rm{B}}}}} + K\frac{{2\left[ {A\left( t \right){H_{\rm{B}}} - \sqrt {A\left( t \right)B} H\left( t \right)} \right]}}{{{H_{\rm{B}}} \cdot M}}, H\left( t \right) > {H_{\rm{A}}}}\\ {\frac{{2\sqrt {A\left( t \right)B} \left( {P\left( t \right) - E\left( t \right)} \right)}}{{{H_{\rm{B}}}}} - K\frac{{2\left[ {A\left( t \right){H_{\rm{B}}} - \sqrt {A\left( t \right)B} H\left( t \right)} \right]}}{{{H_{\rm{B}}} \cdot M}}, {H_{\rm{B}}} < H\left( t \right) \le {H_{\rm{A}}}}\\ {\frac{{2\sqrt {A\left( t \right)B} \left( {P\left( t \right) - E\left( t \right)} \right)}}{{{H_{\rm{B}}}}} - \frac{{2\left[ {A\left( t \right){H_{\rm{B}}}K - \sqrt {A\left( t \right)B} KM} \right]}}{{{H_{\rm{B}}} \cdot M}}, H\left( t \right) < {H_{\rm{B}}}} \end{array}} \right.\;\; $

(8)

公式(8)基本涵盖了草原圆台型湖泊面积与地下水位的响应关系.只要知道湖泊底面积(B)、某一时刻的湖水面积(A)、水深等初始条件后，根据测得的蒸发、降雨、地下水等常观资料，就可以求出任何时刻的湖泊面积.

4 结果 4.1 矿区地下水流场变化

伊敏露天煤矿于1983年开始投产，目前年产原煤2300万t，预计2045年闭矿^[24].根据历史地下水观测孔实测的地下水位数据，采用GIS中的克里金插值法得到1982、2000、2012、2017年8月份的地下水等水位线(图 4).由图 4可知，煤矿开采前研究区内的地下水位在660~674 m之间，波动范围较小，除靠近丘陵的部分西南区域外，研究区地下水位仅相差5 m左右，湖泊周围地下水位最浅，在0~2 m之间. 2000年以后，研究区南部以采坑为中心的地下水位出现大幅下降，地下水位降落111 m，而北部地区的地下水位变化并不大，这与冯海波等的研究结果一致^[21].从2012年开始，采坑周围不再有潜水，相较2000年，研究区北部的地下水位基本无变化；2012- 2017年，研究区的地下水位趋于稳定，不再有大幅的波动.

这是因为煤矿开采前(即1982年)伊敏盆地地下水流场分布主要由地形地貌、地质构造等决定，地下水流向伊加诺尔、巴嘎诺尔湖群为研究区南部局部的地下水排泄区，区域地下水总体由南部和西部台地丘陵向伊敏河排泄.煤矿开采后，受露天采坑和疏干水的影响，地下水位的波动幅度较为剧烈，区域地下水流系统的补径排特征发生明显改变.

4.2 矿区湖泊面积变化

随着开采年限的延长，湖群数量及面积在不断发生变化(表 2、图 5).煤矿开发前(1982年)，研究区有骆驼脖子草库伦、柴达敏诺尔、哈尔呼吉诺尔、哈沙廷布拉格、伊加诺尔和巴嘎诺尔5处湖群，自2000年开始，研究区内的湖群只有骆驼脖子草库伦和柴达敏诺尔两处，其水域面积也在不断变小.截止到2017年，湖泊总面积由采矿前的6.94 km²缩小为1.12 km²，面积减少率达84 %.其中，目前最大的柴达敏诺尔湖泊面积仅存有0.97 km².伊敏煤矿开采35年来，研究区湖泊水面面积一直呈持续萎缩趋势.

本研究区湖泊属于封闭型内流湖，主要接受大气降水和地下水的补给，无地表径流出入，排泄量较少，因此煤矿开采前的湖泊数量多且面积较稳定.然而，露天矿区开发的影响下，研究区南部的伊加诺尔和巴嘎诺尔两个湖群已经被露天采坑所取代，并随着采矿规模和开采年限的不断增大，研究区以北的骆驼脖子草库伦和柴达敏诺尔湖群也出现逐年减少的趋势.

5 讨论 5.1 区域地下水-湖泊补排关系演化

人类活动(如抽水)、气象、地下水等因素均可能导致区域湖泊面积(水量)发生变化.研究区地处草原，区内湖泊均属于咸水湖，目前没有人工利用的途径.根据气象资料，1982-2017年间该区降水和蒸发量波动都比较平稳，而区内湖泊面积总体却呈显著下降趋势(图 6).因而，煤矿开采导致的地下水位变化可能是引起区内湖泊面积(水量)发生萎缩的最主要原因.

伊敏露天煤矿开采前，矿区地下水位高于湖泊水位，区内5个湖群均依赖地下水补给，属于地下水补给型湖泊(图 7a). 1982年以后，受矿山开发强烈抽排地下水的干扰，使地下水位持续下降，进而导致地下水对湖泊的补给逐渐减少，部分湖泊水体开始透过湖底淤泥质弱透水层补给地下水，水面面积不断萎缩.目前哈尔呼吉诺尔、哈沙廷布拉格等湖群正由地下水补给湖泊型向湖泊补给地下水型过渡.在这期间，由于地下水位的不断下降，且湖泊底部淤泥质弱透水层的阻滞作用，地下水补给湖泊越来越困难(图 7b).最终，地下水位下降至湖泊底部淤泥层以下，地下水不再补给湖泊，而湖泊却处于渗漏补给地下水状态，导致湖泊面积不断萎缩甚至消失(图 7c).由此可见，地下水是影响区内地下水-湖泊演变的关键因素.

5.2 矿区湖泊面积变化预测分析

为了进一步分析矿区开发对未来湖泊面积的影响，利用构建的地下水湖泊关系方程，结合实地调查及搜集的资料(表 3)预测了矿区湖泊总面积未来的演化趋势.本文选取表 3给出的相关数据检验模型的可靠性，由于在2000年以后，矿区只有以柴达敏诺尔湖泊骆驼脖子草库伦湖为首的两处湖群，因此，本文选取柴达敏诺尔(图 8a)和骆驼脖子草库伦湖泊的实测面积和模型计算面积进行检验(图 8b)，从图上看计算值与实测值基本吻合，其拟合优度均达到了0.8以上.并且随着地下水流场趋于稳定，地下水与湖泊之间的补排关系发生改变由地下水补给湖泊转变为湖泊补给地下水，因此模型后期不再考虑地下水这一变量，最终使得模型的拟合程度呈现出随时间变化不断增强的趋势.该模型能够应用于预测圆台型湖泊面积受地下水位的影响变化情况.另外，从图 4可以看出，2012年以后，地下水位的波动幅度很小，在假设未来地下水位稳定、气候因素波动不大、矿山开采正常运行的情况下，利用该模型可以预测出2045年(闭矿时)研究区湖泊总面积为0.56 km².

本研究中构建的为圆台型草原湖泊-地下水交互关系解析模型，在自然界中除了圆台型湖泊外，还存在诸如矩形、椭圆型及其他不规则几何体的湖泊类型.虽然它们形状及水文地质条件不同，但湖泊与地下水之间的相互作用机制是一样的.只要根据湖泊的几何形状进行合理的概化，并能够应用数学公式对进行描述，则根据水均衡原理、达西定律、几何学知识，湖泊面积和地下水位的关系模型就能够建立起来.通过调查掌握了某时刻的湖泊面积、湖泊水深、地下水位高程等初始条件，便可在预设气象、水文、人工开采等因素后，利用模型预测未来任意时刻的湖泊水面面积.该模型计算准确度除了与掌握的气象、水文、水文地质等资料的精度有关外，还与时间步长有关.理论上，在迭代计算过程中，Δt取值越小，则方程计算精度越高，在实际应用中，应根据需要选取合适的时间步长.该类地下水位湖泊面积解析模型，参数简单且容易获取，计算方便，实用范围较广.

6 结论

本文分析了伊敏煤矿开采35年来区域地下水系统演化特征，通过建立地下水位与湖泊面积关系模型定量揭示了煤矿开采引发的地下水位变化对湖泊面积的影响.主要结论如下：

1) 煤矿开采前研究区有5处湖泊群，开采至2000年时只存有柴达敏诺尔和骆驼脖子草库伦两处湖泊群，其他皆消失；截至2017年，湖泊总面积由采矿前的6.94 km²萎缩到1.12 km²，减少率达84 %.伊敏煤矿开采35年来，研究区湖泊数量及水体面积呈持续萎缩趋势.

2) 煤矿开采前，矿区地下水位高于湖泊水位，区内5个湖群均依赖地下水补给，属于地下水补给型湖泊；随着煤矿开采年限的增长及采矿规模的不断扩大，矿区地下水位下降至湖泊底部淤泥层以下地下水将不再补给湖泊，湖泊水面面积持续萎缩甚至消失；矿区地下水湖泊补给类型总体上由地下水补给湖泊型向湖泊补给地下水型演化.

3) 基于水均衡原理构建了草原湖泊地下水耦合模型，经检验拟合优度均达到了0.80以上，该模型能够应用于预测圆台型湖泊面积随地下水位变化情况.在气候因素波动不大，煤矿开发稳定的情况下，利用该模型预测发现，伊敏煤矿闭矿时(2045年)矿区湖泊总面积将萎缩至0.56 km².

致谢: 感谢中国地质大学(武汉)王全荣教授对文中构建湖泊面积与地下水位关系模型的指导.