含水层富水性分区及工作面疏放水后涌水量分段预测

2023-11-06雷方超侯恩科赵凯兴

煤田地质与勘探 2023年10期

杨磊，雷方超，侯恩科，路波，冯栋，赵凯兴

(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煤层烧变岩和风化基岩弱至强含水层、延安组正常基岩弱含水层。

根据榆神府矿区浅埋煤层覆岩“三带”发育高度经验值[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。随着工作面的不断开采，涌水量整体呈现“先增后减”趋势，在区段③开采时达到峰值，接近终采线时涌水量相对减小。目前，区段①-区段③已开采完毕，预测涌水量与工作面实测涌水量相比误差较小，吻合度高。