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

邱占林，曾東富，吳超凡，郭玉森，劉鑫堯，孫軼群，黃俊南

（1.龍巖學院資源工程學院，福建龍巖 364012；2.福建師范大學地理科學學院，福建福州 350007；3.福建紫金銅業(yè)有限公司，福建龍巖 364204；4.福建省 197 地質(zhì)大隊，福建泉州 362018）

一直以來，福建小煤礦安全、高效生產(chǎn)受到礦井水害的嚴重制約[1]，其主要問題在于采掘之前所采用的鉆探等傳統(tǒng)探水手段的探查能力、精度和范圍有限[2]。物探技術(shù)近年來雖然在福建小煤礦水害探測中得到越來越廣泛的應用，并取得了一定的成效[3-4]，但還是局限于運用常規(guī)直流電法、高密度電阻率法、礦井瞬變電磁法等[5-6]進行巷道頂?shù)装寮坝^前方探測，且一般在單條巷道中采集易受制于巖性等綜合影響的電磁性參數(shù)，不能快速、全面、準確定位，不能有效圈定富水區(qū)的二維平面（2D）或三維空間（3D）形態(tài)，尤其是受底板水害威脅嚴重的地區(qū)，其應用效果更差。傳統(tǒng)直流電法一般采用“跑極”的方式，效率低下，精度不高；高密度電阻率法采用串行單向循環(huán)數(shù)據(jù)采集方式，且只能探測巷道單條電剖面測線的2D 地電場信息[7]，或者精度不高的3D地球物理場信息，局限性較為顯著；而施工便捷且不接地的礦井瞬變電磁法由于對探測工作環(huán)境條件要求較高，其結(jié)果易受到井下鋼軌、金屬支護材料或地電磁場變化等影響[8-9]，數(shù)據(jù)可靠性較差。而且，上述方法還受到數(shù)據(jù)采集時間差或延遲時間的干擾，其相應的校正方法也尚不成熟；在電阻率數(shù)據(jù)反演方面通常也采用阻尼最小二乘法進行迭代計算，但其反演涉及到非線性函數(shù)（或方程）的求解問題，過參數(shù)化、冗余參數(shù)取舍等問題，易導致收斂效果差，并出現(xiàn)相鄰網(wǎng)格單元之間電阻率參數(shù)反演值畸變現(xiàn)象，往往也會造成迭代計算結(jié)果表現(xiàn)出非唯一性（多解性）[10]。相應的存在搜索線性求解方向中止或使損失函數(shù)只得局部極小值的估算，從而引起反演精度降低。

網(wǎng)絡并行電法自運用于福建小煤礦水害探測以來，其快速、瞬時、同步、連續(xù)等數(shù)據(jù)采集優(yōu)點，對于削弱游散電流等干擾因素，提高數(shù)據(jù)信噪比及可信度方面具有較大的優(yōu)勢[3-4,11-12]，但針對福建小礦井富水異常體的位置、分布形態(tài)等參數(shù)擬合反演還存在一定的缺陷。LMM（Levenberg Marquardt Method）是由Levenberg 提出并經(jīng)Marquardt 發(fā)展完善的一種帶阻尼且介于梯度下降法與高斯-牛頓法的非線性或病態(tài)優(yōu)化算法[13-14]。主要利用近似于二階導數(shù)值信息，快速收斂逼近最優(yōu)解，可有效地擬合地質(zhì)異常體。克服了雅克比矩陣奇異或梯度下降法、高斯-牛頓法單一求解難以收斂或偏離局部極小值等問題，使反演逼真響應效果和計算精度得到極大提升，可有效解決網(wǎng)絡并行電法數(shù)據(jù)反演中出現(xiàn)的上述問題，使得LMM 在電法數(shù)據(jù)反演領(lǐng)域得到越來越廣泛的應用。因此，運用LMM 反演的網(wǎng)絡并行電法開展礦井底板富水性探測，獲取待測區(qū)域全空間地電場信息，進而查明其富水性的基本狀況，可為礦井回采工作面底板水害防治提供一定的水文地質(zhì)依據(jù)。

1 LMM 反演理論及模型構(gòu)建

LMM 反演是一種基于高斯-牛頓法和梯度下降法改進并用于解決非線性或病態(tài)方程組的最優(yōu)化方法[13]。是通過引入信賴域的迭代參數(shù)，采用線性函數(shù)來近似（逼近）非線性或病態(tài)函數(shù)，再模仿阻尼最小二乘法求解目標函數(shù)G（m），并在每一步迭代計算過程中尋求合適的變阻尼因子參數(shù)λ 以適用自適應求解的最優(yōu)化，確保目標函數(shù)線性搜索及其數(shù)值是下降的，以此克服雅可比矩陣幾乎全是奇異矩陣或彌補不良條件下高斯-牛頓法不存在逆矩陣、梯度下降法初始值偏離局部極小值太遠等缺陷。而在LMM 反演計算中，每次迭代均要對阻尼因子λ 進行自適應調(diào)整，而后將上述2 種方法混合改進使用，一般為λ 非常小時，反演算法變?yōu)楦倪M的高斯-牛頓法之最優(yōu)步長計算式；反之，則變?yōu)樘荻认陆捣ㄖ顑?yōu)步長計算式，或兩者交叉使用之最優(yōu)步長計算式。只是在局部迭代步長上稍有差別，但均可求取目標函數(shù)極小值時的參數(shù)向量并進行最優(yōu)化處理，可使并行電法電導率σ（電阻率ρ 之倒數(shù)）反演賦值以及優(yōu)化時過參數(shù)化問題等得以合理解決。

常規(guī)網(wǎng)絡并行電法數(shù)據(jù)反演一般屬于多元非線性病態(tài)函數(shù)求解極值范疇[14]，其基本思想是將其轉(zhuǎn)化為大型線性方程組進行數(shù)值迭代計算，通常采用不動點迭代分析法、梯度下降法、（準）牛頓法、高斯-牛頓法等方法，并引入雅克比矩陣J 對單位矩陣實施對角分解，大多是基于均勻半空間地球物理場構(gòu)建目標函數(shù)G（m）。

式中：ρgc為觀測視電阻率值；ρ（m）為正演理論視電阻率值；Wρ為ρ 的協(xié)方差權(quán)重系數(shù)矩陣；α 為光滑約束模型線性搜索因子；m 為正演理論值初始模型向量；△m 為m 基于迭代計算步長的模型增量（修改）向量；H 為粗糙度矩陣；J為雅克比Jacobi 矩陣。

同時依據(jù)LMM 法正則化，可實現(xiàn)反演迭代計算過程中λ 的自適應化。現(xiàn)將式（1）改為：

在信賴域內(nèi)，光滑約束模型線性搜索因子α 一般會對搜索方向及迭代收斂產(chǎn)生梯度下降法和高斯-牛頓法最優(yōu)化處理效果。當α＞0 且系數(shù)矩陣正定時，保證了最速下降效應，從而使收斂速度得到很大提升。當α=λ 時，式（2）則可轉(zhuǎn)換為：

此時光滑約束模型矩陣S 可由Cole-Cole 模型[15]變換替代，即：

式中：ρ0（m）為正演電阻率理論初始值；η 為極化率；ω 為角頻率；τ 為弛豫時間常數(shù)；γ 為頻率因子。

根據(jù)式（3）求取模型增量（修改）向量△m 后，在給定正演理論值初始模型向量m 的基礎(chǔ)上，即可通過式（5）求出新的數(shù)學模型向量。

2 LMM 反演的三維探測技術(shù)原理

2.1 三維探測技術(shù)原理

基于三維電阻率快速層析成像理念（CT）的網(wǎng)絡并行電法探測技術(shù)，除了改進傳統(tǒng)直流電法觀測系統(tǒng)之外，還加入了擬震式采集算法，實現(xiàn)了快速、高效的探測效果。但在數(shù)據(jù)處理方面，還存在求解多元非線性病態(tài)函數(shù)極值過程中冗余數(shù)據(jù)信息取舍及其影響精度等問題。通過引入LMM 反演數(shù)學模型則可有效解決上述問題，形成針對全空間地電場變化的勘探技術(shù)。探測礦井底板富水性的觀測系統(tǒng)通常布置在運輸-回風巷道中，開切眼處由于巷道空間所限不敷設(shè)電極僅用作輔助觀測，工作站則視探測長度而定，由此可產(chǎn)生地下全空間三維CT 透視快速掃描系統(tǒng)并實施數(shù)據(jù)采集與處理。具體為：在雙巷中布設(shè)供電-測量并行電法系統(tǒng)，并將無窮遠極B 安裝在與供電電極A 不在同一巷道的另外1條巷道中，以保證形成類似于坑透法所掃描的扇形透視觀測區(qū)[16]，并行電法三維探測施工布置如圖1。

圖1 并行電法三維探測施工布置Fig.1 Construction layout of three-dimensional detection by parallel electrical method

測量時，供電點A 極（1#～64#）逐點順次供電，而其余電極作為測量電極則通過程控開關(guān)自動實現(xiàn)AM 法或ABM[5,7]法對全空間地電場的實時、同步、快速測量，完成對雙巷空間的全電場3D 電透視觀測。當1 條巷道采集數(shù)據(jù)結(jié)束時，將其中的觀測裝置與B 極對換，再次進行采樣，即可獲得海量擬震式全空間地電場參數(shù)，并通過數(shù)據(jù)拼接以達到對礦井底板的多次覆蓋掃描測量，從中可提取有效地質(zhì)異常體以供后續(xù)解譯。當測量距離較大時，工作站可布設(shè)多個，將測線號（i-j）作為巷道號-工作站號的現(xiàn)場記錄標識。施工中，前1 個工作站最末尾電極跟后1個工作站起始電極位置重合，并依次同向移動另一巷道中的無窮遠極B（i-j），形成實時、同步測量以保證探測精度。

2.2 三維探測數(shù)據(jù)反演

基于LMM 的電阻率三維反演的一般思路則為全空間地電場并行電法數(shù)據(jù)集的迭代數(shù)值計算，通常將全部實測數(shù)據(jù)集轉(zhuǎn)換、拼接、融合，在給定視電阻率矩陣擬合精度ε（閾值）的前提條件下，構(gòu)建地電場正演理論值初始模型m，并將其剖分為可用于賦值的3D 空間網(wǎng)格，各網(wǎng)格單元寬度通常為0.5～1個電極間距，同時求取基于迭代計算步長的模型增量（修改）向量△m，構(gòu)建全空間地電場反演基礎(chǔ)的目標函數(shù)G（m）并對模型修改量求偏導。在此過程中還要考慮粗糙度矩陣H 的3D 分量系數(shù)。同時考慮到地下探測空間地電場介質(zhì)的不均勻性，一般將該空間網(wǎng)格同步劃分轉(zhuǎn)換成1 個大型矩陣組，通過引入單位矩陣I，并將光滑約束模型線性搜索因子α變換為基于LMM 正則化的變阻尼因子λ，然后采用迭代求解對角轉(zhuǎn)換后的高次方程組，獲得上述網(wǎng)格單元反演所要求的正演電導率σ（即電阻率倒數(shù)值ρ（m））等參數(shù)。需要注意的是，在數(shù)據(jù)LMM 反演迭代處理時，需要把笛卡爾坐標系限定在井下3D 空間地電場范圍內(nèi)，從單位矩陣推導出光滑約束模型S對迭代反演計算進行信賴域內(nèi)的校正，并選取ρ0（m）、η、ω、τ、γ 等Cole-Cole 模型參數(shù)來表征反演三維模型，據(jù)此求得數(shù)學模型向量m(n+1)。重復上述迭代計算過程即可完成3D 視電阻率數(shù)據(jù)的反演。而且，由于LMM 充分利用了近似二階導數(shù)項信息以及變阻尼因子的自適應化使反演計算精度得到提升。

3 現(xiàn)場應用

3.1 礦井概況

黃土坑煤礦位于龍巖市東北側(cè)巖山鄉(xiāng)佳山村，地處政和-大埔深大斷裂西側(cè)閩西南聚煤盆地中部含煤條帶北端。以中二疊統(tǒng)童子巖組（P2t）一、三段含煤[17]，其中37、38#、39#煤層主采且穩(wěn)定。開拓方式為平硐+暗斜井，主要采取走向長壁后退式采煤法[18]，采用單體液壓+木支柱支護，年產(chǎn)量15 萬t/a，屬于福建省小型礦井規(guī)模，采掘工作面為“2 采4 掘”，其中+120-37#N 采煤工作面為本次試驗場地。該礦區(qū)構(gòu)造形態(tài)總體呈軸向為NWW-SEE 的背斜構(gòu)造，同時伴生發(fā)育有正、逆斷層共4 條，其中F6逆斷層和F8正斷層對采掘生產(chǎn)有一定的影響。而礦井水文地質(zhì)類型總體為中等-復雜型，區(qū)內(nèi)地下水主要為大氣降水補給，充水水源則以37#煤層頂?shù)装迳皫r裂隙水為主，其次為中二疊統(tǒng)棲霞組（P2q）底板巖溶裂隙水，而37#煤層底板下至棲霞灰?guī)r頂面距離約為70 m。巖石中的各種構(gòu)造裂隙可構(gòu)成礦井的充水通道，并對礦井充水具有普遍意義。

3.2 數(shù)據(jù)采集與處理

為了查明黃土坑煤礦+120-37#N 采煤工作面的底板富水性分布情況，2021 年1 月15 日在礦方地測科的配合下，運用網(wǎng)絡并行電法勘探技術(shù)對該工作面底板視電阻率ρ 展布特征實施三維全空間探測?，F(xiàn)場數(shù)據(jù)采集規(guī)定以開切眼與運輸巷交點位置為笛卡爾坐標系原點（0，0），沿運輸巷指向F8斷層揭露點（173，0）作為x 軸正向，沿開切眼指向F6斷層揭露點（0，132）作為y 軸正向。井下現(xiàn)場采集選用NPEI-DHZI-1 型網(wǎng)絡并行電法儀，應用AM 方式采集，各站依次使用0.5、2.0 s 恒流供電方波測試1次[7]，工作面雙巷則根據(jù)圖1 觀測系統(tǒng)各布置3 個工作站，運輸巷布設(shè)第1～第3 觀測站，回風巷則為第4～第6 觀測站。第3、第6 觀測站敷設(shè)55 各電極（測線長度270 m），其余各站均為64 個電極（共315 m），電極距均為5 m。同時做好復測以保證ρ 數(shù)據(jù)實時、同步采集的準確性與可靠性。

在數(shù)據(jù)處理階段，首先，要做好現(xiàn)場觀測系統(tǒng)坐標系的設(shè)置和全空間地電場數(shù)據(jù)類型的選?。黄浯?，基于LMM 進行電阻率3D 反演以及電流-電壓模塊的轉(zhuǎn)換，同時需要獲取初始值的正演理論模型構(gòu)建與解析、3D 空間網(wǎng)格剖分以及二次場參數(shù)替代等，而光滑約束模型處理及其線性搜索因子轉(zhuǎn)換以尋求變阻尼因子及其自適應化則是該流程的關(guān)鍵，干擾因素分析與選擇性剔除則是數(shù)據(jù)體整體質(zhì)量的保證；最后，從反演模型結(jié)果中提取自適應迭代電阻率數(shù)據(jù)體水平切片并立體構(gòu)建成圖。

3.3 測試結(jié)果分析

對采煤工作面雙巷所采集的全電場數(shù)據(jù)進行全空間電阻率LMM 反演切片提取成像，可較直觀反映采煤工作面底板下不同層位巖層的電性分布情況與變化規(guī)律，可三維呈現(xiàn)相對低阻異常區(qū)的空間分布與總體連通狀況，也可通過提取工作面底板下不同位置視電阻率異常的水平切片及其組合進行表征，據(jù)此再結(jié)合巷道工程揭露的實際地質(zhì)情況，對該工作面底板下的富水性分布規(guī)律及連通性情況進行全盤分析。但底板下通常因巖性的不均一，表現(xiàn)出不同層位視電阻率閾值響應不一致。因此，常采取同一深度視電阻率值的相對大小來判定與劃分，阻值顯著減小或與背景值差距較大的區(qū)域則為相對富水區(qū)[11]。由此可以得出：該采煤工作面底板下105 m 探測范圍內(nèi)有6 處相對低阻異常區(qū)，其電阻率值均在90 Ω·m 以下，甚至可低至30 Ω·m。采煤工作面底板下不同深度探測異常區(qū)視電阻率空間分布情況如圖2。

圖2 采煤工作面底板下不同深度探測異常區(qū)視電阻率空間分布情況Fig.2 Spatial distribution of apparent resistivity in different depth detection abnormal areas under the floor of coal mining face

依次標定為Dz1～Dz6 三維并行電法相對低阻地質(zhì)異常區(qū)，具體如下：

1）Dz1 低阻異常區(qū)。位于采煤工作面靠近運輸巷x=0～110 m，y=0～70 m 范圍內(nèi)，形態(tài)呈“漏斗”狀，向工作面底板下延展較深，深度可達75 m，且異常區(qū)之間自上而下有所縮小，但電阻率值由淺入深有變大的趨勢，變化幅度也相對較大。

2）Dz2 低阻異常區(qū)。位于工作面靠近回風巷x=60～230 m，y=80～150 m 范圍內(nèi)，沿工作面走向展布約170 m，傾向展布約70 m，該相對低阻異常區(qū)水平展布較大，向工作面底板下延展至87 m 左右，呈“滴管”狀發(fā)育，電阻率值向下變高。

3）Dz3 低阻異常區(qū)?？拷\輸巷位于x=225～315 m，y=0～80 m 范圍內(nèi)，水平展布較大，且向底板下延展很深，呈“酒瓶”狀發(fā)育，電阻率值變化穩(wěn)定，介于80～90 Ω·m 之間。

4）Dz4 低阻異常區(qū)。處在運輸巷x=405～470 m，y=0～80 m 范圍內(nèi)，在工作面底板下方56 m 才開始顯現(xiàn)，且呈“鋼筆”狀發(fā)育，大小形態(tài)相對較穩(wěn)定，往底板負法向延伸很大，電阻率值變化較小。

5）Dz5 低阻異常區(qū)?？拷仫L巷位于采煤工作面x=780～880 m，y=80～150 m 范圍內(nèi)，向下水平橫向展布范圍增大，變成“梯形”狀發(fā)育；水平縱向卻逐漸減小，且向底板下延展超過105 m，電阻率值變化不大，與Dz3 相似。

6）Dz6 低阻異常區(qū)。處于工作面靠近回風巷x=325～400 m，y=105～150 m 范圍內(nèi)，形態(tài)呈“倒錐形”狀朝下尖滅發(fā)育，在z=-56 m 電阻率水平切片處已不可見，電阻率值已變?yōu)楸尘皥鲋担s為95 Ω·m，該異常區(qū)低阻值較為穩(wěn)定。

同時，工作面底板下存在多處相對高阻異常區(qū)，電阻率值一般都大于105 Ω·m，甚至高達130 Ω·m左右，也呈現(xiàn)出與低阻異常向下延伸發(fā)育的特征。但根據(jù)該礦已有的地質(zhì)資料綜合分析，主要存在2個高阻異常區(qū)與所分析底板下低阻異常及其推斷巖層富水性之間關(guān)系緊密，且與該探測范圍內(nèi)的構(gòu)造有關(guān)，分別標記為Gz1 和Gz2。其中Gz1 處于工作面運輸巷一側(cè)x=95～205 m，y=0～90 m 之間，發(fā)育形態(tài)變化不大，面內(nèi)延展約110 m，主要受F8斷層影響；Gz2 高阻區(qū)則位于開切眼附近，面內(nèi)回采方向延展約40 m，向下逐漸變小直至消失，受F6斷層影響。此外，在低阻區(qū)Dz6 附近出現(xiàn)多條引起相對高阻異常區(qū)的分支斷層及其伴生小構(gòu)造，其充、導水條件直接影響到工作面的富水性強弱，工作面底板下24 m處視電阻率水平分布情況如圖3。

圖3 工作面底板下24 m 處視電阻率水平分布情況Fig.3 Horizontal distribution of apparent resistivity at 24 m below the working face floor

通過分析該工作面底板下不同深度的低阻異?？臻g豎直展布及連通情況，認為Dz1、Dz2、Dz6 低阻區(qū)與下伏棲霞灰?guī)r水力聯(lián)系不大，其視電阻率值均有往下升高并逐漸消失，推測Dz1、Dz2 與周圍發(fā)育的斷層構(gòu)造局部導水有關(guān)，且開切眼附近的金屬支柱也對Dz1 異常有一定的影響；而巷道采動效應與邊界效應也會產(chǎn)生伴生裂隙，在充水條件下造成Dz2 和Dz6 呈現(xiàn)低阻異常。Dz3、Dz4、Dz5 相對低阻異常區(qū)在z=-69 m 處平面展布范圍突變?yōu)樽畲?，與棲霞灰?guī)r頂界面發(fā)育位置基本一致（-70 m 左右），判斷Dz3 是由于采動裂隙與該處之下的溶隙導通，在不考慮受巖性變化的影響下，可能與底板裂隙水存在一定的水力聯(lián)系；Dz4 則為溶蝕裂隙朝下進一步發(fā)展演化轉(zhuǎn)為巖溶空洞，阻值集中于85 Ω·m，與導通下部的巖溶水具有一定的關(guān)系；而Dz5 低阻形態(tài)表明底板裂隙連通下伏灰?guī)r裂隙直至尖滅，尤其是在巖性變化的位置，灰?guī)r溶隙的大量發(fā)育導致富水性增強，范圍擴大。因此，Dz1、Dz2、Dz6 低阻區(qū)富水性總體上較弱，而Dz3、Dz4、Dz5 較強，是需要重點防范的區(qū)域。

3.4 可信度分析與鉆探驗證

根據(jù)LMM 條件下光滑約束模型反演迭代計算的最優(yōu)化算法，經(jīng)過10 次循環(huán)迭代運算，計算誤差逐步減小，相應的擬合相對均方根值（RMSV）下降至10.33 ‰，可決系數(shù)R2為0.987 8；有效系數(shù)上升至0.97，R2達0.993，擬合程度高、可信度高。據(jù)此礦方有針對性地選取重點防范區(qū)域Dz2～Dz5 作為鉆孔驗證的靶區(qū)，各實施1 個探水鉆孔，其中在Dz2相對低阻異常區(qū)所施工的鉆孔出水量隨深度逐漸減小，與底板裂隙實際發(fā)育程度及其含水量有關(guān)；Dz3靶區(qū)孔深80 m，其涌水量隨著鉆孔進尺的增加而增大，在底板下70 m 附近水量出現(xiàn)突然增大的現(xiàn)象，由3.42 m3/h 上升到12.24 m3/h，由此表明了該低阻區(qū)富水性較強，主要受底板裂隙水的影響，水量陡變點則與棲霞灰?guī)r水有一定的關(guān)系；Dz4 低阻區(qū)鉆孔在鉆進至z=-66 m 時，鉆桿出現(xiàn)3～4 m 掉鉆現(xiàn)象，其出水量可達20.52 m3/h，水質(zhì)迅速變渾濁，顏色為渾黃色，且伴隨有外來泥沙等物質(zhì)的快速涌入，水壓為2.25 MPa，較該工作面底板正常水壓高0.25 MPa，據(jù)此說明了該相對低阻異常區(qū)是隱伏在工作面底板下66 m 的巖溶空洞，與實際探測判定結(jié)果較為吻合；Dz5 范圍內(nèi)探水鉆孔出水量表現(xiàn)出緩慢增加→快速增大→逐漸穩(wěn)定等階段，尤其是穿越灰?guī)r界面時，涌水量陡增直至穩(wěn)定。測試結(jié)果與鉆探驗證、可信度分析情況基本一致，再次證明了探測分析結(jié)果的可靠性較高。

4 結(jié) 語

1）該礦工作面底板下存在6 處低阻異常區(qū)，Dz3、Dz4、Dz5 富水性較強，Dz1、Dz2、Dz6 則相對較弱；總體上具有北部較南部強，淺部較深部強的特征，與鉆探驗證結(jié)果較為吻合。證明網(wǎng)絡并行電法的LMM反演對于礦井底板富水性全空間觀測具有較強的適用性以及較好的地球物理響應。

2）福建小煤礦由于地質(zhì)條件很差，底板下不同層位視電阻率閾值與其富水性之間的定量關(guān)系變得愈加復雜，針對電阻率反演過程中的迭代計算冗余值的合理利用，作為迭代運算終止的擬合精度值定量化，低阻值的高變質(zhì)程度無煙煤對礦井底板富水性異常圈定的干擾等問題，都需要做進一步的分析。