夏蔓宏，董少剛**，劉白薇，李 銥，李政葵，王 超，周雨澤

(1：內(nèi)蒙古大學(xué)生態(tài)與環(huán)境學(xué)院，呼和浩特 010021) (2：內(nèi)蒙古大學(xué)社會(huì)科學(xué)處，呼和浩特 010021)

強(qiáng)烈人工干預(yù)下的地下水流動(dòng)系統(tǒng)變化，已經(jīng)成為影響局部乃至區(qū)域生態(tài)水文平衡的主要驅(qū)動(dòng)力[1]. 干旱-半干旱地區(qū)的煤炭資源開(kāi)發(fā)往往會(huì)改變區(qū)域地下水循環(huán)模式，破壞地下水與湖泊、河流等地表水體之間的水力聯(lián)系[2-5]，而湖泊對(duì)維持生態(tài)平衡、保護(hù)區(qū)域棲息物種等具有重要意義[6-11]. 國(guó)內(nèi)外學(xué)者們對(duì)于煤碳資源開(kāi)采對(duì)湖泊的影響進(jìn)行了大量研究. 趙寧等和楊立彬等[12-13]研究發(fā)現(xiàn)鄂爾多斯盆地的風(fēng)沙區(qū)依賴地下水補(bǔ)給的湖泊多達(dá)500個(gè)，由于建壩、打井抽取地下水和煤炭資源開(kāi)采等人為因素的影響，引發(fā)區(qū)域地下水位大幅下降、部分湖泊濕地變干，如中國(guó)第一大沙漠湖泊-紅堿淖的水域面積減少了26.66 km2；王文科等[14]研究發(fā)現(xiàn)荒漠化地帶的陜北礦區(qū)內(nèi)河流及湖泊40%～70%的流量來(lái)自地下水的補(bǔ)給，而在大規(guī)模煤礦山開(kāi)發(fā)過(guò)程中，地下水普遍下降10～12 m，導(dǎo)致區(qū)域河水流量及湖泊水面大幅減少；馬雄德等[15]研究發(fā)現(xiàn)，在1990-2011年期間，沙漠灘地的陜西榆神府礦區(qū)內(nèi)水體總體呈現(xiàn)減少的趨勢(shì)，境內(nèi)水體面積從162.63 km2降低到103.84 km2；P?schke等[16]對(duì)斯坦克林湖的研究中發(fā)現(xiàn)只有一部分集水區(qū)的地下水補(bǔ)給湖泊，而其余的則向更深的含水層泄漏；Rudnick等[17]研究發(fā)現(xiàn)湖泊的地下水排放可以在湖泊的水平衡中發(fā)揮主要作用；Hajati等[18]對(duì)漢普頓湖的研究中發(fā)現(xiàn)在低于平均降水量的情況下，該湖經(jīng)歷了一次倒流，在此期間該湖主要轉(zhuǎn)變?yōu)檠a(bǔ)給湖，然后于2012年初返回一個(gè)過(guò)流湖.

綜上所述，學(xué)者們對(duì)于湖泊與地下水之間的演化關(guān)系進(jìn)行了一定的研究，但定量[16,19-20]評(píng)價(jià)地下水對(duì)湖泊的補(bǔ)給作用是很困難的，現(xiàn)有研究主要集中于定性描述煤礦開(kāi)采引起的地下水位變化對(duì)湖泊水量的影響，有鑒于此，本研究以內(nèi)蒙古呼倫貝爾草原伊敏露天煤礦為研究對(duì)象，在對(duì)區(qū)域地下水流系統(tǒng)調(diào)查的基礎(chǔ)上，利用遙感技術(shù)、地下水流系統(tǒng)理論及水均衡原理，探討煤礦開(kāi)采活動(dòng)對(duì)地下水-湖泊系統(tǒng)的演變影響，構(gòu)建地下水位與湖泊面積之間的關(guān)系模型，定量揭示采礦引起的地下水位變化對(duì)湖泊面積(水量)的影響.

1 研究區(qū)概況

1.1 地理地貌

研究區(qū)位于大興安嶺西坡，伊敏盆地南部，伊敏河中下游地區(qū)，隸屬于內(nèi)蒙古自治區(qū)呼倫貝爾市鄂溫克自治旗，面積約100 km2. 伊敏煤田地貌呈一半封閉型的盆地，東西兩側(cè)為丘陵，南部為臺(tái)地，盆地內(nèi)為沖積平原，海拔在651～673 m，南高北低，東西兩側(cè)稍高，中心低洼[10](圖1).

1.2 氣象水文

研究區(qū)屬于典型的中溫帶大陸性季風(fēng)氣候，四季分明，冬冷夏涼. 年均氣溫-1.9℃，境內(nèi)降水自東南向西北遞減，多年年均降水量為354.73 mm，年均蒸發(fā)量為1318.3 mm.

區(qū)內(nèi)最主要的河流為伊敏河，葦子坑河等其余河流均為其支流，伊敏河由南向北從盆地東部穿過(guò)，全長(zhǎng)2000 km，流域面積可達(dá)9000 km2，河床寬約60 m，水深0.50～2.50 m. 目前伊敏露天煤礦區(qū)內(nèi)地下水補(bǔ)給主要依賴于大氣降水、伊敏河以及盆地內(nèi)零星分布的湖泊. 天然狀況下，伊敏河水與地下水的關(guān)系為豐水季節(jié)伊敏河水補(bǔ)給地下水，枯水季節(jié)地下水補(bǔ)給河水. 在煤礦開(kāi)采的影響下，伊敏河受影響河段與地下水的關(guān)系是伊敏河水補(bǔ)給地下水. 由于伊敏礦區(qū)疏排地下水，改變了區(qū)域地下水循環(huán)特征. 目前地下水排泄方式主要是向露天礦疏干區(qū)排泄. 地表水則以徑流的方式，匯集于盆地低洼處，或滲透于地層中或以蒸發(fā)的形式排泄[21].

圖1 研究區(qū)地理地貌Fig.1 Geographical and topographic map of the research area

2 數(shù)據(jù)來(lái)源與處理方法

遙感影像采用Landsat MSS、TM和 OLI 數(shù)據(jù) (表1)；水文地質(zhì)資料來(lái)源于黑龍江省伊敏煤田精查報(bào)告、伊敏煤礦水文地質(zhì)調(diào)查報(bào)告、伊敏煤礦礦區(qū)地下水位動(dòng)態(tài)監(jiān)測(cè)報(bào)告及課題組2016-2018年的實(shí)地調(diào)查；氣象數(shù)據(jù)來(lái)源于當(dāng)?shù)貧庀缶?

表 1 研究采用的 Landsat 影像數(shù)據(jù)集

水體指數(shù)法是實(shí)現(xiàn)水體快速提取的簡(jiǎn)單便捷方法之一[22]. 伊敏盆地的湖泊數(shù)量較多，大小不一，故采用水體指數(shù)法并結(jié)合目視解譯的方式進(jìn)行提取. 本文利用在 ENVI5.1+IDL平臺(tái)開(kāi)發(fā)出的歸一化水體指數(shù)(Normal difference water index，NDWI)從上述 4 期Landsat 系列遙感數(shù)據(jù)中提取伊敏盆地水體信息并繪制湖泊水體圖，然后結(jié)合實(shí)地調(diào)查和收集的數(shù)據(jù)，分析煤礦開(kāi)采35年來(lái)地下水-湖泊系統(tǒng)的演化規(guī)律.

3 區(qū)域地下水-湖泊關(guān)系模型

3.1 湖泊水均衡方程

湖水的補(bǔ)給包括降雨、可能的地表徑流及地下水補(bǔ)給，湖泊的排泄主要是蒸發(fā)、可能的地表徑流排泄以及工農(nóng)業(yè)用水. 因此，一般湖泊的水均衡方程為：

(1)

式中，V為湖泊的湖水體積，m3；?V/?t為單位時(shí)間內(nèi)的湖泊水體體積變化率，m3/a；A為t年份的湖泊水面面積，m2；P為t年份的年降水量，m；f(t)為產(chǎn)匯流函數(shù)；E為t年份的年蒸發(fā)量，m；Q為t年份的湖泊人工抽取水量及地表徑流排泄量，m3/a；Wq為單位時(shí)間內(nèi)地下水與湖泊的交換量，m3/a.

圖2 湖泊水體體積示意圖Fig.2 Diagram of lake water volume

3.2 草原圓臺(tái)型湖泊面積與地下水位關(guān)系方程

前人研究表明直接測(cè)量地下水與湖泊的交換量不僅耗時(shí)、耗力還容易出錯(cuò)[23]，本次研究利用水均衡原理、達(dá)西定律及幾何學(xué)知識(shí)建立地下水位-湖泊面積之間的關(guān)系方程. 由于研究區(qū)湖泊四周地勢(shì)平坦，形態(tài)大多近似圓臺(tái)體，故以圓臺(tái)體湖泊類型為例，湖泊水體體積示意見(jiàn)圖2.

湖泊水體體積可表示為：

(2)

則HA為：

(3)

式中，V為湖泊水體體積，m3；R1為湖泊水面半徑，m；A為湖泊水面積，m2；R2為湖底面半徑，m；HA為圓錐體高度，m；B為湖底面積，m2；HB為湖底圓錐體的高，m.HB和B為定值，則知道了A就可以求出HA，反之亦然.

由于圓臺(tái)體湖泊類型一般為圖2所示，故以湖泊水面和湖底所構(gòu)成的圓錐體頂部為(0，0)，垂直于湖泊水面的方向?yàn)榭v坐標(biāo)，平行于湖泊水面的方向?yàn)闄M坐標(biāo)建立直角坐標(biāo)系(圖3)，則單位時(shí)間內(nèi)地下水與湖泊的交換量為：

(4)

式中，H為t年份湖泊周圍地下水位平均高程，m；M為湖底淤泥層平均厚度，m；HA為t年份的湖泊水位高程(即等于圓錐體高度)，m；K為湖泊淤泥層垂向滲透系數(shù)，m/d.

圖3 地下水-湖泊交互關(guān)系概念示意圖Fig.3 Schematic diagram of groundwater-lake interaction

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

(5)

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

(6)

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

研究區(qū)的湖泊均屬于封閉型湖泊，主要由降水和地下水補(bǔ)給，無(wú)地表徑流出入，排泄方式以蒸發(fā)和滲漏為主. 根據(jù)實(shí)地調(diào)查資料可知，湖泊周圍為牧區(qū)，無(wú)工農(nóng)業(yè)用水，湖泊均為咸水湖不適合牲畜飲用，因此公式(6)可簡(jiǎn)化為：

(7)

根據(jù)湖泊水位與地下水位高程的關(guān)系，將研究區(qū)的地下水-湖泊系統(tǒng)分為地下水位高于湖泊水位，地下水補(bǔ)給湖泊；地下水位低于湖泊水位，高于湖泊底部淤泥質(zhì)底板高程，湖泊補(bǔ)給地下水且補(bǔ)給量與地下水位有關(guān)；地下水位低于湖泊底部淤泥質(zhì)隔水層底板高程，湖泊補(bǔ)給地下水但補(bǔ)給量與地下水位無(wú)關(guān)3種類型，對(duì)應(yīng)的關(guān)系方程為：

(8)

公式(8)基本涵蓋了草原圓臺(tái)型湖泊面積與地下水位的響應(yīng)關(guān)系. 只要知道湖泊底面積(B)、某一時(shí)刻的湖水面積(A)、水深等初始條件后，根據(jù)測(cè)得的蒸發(fā)、降雨、地下水等常觀資料，就可以求出任何時(shí)刻的湖泊面積.

4 結(jié)果

4.1 礦區(qū)地下水流場(chǎng)變化

伊敏露天煤礦于1983年開(kāi)始投產(chǎn)，目前年產(chǎn)原煤2300萬(wàn)t，預(yù)計(jì)2045年閉礦[24]. 根據(jù)歷史地下水觀測(cè)孔實(shí)測(cè)的地下水位數(shù)據(jù)，采用GIS中的克里金插值法得到1982、2000、2012、2017年8月份的地下水等水位線(圖4). 由圖4可知，煤礦開(kāi)采前研究區(qū)內(nèi)的地下水位在660～674 m之間，波動(dòng)范圍較小，除靠近丘陵的部分西南區(qū)域外，研究區(qū)地下水位僅相差5 m左右，湖泊周圍地下水位最淺，在0～2 m之間. 2000年以后，研究區(qū)南部以采坑為中心的地下水位出現(xiàn)大幅下降，地下水位降落111 m，而北部地區(qū)的地下水位變化并不大，這與馮海波等的研究結(jié)果一致[21]. 從2012年開(kāi)始，采坑周圍不再有潛水，相較2000年，研究區(qū)北部的地下水位基本無(wú)變化；2012-2017年，研究區(qū)的地下水位趨于穩(wěn)定，不再有大幅的波動(dòng).

這是因?yàn)槊旱V開(kāi)采前(即1982年)伊敏盆地地下水流場(chǎng)分布主要由地形地貌、地質(zhì)構(gòu)造等決定，地下水流向伊加諾爾、巴嘎諾爾湖群為研究區(qū)南部局部的地下水排泄區(qū)，區(qū)域地下水總體由南部和西部臺(tái)地丘陵向伊敏河排泄. 煤礦開(kāi)采后，受露天采坑和疏干水的影響，地下水位的波動(dòng)幅度較為劇烈，區(qū)域地下水流系統(tǒng)的補(bǔ)徑排特征發(fā)生明顯改變.

圖4 1982-2017年伊敏礦區(qū)地下水等水位線變化Fig.4 Changes of groundwater levels in the Yimin mining area from 1982 to 2017

4.2 礦區(qū)湖泊面積變化

隨著開(kāi)采年限的延長(zhǎng)，湖群數(shù)量及面積在不斷發(fā)生變化(表2、圖5). 煤礦開(kāi)發(fā)前(1982年)，研究區(qū)有駱駝脖子草庫(kù)倫、柴達(dá)敏諾爾、哈爾呼吉諾爾、哈沙廷布拉格、伊加諾爾和巴嘎諾爾5處湖群，自2000年開(kāi)始，研究區(qū)內(nèi)的湖群只有駱駝脖子草庫(kù)倫和柴達(dá)敏諾爾兩處，其水域面積也在不斷變小. 截止到2017年，湖泊總面積由采礦前的6.94 km2縮小為1.12 km2，面積減少率達(dá)84%. 其中，目前最大的柴達(dá)敏諾爾湖泊面積僅存有0.97 km2. 伊敏煤礦開(kāi)采35年來(lái)，研究區(qū)湖泊水面面積一直呈持續(xù)萎縮趨勢(shì).

本研究區(qū)湖泊屬于封閉型內(nèi)流湖，主要接受大氣降水和地下水的補(bǔ)給，無(wú)地表徑流出入，排泄量較少，因此煤礦開(kāi)采前的湖泊數(shù)量多且面積較穩(wěn)定. 然而，露天礦區(qū)開(kāi)發(fā)的影響下，研究區(qū)南部的伊加諾爾和巴嘎諾爾兩個(gè)湖群已經(jīng)被露天采坑所取代，并隨著采礦規(guī)模和開(kāi)采年限的不斷增大，研究區(qū)以北的駱駝脖子草庫(kù)倫和柴達(dá)敏諾爾湖群也出現(xiàn)逐年減少的趨勢(shì).

表2 1982-2017年伊敏礦區(qū)湖泊面積

圖5 1982-2017年伊敏礦區(qū)湖泊示意圖Fig.5 Schematic diagrams of the lakes in Yimin mining area from 1982 to 2017

圖6 氣象因素(降水量和蒸發(fā)量)與湖泊面積變化Fig.6 Changes of meteorological factors and lake area

5 討論

5.1 區(qū)域地下水-湖泊補(bǔ)排關(guān)系演化

人類活動(dòng)(如抽水)、氣象、地下水等因素均可能導(dǎo)致區(qū)域湖泊面積(水量)發(fā)生變化. 研究區(qū)地處草原，區(qū)內(nèi)湖泊均屬于咸水湖，目前沒(méi)有人工利用的途徑. 根據(jù)氣象資料，1982-2017年間該區(qū)降水和蒸發(fā)量波動(dòng)都比較平穩(wěn)，而區(qū)內(nèi)湖泊面積總體卻呈顯著下降趨勢(shì)(圖6). 因而，煤礦開(kāi)采導(dǎo)致的地下水位變化可能是引起區(qū)內(nèi)湖泊面積(水量)發(fā)生萎縮的最主要原因.

圖7 礦區(qū)地下水-湖泊補(bǔ)給類型演化Fig.7 Evolution of groundwater-lake recharge type in mining area

伊敏露天煤礦開(kāi)采前，礦區(qū)地下水位高于湖泊水位，區(qū)內(nèi)5個(gè)湖群均依賴地下水補(bǔ)給，屬于地下水補(bǔ)給型湖泊(圖7a). 1982年以后，受礦山開(kāi)發(fā)強(qiáng)烈抽排地下水的干擾，使地下水位持續(xù)下降，進(jìn)而導(dǎo)致地下水對(duì)湖泊的補(bǔ)給逐漸減少，部分湖泊水體開(kāi)始透過(guò)湖底淤泥質(zhì)弱透水層補(bǔ)給地下水，水面面積不斷萎縮. 目前哈爾呼吉諾爾、哈沙廷布拉格等湖群正由地下水補(bǔ)給湖泊型向湖泊補(bǔ)給地下水型過(guò)渡. 在這期間，由于地下水位的不斷下降，且湖泊底部淤泥質(zhì)弱透水層的阻滯作用，地下水補(bǔ)給湖泊越來(lái)越困難(圖7b). 最終，地下水位下降至湖泊底部淤泥層以下，地下水不再補(bǔ)給湖泊，而湖泊卻處于滲漏補(bǔ)給地下水狀態(tài)，導(dǎo)致湖泊面積不斷萎縮甚至消失 (圖7c). 由此可見(jiàn)，地下水是影響區(qū)內(nèi)地下水-湖泊演變的關(guān)鍵因素.

5.2 礦區(qū)湖泊面積變化預(yù)測(cè)分析

為了進(jìn)一步分析礦區(qū)開(kāi)發(fā)對(duì)未來(lái)湖泊面積的影響，利用構(gòu)建的地下水-湖泊關(guān)系方程，結(jié)合實(shí)地調(diào)查及搜集的資料(表3)預(yù)測(cè)了礦區(qū)湖泊總面積未來(lái)的演化趨勢(shì). 本文選取表3給出的相關(guān)數(shù)據(jù)檢驗(yàn)?zāi)Ｐ偷目煽啃?，由于?000年以后，礦區(qū)只有以柴達(dá)敏諾爾湖泊駱駝脖子草庫(kù)倫湖為首的兩處湖群，因此，本文選取柴達(dá)敏諾爾(圖8a)和駱駝脖子草庫(kù)倫湖泊的實(shí)測(cè)面積和模型計(jì)算面積進(jìn)行檢驗(yàn)(圖8b)，從圖上看計(jì)算值與實(shí)測(cè)值基本吻合，其擬合優(yōu)度均達(dá)到了0.8以上. 并且隨著地下水流場(chǎng)趨于穩(wěn)定，地下水與湖泊之間的補(bǔ)排關(guān)系發(fā)生改變由地下水補(bǔ)給湖泊轉(zhuǎn)變?yōu)楹囱a(bǔ)給地下水，因此模型后期不再考慮地下水這一變量，最終使得模型的擬合程度呈現(xiàn)出隨時(shí)間變化不斷增強(qiáng)的趨勢(shì). 該模型能夠應(yīng)用于預(yù)測(cè)圓臺(tái)型湖泊面積受地下水位的影響變化情況. 另外，從圖4可以看出，2012年以后，地下水位的波動(dòng)幅度很小，在假設(shè)未來(lái)地下水位穩(wěn)定、氣候因素波動(dòng)不大、礦山開(kāi)采正常運(yùn)行的情況下，利用該模型可以預(yù)測(cè)出2045年(閉礦時(shí))研究區(qū)湖泊總面積為0.56 km2.

表3 地下水-湖泊解析模型參數(shù)取值

圖8 1982-2017年間柴達(dá)敏諾爾湖泊(a)和駱駝脖子草庫(kù)倫湖泊(b)面積變化Fig.8 Changes in area of Lake Chedamenor (a) and Lake Luotuobozicaoykulun (b) from 1982 to 2017

本研究中構(gòu)建的為圓臺(tái)型草原湖泊-地下水交互關(guān)系解析模型，在自然界中除了圓臺(tái)型湖泊外，還存在諸如矩形、橢圓型及其他不規(guī)則幾何體的湖泊類型. 雖然它們形狀及水文地質(zhì)條件不同，但湖泊與地下水之間的相互作用機(jī)制是一樣的. 只要根據(jù)湖泊的幾何形狀進(jìn)行合理的概化，并能夠應(yīng)用數(shù)學(xué)公式對(duì)進(jìn)行描述，則根據(jù)水均衡原理、達(dá)西定律、幾何學(xué)知識(shí)，湖泊面積和地下水位的關(guān)系模型就能夠建立起來(lái). 通過(guò)調(diào)查掌握了某時(shí)刻的湖泊面積、湖泊水深、地下水位高程等初始條件，便可在預(yù)設(shè)氣象、水文、人工開(kāi)采等因素后，利用模型預(yù)測(cè)未來(lái)任意時(shí)刻的湖泊水面面積. 該模型計(jì)算準(zhǔn)確度除了與掌握的氣象、水文、水文地質(zhì)等資料的精度有關(guān)外，還與時(shí)間步長(zhǎng)有關(guān). 理論上，在迭代計(jì)算過(guò)程中，Δt取值越小，則方程計(jì)算精度越高，在實(shí)際應(yīng)用中，應(yīng)根據(jù)需要選取合適的時(shí)間步長(zhǎng). 該類地下水位-湖泊面積解析模型，參數(shù)簡(jiǎn)單且容易獲取，計(jì)算方便，實(shí)用范圍較廣.

6 結(jié)論

本文分析了伊敏煤礦開(kāi)采35年來(lái)區(qū)域地下水系統(tǒng)演化特征， 通過(guò)建立地下水位與湖泊面積關(guān)系模型定量揭示了煤礦開(kāi)采引發(fā)的地下水位變化對(duì)湖泊面積的影響. 主要結(jié)論如下：

1) 煤礦開(kāi)采前研究區(qū)有5處湖泊群，開(kāi)采至2000年時(shí)只存有柴達(dá)敏諾爾和駱駝脖子草庫(kù)倫兩處湖泊群，其他皆消失；截至2017年，湖泊總面積由采礦前的6.94 km2萎縮到1.12 km2，減少率達(dá)84%. 伊敏煤礦開(kāi)采35年來(lái)，研究區(qū)湖泊數(shù)量及水體面積呈持續(xù)萎縮趨勢(shì).

2) 煤礦開(kāi)采前，礦區(qū)地下水位高于湖泊水位，區(qū)內(nèi)5個(gè)湖群均依賴地下水補(bǔ)給，屬于地下水補(bǔ)給型湖泊；隨著煤礦開(kāi)采年限的增長(zhǎng)及采礦規(guī)模的不斷擴(kuò)大，礦區(qū)地下水位下降至湖泊底部淤泥層以下地下水將不再補(bǔ)給湖泊，湖泊水面面積持續(xù)萎縮甚至消失；礦區(qū)地下水-湖泊補(bǔ)給類型總體上由地下水補(bǔ)給湖泊型向湖泊補(bǔ)給地下水型演化.

3) 基于水均衡原理構(gòu)建了草原湖泊-地下水耦合模型，經(jīng)檢驗(yàn)擬合優(yōu)度均達(dá)到了0.80以上，該模型能夠應(yīng)用于預(yù)測(cè)圓臺(tái)型湖泊面積隨地下水位變化情況. 在氣候因素波動(dòng)不大，煤礦開(kāi)發(fā)穩(wěn)定的情況下，利用該模型預(yù)測(cè)發(fā)現(xiàn)，伊敏煤礦閉礦時(shí)(2045年)礦區(qū)湖泊總面積將萎縮至0.56 km2.

致謝：感謝中國(guó)地質(zhì)大學(xué)(武漢)王全榮教授對(duì)文中構(gòu)建湖泊面積與地下水位關(guān)系模型的指導(dǎo).