基于Levenberg Marquardt 反演的福建小煤矿底板富水性分析

2022-02-10邱占林曾东富吴超凡郭玉森刘鑫尧孙轶群黄俊南

煤矿安全 2022年12期

邱占林，曾东富，吴超凡，郭玉森，刘鑫尧，孙轶群，黄俊南

（1.龙岩学院资源工程学院，福建龙岩 364012；2.福建师范大学地理科学学院，福建福州 350007；3.福建紫金铜业有限公司，福建龙岩 364204；4.福建省 197 地质大队，福建泉州 362018）

一直以来，福建小煤矿安全、高效生产受到矿井水害的严重制约[1]，其主要问题在于采掘之前所采用的钻探等传统探水手段的探查能力、精度和范围有限[2]。物探技术近年来虽然在福建小煤矿水害探测中得到越来越广泛的应用，并取得了一定的成效[3-4]，但还是局限于运用常规直流电法、高密度电阻率法、矿井瞬变电磁法等[5-6]进行巷道顶底板及迎头前方探测，且一般在单条巷道中采集易受制于岩性等综合影响的电磁性参数，不能快速、全面、准确定位，不能有效圈定富水区的二维平面（2D）或三维空间（3D）形态，尤其是受底板水害威胁严重的地区，其应用效果更差。传统直流电法一般采用“跑极”的方式，效率低下，精度不高；高密度电阻率法采用串行单向循环数据采集方式，且只能探测巷道单条电剖面测线的2D 地电场信息[7]，或者精度不高的3D地球物理场信息，局限性较为显著；而施工便捷且不接地的矿井瞬变电磁法由于对探测工作环境条件要求较高，其结果易受到井下钢轨、金属支护材料或地电磁场变化等影响[8-9]，数据可靠性较差。而且，上述方法还受到数据采集时间差或延迟时间的干扰，其相应的校正方法也尚不成熟；在电阻率数据反演方面通常也采用阻尼最小二乘法进行迭代计算，但其反演涉及到非线性函数（或方程）的求解问题，过参数化、冗余参数取舍等问题，易导致收敛效果差，并出现相邻网格单元之间电阻率参数反演值畸变现象，往往也会造成迭代计算结果表现出非唯一性（多解性）[10]。相应的存在搜索线性求解方向中止或使损失函数只得局部极小值的估算，从而引起反演精度降低。

网络并行电法自运用于福建小煤矿水害探测以来，其快速、瞬时、同步、连续等数据采集优点，对于削弱游散电流等干扰因素，提高数据信噪比及可信度方面具有较大的优势[3-4,11-12]，但针对福建小矿井富水异常体的位置、分布形态等参数拟合反演还存在一定的缺陷。LMM（Levenberg Marquardt Method）是由Levenberg 提出并经Marquardt 发展完善的一种带阻尼且介于梯度下降法与高斯-牛顿法的非线性或病态优化算法[13-14]。主要利用近似于二阶导数值信息，快速收敛逼近最优解，可有效地拟合地质异常体。克服了雅克比矩阵奇异或梯度下降法、高斯-牛顿法单一求解难以收敛或偏离局部极小值等问题，使反演逼真响应效果和计算精度得到极大提升，可有效解决网络并行电法数据反演中出现的上述问题，使得LMM 在电法数据反演领域得到越来越广泛的应用。因此，运用LMM 反演的网络并行电法开展矿井底板富水性探测，获取待测区域全空间地电场信息，进而查明其富水性的基本状况，可为矿井回采工作面底板水害防治提供一定的水文地质依据。

1 LMM 反演理论及模型构建

LMM 反演是一种基于高斯-牛顿法和梯度下降法改进并用于解决非线性或病态方程组的最优化方法[13]。是通过引入信赖域的迭代参数，采用线性函数来近似（逼近）非线性或病态函数，再模仿阻尼最小二乘法求解目标函数G（m），并在每一步迭代计算过程中寻求合适的变阻尼因子参数λ 以适用自适应求解的最优化，确保目标函数线性搜索及其数值是下降的，以此克服雅可比矩阵几乎全是奇异矩阵或弥补不良条件下高斯-牛顿法不存在逆矩阵、梯度下降法初始值偏离局部极小值太远等缺陷。而在LMM 反演计算中，每次迭代均要对阻尼因子λ 进行自适应调整，而后将上述2 种方法混合改进使用，一般为λ 非常小时，反演算法变为改进的高斯-牛顿法之最优步长计算式；反之，则变为梯度下降法之最优步长计算式，或两者交叉使用之最优步长计算式。只是在局部迭代步长上稍有差别，但均可求取目标函数极小值时的参数向量并进行最优化处理，可使并行电法电导率σ（电阻率ρ 之倒数）反演赋值以及优化时过参数化问题等得以合理解决。

常规网络并行电法数据反演一般属于多元非线性病态函数求解极值范畴[14]，其基本思想是将其转化为大型线性方程组进行数值迭代计算，通常采用不动点迭代分析法、梯度下降法、（准）牛顿法、高斯-牛顿法等方法，并引入雅克比矩阵J 对单位矩阵实施对角分解，大多是基于均匀半空间地球物理场构建目标函数G（m）。

式中：ρgc为观测视电阻率值；ρ（m）为正演理论视电阻率值；Wρ为ρ 的协方差权重系数矩阵；α 为光滑约束模型线性搜索因子；m 为正演理论值初始模型向量；△m 为m 基于迭代计算步长的模型增量（修改）向量；H 为粗糙度矩阵；J为雅克比Jacobi 矩阵。

同时依据LMM 法正则化，可实现反演迭代计算过程中λ 的自适应化。现将式（1）改为：

在信赖域内，光滑约束模型线性搜索因子α 一般会对搜索方向及迭代收敛产生梯度下降法和高斯-牛顿法最优化处理效果。当α＞0 且系数矩阵正定时，保证了最速下降效应，从而使收敛速度得到很大提升。当α=λ 时，式（2）则可转换为：

此时光滑约束模型矩阵S 可由Cole-Cole 模型[15]变换替代，即：

式中：ρ0（m）为正演电阻率理论初始值；η 为极化率；ω 为角频率；τ 为弛豫时间常数；γ 为频率因子。

根据式（3）求取模型增量（修改）向量△m 后，在给定正演理论值初始模型向量m 的基础上，即可通过式（5）求出新的数学模型向量。

2 LMM 反演的三维探测技术原理

2.1 三维探测技术原理

基于三维电阻率快速层析成像理念（CT）的网络并行电法探测技术，除了改进传统直流电法观测系统之外，还加入了拟震式采集算法，实现了快速、高效的探测效果。但在数据处理方面，还存在求解多元非线性病态函数极值过程中冗余数据信息取舍及其影响精度等问题。通过引入LMM 反演数学模型则可有效解决上述问题，形成针对全空间地电场变化的勘探技术。探测矿井底板富水性的观测系统通常布置在运输-回风巷道中，开切眼处由于巷道空间所限不敷设电极仅用作辅助观测，工作站则视探测长度而定，由此可产生地下全空间三维CT 透视快速扫描系统并实施数据采集与处理。具体为：在双巷中布设供电-测量并行电法系统，并将无穷远极B 安装在与供电电极A 不在同一巷道的另外1条巷道中，以保证形成类似于坑透法所扫描的扇形透视观测区[16]，并行电法三维探测施工布置如图1。

图1 并行电法三维探测施工布置Fig.1 Construction layout of three-dimensional detection by parallel electrical method

测量时，供电点A 极（1#～64#）逐点顺次供电，而其余电极作为测量电极则通过程控开关自动实现AM 法或ABM[5,7]法对全空间地电场的实时、同步、快速测量，完成对双巷空间的全电场3D 电透视观测。当1 条巷道采集数据结束时，将其中的观测装置与B 极对换，再次进行采样，即可获得海量拟震式全空间地电场参数，并通过数据拼接以达到对矿井底板的多次覆盖扫描测量，从中可提取有效地质异常体以供后续解译。当测量距离较大时，工作站可布设多个，将测线号（i-j）作为巷道号-工作站号的现场记录标识。施工中，前1 个工作站最末尾电极跟后1个工作站起始电极位置重合，并依次同向移动另一巷道中的无穷远极B（i-j），形成实时、同步测量以保证探测精度。

2.2 三维探测数据反演

基于LMM 的电阻率三维反演的一般思路则为全空间地电场并行电法数据集的迭代数值计算，通常将全部实测数据集转换、拼接、融合，在给定视电阻率矩阵拟合精度ε（阈值）的前提条件下，构建地电场正演理论值初始模型m，并将其剖分为可用于赋值的3D 空间网格，各网格单元宽度通常为0.5～1个电极间距，同时求取基于迭代计算步长的模型增量（修改）向量△m，构建全空间地电场反演基础的目标函数G（m）并对模型修改量求偏导。在此过程中还要考虑粗糙度矩阵H 的3D 分量系数。同时考虑到地下探测空间地电场介质的不均匀性，一般将该空间网格同步划分转换成1 个大型矩阵组，通过引入单位矩阵I，并将光滑约束模型线性搜索因子α变换为基于LMM 正则化的变阻尼因子λ，然后采用迭代求解对角转换后的高次方程组，获得上述网格单元反演所要求的正演电导率σ（即电阻率倒数值ρ（m））等参数。需要注意的是，在数据LMM 反演迭代处理时，需要把笛卡尔坐标系限定在井下3D 空间地电场范围内，从单位矩阵推导出光滑约束模型S对迭代反演计算进行信赖域内的校正，并选取ρ0（m）、η、ω、τ、γ 等Cole-Cole 模型参数来表征反演三维模型，据此求得数学模型向量m(n+1)。重复上述迭代计算过程即可完成3D 视电阻率数据的反演。而且，由于LMM 充分利用了近似二阶导数项信息以及变阻尼因子的自适应化使反演计算精度得到提升。

3 现场应用

3.1 矿井概况

黄土坑煤矿位于龙岩市东北侧岩山乡佳山村，地处政和-大埔深大断裂西侧闽西南聚煤盆地中部含煤条带北端。以中二叠统童子岩组（P2t）一、三段含煤[17]，其中37、38#、39#煤层主采且稳定。开拓方式为平硐+暗斜井，主要采取走向长壁后退式采煤法[18]，采用单体液压+木支柱支护，年产量15 万t/a，属于福建省小型矿井规模，采掘工作面为“2 采4 掘”，其中+120-37#N 采煤工作面为本次试验场地。该矿区构造形态总体呈轴向为NWW-SEE 的背斜构造，同时伴生发育有正、逆断层共4 条，其中F6逆断层和F8正断层对采掘生产有一定的影响。而矿井水文地质类型总体为中等-复杂型，区内地下水主要为大气降水补给，充水水源则以37#煤层顶底板砂岩裂隙水为主，其次为中二叠统栖霞组（P2q）底板岩溶裂隙水，而37#煤层底板下至栖霞灰岩顶面距离约为70 m。岩石中的各种构造裂隙可构成矿井的充水通道，并对矿井充水具有普遍意义。

3.2 数据采集与处理

为了查明黄土坑煤矿+120-37#N 采煤工作面的底板富水性分布情况，2021 年1 月15 日在矿方地测科的配合下，运用网络并行电法勘探技术对该工作面底板视电阻率ρ 展布特征实施三维全空间探测。现场数据采集规定以开切眼与运输巷交点位置为笛卡尔坐标系原点（0，0），沿运输巷指向F8断层揭露点（173，0）作为x 轴正向，沿开切眼指向F6断层揭露点（0，132）作为y 轴正向。井下现场采集选用NPEI-DHZI-1 型网络并行电法仪，应用AM 方式采集，各站依次使用0.5、2.0 s 恒流供电方波测试1次[7]，工作面双巷则根据图1 观测系统各布置3 个工作站，运输巷布设第1～第3 观测站，回风巷则为第4～第6 观测站。第3、第6 观测站敷设55 各电极（测线长度270 m），其余各站均为64 个电极（共315 m），电极距均为5 m。同时做好复测以保证ρ 数据实时、同步采集的准确性与可靠性。

在数据处理阶段，首先，要做好现场观测系统坐标系的设置和全空间地电场数据类型的选取；其次，基于LMM 进行电阻率3D 反演以及电流-电压模块的转换，同时需要获取初始值的正演理论模型构建与解析、3D 空间网格剖分以及二次场参数替代等，而光滑约束模型处理及其线性搜索因子转换以寻求变阻尼因子及其自适应化则是该流程的关键，干扰因素分析与选择性剔除则是数据体整体质量的保证；最后，从反演模型结果中提取自适应迭代电阻率数据体水平切片并立体构建成图。

3.3 测试结果分析

对采煤工作面双巷所采集的全电场数据进行全空间电阻率LMM 反演切片提取成像，可较直观反映采煤工作面底板下不同层位岩层的电性分布情况与变化规律，可三维呈现相对低阻异常区的空间分布与总体连通状况，也可通过提取工作面底板下不同位置视电阻率异常的水平切片及其组合进行表征，据此再结合巷道工程揭露的实际地质情况，对该工作面底板下的富水性分布规律及连通性情况进行全盘分析。但底板下通常因岩性的不均一，表现出不同层位视电阻率阈值响应不一致。因此，常采取同一深度视电阻率值的相对大小来判定与划分，阻值显著减小或与背景值差距较大的区域则为相对富水区[11]。由此可以得出：该采煤工作面底板下105 m 探测范围内有6 处相对低阻异常区，其电阻率值均在90 Ω·m 以下，甚至可低至30 Ω·m。采煤工作面底板下不同深度探测异常区视电阻率空间分布情况如图2。

图2 采煤工作面底板下不同深度探测异常区视电阻率空间分布情况Fig.2 Spatial distribution of apparent resistivity in different depth detection abnormal areas under the floor of coal mining face

依次标定为Dz1～Dz6 三维并行电法相对低阻地质异常区，具体如下：

1）Dz1 低阻异常区。位于采煤工作面靠近运输巷x=0～110 m，y=0～70 m 范围内，形态呈“漏斗”状，向工作面底板下延展较深，深度可达75 m，且异常区之间自上而下有所缩小，但电阻率值由浅入深有变大的趋势，变化幅度也相对较大。

2）Dz2 低阻异常区。位于工作面靠近回风巷x=60～230 m，y=80～150 m 范围内，沿工作面走向展布约170 m，倾向展布约70 m，该相对低阻异常区水平展布较大，向工作面底板下延展至87 m 左右，呈“滴管”状发育，电阻率值向下变高。

3）Dz3 低阻异常区。靠近运输巷位于x=225～315 m，y=0～80 m 范围内，水平展布较大，且向底板下延展很深，呈“酒瓶”状发育，电阻率值变化稳定，介于80～90 Ω·m 之间。

4）Dz4 低阻异常区。处在运输巷x=405～470 m，y=0～80 m 范围内，在工作面底板下方56 m 才开始显现，且呈“钢笔”状发育，大小形态相对较稳定，往底板负法向延伸很大，电阻率值变化较小。

5）Dz5 低阻异常区。靠近回风巷位于采煤工作面x=780～880 m，y=80～150 m 范围内，向下水平横向展布范围增大，变成“梯形”状发育；水平纵向却逐渐减小，且向底板下延展超过105 m，电阻率值变化不大，与Dz3 相似。

6）Dz6 低阻异常区。处于工作面靠近回风巷x=325～400 m，y=105～150 m 范围内，形态呈“倒锥形”状朝下尖灭发育，在z=-56 m 电阻率水平切片处已不可见，电阻率值已变为背景场值，约为95 Ω·m，该异常区低阻值较为稳定。

同时，工作面底板下存在多处相对高阻异常区，电阻率值一般都大于105 Ω·m，甚至高达130 Ω·m左右，也呈现出与低阻异常向下延伸发育的特征。但根据该矿已有的地质资料综合分析，主要存在2个高阻异常区与所分析底板下低阻异常及其推断岩层富水性之间关系紧密，且与该探测范围内的构造有关，分别标记为Gz1 和Gz2。其中Gz1 处于工作面运输巷一侧x=95～205 m，y=0～90 m 之间，发育形态变化不大，面内延展约110 m，主要受F8断层影响；Gz2 高阻区则位于开切眼附近，面内回采方向延展约40 m，向下逐渐变小直至消失，受F6断层影响。此外，在低阻区Dz6 附近出现多条引起相对高阻异常区的分支断层及其伴生小构造，其充、导水条件直接影响到工作面的富水性强弱，工作面底板下24 m处视电阻率水平分布情况如图3。

图3 工作面底板下24 m 处视电阻率水平分布情况Fig.3 Horizontal distribution of apparent resistivity at 24 m below the working face floor

通过分析该工作面底板下不同深度的低阻异常空间竖直展布及连通情况，认为Dz1、Dz2、Dz6 低阻区与下伏栖霞灰岩水力联系不大，其视电阻率值均有往下升高并逐渐消失，推测Dz1、Dz2 与周围发育的断层构造局部导水有关，且开切眼附近的金属支柱也对Dz1 异常有一定的影响；而巷道采动效应与边界效应也会产生伴生裂隙，在充水条件下造成Dz2 和Dz6 呈现低阻异常。Dz3、Dz4、Dz5 相对低阻异常区在z=-69 m 处平面展布范围突变为最大，与栖霞灰岩顶界面发育位置基本一致（-70 m 左右），判断Dz3 是由于采动裂隙与该处之下的溶隙导通，在不考虑受岩性变化的影响下，可能与底板裂隙水存在一定的水力联系；Dz4 则为溶蚀裂隙朝下进一步发展演化转为岩溶空洞，阻值集中于85 Ω·m，与导通下部的岩溶水具有一定的关系；而Dz5 低阻形态表明底板裂隙连通下伏灰岩裂隙直至尖灭，尤其是在岩性变化的位置，灰岩溶隙的大量发育导致富水性增强，范围扩大。因此，Dz1、Dz2、Dz6 低阻区富水性总体上较弱，而Dz3、Dz4、Dz5 较强，是需要重点防范的区域。

3.4 可信度分析与钻探验证

根据LMM 条件下光滑约束模型反演迭代计算的最优化算法，经过10 次循环迭代运算，计算误差逐步减小，相应的拟合相对均方根值（RMSV）下降至10.33 ‰，可决系数R2为0.987 8；有效系数上升至0.97，R2达0.993，拟合程度高、可信度高。据此矿方有针对性地选取重点防范区域Dz2～Dz5 作为钻孔验证的靶区，各实施1 个探水钻孔，其中在Dz2相对低阻异常区所施工的钻孔出水量随深度逐渐减小，与底板裂隙实际发育程度及其含水量有关；Dz3靶区孔深80 m，其涌水量随着钻孔进尺的增加而增大，在底板下70 m 附近水量出现突然增大的现象，由3.42 m3/h 上升到12.24 m3/h，由此表明了该低阻区富水性较强，主要受底板裂隙水的影响，水量陡变点则与栖霞灰岩水有一定的关系；Dz4 低阻区钻孔在钻进至z=-66 m 时，钻杆出现3～4 m 掉钻现象，其出水量可达20.52 m3/h，水质迅速变浑浊，颜色为浑黄色，且伴随有外来泥沙等物质的快速涌入，水压为2.25 MPa，较该工作面底板正常水压高0.25 MPa，据此说明了该相对低阻异常区是隐伏在工作面底板下66 m 的岩溶空洞，与实际探测判定结果较为吻合；Dz5 范围内探水钻孔出水量表现出缓慢增加→快速增大→逐渐稳定等阶段，尤其是穿越灰岩界面时，涌水量陡增直至稳定。测试结果与钻探验证、可信度分析情况基本一致，再次证明了探测分析结果的可靠性较高。