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

1 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院， 合肥 230026

0 引言

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

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

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

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

1 方法原理

1.1 經(jīng)典PSO優(yōu)化算法及其改進(jìn)型算法

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

(1)

(2)

其中，ω是慣性權(quán)重；rand()是(0,1)區(qū)間內(nèi)的隨機(jī)數(shù)；c1,c2為學(xué)習(xí)因子，取值在區(qū)間[0,4]，一般取c1=c2=2；t為迭代次數(shù).

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

林衛(wèi)星和陳炎海(2011)提出的改進(jìn)粒子群(Modified Particle Swarm Optimization, MPSO)算法，針對(duì)慣性權(quán)重系數(shù)早期較大、晚期較小的特點(diǎn)，根據(jù)最佳阻尼比的原理，對(duì)標(biāo)準(zhǔn)PSO算法的慣性權(quán)重系數(shù)做出修改：

(3)

(4)

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

趙玉靜(2011)提出的MPSO，是在公式(3)、(4)的基礎(chǔ)上對(duì)速度更新方程做出了調(diào)整，將第二項(xiàng)當(dāng)前迭代的位置用當(dāng)前迭代的局部最優(yōu)解代替，公式如下：

(5)

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

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

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

(6)

其中α為常數(shù)，取值為[0.1,0.5]區(qū)間內(nèi)，α常取0.1.

1.2 Lévy flight基本理論

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

(7)

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

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

(8)

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

(9)

L(λ)為(8)式的Lévy(λ)，參數(shù)μ和v都服從正態(tài)分布，其中，v服從標(biāo)準(zhǔn)正態(tài)分布：N(0,1)；μ服從：N(0,σ2)；σ由(10)式確定：

(10)

其中，Γ(λ)是參數(shù)為λ的標(biāo)準(zhǔn)Gamma函數(shù).λ一般取值為1.5.

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

(11)

1.3 L-PSO算法基本理論

PSO算法在迭代初期由于粒子的多樣性，使得在早期具有較強(qiáng)的全局搜索能力，但當(dāng)?shù)M(jìn)行，粒子會(huì)向著當(dāng)前最優(yōu)解聚集，使得算法出現(xiàn)早熟現(xiàn)象，尤其是當(dāng)粒子數(shù)不足時(shí)，存在著忽略全局最優(yōu)的小概率事件，一旦迭代進(jìn)行，容易陷入局部最優(yōu).而Lévy flight重尾分布的特點(diǎn)，使得粒子即使在局部最優(yōu)的情況下，也能夠有機(jī)會(huì)進(jìn)行長距離跳躍，正好彌補(bǔ)了PSO算法容易早熟的不足.

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

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

(1) 初始化粒子群：設(shè)定粒子總數(shù)，并對(duì)每個(gè)粒子初始位置xi=(xi1,xi2,…,xin)和初始的速度vi=(vi1,vi2,…,vin)隨機(jī)選擇.設(shè)定最高迭代次數(shù)tmax以及截?cái)嗾`差RMS，每次迭代時(shí)速度v的取值范圍，c1、c2的取值以及慣性權(quán)重系數(shù)ω的取值策略；

(2)計(jì)算粒子群的適應(yīng)度值：根據(jù)適應(yīng)度函數(shù)(即目標(biāo)函數(shù))，計(jì)算每個(gè)粒子的適應(yīng)度值，并根據(jù)適應(yīng)度值更新每個(gè)粒子的個(gè)體最優(yōu)位置pbesti=(pbesti1,pbesti2,…,pbestin)和群體最優(yōu)位置gbest=(gbest1,gbest2,…,gbestn)，在第一次迭代時(shí)，個(gè)體最優(yōu)位置可設(shè)為初始位置，群體最優(yōu)位置是適應(yīng)度最好的粒子位置；

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

圖1展示了L-PSO優(yōu)化算法的基本流程.

圖1 L-PSO優(yōu)化算法流程圖

2 數(shù)值模擬

為了驗(yàn)證L-PSO算法的效果，本文在Matlab2014a平臺(tái)上編程實(shí)現(xiàn)了不同的慣性權(quán)重的取值策略的算法，并設(shè)計(jì)了數(shù)值模擬實(shí)驗(yàn)，取得了理想的結(jié)果.

2.1 常用多極值測試函數(shù)

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

Aekley函數(shù)的表達(dá)式如下：

xi∈[-5.12,+5.12].

(12)

Aekley函數(shù)在定義域內(nèi)具有極密集的極值點(diǎn)，經(jīng)常被用來測試優(yōu)化算法的全局搜索能力.該函數(shù)只存在一個(gè)全局極值點(diǎn)，即(x,y)=(0,0)時(shí)，f(0,0)=0(如圖2).

圖2 Aekley函數(shù)圖像

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

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

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

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

xi∈[-5.12,+5.12]

(13)

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

圖4 Schaffer函數(shù)圖像

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

圖5 6種PSO算法的計(jì)算結(jié)果

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

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

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

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

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

ρ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六種粒子群算法進(jìn)行反演.其中DPSO采用ω=0.99k·rand()/2+0.1的策略，其他參數(shù)統(tǒng)一設(shè)置為：c1=c2=2，粒子數(shù)W=30.迭代次數(shù)上限50次，誤差RMS設(shè)置為3%.為了與真實(shí)模型進(jìn)行更準(zhǔn)確的對(duì)比、方便計(jì)算擬合誤差，反演模型層數(shù)與真實(shí)模型層數(shù)相同，并且其初始模型范圍相同，經(jīng)過5次試驗(yàn)反演之后，平均結(jié)果如表1及圖6所示.

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

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

表1 六種PSO算法反演對(duì)比結(jié)果

對(duì)5次結(jié)果的統(tǒng)計(jì)表明LDPSO算法對(duì)理論模擬地層的反演精度明顯優(yōu)于其他PSO算法，尤其是對(duì)薄層的反演上具有更好的反演效果.值得提及的是，在這10次反演中，LDPSO算法都只用數(shù)次迭代就達(dá)到了誤差精度內(nèi)，展現(xiàn)了LDPSO算法高效的特點(diǎn).

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

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

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.

兩種算法的反演結(jié)果如圖7.

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

對(duì)于KH模型，為更合理的評(píng)估兩種算法和真實(shí)模型的擬合誤差，反演模型層數(shù)與真實(shí)模型采用相同層數(shù)，兩種算法采用相同的最大迭代次數(shù)，并且初始模型相同，反演模型參數(shù)ρ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.由上述反演擬合結(jié)果來看，在相同迭代次數(shù)的情況下，LDPSO算法效果明顯優(yōu)于經(jīng)典SA算法.

圖8是低阻薄層的反演結(jié)果.兩種算法同樣具有相同的最大迭代次數(shù)，迭代次數(shù)上限為100次，反演模型參數(shù)設(shè)置如下：ρ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可以看出，雖然兩種算法對(duì)低阻薄層都有所反映, 但LDPSO算法反演的薄層厚度比SA算法效果更為精確.

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

3 實(shí)測數(shù)據(jù)試驗(yàn)

為驗(yàn)證L-PSO算法的實(shí)際應(yīng)用效果，結(jié)合陜西韓城某煤礦井采空區(qū)瞬變電磁調(diào)查數(shù)據(jù)進(jìn)行測試，并與常用的瞬變電磁反演方法進(jìn)行對(duì)比，分析反演效果.

3.1 地質(zhì)背景

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

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

現(xiàn)場數(shù)據(jù)采集儀器使用重慶璀陸探測技術(shù)有限公司的FCTEM-60瞬變電磁探測系統(tǒng)，該套采集系統(tǒng)采用重疊小線圈收發(fā)裝置，高密度時(shí)序序列數(shù)據(jù)采樣方式，橫向分辨率高.現(xiàn)場施工嚴(yán)格按照地面瞬變電磁施工規(guī)范執(zhí)行.瞬變電磁發(fā)射線框?yàn)榘霃綖?.55 m的多匝圓形線框，發(fā)射電壓為21 V，點(diǎn)距2 m.本地區(qū)實(shí)際完成45條測線，實(shí)際探測瞬變電磁物理點(diǎn)1866個(gè)，測線分布見圖9.

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

3.3 反演

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

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

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

3.4 資料解釋

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

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

圖13 測區(qū)異常體分布與塌陷規(guī)律對(duì)比圖，黃色線段表示推測異常區(qū)在空間上的分布聯(lián)系

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

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

4 結(jié)論與建議

數(shù)值模擬及實(shí)測數(shù)據(jù)的反演案例表明常用的粒子群類算法在引入Lévy flight之后，在全局搜索能力上具有更出色的能力.其具有兩個(gè)優(yōu)點(diǎn)：首先，在保留了經(jīng)典PSO算法收斂的快速性的同時(shí)，提高了算法游走全局的能力，大大增加了算法搜索到全局最優(yōu)的概率；其次，由于Lévy flight重尾的特點(diǎn)，在早期加快了PSO算法的收斂速度，而在極值附近保持了PSO算法的性能.相比蒙特卡羅類算法和模擬退火等優(yōu)化算法，L-PSO算法調(diào)用參數(shù)較少.但值得注意，雖然Lévy flight增加了PSO算法在迭代的任何時(shí)刻都具備探索全局的能力，但面對(duì)單峰函數(shù)，Lévy flight不僅不會(huì)提高PSO算法的性能，反而會(huì)因?yàn)槎嘁徊剿阉鞫晕⒔档褪諗克俣?因此，我們認(rèn)為Lévy flight更適用于多極值的非線性問題.

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

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