基于改進粒子群算法的大地電磁反演

李麗麗,, 李長偉,2*, 程勃,2, 陳漢波, 呂玉增,2, 熊彬,2, 張媛, 黃楊

(1. 桂林理工大學(xué)地球科學(xué)學(xué)院, 桂林 541000; 2. 廣西隱伏金屬礦產(chǎn)勘查重點實驗室, 桂林 541000)

大地電磁(magnetotelluric sound, MT)反演方法大體可分為線性和非線性反演[1]。線性反演迭代收斂快,主要有高斯牛頓法[2]、共軛梯度法[3]、馬夸特法[4]等,但是這些方法存在反演結(jié)果依賴于初始模型,易陷入局部最優(yōu)解的缺點。非線性反演方法既能克服線性反演的缺點,又能有效地避免算法陷入局部最優(yōu)[5]。因此,非線性反演方法一直是學(xué)者們研究的熱點。常見的非線性反演方法包括粒子群算法(particle swarm optimization, PSO)、蒙特卡洛、遺傳算法,神經(jīng)網(wǎng)絡(luò)算法等。由于粒子群算法原理簡明、不依賴于初始模型、具有全局尋優(yōu)、收斂速度快等優(yōu)點,因此廣泛應(yīng)用于電力、地球物理、機器學(xué)習(xí)等各種領(lǐng)域中[6-9]。

PSO算法是由Shi等[10]在1995年通過模擬鳥覓食行為而發(fā)展起來的全空間搜索的算法,由Shaw等[11]最先應(yīng)用在MT數(shù)據(jù)的反演中,從可行性分析來看,該方法可以有效解決電磁數(shù)據(jù)的反演問題。PSO算法目前主要應(yīng)用于磁法、地震、測井、電法等方面。

趙平起等[12]在效差分算子的基礎(chǔ)上,將BFP算法中的趨化、復(fù)制和擴散3個步驟加入PSO算法中進行反演,彈性波在不同模型上的數(shù)值模擬結(jié)果表明,改進算法不僅具有更高的精度,而且能夠有效壓制數(shù)值頻散。陳杰等[13]將GA算法和PSO算法同時用于鐵路隧道實測數(shù)據(jù)進行反演,結(jié)果表明兩種方法對理論模型模擬的觀測數(shù)據(jù)均具有較高的反演精度,且PSO算法的迭代收斂速度更快,但在加入噪聲后的反演中GA算法卻具有更強的抗干擾能力。王書明等[14]采用Levy flight 飛行策略結(jié)合了短距離搜索與偶爾長距離游走的特點,用于求解瞬變電磁二維反演問題,反演結(jié)果表明改進算法在保證算法收斂效率的同時提高了算法收斂全局的能力。李曦等[15]在優(yōu)化支持向量機的基礎(chǔ)上采用PSO算法進行巖性識別,不僅正確識別了巖性,而且提高了測井儲層解釋精度。Abril等[16]將帶有遺傳算子的PSO算法(EMPSO)用來解決二維電阻率成像問題,并通過將純MPI和混合MPI-OpenMPI并行來解決時間問題,真實數(shù)據(jù)反演效果良好,與沒有并行的相比,提高了計算時間。PSO算法經(jīng)過多年地發(fā)展,已經(jīng)衍生出很多優(yōu)化的算法。但是大多數(shù)學(xué)者都是針對算法參數(shù)學(xué)習(xí)因子, 和慣性權(quán)重的優(yōu)化,并未從根本上改善算法前期容易陷入局部極值的情況。

針對PSO算法維度越高,越容易陷入局部極值,大地電磁傳統(tǒng)反演方法依賴初始模型的缺點,現(xiàn)提出一種改進的粒子群算法(cooperation adaptivity particle swarm optimization, CA-PSO),來進行大地電磁反演模擬,并在幾個測試函數(shù)和物理模型中進行了測試,與其他改進的粒子群算法進行了對比,討論該算法的性能。CA-PSO算法一維反演的目標(biāo)函數(shù)是利用觀測視電阻率和正演得到的視電阻率構(gòu)建,二維反演的目標(biāo)函數(shù)在一維的基礎(chǔ)上添加核函數(shù)以及先驗信息的限制。為了改善PSO算法容易陷入局部極值,增加粒子的活性,與標(biāo)準(zhǔn)粒子群算法相比,改進算法的迭代方式多了局部進化,將靜態(tài)慣性權(quán)重替換為動態(tài)慣性權(quán)重,并且增加了收縮因子參數(shù)。

1 改進的粒子群算法(CA-PSO)

1.1 粒子群算法

粒子群算法起源于對鳥群運動行為的研究。每只鳥在某處位置能夠找到食物的可能通過適應(yīng)度來刻畫,每只鳥都能記著它的覓食地點,并找到最好的(局部最優(yōu),相當(dāng)于極值點),鳥群中所有個體的最佳位置就可以看作整個鳥群的最佳覓食點(全局最優(yōu),相當(dāng)于最值點)。整個鳥類群體的覓食活動必然會移動到整個搜索范圍的最佳覓食區(qū),通過不斷改變鳥群覓食位置,即不斷迭代更新位置,速度也不斷更新,鳥群向最優(yōu)位置靠攏。

原理如下:在一個D維解空間中,有N個粒子構(gòu)成一個群體,每個粒子均以一定的速度移動。其中第i個粒子表示為一個D維向量Xi=[Xi1,Xi2,…,Xid]T(i=1,2,…,N,1≤d≤D),即第i個粒子在D維空間中的位置為Xi,將Xi代入目標(biāo)函數(shù)就可以計算出粒子i的適應(yīng)度,依靠適應(yīng)度值的大小來衡量粒子Xi的好壞;第i個粒子的速度表示為Vi=[Vi1,Vi2,…,Vid]T;第i個粒子自身的最優(yōu)值為Pbesti=[Pbesti1,Pbesti2,…,Pbestid]T;種群粒子最優(yōu)值為Gbesti=[Gbesti1,Gbesti2,…,Gbestid]T。粒子速度和位置每次更新公式為

(1)

(2)

式中:w為慣性權(quán)重,是用來均衡全局以及局部搜索的能力;c1和c2為學(xué)習(xí)因子,代表粒子向自己和群體中有優(yōu)秀粒子學(xué)習(xí)的能力;然后朝著歷史最優(yōu)點Pbest靠攏,一般取c1=c2=2;r1和r2為[0.1]的偽隨機數(shù)。粒子通過多次迭代更新來貼近最優(yōu)解,最終得到的Gbest就是算法找到的最優(yōu)解。

由式(1)可知粒子下一步運動速度由三部分決定:慣性權(quán)重和上一步的速度項,自我學(xué)習(xí)因子和個體最佳位置項,群體學(xué)習(xí)因子和群體最佳位置項。由式(2)可知下一步粒子,移動的位置由上一步的位置和上一步的速度來決定,具體流程如圖1所示。

在CA-PSO算法中,用局部進化替代流程圖中紅框圈出來的步驟

1.2 改進粒子群算法

本文提出的CA-PSO優(yōu)化算法的基本原理是:在原標(biāo)準(zhǔn)PSO算法的基礎(chǔ)上,增加局部進化策略來提高粒子迭代的效率,引入線性遞變的慣性權(quán)重來確保粒子速度過大跳出搜索區(qū)域,過小陷入局部極值。同時也引入收縮因子[17]來均衡算法的收斂性,同時在算法迭代時,采用輪換法則進行計算。

算法局部進化策略:假設(shè)初始種群的個數(shù)為N,小組個數(shù)為Nm,小組內(nèi)的粒子數(shù)為Nn=N/Nm,將種群粒子代入適應(yīng)度函數(shù)獲得適應(yīng)度,這里為了編程方便,按適應(yīng)度大小進行降序排序,按照轉(zhuǎn)輪法則分為Nmm。例如:分為10組,每組粒子數(shù)就有N/10個,則Nn=N/10,所有粒子就被表示為一個Nn行、10列的矩陣,為了便于計算,第一行數(shù)編號為(1, 1),(1, 2),……,(1,Nn);第二行編號(2, 1),(2, 2),……,(2,Nn), 其他行依次類推,這里采用輪換法則是為了便于編程計算。

收縮因子計算公式為

(3)

式(3)中:Nn為小組內(nèi)粒子個數(shù);i為當(dāng)前粒子數(shù)。

慣性權(quán)重的計算公式為

(4)

式(4)中:wmax為最大慣性權(quán)重,一般取0.9;wmin為最小慣性權(quán)重,一般取0.4;T為最大迭代次數(shù);t為當(dāng)前迭代次數(shù)。

每組中有一個最優(yōu)粒子,組內(nèi)迭代Nger次,適應(yīng)度越好收縮因子S越小,適應(yīng)度越差,收縮因子S越大;當(dāng)某個粒子連續(xù)Nger次都是適應(yīng)度比較差的,則將這個粒子初始化;組內(nèi)迭代完后再次按適應(yīng)度分組進行組內(nèi)迭代,直到迭代次數(shù)足夠或者最佳適應(yīng)度滿足所需條件時迭代停止。

小組內(nèi)更新公式為

(5)

(6)

1.3 算法流程

基于以上理論得到本文算法(CA-PSO),算法的實現(xiàn)步驟分為全局進化步驟和局部進化步驟兩部分。

全局進化步驟如下:①初始化種群(粒子群共有N個粒子)中的粒子數(shù)量,并給每個粒子一個隨機的初始位置和初始速度;②根據(jù)個體粒子的適應(yīng)度,將它們按降序排列,得到種群的全局最優(yōu)粒子Pbest;③按照轉(zhuǎn)輪法則將粒子劃分為若干個小組;④在每個小組內(nèi)部實行CA-PSO的局部進化策略;⑤判斷組內(nèi)局部進化步驟是否結(jié)束,是轉(zhuǎn)入步驟⑥,否轉(zhuǎn)入步驟④;⑥判斷種群迭代有沒有結(jié)束,是轉(zhuǎn)入步驟⑦,若沒有,則轉(zhuǎn)入步驟②。⑦最佳適應(yīng)度值所對應(yīng)的全局最優(yōu)粒子,就是最終結(jié)果,算法結(jié)束。

局部進化步驟如下:①初始化小組計數(shù)器ik和小組內(nèi)粒子迭代次數(shù)Nger;②根據(jù)小組計數(shù)器的值選擇局部進化的小組;③計算出小組內(nèi)的最好個體xmi;④依據(jù)式(5)和式(6)更新粒子位置和速度;⑤判斷粒子位置是否越過邊界,若是,則按無形/吸收的方式重新生成粒子;⑥計算組內(nèi)粒子的適應(yīng)度,并排序,按適應(yīng)度更新收縮因子;⑦確定組內(nèi)的迭代次數(shù)是否大于最大迭代次數(shù),是繼續(xù)步驟⑧,否繼續(xù)步驟③;⑧令迭代計數(shù)器ik=ik+1,判斷ik計數(shù)器是否大于種群數(shù),否轉(zhuǎn)入步驟②,是判斷組內(nèi)每個粒子的歷史最優(yōu)位置是否發(fā)生變化,若沒有,且適應(yīng)度排序靠后,否則重找此粒子,跳出局部進化。

2 數(shù)值實驗

2.1 測試函數(shù)

基于上述分析,利用MATLAB實現(xiàn)了CA-PSO算法,選取6個經(jīng)典測試函數(shù)[18]來檢驗CA-PSO算法的綜合性能,函數(shù)詳細信息如表1所示,最佳值是測試函數(shù)在搜索空間中找到的最佳結(jié)果。

表1 測試函數(shù)

2.2 對比算法

為進一步驗證CA-PSO算法的求解精度和收斂速度,將CA-PSO與其他PSO改進算法進行對比,分別為:粒子群算法PSO,包含變異的慣性和經(jīng)驗相互影響的算法(BEPSO)[19],重構(gòu)了模型收縮因子的算法(CIPSO)[20]、自適應(yīng)粒子群算法(IPSO)[21]。實驗在相同的硬件環(huán)境下進行,種群大小N為500;維數(shù)D為500;最大的迭代次數(shù)為1 000;算法運行次數(shù)是20,上述對比算法的詳細參數(shù)取值如表2所示,每種算法運行20次后取其平均最優(yōu)值,結(jié)果如表3所示。對比表3中的數(shù)據(jù)可以看出,與其他4種算法相比,CA-PSO在函數(shù)f1～f6上的均值都為0,這說明算法的尋優(yōu)精度比其他4種算法都高。同時BEPSO算法在f1偶次多項式函數(shù)和f4具有較高尋優(yōu)難度上的均值不為0,說明該算法收斂度和收斂精度較低,陷入了局部誤差。

表2 對比算法的參數(shù)列表

表3 算法結(jié)果對比

2.3 收斂曲線分析

圖2和圖3顯示了4種算法在6個測試函數(shù)上的適應(yīng)度平均值的變化曲線。可以很直觀地看出,隨著迭代次數(shù)的增多,CA-PSO算法在6種測試函數(shù)上都能快速收斂到最優(yōu)解。f2、f3、f6中CIPSO是最慢找到最優(yōu)解的,f1、f4中BEPSO是最慢收斂到最優(yōu)解的。結(jié)果表明,CA-PSO算法比其他4種對比算法有更好的全局搜索能力和更快的收斂速度。

圖2 單峰測試函數(shù)收斂圖

圖3 多峰測試函數(shù)收斂圖

3 MT反演

3.1 MT一維反演

一維定義理論模型和實際模型的相對均方差作為目標(biāo)函數(shù)φ,即

(7)

構(gòu)建了2個三層模型:K和H型,模型參數(shù)如表4所示。取初始種群數(shù)50,空間維度為5(3個電阻率,2個層厚),最大迭代次數(shù)1 000,最大誤差界限常數(shù)ε=1×10-4,反演結(jié)果如表5所示。

表4 模型參數(shù)

表5 反演結(jié)果

一維大地電磁改進粒子群算法反演具體步驟如下。

(1)初始化:根據(jù)模型層數(shù)和反演參數(shù)的多少及問題復(fù)雜程度確定相關(guān)參數(shù),粒子維度D與地層參數(shù)(2n-1)個數(shù)一致。

(3)擇優(yōu):根據(jù)評價的情況,選擇當(dāng)前最優(yōu)適應(yīng)度值對應(yīng)的粒子,并判斷是否滿足中止條件,如滿足,則停止搜素;否則,進入下一步種群更新。

(4)更新:對粒子群按照式(5)和式(6)更新當(dāng)前種群(小組內(nèi))個體對應(yīng)的參數(shù)值,再更新個體粒子適應(yīng)度,從而完成當(dāng)前整個種群的更新。

(5)進入下一代循環(huán):算法轉(zhuǎn)至步驟(2)繼續(xù)執(zhí)行,直到:①不滿足迭代次數(shù)t>2且φ<ε(最大誤差界限常數(shù))時;②最大迭代次數(shù)為100;③φ<1×10-5中的任何一個條件時迭代中止。

CA-PSO算法反演與PSO算法反演流程區(qū)別在于步驟(4)中多了一步添加慣性權(quán)重和收縮因子的小組迭代操作。

由表5可以看出,H型曲線的反演電阻率均比理論模型電阻率小,最大60 Ω·m,最小0.1 Ω·m,層厚但也在40 m以內(nèi)有略微差異;K型曲線反演電阻率均比理論模型大,范圍在0.02～42 Ω·m,層厚更接近理論模型,差距在10 m以內(nèi),K型反演效果更好,總體來說該方法可用。

一般來說,H型測深曲線模型在反演中最難解決,所以這里只給出H型曲線模型的反演結(jié)果圖,如圖4所示。從圖4可以看出,LS與模型的一致性最高,而粒子群算法的一致性最差。CA-PSO方法與理論模型在650 m左右的深度范圍內(nèi)曲線基本一致,表明CA-PSO方法適用于一維MT是可行的。

圖4 H型模型結(jié)果圖

3.2 MT二維反演

CA-PSO用在MT二維反演上,由于維度增加,計算量大且計算時間長。為了改善計算效果,二維反演的目標(biāo)函數(shù)跟一維相比加入核函數(shù)[22-23]。例如:一組輸入樣本{xi,yi}i=1,2,…,m,輸入量為xi∈Rn,輸出量為yi∈R,通過核函數(shù)k(·,·)將在低維映射到更高的維度。目標(biāo)函數(shù)為

φ=k(xi,yi)w+b

(8)

式(8)中:w為慣性權(quán)重;b為偏置量,增加平移的能力,一般取值為2;{xi,yi}i=1,2,…,m類比為

式(8)的求解轉(zhuǎn)化為規(guī)劃問題的最優(yōu)解。

核函數(shù)取為Gaussian核函數(shù),即

(9)

式(9)中:σ為輸入?yún)?shù),取值為1×10-3。

(10)

為驗證算法在二維反演的適用性,設(shè)計了從簡單到復(fù)雜的地電模型進行測試,地表布設(shè)一條長6 km的測線,13個測點,頻率取0.01、1、100 Hz,計算這些測點處的TE(橫電波)和TM(橫磁波)極化模式下的視電阻率,來進行擬合反演尋找與所設(shè)異常體模型相對應(yīng)的在背景場中反演異常體的位置。設(shè)置初始種群數(shù)N=100, 搜索空間維度D=10,最大迭代次Nger=1 000,學(xué)習(xí)因子c1=c3=2,c2=0.8,搜索空間的其他參數(shù)跟一維類似。背景電阻率均為100 Ω·m。

(1)模型一:該模型(圖5)反映的是1個低阻異常體,埋深1 km,規(guī)模為2 km×1 km,異常體電阻率10 Ω·m。由圖6可以清楚地看出,異常體位置與實際模型相符合,表明該方法反演單個規(guī)則異常體效果很好。

圖5 單異常低阻模型

圖6 單異常低阻反演結(jié)果

(2)模型二:該模型(圖7)反映的是2個一樣大小的低阻異常體,規(guī)模為2 km×2 km,埋深 1 km,異常體電阻率均為30 Ω·m。由圖8可以明顯地看出,異常體與實際模型的位置非常相似,表明該方法對于反演2個同樣大小的低阻異常體是可行的。

圖7 雙異常低阻模型

圖8 雙異常低阻反演結(jié)果

(3)模型三:該模型(圖9)反映的是一個電阻率為30 Ω·m,規(guī)模3 km×0.5 km,埋深 0.5 km的低阻異體,和一個電阻率為200 Ω·m,規(guī)模3 km×3 km,埋深 2 km的高阻異常體。 由圖10可以看出,圖10(b)比圖10(a)能更好地反映實際模型,這表明在目標(biāo)函數(shù)中只引入核函數(shù)和又加入模型修改量相比,后者反演效果更好,而且可以很清楚地區(qū)分高低阻。對于圖10(a)反演結(jié)果較差,是因為少了模型修改量和一些限制。綜上,改進算法CA-PSO應(yīng)用在MT二維反演上可以很好地定位異常體位置,而且對于多個異常體也可以明顯地反演出來,該改進方法用在二維反演是可行的。

圖9 高低阻雙異常模型

圖10 高低阻雙異常反演結(jié)果

4 結(jié)論

基于以上內(nèi)容,得到以下結(jié)論。

(1)提出了一種基于增加局部進化策略和加入收縮因子的CA-PSO優(yōu)化算法。算法加入了含有遞變的慣性權(quán)重和收縮因子的局部進化步驟,將原來需要每個粒子之間相互對比的方式,變成了只需每個組的最優(yōu)粒子進行對比,不僅兼顧了粒子集的多樣性和收斂性而且提高了粒子進化的效率,也較好地避免了算法過早收斂和陷入局部誤差的問題;在測試函數(shù)上的對比結(jié)果表明,該算法不依賴于初始模型,計算簡單,在一定程度上改善了PSO算法容易陷入局部極值的優(yōu)點。

(2)通過設(shè)置目標(biāo)函數(shù),在目標(biāo)函數(shù)中加入先驗信息限制,利用CA-PSO算法進行一、二維MT反演,能夠準(zhǔn)確地進行異常體的定位,反演結(jié)果與實際模型吻合,表明了算法應(yīng)用于MT反演的可行性。

(3)利用CA-PSO算法實現(xiàn)二維MT反演,在目標(biāo)函數(shù)中引入核函數(shù),將傳統(tǒng)擬合函數(shù)所使用的二范數(shù)改為一范數(shù),能夠減少計算量,得到準(zhǔn)確的反演結(jié)果。

(4)由于粒子群算法是全局搜索,維度越高,計算量越大,計算時間也越長,所以目前很少有人用改進的PSO算法來做MT的二維反演,大都做的是一維反演。所以本文研究做二維反演只是一種嘗試,而且每計算一次時間較長,計算結(jié)果也是理想化的,沒有考慮實際的因素。對于減少二維反演計算時間和實際問題的研究,是下一步的研究計劃。