瞬變電磁數據L-PSO反演方法

王書明， 底青云， 夏彤， 任子乾， 宋江濤， 鄒貴安

1 中國地質大學(武漢)地球物理與空間信息學院， 地球內部多尺度成像湖北省重點實驗室， 武漢 430074 2 中國科學院地質與地球物理研究所， 北京 100029 3 中國科學技術大學地球和空間科學學院， 合肥 230026

0 引言

反演問題是地球物理學中的核心問題，是從觀測數據求解地質模型的過程.由于實際地層理論上屬于無限維的空間，然而觀測數據往往是離散且有限的，因此，反演問題的解必定是非唯一的(楊文采，1997).除了觀測資料的局限性之外，場的等效性和觀測數據的誤差也會使得解的非唯一性更加嚴重.

基本的反演方法可分為線性反演和非線性反演兩大類方法.由于大多數地球物理反問題都是非線性的，用傳統的線性反演方法很容易陷入到局部最優解當中，并且線性反演方法非常依賴初始設置的模型.因此，20世紀中期之后，出現了地球物理非線性反演方法，PSO(Particle Swarm Optimization)算法就是其中運用較為廣泛的一種反演方法.

PSO算法在地球物理領域的發展可追溯到20世紀，由印度的Shaw和Srivastava(2007)首先應用于大地電磁(MT)數據的反演中，并做了可行性分析，結果表明該方法可以解決電磁數據的反演問題.PSO方法被引進國內后，主要應用于地震和磁法勘探反演解釋工作中，如在地震波阻抗反演(易遠元等，2007)、磁資料井地聯合反演(張大蓮等，2009)、板狀體磁異常數據反演(吳招才和劉天佑，2009)；在電磁法及瞬變電磁法方面，Olalekan和Di(2018)、易遠元和王家映(2009)、李明星等(2014)、肖敏(2010)、師學明等(2009)、程久龍等(2014)、底青云等(Di et al., 2018)學者發表了一些相關文章.PSO優化算法公式簡單實用，經過多年發展，衍生出很多改進型的算法.但是以上多數工作都是針對于參數c1，c2以及慣性權重ω的改進，并沒有從根本上改善算法前期容易陷入局部最優解的境況.

Levy飛行軌跡是Levy提出，Benoist-Madelbrot進行描述的一種穩定分布(Mandelbrot et al., 1983).飛行是一種隨機游走過程，它的步長是一種連續的重尾分布.Yang和Deb(2009)首次提出了一種稱為布谷鳥算法(CS)的優化算法，在該算法中引入了Lévy flight策略，以此提高長距離搜索能力，避免搜索陷入局部解的情況，取得了良好的效果.侯慧超(2014)針對粒子群算法容易陷入局部最優解的問題，引入萊維飛行的策略，增強了種群的多樣性，也提高了粒子群算法的收斂速度和求解能力.牛海帆等(2016)完成了萊維飛行與粒子群的混合搜索算法.李焰華等(2018)將Lévy flight與烏鴉算法(CSA)相結合，實驗結果表明LCSA算法(Levy Crow Search Algorithm)有效地降低了“烏鴉群”在沒有“領導者”情況下搜索的盲目性.王學武等(2017)利用萊維飛行粒子群算法求解焊接機器人路徑優化問題.因此本文將此策略引入PSO算法中，通過與目前經典的PSO算法及其他優化算法做對比，探究Lévy flight 對PSO算法的影響.

1 方法原理

1.1 經典PSO優化算法及其改進型算法

粒子群優化算法，是受自然中群體行為的啟發發展而來的優化算法，在一個種群中，每個個體都是一個“粒子”，種群(如鳥群)在尋找食物過程中，每個個體都在自己周圍的一個范圍內尋找，通過行為的信息交換，群體會整體向著食物源的方向移動，最后尋找到食物源.在PSO算法中，每個粒子都有個體的位置和速度，對于第i個粒子，定義自己的位置：xi=(xi1,xi2,…,xin)和速度：vi=(vi1,vi2,…,vin).通過三個量確定和更新自身的位置和速度，分別是自身慣性速度vi=(vi1,vi2,…,vin)、個體最優位置pbesti=(pbesti1,pbesti2,…,pbestin)和群體最優位置gbest=(gbest1,gbest2,…,gbestn).并且每個粒子通過(1) (2)式更新自身的速度和位置：

(1)

(2)

其中，ω是慣性權重；rand()是(0,1)區間內的隨機數；c1,c2為學習因子，取值在區間[0,4]，一般取c1=c2=2；t為迭代次數.

由(1)式可以看出，慣性權重系數ω決定了個體的慣性速度vi，較大的慣性權重系數使得迭代時的速度變大，使得粒子“飛”的更遠，有利于算法跳出局部極小，但是較大的ω會不利于算法的收斂，因此為了算法能夠在早期迅速收斂到全局極值附近，應選取較大的ω，在迭代次數增加之后，為了算法能夠迅速收斂，應選取較小的ω.對于參數選取策略，有很多學者做了改進工作.

林衛星和陳炎海(2011)提出的改進粒子群(Modified Particle Swarm Optimization, MPSO)算法，針對慣性權重系數早期較大、晚期較小的特點，根據最佳阻尼比的原理，對標準PSO算法的慣性權重系數做出修改：

(3)

(4)

其中，φ1=c1rand()，φ2=c2rand()，rand()為符合均勻分布的區間(0,1)上的隨機數.

趙玉靜(2011)提出的MPSO，是在公式(3)、(4)的基礎上對速度更新方程做出了調整，將第二項當前迭代的位置用當前迭代的局部最優解代替，公式如下：

(5)

其中，慣性權重系數ω的取值策略與(4)式相同.該方法使得全局最優解與局部最優解影響更大，加快算法收斂速度.

師學明等(2009)的離散粒子群(Discrete Particle Swarm Optimization,DPSO)算法，對慣性因子采取非線性波動遞減的取值策略.使得慣性權重系數ω的振蕩遞減呈現出與模擬退火法類似的退火過程，其速度更新方程等同于(3)式，慣性權重系數取值策略如下：

ω=0.99k·rand()/2+α,

(6)

其中α為常數，取值為[0.1,0.5]區間內，α常取0.1.

1.2 Lévy flight基本理論

Lévy flight是一種非高斯隨機過程(李煜和馬良，2012)，其步長服從Lévy分布.Lévy分布是一種“重尾”的穩定分布，用以描述自然界中多數種群的覓食游走行為特征.同時，Lévy flight也是以偶爾發生長步長跳躍為特點的一類馬爾科夫性質的隨機游走.其位置的更新方程為：

(7)

其中，α為步長因子，常取α=0.01，⊕算子表示點對點的乘積，Lévy(λ)表示服從參數為λ的連續的Lévy分布，其中1<λ<3：

Lévy(λ)～μ=t-λ,

(8)

地球物理反演模型參數通常是有限的，即xi=(xi1,xi2,…,xin)為有限且離散的，所以一般采用Mantegna(1994)提出的離散化公式模擬Lévy分布函數：

(9)

L(λ)為(8)式的Lévy(λ)，參數μ和v都服從正態分布，其中，v服從標準正態分布：N(0,1)；μ服從：N(0,σ2)；σ由(10)式確定：

(10)

其中，Γ(λ)是參數為λ的標準Gamma函數.λ一般取值為1.5.

由以上可得Lévy flight的位置更新方程為：

(11)

1.3 L-PSO算法基本理論

PSO算法在迭代初期由于粒子的多樣性，使得在早期具有較強的全局搜索能力，但當迭代進行，粒子會向著當前最優解聚集，使得算法出現早熟現象，尤其是當粒子數不足時，存在著忽略全局最優的小概率事件，一旦迭代進行，容易陷入局部最優.而Lévy flight重尾分布的特點，使得粒子即使在局部最優的情況下，也能夠有機會進行長距離跳躍，正好彌補了PSO算法容易早熟的不足.

粒子群算法根據不同的參數策略形成了不同變種的算法，其統稱為PSO算法，Lévy flight與具體的算法并無關系，因此可以引入任何一種變體的PSO算法之中，本文將引入Lévy flight之后的任何PSO算法統稱為L-PSO算法，涉及具體算法時，例如：MPSO引入Lévy flight之后稱之為LMPSO；經典PSO算法引入Lévy flight之后簡稱為LPSO算法.

L-PSO算法的具體步驟，實際上是在粒子群進行完一次位置更新后，再進行一步Lévy flight隨機游走，因為Lévy flight偶爾長距離的跳躍的特點，使得部分粒子能夠跳出局部進行搜索.從而增強粒子群搜索全局極值的能力.具體實施如下：

(1) 初始化粒子群：設定粒子總數，并對每個粒子初始位置xi=(xi1,xi2,…,xin)和初始的速度vi=(vi1,vi2,…,vin)隨機選擇.設定最高迭代次數tmax以及截斷誤差RMS，每次迭代時速度v的取值范圍，c1、c2的取值以及慣性權重系數ω的取值策略；

(2)計算粒子群的適應度值：根據適應度函數(即目標函數)，計算每個粒子的適應度值，并根據適應度值更新每個粒子的個體最優位置pbesti=(pbesti1,pbesti2,…,pbestin)和群體最優位置gbest=(gbest1,gbest2,…,gbestn)，在第一次迭代時，個體最優位置可設為初始位置，群體最優位置是適應度最好的粒子位置；

(4)判斷：是否達到誤差要求范圍內或達到最大迭代次數，如否，進行步驟(3)，依此進行，直到滿足終止條件，停止迭代，輸出結果并保存.

圖1展示了L-PSO優化算法的基本流程.

圖1 L-PSO優化算法流程圖

2 數值模擬

為了驗證L-PSO算法的效果，本文在Matlab2014a平臺上編程實現了不同的慣性權重的取值策略的算法，并設計了數值模擬實驗，取得了理想的結果.

2.1 常用多極值測試函數

全局優化算法主要有兩個指標決定其優化能力，一個是在多極值函數中搜索到全局最優即全局搜索能力；第二是收斂能力，即收斂于全局最優解的精度和速度；因此我們選取了兩個典型的函數：Aekley函數與Schaffer函數.

Aekley函數的表達式如下：

xi∈[-5.12,+5.12].

(12)

Aekley函數在定義域內具有極密集的極值點，經常被用來測試優化算法的全局搜索能力.該函數只存在一個全局極值點，即(x,y)=(0,0)時，f(0,0)=0(如圖2).

圖2 Aekley函數圖像

為對比引入Lévy flight后的PSO算法收斂效果，記加入Lévy flight后的經典PSO算法、MPSO(林衛星和陳炎海，2011)以及DPSO算法(師學明等，2009)分別為：LPSO、LMPSO和LDPSO.因為PSO算法也是一種隨機算法，因此我們對每種(引入Lévy flight前和引入Lévy flight后，共6種)算法都進行多次試驗，不設精度，以迭代50次結果平均值作為參考值.

經過10次重復試驗，6種算法都搜索到了全局極值點，說明PSO算法的優良特性.由圖3a對比來看，凡是引入Lévy flight游走的之后的算法，如：PSO與LPSO、MPSO與LMPSO以及DPSO與LDPSO都比未經過Lévy flight的算法得到了更好的結果.其中，師學明(2009)的慣性因子阻尼振蕩遞減取值策略，在加入Lévy flight之后取得了最為理想的收斂效果，其迭代誤差最小，并且也是最快達到收斂的算法，僅用了不到10次就搜索到全局極值點.原因在于，在引入Lévy flight之后，算法在每次迭代都有機會進行“大步長”的搜索，如果搜索效果不好則待在原位置，如果“有幸”得到更好的解，則會加快算法的收斂速度，因此可以看出，Lévy flight具有提高算法收斂速度的特性.

圖3 (a)表示迭代50次，絕對誤差隨迭代次數的變化； (b)可以看出6種算法最終都收斂到全局極值

Schaffer函數是一種常見的用來檢測優化算法全局收斂能力的一類特殊函數，其公式如下：

xi∈[-5.12,+5.12]

(13)

其圖像如圖4所示.該函數特點就是在中心點(x1=x2=0處)存在全局最大值f(x)=1，而在其周圍存在一圈次極值，并且次極值約為0.99，與極值相差甚小，且全局最優解周圍斜率較大，反演結果很容易越過最優極值，落在次級值上.因此該函數適合測試算法的收斂能力以及全局搜索能力.依舊以三種PSO算法對比加入Lévy flight之后的效果.由于其隨機性，進行了大量的測試.

圖4 Schaffer函數圖像

圖5a展示的是6種算法的迭代誤差的遞減，可以看到LDPSO最終誤差最小，對應圖5b，前5種算法最終都收斂到了次極值的圈脊上(兩個綠色的點代表PSO與LPSO，兩個藍色的點代表MPSO與LMPSO，黃色的點代表DPSO，黑色的點代表LDPSO)，只有LDPSO收斂到了全局最優點，精度達到了足以跳出次極值點的精度.這歸功于加入Lévy flight之后，使得算法在搜索到局部解之后，仍然具備一定全局搜索的能力，從而最終找到最優的全局解.可見LDPSO算法具有更好的全局搜索能力.

圖5 6種PSO算法的計算結果

2.2 瞬變電磁L-PSO算法反演

為探究引入Lévy flight策略之后的粒子群算法在瞬變電磁數據反演方面的效果，分別對三種不同慣性參數策略的粒子群算法：PSO、MPSO、DPSO以及其引入Lévy flight之后的改進算法：LPSO、LMPSO、LDPSO，共6種算法進行同類對比.并且與另一種常用的隨機算法——模擬退火算法(SA)進行對比，比較L-PSO算法與經典隨機優化算法的區別.

反演數據由正演模擬得到，模擬裝置采用40 m×40 m的單匝方形中心回線裝置，發射電流為10 A，根據一維層狀模型的解析解(Nabighian, 1992)，首先計算得到頻率域的解(牛之璉, 2007)，然后通過Euler方法進行時-頻轉化，得到時間域瞬變電磁響應(Li et al., 2016; 王鵬飛等, 2019).在Matlab2014a平臺下編程實現一維正演與本文提及的反演方法，處理器為Inter Core i5-6300HQ，主頻為2.3 Hz(雙核四線程).

(1)PSO及其改進型算法的同類對比

設置含有低阻與高阻夾層的7層復雜地層作為理論模型，其地層參數如下：

ρ1=50 Ωm，ρ2=10 Ωm，ρ3=100 Ωm，ρ4=200 Ωm，ρ5=100 Ωm，ρ6=600 Ωm，ρ7=200 Ωm；h1=50 m，h2=5 m，h3=50 m，h4=50 m，h5=100 m，h6=10 m.

分別采用：PSO、LPSO、MPSO、LMPSO、DPSO、LDPSO六種粒子群算法進行反演.其中DPSO采用ω=0.99k·rand()/2+0.1的策略，其他參數統一設置為：c1=c2=2，粒子數W=30.迭代次數上限50次，誤差RMS設置為3%.為了與真實模型進行更準確的對比、方便計算擬合誤差，反演模型層數與真實模型層數相同，并且其初始模型范圍相同，經過5次試驗反演之后，平均結果如表1及圖6所示.

圖6第一行的三張圖為采用不同慣性參數的PSO算法，第二行分別對應其引入Lévy flight之后的改進型算法.從表1的統計結果來看，所有方法最后都得到了收斂，誤差在3%以下.首先從慣性參數的取值策略來看，DPSO的慣性參數策略好于PSO和MPSO.根據圖6的結果，在6種算法中，明顯可以看到LDPSO算法的擬合結果最優，并且相對于3種未采用Lévy flight的PSO算法，在其引入Lévy flight之后，擬合結果都得到了相應的提升，擬合誤差進一步減小.

圖6 六種PSO算法反演得出的模型擬合圖

表1 六種PSO算法反演對比結果

對5次結果的統計表明LDPSO算法對理論模擬地層的反演精度明顯優于其他PSO算法，尤其是對薄層的反演上具有更好的反演效果.值得提及的是，在這10次反演中，LDPSO算法都只用數次迭代就達到了誤差精度內，展現了LDPSO算法高效的特點.

(2)與模擬退火算法對比

粒子群算法的同類算法對比表明LDPSO算法具有更優秀的性能，因此將LDPSO算法與經典模擬退火(SA)算法作對比，以四層地層模型：典型KH型和低阻薄層地層模型為例.設置參數如下：

KH型：ρ1=10 Ωm，ρ2=100 Ωm，ρ3=20 Ωm，ρ4=400 Ωm，h1=50 m，h2=100 m，h3=100 m；

低阻薄層：ρ1=50 Ωm，ρ2=10 Ωm，ρ3=200 Ωm，ρ4=500 Ωm，h1=50 m，h2=10 m，h3=100 m.

兩種算法的反演結果如圖7.

圖7 (a) LDPSO算法反演KH型地層模型的擬合圖，(b) SA算法反演KH型地層模型的擬合圖

對于KH模型，為更合理的評估兩種算法和真實模型的擬合誤差，反演模型層數與真實模型采用相同層數，兩種算法采用相同的最大迭代次數，并且初始模型相同，反演模型參數ρ1為0～100 Ωm，ρ2為10～300 Ωm，ρ3為0～100 Ωm，ρ4為100～1000 Ωm；h1為0～100 Ωm，h2為10～200 Ωm，h3為10～200 Ωm.由上述反演擬合結果來看，在相同迭代次數的情況下，LDPSO算法效果明顯優于經典SA算法.

圖8是低阻薄層的反演結果.兩種算法同樣具有相同的最大迭代次數，迭代次數上限為100次，反演模型參數設置如下：ρ1為0～150 Ωm，ρ2為0～100 Ωm，ρ3為50～600 Ωm，ρ4為100～1000 Ωm；h1為0～100 Ωm，h2為0～100 Ωm，h3為10～200 Ωm.由圖8可以看出，雖然兩種算法對低阻薄層都有所反映, 但LDPSO算法反演的薄層厚度比SA算法效果更為精確.

圖8 (a) LDPSO算法反演低阻薄層地層模型的擬合圖，(b) SA算法反演低阻薄層地層模型的擬合圖

3 實測數據試驗

為驗證L-PSO算法的實際應用效果，結合陜西韓城某煤礦井采空區瞬變電磁調查數據進行測試，并與常用的瞬變電磁反演方法進行對比，分析反演效果.

3.1 地質背景

礦井位于陜西省韓城市東北部，礦井已經開采完畢，主采1#煤層，煤層平均厚度6.1 m，采煤高度2.3～2.5 m，放煤高度3.6～3.8 m，主采煤層海拔約為+230至+240 m.礦區地表為典型黃土高原地貌，溝壑縱橫，平均海拔約為+450 m，地表平均高差約為150 m.淺層以黃土層覆蓋，經受強烈的剝蝕和地表水的長期沖蝕切割，下切很深.中深部因構造應力較弱，地質構造相對比較簡單，除寬緩的褶曲外，巖層傾角平緩，一般10°以下.礦井工作面處于地層走向方向上的寬緩褶曲構造的背斜區域，同時距離邊淺部復雜地質構造帶較遠.在該井邊界內及其附近煤礦，發育中、小型斷層和褶曲構造9條，落差3～9 m不等，造成該地區采空區容易充水，為探明該礦井積水范圍和采空區分布特征，在地面進行了瞬變電磁法探測.

3.2 數據采集與預處理

現場數據采集儀器使用重慶璀陸探測技術有限公司的FCTEM-60瞬變電磁探測系統，該套采集系統采用重疊小線圈收發裝置，高密度時序序列數據采樣方式，橫向分辨率高.現場施工嚴格按照地面瞬變電磁施工規范執行.瞬變電磁發射線框為半徑為1.55 m的多匝圓形線框，發射電壓為21 V，點距2 m.本地區實際完成45條測線，實際探測瞬變電磁物理點1866個，測線分布見圖9.

圖9 工區測線分布圖.由于地表溝壑縱橫，山體起伏較大且陡峭,很難布設貫穿測區的南-北走向的較長測線

3.3 反演

預處理之后，根據瞬變電磁法的電磁場擴散規律，用“煙圈法”進行半定性的解釋.根據“煙圈法”和已掌握的地質資料設置反演初始模型，在反演過程中設置30個粒子數，最大迭代次數為50次，反演采用10層模型，為顯示反演效果，選取了常用的阻尼最小二乘算法進行對比，圖10將阻尼最小二乘算法和LDPSO算法反演的模型參數響應與實測數據的擬合情況同時展示.

圖10 (a) 中左圖為阻尼最小二乘算法反演擬合結果，右圖是LDPSO優化算法反演擬合結果； (b) 是兩種算法在50次迭代中的誤差變化曲線

由擬合結果圖10a和誤差變化折線圖10b的對比可以看出，對于瞬變電磁法早期和晚期數據差異較大的數據，LDPSO算法的反演效果明顯優于阻尼最小二乘算法，并且在50次的迭代過程中，盡管LDPSO在初始位置的解相較于最小二乘方法更差，但其在初始的幾步迭代中得到迅速收斂，到迭代晚期搜尋到了更好的解，再次印證了該方法在瞬變電磁實測資料的一維反演中的適用性.

3.4 資料解釋

我們采用LDPSO算法，對所有測線數據進行了反演，根據測線的實際坐標和海拔標高，繪制了幾條主要的測線剖面在測區的實際位置(如圖11)；并根據反演電阻率的相對值，估計采空區電阻率范圍勾勒出測區的三維異常圖(如圖12)，以及反演目標層范圍內的異常體平面分布(如圖13).

圖11 左圖表示工區內穿過異常體的幾條測線反演的電阻率剖面在空間的位置，右圖為這幾條測線的俯視圖

圖13 測區異常體分布與塌陷規律對比圖，黃色線段表示推測異常區在空間上的分布聯系

從反演電阻率剖面切片圖11可以看出，通過異常體的幾條主要測線，反演的異常體均在煤層所在深度，并且各個測線之間的反演結果在平面上能夠相互對應，具有一致性.從圖12的異常體分布范圍可以看出，低阻異常體呈現條帶狀分布，符合采空區的空間分布，且異常體的高度約為15 m，根據煤礦采空區一般的垮落規律，采空區冒落帶高度一般為采煤高度的3倍左右.因此推測低阻異常體是由于采空區積水造成的，而反演目標深度范圍內的高阻異常體，推斷為由塌陷引起的主要空洞、裂隙、破碎帶非積水區.

通過掌握的已有地質資料結合瞬變電磁反演結果，我們認為，由于該工作面已經停產了相當的一段時間，地表溝壑縱橫，且土質疏松，由于采空區塌陷影響和原本的裂隙發育，形成滲水通道，易于大氣降水及上部含水層中水的補給和儲存，造成了采空區積水.其次，根據圖13的異常體分布范圍，可以看出，工區東西兩側的采空區范圍明顯比中間區域的采空區范圍大，這是由于單一工作面塌陷時，靠近中心位置的采空區越塌實，不易形成積水空間，而越靠近兩側時采空區塌陷情況復雜，不能完全塌實，因此形成了可能的積水空間.通過本次探測的結果，該區的塌實為“畸變”的OX型結構，從整體上符合采空區塌陷的一般規律.

4 結論與建議

數值模擬及實測數據的反演案例表明常用的粒子群類算法在引入Lévy flight之后，在全局搜索能力上具有更出色的能力.其具有兩個優點：首先，在保留了經典PSO算法收斂的快速性的同時，提高了算法游走全局的能力，大大增加了算法搜索到全局最優的概率；其次，由于Lévy flight重尾的特點，在早期加快了PSO算法的收斂速度，而在極值附近保持了PSO算法的性能.相比蒙特卡羅類算法和模擬退火等優化算法，L-PSO算法調用參數較少.但值得注意，雖然Lévy flight增加了PSO算法在迭代的任何時刻都具備探索全局的能力，但面對單峰函數，Lévy flight不僅不會提高PSO算法的性能，反而會因為多一步搜索而略微降低收斂速度.因此，我們認為Lévy flight更適用于多極值的非線性問題.

L-PSO算法相較于其他改進型的PSO算法，其多出了一步游走的步驟，因此多調用了一次正演，使得速度方面有所降低，對于瞬變一維反演來說，時間增加并不明顯，但是對于二、三維反演，瞬變正演過程耗時巨大，每多一次正演對反演的效率都造成很大影響，因此，如何將L-PSO算法應用于高維瞬變電磁反演是下一步需要研究的問題.

致謝感謝重慶璀陸探測技術有限公司提供的儀器，感謝審稿人給出的建設性意見和建議.