曾一凡，梅傲霜，武 強(qiáng)，華照來(lái)，趙 頔，杜 鑫，王 路，呂 揚(yáng)，潘 旭

(1.中國(guó)礦業(yè)大學(xué)(北京) 國(guó)家煤礦水害防治工程技術(shù)研究中心，北京 100083；2.陜西陜煤曹家灘礦業(yè)有限公司，陜西 榆林 719000；3. 北礦大(南京)新能源環(huán)保技術(shù)研究院，江蘇 南京 210005)

煤炭是我國(guó)重要的基礎(chǔ)能源，受制于我國(guó)缺油、少氣、相對(duì)富煤的能源稟賦影響，在今后相當(dāng)長(zhǎng)的一段時(shí)間內(nèi)，以煤炭為主體的能源格局將在我國(guó)長(zhǎng)期存在[1-2]。隨著煤炭資源開(kāi)采強(qiáng)度、深度、規(guī)模和資源量的日益增大，與煤礦開(kāi)采有關(guān)的水害問(wèn)題愈發(fā)突出，由此引發(fā)的水害事故也愈發(fā)嚴(yán)重[3]。在深入系統(tǒng)地查明礦區(qū)水文地質(zhì)條件的基礎(chǔ)上，利用多種方法對(duì)礦井涌(突)水進(jìn)行準(zhǔn)確判識(shí)是確保煤礦防治水工作有效開(kāi)展、降低煤礦水害事故損失的前提和基礎(chǔ)，對(duì)保障煤礦安全生產(chǎn)具有重要意義。

目前，基于不同含水層(水體)水化學(xué)特征上的差異性，運(yùn)用模糊數(shù)學(xué)[4]、Fisher判別分析[2,5-6]、Piper三線(xiàn)圖[7]、距離判別分析[8-9]等多種方法對(duì)涌(突)水及其可能來(lái)源含水層(水體)的水化學(xué)數(shù)據(jù)進(jìn)行處理對(duì)比分析，是較為常用的方法。隨著計(jì)算機(jī)技術(shù)的發(fā)展，近年來(lái)機(jī)器學(xué)習(xí)，特別是支持向量機(jī)(support vector machine,SVM)被廣泛應(yīng)用于礦井涌(突)水的水源判別中[10-12]。另一方面，隨機(jī)森林(random forest,RF)由于其良好的穩(wěn)定性和魯棒性，近年來(lái)也被逐漸應(yīng)用采空區(qū)自燃預(yù)測(cè)、礦井涌(突)水水源判別等方面[13-14]。上述方法本質(zhì)都是利用數(shù)學(xué)方法對(duì)涌(突)水及其可能來(lái)源含水層(水體)的水化學(xué)數(shù)據(jù)進(jìn)行處理，忽略了涌(突)水可能來(lái)源含水層(水體)間的循環(huán)演化過(guò)程，使判識(shí)結(jié)果缺少實(shí)際水循環(huán)情況的支撐與驗(yàn)證；另一方面，上述方法僅對(duì)來(lái)源進(jìn)行判識(shí)，缺少對(duì)實(shí)際采礦過(guò)程中涌(突)水現(xiàn)象與礦井立體水文地質(zhì)模型等結(jié)合。因此，筆者提出一種基于水化學(xué)場(chǎng)機(jī)器學(xué)習(xí)分析與水動(dòng)力場(chǎng)反向示蹤模擬耦合的礦井涌(突)水水源綜合判識(shí)技術(shù)，結(jié)合陜西榆林曹家灘煤礦工程背景，在已有數(shù)據(jù)資料的支撐下，利用水文地球化學(xué)的原理和方法對(duì)礦井涌(突)水及其可能來(lái)源含水層(水體)的水化學(xué)特征進(jìn)行分析，一方面利用特征上的相似性，定性對(duì)涌(突)水來(lái)源進(jìn)行判識(shí)，另一方面水化學(xué)特征方面的分析賦予了機(jī)器學(xué)習(xí)定量判識(shí)結(jié)果的實(shí)際意義；此外，結(jié)合實(shí)際情況，建立地下水滲流場(chǎng)，利用GMS軟件中MODPATH反向示蹤涌(突)水水源，實(shí)現(xiàn)了涌(突)水路徑的可視化及對(duì)水化學(xué)機(jī)器學(xué)習(xí)判識(shí)結(jié)果的驗(yàn)證，為礦井涌(突)水來(lái)源判識(shí)提供了一種新的思路。

1 研究區(qū)概況

陜西省榆林市曹家灘煤礦，地處陜北黃土高原北部，鄂爾多斯高原東北部，毛烏素沙漠東南緣，為沙丘沙地和風(fēng)沙灘地、黃土梁峁地貌。曹家灘煤礦位于我國(guó)西北內(nèi)陸，為典型的溫帶干旱、半干旱大陸性季風(fēng)氣候，多年平均降水量535.51 mm，多年平均蒸發(fā)量1 916.1 mm，蒸發(fā)量遠(yuǎn)大于降水量[15]。曹家灘煤礦內(nèi)NW—SE向的分水嶺大致將礦區(qū)潛水劃分為東西2個(gè)面積大致相等且相對(duì)獨(dú)立的水文地質(zhì)單元，西南部屬榆溪河流域，東北部屬禿尾河流域。煤礦東南有屬禿尾河支流的野雞河、高羔兔溝兩小溝流，為季節(jié)性溝流，在溝流的下游因老鄉(xiāng)農(nóng)灌截流常出現(xiàn)斷流，區(qū)內(nèi)原有一些海子多數(shù)已經(jīng)干枯，現(xiàn)存的一些海子水位不深，蓄水量不大。

按照含水介質(zhì)的不同，曹家灘井田含水層自2-2煤層向上可分為侏羅系中統(tǒng)延安組第5段孔隙裂隙承壓含水層、侏羅系中統(tǒng)直羅組孔隙裂隙承壓含水層、侏羅系中統(tǒng)安定組孔隙裂隙承壓含水層、風(fēng)化巖基巖裂隙承壓水含水層、第四系中更新統(tǒng)離石黃土弱含水層和第四系上更新統(tǒng)薩拉烏蘇組孔隙潛水含水層。隔水層主要是新近系保德組紅土，其次為基巖中的泥巖和粉砂巖。其中，保德組紅土為研究區(qū)隔水性能較好且較穩(wěn)定的隔水層，其沿分水嶺兩側(cè)最厚，自東向西、西南逐漸變薄。位于井田東翼的122108工作面和西翼的122109工作面所開(kāi)采的2-2煤煤厚均為10 m左右，同為綜采放頂煤開(kāi)采工藝且已完成回采，2個(gè)工作面面積相近，分別為1.67，1.56 km2。但根據(jù)涌水量數(shù)據(jù)統(tǒng)計(jì)顯示，井田東翼122108工作面平均涌水量約為341.21 m3/h，而西翼122109工作面平均涌水量約為758.33 m3/h，涌水量差異較大。筆者嘗試通過(guò)基于水化學(xué)場(chǎng)機(jī)器學(xué)習(xí)分析與水動(dòng)力場(chǎng)反向示蹤模擬耦合的礦井涌(突)水水源綜合判識(shí)技術(shù)，對(duì)兩工作面的涌水來(lái)源進(jìn)行判識(shí)。其中，122109工作面走向方向剖面示意如圖1所示。

圖1 122109工作面走向水文地質(zhì)剖面示意Fig.1 Hydrogeologic profile of No.122109 working face advancing direction

2 研究區(qū)地下水化學(xué)特征成因分析

地下水中各離子的質(zhì)量濃度是不同含水層信息表征的載體，具有化學(xué)指紋識(shí)別的功能[2]。不同含水層地下水與含水介質(zhì)發(fā)生著不同的水文地球化學(xué)作用，使其具有不同的水化學(xué)特征，這些差異是水化學(xué)特征可作為水源判別依據(jù)的根本原因。涌(突)水水源判別中，掌握當(dāng)?shù)厮瘜W(xué)特征的成因，了解不同含水層之間的演化關(guān)系，是保證判別結(jié)果準(zhǔn)確的基礎(chǔ)。

2.1 水化學(xué)參數(shù)統(tǒng)計(jì)特征

運(yùn)用SPSS軟件對(duì)研究區(qū)地表水和各含水層地下水水樣的水化學(xué)成分進(jìn)行統(tǒng)計(jì)分析，結(jié)果見(jiàn)表1。

表1 水化學(xué)參數(shù)統(tǒng)計(jì)特征值Table 1 Statistical summary of hydrochemical parameters

2.2 水化學(xué)參數(shù)相關(guān)性分析

表2 研究區(qū)地下水水化學(xué)參數(shù)相關(guān)性系數(shù)矩陣Table 2 Pearson’s correlation coefficients of groundwater chemical parameters in study area

2.3 水化學(xué)類(lèi)型及其空間分布

研究區(qū)各含水層地下水的Piper三線(xiàn)圖(圖2)顯示，從空間上來(lái)看，研究區(qū)內(nèi)地表水和各含水層地下水中陽(yáng)離子主要成分變化較大，地表水和第四系地下水以Ca2+為主導(dǎo)，風(fēng)化層水無(wú)主導(dǎo)陽(yáng)離子，直羅組水分布覆蓋Ca2+主導(dǎo)、無(wú)主導(dǎo)和Na++K+主導(dǎo)3種情況，延安組水分布覆蓋無(wú)主導(dǎo)和Na++K+主導(dǎo)2種情況，礦井水以Na++K+為主導(dǎo)，隨著水樣深度的增加，陽(yáng)離子含量呈現(xiàn)由Ca2+為主導(dǎo)過(guò)渡到Na++K+為主導(dǎo)的趨勢(shì)。

從圖2中讀取了它們的水化學(xué)類(lèi)型見(jiàn)表3(M801，M802和M901，M902分別為122108和122109工作面涌水)。地表水和第四系水的水化學(xué)主要類(lèi)型均為HCO3-Ca型(75%和69%)，代表這2類(lèi)水之間水力聯(lián)系較強(qiáng)。風(fēng)化層水的水化學(xué)類(lèi)型比較復(fù)雜，沒(méi)有占明顯主導(dǎo)作用的水化學(xué)類(lèi)型，部分陽(yáng)離子以Ca2+(Mg2+)為主導(dǎo)，部分陽(yáng)離子以Na+(K+)為主導(dǎo)，也可能是由于部分地區(qū)存在紅土薄弱區(qū)，使風(fēng)化層水與第四系水存在水力聯(lián)系。直羅組水的水化學(xué)類(lèi)型以HCO3-Na(50%)為主，可能和陽(yáng)離子交替吸附作用及硅酸鹽礦物的溶解有關(guān)。

表3 研究區(qū)水樣水化學(xué)類(lèi)型統(tǒng)計(jì)Table 3 Statistical table of hydrochemistry types in study area

圖2 研究區(qū)地表水和地下水Piper三線(xiàn)圖Fig.2 Piper diagram of surface water and groundwater in study area

延安組的水化學(xué)類(lèi)型主要為HCO3·SO4-Na(40%)和HCO3-Na·Mg·Ca(20%)，可能和含水層中石膏的溶解以及Ca2+,Na+發(fā)生陽(yáng)離子交替吸附作用有關(guān)。礦井水中122108，122109工作面水化學(xué)類(lèi)型主要為HCO3-Na，礦井水的水化學(xué)類(lèi)型與直羅組和延安組地下水的水化學(xué)類(lèi)型相近，與地表水和第四系水的水化學(xué)類(lèi)型不同，與部分風(fēng)化層水的水化學(xué)類(lèi)型相似。初步判斷，礦井水來(lái)源于直羅組和延安組裂隙含水層的補(bǔ)給，可能有少量風(fēng)化層水的來(lái)源。

3 機(jī)器學(xué)習(xí)判別模型

SVM基于統(tǒng)計(jì)學(xué)習(xí)理論基礎(chǔ)和結(jié)構(gòu)風(fēng)險(xiǎn)最小化原理實(shí)現(xiàn)模式分類(lèi)和非線(xiàn)性回歸[18]。對(duì)于模式分類(lèi)問(wèn)題，其主要思想是建立一個(gè)超平面作為分類(lèi)的決策曲面，使類(lèi)域邊界之間的隔離距離最大，進(jìn)而將不同類(lèi)別數(shù)據(jù)分開(kāi)。SVM的結(jié)構(gòu)如圖3所示，其中K為核函數(shù)，種類(lèi)主要有線(xiàn)性核函數(shù)、多項(xiàng)式核函數(shù)、徑向基核函數(shù)和sigmoid核函數(shù)等。本文中SVM核函數(shù)選取最為常見(jiàn)且已被證實(shí)具有最優(yōu)水源判別能力的徑向基核函數(shù)[11]，其形式為

圖3 SVM結(jié)構(gòu)示意Fig. 3 SVM structure diagram

(1)

式中，σ為核寬度。

RF是一種包含若干個(gè)決策樹(shù)的模型，這些隨機(jī)樹(shù)的形成采用的是隨機(jī)的方法。其基本原理是輸入測(cè)試數(shù)據(jù)后，首先讓每一個(gè)決策樹(shù)進(jìn)行單獨(dú)分類(lèi)，最后選取分類(lèi)結(jié)果最多的那個(gè)類(lèi)別作為最終輸出結(jié)果。RF結(jié)構(gòu)如圖4所示，具體流程：

圖4 RF模型計(jì)算流程Fig. 4 RF model calculation process

(1)重采樣。利用bootstrap方法對(duì)訓(xùn)練數(shù)據(jù)進(jìn)行重采樣，隨機(jī)產(chǎn)生T個(gè)訓(xùn)練集，每次沒(méi)有被抽到的訓(xùn)練數(shù)據(jù)組成T個(gè)袋外數(shù)據(jù)(out-of-bag,OOB)，抽到的訓(xùn)練數(shù)據(jù)為袋內(nèi)數(shù)據(jù)(In bag)；

(2)生成決策樹(shù)。從每個(gè)訓(xùn)練集樣本中的M個(gè)屬性中隨機(jī)挑選m個(gè)屬性，作為節(jié)點(diǎn)分裂屬性集，然后從屬性集中選出最優(yōu)屬性進(jìn)行節(jié)點(diǎn)分裂，構(gòu)建出每個(gè)訓(xùn)練集對(duì)應(yīng)的CART樹(shù)，且每棵樹(shù)都保持完整成長(zhǎng)，不進(jìn)行任何裁枝，即m保持不變；

(3)決策。由于每個(gè)CART樹(shù)在訓(xùn)練集選擇和屬性選擇上都是隨機(jī)的，因此，這T個(gè)決策樹(shù)是獨(dú)立的，將測(cè)試集輸入每個(gè)決策樹(shù)，得到T個(gè)預(yù)測(cè)結(jié)果，對(duì)于分類(lèi)問(wèn)題，利用少數(shù)服從多數(shù)的原則，選擇輸出結(jié)果最多的類(lèi)別為測(cè)試集所屬的類(lèi)別。

3.1 數(shù)據(jù)預(yù)處理

數(shù)據(jù)是否采用歸一化處理，以及不同的歸一化方式對(duì)模型的準(zhǔn)確性有著不同的影響[19]。交叉驗(yàn)證(Cross Validation,CV)是一種測(cè)試分類(lèi)器性能的統(tǒng)計(jì)分析方法，其基本原理是將原始數(shù)據(jù)在某種意義下分組，一部分作為訓(xùn)練集，一部分作為驗(yàn)證集，然后用訓(xùn)練集訓(xùn)練出分類(lèi)模型并用驗(yàn)證集檢驗(yàn)?zāi)Ｐ偷臏?zhǔn)確率。其中K折交叉驗(yàn)證(K-fold Cross Validation)的方法可以有效避免欠學(xué)習(xí)或過(guò)學(xué)習(xí)狀態(tài)的發(fā)生，能得到具有說(shuō)服力的結(jié)果。筆者利用K折交叉驗(yàn)證對(duì)不同數(shù)據(jù)預(yù)處理方式得到的模型驗(yàn)證集準(zhǔn)確率進(jìn)行計(jì)算，進(jìn)而對(duì)不同歸一化方式進(jìn)行選取。結(jié)果(表4)表明，RF和SVM采用[-1,1]歸一化方式對(duì)本文的數(shù)據(jù)進(jìn)行預(yù)處理，可以有效提高模型的精度。

表4 采用不同歸一化方式準(zhǔn)確率對(duì)比Table 4 Accuracy comparison of different normalization methods

3.2 模型參數(shù)的選取優(yōu)化

在使用SVM做分類(lèi)預(yù)測(cè)時(shí)，需要調(diào)整懲罰系數(shù)c和核函數(shù)參數(shù)g，以得到比較理想的分類(lèi)準(zhǔn)確率。與傳統(tǒng)的網(wǎng)格搜索(Grid search)算法尋優(yōu)相比，啟發(fā)式的遺傳算法(Genetic Algorithm,GA)尋優(yōu)不必遍歷網(wǎng)格內(nèi)的所有參數(shù)點(diǎn)，具有較強(qiáng)的魯棒性，且已被證實(shí)可以在復(fù)雜的參數(shù)空間中快速選取到最佳c,g值[10,12]。因此，筆者利用GA算法對(duì)SVM的參數(shù)進(jìn)行優(yōu)化，建立參數(shù)優(yōu)化的GA-SVM模型。利用LIBSVM工具箱[20]、Sheffield遺傳算法工具箱和gaSVMcgForClass函數(shù)[21]實(shí)現(xiàn)懲罰系數(shù)c和核函數(shù)參數(shù)g的尋優(yōu)，計(jì)算得到最優(yōu)懲罰系數(shù)c=1.996 8，最優(yōu)核函數(shù)參數(shù)g=4.093 6。

由RF模型的計(jì)算流程可知，每棵決策樹(shù)在生成時(shí)，沒(méi)有被抽到的樣本數(shù)據(jù)被保留在了OOB子集中，利用OOB子集數(shù)據(jù)中的樣本可以計(jì)算得到每棵樹(shù)的OOB錯(cuò)誤率，則模型中所有決策樹(shù)的OOB錯(cuò)誤率均值可以對(duì)RF模型的性能進(jìn)行評(píng)價(jià)。此外，在決策樹(shù)生長(zhǎng)過(guò)程中，節(jié)點(diǎn)預(yù)選變量m過(guò)多可能導(dǎo)致模型過(guò)擬合，預(yù)測(cè)精度下降[13]。筆者利用MATLAB對(duì)不同決策樹(shù)棵樹(shù)T和節(jié)點(diǎn)預(yù)選變量m對(duì)應(yīng)的OOB錯(cuò)誤率進(jìn)行計(jì)算，結(jié)果如圖5所示。由圖5可以看出，隨著決策樹(shù)棵數(shù)的增加，運(yùn)算量增大，但是OOB錯(cuò)誤率的降低并不明顯，因此，本文所使用的RF模型中選取決策樹(shù)棵數(shù)T=150，節(jié)點(diǎn)預(yù)選變量m=1。

圖5 隨機(jī)森林OOB錯(cuò)誤率Fig.5 OOB error rate of random forest

3.3 模型選取與水源判別

將原始數(shù)據(jù)輸入模型進(jìn)行10次十折交叉驗(yàn)證，得到GA-SVM模型和RF模型的訓(xùn)練集與驗(yàn)證集準(zhǔn)確率如圖6和表5所示。

從圖6可以看出，GA-SVM和RF的訓(xùn)練集和驗(yàn)證集準(zhǔn)確率均處于高位，說(shuō)明2種模型的預(yù)測(cè)分類(lèi)結(jié)果均具有一定的可信度。GA-SVM的訓(xùn)練集正確率在接近100%處小范圍波動(dòng)，驗(yàn)證集的正確率在接近100%處波動(dòng)較大，產(chǎn)生了較大的預(yù)測(cè)誤差，RF訓(xùn)練集正確率均為100%，驗(yàn)證集正確率在接近100%處小范圍波動(dòng)，說(shuō)明與GA-SVM相比，RF可以更好地降低噪聲對(duì)預(yù)測(cè)分類(lèi)結(jié)果的影響，具有更好的魯棒性，也說(shuō)明了RF具有優(yōu)于GA-SVM的穩(wěn)定性。從表5可以看出，RF的訓(xùn)練集準(zhǔn)確率和驗(yàn)證集準(zhǔn)確率分別為100%和98.00%，分別高出GA-SVM訓(xùn)練集和驗(yàn)證集準(zhǔn)確率1.37%和7.83%。此外，SVM在進(jìn)行模型訓(xùn)練時(shí)，參數(shù)c和g對(duì)模型分類(lèi)準(zhǔn)確率影響較大，需要額外利用GA算法等方法對(duì)其進(jìn)行優(yōu)選，增加了模型的復(fù)雜度，而RF無(wú)需復(fù)雜的參數(shù)設(shè)置和優(yōu)化便能得到較為滿(mǎn)意的性能，且決策樹(shù)為非線(xiàn)性處理器，可將RF視為若干個(gè)非線(xiàn)性關(guān)系組合形成的更為復(fù)雜的非線(xiàn)性關(guān)系處理器?？傊?，[-1,1]歸一化后的RF模型在研究區(qū)礦井水水源判別方面與GA-SVM相比具有更優(yōu)的性能。

圖6 GA-SVM和RF10次十折交叉驗(yàn)證結(jié)果對(duì)比Fig.6 Comparison of GA-SVM and RF 10-fold cross validation results

表5 10次十折交叉驗(yàn)證正確率對(duì)比Table 5 Comparison of 10-fold cross validation accuracy

將待判別的122108和122109工作面的礦井水水化學(xué)數(shù)據(jù)進(jìn)行[-1,1]歸一化后輸入訓(xùn)練好的RF模型，即可快速得到預(yù)測(cè)結(jié)果。RF利用各個(gè)決策樹(shù)對(duì)數(shù)據(jù)進(jìn)行預(yù)測(cè)分類(lèi)，且每個(gè)決策樹(shù)給出一個(gè)預(yù)測(cè)分類(lèi)結(jié)果，將每個(gè)預(yù)測(cè)分類(lèi)結(jié)果與所有預(yù)測(cè)結(jié)果的比值稱(chēng)為得分，得分最高的那個(gè)預(yù)測(cè)分類(lèi)結(jié)果即為最終預(yù)測(cè)類(lèi)別。預(yù)測(cè)結(jié)果見(jiàn)表6，可以看出，122108工作面礦井水主要來(lái)源于直羅組水和延安組水，122109工作面礦井水主要來(lái)源于第四系水。

表6 RF預(yù)測(cè)結(jié)果Table 6 RF prediction results

4 剖面數(shù)值模擬

水文地球化學(xué)和機(jī)器學(xué)習(xí)的方法，分別對(duì)研究區(qū)涌(突)水水源進(jìn)行判別，然而2種方法在122109工作面M901和M902水樣的判別結(jié)果上出現(xiàn)了差異。下面將利用滲流場(chǎng)模型可視化示蹤的方法，利用GMS軟件建立研究區(qū)典型二維水文地質(zhì)剖面滲流場(chǎng)數(shù)值模型，進(jìn)一步對(duì)M901和M902的涌水來(lái)源及通道進(jìn)行判別與可視化輸出。

4.1 典型剖面選取

研究區(qū)分水嶺西部潛水總體向西南方向徑流，研究區(qū)紅土層下部裂隙承壓水徑流方向基本順巖層傾向由東向西南方向運(yùn)移。另一方面，根據(jù)勘探鉆孔揭露，研究區(qū)紅土隔水層在西部存在薄弱區(qū)。典型剖面的選取主要目的是從地下水動(dòng)力學(xué)方面分析122109工作面礦井水垂向來(lái)源于第四系地下水補(bǔ)給的可能性，并嘗試給出補(bǔ)給路徑，因此所選取的剖面應(yīng)大致沿地下水徑流西南走向，且穿過(guò)研究區(qū)西部紅土薄弱區(qū)和122109工作面并靠近M901和M902取樣點(diǎn)位置。為了盡可能滿(mǎn)足研究目的，結(jié)合研究區(qū)鉆孔的分布情況，選取SW-NE方向22.1°的剖面為研究剖面(圖7)。

圖7 第四系與風(fēng)化基巖水水化學(xué)類(lèi)型與紅土厚度分布示意Fig. 7 Indications of hydrochemical types of Quaternary and weathered bedrock and thickness distribution of laterite

4.2 剖面地下水?dāng)?shù)值模型建立

由于本次剖面模擬的目的是分析礦井涌水來(lái)源和路徑，不考慮地下水隨時(shí)間改變而變化的情況，且所選剖面上含水層在水平上無(wú)明顯分帶，在垂向上不同深度的巖性有一定的分層性，因此將剖面結(jié)構(gòu)概化為非均質(zhì)各向異性穩(wěn)定流地下水系統(tǒng)。垂向上概化為6層，即第四系含水層、保德組紅土隔水層、風(fēng)化基巖含水層、直羅組含水層、延安組第五段含水層和煤層。模型的上部邊界，由于薩拉烏蘇組含水層水量豐富且剖面附近長(zhǎng)觀(guān)孔監(jiān)測(cè)數(shù)據(jù)顯示水位半年波動(dòng)僅為0.19 m，因此，將剖面上部邊界視為定水頭邊界。模型的下部定義為隔水邊界，東部和西部為控制整個(gè)剖面內(nèi)地下水的流向，設(shè)置為定水頭邊界。由于穩(wěn)定流中地下水水位不隨時(shí)間發(fā)生變化，因此，初始水頭對(duì)模型流場(chǎng)無(wú)影響，設(shè)置為與地表高程一致。模型水平長(zhǎng)5.8 km，垂向上標(biāo)高在948.66～1 288.32 m，因此，將模型Y軸(水平方向)剖分為583列，X軸(厚度)剖分為2列(11.23 m)，Z軸(垂直方向)剖分為31層，共計(jì)36 146個(gè)活動(dòng)單元格。通過(guò)反復(fù)調(diào)參，使剖面滲透系數(shù)符合附近歷次抽水試驗(yàn)成果并與巖性變化規(guī)律相對(duì)應(yīng)，且剖面流場(chǎng)形態(tài)與實(shí)際流場(chǎng)形態(tài)相一致，最終利用GMS軟件運(yùn)行后，得到煤層開(kāi)采前等水位線(xiàn)分布如圖8(a)所示。

圖8 開(kāi)采前后典型剖面水頭分布情況Fig. 8 Water head distribution of typical section before and after mining

4.3 采動(dòng)后涌水來(lái)源情況模擬

煤層開(kāi)采影響下，頂板巖體中，上部巖體的移動(dòng)變形小于下部巖體，因而采空區(qū)上方的巖層會(huì)自下而上發(fā)育導(dǎo)水裂隙帶和彎曲下沉帶。其中導(dǎo)水裂隙帶自下而上可分為垮落帶和斷裂帶。根據(jù)巖層的斷裂、開(kāi)裂及離層發(fā)育程度及導(dǎo)水能力，可以將斷裂帶區(qū)域自下而上進(jìn)一步分為嚴(yán)重?cái)嗔褏^(qū)、一般開(kāi)裂區(qū)和微小開(kāi)裂區(qū)(表7)，其中垮落帶高度約為導(dǎo)水裂隙帶高度的1/4，嚴(yán)重?cái)嗔褏^(qū)和一般開(kāi)裂區(qū)高度約為導(dǎo)水裂隙帶高度的1/2，微小開(kāi)裂區(qū)高度約為導(dǎo)水裂隙帶高度的1/4[22]。煤層上覆巖體的斷裂變形也會(huì)造成其滲透率發(fā)生變化，隨著與開(kāi)采煤層垂直距離的增加，導(dǎo)水裂隙帶的滲透系數(shù)增加幅度逐漸減小，變化幅度由上至下均為開(kāi)采前的1～10倍[22-26]。

表7 斷裂帶分區(qū)情況[22]Table 7 Fractured zone division[22]

通過(guò)在122109工作面施工的“兩帶”鉆孔的手段，以鉆進(jìn)過(guò)程中沖洗液漏失量、鉆孔內(nèi)水位變化、巖芯鑒定、鉆孔電視等多種方法，綜合確定122109工作面垮落帶最大發(fā)育高度為53.99 m，波及到延安組含水層；斷裂帶最大發(fā)育高度為215.01 m，波及到風(fēng)化基巖含水層。因此，結(jié)合導(dǎo)水裂隙帶分區(qū)情況及導(dǎo)水裂隙帶滲透系數(shù)變化情況，將模型對(duì)應(yīng)導(dǎo)水裂隙帶發(fā)育高度下部1/4位置、中部1/2位置和上部1/4位置滲透系數(shù)分別設(shè)置為原有10倍、5倍和2倍。此外，根據(jù)研究區(qū)附近已有紅土層采動(dòng)前后壓水試驗(yàn)成果，雖然導(dǎo)水裂隙帶未波及到紅土層，但應(yīng)力的改變也導(dǎo)致了其滲透系數(shù)變?yōu)殚_(kāi)采前的10倍[27]，將模型對(duì)應(yīng)層位滲透系數(shù)也進(jìn)行調(diào)整。而紅土層上覆第四系松散沉積物主要成分為細(xì)沙和少量亞沙土、亞黏土，位于彎曲下沉帶，研究區(qū)已有成果表明在采動(dòng)影響下，其滲透性變異程度不大，模型中對(duì)應(yīng)滲透系數(shù)不做調(diào)整[28-29]。將調(diào)整后的參數(shù)輸入至模型中，使用GMS軟件中的Drain模塊模擬煤層開(kāi)采后的情況，利用MODPATH模塊的反向示蹤功能對(duì)煤層開(kāi)采后的涌水來(lái)源進(jìn)行模擬，其中開(kāi)采后水頭分布情況如圖8(b)所示，122109工作面礦井水來(lái)源路徑情況如圖9所示。

圖9 開(kāi)采后涌水來(lái)源判識(shí)示意Fig.9 Source identification diagram of water inflow after mining

從模型計(jì)算結(jié)果可以看出，開(kāi)采后，典型剖面流場(chǎng)較開(kāi)采前發(fā)生了明顯變化，受煤層開(kāi)采影響，采空區(qū)附近含水層水頭明顯降低。MODPATH模塊反向示蹤計(jì)算結(jié)果顯示，在紅土薄弱區(qū)存在第四系水越流并通過(guò)導(dǎo)水裂隙帶補(bǔ)給礦井水的情況。

5 討 論

根據(jù)水文地球化學(xué)特征相似性和隨機(jī)森林方法的判識(shí)結(jié)果，122108工作面礦井水的主要來(lái)源為直羅組和延安組含水層地下水。此工作面北部緊鄰的122106工作面“三帶”發(fā)育規(guī)律研究項(xiàng)目，通過(guò)理論分析、現(xiàn)場(chǎng)實(shí)測(cè)、物理相似模擬及數(shù)值計(jì)算等方法，綜合確定導(dǎo)水裂隙帶最大發(fā)育高度為162.0 m，波及到直羅組地層，與水文地球化學(xué)和隨機(jī)森林方法判識(shí)結(jié)果吻合?？梢源_定122108工作面礦井涌水的主要來(lái)源為直羅組和延安組含水層。

水文地球化學(xué)特征相似性判識(shí)122109工作面礦井水主要來(lái)源為深層直羅組或延安組含水層地下水，而隨機(jī)森林方法判識(shí)結(jié)果顯示其主要來(lái)源為第四系含水層地下水。通過(guò)引入典型剖面二維地下水?dāng)?shù)值模擬的方法，發(fā)現(xiàn)存在第四系含水層地下水通過(guò)紅土薄弱區(qū)和缺失區(qū)及導(dǎo)水裂隙帶涌入122109工作面的情況。研究區(qū)內(nèi)保德組紅土結(jié)構(gòu)致密，半堅(jiān)硬狀，透水性差，分布連續(xù)，是較為穩(wěn)定的隔水層。因此，當(dāng)導(dǎo)水裂隙帶發(fā)育至安定組風(fēng)化基巖含水層時(shí)，紅土層是否可以有效阻隔第四系水越流補(bǔ)給安定組風(fēng)化基巖含水層，直接影響到122109工作面礦井水是否存在第四系含水層地下水作為來(lái)源。收集整理研究區(qū)范圍內(nèi)各類(lèi)鉆孔共129個(gè)，插值生成煤礦礦權(quán)范圍內(nèi)紅土厚度分布圖，并將收集整理的水化學(xué)數(shù)據(jù)中，存在坐標(biāo)的14個(gè)第四系含水層取樣點(diǎn)和14個(gè)風(fēng)化基巖含水層取樣點(diǎn)投射至紅土厚度分布圖，最終成圖如圖7所示。

從圖7可以看出，研究區(qū)紅土分布情況為中部及東部較厚，厚度可達(dá)50～134 m；而研究區(qū)西部紅土較薄，厚度普遍在50 m以下，且存在紅土缺失區(qū)，經(jīng)計(jì)算，缺失區(qū)面積為1.13 km2。由前文分析可知，研究區(qū)第四系地下水中陽(yáng)離子為Ca2+，而風(fēng)化層水無(wú)主導(dǎo)離子。結(jié)合圖水化學(xué)類(lèi)型分布情況來(lái)看，在中部和東部紅土較厚區(qū)，風(fēng)化層水中Ca2+質(zhì)量濃度不大，其與臨近的第四系水水化學(xué)類(lèi)型存在明顯差異，而在紅土薄弱區(qū)，特別是在紅土缺失區(qū)，風(fēng)化層水中Ca2+質(zhì)量濃度顯著增加。經(jīng)進(jìn)一步計(jì)算發(fā)現(xiàn)，研究區(qū)第四系水中Ca2+毫克當(dāng)量百分比平均為61.60%，紅土厚度15～134 m處Ca2+毫克當(dāng)量百分比平均為31.77%，而紅土厚度0～15 m處Ca2+毫克當(dāng)量百分比平均為40.82%。對(duì)比分析判斷，在西部紅土缺失區(qū)和薄弱區(qū)，第四系含水層與風(fēng)化層水之間存在水力聯(lián)系?？紤]到紅土缺失區(qū)域與122109工作面存在重疊，重疊面積為0.17 km2，因此，在紅土缺失區(qū)域內(nèi)，存在第四系含水層地下水通過(guò)導(dǎo)水裂隙帶涌入122109工作面的情況，這與隨機(jī)森林判識(shí)和典型剖面二維數(shù)值模擬判識(shí)結(jié)果相吻合。水文地球化學(xué)特征相似性判識(shí)122109工作面涌水來(lái)源為直羅組或延安組含水層地下水，可能是因?yàn)榕凶R(shí)過(guò)程僅參照了K+，Na+，Ca2+，Mg2+等常規(guī)離子質(zhì)量濃度，忽略了pH、TDS、總硬度和總堿度等理化指標(biāo)，導(dǎo)致判識(shí)結(jié)果出現(xiàn)誤差。此外，在實(shí)際生產(chǎn)過(guò)程中，122108工作面平均涌水量約為341.21 m3/h，而122109工作面平均涌水量約為758.33 m3/h，達(dá)到122108工作面涌水量的2.22倍，結(jié)合122109工作面上部存在0.17 km2紅土缺失區(qū)，且第四系薩拉烏蘇組含水層水量大的特點(diǎn)，綜合判斷122109工作面涌水主要來(lái)源于第四系含水層地下水。

6 結(jié) 論

(2)通過(guò)十折交叉驗(yàn)證方法，得出RF方法正確率為98.00%，優(yōu)于GA-SVM的90.17%。此外，SVM需要額外利用GA算法等方法對(duì)參數(shù)c和g進(jìn)行優(yōu)選，增加了模型的復(fù)雜度，而RF無(wú)需復(fù)雜的參數(shù)設(shè)置和優(yōu)化便能得到較為滿(mǎn)意的性能。因此，利用RF得出122108工作面涌水來(lái)源于直羅組和延安組含水層地下水，而122109工作面涌水來(lái)源于第四系含水層地下水。

(3)選取122109工作面取樣點(diǎn)附近典型剖面，利用反向示蹤原理，結(jié)合曹家灘礦水文地質(zhì)條件、工程地質(zhì)條件和抽水試驗(yàn)滲透系數(shù)等信息，建立MODPATH礦井涌(突)水滲流場(chǎng)可視化模型，實(shí)現(xiàn)涌水路徑的可視化輸出，并對(duì)RF判識(shí)結(jié)果進(jìn)行驗(yàn)證。結(jié)果顯示，122109工作面在紅土薄弱區(qū)和缺失區(qū)位置，存在第四系含水層地下水通過(guò)導(dǎo)水裂隙帶涌入工作面的情況。

(4)基于水化學(xué)場(chǎng)機(jī)器學(xué)習(xí)分析與水動(dòng)力場(chǎng)反向示蹤模擬耦合的礦井涌(突)水水源綜合判識(shí)結(jié)果，結(jié)合工程實(shí)例進(jìn)行討論驗(yàn)證分析，發(fā)現(xiàn)判識(shí)結(jié)果與工程實(shí)際情況相吻合，證明該方法具有較高的準(zhǔn)確性和可信程度。最終判斷122108工作面涌水來(lái)源于直羅組和延安組含水層地下水，122109工作面涌水來(lái)源于第四系含水層地下水。判識(shí)結(jié)果可為曹家灘煤礦以后防治水工作的開(kāi)展提供參考依據(jù)。