基于粒子群算法的GPR與電測深共同反演及等值現(xiàn)象消除

戴前偉，張 豪，張 彬

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2.中南大學(xué) 有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室，長沙 410083)

0 引言

在地球物理勘探中，以巖礦石間導(dǎo)電性的差異為基礎(chǔ)，通過觀測和研究人工電場的地下分布特點(diǎn)，來解決實(shí)際問題的電阻率法是電法勘探的一個(gè)重要分支，電阻率ρ是描述導(dǎo)電性強(qiáng)弱的電性參數(shù)。在實(shí)際的工作中，不可避免的會(huì)存在著一定的觀測誤差，從而使得不同的地電斷面所對(duì)應(yīng)的電測深曲線之間出現(xiàn)等值現(xiàn)象，造成錯(cuò)誤的解釋。而探地雷達(dá)則是根據(jù)高頻電磁波在地下介質(zhì)中的傳播特性，探測地下結(jié)構(gòu)和特性的一種地球物理方法。依照電磁場與電磁波理論，介質(zhì)的介電常數(shù)、電導(dǎo)率和磁導(dǎo)率的差異決定了電磁波的穿透深度、傳播路徑、波形以及發(fā)生反射、繞射現(xiàn)象等。因此，探地雷達(dá)的勘探結(jié)果不會(huì)出現(xiàn)等值現(xiàn)象。

絕大多數(shù)的地球物理反演問題為非線性問題[1-2]，完全非線性反演的各類方法仍在不斷地更新和發(fā)展。粒子群優(yōu)化算法(Particle Swarm Optimization, PSO)是由Eberhart和Kennedy[3-5]通過研究鳥群社會(huì)系統(tǒng)，提出的一種基于迭代的群智能隨機(jī)優(yōu)化算法。算法的基本概念源于對(duì)鳥群覓食過程中的遷徙和集聚行為的研究。近年來粒子群反演越來越多的應(yīng)用于地球物理領(lǐng)域。韓家興[6]提出了基于粒子群優(yōu)化算法的多元線性擬合法；孫中科[7]進(jìn)行了重力張量的粒子群反演研究；張進(jìn)[8]提出了加入模擬退火算子的隨機(jī)粒子群算法用于反演彈性阻抗；羅德相[9]提出了多種群粒子群優(yōu)化算法；易遠(yuǎn)元等[10]以地震波阻抗數(shù)值模擬和實(shí)際資料為例，說明了粒子群反演方法的可行性及有效性。

一條n層電測深曲線，可能對(duì)應(yīng)一組不同參數(shù)(ρ,h)的n層地電斷面，稱為同層等值現(xiàn)象。它包括S等值和T等值兩類[11]。為了解決電法勘探的S等值現(xiàn)象，筆者提出一種改進(jìn)的多種群粒子群反演算法，在原有的速度更新公式的基礎(chǔ)上，增加一個(gè)新的全局極值點(diǎn)，實(shí)現(xiàn)電法參數(shù)和探地雷達(dá)參數(shù)的共同反演，得到更加準(zhǔn)確的地電模型。

1 同層等值現(xiàn)象

根據(jù)勢場問題解的唯一性定理，一定的地電斷面對(duì)應(yīng)唯一的電測深曲線[12]。而實(shí)際工作中，當(dāng)有些不同地電斷面所對(duì)應(yīng)的電測深曲線間的差別在觀測誤差范圍以內(nèi)，常將其看成為“同一條”電測深曲線，這種情況稱為電測深曲線的等值現(xiàn)象。等值現(xiàn)象的存在，使得一條實(shí)測電測深曲線可對(duì)應(yīng)一組不同的地電斷面，從而造成錯(cuò)誤的解釋。

電法勘探的同層等值現(xiàn)象包括S等值和T等值兩類。本次實(shí)驗(yàn)的重點(diǎn)在于解決S等值現(xiàn)象。因?yàn)殡娦詫訁?shù)決定著電阻率轉(zhuǎn)換函數(shù)，相同轉(zhuǎn)換函數(shù)的地電斷面對(duì)應(yīng)著相同的電測深曲線，所以可以對(duì)轉(zhuǎn)換函數(shù)的等值性進(jìn)行分析，研究S等值現(xiàn)象的原因。

以三層曲線為例，對(duì)于H型和A型三層介質(zhì)，各層電阻率為ρ1、ρ2、ρ3,層厚度為h1、h2、h3，電阻率轉(zhuǎn)換函數(shù)T1(m)可表示為

T1(m)=ρ1cth{mh1+cth-1[μ2·

cth(mh1v2+cth-1μ3)]}

(1)

v2=h2/h1

(2)

μ2=ρ2/ρ1

(3)

μ3=ρ3/ρ2

(4)

其中：cth(x)為雙曲函數(shù)中的余切雙曲函數(shù)。當(dāng)v2?1，且μ3?1，即第二層的層厚度小且第三層電阻率值遠(yuǎn)大于第二層電阻率時(shí)

(5)

(6)

故

(7)

由此可見，第二層厚度或電阻率發(fā)生改變時(shí)，只要S2=h2/ρ2保持不變，則T1(m)不變，故稱為S等值現(xiàn)象。從前面分析可知，發(fā)生S等值現(xiàn)象的條件是：v2?1即第二層薄，且μ3?1。中間層越薄，ρ2越小等值范圍越寬。

2 粒子群算法

2.1 標(biāo)準(zhǔn)PSO算法

(8)

(9)

解空間規(guī)定了粒子飛行的范圍，在每次迭代中，粒子通過速度與位置公式不斷進(jìn)行著更新。粒子的速度更新公式分為三部分，①表示粒子當(dāng)前的飛行速度，提供粒子在搜索空間中飛行的動(dòng)力；②被稱為“個(gè)體認(rèn)知”部分，代表了粒子的個(gè)體經(jīng)驗(yàn)，促使粒子朝著自身所經(jīng)歷過的最優(yōu)位置移動(dòng)；③被稱為“群體認(rèn)知”部分，代表了群體經(jīng)驗(yàn)對(duì)粒子飛行軌跡的影響，促使粒子朝著目前群體發(fā)現(xiàn)的最優(yōu)位置移動(dòng)。

2.2 線性調(diào)整w的策略

線性調(diào)整w的算法，又稱為慣性權(quán)重線性遞減算法[14]。使用該算法對(duì)問題進(jìn)行優(yōu)化處理時(shí)，初始采用全局搜索，使搜索空間快速收斂于某一區(qū)域，之后采用局部精細(xì)搜索，可以求得更高精度的解。在vi前乘以慣性權(quán)重w，w較大時(shí)，算法具有較強(qiáng)的全局搜索能力；w較小，算法傾向于局部搜索。w調(diào)整公式為式(10)。

(10)

其中：wmax、wmin分別是w的最大值和最小值；iter和itermax分別是當(dāng)前迭代次數(shù)和最大迭代次數(shù)。

2.3 多種群粒子群優(yōu)化算法

多種群粒子群優(yōu)化算法與標(biāo)準(zhǔn)粒子群優(yōu)化算法的區(qū)別在于，速度更新公式中，除了當(dāng)前速度、“個(gè)體認(rèn)知”和“群體認(rèn)知”外，再增加一個(gè)不同約束條件下的“群體認(rèn)知”，以獲得問題的最優(yōu)解。速度更新公式為式(11)。

(11)

對(duì)于標(biāo)準(zhǔn)PSO易陷入局部極值的問題，多種群粒子群優(yōu)化算法可以通過新的約束條件得到相同解空間下的新的群體最佳，跳出局部極值。以本次實(shí)驗(yàn)為例，速度公式中的gbest1是以電法勘探數(shù)據(jù)通過多種群粒子群反演得到的全局最佳，gbest2是以探地雷達(dá)數(shù)據(jù)通過標(biāo)準(zhǔn)粒子群反演得到的全局最佳。由于電法勘探存在著等值現(xiàn)象，而探地雷達(dá)反演并不受其影響，因此gbest2可以準(zhǔn)確地反演出層厚度參數(shù)，在迭代過程中可以有效的避免因等值現(xiàn)象而使結(jié)果陷入局部極值的情況。

3 實(shí)驗(yàn)算法流程

圖1 算法流程圖

本次實(shí)驗(yàn)過程中，電法勘探參數(shù)下的多種群粒子群反演作為反演主程序，通過探地雷達(dá)參數(shù)下的標(biāo)準(zhǔn)粒子群反演求得gbest2id，參與主程序反演。根據(jù)模擬探測地層的深度，探地雷達(dá)正演模擬采用一階Mur吸收邊界的FDTD算法[15-19]，方式為位于地表的自激自收，激勵(lì)源采用主頻為100 MHz的布萊克曼-哈里斯脈沖。由于參數(shù)維度的不同，我們對(duì)電法勘探模型參數(shù)和探地雷達(dá)模型參數(shù)進(jìn)行處理，將兩種勘探方法所得到的參數(shù)分為電性參數(shù)和幾何參數(shù)。以n層地電斷面為例，對(duì)電法勘探參數(shù)和探地雷達(dá)參數(shù)進(jìn)行處理：電法參數(shù)model=(ρ1，ρ2，…，ρn與h1，h2，…，hn-1),ρ1、ρ2、…、ρn稱為電性參數(shù)ρs，h1、h2、…、hn-1稱為層厚度參數(shù)h1。同理，探地雷達(dá)參數(shù)model=(ε1,ε2,…εn,σ1,σ2,…σn,h1,h2,…h(huán)n-1)，ε1、ε2、…、εn與σ1、σ2、…、σn稱為電性參數(shù)εσ，h1、h2、…、hn-1稱為h2。

運(yùn)行反演主程序，當(dāng)多種群粒子群反演算法進(jìn)入迭代后，程序首先調(diào)用探地雷達(dá)參數(shù)下的粒子群反演程序，得到一組探地雷達(dá)模型參數(shù)下的速度v2id、位置x2id、個(gè)體極值點(diǎn)pbestid、全局極值點(diǎn)gbestid并返回給電法反演主程序。取此全局極值點(diǎn)中的層厚度h2，介電常數(shù)與電導(dǎo)率εσ，并用h2替換電法反演程序上次迭代(第一次為粒子初始化)所得的全局極值點(diǎn)中的層厚度h1，組成一個(gè)新的與電法勘探參數(shù)維度相同的全局極值點(diǎn)gbest2，并進(jìn)行反演。一次迭代過后，得到新的vid、xid、gbest1id。取各個(gè)數(shù)組中的層厚度參數(shù)h1，并用h1替換探地雷達(dá)程序上次迭代所得對(duì)應(yīng)數(shù)組中的層厚度h2，代入探地雷達(dá)粒子群反演程序，以此循環(huán)。

通過這種方式，層厚度參數(shù)在整個(gè)反演過程中不斷進(jìn)行著更新，矯正電法勘探因等值現(xiàn)象而出現(xiàn)的誤差，得到更準(zhǔn)確的地電模型。

表1 H型模型反演結(jié)果

4 算例分析

本次實(shí)驗(yàn)分析了水平地層上的H型和A型三層電測深曲線，曲線的基本形態(tài)由ρ1、ρ2、ρ3之間的大小關(guān)系決定。前兩層厚度為h1、h2，第三層為無限大。

4.1 H型三層曲線

H型三層曲線的電阻率關(guān)系為ρ1>ρ2<ρ3，且要滿足S等值現(xiàn)象的條件，因此給定理論模型model=(1 500，200，2 000，10，2)，分別進(jìn)行標(biāo)準(zhǔn)粒子群反演和多種群粒子群反演，對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行分析。

用標(biāo)準(zhǔn)粒子群反演程序?qū)碚撃Ｐ瓦M(jìn)行反演，取連續(xù)三次的模擬結(jié)果，記為model1、model2和model3。反演過程中，粒子群個(gè)數(shù)M=50，學(xué)習(xí)因子c1=c2=1.496 2，wmax=0.9，wmin=0.4。在最大迭代次數(shù)itermax=50的情況下均可達(dá)到最大誤差界限常數(shù)，記錄反演結(jié)果

給定與理論模型相同的地電斷面的探地雷達(dá)參數(shù)為model=(3，7，10，1/1500，1/200，1/2000，10，2)，其中第一層的相對(duì)介電常數(shù)ε1為3，電導(dǎo)率ρ1為1/1500；第二層的相對(duì)介電常數(shù)ε2為7，電導(dǎo)率ρ2為1/200；第三層的相對(duì)介電常數(shù)ε3為10，電導(dǎo)率ρ3為1/2000。用多種群粒子群反演程序進(jìn)行反演，結(jié)果記為model4。兩種反演算法的反演結(jié)果見表1，電法模型參數(shù)(即給定的理論模型)與反演結(jié)果對(duì)比見圖2，各模型參數(shù)深度與視電阻率對(duì)比結(jié)果見圖3。最后一次迭代后，探地雷達(dá)參數(shù)的模型參數(shù)與反演結(jié)果對(duì)比見圖4。

圖2 H型模型中模型參數(shù)與反演結(jié)果對(duì)比

圖3 深度與視電阻率對(duì)比

由圖4可見，從初始發(fā)射約116 ns第一層下底面反射波到達(dá)，151 ns第二層下底面反射波到達(dá)。其中，第一層單程走時(shí)為58 ns，電磁波速度為1.730 9×108m/s,相乘可得第一層層厚度為10.03 m; 第二層單程走時(shí)為18 ns，電磁波速度為1.133 1×108m/s,相乘可得第二層層厚度為19.83 m，與理論模型吻合，可以得到準(zhǔn)確的層厚度參數(shù)。

對(duì)比表1和圖1可以看出，在所給定的參數(shù)情況下，兩種反演方法求得模型參數(shù)的正演響應(yīng)曲線均與理論模型的正演曲線擬合，且都可以很好地求出ρ1、ρ3、h1的值。對(duì)比圖2可以看出model1、model2和model3均出現(xiàn)了S等值現(xiàn)象，縱向電導(dǎo)S2=h2/ρ2≈100，第二層的電阻率ρ2和層厚度h2與理論模型model數(shù)據(jù)相比同倍數(shù)增大；model4與理論模型相比，未出現(xiàn)等值現(xiàn)象，反演結(jié)果較為理想。

表2 A型模型反演結(jié)果

圖4 探地雷達(dá)參數(shù)的模型參數(shù)與反演結(jié)果對(duì)比

圖5 A型模型中模型參數(shù)與反演結(jié)果對(duì)比

圖6 深度與視電阻率對(duì)比圖

4.2 A型三層曲線

A型三層曲線的電阻率關(guān)系ρ1<ρ2<ρ3，給定理論模型model=(1000，2000，10000，10，2)，用標(biāo)準(zhǔn)粒子群反演程序?qū)碚撃Ｐ瓦M(jìn)行反演，取連續(xù)三次的模擬結(jié)果，記為model1、model2和model3。粒子群反演各參數(shù)設(shè)置均與H型三層相同。給定與理論模型相同的地電斷面的探地雷達(dá)參數(shù)為model=(3，7，10，1/1000，1/2000，1/10000，10，2)，用多種群粒子群反演程序進(jìn)行反演，結(jié)果記為model4。兩種反演算法的反演結(jié)果見表2，模型參數(shù)(即給定的理論模型)與反演結(jié)果對(duì)比見圖5，各模型參數(shù)深度與視電阻率對(duì)比結(jié)果見圖6。最后一次迭代后，探地雷達(dá)參數(shù)的模型參數(shù)與反演結(jié)果對(duì)比見圖7。

圖7 探地雷達(dá)參數(shù)的模型參數(shù)與反演結(jié)果對(duì)比

由圖7可知，由于給出的H型和A型三層地電模型的層厚度、相對(duì)介電常數(shù)、探地雷達(dá)正反演程序的各參數(shù)均相同，所以A型模型的單程走時(shí)、電磁波速度等均與H型基本相同，探地雷達(dá)參數(shù)的最后一次粒子群反演迭代結(jié)果與理論模型吻合，可以得到準(zhǔn)確的層厚度參數(shù)。

從表2和圖5可以看出model1、model2和model3均出現(xiàn)了S等值現(xiàn)象，縱向電導(dǎo)S2=h2/ρ2≈300，第二層的電阻率和層厚度與理論模型model數(shù)據(jù)相比同倍數(shù)增大；model4與理論模型相比，未出現(xiàn)等值現(xiàn)象，反演結(jié)果較為理想。

5 結(jié)論

筆者改進(jìn)了多種群粒子群反演算法，并將包含探地雷達(dá)信息的幾何參數(shù)加載到算法中，為原算法增加一個(gè)新的全局極值點(diǎn)，實(shí)現(xiàn)電性參數(shù)和幾何參數(shù)的共同反演，得到以下幾點(diǎn)結(jié)論：

1)當(dāng)理論模型或勘探結(jié)果滿足S等值現(xiàn)象的發(fā)生條件時(shí)，用標(biāo)準(zhǔn)粒子群算法進(jìn)行一維層狀電法反演，容易陷入局部極值，出現(xiàn)等值現(xiàn)象，無法得到準(zhǔn)確的地電模型。

2)多種群粒子群反演與標(biāo)準(zhǔn)粒子群反演相比，可以增加新的約束條件，不易陷入局部極值，取得更好的優(yōu)化結(jié)果。探地雷達(dá)參數(shù)反演求得的第二個(gè)全局極值點(diǎn)，可以更好地求出層厚度參數(shù)。相比于電法參數(shù)的標(biāo)準(zhǔn)粒子群反演，在整個(gè)迭代過程中層厚度參數(shù)不斷的進(jìn)行著更新，在解空間內(nèi)向著理論模型參數(shù)逼近，從而跳出局部極值點(diǎn)，避免了電法勘探的S等值現(xiàn)象。

3)多種群粒子群反演還有很多可以優(yōu)化改進(jìn)的地方，如c3的取值問題，c1、c2和c3之間的關(guān)系等還需進(jìn)一步完善，如何解決T等值現(xiàn)象還需進(jìn)一步改進(jìn)。