潘國營， 南立超， 郭二濤， 林 云

(河南理工大學(xué) 資源環(huán)境學(xué)院， 河南 焦作 454000)

焦作礦區(qū)巖溶水系統(tǒng)是我國北方巖溶區(qū)發(fā)育比較典型的系統(tǒng)之一。區(qū)內(nèi)奧陶系灰?guī)r厚度400 m，廣泛出露于礦區(qū)北部的太行山區(qū)，地表及地下巖溶非常發(fā)育，降水和河水入滲補(bǔ)給條件好。焦作山前地區(qū)是巖溶水集中排泄區(qū)，有九里山泉群和三姑泉群，九里山泉群歷史上最大流量達(dá)到12 m3/s，三姑泉目前流量在4 m3/s左右[1]。

焦作豐富的巖溶水資源為焦作市工業(yè)和居民生活用水提供了寶貴的水源。圍繞焦作礦區(qū)巖溶水的開發(fā)利用，河南省地礦局第一水文地質(zhì)隊、中國地質(zhì)科學(xué)院巖溶地質(zhì)研究所等單位開展了大量勘探和研究工作，查明了水文地質(zhì)條件，確定了巖溶水資源量為7.732 m3/s，據(jù)此，焦作市相繼建設(shè)了崗莊、東小莊、周莊等幾個大型開采巖溶水的水廠[2]。礦井排水和城市開采是巖溶水的兩種排泄方式，20世紀(jì)80年代，二者之和超過10 m3/s，超過巖溶水多年平均補(bǔ)給量。巖溶水位呈階梯狀下降，引發(fā)了九里山泉斷流及水質(zhì)惡化等環(huán)境水文地質(zhì)問題，并引起人們的高度關(guān)注，例如，劉白宙[3]研究了山區(qū)降雨對地下水位的影響，黃平華等[4]建立了降雨與水位間的統(tǒng)計模型，潘國營等[5]分析了焦作礦區(qū)巖溶水的氯化物污染原因。歷經(jīng)多年的采煤，焦作礦區(qū)的煤炭資源接近枯竭，一些煤礦相繼關(guān)閉，礦井排水量將大幅度減少。與此同時，焦作電廠、化工廠等高耗水企業(yè)或遷至外地或?qū)崿F(xiàn)了轉(zhuǎn)型，且開采巖溶水的水廠因有了南水北調(diào)工程的水源，將逐漸停止開采巖溶水，城市開采巖溶水量和礦井排水量都將繼續(xù)減少。在閉礦和城市減采雙向因素影響下，巖溶水位如何變化以及九里山一帶能否出現(xiàn)歷史上“百泉噴涌、湖水碧綠”的美景，是焦作市打造山水園林城市迫切想要知道的問題。

因此，本文基于50年來焦作礦區(qū)巖溶水水位、降水量、開采量及排水量等長時間序列數(shù)據(jù)，利用相關(guān)分析和關(guān)聯(lián)度分析等方法，來分析降水、人工開采地下水對巖溶水位的影響，并預(yù)測未來巖溶水位的變化。

1 焦作礦區(qū)水文地質(zhì)背景

焦作礦區(qū)奧陶系灰?guī)r厚400 m，是區(qū)域性的強(qiáng)富水層，巖溶裂隙發(fā)育。在北部山區(qū)呈裸露型，山前傾斜平原區(qū)則掩埋于石炭-二疊和新生界地層之下。奧陶系灰?guī)r在太行山區(qū)露頭面積近2 000 km2，地表及地下巖溶發(fā)育，山區(qū)大氣降水豐富，大氣降水入滲是焦作巖溶水重要的補(bǔ)給來源之一[6]。丹河常年有水，流經(jīng)灰?guī)r分布區(qū)，河床滲漏嚴(yán)重，多年平均滲漏量為1.60 m3/s。西石河、山門河和子房溝河屬季節(jié)性河流，經(jīng)碳酸鹽巖分布區(qū)，河水在距出山口5～10 km地段全部漏失補(bǔ)給地下水。地表水沿河滲漏也是焦作巖溶水的重要補(bǔ)給來源之一[7]。巖溶水在焦作北部和西部接受補(bǔ)給后，由北向南、自西向東南流向太行山前的焦作礦區(qū)。焦作礦區(qū)發(fā)育了十幾條縱橫交錯的隱伏斷層，來自山區(qū)的巖溶水沿鳳凰嶺斷層、九里山斷層、朱村斷層等強(qiáng)導(dǎo)水?dāng)嗔堰\動、富集，并形成巖溶水強(qiáng)徑流帶。

天然條件下，巖溶水在九里山奧陶系灰?guī)r殘丘一帶的“天窗”處以泉群形式集中排泄，泉水最高流量達(dá)到12 m3/s，目前，人工開采和礦井排水是巖溶水主要的排泄方式，其中城市開采量為2.6 m3/s，礦井排水量為5.2 m3/s[8-10]。

2 巖溶水位動態(tài)變化特征

巖溶水位動態(tài)變化主要受山區(qū)大氣降水和人工開采(包括礦井排水)雙重因素的影響。20世紀(jì)60年代，巖溶水開采量很小，巖溶水補(bǔ)給與排泄處于天然動態(tài)均衡狀態(tài)，水位在110 m(海拔高程，下同)上下波動，高出九里山泉群排泄極限標(biāo)高(95 m)，巖溶水在九里山奧灰露頭周圍以泉群形式排泄。隨著礦井排水量的逐步增加，巖溶水位逐步下降，至20世紀(jì)80年代，礦井排水量達(dá)到峰值，與此同時，隨著崗莊水源地的建設(shè)，城市進(jìn)入了大規(guī)模開采利用巖溶水的時代，礦井排水和城市開采巖溶水水量超過10 m3/s，至20世紀(jì)80年代末巖溶水位降至75 m左右。自1994年開始，王封礦、焦西礦和焦東礦相繼被關(guān)閉，礦井排水量逐漸下降，但隨著城市發(fā)展和人口的增加，巖溶水開采量卻逐年上升，巖溶水開采總量仍然保持在較高的水平，在8.5～9.5 m3/s之間，巖溶水水位在75～90 m之間波動。2008年以來，焦作礦區(qū)各礦普遍采用了突水點注漿和采煤工作面底板預(yù)注漿等防治水措施，礦井排水量明顯下降，巖溶水水位有所回升，尤其是2012年焦作電廠異地重建后，其崗莊水源地開采量大幅度下降，水位回升到90 m。本文以5 a為滑動長度對焦作巖溶水年平均水位進(jìn)行計算，得出5 a平均水位滑動平均值做出統(tǒng)計圖，見圖1。

圖1 巖溶水位5 a滑動平均統(tǒng)計圖

經(jīng)過滑動平均后，序列中短于滑動長度的周期大大削弱，顯現(xiàn)出變化趨勢。由圖1可以看出,從1965-1993年焦作巖溶水位呈逐漸下降趨勢；1994-2007年巖溶水位趨于穩(wěn)定并有些回升的趨勢；2008-2015年巖溶水位呈上升趨勢。因此，本文將焦作市巖溶水位變化分為1965-1993年,1994-2007年和2008-2015年3個階段。

3 巖溶水位影響因素關(guān)聯(lián)性分析

在開發(fā)利用地下水資源和地下水預(yù)測工作中，一般都要對該地區(qū)地下水的動態(tài)變化進(jìn)行分析研究，尋找動態(tài)變化的影響因素以及其各因素的主次關(guān)系，以提高評估和預(yù)測工作的質(zhì)量[11]。

為了定量分析降水和地下水開采量對地下水動態(tài)變化的影響，采用斜率灰色關(guān)聯(lián)度法計算二者的影響程度。斜率灰色關(guān)聯(lián)度分析法，是基于系統(tǒng)各因素時間序列曲線間相似相異程度來衡量其關(guān)聯(lián)度大小的量化方法。對于離散數(shù)據(jù)列，若兩曲線在各時段上曲線斜率相等或相差較小，則二者的關(guān)聯(lián)系數(shù)就大；反之，關(guān)聯(lián)系數(shù)就越小[12-13]。

設(shè)xi為系統(tǒng)因素，在時間點上的觀測數(shù)據(jù)為xi(k)(k=1，2，…，n)，則稱：

xi=(xi(1),xi(2),…,xi(n))

(1)

(i=0,1,2,…,m)

為系統(tǒng)因素xi的時間序列。定義x0為母序列；xi為子序列，i=1，2，…，m。為了方便關(guān)聯(lián)度的計算，對各序列進(jìn)行無量綱化處理，得到新的數(shù)列：

yi=xid=(yi(1),yi(2),…,yi(n))

(2)

(i=0,1,2,…,m)

式中：d為序列算子。對新的序列求一階差商，找到曲線在各時點的斜率：

Δyi(k)=yi(k+1)-yi(k)

(3)

(k=1,2,…,n-1;i=0,1,…,m)

求灰色關(guān)聯(lián)系數(shù)：

(4)

(k=1,2,…,n-1;i=0,1,…,m)

綜合各點的關(guān)聯(lián)系數(shù)，得整個xi曲線與參考曲線x0的關(guān)聯(lián)度：

(5)

(i=1,2,…,m)

根據(jù)以上所述，以焦作市1965-2015年降水量、地下水開采量以及地下水位資料為依據(jù)，計算得出焦作地區(qū)降水量與巖溶水位關(guān)聯(lián)度為0.78，地下水開采量與巖溶水位關(guān)聯(lián)度為0.94。計算結(jié)果表明，降水量與地下水開采量都對焦作市地下水動態(tài)變化影響較大。

4 巖溶水動態(tài)影響因素的分析

由于降水入滲是焦作礦區(qū)巖溶水的重要補(bǔ)給來源，因此人們將水位下降的原因歸結(jié)為降水量的減少。為了分析降水量對巖溶水位的影響，根據(jù)焦作1955-2015年年均水位、年均降水量和人工開采量(包括城市開采和礦井排水)數(shù)據(jù)，運用Mann-Kendall非參數(shù)統(tǒng)計法進(jìn)行研究。

4.1 降水量對巖溶水位的影響

圖2為1965-2015年降水量和地下水位的統(tǒng)計圖。從圖2中可以看出，各序列變化都有一定的波動性，具有相似的“峰”、“谷”特征，年降水量出現(xiàn)豐枯變化時，地下水位也隨之產(chǎn)生高低變化，巖溶水位動態(tài)變化與降水量關(guān)系密切，但巖溶水位的多年波動下降是否是降水量減少所致，需要進(jìn)一步研究[14]。因此采用Mann-Kendall非參數(shù)檢驗方法[15]來檢測焦作巖溶水位與降水量長期變化趨勢。

圖2 1965-2015年降水量和地下水位統(tǒng)計圖

在雙邊趨勢檢驗中，對于給定的置信水平α，若對統(tǒng)計量Z，有∣Z∣≥Z1-α/2，則時間序列數(shù)據(jù)存在明顯的上升或下降趨勢。Z為正值表示增加趨勢，負(fù)值表示減少趨勢。Z的絕對值在大于等于1.28、1.64、2.32時表示分別通過了信度90%、95%、99%顯著性檢驗。由1965-2015年焦作巖溶水位和降水量資料計算得出Z值見表1。

表1 不同時段巖溶水位與降水量Z值

由計算結(jié)果可知巖溶水位在1965-1993年是下降趨勢，并通過了99%的顯著性檢驗；1994-2007年，巖溶水位有上升趨勢，但不顯著；2008-2015年巖溶水位有上升趨勢，并較為顯著，通過了95%的顯著性檢驗。降水量雖在1965-1993年有增大趨勢，但很不顯著，并且在1994-2007年和2008-2015年兩個階段Z值均為0，均沒有增大或減小的趨勢。在多年降水量沒有增大或減小趨勢的情況下，年均巖溶水水位下降或上升的趨勢顯著，由此可知焦作巖溶水位的多年波動下降并非降雨入滲補(bǔ)給量減少所致。

4.2 人工開采對巖溶水位的影響

焦作市巖溶水開采主要是城市巖溶水開采和礦井排水。焦作市1965年以前地下水開采量小于1.5 m3/s，水位在100～110 m間波動，最高達(dá)119 m，水位不呈下降趨勢；1965-1993年，隨著煤炭工業(yè)的大力發(fā)展，焦作電廠崗莊水源地建成使用，地下水開采總量升至10.46 m3/s，其中礦井排水量占87%，地下水位持續(xù)下降，最低降至70 m。1994-2008年，王封礦、焦西礦和焦東礦相繼關(guān)閉停產(chǎn)，整個礦區(qū)排水總量開始呈逐年遞減趨勢，但城市建設(shè)和發(fā)展速度加快，城市開采地下水量大幅度增加，地下水總開采量仍然保持在8.5～9.5 m3/s，地下水位在70～90 m之間波動。2008年以后，焦作礦區(qū)各礦普遍采用了注漿堵水技術(shù)，再加上部分礦井的關(guān)停停產(chǎn)，礦井排水量降至4.5 m3/s以下，并逐年減少，由于農(nóng)業(yè)節(jié)水的工程措施和城市工業(yè)用水節(jié)水的發(fā)展，城市地下水開采量也進(jìn)一步下降，總排水量穩(wěn)定在6 m3/s，地下水位有所回升，尤其是2012年焦作電廠異地重建后，其崗莊水源地開采量大幅度下降，開采量與地下水位統(tǒng)計圖見圖3。

圖3 開采量與地下水位統(tǒng)計圖

運用Mann-Kendall非參數(shù)檢驗方法對焦作巖溶水位與巖溶水開采量長期變化趨勢進(jìn)行檢測。計算結(jié)果如下：

由計算結(jié)果(表2)可以看出巖溶水開采量的變化趨勢。在1965-1993年巖溶水開采量增大趨勢顯著，并通過了99%的顯著性檢驗； 1994-2007年與2008-2015年兩個階段巖溶水開采量均為顯著下降趨勢。開采量呈上升趨勢時巖溶水位呈下降趨勢，開采量呈減少趨勢時巖溶水位呈上升趨勢。由此知焦作市巖溶水位的下降主要是由巖溶水開采量的增加造成的。

表2 巖溶水位與開采量Z值

5 巖溶水位預(yù)測

5.1 預(yù)測巖溶水位的多元回歸模型

前1年集中反映了前諸年的降水量、開采量對巖溶水位的綜合影響，故本年以前影響地下水位的各項因素可用前1年水位這一綜合因素來表示，但由于降水對巖溶水位的影響具有滯后性，分析巖溶水位與前1年降水量相關(guān)關(guān)系，其相關(guān)性比當(dāng)年更密切，故前1年降水量也是一項重要指標(biāo)。同時，當(dāng)年的降水量、開采量影響巖溶水位的因素對回歸計算的影響也有了綜合的體現(xiàn)[16]。

綜上分析可知，本年巖溶水位ht(m)受到前1年降水量Pt-1(mm)、本年降水量Pt(mm)和本年開采量Qt(m3/s)的影響，模型結(jié)果為：

ht=31.194+0.661ht-1+0.009Pt-1+

0.002Pt-0.951Qt

(6)

復(fù)相關(guān)系數(shù)R=0.901，決定系數(shù)為R2=0.812，模型有較好的線性關(guān)系，且回歸結(jié)果較好。如圖4所示。

圖4 巖溶水位多元回歸模型擬合圖

利用1965-2015年歷史水位、降水量和開采量數(shù)據(jù)對上述水位預(yù)測模型進(jìn)行了誤差檢驗，檢驗表明，模型顯著水平符合統(tǒng)計學(xué)要求，計算值與實際值的平均誤差為2.64%，擬合值與真實值差異較小，模型檢驗較好。

5.2 巖溶水位預(yù)測

預(yù)報未來礦井關(guān)停后的水位變化，設(shè)計了如下預(yù)報條件：

(1)根據(jù)1565-2016年焦作市降水量數(shù)據(jù)，采用皮爾遜模型計算出了不同降水保證率對應(yīng)的降水量，保證率95%(枯水年)、降水保證率75%(平水年)、降水保證率50%(豐水年)所對應(yīng)的降水量分別為378、532和612 mm。

(2)對于城市開采巖溶水量，考慮到焦作電廠已經(jīng)搬遷，向電廠供水的崗莊水源地開采巖溶水量已經(jīng)削減，焦作市利用南水北調(diào)工程的供水廠尚在建設(shè)，且用于滿足新建高新技術(shù)開發(fā)區(qū)的供水，現(xiàn)有開采巖溶水的水源地將維持現(xiàn)狀開采強(qiáng)度，因此，未來兩年城市開采量仍按2016年取值。

(3)2016年朱村礦和韓王礦已經(jīng)關(guān)閉，整個礦區(qū)剩余礦井排水量為3.09 m3/s。2017年計劃關(guān)停演馬礦，減少礦井排水0.80 m3/s，2018年將關(guān)閉方莊礦，再減少礦井排水0.35 m3/s，假如其他礦井的排水量仍維持2016年的水平不變，2017和2018年礦區(qū)總排水量將降低至2.29 m3/s和1.94 m3/s。

預(yù)測結(jié)果為見表3。

表3 2017、2018年巖溶水位預(yù)測值

預(yù)測結(jié)果顯示，2017年演馬礦停止排水后，礦區(qū)巖溶水位將繼續(xù)上升，枯水年水位為92.6 m，平水年為94.4 m，豐水年為95.5 m，2018年再關(guān)停方莊礦后，枯水年、平水年和豐水年對應(yīng)的巖溶水水位分別為92.8、95.6和97.2 m，即在降水量較為充沛的平水年和豐水年，焦作礦區(qū)巖溶水水位將超過九里山泉群最低排泄標(biāo)高，九里山泉群將恢復(fù)出流。

6 結(jié) 論

(1)1965年以來，焦作礦區(qū)年際降水量沒有明顯的增大或減小的趨勢，該區(qū)域巖溶水位的多年波動下降并非降雨入滲補(bǔ)給量減少所致。

(2)巖溶水人工開采量對焦作市巖溶水位動態(tài)變化影響很顯著，關(guān)聯(lián)度為0.94，呈負(fù)相關(guān)關(guān)系，是焦作溶水位多年波動下降的主要影響因素。

(3)繼2016年底關(guān)閉朱村礦與馮營礦后，2018年再關(guān)閉演馬礦和方莊礦，即開采總量降至4.05 m3/s時，巖溶水位將明顯上升，上升幅度在5～7 m，甚至將超過九里山泉群排泄標(biāo)高(95 m)，泉群將會復(fù)涌，這將為焦作市九里山白鷺濕地公園的恢復(fù)及建設(shè)提供水源條件。

參考文獻(xiàn)：

[1] 武 強(qiáng),潘國營,管恩太,等.焦作礦區(qū)突水災(zāi)害研究綜述[J].中國地質(zhì)災(zāi)害與防治學(xué)報，1995,6(4)：44-49.

[2] 管恩太,武 強(qiáng),李 鐸.焦作市地下水資源保護(hù)與利用研究[J].煤田地質(zhì)與勘探，2005,33(1):38-40.

[3] 劉白宙.焦作礦區(qū)地下水位動態(tài)變化特征研究[J].工程勘察，2007(11)：40-43.

[4] 黃平華,白萬備,鄧 勇.焦作煤礦區(qū)巖溶水水位統(tǒng)計模型研究[J].中國巖溶，2013,32(3)：299-304.

[5] 潘國營,賈軍輝.焦作礦區(qū)巖溶水Cl-污染原因初探[J].水文地質(zhì)工程地質(zhì)，2002，29(5)：50-51+64.

[6] 夏大平.焦作礦區(qū)礦井水資源化研究[D].焦作:河南理工大學(xué),2009.

[7] 潘國營,秦永泰,馬亞芬,等.基于R/S和Morlet小波分析的丹河徑流變化特征研究[J].水資源與水工程學(xué)報，2015,26(3)：41-45+50.

[8] 侯玉松,馬振民,雒蕓蕓,等.焦作地區(qū)水文地質(zhì)條件對淺層地下水污染的控制作用研究[J].中國農(nóng)村水利水電，2013(4)：40-44.

[9] 黃平華,陳建生,寧 超.焦作礦區(qū)地下水中氫氧同位素分析[J].煤炭學(xué)報，2012,37(5)：770-775.

[10] 潘國營,武 強(qiáng),董東林,等.焦作礦區(qū)巖溶裂隙網(wǎng)絡(luò)滲流特征及研究方法[J].煤炭學(xué)報，1998,23(6)：566-570.

[11] 馬而超.地下水位動態(tài)分析及預(yù)報研究[D].西安:長安大學(xué),2009.

[12] 周秀文.灰色關(guān)聯(lián)度的研究與應(yīng)用[D].長春:吉林大學(xué),2007.

[13] 李明涼.灰色關(guān)聯(lián)度新判別準(zhǔn)則及其計算公式[J].系統(tǒng)工程，1998,16(1)：68-70+61.

[14] 楊耀棟,李曉華,王蘭化,等.天津平原區(qū)地下水位動態(tài)特征與影響因素分析[J].地質(zhì)調(diào)查與研究，2011,34(4)：313-320.

[15] 劉 娟,陳濤濤,遲道才.基于Daniel及Mann-kendall檢驗的遼西北地區(qū)降雨量趨勢分析[J].沈陽農(nóng)業(yè)大學(xué)學(xué)報,2014,45(5)：599-603

[16] 束龍倉.地下水動態(tài)預(yù)測方法及其應(yīng)用[M].北京:中國水利水電出版社, 2010.