基于混合智能算法優(yōu)化的大地電磁反演法研究

王 新，白寧波，周君君，王玉喜

(1.河南理工大學(xué) 物理與電子信息學(xué)院，河南 焦作 454000；2.中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北 武漢 430074； 3.河南景惠科技有限公司，河南 鄭州 450002)

0 引 言

大地電磁測(cè)深(Magnetotelluric，MT)是一種發(fā)展較早的地球物理勘探方法，該方法利用天然的電磁場(chǎng)作為場(chǎng)源進(jìn)行研究地球內(nèi)部的電性結(jié)構(gòu)，如今已廣泛應(yīng)用于固體礦產(chǎn)的勘察、石油天然氣的勘探等各個(gè)領(lǐng)域[1-2]。大地電磁反演是大地電磁測(cè)深的關(guān)鍵步驟之一，而反演算法的優(yōu)劣將直接影響最后的解釋結(jié)果。目前，大地電磁反演算法可分為線性算法和非線性算法，這2種算法的區(qū)別在于是否對(duì)目標(biāo)函數(shù)進(jìn)行線性化處理[3]。一維典型的線性反演算法包含Occam法，馬夸特法(Marquart)，正則化方法等[4-6]。但是反演問(wèn)題本身存在非線性、不適定性和解的非一性等問(wèn)題，因此線性化反演算法容易使反演陷入局部極小值，并且線性算法反演是否成功還依賴于對(duì)初始模型的選擇。為了克服線性算法的缺陷，一些學(xué)者將遺傳算法[7]、模擬退火法[8]、粒子群算法[9]、差分進(jìn)化算法[10-11]、人工魚(yú)群算法[12]、神經(jīng)網(wǎng)絡(luò)[13]和一些改進(jìn)的群智能算法[14-17]引入大地電磁的反演之中，并取得較好的效果。

目前，單一非線性算法針對(duì)大地電磁進(jìn)行反演時(shí)仍然存在收斂速度慢，求值精度低，尤其是隨著反演層數(shù)的增多，會(huì)導(dǎo)致反演的維度增加，這使得單一的非線性算法很難滿足反演精度和效率的需求。相關(guān)學(xué)者已經(jīng)結(jié)合其他算法對(duì)單一的非線性智能算法進(jìn)行改進(jìn)并應(yīng)用到反演中[18-20]。在綜合研究非線性算法發(fā)現(xiàn)，粒子群算法(Particle Swarm Optimization，PSO)和差分進(jìn)化算法(Differential Evolution，DE)盡管都有著良好的全局搜索能力，但是這2種算法也都有著自身的缺陷。如PSO算法極易陷入局部極小值而且尋優(yōu)結(jié)果不穩(wěn)定；DE算法容易陷入停滯和早熟狀態(tài)，而且在搜索速度上較PSO慢，同時(shí)這2種算法對(duì)局部的尋優(yōu)能力也相對(duì)較差。因此本文針對(duì)算法問(wèn)題，提出了一種混合智能算法(Hybrid Intelligent Algorithm，HIA)。HIA算法融合了PSO算法和DE算法對(duì)全局進(jìn)行迭代尋優(yōu)的搜索能力，讓兩種算法進(jìn)行優(yōu)勢(shì)互補(bǔ)，同時(shí)又利用Nelder-Mead優(yōu)化算法對(duì)尋優(yōu)結(jié)果進(jìn)行局部尋優(yōu)，使得算法能夠快速的收斂到最優(yōu)值，而且在進(jìn)行大地電磁反演時(shí)，引入正則化因子，可以有效地避免反演問(wèn)題的非唯一性。筆者詳細(xì)地?cái)⑹隽薍IA算法的工作原理，通過(guò)多維函數(shù)、地電模型以及實(shí)測(cè)數(shù)據(jù)對(duì)HIA算法進(jìn)行驗(yàn)證。

1 HIA算法的原理和實(shí)現(xiàn)

1.1 PSO算法原理

PSO算法自1995年由Kennedy和Eberhart提出以來(lái)，被廣泛的應(yīng)用于各領(lǐng)域。PSO算法在進(jìn)行迭代時(shí)，主要是利用速度和位置的更新來(lái)完成對(duì)適應(yīng)值函數(shù)的尋優(yōu)。

粒子群的速度和位置的更新公式如下：

(1)

(2)

PSO算法中，ω越大，粒子飛行的速度越大，粒子會(huì)以較大的步長(zhǎng)對(duì)全局進(jìn)行搜索，但是不利于收斂；反之，ω較小時(shí)，粒子的速度變小，位移變小，粒子相同迭代次數(shù)搜索范圍縮小，體現(xiàn)了粒子的局部搜索能力。c1和c2分別表示粒子自身的學(xué)習(xí)能力和粒子間的相互學(xué)習(xí)能力，因此慣性權(quán)重ω、c1和c2共同維護(hù)著粒子對(duì)全局和局部的平衡。為了保持算法的穩(wěn)定，文獻(xiàn)[21]利用式(3)來(lái)確定三者的關(guān)系：

(3)

1.2 DE算法原理

DE算法由Rainer Storn和Kenneth Price兩位學(xué)者在遺傳算法的基礎(chǔ)上提出來(lái)，基本的流程包含變異、交叉和選擇。關(guān)于差分進(jìn)化的具體描述如下：

1)初始化。在D維的空間中，對(duì)Np個(gè)個(gè)體進(jìn)行初始化，具體產(chǎn)生方式如下：

(4)

2)變異。從種群中隨機(jī)選取3個(gè)不同的個(gè)體，xr1、xr2、xr3(r1≠r2≠r3)，式(4)產(chǎn)生變異個(gè)體vi，具體如下：

(5)

式中：F為變異因子。

3)交叉。具體公式如下：

(6)

式中：CR∈(0,1)為交叉因子；rang(j)為[1,2,…，D]的隨機(jī)整數(shù)。

(7)

5)迭代停止。若種迭代次數(shù)達(dá)到最大迭代數(shù)K或者停機(jī)條件則尋優(yōu)結(jié)束，否則轉(zhuǎn)入步驟2)，開(kāi)始下一輪的迭代。

1.3 Nelder-Mead優(yōu)化算法

Nelder-Mead優(yōu)化算法是一種面向局部的搜索方法，而且有較高的搜索速率。其核心思想是：在n維的搜索空間中，隨機(jī)的生成有n+1個(gè)頂點(diǎn)相互連接生成簡(jiǎn)單的多面體。記xl為多面體的頂點(diǎn)，l=0,1,…,n，根據(jù)適應(yīng)值函數(shù)計(jì)算f(xl)，假設(shè)f(x0)≤f(x1)≤…≤f(xn)。其中，x0為最優(yōu)點(diǎn)，xn為最差點(diǎn)，然后經(jīng)過(guò)反射、擴(kuò)張、壓縮等操作尋找極小值點(diǎn)，關(guān)于頂點(diǎn)的更新的步驟如下：

5)令xq=x0+(xq-x0)，其中,q=0,1,…,n，形成新的多面體，重新計(jì)算f(xl)轉(zhuǎn)步驟1)。

1.4 HIA工作步驟

HIA算法是結(jié)合了PSO算法和DE算法的全局搜索能力和Nelder-Mead優(yōu)化算法的局部搜索能力而形成的一種的混合算法。HIA算法融合了2種算法對(duì)全局進(jìn)行尋優(yōu)的優(yōu)點(diǎn)，既能夠很大程度上避免了PSO算法容易陷入局部最小的缺點(diǎn)，又彌補(bǔ)了DE算法尋優(yōu)效率和容易早熟的不足，同時(shí)引入Nelder-Mead優(yōu)化算法進(jìn)行局部搜索，提高了算法整體的穩(wěn)定性和魯棒性，進(jìn)一步保證了算法尋優(yōu)的效率和結(jié)果的正確性。HIA算法工作步驟如下：

1)首先固定種群Np、最大迭代次數(shù)K、交換周期τ和搜索條件算子ps(0

2)利用目標(biāo)函數(shù)計(jì)算每個(gè)種群的適應(yīng)度函數(shù)值，判斷是否滿足停機(jī)條件。若滿足，保存最優(yōu)值，循環(huán)結(jié)束；若不滿足停機(jī)條件，則進(jìn)入下一步驟3)。

3)將Np/2種群利用PSO算法進(jìn)行迭代尋優(yōu)，判斷是否滿足r

4)剩余Np/2種群利用DE算法進(jìn)行迭代尋優(yōu)，然后判斷是否滿足r

5)判斷算法是否迭代了τ次，若滿足，將種群Np進(jìn)行重組，然后再將Np隨機(jī)均勻分配成2組，轉(zhuǎn)到步驟2)；若不滿足，直接轉(zhuǎn)步驟2)。

1.5 HIA算法測(cè)試

為了測(cè)試HIA算法的算法性能，用Rastrigin函數(shù)對(duì)算法進(jìn)行驗(yàn)證。Rastrigin函數(shù)常用來(lái)評(píng)價(jià)優(yōu)化算法處理復(fù)雜問(wèn)題的能力，此外，函數(shù)會(huì)因?yàn)榫S度不同而存在許多局部極小值點(diǎn)，但其全局極小值點(diǎn)只有0。函數(shù)的表達(dá)式如下：

(8)

式中：j=1,2,…D為函數(shù)的維度；對(duì)于R(X)，xj∈(-5.12,5.12)。

為了彰顯混合智能算法的優(yōu)勢(shì)，將HIA算法尋優(yōu)得到的結(jié)果與文獻(xiàn)[9]中提出的阻尼PSO算法和DE算法的尋優(yōu)結(jié)果進(jìn)行對(duì)比。為保證算法結(jié)果對(duì)比的可靠性，3種算法設(shè)置共同參數(shù)，最大迭代次數(shù)K=1 000,種群Np=60，ps=0.1，停止條件為1×10-20。在進(jìn)行HIA迭代時(shí)，交換周期τ=250,函數(shù)的維度D=15。考慮到算法的隨機(jī)性，該算法獨(dú)立運(yùn)行80次的平均值作為最終輸出結(jié)果，具體見(jiàn)表1。

表1 Rastrigin函數(shù)的尋優(yōu)結(jié)果對(duì)比Table 1 Comparison of optimization results of Rastrigin function

由表1可得，對(duì)于Rastrigin函數(shù)的尋優(yōu)，HIA算法可以準(zhǔn)確的求得函數(shù)的全局極小值，而阻尼PSO算法和DE算法都沒(méi)有收斂到全局極小值點(diǎn)。圖1給出了3種算法的迭代下降曲線。由圖1可得，HIA算法會(huì)在迭代800多次時(shí)就達(dá)到停機(jī)條件，而阻尼PSO算法，直接陷入了局部極小值；DE算法的迭代曲線雖然仍在下降，但是算法的收斂速度和尋優(yōu)結(jié)果都不理想。因此，HIA算法無(wú)論是尋優(yōu)效率還是尋優(yōu)結(jié)果都明顯優(yōu)于阻尼PSO算法和DE算法。

圖1 Rastrigin函數(shù)的迭代下降曲線Fig.1 Iterative descent curves of Rastrigin function

2 理論模型的演算與分析

2.1 反演問(wèn)題的實(shí)現(xiàn)

在地球物理反演中，目標(biāo)函數(shù)通常采用以下形式：

Φ=Φd+λΦm

(9)

式中：Φd是數(shù)據(jù)的擬合差；Φm是模型的擬合差；λ為正則化因子；m是模型向量，引入模型擬合差可以有效避免反演的多解性。

對(duì)于一維r層地電斷面而言，模型向量m表示為

m=(ρ1，ρ2，…ρr，h1,h2,…h(huán)r-1)

(10)

一維反演數(shù)據(jù)擬合差Φd表達(dá)式為

(11)

根據(jù)不同的約束條件，模型擬合差Φm的約束有很多種，常用的有3種，模型參數(shù)的平方和最小，模型參數(shù)導(dǎo)數(shù)的平方和最小和模型參數(shù)二階導(dǎo)數(shù)的平方和最小。在此選取最小模型的約束函數(shù)，式中，mref為背景模型。

Φm=‖m-mref‖2

(12)

對(duì)于正則化因子的選取，本文采取文獻(xiàn)[5]提出的自適應(yīng)調(diào)節(jié)正則化因子的調(diào)節(jié)方法，該方法在一定程度上同時(shí)滿足了使Φd和Φm極小的條件，而且能很好的擬合數(shù)據(jù)和壓制模型的冗余構(gòu)造。具體公式如下：

(13)

2.2 四層地電模型

文中數(shù)值模型反演程序在Intel(R)Core(TM)i5-7200 2.5GHz(四核)上運(yùn)行。

1)HIA反演分析。本文采用HIA算法對(duì)一維四層地電模型(KH)進(jìn)行反演。首先對(duì)不含噪聲的模型進(jìn)行反演，并把正演理論數(shù)據(jù)作為反演的觀察數(shù)據(jù)。為了模擬真實(shí)的觀測(cè)數(shù)據(jù)和驗(yàn)證算法的魯棒性，筆者對(duì)理論數(shù)據(jù)分別添加5%和10%的隨機(jī)的噪聲進(jìn)行反演，反演結(jié)果見(jiàn)表2。在HIA算法進(jìn)行反演時(shí)，搜索范圍的上下限與真實(shí)值相差2倍或2倍以上。反演中選取的種群個(gè)數(shù)Np=60，K=100，交換次數(shù)τ=30，搜索算子ps=0.1，初始正則化因子λ=0.5。在進(jìn)行無(wú)噪聲反演時(shí)，設(shè)定目標(biāo)函數(shù)的停止條件為1×10-10。為了保證反演結(jié)果的可靠性，文中對(duì)觀測(cè)數(shù)據(jù)進(jìn)行了80次反演，并把80次反演結(jié)果的平均值作為最終的結(jié)果。具體結(jié)果見(jiàn)表2。由表2的反演結(jié)果可知，HIA算法在不含噪的條件下能夠完全的反演出四層的地電模型，而且反演所需要的時(shí)間也很短。在含5%和10%的噪聲下，反演結(jié)果的最大誤差率分別為2.66%和5.16%。由此可以得出，噪聲對(duì)反演結(jié)果雖然有一定的影響，但是HIA算法也能夠反演出KH模型,可見(jiàn)HIA具有良好的魯棒性。圖2給出了HIA算法反演在添加10%的噪聲的視電阻率擬合曲線、相位擬合曲線。

表2 四層地電模型HIA算法反演結(jié)果Table 2 Inversion results of HIA algorithm for four-layer geoelectric model

圖2 四層含10%噪聲地電模型HIA算法反演擬合結(jié)果Fig.2 HIA algorithm for inverse fitting of four-layer geoelectric model with 10% noise

2)HIA與PSO和DE對(duì)比分析。為了彰顯改進(jìn)算法的優(yōu)勢(shì)，分別將四層地電模型在不含噪聲的條件下與阻尼PSO算法和DE算法的反演結(jié)果進(jìn)行了對(duì)比。同時(shí)，筆者為保證反演結(jié)果的可靠性，三種算法設(shè)置相同的參數(shù)，具體和式(1)中一致。3種算法分別獨(dú)立進(jìn)行80次反演，并把反演結(jié)果的平均值作為輸出結(jié)果。表3給出了3種算法反演的對(duì)比結(jié)果，圖3分別給出了3種算法的反演迭代收斂對(duì)比曲線和反演結(jié)果對(duì)比曲線。由表3的反演結(jié)果可知，相同的迭代次數(shù)下，HIA算法幾乎可以恢復(fù)地下的模型，時(shí)間上也比DE反演和阻尼PSO反演較快。由圖3的反演結(jié)果可以得到，HIA算法在進(jìn)行迭代17次就達(dá)到了停止條件，反演結(jié)果幾乎和理論模型重合，而阻尼PSO算法和DE算法反演的結(jié)果相對(duì)較差。由此可見(jiàn)，HIA算法比單獨(dú)的DE算法和阻尼PSO算法在大地電磁反演方面具有優(yōu)勢(shì)。

表3 四層地電模型不同算法反演結(jié)果對(duì)比Table 3 Comparison of inversion results of different algorithms for four-layer geoelectric model

圖3 四層地電模型不同算法反演的對(duì)比曲線Fig.3 Contrast curves inverted by different algorithms for four-layer geoelectric model

2.3 實(shí)測(cè)數(shù)據(jù)反演

為進(jìn)一步驗(yàn)證算法的可行性，利用HIA算法對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行反演。實(shí)測(cè)數(shù)據(jù)來(lái)自文獻(xiàn)[5],取自吉林樺皮廠地?zé)崽锏囊粋€(gè)實(shí)測(cè)點(diǎn)數(shù)據(jù)，由鉆井資料可知，地下電阻率大概呈現(xiàn)“高低高低”的縱向分布。采用HIA算法反演的參數(shù)設(shè)置為Np=60，K=400，交換次數(shù)τ=100，單純形搜索算子ps=0.1，初始正則化因子λ=0.5。在算法進(jìn)行反演時(shí)，采用一個(gè)六層的地電模型進(jìn)行擬合反演，并將反演結(jié)果和目前常用的實(shí)測(cè)反演算法Bostick反演以及Occam的反演的結(jié)果進(jìn)行了對(duì)比，圖4給出了3種算法的視電阻率和相位的擬合曲線。由圖4a的視電阻率曲線的擬合結(jié)果可知，Occam反演結(jié)果和本文HIA算法的反演結(jié)果效果較佳，但是Occam的高頻部分?jǐn)M合較差，而本文算法高低頻的擬合效果都很好；由圖4b相位曲線的擬合結(jié)果可以得出，Bostick和本文算法反演結(jié)果較好，但是Bostick在高低頻的擬合效果較好，中頻的擬合效果較差，而HIA算法雖然低頻效果較Bosick的結(jié)果差，總體效果卻是最佳的。采用式(11)計(jì)算3種算法的擬合差，Occam的結(jié)果為8.78，Bostick的結(jié)果為2.36，而HIA算法結(jié)果為0.781 5。由此可得，本文算法用于實(shí)測(cè)大地電磁的反演也是十分有效的。

圖4 實(shí)測(cè)數(shù)據(jù)反演結(jié)果對(duì)比Fig.4 Comparison of measured data inversion results

3 結(jié) 論

1)提出了混合HIA算法解決了PSO算法和DE算法尋優(yōu)效率低和易陷入局部極小值的問(wèn)題。通過(guò)對(duì)多維函數(shù)和理論模型進(jìn)行對(duì)比試驗(yàn)，證明了HIA算法比這2種算法反演的尋優(yōu)效率更高。

2)利用HIA算法對(duì)實(shí)測(cè)點(diǎn)數(shù)據(jù)進(jìn)行反演，反演結(jié)果表明，HIA反演結(jié)果與實(shí)測(cè)數(shù)據(jù)更加吻合。

3)考慮到二維或三維更符合實(shí)際的地質(zhì)情況，將來(lái)可以對(duì)HIA算法進(jìn)一步研究，以便能夠拓展到高維反演。