基于混沌天牛群算法的大地電磁反演

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

(成都理工大學 地球物理學院，成都 610059)

0 引言

大地電磁測深法(MT)作為地球物理勘探最主要的手段之一，其反演問題一直是研究熱點[1-2]。目前MT一維線性反演方法很多，如梯度法、差商牛頓法、馬奎特法和廣義逆矩陣法等，此類方法依賴于初始模型的選擇，全局尋優能力比較差，容易陷入局部最優解，魯棒性較差，有時甚至無法收斂[3-4]。鑒于此原因，國內、外學者對非線性反演展開了一系列研究，如模擬退火算法、重力搜索算法、遺傳算法、人工神經網絡及各種啟發式智能仿生算法等。柳建新等[5]將實數編碼遺傳算法應用于大地電磁二維理論模型的反演；師學明等[6]將阻尼粒子群算法應用于大地電磁的優化反演；胡祖志等[7]將人工魚群算法應用于大地電磁的最優化約束反演；王鶴等[8-10]基于神經網絡對大地電磁的非線性反演做了大量研究；Godio A等[11]再次將粒子群算法引入大地電磁反演，并與Occam反演方法進行了對比；王鵬飛等[12]在大地電磁的理論模型反演中引入改進布谷鳥算法。大量的研究表明，非線性反演方法在地球物理反演中應用前景十分廣闊。天牛須搜索算法(Beetle Antennae Search Algorithm,BAS)是理想化甲蟲類生物覓食和求偶的活動，并將其建立起數學模型的一種智能算法[13]。算法主要以天牛為例，天牛的頭部左右兩側具有兩根觸須，當天牛附近區域存在目標體(食物或者異性)時，其頭部的觸須檢測目標體所釋放的氣味，若右觸須感知到的氣味更強，天牛接下來往右側運動，反之，就往左側運動。根據天牛找尋目標體的這一生物活動，建立理想化數學模型，目標體或者說目標體的氣味相當于目標函數，然后利用此數學模型進行搜索尋優。

針對標準的BAS算法，一些學者提出天牛須搜索算法雖然能有目標地進行搜索，但是當搜索目標增多，尤其大于4個目標體時，收斂結果容易陷入局部極值，而且在多峰值函數優化時，陷入局部更容易陷入局部最優。趙玉強等[14]提出一類帶學習與競爭策略的混沌天牛群搜索算法(Learning and Competing Chaos Beetle Swarm Algorithm,LCCBSA)，這種改進的算法將一個天牛的搜索變成群體天牛的搜索，之后引入參考粒子群算法的指導性學習策略及競爭場策略。引入的策略使得天牛有目的的搜索且避免了陷入局部搜索的問題。這里利用LCCBSA算法對四個測試函數試算并與粒子群算法和遺傳算法的結果進行對比，該算法搜索效率和收斂速度相比其他算法有一定的優勢，于是進一步應用于大地電磁反演問題中。

1 大地電磁測深反演理論

大地電磁一維層狀模型視電阻率可用以下函數表示，即

ρa(f)=F(ρ1,ρ2…,ρn,h1,h2…,hn-1,f)

(1)

其中：ρi(ρ1,ρ2,…,ρn)表示模型電阻率參數;hi(h1,h2,…,hn-1)表示厚度參數;f表示頻率。大地電磁的反演的本質是根據實際觀測的ρa(f)求解模型參數ρ1、ρ2、…、ρn、h1、h2、…、hn-1。

將混沌天牛群算法應用于一維大地電磁的n層模型反演中，假設有m個天牛，相當于從m個初始模型出發，利用Tent映射將隨機向量(ρi,1,ρi,2,…,ρi,n,hi,1,hi,2,hi,n-1)T(i=1,2,…,m)對天牛群進行初始賦值，參數個數D=2n-1,并給出解空間每一個參數的數值大小范圍，然后利用混沌天牛群算法在解空間內搜索最優解。選取胡祖志等[15]所提出的大地電磁反演目標函數，如式(2)所示。

(2)

其中：N表示觀測頻點個數；前一項表示觀測視電阻率和計算視電阻率的誤差；后一項表示觀測相位和計算相位的誤差。

2 混沌天牛群搜索算法(LCCBSA)

2.1 標準BAS算法介紹

單個天牛個體朝著觸須感知到的適應度優勢的方向運動，最終搜索到全局的最優解。假設在解空間為D維，天牛的當前位置為X=(x1,x2,…,xD),則天牛左側和右側兩只觸須的搜索行為可用式(3)定義。

(3)

(4)

2.2 初始種群混沌化

在對群智能算法的研究過程中，很多文獻表明算法的收斂性與群體初始值的分布相關，初始群體分布如果具有更好的均勻性和遍歷性，算法的收斂性越好。筆者采用Tent映射得到天牛初始種群，Tent映射公式如式(5)。

(5)

式中：i=1、2、…、D表示混沌變量；n=1、2、…、N表示種群個體；μ表示混沌形態控制參數，一般情況下，其值在接近2時使得混亂狀態較好。為了得到均勻性和遍歷性較好的混沌初始種群，首先，種群中第一個天牛個體隨機生成，取值范圍為[0，1]，然后從第二個天牛個體開始初始值依次迭代產生，直到第N個天牛，最后再映射到解空間中。

2.3 指導性學習策略

在標準BAS算法中，天牛個體的搜索能力取決于觸須的長度，從當前位置向相鄰位置的移動并不容易，個體之間信息的交流限制較大。為了提升尋優能力，筆者參考粒子群算法個體位置更新的機制，引入指導性學習策略更新天牛個體位置。思路如下：群體中的天牛個體在當前位置會對目標體有一個感知方向和歷史群體最優個體的感知方向，當前個體感知目標的同時學習歷史群體最優值，向群體最優方向靠攏，到達圖中右側天牛位置(圖1)。策略將感知目標這一活動定義為認知部分，學習歷史群體最優這一活動定義為社會經驗部分。這種兼顧學習和自身判斷的策略，由于受到了認知部分和社會經驗部分的雙重影響，使得算法在多目標尋優中不會陷入局部最優。

圖1 天牛尋優示意圖

(6)

個體位置更新公式為式(7)。

(7)

(8)

2.4 競技場策略

在實際實驗中，群智能算法尋優在后期搜索時往往會陷入局部最優，難以收斂到全局最優解。Cheng R[16]提出優化大規模問題時可以采用競技的方法，蔣瑩瑩等[17]提高算法精度和收斂速度時采取了個體之間相互競爭的方式。筆者提出競技場策略，思路為將群體中個體的聚集度進行比較，進而競爭產生適應度高的個體，并對適應度小的個體進行位置更新。某代群體中個體聚集度的定義如式(9)所示。

(9)

其中，S表示兩只天牛個體的相似度，定義如下：

(10)

從式(9)中發現，天牛群體中任一個體都需要計算與其他個體的歐幾里得范數。閾值δ的值本文中取所得歐式范數的中位數，若其大于閾值δ，則兩個天牛個體之間無相似性；否則就表示兩者相似，相似度增加1，相似度累加值表示個體聚集度。顯然，聚集度高的個體越多，個體之間越相似，群體的多樣性越少，對算法的全局搜索不利，算法最終容易陷入局部極值。

為了解決群體中上述問題，將某代群體中聚集度的個體降序排序，選取前m個聚集度高個體參加競爭，通過計算它們的適應度，對其中n個適應度較小的個體位置更新。這一部分個體的位置更新采用偏好隨機的方式如式(11)所示。

(11)

2.5 算法流程

LCCBSA算法流程如圖2所示。

圖2 LCCBSA算法簡要流程圖

3 測試函數及大地電磁應用

3.1 測試函數試算

為了驗證本文提出的LCCBSA算法的性能，利用測試函數進行試算是最常用的手段。筆者選取4個多維基準測試函數，函數的基本信息如表1所示，函數表達式中D表示變量的維數。

表1 測試函數

在數值模擬實驗中，選擇已經較為成熟的粒子群算法(PSO)和遺傳算法(GA)與本文的LCCBSA算法進行對比，三種算法種群規模設置N=100，最大迭代次數200次，LCCBSA算法參數設置為：Cmax=2.5，Cmin=1.5，s=1.5，l=2.0，m=N/3，n=m/2；PSO算法參數設置為：慣性權重w=0.9，學習因子C1=C2=1.5；實數GA算法參數設置為：交叉概率Pc=0.8，變異概率Pm=0.1。表2為利用這三種算法對4個基準函數的分別重復進行30次實驗最優解的平均值對比。

表2 平均最優解對比

在相同測試函數和測試維度條件下，LCCBSA算法除F2測試函數外，其他三個測試函數較PSO算法和GA算法都有更好的尋優能力，且具有更快的收斂速度。

圖3 測試函數目標函數值對比

3.2 MT理論模型反演

3.2.1 三層(H型)理論模型反演

構建H型反演模型，反演過程中，迭代次數設置為2 000次，混沌天牛群其他參數設置與前面的測試函數一致，最小誤差ε=1×10-3，表3給出三層H型地電模型參數及LCCBSA算法反演結果，未添加噪聲正演響應曲線與反演計算曲線見圖4，兩者具有很好的擬合效果。

表3 三層H型地電模型LCCBSA算法反演結果

圖4 三層地電模型觀測數據與反演結果計算值對比

圖5 三層地電模型含噪10%觀測數據與反演結果計算值對比

3.2.2 六層理論模型反演

設計一個六層地電模型，模型參數、模型值和反演結果見表4。該模型為參考文獻[18]中計算大地電磁多尺度反演(MI)時所設計的模型，其主要是為了檢驗反演方法恢復地下深部低阻層的能力(如表4中9 450 m處的50 Ω·m的低阻地層)。LCCBSA算法反演參數設置同三層理論模型，未添加噪聲正演響應曲線與反演計算視電阻率曲線見圖6，由表4可以看出，MI反演結果的最大相對誤差達到44.6%，遠大于LCCBSA的最大相對誤差，說明LCCBSA反演結果較MI反演結果更好，而且對深部的低阻層電阻率恢復也更加令人滿意。

圖6 六層地電模型觀測數據與反演結果計算值對比

3.2.3 含噪理論數據分析

由于MT野外采集數據受到各種噪聲的影響，造成反演算法的不穩定性，為了檢驗本文算法的抗噪能力，對上述三層和六層理論正演計算結果數據添加10%的高斯噪聲，之后使用LCCBSA算法對含噪數據進行反演計算。三層模型和六層模型含噪數據反演結果分別見表3和表4，正演響應曲線與反演計算曲線分別見圖5和圖7。

從表3和表4中可以看出，理論模型含噪數據相對于無噪數據的反演結果，與真實模型值的差異增大，可見噪聲會對反演結果的準確性產生一定的影響，三層模型和六層模型參數的最大相對誤差分別為4.72%和14.85%，說明隨著模型的復雜度增加，添加噪聲后的數據反演結果差異也會變大。雖然噪聲一定程度上會影響反演結果的準確性，但是反演結果的相對誤差較小，與真實值相差不大，且圖5和圖7可以看出反演計算視電阻率曲線和相位曲線能與正演理論曲線擬合效果較好，說明含噪反演結果準確可靠，LCCBSA算法的抗噪性較好。

表4 六層地電模型LCCBSA算法反演結果對比表

圖7 六層地電模型含噪10%觀測數據與反演結果計算值對比

3.3 MT實測數據反演

應用新疆塔里木盆地布設的某條MT實測數據，驗證LCCBSA的反演效果。工區內第四系和新近系沉積大套礫石層，但其橫縱向分布范圍變化劇烈，造成地震資料反射凌亂，影響圈閉的準確落實。在實際生產工作中發現這套礫石層在電測井資料顯示為高阻特征，利用大地電磁刻畫可以取得較好的效果，且Bostick反演結果與地質規律吻合率較高。本文中選取的測線長為24.075 km，點距為500 m，頻點38個。對該測線進行分別進行LCCBSA反演和Bostick反演。

LCCBSA算法反演實測數據時，建立了首層厚度為50 m，遞增因子為1.1，總計37層的初始層厚模型，反演的初始電阻率模型為經過預處理后的電阻率，反演時保持層厚度不變，只對38個電阻率參數進行反演，電阻率搜索空間為經預處理后電阻率的50%～150%。LCCBSA算法計算所需要的天牛群數100個，迭代次數2 000次，認知因子Cmax=2.5，Cmin=1.5，學習因子s=1.5，天牛須長度l=2.0，迭代終止條件ε=1×10-3。反演時對測線上的每個測點獨立反演20次，求取20次結果的平均值作為每個測點反演結果。

圖8為bz1井旁測點的反演結果，圖8(a)顯示實際觀測數據與LCCBSA反演結果計算的視電阻率值變化趨勢一致，且擬合效果較好;圖8(b)為LCCBSA反演結果與Bostick反演結果和測井曲線的對比，可以發現LCCBSA和Bostick反演的電阻率與測井曲線具有很好地吻合，但在深部6 000 m處LCCBSA相對于Bostick反演電阻率偏大。

圖8 bz1井旁測點實測數據反演結果

圖9為實測電阻率和LCCBSA反演計算視電阻率擬斷面圖的對比，可以看出反演計算的視電阻率和實測電阻率的一致性較好，說明反演的電阻率較為可靠。圖10是LCCBSA和Botick擬二維反演解釋剖面圖的對比，兩種方法反演剖面具有一定的相似性，推斷解釋存在三套地層，第一層為高電阻率的第四系地層(Q)，電阻率一般在1 000 Ω·m左右，主要分布未成巖礫石；第二層為次高阻的新近系地層(N)，一般電阻率為幾十至幾百歐姆米，主要分布準成巖和成巖礫石；第三層低電阻率的古近系地層(E)和白堊系地層(K)，主要巖性為砂巖、泥巖。同時，LCCBSA算法反演剖面深部低阻范圍較Bostick反演剖面低阻范圍大，是由于初始層厚模型遞增變化導致深部層厚大，可以利用測井和地震解釋層位約束，反演出更為準確的低阻層。

圖9 實測視電阻率和反演計算視電阻率擬斷面圖對比

圖10 Bostick和LCCBSA算法反演剖面對比

4 結論

1)測試函數的試算表明，本文介紹的混沌天牛群搜索算法(LCCBSA)，在大部分情況下相對于粒子群算法和遺傳算法收斂速度較快，具有一定的優勢性。

2)MT理論模型和實測數據的試算結果表明，LCCBSA算法能夠較好地恢復地下介質的電阻率，是一種比較可行的非線性反演方法，避免了傳統線性反演方法初始模型的選取、陷入局部最優和計算雅可比偏導數矩陣的問題，并且可以用于航空電磁、瞬變電磁、地震等地球物理方法的反演。

3)LCCBSA算法跟其他非線性反演方法一樣，模型空間的選取、參數的改變都對反演效率有較大的影響，當前LCCBSA算法尚不成熟，應用于大地電磁的實驗是一種嘗試，并且其計算時間長的缺陷，極大地影響其向二、三維反演方向發展，但是將其和神經網絡結合是一個具有研究意義的方向。