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

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

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

0 引言

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

針對(duì)標(biāo)準(zhǔn)的BAS算法，一些學(xué)者提出天牛須搜索算法雖然能有目標(biāo)地進(jìn)行搜索，但是當(dāng)搜索目標(biāo)增多，尤其大于4個(gè)目標(biāo)體時(shí)，收斂結(jié)果容易陷入局部極值，而且在多峰值函數(shù)優(yōu)化時(shí)，陷入局部更容易陷入局部最優(yōu)。趙玉強(qiáng)等[14]提出一類帶學(xué)習(xí)與競爭策略的混沌天牛群搜索算法(Learning and Competing Chaos Beetle Swarm Algorithm,LCCBSA)，這種改進(jìn)的算法將一個(gè)天牛的搜索變成群體天牛的搜索，之后引入?yún)⒖剂Ｗ尤核惴ǖ闹笇?dǎo)性學(xué)習(xí)策略及競爭場策略。引入的策略使得天牛有目的的搜索且避免了陷入局部搜索的問題。這里利用LCCBSA算法對(duì)四個(gè)測試函數(shù)試算并與粒子群算法和遺傳算法的結(jié)果進(jìn)行對(duì)比，該算法搜索效率和收斂速度相比其他算法有一定的優(yōu)勢，于是進(jìn)一步應(yīng)用于大地電磁反演問題中。

1 大地電磁測深反演理論

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

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

(1)

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

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

(2)

其中：N表示觀測頻點(diǎn)個(gè)數(shù)；前一項(xiàng)表示觀測視電阻率和計(jì)算視電阻率的誤差；后一項(xiàng)表示觀測相位和計(jì)算相位的誤差。

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

2.1 標(biāo)準(zhǔn)BAS算法介紹

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

(3)

(4)

2.2 初始種群混沌化

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

(5)

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

2.3 指導(dǎo)性學(xué)習(xí)策略

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

圖1 天牛尋優(yōu)示意圖

(6)

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

(7)

(8)

2.4 競技場策略

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

(9)

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

(10)

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

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

(11)

2.5 算法流程

LCCBSA算法流程如圖2所示。

圖2 LCCBSA算法簡要流程圖

3 測試函數(shù)及大地電磁應(yīng)用

3.1 測試函數(shù)試算

為了驗(yàn)證本文提出的LCCBSA算法的性能，利用測試函數(shù)進(jìn)行試算是最常用的手段。筆者選取4個(gè)多維基準(zhǔn)測試函數(shù)，函數(shù)的基本信息如表1所示，函數(shù)表達(dá)式中D表示變量的維數(shù)。

表1 測試函數(shù)

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

表2 平均最優(yōu)解對(duì)比

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

圖3 測試函數(shù)目標(biāo)函數(shù)值對(duì)比

3.2 MT理論模型反演

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

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

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

圖4 三層地電模型觀測數(shù)據(jù)與反演結(jié)果計(jì)算值對(duì)比

圖5 三層地電模型含噪10%觀測數(shù)據(jù)與反演結(jié)果計(jì)算值對(duì)比

3.2.2 六層理論模型反演

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

圖6 六層地電模型觀測數(shù)據(jù)與反演結(jié)果計(jì)算值對(duì)比

3.2.3 含噪理論數(shù)據(jù)分析

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

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

表4 六層地電模型LCCBSA算法反演結(jié)果對(duì)比表

圖7 六層地電模型含噪10%觀測數(shù)據(jù)與反演結(jié)果計(jì)算值對(duì)比

3.3 MT實(shí)測數(shù)據(jù)反演

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

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

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

圖8 bz1井旁測點(diǎn)實(shí)測數(shù)據(jù)反演結(jié)果

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

圖9 實(shí)測視電阻率和反演計(jì)算視電阻率擬斷面圖對(duì)比

圖10 Bostick和LCCBSA算法反演剖面對(duì)比

4 結(jié)論

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

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

3)LCCBSA算法跟其他非線性反演方法一樣，模型空間的選取、參數(shù)的改變都對(duì)反演效率有較大的影響，當(dāng)前LCCBSA算法尚不成熟，應(yīng)用于大地電磁的實(shí)驗(yàn)是一種嘗試，并且其計(jì)算時(shí)間長的缺陷，極大地影響其向二、三維反演方向發(fā)展，但是將其和神經(jīng)網(wǎng)絡(luò)結(jié)合是一個(gè)具有研究意義的方向。