基于粒子群的全空間瞬變電磁二維反演方法研究

黃煒偉

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

0 引言

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

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

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

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

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

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

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

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

具體步驟如下：

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

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

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

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

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

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

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

(1)

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

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

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

(2)

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

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

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

(3)

(4)

3 模型非線(xiàn)性反演

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

以煤礦井下(含巷道)全空間水平層狀地層為基礎(chǔ)建立典型的3層模型，進(jìn)行二維反演模擬，驗(yàn)證反演算法的準(zhǔn)確性，從而為實(shí)際資料的反演提供參考。Q型地電模型如圖2所示，l軸表示模型長(zhǎng)度，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ā)射線(xiàn)框和接收線(xiàn)框以重疊回線(xiàn)的方式水平放置于巷道中央；第3層是底板灰?guī)r層，設(shè)電阻率ρ3為200 Ω·m，層厚h3為100 m[14-16]。模型的長(zhǎng)度和寬度都設(shè)置為200 m，網(wǎng)格剖分時(shí)1個(gè)單位代表1 m。

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

(a) 時(shí)間步長(zhǎng)為800

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

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

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

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

4 工程應(yīng)用

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

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

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

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

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

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

4.2 反演誤差分析

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

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

4.3 地質(zhì)資料解釋

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

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

4.4 驗(yàn)證情況

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

5 結(jié)論

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

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

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