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

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

0 引言

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

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

1 算法理論

1.1 GB分解原理

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

Z=RCZ2DRT

C=gTSA

(1)

其中：

(2)

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

(3)

1.2 標(biāo)準(zhǔn)粒子群優(yōu)化算法

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

粒子群算法通常的數(shù)學(xué)描述為設(shè)在一個(gè)N維空間中，由m個(gè)粒子組成的種群X=(x1,…,xi,…,xN)，其中第i個(gè)粒子的位置為xi=(xi1,xi2,…,xiN)T, 其中速度為Vi=(vi1,vi2,…,vid,…,viN)T。它的個(gè)體極值為pi=(pi1,pi2,…,piN)T,種群的全局極值為pg=(pg1,pg2,…,pgN)T,按照追隨當(dāng)前最優(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)

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

圖1 粒子群算法流程圖

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

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

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

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)

慣性權(quán)重對(duì)粒子群算法來(lái)說(shuō)是一個(gè)非常重要的參數(shù)，它的物理意義是對(duì)當(dāng)前粒子速度的繼承，標(biāo)準(zhǔn)粒子群算法的慣性權(quán)重，通常是一個(gè)常數(shù)或者隨著迭代次數(shù)的增加成線性遞減，但在實(shí)際應(yīng)用中，有其不適用之處，容易導(dǎo)致收斂過(guò)慢或者陷入局部極值中。若將慣性權(quán)重設(shè)定為服從某種隨機(jī)分布的隨機(jī)數(shù)，這樣就可以在一定程度上克服傳統(tǒng)慣性權(quán)重的不足。如設(shè)ω'=ω-c1r1-c2r2，因?yàn)閞1、r2是(0,1)之間的隨機(jī)數(shù)，所以ω'實(shí)際上也是一個(gè)隨機(jī)數(shù)。因此，式(6)和式(7)就可以改寫(xiě)為式(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)

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

2 算法驗(yàn)證

2.1 仿真實(shí)驗(yàn)

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

Matyas函數(shù)：

Sphere函數(shù)：

測(cè)試函數(shù)及算法參數(shù)如表1所示，Matyas函數(shù)維數(shù)為2，Sphere函數(shù)維度為50，除慣性權(quán)重，其余參數(shù)均一致。通過(guò)高低兩個(gè)維度的測(cè)試函數(shù)來(lái)測(cè)試在各算法的性能表現(xiàn)，實(shí)驗(yàn)結(jié)果如表2所示，算法1在低維度函數(shù)中勉強(qiáng)達(dá)到求解要求，但在高維度函數(shù)中計(jì)算結(jié)果出現(xiàn)較大偏差；算法2的計(jì)算表現(xiàn)稍強(qiáng)于算法1，但不及算法3和算法4，平均耗時(shí)卻少于算法3；算法3是計(jì)算結(jié)果最好的，但平均耗時(shí)不如算法2和算法4；算法4計(jì)算表現(xiàn)雖略遜于算法3，但平均耗時(shí)卻強(qiáng)于算法3?？偟膩?lái)說(shuō)，除算法1以外，其余算法均能達(dá)到測(cè)試計(jì)算要求。

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

表2 實(shí)驗(yàn)結(jié)果

2.2 理論模型

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

2.2.1 單頻點(diǎn)數(shù)據(jù)

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

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

圖2 單頻點(diǎn)理論模型GB分解結(jié)果

2.2.2 三維正演模型

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

圖3 三維/二維構(gòu)造模型示意圖

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

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

圖5 三維/二維模型GB分解結(jié)果

2.3 實(shí)際資料分解

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

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

圖7 實(shí)際資料分解結(jié)果

3 結(jié)論

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

1)單純使用標(biāo)準(zhǔn)粒子算法無(wú)法滿足計(jì)算要求，而該進(jìn)后的隨機(jī)慣性權(quán)重二階振蕩粒子群算法，無(wú)論在計(jì)算精度還是效率上都能滿足需求。

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

3)改進(jìn)的粒子群算法，具有較強(qiáng)的隨機(jī)性，控制參數(shù)較多，為達(dá)到滿意的計(jì)算結(jié)果，可能需要反復(fù)試驗(yàn)，多次計(jì)算，但過(guò)程簡(jiǎn)單方便，易于實(shí)現(xiàn)。