改進的布谷鳥MT反演算法

王鵬飛 王書明

(中國地質大學(武漢)地球物理與空間信息學院，湖北武漢 430074)

0 引言

大地電磁(Magnetotelluric，MT)勘探以天然電磁場為場源，通過測量地表相互正交的不同頻率的電磁場分量，研究地球內部電性結構[1-3]。該方法具有探測深度大(可從地表至地下上千公里深度范圍)、不需要人工場源、對低阻體敏感等特點，因而被廣泛應用于深部地殼、地幔電性結構研究、油氣資源普查、工程勘察、環境監測等領域[4-5]。

MT數據的反演研究經歷了三個發展階段。最早使用近似反演類方法，如Bostick反演[6]、大地電磁擬地震反演[7]、曲線對比法[8]等； 隨后發展了各種線性反演方法，如馬奎特法[9]、廣義逆矩陣法[10]、奧康姆法[11-12]、正則化法[13]等。由于MT反演是高度非線性的，線性反演方法容易陷入局部極小值，且初始模型的選擇對反演結果影響很大[14]，目前的反演研究轉向非線性方法和仿生算法，如模擬退火法[15]、遺傳算法[16]、人工神經網絡[17]、粒子群算法[18-19]、人工魚群算法[20]等。

布谷鳥搜索(Cuckoo search algorithm，CS)算法是受自然界中布谷鳥寄生育雛習性的啟發而提出的，該算法采用萊維飛行模式更新個體，具有很強的全局搜索能力[21-23]。該算法一經提出，備受國內外學者的關注，并被不斷改進以獲得更好的優化效果。鄭洪清[24]采用自適應的發現概率和步長，加快算法的收斂速度； 李明等[25]提出一種結合差分進化算法的混合布谷鳥算法。

本文將粒子群算法中的全局最優解思想引入布谷鳥算法，加強算法的局部搜索能力；同時結合單純形法，在每次迭代后，利用反射、擴張、壓縮等手段改進較差的鳥巢，進一步提高算法的尋優精度。為了評價改進后算法對MT反演的有效性和穩定性，首先對三個典型模型的視電阻率和相位曲線進行聯合反演試算，然后加入隨機噪聲分析算法的抗噪能力，最后反演了來自中國內蒙古的實測數據。

1 改進布谷鳥算法原理

1.1 基本布谷鳥算法

自然界中布谷鳥不會筑巢，也不會孵化自己的卵，而是將卵產在其他鳥類的鳥巢中，由它們代為孵化。而其他鳥類不愿意孵化外來的卵，當它們發現外來的卵就會將其扔掉，有時甚至拋棄整個鳥巢，另筑新巢。為了不被發現，布谷鳥會根據自己的卵的形狀和顏色選擇宿主[26-28]。

布谷鳥算法規定一個布谷鳥只在一個鳥巢中產一個蛋，設總的鳥巢個數為n，用模型xi(i=1,2,…，n)表示每個鳥巢的位置，其目標函數值表示鳥巢的優劣，每個模型xi中包括d個參數xij(j=1,2,…,d)，在這些參數范圍內隨機產生n個初始模型。

布谷鳥算法有兩種方式更新鳥巢的位置：一是萊維飛行模式，另一種是隨機搜索。

萊維飛行模型是模擬飛行動物覓食的過程，是一種隨機游走的方式，其軌跡如圖1所示。在飛行過程中，大量的短距離搜索與偶爾的長距離跳躍交替出現，短距離搜索有利于提高算法的局部搜索能力，而偶爾的長距離跳躍使得優化算法不易陷入局部極值[29-30]。

圖1 萊維飛行軌跡示意圖

萊維飛行公式為

(1)

(2)

式中：u、v是服從標準正態分布的隨機變量；φ的取值為

(3)

其中Γ(·)是伽馬函數。

另一種更新方式是對每個鳥巢都產生一個[0,1]的隨機數ε，與宿主發現布谷鳥鳥蛋的概率P做比較，通過下式判斷是否產生新鳥巢

(4)

1.2 改進的布谷鳥(ICS)算法

由于采用萊維飛行更新個體，布谷鳥算法具有很強的探索能力，但開發能力稍顯不足。本文借鑒粒子群中的全局最優解思想，對鳥巢位置的更新(式(1))進行改進

(5)

反射操作：xr=xc+δr(xc-xs),其中xr為反射點，反射系數δr通常取1.0。

擴張操作：xe=xc+δe(xr-xc),其中xe為擴張點，擴張系數δe通常取2.0。

收縮操作：xw=xc+δw(xs-xc),其中xw為收縮點，收縮系數δw的取值一般等于壓縮系數δt。

單純形法(圖2)的計算步驟如下：

(1)通過比較各個鳥巢的目標函數值，得到最優鳥巢xg和次優鳥巢xb的目標函數值f(xg)和f(xb)。

(2)找到最差鳥巢位置的目標函數值f(xs)，對xs執行反射操作，得到反射點xr。

(3)若f(xr)

(4)若f(xr)>f(xs)，說明位置變得更差，則執行壓縮操作，得到壓縮點xt;若f(xt)

信息化教學是目前各級各類學校普遍應用的教學方式，因為通過信息化的手段可以使師生輕松獲得更多優質教學資源，改善教學方法、提高教學效率、拓展知識視野，真正實現人人處處時時皆可學習。在信息化教學中，考試也逐漸從單一的筆試慢慢過渡到筆試+機試等多樣化的考試方式，通過計算機考試、評閱可以最大限度保證閱卷的公平性，降低誤判、錯判的幾率，節省人力、物力，克服了人為差錯的可能性。目前，很多考試、測驗系統中客觀題均實現了自動評閱，但主觀題的自動評閱一直是瓶頸問題，本文重點研究了主觀題自動評閱在信息化教學中的重要意義及可實施性。

(5)若f(xg)

圖2 單純形法結構示意圖

1.3 算法測試

為了檢驗改進算法的尋優性能，選擇Ackely函數進行實驗，主要測試算法的優化精度，并與基本布谷鳥算法進行對比。Ackely函數在(x1,x2,…,xd)=(0,0,…,0)處取最小值0，其中d為函數自變量的個數，本次實驗d取6。此函數是典型的多極值函數[32]，其表達式為

xi∈[-32,32]

(6)

分別應用ICS和CS算法對上述函數求取最小值，迭代次數為100次。表1列出了兩種算法的優化結果，每次結果均為6次計算結果的平均值。從表1可以看出，改進后的算法在優化精度上有了明顯的提高。

表1 兩種算法性能比較

2 理論模型試算

為了檢驗ICS算法的有效性和穩定性，對三個典型地電模型進行反演試算。模型A為二層G型，模型B為三層H型，模型C為四層KH型，模型參數分別如表2所示。在實際應用中，通常沒有足夠的信息估計地層的電阻率和厚度。為了更符合實際情況，使用較寬的模型搜索范圍，規定搜索范圍的下限和上限與真實值相差2倍以上。為了避免單次反演的偶然性，對每個模型都進行10次獨立反演，并用這10次反演的平均值作為最終的反演結果。反演參數設置為：鳥巢個數為20，發現外來鳥蛋的最小概率為Pmin=0.1，最大概率為Pmax=0.8，最大迭代次數為200。

MT反演是求解目標函數最小值的優化問題，由于視電阻率可以直觀地反映介質的電阻率，因此對視電阻率數據進行反演的情況居多，而相位數據的反演則相對較少。本文同時對視電阻率和相位進行聯合反演，目標函數公式為

(7)

式中：N為反演的頻點數；ρcal,i和ρobs,i分別為第i個頻率的視電阻率計算值和觀測值；θcal,i和θobs,i分別為第i個頻率的相位計算值和觀測值。

表2 不同地電模型參數及反演搜索范圍

2.1 無噪聲理論模型

首先，應用ICS算法對無噪聲的理論數據進行反演，三個模型的反演結果分別如表3和圖3～圖5所示。由反演結果可以看出，使用較寬模型搜索范圍時，反演的視電阻率和相位曲線與理論曲線基本重合。三個模型所有參數反演結果的最大相對誤差分別為0、0、0.66%，即反演結果與模型十分吻合，各個模型的真實值均被ICS算法較準確地反演和重建?？梢奍CS算法可以對理論大地磁場數據進行精確反演。

圖3 模型A不含噪聲理論數據與反演結果對比

圖4 模型B不含噪聲理論數據與反演結果對比

圖5 模型C不含噪聲理論數據與反演結果對比

2.2 含噪聲理論數據

在野外采集MT數據時，不可避免地會接收到各種噪聲干擾，這會導致反演算法的不穩定，可能會得到錯誤的反演結果，因此有必要對本文ICS算法的抗噪能力進行測試。在前面三種模型的理論視電阻率和相位數據中都加入10%的高斯隨機噪聲，使用ICS算法對含噪數據進行反演計算，反演結果見圖6～圖8?？梢钥闯?，與不含噪聲相比，含噪數據的反演結果與真實模型之間存在一定偏差，三個模型中參數(電阻率或者厚度)最大相對誤差分別增至2.05%、3.44%、5.33%，因此噪聲對反演結果的準確性有一定影響。但是反演視電阻率和相位曲線依然能夠很好地擬合理論曲線，計算的地層電阻率和厚度與實際模型基本吻合，說明ICS具有較好的抗噪性。

2.3 反演結果分析

反演結果的穩定性和收斂性是評價算法的重要標準。利用ICS算法對每個模型都行了10次獨立反演，反演結果中沒有出現與真實模型相差很大的情況，即反演成功率為100%。

圖6 模型A含10%白噪聲理論數據與反演結果對比

圖7 模型B含10%白噪聲理論數據與反演結果對比

圖8 模型C含10%白噪聲理論數據與反演結果對比

為了定量分析ICS算法的穩定性，統計了反演參數的相對誤差和標準差。結果顯示，對于無噪數據，模型的反演結果與真實值幾乎完全一致，相對誤差和標準差都趨近于零，說明反演結果準確可靠。對于含噪數據，地層厚度的相對誤差最大不超過5.33%，標準差最大不超過32.08m，說明噪聲在一定程度上影響了反演結果的精度，但所得解依然穩定，且與真實值相差不大，反演結果仍令人滿意。

為了評價ICS算法的收斂性，以模型C為例，對反演過目標函數的收斂曲線進行分析?？梢钥闯?，無論含噪數據還是不含噪數據，目標函數在前20次迭代中快速收斂，隨后逐漸收斂到零值(黑線)或某一常數(紅線)，這表明ICS算法完成了對最優解的搜索。由此可見，在無噪聲和含10% 噪聲的情況下，ICS算法對MT數據的反演均具有良好的收斂性。

2.4 ICS與CS算法的對比分析

為了進一步分析改進后算法的收斂性能，分別應用ICS和CS算法對模型C不含噪聲的理論數據進行反演試算。使用表2中的模型搜索范圍，布谷鳥群個數為20，最大迭代次數為200。每種算法都進行10次獨立反演，取其平均值作為最終的反演結果，反演結果如表3和表4所示。目標函數收斂曲線如圖9b所示。

表3 模型A～模型C含噪與不含噪數據ICS法反演結果統計表

從圖9b可以看出，在200次迭代內，ICS比CS算法具有更好的收斂能力，且目標函數更趨近于零。為了定量分析，對ICS和CS算法反演結果的相對誤差進行了統計(表3和表4)，結果顯示兩種算法反演的最大相對誤差分別為0.06%和1.13%，最大標準差分別為1.47m和16.82m，說明ICS反演更加穩定、反演結果的精度更高。在耗時方面，迭代反演200次CS算法用時21.227s，而ICS算法為22.674s，用時稍多，這是由于ICS算法中加入了單純形法的反射、壓縮、擴張操作，增加了反演時間，但增加時間并不太長，而在收斂性能和尋優精度上的提升效果更明顯。

表4 模型C不含噪聲數據CS法反演結果統計表

圖9 模型C反演目標函數收斂曲線

3 實測數據反演

應用中國內蒙古某地的實測MT數據驗證ICS算法的反演效果。測區內布設了8條測線，共130個測點，點距為500m。分別使用WinGlink軟件中的OCCAM方法和本文的ICS算法對4號線第15個測點的數據進行反演。

圖10a為ICS和OCCAM方法反演得到的該測點地層電阻率。由圖中可見，兩種算法反演的電阻率變化趨勢大致相同。從地表向下大致可以分為三個電性層：第一層為低阻層，電阻率僅十幾歐姆米，厚度約為210m；第二層為次高阻層，電阻率為數十到上千歐姆米，厚度為2700m；第三層為高阻層，電阻率最大值達17600Ω·m。

圖10 內蒙古MT數據兩種方法反演結果對比

由圖10b可以看出，OCCAM和ICS反演模型的視電阻率曲線與實測數據都基本吻合，但是ICS反演結果與真實地層更加接近，即ICS算法反演結果擬合得更好。

4 結論

本文將布谷鳥算法應用到MT數據反演中。通過引入粒子群中的全局最優解引導思想，提高算法的局部搜索能力，并結合單純形法，進一步提高算法的尋優精度。首先采用無噪聲和含10%噪聲的理論數據對ICS算法的有效性和穩定性進行了試驗；然后將ICS和CS算法進行比較，分析了改進前后布谷鳥算法的收斂性能；最后反演了來自內蒙古某地的實測MT數據。

理論和實測數據反演結果表明：

(1)ICS算法是一種性能優異的群智能優化算法，可以有效地對MT數據進行反演，且反演結果穩定、精度較高;

(2)與傳統的布谷鳥算法相比，ICS算法收斂能力更強、尋優精度更高，但反演用時稍多。