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

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

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

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

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

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

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

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

1.1 粒子群算法

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

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

(1)

(2)

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

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

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

1.2 改進粒子群算法

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

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

收縮因子計算公式為

(3)

式(3)中:Nn為小組內粒子個數;i為當前粒子數。

慣性權重的計算公式為

(4)

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

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

小組內更新公式為

(5)

(6)

1.3 算法流程

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

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

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

2 數值實驗

2.1 測試函數

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

表1 測試函數

2.2 對比算法

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

表2 對比算法的參數列表

表3 算法結果對比

2.3 收斂曲線分析

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

圖2 單峰測試函數收斂圖

圖3 多峰測試函數收斂圖

3 MT反演

3.1 MT一維反演

一維定義理論模型和實際模型的相對均方差作為目標函數φ,即

(7)

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

表4 模型參數

表5 反演結果

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

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

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

(4)更新:對粒子群按照式(5)和式(6)更新當前種群(小組內)個體對應的參數值,再更新個體粒子適應度,從而完成當前整個種群的更新。

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

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

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

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

圖4 H型模型結果圖

3.2 MT二維反演

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

φ=k(xi,yi)w+b

(8)

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

式(8)的求解轉化為規劃問題的最優解。

核函數取為Gaussian核函數,即

(9)

式(9)中:σ為輸入參數,取值為1×10-3。

(10)

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

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

圖5 單異常低阻模型

圖6 單異常低阻反演結果

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

圖7 雙異常低阻模型

圖8 雙異常低阻反演結果

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

圖9 高低阻雙異常模型

圖10 高低阻雙異常反演結果

4 結論

基于以上內容,得到以下結論。

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

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

(3)利用CA-PSO算法實現二維MT反演,在目標函數中引入核函數,將傳統擬合函數所使用的二范數改為一范數,能夠減少計算量,得到準確的反演結果。

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