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

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

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

0 引言

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

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

1 算法理論

1.1 GB分解原理

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

Z=RCZ2DRT

C=gTSA

(1)

其中：

(2)

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

(3)

1.2 標準粒子群優化算法

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

粒子群算法通常的數學描述為設在一個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,按照追隨當前最優粒子的原理，粒子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)

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

圖1 粒子群算法流程圖

1.3 改進的粒子群優化算法

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

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

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)

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

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

Matyas函數：

Sphere函數：

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

表1 測試函數及算法參數

表2 實驗結果

2.2 理論模型

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

2.2.1 單頻點數據

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

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

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

2.2.2 三維正演模型

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

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

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

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

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

2.3 實際資料分解

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

圖6 工區地質略圖[20]

圖7 實際資料分解結果

3 結論

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

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

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

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