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

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

0 引 言

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

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

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

1.1 PSO算法原理

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

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

(1)

(2)

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

(3)

1.2 DE算法原理

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

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

(4)

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

(5)

式中：F為變異因子。

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

(6)

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

(7)

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

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

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

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

1.4 HIA工作步驟

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

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

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

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

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

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

1.5 HIA算法測試

為了測試HIA算法的算法性能，用Rastrigin函數(shù)對算法進行驗證。Rastrigin函數(shù)常用來評價優(yōu)化算法處理復雜問題的能力，此外，函數(shù)會因為維度不同而存在許多局部極小值點，但其全局極小值點只有0。函數(shù)的表達式如下：

(8)

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

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

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

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

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

2 理論模型的演算與分析

2.1 反演問題的實現(xiàn)

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

Φ=Φd+λΦm

(9)

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

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

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

(10)

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

(11)

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

Φm=‖m-mref‖2

(12)

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

(13)

2.2 四層地電模型

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

1)HIA反演分析。本文采用HIA算法對一維四層地電模型(KH)進行反演。首先對不含噪聲的模型進行反演，并把正演理論數(shù)據(jù)作為反演的觀察數(shù)據(jù)。為了模擬真實的觀測數(shù)據(jù)和驗證算法的魯棒性，筆者對理論數(shù)據(jù)分別添加5%和10%的隨機的噪聲進行反演，反演結(jié)果見表2。在HIA算法進行反演時，搜索范圍的上下限與真實值相差2倍或2倍以上。反演中選取的種群個數(shù)Np=60，K=100，交換次數(shù)τ=30，搜索算子ps=0.1，初始正則化因子λ=0.5。在進行無噪聲反演時，設(shè)定目標函數(shù)的停止條件為1×10-10。為了保證反演結(jié)果的可靠性，文中對觀測數(shù)據(jù)進行了80次反演，并把80次反演結(jié)果的平均值作為最終的結(jié)果。具體結(jié)果見表2。由表2的反演結(jié)果可知，HIA算法在不含噪的條件下能夠完全的反演出四層的地電模型，而且反演所需要的時間也很短。在含5%和10%的噪聲下，反演結(jié)果的最大誤差率分別為2.66%和5.16%。由此可以得出，噪聲對反演結(jié)果雖然有一定的影響，但是HIA算法也能夠反演出KH模型,可見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對比分析。為了彰顯改進算法的優(yōu)勢，分別將四層地電模型在不含噪聲的條件下與阻尼PSO算法和DE算法的反演結(jié)果進行了對比。同時，筆者為保證反演結(jié)果的可靠性，三種算法設(shè)置相同的參數(shù)，具體和式(1)中一致。3種算法分別獨立進行80次反演，并把反演結(jié)果的平均值作為輸出結(jié)果。表3給出了3種算法反演的對比結(jié)果，圖3分別給出了3種算法的反演迭代收斂對比曲線和反演結(jié)果對比曲線。由表3的反演結(jié)果可知，相同的迭代次數(shù)下，HIA算法幾乎可以恢復地下的模型，時間上也比DE反演和阻尼PSO反演較快。由圖3的反演結(jié)果可以得到，HIA算法在進行迭代17次就達到了停止條件，反演結(jié)果幾乎和理論模型重合，而阻尼PSO算法和DE算法反演的結(jié)果相對較差。由此可見，HIA算法比單獨的DE算法和阻尼PSO算法在大地電磁反演方面具有優(yōu)勢。

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

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

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

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

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

3 結(jié) 論

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

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

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