含水層富水性分區及工作面疏放水后涌水量分段預測

楊 磊，雷方超，侯恩科，路 波，馮 棟，趙凱興

(1.陜煤集團神木紅柳林礦業公司，陜西 神木 719300；2.西安科技大學 地質與環境學院，陜西 西安 710054)

礦井頂板水害是制約我國陜北侏羅紀煤田礦井安全生產的主要因素，預防頂板水害事故的發生成為防治水工作的重中之重，尤其是頂板充水含水層預處理前、后工作面涌水量的精準預測，更是成為了近年來防治水工作的重點討論對象[1]。

目前，針對工作面涌水量的預測方法主要包括解析法[2-3]、水文地質比擬法[4]、數值法[5-6]、相關分析法[7]、灰色系統理論[8]、時間序列[9]等方法。經過長期研究，以往工作面涌水量預測大多是將工作面作為一個整體進行，未考慮工作面內水文地質條件差異對涌水量的影響，預測結果與實際相比誤差較大，因此，許多學者逐漸從不同角度、因地制宜地對工作面涌水量進行更加精細的時空動態預測。陳思佳等[10]根據月掘進進度將工作面分段，提出了基于“分段大井法”的工作面涌水量預測方法；虎維岳[11]分析了工作面頂板涌水模式及涌水量的時空組成特點，預測了多個回采工作面同時生產時不同年份礦井涌水量及其動態變化過程；侯恩科等[12]以陜西檸條塔井田1-2煤層隱伏火燒區下3 個連續開采的工作面為研究對象，利用GMS 軟件建立地下水數值模型，模擬了采動過程中涌水量隨開采時間的動態變化。

陜北侏羅紀煤田煤質優良、地質構造簡單、開采技術條件優越，往往選擇布置超長走向的工作面以提升生產效率，而工作面上覆含水層富水性在空間上存在不均一性，從而引起了工作面在不同采動位置時涌水量的變化，且現階段依據工作面主要充水含水層富水性分區預測工作面涌水量的研究較少。因此，筆者以陜西紅柳林煤礦15217 工作面為研究對象，建立風化基巖及燒變巖含水層的富水性分區模型，依據富水性分區結果對工作面劃分不同富水等級區段，并利用動靜儲量法對井下疏放水后不同區段的涌水量進行分區預測，以期為風化基巖與燒變巖復合含水層充水模式下的工作面涌水量預測提供新思路。

1 工作面概況

15217 工作面位于紅柳林井田東部，工作面傾向寬348 m，走向長2 400 m，5-2煤層厚3.8～5.7 m，傾角0°～3°，工作面地表以波狀沙丘地貌為主，受雨水沖刷形成多條沖溝，匯入毛驢灘溝后沿南偏東方向流經工作面，如圖1 所示。該區域5-2煤層賦存整體為南高北低、西高東低。15217 工作面呈南北向分布，其西側為尚未開采的15218 工作面，東側為已開采的15216工作面采空區。15216 工作面5-2煤層底板整體低于15217 工作面，采空區積水整體向東排泄，對15217 工作面安全生產不構成威脅。

圖1 15217工作面井上下對照及A-A’水文剖面Fig.1 Ground and underground conditions contrast and A-A' hydrographic profile of working face 15217

侏羅系延安組5-2煤層上覆含隔水層自上而下依次為富水性極弱的第四系沙層與黃土、新近系保德組紅土隔水層、延安組4-2煤層燒變巖和風化基巖弱至強含水層、延安組正?；鶐r弱含水層。

根據榆神府礦區淺埋煤層覆巖“三帶”發育高度經驗值[13]與相鄰礦井實測數據，采用28 倍裂采比計算導水裂隙帶最大發育高度，得出15217 工作面導水裂隙帶發育至紅土隔水層以上，局部已溝通地表。

根據毛驢灘溝附近松散層以往水文探查結果，15217 工作面及附近區域松散層不含水。工作面采前對地表毛驢灘溝流水進行了截流–導流處理，對低洼積水區進行了疏排，采中對地表裂縫進行填埋治理，保證大氣降水與地表水不會涌入工作面。根據以往開采工作面涌水量資料，雨季工作面涌水量未出現明顯增大。因此，松散層水、地表水及大氣降水對15217 工作面涌水量影響較小。

風化基巖與燒變巖含水層是15217 工作面主要頂板充水含水層。風化基巖厚度4.50～23.95 m，巖性以風化砂巖為主，其富水性受巖性、厚度及裂隙發育程度影響，具有明顯的分區特征。5-2煤層上覆的4-2煤層沿煤層露頭邊界向賦存區內自燃而形成較大范圍燒變巖，燒變巖厚度0～11.1 m，受燒變程度及補給條件影響，富水性在平面上表現出不均一性。在4-2煤火燒區內部，燒變巖與風化基巖直接接觸，兩者水力聯系密切。因此，15217 工作面受4-2煤燒變巖與風化基巖復合含水層水害威脅。

2 含水層富水性影響因素及等級

2.1 影響因素

通過對15217 工作面所在區域風化基巖與燒變巖含水層發育特征的研究，結合抽水試驗成果對比分析，綜合選定含水層厚度、含水層巖性組合指數、巖石燒變及風化程度指數、巖心采取率作為含水層富水性評價的主控因素[14-16]，建立量化指標進行富水性精確預測。

2.1.1含水層厚度(D)

風化基巖與燒變巖含水層因裂隙、孔隙發育導致儲水性與透水性增強，一般而言，當其他影響含水層富水性的因素不發生變化時，含水層厚度越大，富水性越強。

2.1.2巖性組合指數(O)

研究區風化基巖及燒變巖均由單一巖性或多層不同巖性的巖層組成，巖性以砂巖與泥巖為主。抽水鉆孔實測單位涌水量資料表明：對于不同巖性巖層形成的風化基巖，脆性砂巖的單位涌水量明顯大于塑性泥巖，其根本原因在于砂巖經風化后孔隙率增大且裂隙發育，形成連通性好的風化裂隙，而泥巖風化后松軟且遇水易泥化膨脹，對風化裂隙具有一定壓實彌合作用；對于同一巖性的風化基巖，碎屑粒徑越大的巖層單位涌水量也越大，如風化后粗砂巖較細砂巖的單位涌水量明顯增強。

因此，根據巖體巖性對富水性的影響程度將其量化，結合其在含水層厚度中的占比，建立巖性組合指數O，O與富水性呈正相關關系，公式如下：

式中：n為依據巖性將巖體劃分的層數；u為第u層，1≤u≤n；Du為第u層巖層厚度，m；λu為第u層巖層巖性量化值，按巖體巖性不同，將泥巖、砂質泥巖、粉砂巖、細砂巖、中砂巖、粗砂巖分別用數值1～6 表示。

2.1.3巖石燒變及風化程度指數(G)

基巖受風化程度不同，其巖體結構、裂隙發育情況、巖體力學性質和含隔水性均存在顯著差異。巖體受風化程度越強，內部孔隙、裂隙越發育，單位涌水量越大。燒變巖由于內部發育連通性裂隙孔洞，構成了良好的儲水場所，在受補給強的區域，其富水性相較于風化基巖往往更強。據此特征建立巖石燒變及風化程度指數G，公式如下：

式中：δu為第u層巖石燒變及風化程度量化值，將弱風化、中等風化、強風化、燒變巖分別用數值1～4 表示。

2.1.4巖心采取率指數(P)

在相同的鉆探工藝下，巖心采取率與巖石裂隙發育程度息息相關。巖心采取率越低的巖層，巖體結構越破碎，導水能力相對越好；根據巖心采取率與含水層富水性的關系建立巖心采取率指數P，以此表征基巖含水層富水性強弱，公式如下：

式中：τu為巖心采取率，%。

2.2 富水性等級劃分

根據《煤礦防治水細則》規定，含水層富水性按照抽水鉆孔單位涌水量q劃分為4 個等級：弱富水性(q≤0.1 L/(s·m))、中等富水性(0.1 L/(s·m)5.0 L/(s·m))。

本次選取15217 工作面所在區域51 組風化基巖與燒變巖含水層抽水鉆孔資料組成訓練樣本集，其中未發現極強富水性鉆孔，弱富水性鉆孔占比56.86%。富水性等級劃分在《煤礦防治水細則》規定的基礎上，將原有弱富水性等級(q≤0.10 L/(s·m))細化為弱富水性等級(0.01 L/(s·m)

表1 富水性等級量化值及樣本容量Table 1 Quantitative values of water richness grade and sample size

表2 訓練集樣本指標及實測富水性數據Table 2 Sample indexes and measured water richness data of training set

3 基于蝠鲼覓食算法優化支持向量機的含水層富水性預測

3.1 支持向量機原理

支持向量機(Support Vector Machine，SVM)是監督學習算法中應用廣泛、效果良好的一類算法，針對小樣本、非線性問題有優秀的泛化與學習能力，其解釋結果具有良好的推廣性。原始的SVM 致力于尋求可將特征空間上二元數據間隔最大化的超平面(決策邊界)，后續通過選取核函數將二維向量映射至高維特征空間中尋找最優的超平面，從而解決線性不可分的問題[17-19]。

給定訓練樣本集W：

SVM 解決線性不可分問題的目標函數如下：

約束條件為：

式中：i=1,2,···,m，m為樣本個數；ω為權值向量，決定超平面方向；b為偏置參數，決定超平面與原點距離；ξi為松弛變量；C為懲罰因子。

通常利用拉格朗日函數將該問題轉化為對偶問題繼續求解，其對偶問題的目標函數為：

約束條件為：

式中：α為拉格朗日乘子；K(xi,xj)為核函數，本次計算選取徑向基核函數：

式中：σ為核參數。

在支持向量機模型中，懲罰因子C及核參數σ分別影響著誤差的容忍程度和模型的復雜程度，當C值過大時，預測結果更容易出現過擬合現象，反之，容易導致欠擬合；當σ值越小時，模型越簡單，對預測集的泛化程度越高，反之，模型越復雜且泛化程度越低。懲罰因子C及核參數σ的取值決定了支持向量機預測結果的準確性，選取合適的參數取值是提高預測結果準確性的前提。因此，筆者擬采用蝠鲼覓食算法優化支持向量機參數以提高預測結果的準確性。

3.2 蝠鲼覓食優化算法

蝠鲼覓食優化算法(Manta Ray Foraging Optimization，MRFO)是Zhao Weiguo 等[20]在2020 年提出的一種基于群體優化的算法，其核心是模擬蝠鲼的鏈式、螺旋及翻轉3 種覓食策略來促進數學模型達到全局最優，以實現簡便、自適應強、參數少、收斂速度快的特點而備受關注[21-22]。

(1) 鏈式覓食：蝠鲼群體覓食時，各蝠鲼個體從頭到尾有序排列形成一條覓食鏈，該覓食鏈中的每個蝠鲼個體不僅向著食物源前進，也向著前一個個體前進，這代表著MRFO 在每次迭代時，都會更新最佳解決方案及其前面的解決方案。

(2) 螺旋覓食：當排列成一條覓食鏈的蝠鲼群體發現食物源時，排列有序的蝠鲼個體以食物源為源點，首尾相連呈螺旋式行進軌跡靠近食物源，使得該算法具有良好的局部搜索功能。

(3) 翻轉覓食：各蝠鲼個體會以食物源位置為原點，通過圍繞原點的翻轉運動使個體位置與最佳位置的距離逐漸減小，迫使所有個體在搜索空間中逐漸逼近最優解。

3.3 富水性評價模型構建

本文選取MRFO 優化SVM 模型預測含水層富水性，核函數選用徑向基核函數K(xi,xj)，具體計算流程如圖2 所示。為了消除各富水性影響因素之間的量綱影響，方便數據指標之間的可比性，在模型初始階段對51 組訓練樣本數據(表2)進行歸一化處理，數據歸一化結果見表3。

表3 訓練集樣本指標歸一化結果Table 3 Indicator normalization results of training set samples

圖2 MRFO 優化SVM 計算流程Fig.2 Calculation process of SVM optimization by MRFO

在數據歸一化的基礎上，選擇合適的初始參數是模型尋優的前提，種群規模關系著模型全局尋優的準確性與效率。較大的種群有助于模型尋找全局最優解，但同時伴隨著收斂速度慢的特點，反之，較小的種群規模使模型快速收斂，但所求結果可能為局部最優。且一般而言，當模型迭代次數為種群規模的1～10 倍時，更利于算法尋優。由于訓練樣本數據量較小，本次MRFO 的初始參數設定如下：種群規模m取20，最大迭代次數T取100。

在完成模型初始參數設定的基礎上，將歸一化后的富水性主控因素數據隨機生成初始化種群，利用SVM 計算訓練集的適應度，結合MRFO 中鏈式、螺旋、翻轉3 種覓食優化方式完成現階段最優適應度個體位置的迭代更新，并源源不斷地將最新一次的個體位置回代入SVM 計算訓練集的適應度，當算法滿足迭代終止條件時，利用最優適應度對應的參數計算驗證集適應度，并輸出預測集富水性結果。

3.4 富水性分區結果

本次選取工作面內及周邊14 組水文孔數據作為驗證集樣本，利用Python 語言下sklearn 工具構建MFRO 優化SVM 模型對樣本集進行訓練，同時選取鯨魚優化(Whale Optimization Algorithm，WOA)、粒子群優化(Particle Swarm Optimization，PSO)、麻雀搜索算法優化(Sparrow Search Algorithm，SSA)及獅群優化(Loin Swarm Optimization，LSO)作為對比訓練模型。

從基于各種優化算法改進的支持向量機模型預測所得鉆孔富水性結果來看，MRFO 與WOA 優化算法適應度均可達0.86，預測效果最好，其余優化算法所預測富水性適應度均低于該水平，見表4。

表4 驗證樣本不同優化算法預測結果及適應度Table 4 The fitness and prediction results of the verification samples under different optimization algorithms

就各優化算法達到最優結果的效率而言，具有最優適應度的MRFO 與WOA 相比較，MRFO 于第3 次迭代已達最優，WOA 于第15 次迭代達到最優，如圖3所示。因此，基于MRFO 優化SVM 的含水層富水性預測模型效果最優。

圖3 適應度-迭代次數關系Fig.3 Fitness-number of iterations relationship

本次選取MRFO 優化SVM 模型對15217 工作面內及周邊17 組探煤孔數據組成的預測樣本集的風化基巖與燒變巖含水層富水性進行預測，按照探煤孔富水性預測結果，結合已有抽水試驗鉆孔對風化基巖與燒變巖含水層富水性進行分區，分區結果如圖4所示。

圖4 基于MRFO 優化SVM 的富水性分區Fig.4 Waterrichness partition based on SVM optimized by MRFO

4 基于井下疏放水的富水性分區結果驗證

4.1 工作面富水性分段

在風化基巖與燒變巖含水層富水性分區預測的基礎上，基于富水性分區結果將15217 工作面劃分為4 個區段(圖4)：區段①整體屬于弱富水區，受風化基巖含水層影響；區段②整體位于弱–中等富水區；區段④位于中等–極弱富水區，兩段均受風化基巖含水層及局部賦存的燒變巖含水層影響；區段③主要位于強富水區，整體受到大范圍燒變巖及風化基巖含水層影響。

4.2 分區結果驗證

為驗證富水性分區結果的準確性，結合15217 工作面采前井下疏放水鉆孔實測資料，通過對井下疏放水鉆孔初始涌水量、各鉆場疏放水強度與時間、各區段總疏放水量占比的分析，驗證富水性分區結果的準確性。

初始涌水量Q0>10 m3/h 的鉆孔主要位于區段③，集中分布于強富水區，且該區段為燒變巖主要賦存區域；區段④也分布少量Q0>10 m3/h 的鉆孔，其中T22 鉆場位于燒變巖區域，區段①與區段②探放水鉆孔均為Q0≤10 m3/h，如圖5 所示。

圖5 工作面富水性分段及井下疏放水鉆孔布置Fig.5 Layout of waterrichness section on working face and boreholes for drainage

從各鉆場的疏放水強度與時長來看，位于強富水區的大部分疏放水鉆場涌水量整體較大且持續時間較長，位于中等富水區和弱富水區疏放水鉆場涌水量較小且持續時間短。就Q0>10 m3/h 的疏放水鉆場而言，大部分鉆場位于燒變巖強富水區(T11、T12、T15、T16、T17、T19、T20 鉆場)，涌水量大、持續時間長；少量鉆場位于燒變巖中等富水區(T21、T22 鉆場)，涌水量較大、持續時間較短；僅T26 鉆場位于風化基巖弱富水區，但其涌水量很小且持續時間較短，如圖6 所示。

圖6 Q0>10 m3/h 疏放水鉆場放水量歷時曲線(2022 年)Fig.6 Q0>10 m3/h drainage drill sites water release time curves in 2022

工作面各區段疏放水鉆孔總放水量如圖7 所示，各區段總放水量的特征為：區段③>區段④>區段②>區段①。這一結果表明，富水性較強區域的放水量大于較弱區域，其根本原因在于該區域燒變巖含水層富水性強于風化基巖含水層，從而使得溝通燒變巖含水層的疏放水鉆孔疏放水量往往大于僅溝通風化基巖含水層的鉆孔。

圖7 工作面不同區段總放水量及占比Fig.7 Total discharge water and its proportion in different sections of working face

經過井下疏放水實測資料驗證，本次富水性分區結果較為準確，因此，基于MRFO 優化SVM 的富水性預測方法準確度較高。

5 基于富水性分區的工作面涌水量分段預測

5.1 疏放水后水文地質條件分析

15217 工作面采前經過長時間井下疏放水，上覆基巖含水層已經由承壓水轉為無壓水。如圖8 所示，工作面上覆風化基巖含水層水位高程整體高于含水層底面高程，仍存在一定殘余水頭；燒變巖含水層底面高程在區段②與區段③區域基本低于基巖含水層水位高程，該區域燒變巖含水層已基本疏干，但區段④區域仍有一定殘余水頭，未疏放完畢。在工作面開采階段，導水裂隙帶發育至風化基巖與燒變巖含水層時，存在的殘余水頭仍會影響工作面的安全回采。

圖8 風化基巖底面–燒變巖底面–基巖含水層水位等值線Fig.8 Isoline of water level of weathered bedrock bottom-burned rock bottom-bedrock aquifer

5.2 工作面疏放水后涌水量分段預測

目前，大井法是計算工作面涌水量最常用的解析方法，該方法將工作面采后形成的采空區近似看成一個理想化的具有一定影響半徑的“大井”，含水層水源源不斷流入該“大井”的水量即為工作面的涌水量。但是由于傳統大井法僅考慮了含水層水流入采空區的側向補給水量，忽略了含水層自身賦存的水在采動裂隙溝通下向采空區的釋放水量，導致所計算的涌水量比實際涌水量普遍較小。鑒于此問題，本文采用動靜儲量法[23]計算井下疏放水后的工作面涌水量，將傳統大井法中采空區側向補給水量作為動態補給量，加入工作面采后垮落區內流入的含水層自身賦存的水量作為靜態釋放水，兩者之和即為動靜儲量法所計算的工作面涌水量。

由于工作面采前經過長時間的井下疏放水，工作面開采前基巖含水層水已由承壓水轉化為無壓水，因此利用大井法潛水公式計算動態補給水量，計算公式如下：

式中：Qd為動態水補給量，m3/h；H為初始水位高程，m；S為含水層水位降深，m；R0為引用影響半徑，m；r0為引用半徑，m；K為含水層滲透系數，m/d，選取各區段內及附近各水文孔抽水試驗所得滲透系數的算術平均值作為該區段計算動態補給水量所用K值。

15217 工作面疏放水后各區段動態水量計算結果見表5。動態補給量在區段①-區段③呈正增長趨勢，主要原因是采空區面積不斷擴大引起影響半徑的增大，當開采至區段④時，該區域附近為風化基巖底面抬升區，底面最低高程相較于區段③增大12 m，導致了該區域風化基巖含水層水位降深減小，從而使動態補給量表現出下降的趨勢。

表5 工作面各區段動態水量預計結果Table 5 Dynamic water yield calculation results for each section of the working face

靜態釋放水量計算公式如下：

式中：Qj為靜態釋放水量，m3/h；F為采空區面積，m2；L為回采長度，m；B為工作面寬度，m；M為承壓含水層厚度，m；μ為含水層貯水系數；t為預計回采時間，h。

由于采前工作面上覆風化基巖與燒變巖承壓含水層經長期疏放水已轉變為無壓含水層，當工作面采后導水裂隙帶溝通上覆含水層時，僅有殘余水頭范圍內賦存水量流入工作面，當殘余水頭完全釋放入工作面時，殘余水頭高度等于含水層水位降深(S)。潛水含水層靜態釋放水量計算公式如下：

15127 工作面于2022 年6 月下旬開始回采，區段①于7 月末回采完畢，區段②于9 月上旬回采完畢，區段③于11 月末回采完畢，區段④尚未回采，在本次靜態釋放水量計算過程中，根據以往工作面回采速度及15217 工作面回采接續規劃，選取工作面回采速度10 m/d 預計回采時間。

需要說明的是，在計算后續區段的動態水量時，應考慮已采多個區段已形成的采空區對概化“大井”影響半徑的影響，同時，已采多個區段已經完成自身賦存水量的釋放，因此，后續區段的靜態釋放水量的計算可以忽略已采區段靜態水的影響。15217 工作面疏放水后各區段靜態水量計算結果見表6。

表6 工作面各區段靜態水量計算結果Table 6 Static water yield calculation results for each section of the working face

工作面涌水量(Q)應為上覆含水層動態補給水量(Qd)與靜態釋放水量(Qj)之和，工作面各區段涌水量計算結果見表7。隨著工作面的不斷開采，涌水量整體呈現“先增后減”趨勢，在區段③開采時達到峰值，接近終采線時涌水量相對減小。目前，區段①-區段③已開采完畢，預測涌水量與工作面實測涌水量相比誤差較小，吻合度高。

表7 工作面各區段涌水量計算結果Table 7 Water inflow calculation results for each section of the working face

6 結論

a.針對具有密切水力聯系的風化基巖與燒變巖復合含水層，提出了基于蝠鲼覓食算法優化支持向量機富水性預測方法。

b.通過對井下疏放水鉆孔初始涌水量、疏放水鉆場放水強度與時間、各區段總疏放水量的分析，驗證了富水性分區結果的準確性，證明了基于蝠鲼覓食算法優化支持向量機的富水性預測方法準確性較高。

c.按照富水性分區結果對工作面劃分了不同富水等級的區段，在綜合分析工作面疏放水后水文地質條件的基礎上，利用動靜儲量法對工作面疏放水后的涌水量進行了分段預測，預測結果與涌水量實測結果較為吻合。