近地表空間梯度瞬變電磁粒子群優(yōu)化反演研究

賴虔林，朱德兵，高 堤，朱 涵

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083；2.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室，湖南 長沙 410083)

1 引 言

瞬變電磁法(Transient Electromagnetic Methods，TEM)是利用不接地回線或接地線源向地下發(fā)送脈沖式一次電磁場，采用回線或接地電極觀測該脈沖電磁場感應(yīng)的地下渦流產(chǎn)生的二次電磁場時空分布，以達(dá)到探測地下地質(zhì)目標(biāo)體的一種電磁感應(yīng)類方法[1,2]。在近地表瞬變電磁勘探中，常規(guī)瞬變電磁法受關(guān)斷時間延遲、線圈自感及互感、一次場及大地二次場混合疊加等因素影響，瞬變早期階段難以觀測到有效信號，造成瞬變電磁法探測的“盲區(qū)”問題[3-7]。為有效解決此問題，朱德兵[8]提出一種瞬變電磁響應(yīng)信號梯度測量方案，該方案能有效消除瞬變電磁一次場互感影響，使瞬變電磁早時段高精度觀測成為可能。周光建[9]基于此方案深入研究空間梯度瞬變電磁法正演理論，得出這種觀測模式顯著提高信噪比和異常分辨能力的結(jié)論?；谏鲜鲅芯炕A(chǔ)，本文開展空間梯度瞬變電磁法一維反演研究。

一維瞬變電磁反演解釋主要有兩種方法。一種是煙圈反演方法，Nabighian[10]于1979年提出把向地下擴散的電流環(huán)近似等效為瞬變電磁響應(yīng)，把隨時間變化的瞬變電磁響應(yīng)曲線轉(zhuǎn)換為電阻率隨深度變化的曲線以實現(xiàn)瞬變電磁的反演；另一種是基于正演公式的最優(yōu)化方法，主要有高斯法、負(fù)梯度法、馬奎特法等。具體過程是先給定初始模型參數(shù)，利用最優(yōu)化方法對瞬變電磁響應(yīng)數(shù)據(jù)進行擬合，進而不斷修正模型參數(shù)，最終達(dá)到實測數(shù)據(jù)和正演數(shù)據(jù)某種意義下的最佳擬合，從而得到最佳反演結(jié)果。Nekut[11]于1987年提出一種中心回線瞬變電磁法的直接反演方法，其運算速度比最小二乘反演更快，能夠估算電導(dǎo)率剖面，但精度不高。Christensen[12]于2002年提出基于自適應(yīng)玻恩近似的快速近似一維反演算法，比普通非線性反演算法快約100倍。張繼令等[13]于2007年將Occam反演技術(shù)應(yīng)用到中心回線瞬變電磁測深數(shù)據(jù)的反演中，實現(xiàn)相對視電阻率的定性或半定量解釋。薛國強等[14]于2008年總結(jié)瞬變電磁反演問題，敘述成像類反演方法主要有時頻等效轉(zhuǎn)換和波場轉(zhuǎn)換兩種方法。李永興等[15]于2010年提出模型交替調(diào)整反演算法，并應(yīng)用于航空瞬變電磁一維反演中，與Zohdy法相比，該方法具有更高精度。李明星等[16]于2014年將粒子群算法應(yīng)用于瞬變電磁反演計算中，取得一定的反演效果。李展輝等[17]于2018年采用正則化Gauss-Newton迭代方法反演瞬變電磁數(shù)據(jù)。楊海燕等[18]于2018年針對圓錐形場源激發(fā)的瞬變電磁響應(yīng)，采用阻尼最小二乘法實現(xiàn)礦井瞬變電磁反演。孫懷鳳等[19]于2019年基于觀測數(shù)據(jù)與模擬數(shù)據(jù)的L1范數(shù)建立目標(biāo)函數(shù)，采用模擬退火非線性全局最優(yōu)化方法實現(xiàn)瞬變電磁一維反演。Guo等[20]于2019年將SDM法應(yīng)用于瞬變電磁資料的反演，通過引入先驗信息，提高了反演精度，加快了反演過程。

本文首先開展空間梯度瞬變電磁法一維正演理論研究，基于Nabighian頻率域響應(yīng)公式并結(jié)合空間梯度觀測模式，采用Gaver-Stehfest數(shù)字濾波算法、正余弦變換，推導(dǎo)出空間梯度瞬變電磁響應(yīng)公式。隨后針對近地表介質(zhì)電性結(jié)構(gòu)復(fù)雜、埋深小等特點，提出一種同時提高全局和局部搜索能力的基于反向?qū)W習(xí)策略的自適應(yīng)混沌粒子群算法(Opposition-Based Learning Adaptive Chaotic Particle Swarm Optimization,OBL-ACPSO)。最后通過Rastrigin、Griewank兩種測試函數(shù)，測試所提方案的性能。利用PSO算法和OBL-ACPSO算法分別對典型理論模型的無噪和含噪瞬變電磁數(shù)據(jù)進行反演，最終評價所提方案的適用性。

表1 模型參數(shù)-中間層電阻率變化

2 一維正演理論

根據(jù)瞬變電磁一次場空間分布的同軸等值性，空間梯度瞬變電磁測量模式利用以發(fā)射線圈為對稱平面的上下平行共軸的兩個相同接收線圈所得信號相減，消除瞬變電磁一次場及收發(fā)線圈互感的影響，實現(xiàn)TEM純二次場響應(yīng)觀測，為近地表電性結(jié)構(gòu)層高精度探測提供了技術(shù)支持。為了分析近地表空間梯度瞬變電磁觀測模式的特點，建立一維層狀介質(zhì)模型。設(shè)地下半空間為n層介質(zhì)組成，1～n層電阻率、厚度分別為ρi、Hi。建立柱坐標(biāo)系，垂直向下Z為正，半徑為a的圓回線位于地表以上h高度，兩個同規(guī)格接收回線關(guān)于發(fā)射回線對稱布置，上、下接收回線分別置于z1,z2位置，空間梯度回線裝置見圖1。

圖1 空間梯度回線裝置示意Fig.1 Schematic diagram of space gradient loop device

由Nabighian[21]研究TEM正演理論可知，當(dāng)圓形發(fā)射回線位于距地面h高度、接收回線置于z位置組成中心回線裝置時，所產(chǎn)生的垂向磁場響應(yīng)為：

(1)

根據(jù)式(1)并結(jié)合空間梯度觀測模式，可推導(dǎo)出空間梯度回線裝置磁場垂直分量頻率域響應(yīng)為：

(2)

式(2)涉及到貝塞爾函數(shù)計算，采用Gaver-Stehfest線性數(shù)字濾波算法[22]求解，對于J1貝塞爾函數(shù)計算，本文選用一階140點的濾波系數(shù)，相關(guān)參數(shù)可見文獻[22]，本文不予具體列出。

計算出磁場垂直分量頻率域響應(yīng)，進行正余弦變換[23]，可求得空間梯度瞬變電磁時間域響應(yīng)，二次場磁場分量對時間的導(dǎo)數(shù)及二次場感應(yīng)電動勢可表示為:

(3)

(4)

為了研究空間梯度回線裝置和重疊回線裝置對近地表介質(zhì)的垂向探測能力，以H型、K型的三層介質(zhì)地電模型為例，通過改變介質(zhì)物性參數(shù)(表1、表2)，正演模擬兩種裝置模式的瞬變電磁響應(yīng)，評價這兩種裝置的探測能力。表1與表2中H1、ρ1分別表示第一層的厚度和電阻率，H2、ρ2分別表示中間層的厚度和電阻率，ρ3表示第三層電阻率?？臻g梯度瞬變電磁裝置參數(shù)設(shè)置為：發(fā)射回線半徑1 m；距地面高度0.5 m；發(fā)射電流10 A；接收回線半徑1 m，下方接收回線置于地面，上方接收回線置于距地面高度1 m位置；發(fā)射及接收線圈的匝數(shù)均為1。裝置如圖1所示。重疊回線裝置置于地面，其它參數(shù)與空間梯度回線裝置一致。

表2 模型參數(shù)——中間層厚度變化

圖2 中間層電阻率變化條件下重疊回線和空間梯度回線異常曲線對比 (cd：重疊回線裝置，td：空間梯度回線裝置)Fig.2 Comparison of abnormal curves of overlapping loop and spatial gradientloop under the condition of resistivity variation of interlayer (cd: Overlapping loop device, td: spatial gradient loop device)

以均勻半空間下的瞬變電磁響應(yīng)作為背景，計算兩組模型在重疊回線裝置和空間梯度回線裝置的正演模擬結(jié)果，研究模型參數(shù)改變時的絕對異常變化和相對異常變化情況，結(jié)果見圖2和圖3，其相對異常計算表達(dá)式如式(5)。

(5)

圖3 中間層厚度變化條件下重疊回線和空間梯度回線異常曲線對比 (cd：重疊回線裝置，td：空間梯度回線裝置)Fig.3 Comparison of abnormal curves of overlapping loop and spatial gradient loop under the condition of Thickness variation of interlayer (cd: overlapping loop device, td: spatial gradient loop device)

式中：δi表示相對異常；Ai表示在物性參數(shù)變化時的瞬變電磁響應(yīng)；A0表示ρ=100 Ω·m時的均勻半空間瞬變電磁響應(yīng)。

由圖2(a)、圖2(b)和圖3(a)、圖3(b)可以看出，不論地電模型是H型還是K型，重疊回線裝置的觀測信號遠(yuǎn)大于梯度回線裝置觀測信號，表明前者在信號幅值上具有絕對優(yōu)勢。由圖2(c)、圖2(d)和圖3(c)、圖3(d)則可以看出，梯度回線裝置測量的相對異常具有明顯優(yōu)勢，尤其是對于1 ms內(nèi)的早時段相對異常，重疊回線裝置相對異常遠(yuǎn)不及梯度觀測模式相對異常。這是因為瞬變早期階段由于梯度回線裝置壓制一次場及收發(fā)線圈互感，而重疊回線裝置信號采集包含瞬變電磁一次場信息，故相對異常信號提取中梯度回線裝置更有優(yōu)勢；瞬變晚期階段由于重疊回線裝置測量一次場衰減完后的總二次場，而空間梯度裝置測量純二次場的差值，其幅值很小，所以重疊回線裝置絕對值測量結(jié)果有優(yōu)勢。

H型和K型地電模型計算結(jié)果表明，空間梯度瞬變電磁觀測模式對低阻介質(zhì)體具有明顯的異常分辨優(yōu)勢，且這種觀測模式尤其適合在早時段(1 ms內(nèi))近地表測量，對電阻率變化反映靈敏度更高，有利于近地表快速、精確探測。同時考慮到近地表介質(zhì)幾何尺度變化大、結(jié)構(gòu)復(fù)雜，其反演難度較大，故需要采取精確穩(wěn)定的反演方案以實現(xiàn)近地表小尺度介質(zhì)梯度瞬變電磁反演。

3 粒子群算法原理

3.1 基本原理

1995年Kennedy和 Eberhart提出了粒子群優(yōu)化算法(Particle Swarm Optimization, PSO)[24]。在PSO算法中，每個“粒子”都會在搜索空間中飛行，并且存在一個速度決定其飛行方向和距離。算法初始化為一群隨機粒子，在每一次迭代中，每個粒子根據(jù)個體歷史最優(yōu)位置pBest和當(dāng)前種群粒子全局最優(yōu)位置gBest來調(diào)整自身的速度和位置，以搜尋全局最優(yōu)位置。根據(jù)式(6)和式(7)粒子更新自身速度和位置。

3.2 改進的粒子群算法

3.2.1 自適應(yīng)慣性權(quán)重系數(shù)

在PSO算法尋優(yōu)中，平衡全局搜索能力和局部搜索能力是有效求解最優(yōu)解的關(guān)鍵。慣性權(quán)重w是PSO算法參數(shù)設(shè)置中至關(guān)重要的參數(shù)，w決定算法的全局與局部尋優(yōu)能力。因此，在PSO算法中須選擇合適w以精確高效地尋找全局最優(yōu)解，文中根據(jù)粒子更新過程中變化的適應(yīng)度值設(shè)置w。具體而言，根據(jù)式(8)確定自適應(yīng)慣性權(quán)重因子(AIWF)。

(8)

式中：wmax、wmin分別表示慣性權(quán)重的最大值、最小值，本文取0.9、0.4；f表示粒子當(dāng)前目標(biāo)函數(shù)值；favg、fmin分別表示當(dāng)前種群平均目標(biāo)函數(shù)值和最小目標(biāo)函數(shù)值。

根據(jù)式(8)，w隨粒子當(dāng)前目標(biāo)函數(shù)值變化而變化。目標(biāo)函數(shù)值小的粒子其慣性權(quán)重變化小，利于在可能的最優(yōu)解附近精細(xì)搜索；目標(biāo)函數(shù)值大的粒子其慣性權(quán)重變化相對大，其優(yōu)勢在于跳出局部極小值進而搜索新的區(qū)域。因此，自適應(yīng)慣性權(quán)重因子為維持種群多樣性和全局收斂能力提供一種良好途徑。

3.2.2 反向?qū)W習(xí)策略

在算法進化后期，粒子多樣性減少，導(dǎo)致收斂速度過快且過于早熟。為了克服這一問題，引入反向?qū)W習(xí)策略。即對PSO算法進化過程中產(chǎn)生的當(dāng)前種群進行性能評價，選擇目標(biāo)函數(shù)值較差的前20 %粒子X，考慮其在搜索空間中對稱位置的粒子X′性能。在多維搜索空間[ai,bi]中，假設(shè)粒子X(x1,x2,...,xn)為n維空間的粒子，則其反向位置X′(x′1,x′2,...,x′n)滿足：

x′i=ai+bi-xi，i=1,2,...,n

(9)

式中：x′i表示某一維粒子xi在搜索空間中對稱位置的粒子；[ai,bi]表示多維搜索空間。

反向?qū)W習(xí)策略同時搜索當(dāng)前位置和反向位置粒子，評估兩者性能，從中選擇較優(yōu)的粒子繼續(xù)更新進化。這有效拓寬了粒子搜索范圍，保持了種群多樣性，增強了PSO算法全局尋優(yōu)能力。

3.2.3 單維混沌局部搜索

式中：x1∈(0,1)且x1?{0.25,0.5,0.75};k表示當(dāng)前混沌搜索次數(shù);M表示最大混沌搜索次數(shù);tpk,d表示第k次搜索第d維局部混沌產(chǎn)生的新解;t表示當(dāng)前粒子群的迭代次數(shù);Tmax表示粒子群的最大迭代次數(shù)。

采用單維混沌局部搜索方式，在PSO算法搜索的當(dāng)前種群最優(yōu)位置gBest附近進行局部搜索，當(dāng)混沌局部搜索的最優(yōu)新解優(yōu)于當(dāng)前位置最優(yōu)解時，用最優(yōu)新解代替當(dāng)前位置最優(yōu)解。通過這種方式可間接提高PSO算法的局部搜索能力，實現(xiàn)小范圍內(nèi)精細(xì)搜索。

3.3 算法測試

為檢驗OBL-ACPSO算法的尋優(yōu)能力，選擇Rastrigin、Griewank函數(shù)作為測試函數(shù)。這兩個函數(shù)相似，有廣泛分布的局部極值，屬于典型的非線性復(fù)雜問題。利用上述兩種測試函數(shù)檢驗OBL-ACPSO算法的性能，兩種函數(shù)的二維表達(dá)式分別是：

上述兩個測試函數(shù)存在全局最小值f(0,0)=0。圖4為兩個測試函數(shù)的三維圖像，可以看出均存在大量局部極小值。

圖4 測試函數(shù)三維圖像Fig.4 3D image of test function

同時選擇GA、PSO、CPSO作為對比算法，其中相關(guān)參數(shù)設(shè)置為：迭代次數(shù)200次；種群數(shù)50；學(xué)習(xí)因子c1=c2=2；慣性權(quán)重系數(shù)wmax=1.2，wmin=0.2。對于GA算法，Pc=0.7，Pm=0.05；對于PSO算法，采用線性遞減的慣性權(quán)重系數(shù)；對于CPSO與OBL-ACPSO兩種算法，采用自適應(yīng)慣性權(quán)重系數(shù)，混沌局部搜索最大步數(shù)為10。為了測試算法的穩(wěn)定性，對上述所有優(yōu)化算法進行20輪運算，并取均值作為最優(yōu)解，其尋優(yōu)結(jié)果對比見表3，相應(yīng)的目標(biāo)函數(shù)曲線見圖5。

表3 優(yōu)化算法對測試函數(shù)的尋優(yōu)結(jié)果對比

圖5 四種優(yōu)化算法對測試函數(shù)的尋優(yōu)能力對比Fig.5 Comparison of four optimization algorithms for test function

由表3可知，對于上述兩種測試函數(shù)，比較四種優(yōu)化算法的最優(yōu)解發(fā)現(xiàn)，OBL-ACPSO算法的尋優(yōu)結(jié)果優(yōu)于GA、PSO和CPSO算法的尋優(yōu)結(jié)果；在20輪計算用時方面，GA算法較其它算法明顯耗時，而對于PSO、CPSO和OBL-ACPSO三種算法的用時比較，T(OBL-ACPSO)>T(CPSO)>T(PSO)，主要原因是CPSO算法在PSO算法的基礎(chǔ)上增加了混沌局部搜索方式，OBL-ACPSO在CPSO算法的基礎(chǔ)上增加了反向?qū)W習(xí)策略。由圖5可知，PSO算法收斂速度最慢，GA算法次之且存在多代未進化現(xiàn)象，CPSO與OBL-ACPSO兩種算法較快，但CPSO算法在40次迭代停止，而OBL-ACPSO算法能夠繼續(xù)迭代尋優(yōu)。經(jīng)過以上對比，OBL-ACPSO算法相比較其它三種算法展現(xiàn)出更好的尋優(yōu)能力，在計算用時相當(dāng)?shù)幕A(chǔ)上，增強了全局和局部搜索能力，提升了全局收斂能力，保證了計算結(jié)果的精度，這為后續(xù)處理復(fù)雜尋優(yōu)問題提供有效解決方案。

4 理論模型實驗

基于上述改進的粒子群優(yōu)化算法，將其應(yīng)用于梯度瞬變電磁一維反演。在反演效果評價中，目標(biāo)函數(shù)的選取至關(guān)重要。因此，針對梯度瞬變電磁數(shù)據(jù)動態(tài)范圍大的特點，本文選取如下目標(biāo)函數(shù)作為梯度瞬變電磁反演的評價準(zhǔn)則：

(14)

設(shè)計實際勘查工作中，常見的三層介質(zhì)地電模型，求解瞬變電磁二次場感應(yīng)電動勢理論值，并將OBL-ACPSO和PSO算法應(yīng)用于空間梯度瞬變電磁一維反演中，通過對比反演結(jié)果，分析OBL-ACPSO的尋優(yōu)能力。在反演過程中，反演模型參數(shù)的搜索范圍設(shè)置為不小于理論參數(shù)值的±60%，其他反演參數(shù)為迭代次數(shù)80次，種群數(shù)50?？臻g梯度瞬變電磁裝置參數(shù)設(shè)置為：發(fā)射回線半徑1 m，距地面高度0.5 m;發(fā)射電流10A;接收回線半徑1 m，下方接收回線置于地面；發(fā)射及接收線圈的匝數(shù)均為1。

4.1 三層介質(zhì)地電模型

圖6～圖9分別為K型、H型、A型、Q型地電模型反演結(jié)果。從圖6～圖9的(a)圖看出：對于三層介質(zhì)地電模型反演，OBL-ACPSO算法比PSO算法效果更好，對層界面刻畫明顯，與理論模型層界面高度擬合；OBL-ACPSO算法感應(yīng)電壓擬合程度較高；圖6～圖9中的(c)圖前60次反演迭代PSO算法下降較為緩慢，而OBL-ACPSO算法在30次迭代迅速下降至趨于穩(wěn)定，迭代效率大幅度提升。

圖6 K型地電模型反演結(jié)果對比Fig.6 Comparison of inversion results of K-type geoelectric model

圖7 H型地電模型反演結(jié)果對比Fig.7 Comparison of inversion results of H-type geoelectric model

圖8 A型地電模型反演結(jié)果對比Fig.8 Comparison of inversion results of A-type geoelectric model

圖9 Q型地電模型反演結(jié)果對比Fig.9 Comparison of inversion results of Q-type geoelectric model

圖10 K型地電模型含5 %高斯噪聲反演結(jié)果對比Fig.10 Comparison of inversion results of K-type geoelectric model with 5 % Gaussian noise

圖11 H型地電模型含5 %高斯噪聲反演結(jié)果對比Fig.11 Comparison of inversion results of H-type geoelectric model with 5 % Gaussian noise

4.2 含噪聲三層介質(zhì)地電模型

瞬變電磁實測數(shù)據(jù)往往含有噪聲干擾，為了測試OBL-ACPSO算法對含噪聲瞬變電磁數(shù)據(jù)的反演性能，對典型三層地電模型正演得到的瞬變電磁感應(yīng)電壓數(shù)據(jù)加入5 %的高斯噪聲，并對含噪聲的瞬變電磁感應(yīng)電壓數(shù)據(jù)進行PSO、OBL-ACPSO算法反演分析。圖10～圖13給出反演結(jié)果，可以看出：圖10～圖13的(a)圖地電模型圖中，OBL-ACPSO反演結(jié)果更接近于理論模型，優(yōu)于PSO反演結(jié)果，表明其抗干擾能力更強。圖10～圖13的(b)圖感應(yīng)電壓擬合曲線中展示了0.07 ms時刻的擬合效果，采用OBL-ACPSO算法反演得到的地電參數(shù)，再通過正演得到感應(yīng)電壓曲線擬合效果優(yōu)于PSO算法；圖10～圖13的(c)圖迭代誤差曲線圖中，兩種算法的最小擬合誤差值趨于一致，但對于達(dá)到該值的最小迭代次數(shù)對比中，OBL-ACPSO算法低于20次，而PSO算法高于60次，說明OBL-ACPSO算法擬合誤差下降速率較快，迭代效率更高。

圖12 A型地電模型含5 %高斯噪聲反演結(jié)果對比Fig.12 Comparison of inversion results of A-type geoelectric model with 5 % Gaussian noise

圖13 Q型地電模型含5 %高斯噪聲反演結(jié)果對比Fig.13 Comparison of inversion results of Q-type geoelectric model with 5 % Gaussian noise

根據(jù)上述無噪與含噪聲瞬變電磁感應(yīng)電壓數(shù)據(jù)反演結(jié)果，OBL-ACPSO算法尋優(yōu)能力較PSO算法好，計算精度更高，反演結(jié)果電阻率曲線形態(tài)與理論模型相符合，尤其表現(xiàn)在對目標(biāo)體層界面的高分辨刻畫，具有一定的分層效果。對于加入一定高斯噪聲情況下的瞬變電磁感應(yīng)電壓數(shù)據(jù)，反演結(jié)果依然和理論模型曲線形態(tài)一致性較好，表現(xiàn)出良好的抗噪聲能力。

5 結(jié) 論

1)空間梯度瞬變電磁法一維正演理論研究表明，空間梯度瞬變電磁觀測模式在1 ms內(nèi)的早時段記錄具有明顯優(yōu)勢，對近地表30 m深度范圍介質(zhì)的地電參數(shù)變化反映較為靈敏。

2)近地表介質(zhì)幾何尺度變化大、結(jié)構(gòu)復(fù)雜，精確反演難度增大。OBL-ACPSO算法明顯增強粒子尋優(yōu)能力，提高了計算精度，抗噪性能突出，反演結(jié)果能夠有效反映地電模型曲線形態(tài)，為近地表空間梯度瞬變電磁數(shù)據(jù)實現(xiàn)小尺度介質(zhì)快速反演解釋提供了新算法。

3)近地表30 m深度范圍內(nèi)介質(zhì)結(jié)構(gòu)地電參數(shù)縱橫向變化大，實際工程中要將算法投入應(yīng)用，還需要對梯度瞬變電磁響應(yīng)接收數(shù)據(jù)進行更精細(xì)地剖析，確認(rèn)1 ms早時段響應(yīng)信號來自于TEM渦流響應(yīng)或僅與介質(zhì)電阻率這一物性參數(shù)關(guān)聯(lián)。