基于改進粒子群算法的大地電磁阻抗張量分解方法

陳先潔, 王緒本, 李德偉, 謝卓良, 乃國茹

(成都理工大學 地球物理學院，成都 610000)

0 引言

在大地電磁測深法中往往存在由于淺地表電性結構不均勻，使得在不均勻體的周圍電流密度發(fā)生畸變，從而使電場分量發(fā)生變化，引起局部的畸變異常[1-4]。為消除此種局部畸變異常帶來的影響，Groom和Bailey假設局部三維異常體覆蓋在二維區(qū)域異常之上，將實際阻抗張量分解為畸變阻抗張量和區(qū)域阻抗張量，以此來恢復未受畸變影響的區(qū)域阻抗張量(GB分解法)。它的本質(zhì)是求解由區(qū)域阻抗張量和畸變方程組成的非線性超定方程組[6-8]。傳統(tǒng)的GB分解法將非線性問題線性化求解，導致求解結果不穩(wěn)定，易陷入局部極值中，很難得到較為準確的分解結果。李洋等[5]提出來一種混合優(yōu)化算法，利用模擬退火法得到的全局最優(yōu)解作為局部非線性最小二乘法的初始值，提高了求解穩(wěn)定性；尹曜田等[6]研究了基于遺傳算法的阻抗張量分解方法，采用遺傳算法求解非線性方程組；蔡軍濤等[7]提出共軛阻抗變換的概念，不受電場畸變影響，無需地下電性維性信息的假設；謝成良等[8]研究了基于相位張量約束下的阻抗張量分解方法，利用相位張量不受電場畸變影響的優(yōu)點，分析得到合理的初始模型，提高了計算穩(wěn)定性和計算效率，但仍然存在對初始值依賴性嚴重等問題。

筆者在求解GB分解法中的非線性超定方程組時,采用了粒子群優(yōu)化算法，但傳統(tǒng)的粒子群算法極易陷入局部極值當中，無法滿足求解要求。為此，筆者在前人研究基礎上，在傳統(tǒng)粒子群算法中引入振蕩環(huán)節(jié)和隨機慣性權重等思想，形成基于隨機慣性權重的二階振蕩粒子群算法，并利用該算法實現(xiàn)大地電磁阻抗張量分解，通過理論數(shù)據(jù)和實測資料的分解處理，驗證了改進算法在阻抗張量分解中的適用性和有效性。

1 算法理論

1.1 GB分解原理

GB分解法針對三維/二維電性結構模型，將觀測阻抗分解為區(qū)域阻抗張量和畸變阻抗張量，以此恢復未受畸變影響的區(qū)域阻抗張量，進而求得區(qū)域構造走向的角度以及電性主軸的特征。在不考慮磁場的感應型畸變的情況下，觀測阻抗張量Z可以分解為：

Z=RCZ2DRT

C=gTSA

(1)

其中：

(2)

上述各矩陣帶入式(2)中得到：

(3)

1.2 標準粒子群優(yōu)化算法

粒子群最早是由Eberhaet等提出[9]，該算法概念源于鳥群 的覓食行為，通過 “鳥群”中 “每只鳥”之間的協(xié)作和信息共享來尋找最優(yōu)解[10]。粒子群算法采用具有速度和位置屬性的粒子來模擬鳥群中的鳥，速度代表移動的快慢，位置代表移動的方向。在N維空間中隨機生成若干粒子，每個粒子在搜索空間中單獨的搜尋最優(yōu)解，并將其記為當前個體極值，并將個體極值與整個粒子群里的其他粒子“共享”，并找到其中最優(yōu)的那個個體極值作為整個粒子群的當前全局最優(yōu)解，粒子群中的所有粒子，根據(jù)自己找到的當前個體極值和整個粒子群共享的當前全局最優(yōu)解，來調(diào)整自己的速度和位置[10-14]。粒子群算法雖然是一個非常簡單的算法，但其應用領域十分廣泛，主要包括函數(shù)優(yōu)化，約束優(yōu)化，工程設計，神經(jīng)網(wǎng)絡訓練等，具有簡單高效且不依賴初始值等特點。

粒子群算法通常的數(shù)學描述為設在一個N維空間中，由m個粒子組成的種群X=(x1,…,xi,…,xN)，其中第i個粒子的位置為xi=(xi1,xi2,…,xiN)T, 其中速度為Vi=(vi1,vi2,…,vid,…,viN)T。它的個體極值為pi=(pi1,pi2,…,piN)T,種群的全局極值為pg=(pg1,pg2,…,pgN)T,按照追隨當前最優(yōu)粒子的原理，粒子xi將按照式(4)，式(5)改變自己的速度和位置。

vij(t+1)=ωvij(t)+c1r1(t)(pij(t)-

xij(t))+c2r2(t)(pgj(t)-xij(t))

(4)

xij(t+1)=xij(t)+vij(t+1)

(5)

其中：ω是慣性權重系數(shù)，其值較大時，全局尋優(yōu)能力強，局部尋優(yōu)能力弱；其值較小時，全局尋優(yōu)能力弱，局部尋優(yōu)能力強;j=1、2、…、N，i=1、2、…、m；m為種群規(guī)模；t為當前進化代數(shù)；r1、r2為分布于[0,1]之間的隨機數(shù)；c1、c2為加速常數(shù)。從式(4)中可知，每個粒子由三部分組成：第一部分為“記憶”部分，表示上次速度大小和方向的影響；第二部分為“認知”部分，是從當前點指向粒子自身最好點的一個矢量，表示粒子自身的“思考”；第三部分為“社會”部分，是一個從當前點指向種群最好點的矢量，表示粒子之間的” 信息共享與協(xié)同合作”[14]。在使用粒子群算法時，需要重點設置待定參數(shù)變化范圍和參數(shù)的每步迭代最大允許值，尤其是參數(shù)的每步迭代最大允許值，一般為變化范圍的10%～20%，該值越小，收斂分辨率越高，但收斂速度慢，反之，該值越大，收斂速度越快，但易跳出全局最優(yōu)值，故需要反復試驗，以達到最佳優(yōu)化效果。圖1為粒子群算法流程圖。

圖1 粒子群算法流程圖

1.3 改進的粒子群優(yōu)化算法

標準的粒子群算法的性能很大程度上依賴于其參數(shù)(慣性權重、隨機因子、加速常數(shù)、最大迭代步長等)，這就導致在應用過程中出現(xiàn)容易喪失種群多樣性，過早陷入局部極值，精度不高等問題。為此國內(nèi)、外學者提出了許多種改進方法，其中包括將其他優(yōu)化算法思想引入到粒子群算法中(如基于模擬退火的粒子群優(yōu)化算法[15])，另外通過改進標準粒子群中的速度和位置公式(式(4)和式(5))，提高算法性能也是行之有效的方法(如在速度公式中增加二次項，粒子位置采用高斯分布[16]等)。不管從什么角度對粒子群算法進行改進，都能一定程度上提高粒子群算法的性能。筆者在二階振蕩粒子群算法的基礎上，將隨機慣性權重引入該算法中，以此改善算法的全局收斂性，提高算法精度。

二階振蕩粒子算法由國內(nèi)學者胡建秀等[17]提出，在速度方程中增加了二次項和振蕩環(huán)節(jié)，使算法在前期具有較強全局收斂能力，呈振蕩收斂，而在后期加強了局部搜索能力，呈漸進收斂。蔣麗等[18]通過改進在振蕩環(huán)節(jié)中的參數(shù)及其取值范圍，進而加強全局搜索能力，改進后的速度和位置公式如下：

vij(t+1)=ωvij(t)+c1r1(t)(pij(t)-

(1+ξ1)xij(t)+ξ2xij(t-1))+

c2r2(t)(pgj(t)-(1+ξ3)xij(t)+

ξ4xij(t-1))

(6)

xij(t+1)=xij(t)+vij(t+1)

(7)

慣性權重對粒子群算法來說是一個非常重要的參數(shù)，它的物理意義是對當前粒子速度的繼承，標準粒子群算法的慣性權重，通常是一個常數(shù)或者隨著迭代次數(shù)的增加成線性遞減，但在實際應用中，有其不適用之處，容易導致收斂過慢或者陷入局部極值中。若將慣性權重設定為服從某種隨機分布的隨機數(shù)，這樣就可以在一定程度上克服傳統(tǒng)慣性權重的不足。如設ω'=ω-c1r1-c2r2，因為r1、r2是(0,1)之間的隨機數(shù)，所以ω'實際上也是一個隨機數(shù)。因此，式(6)和式(7)就可以改寫為式(8)和式(9)。

vij(t+1)=(ω-c1r1(t)-c2r2(t))vij(t)+

c1r1(t)(pij(t)-(1+ξ1)xij(t)+

ξ2xij(t-1))+c2r2(t)(pgj(t)-

(1+ξ3)xij(t)+ξ4xij(t-1))

(8)

xij(t+1)=xij(t)+vij(t+1)

(9)

基于此，筆者將隨機慣性權重引入到改進后的二階振蕩粒子群中，使算法的全局搜索能力和收斂速度得到一定程度地提高。

2 算法驗證

2.1 仿真實驗

利用兩個測試函數(shù)來對上述改進方法進行試驗比較，測試函數(shù)分別是Matyas函數(shù)和Sphere函數(shù)。為體現(xiàn)仿真實驗的有效性，利用標準粒子群算法(算法1)，恒定慣性權重二階振蕩粒子算法(算法2)，線性慣性權重二階振蕩粒子群算法(算法3)，隨機慣性權重二階振蕩粒子群算法(算法4)，對測試函數(shù)分別進行100次仿真計算實驗，并統(tǒng)計分析計算后的數(shù)據(jù)，從最優(yōu)值，平均值，標準差，平均耗時等對算法進行統(tǒng)計評估。

Matyas函數(shù)：

Sphere函數(shù)：

測試函數(shù)及算法參數(shù)如表1所示，Matyas函數(shù)維數(shù)為2，Sphere函數(shù)維度為50，除慣性權重，其余參數(shù)均一致。通過高低兩個維度的測試函數(shù)來測試在各算法的性能表現(xiàn)，實驗結果如表2所示，算法1在低維度函數(shù)中勉強達到求解要求，但在高維度函數(shù)中計算結果出現(xiàn)較大偏差；算法2的計算表現(xiàn)稍強于算法1，但不及算法3和算法4，平均耗時卻少于算法3；算法3是計算結果最好的，但平均耗時不如算法2和算法4；算法4計算表現(xiàn)雖略遜于算法3，但平均耗時卻強于算法3。總的來說，除算法1以外，其余算法均能達到測試計算要求。

表1 測試函數(shù)及算法參數(shù)

表2 實驗結果

2.2 理論模型

從計算精度和計算效率兩方面考慮，算法4相對更加適用于大地電磁阻抗張量分解。為驗證其有效性，利用單個頻點數(shù)據(jù)和三維模型數(shù)據(jù)進行分解試驗。

2.2.1 單頻點數(shù)據(jù)

利用算法4，經(jīng)過反復試驗，其中主要幾個控制參數(shù)設置如下：種群中粒子數(shù)為24個，加速常數(shù)c1、c2均為2，最大迭代次數(shù)為2 000次，閾值為1e-25，最大迭代步長為取值范圍的15%。 對上述畸變數(shù)據(jù)計算100次，計算出的統(tǒng)計結果為：扭曲因子為2.101°，剪切因子為24.68°，構造走向角為0.070°。所得結果如圖2所示，其中橫軸為計算次數(shù)，縱軸為計算結果。

從圖2可看出， 采用改進算法計算出100次結果非常穩(wěn)定，并且每次的結果都十分接近真實值，誤差在允許范圍之內(nèi)，說明改進算法是能夠準確地計算出模型的各個參數(shù)，進而可以達到理想的分解效果。

圖2 單頻點理論模型GB分解結果

2.2.2 三維正演模型

為進一步驗證本文算法的分解效果，依據(jù)GB分解法要求的假設條件，建立一個三維/二維電性結構模型。其中背景區(qū)域電阻率為500 Ω·m，二維電性結構為一個無限延伸的低阻體，電阻率為100 Ω·m，局部三維塊體電阻率為1 Ω·m。計算所用網(wǎng)格為44×45×60，計算頻段為0.01 Hz～100 Hz，共20個頻點。這里選取典型測點進行分析，測點與異常體相對位置如圖3所示。圖4是畸變前、后的視電阻率曲線，畸變前曲線表現(xiàn)出明顯的二維特征，受到局部畸變影響后，視電阻率曲線出現(xiàn)平移，表現(xiàn)出強烈的三維性。

圖3 三維/二維構造模型示意圖

圖4 畸變前后視電阻率曲線

利用改進算法進行分解處理，得到如圖5所示的分解結果。從結果上看，在正演頻段內(nèi)，剪切角分布在0°到20°之間，扭曲角分布在0°到10°之間，電性主軸角分布在5°到33°之間。這表明利用基于粒子群算法的GB分解法所計算出的模型各參數(shù)均十分穩(wěn)定且可靠。

圖5 三維/二維模型GB分解結果

2.3 實際資料分解

BC87數(shù)據(jù)集(數(shù)據(jù)來源：https://www.mtnet.info/data/bc87/bc87.html)位于大不列顛哥倫比亞省東南部，大地電磁剖面長度為150 km，總計27個測點,觀測周期為0.002 s～1 820 s.數(shù)據(jù)質(zhì)量較好，但地質(zhì)構造復雜，受到的畸變效應嚴重，為我們研究阻抗張量分解方法，恢復區(qū)域阻抗響應研究提供了有效的實測數(shù)據(jù)集。該地區(qū)的地質(zhì)略圖如圖6所示，眾多學者對該數(shù)據(jù)集進行了分析研究[17-18]。筆者選擇其中位于Valhalla(瓦爾哈拉片麻巖)附近的lit000號測點(圖6中紅色測點)進行阻抗張量分解研究，并通過商業(yè)軟件(MT-Pioneer)處理，得到分解結果真值，對比結果如圖7所示，其中走向角(Strike)集中分布在0°到60°，扭曲因子(Twist)集中分布在-60°到30°，剪切因子(Shear)集中分布在0°到30°。利用本文算法得到的計算結果與商業(yè)軟件的分解結果基本一致，但由于實測數(shù)據(jù)存在噪聲干擾，對分解結果造成一定干擾，也就導致某些頻點的分解結果與真實值存在差異。

圖6 工區(qū)地質(zhì)略圖[20]

圖7 實際資料分解結果

3 結論

在前人的研究基礎上，將粒子群算法引入到大地電磁阻抗張量分解方法中，并對算法進行改進，提高其算法性能。利用測試函數(shù)對算法進行仿真實驗，并經(jīng)過對理論數(shù)據(jù)和實測數(shù)據(jù)的阻抗張量分解試驗，主要得到以下結論：

1)單純使用標準粒子算法無法滿足計算要求，而該進后的隨機慣性權重二階振蕩粒子群算法，無論在計算精度還是效率上都能滿足需求。

2)GB分解算法對初始模型依賴性較強，初始模型選擇對算法結果影響較大，如選擇不當，會造成計算結果偏差較大，而粒子群算法隨機生成初始值，在一定程度上降低模型初始值對分解效果的影響。再者，若初始值并不理想，改進算法中的隨機慣性權重也有可能產(chǎn)生相對較大的值，加強全局搜索能力，避免陷入局部極值。

3)改進的粒子群算法，具有較強的隨機性，控制參數(shù)較多，為達到滿意的計算結果，可能需要反復試驗，多次計算，但過程簡單方便，易于實現(xiàn)。