黃煒偉

(山西潞安環(huán)保能源開發(fā)股份有限公司 漳村煤礦， 山西 長治 046032)

0 引言

瞬變電磁法是指利用發(fā)射線圈向地下發(fā)射一次脈沖磁場，通過接收線圈接收地下的感應(yīng)電流產(chǎn)生的二次場，從而探測介質(zhì)電阻率的一種方法。瞬變電磁法具有施工效率高、對低阻體敏感等優(yōu)點(diǎn)，適用于井下掘進(jìn)巷道及回采工作面頂?shù)装搴蕴綔y、掘進(jìn)巷道超前探測及隧道地質(zhì)構(gòu)造富水性探測等。

在電磁學(xué)中，正演是指將已知的電荷屬性換算成電場分布，并計(jì)算得到空間某點(diǎn)的電場屬性；反演是指通過電場分布得到一種電荷分布模型或近似模型，得到實(shí)際地層的地電特性。通過一維反演只能得到深度信息，而通過二維反演則可以得到橫向和縱向地電特性，對異常區(qū)的圈定由1個(gè)方向變?yōu)?個(gè)方向，從而提高探測精度。煤礦井下瞬變電磁的探測范圍是全空間的，在全空間瞬變電磁法數(shù)據(jù)處理過程中很少加入反演計(jì)算，常用的線性反演難以滿足生產(chǎn)需求[1-2]。然而，非線性反演問題復(fù)雜，現(xiàn)有研究多停留在一維反演階段，很少涉及二維非線性反演。為了滿足生產(chǎn)需求，提高反演精度，有必要研究二維非線性反演方法。在通過全空間瞬變電磁法采集的原始數(shù)據(jù)處理過程中加入二維反演計(jì)算，對提高探測精度和分辨率具有重要意義。

粒子群算法是一種非線性反演計(jì)算方法，在反演過程中可以適應(yīng)地電特性的變化，因此，學(xué)者們將其應(yīng)用到瞬變電磁反演計(jì)算中并進(jìn)行了改進(jìn)研究。文獻(xiàn)[3-4]對全空間瞬變電磁法進(jìn)行了粒子群優(yōu)化反演研究，改進(jìn)了粒子群尋優(yōu)算法，但只進(jìn)行了一維反演計(jì)算，未解決體積效應(yīng)問題。文獻(xiàn)[5]主要從算法方面對粒子群反演方法進(jìn)行改進(jìn)，未結(jié)合實(shí)際應(yīng)用數(shù)據(jù)進(jìn)行分析驗(yàn)證。文獻(xiàn)[6]最早將一維非線性反演技術(shù)引入到地球物理探測中，開拓了粒子群反演的先河，但是尋優(yōu)算法還有待改進(jìn)。文獻(xiàn)[7-8]對地球物理反演基礎(chǔ)理論進(jìn)行了研究，但在全空間瞬變電磁方面未進(jìn)行理論研究和實(shí)際應(yīng)用。文獻(xiàn)[9-10]引入粒子群算法并進(jìn)行了系統(tǒng)研究，為全空間瞬變電磁反演計(jì)算研究提供了基礎(chǔ)。

隨著煤礦開采深度的加大，水害威脅越來越嚴(yán)重，針對煤礦實(shí)際采掘情況的復(fù)雜性，為了提高探測精度，為鉆探提供有效靶點(diǎn)，有必要開展全空間瞬變電磁二維反演方法的研究，為礦井水害防治提供可靠的地質(zhì)資料。本文在一維反演和二維正演研究的基礎(chǔ)上，結(jié)合煤礦實(shí)際探測精度要求，提出了基于粒子群的全空間瞬變電磁二維反演方法，通過粒子群算法對電阻率和地層厚度參數(shù)進(jìn)行尋優(yōu)，以提高反演精度，為全空間瞬變電磁資料處理解釋奠定基礎(chǔ)。

1 全空間瞬變電磁二維反演方法

在全空間的前提下，二維反演和一維反演都是根據(jù)瞬變電磁儀探測到的歸一化感應(yīng)電位值和與之對應(yīng)的時(shí)間進(jìn)行反演計(jì)算，最終得到探測目標(biāo)體的電性參數(shù)。由于實(shí)際地層的電性結(jié)構(gòu)是復(fù)雜的三維結(jié)構(gòu)，僅僅依靠一維反演結(jié)果來反映地層的地電特性誤差可能較大，需通過二維反演來反映地層的真實(shí)情況。

基于粒子群的全空間瞬變電磁二維反演方法流程如圖1所示。

圖1 基于粒子群的全空間瞬變電磁二維反演方法流程Fig.1 Flow of two-dimensional inversion method of transient electromagnetic in whole-space based on particle swarm

具體步驟如下：

(1) 基于實(shí)測資料，通過粒子群算法進(jìn)行全空間瞬變電磁一維反演。以一維反演結(jié)果中的探測深度值及其對應(yīng)的視電阻率值作為初始模型參數(shù)，建立二維正演模型。

(2) 選用等間距分布的網(wǎng)格對二維正演模型進(jìn)行網(wǎng)格化，每個(gè)網(wǎng)格的大小為1 m×1 m。對網(wǎng)格化的二維正演模型進(jìn)行時(shí)域有限差分模擬，得到不同時(shí)刻各網(wǎng)格節(jié)點(diǎn)處的磁場強(qiáng)度值。

(3) 按照原始數(shù)據(jù)采集時(shí)的時(shí)間門所對應(yīng)的時(shí)間和一維反演結(jié)果中的探測深度來抽取正演模擬后的磁場強(qiáng)度值，使得抽取的磁場強(qiáng)度值與實(shí)測資料中的磁場強(qiáng)度值具有時(shí)間和深度的對應(yīng)關(guān)系。

(4) 基于最小二乘反演計(jì)算方法，對抽取的磁場強(qiáng)度值與實(shí)測資料中的磁場強(qiáng)度值進(jìn)行反演計(jì)算。通過調(diào)整二維正演模型中的地電參數(shù)，使得最小二乘反演后的擬合誤差處在適合二維反演的誤差區(qū)間。

(5) 對經(jīng)過最小二乘法反演后得到的磁場強(qiáng)度值和與之對應(yīng)的時(shí)間值按測點(diǎn)號順序依次進(jìn)行二維反演[11-13]。

模型M(q)的參數(shù)包括K層地層的電阻率和K-2層地層的厚度，給定參數(shù)初始值x0，即一維反演得到的視電阻率和與之對應(yīng)的探測深度，令迭代過程中模型參數(shù)的修正方程為[10-11,13]

AT(f(x0)+AΔx0)=0

(1)

式中：A為雅克比距陣；f(x0)為函數(shù)F泰勒展開的一階項(xiàng)；Δx0為迭代過程中模型參數(shù)的修正值。

由式(1)可求得模型參數(shù)修正值為

Δx0=-(ATA)-1ATf(x0)

(2)

模型的近似解為x1=x0+Δx0，將x1作為下次迭代的初始值，進(jìn)行多次迭代，直至擬合誤差滿足精度要求為止。

反演過程中模型參數(shù)修正對反演效果起著至關(guān)重要的作用。給定一組參數(shù)，可得到一個(gè)對應(yīng)的模型參數(shù)修正值，通過參數(shù)調(diào)整，可得到若干模型參數(shù)修正值。將各模型參數(shù)修正值作為迭代的初始值，經(jīng)過多次迭代后，若結(jié)果擬合誤差符合精度要求，則停止迭代，該結(jié)果為最終的反演結(jié)果；若擬合誤差不滿足精度要求，則繼續(xù)迭代或更換模型參數(shù)修正值后繼續(xù)迭代，直到擬合誤差滿足精度要求為止[11]。

2 基于粒子群算法的參數(shù)尋優(yōu)

(3)

(4)

3 模型非線性反演

3.1 模型瞬變電磁響應(yīng)規(guī)律分析

以煤礦井下(含巷道)全空間水平層狀地層為基礎(chǔ)建立典型的3層模型，進(jìn)行二維反演模擬，驗(yàn)證反演算法的準(zhǔn)確性，從而為實(shí)際資料的反演提供參考。Q型地電模型如圖2所示，l軸表示模型長度，z軸表示模型寬度。

圖2 Q型地電模型Fig.2 Q-type geoelectric model

模型共分為3層：第1層是頂板砂巖層，設(shè)電阻率ρ1為500 Ω·m，層厚h1為80 m；第2層是煤層，設(shè)電阻率ρ2為400 Ω·m，層厚h2為20 m，巷道位于該層最下端，發(fā)射線框和接收線框以重疊回線的方式水平放置于巷道中央；第3層是底板灰?guī)r層，設(shè)電阻率ρ3為200 Ω·m，層厚h3為100 m[14-16]。模型的長度和寬度都設(shè)置為200 m，網(wǎng)格剖分時(shí)1個(gè)單位代表1 m。

激發(fā)源位于l=100 m,z=0處，通過發(fā)射線框發(fā)射一次脈沖磁場，并分析磁場在模型地層中的傳播規(guī)律[15-17]。時(shí)間步長為800，1 600時(shí)Q型地電模型瞬變電磁場擴(kuò)散切片如圖3所示，其中y為與層狀地層走向垂直方向距巷道頂板的距離，在巷道頂板處y=0。從圖3可看出，自y=0斷面向y=80 m斷面，磁場強(qiáng)度值及擴(kuò)散范圍均逐漸變??；z>100 m范圍內(nèi)的磁場擴(kuò)散速率大于z<100 m范圍內(nèi)的磁場擴(kuò)散速率，反映了z>100 m范圍內(nèi)的電阻率高于z<100 m范圍內(nèi)的電阻率，與地層模型中設(shè)置的電阻率相對應(yīng)。綜合分析可知，隨著時(shí)間的增加，磁場擴(kuò)散范圍明顯增大；因?yàn)榇艌鲈跀U(kuò)散過程中不斷衰減，所以靠近激發(fā)源處磁場強(qiáng)度值較大，遠(yuǎn)離激發(fā)源處磁場強(qiáng)度值較小。

(a) 時(shí)間步長為800

(b) 時(shí)間步長為1 600

3.2 模型粒子群反演誤差分析

Q型地電模型瞬變電磁反演計(jì)算擬合誤差如圖4所示。迭代2次時(shí)擬合誤差為1×10-2，迭代10次時(shí)擬合誤差為1×10-3，迭代50次時(shí)擬合誤差為1.3×10-4。對于前20次迭代，隨著迭代次數(shù)的增大，擬合誤差下降速率較大；對于第21次到第50次迭代，隨著迭代次數(shù)的增大，擬合誤差下降速率明顯減小，整體反演精度相對較高。

圖4 Q型地電模型瞬變電磁反演計(jì)算擬合誤差Fig.4 Fitting error of transient electromagnetic inversion calculation of Q-type geoelectric model

4 工程應(yīng)用

選取某礦-550 m水平西翼總回巷道底板瞬變電磁探測資料進(jìn)行反演計(jì)算，由巷道東段向西段依次采集數(shù)據(jù)，探測起點(diǎn)定為0號測點(diǎn)，測點(diǎn)間距為10 m，探測終點(diǎn)為200號測點(diǎn)，測點(diǎn)總數(shù)為21個(gè)，總探測長度為200 m[11]。

4.1 瞬變電磁場擴(kuò)散規(guī)律

時(shí)間步長為2 500時(shí)的瞬變電磁場擴(kuò)散圖如圖5所示，激發(fā)源位于l=100 m,z=0處，z軸坐標(biāo)值的負(fù)號表示底板以下探測深度，正號表示巷道頂板以上探測深度；l軸坐標(biāo)值與實(shí)測數(shù)據(jù)的測點(diǎn)號相對應(yīng)，例如，l=10 m時(shí)，對應(yīng)的測點(diǎn)號為10號。

圖5 時(shí)間步長為2 500時(shí)的瞬變電磁場擴(kuò)散圖Fig.5 Transient electromagnetic field diffusion diagram when the time step is 2 500

從圖5可看出，由于存在巷道空間，瞬變電磁場擴(kuò)散到巷道邊界時(shí)發(fā)生畸變，隨著時(shí)間的推移，巷道邊界對瞬變電磁場擴(kuò)散造成的影響逐漸減弱，瞬變電磁場開始向巷道頂板以上地層和底板以下地層擴(kuò)散。因?yàn)橄锏理敯逡陨系貙記]有異常體，磁場均勻擴(kuò)散，所以瞬變電磁場強(qiáng)度逐漸減弱。底板以下地層中存在低阻異常體，瞬變電磁場在低阻體中擴(kuò)散速度較慢，因此，瞬變電磁場在低阻異常體邊界處形成新的激發(fā)擴(kuò)散源。

根據(jù)圖5可確定低阻異常體的邊界，只是140號測點(diǎn)到200號測點(diǎn)探測深度為底板以下65～100 m時(shí)，低阻體處于模型的右下角，此時(shí)該區(qū)域瞬變電磁場強(qiáng)度較弱，使得低阻異常體的邊界不是很明顯。

4.2 反演誤差分析

磁場強(qiáng)度擬合誤差曲線如圖6所示。在前11次迭代過程中，隨著迭代次數(shù)增加，擬合誤差迅速下降；到第11次迭代時(shí)，擬合誤差為0.8%，小于擬合誤差精度要求，擬合效果較好；之后，隨著迭代次數(shù)的增加，擬合誤差下降趨勢不明顯，曲線逐漸趨于平緩；到第30次迭代時(shí)，擬合誤差僅為0.44%。當(dāng)?shù)螖?shù)增大到一定程度時(shí)，如果繼續(xù)增大迭代次數(shù)，則擬合誤差減小幅度很小，但是反演計(jì)算所用時(shí)間會成倍增加，反而會降低計(jì)算效率。因此，在反演過程中需要選取合適的迭代次數(shù)，使得反演后的擬合誤差滿足精度要求。本次反演計(jì)算所用時(shí)間為76.592 232 s，計(jì)算效率較高。

圖6 磁場強(qiáng)度擬合誤差曲線Fig.6 Fitting error curve of magnetic field intensity

4.3 地質(zhì)資料解釋

全空間瞬變電磁二維反演視電阻率斷面如圖7所示。圖7可以反映巷道頂?shù)装鍘r層的電性特征，巷道頂板以上30～100 m及巷道底板以下距離巷道底部0～50 m探測范圍內(nèi)巖層的視電阻率值相對較高，判斷為巖層不含水，這主要反映了底板以砂泥巖為主的巖性特征。巷道頂板以上20～30 m及巷道底板以下50～100 m范圍內(nèi)巖層的視電阻率值相對較低，低阻異常區(qū)共劃分為3處：① 測點(diǎn)號40～70，115～195，頂板以上20～30 m探測范圍；② 測點(diǎn)號0～145，底板以下55～100 m探測范圍；③ 測點(diǎn)號160～200，底板以下62～100 m探測范圍。上述3個(gè)低阻異常區(qū)呈層狀分布，判斷為巖層相對弱含水，這主要反映了底板深部砂巖含水層的特征。

圖7 全空間瞬變電磁二維反演視電阻率斷面Fig.7 Cross section of apparent resistivity inversion of two-dimensional inversion of transient electromagnetic in whole-space

4.4 驗(yàn)證情況

實(shí)測數(shù)據(jù)處理過程中加入了二維反演計(jì)算，處理結(jié)果解釋后，結(jié)合煤礦井下現(xiàn)有巷道條件，在煤層頂板設(shè)計(jì)測線50，100，130，180 m處各布置1個(gè)鉆孔，在煤層底板設(shè)計(jì)測線20，120，150 m處各布置1個(gè)鉆孔，進(jìn)行鉆探驗(yàn)證，鉆孔深度設(shè)計(jì)為80 m。頂板鉆孔基本無出水情況，只有反演結(jié)果低阻區(qū)域的2個(gè)鉆孔有滴漏水情況；底板鉆孔出水量相對較大，尤其是測線20 m和120 m處的鉆孔出水量為2 m3/h；測線150 m處的鉆孔出水量為0.1 m3/h。驗(yàn)證結(jié)果表明，反演后的結(jié)果明顯提高了探測精度和分辨率，對含水低阻異常區(qū)范圍圈定準(zhǔn)確，二維反演結(jié)果與實(shí)際驗(yàn)證情況基本吻合，反演后精度明顯提高。

5 結(jié)論

(1) 介紹了全空間瞬變電磁二維反演方法的原理、流程及具體步驟，并通過粒子群算法對電阻率和地層厚度參數(shù)進(jìn)行尋優(yōu)。

(2) 采用Q型地電斷面模型分析得出瞬變電磁場擴(kuò)散規(guī)律：隨著時(shí)間的增加，磁場擴(kuò)散范圍明顯增大；因?yàn)榇艌鲈跀U(kuò)散過程中不斷衰減，所以靠近激發(fā)源處磁場強(qiáng)度值較大，遠(yuǎn)離激發(fā)源處磁場強(qiáng)度值較小。

(3) 在理論研究的基礎(chǔ)上進(jìn)行工程應(yīng)用，對實(shí)際資料進(jìn)行二維反演，分析了瞬變電磁場的擴(kuò)散規(guī)律和擬合誤差，解釋了視電阻率探測成果，并進(jìn)行了鉆探驗(yàn)證。應(yīng)用結(jié)果表明，全空間瞬變電磁二維反演方法是可行的，反演結(jié)果中的低阻異常體具有成層特性及一定的連通性，更能反映層狀地層以砂泥巖為主的巖性結(jié)構(gòu)特征及含水特性，提高了探測精度和分辨率。