雙種群人工蜂群算法及其在MT和重力聯(lián)合反演中的應(yīng)用

曾志文 陳 曉* 郭 冬 鄧居智 張志勇 陳 輝

(①東華理工大學(xué)核資源與環(huán)境國家重點(diǎn)實(shí)驗室，江西南昌 330013；②東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌 330013；③安徽省勘查技術(shù)院，安徽合肥 230041)

0 引言

地球物理單一方法的局限性、反演的多解性以及地質(zhì)條件的復(fù)雜性使得聯(lián)合反演成為地球物理領(lǐng)域的必然發(fā)展趨勢[1]。聯(lián)合反演是一種定量的綜合地球物理解釋技術(shù)，可以減少解的非唯一性。隨著地球物理勘探向地球深部發(fā)展，聯(lián)合反演的優(yōu)勢也越加突顯[2-4]。基于線性優(yōu)化算法的聯(lián)合反演收斂速度快，但容易陷入局部極小。基于非線性優(yōu)化算法的聯(lián)合反演具有全局尋優(yōu)、無需求偏導(dǎo)數(shù)矩陣、便于先驗信息融入等特點(diǎn)。目前，很多非線性優(yōu)化算法在地球物理聯(lián)合反演領(lǐng)域得到了廣泛應(yīng)用，如模擬退火算法[5-7]、遺傳算法[8-9]、差分進(jìn)化算法[10-11]、人工魚群算法[12]等。可見基于非線性優(yōu)化算法的聯(lián)合反演是地球物理聯(lián)合反演領(lǐng)域的一個重要發(fā)展方向。

人工蜂群(ABC)算法是Karaboga[13]于2005年提出的一種非線性優(yōu)化算法。在該算法中，存在引領(lǐng)蜂、跟隨蜂和偵察蜂，每種蜜蜂各司其職，并存在獨(dú)特的角色轉(zhuǎn)換機(jī)制。由于該算法是一種較新穎的優(yōu)化算法，相比傳統(tǒng)的非線性優(yōu)化算法(如模擬退火、遺傳算法等)，該算法在地球物理反演領(lǐng)域的應(yīng)用相對偏少。王猛等[14]將ABC算法應(yīng)用于瞬變電磁測深資料反演，并指出ABC算法可以提高瞬變電磁資料的解釋精度；候征等[15]將ABC算法用于一維瑞雷波多階模式聯(lián)合反演，研究結(jié)果表明該算法可提高反演精度；Wen等[16]將該算法應(yīng)用于某煤礦2.5維CSAMT數(shù)據(jù)反演，揭示出疑似陷落柱的地電結(jié)構(gòu)。需要指出的是，ABC算法雖然具有較好的探索能力，但局部搜索能力較弱，收斂相對較慢[17-18]。引領(lǐng)蜂在蜜源附近搜索時，只隨機(jī)地對解的某一個分量進(jìn)行更新，當(dāng)待解參數(shù)較多時，這樣搜索方式效率偏低。此外，優(yōu)質(zhì)蜜源對蜂群影響不足。

在非線性優(yōu)化算法研究領(lǐng)域，雙種群架構(gòu)是增強(qiáng)算法尋優(yōu)能力的重要思路。趙燕偉等[19]實(shí)現(xiàn)了基于雙種群的遺傳算法，指出改進(jìn)算法可以提高算法的全局收斂性能；吳亮紅等[20]提出一種雙種群差分進(jìn)化(DE)算法，不同的種群使用不同的變異策略，仿真實(shí)驗表明該算法全局搜索能力強(qiáng)、收斂速度快；暴勵等[21]采用雙種群架構(gòu)使ABC算法與DE算法并行進(jìn)化，每隔一定進(jìn)化次數(shù)就分別比較兩個種群的最優(yōu)解個體，用最優(yōu)個體替代次優(yōu)個體，以促進(jìn)種群間的相互學(xué)習(xí)；何光杰等[22]提出了一種雙種群的準(zhǔn)粒子群算法，兩個種群分別使用不同的粒子進(jìn)化方式，增強(qiáng)了算法的尋優(yōu)能力；陳亞峰等[23]借助雙種群對螢火蟲算法的種群進(jìn)行了全局和局部種群劃分，一定程度上解決了螢火蟲算法全局搜索和局部搜索不能兼顧、快速收斂和過早停滯相矛盾的問題。

基于前人的研究成果，本文提出一種基于雙種群架構(gòu)的ABC算法，并以大地電磁測深(MT)和重力聯(lián)合反演為例，驗證該算法的適用性和實(shí)用性。

1 基于雙種群架構(gòu)的ABC算法

1.1 方法原理

標(biāo)準(zhǔn)ABC算法的基本流程可以概括為以下四個方面：蜜源的生成、引領(lǐng)蜂尋找更優(yōu)質(zhì)蜜源、跟隨蜂決定是否跟隨及判斷是否達(dá)到局部搜索限制次數(shù)[24]。在標(biāo)準(zhǔn)ABC算法的基礎(chǔ)上，本文提出基于雙種群架構(gòu)的改進(jìn)ABC算法：將種群一分為二，在種群1中，嘗試引入交叉和變異策略以提高ABC算法的搜索效率；在種群2中，嘗試引入最優(yōu)解鄰域搜索以提高ABC算法的局部尋優(yōu)能力。兩個種群每經(jīng)歷一定代數(shù)后就進(jìn)行互相交流。

1.1.1 變異和交叉操作

在標(biāo)準(zhǔn)ABC算法生成蜜源階段，只對模型分量進(jìn)行隨機(jī)的一維搜索，即

(1)

在標(biāo)準(zhǔn)ABC算法搜索策略更新的蜜源U1a劣于上代蜜源時，按照下式進(jìn)行搜索以提高尋優(yōu)效率

(2)

式中：U1b是根據(jù)變異和交叉操作搜索出的新蜜源；F和CR分別表示變異因子和交叉因子，取值范圍均為[0,1];rand (·)表示取隨機(jī)數(shù)。

1.1.2 最優(yōu)解鄰域搜索

上節(jié)介紹的變異和交叉操作加快了對個體的開發(fā)，但是優(yōu)質(zhì)蜜源在蜂群中仍然無法傳播。故此，本文在種群2中加入最優(yōu)解鄰域搜索策略，利用下式對優(yōu)質(zhì)蜜源進(jìn)行充分開采

(3)

式中：U2是種群2更新的蜜源；Xbest表示種群2中的最優(yōu)個體，該參數(shù)的引入旨在提高ABC算法局部尋優(yōu)能力并加快收斂速度。具體而言，式(3)表示將上次迭代中的最優(yōu)個體作為被繼承的父代，在最優(yōu)解鄰域進(jìn)行探索。

兩個種群每隔一定的代數(shù)就分別隨機(jī)選取若干個個體并交換其在種群中的位置，進(jìn)而促進(jìn)優(yōu)良蜜源信息在種群間的傳播。該操作可在平衡全局尋優(yōu)和局部尋優(yōu)能力的同時維護(hù)種群的多樣性，進(jìn)一步提高解的搜索效率。

1.1.3 基本流程

基于雙種群架構(gòu)的ABC算法流程見圖1，具體如下。

圖1 改進(jìn)ABC算法流程圖

(1)生成初始蜜源：設(shè)置最大蜜源數(shù)、局部限制次數(shù)limit、最大迭代次數(shù)T。將蜜源平均分為兩部分，即種群1和種群2。

(2)種群1中的引領(lǐng)蜂搜尋蜜源：在種群1中，引領(lǐng)峰按照標(biāo)準(zhǔn)ABC算法搜索策略(式(1))搜索新蜜源，對比新蜜源適應(yīng)度f(U1a)與上一次迭代發(fā)現(xiàn)蜜源的適應(yīng)度f(Xi)并擇優(yōu)選擇，若發(fā)現(xiàn)新蜜源的適應(yīng)度優(yōu)于上次迭代發(fā)現(xiàn)蜜源，則保留新蜜源，直接進(jìn)入步驟(3)；反之，按照式(2)借鑒交叉和變異的思想搜索新蜜源，再對比新蜜源的適應(yīng)度f(U1b)與上一次迭代發(fā)現(xiàn)蜜源的適應(yīng)度，擇優(yōu)選擇新蜜源。

(3)種群2中的引領(lǐng)峰搜索蜜源：在種群2中，按照式(3)在最優(yōu)解附近搜索蜜源，同樣對比新蜜源適應(yīng)度f(U2)與上一次迭代發(fā)現(xiàn)蜜源的適應(yīng)度，擇優(yōu)選擇新蜜源。每迭代一定次數(shù)，兩個種群隨機(jī)交換種群中若干個蜜源的位置,生成新蜜源。

(4)跟隨蜂的局部尋優(yōu)：完成上述過程后，引領(lǐng)蜂飛回，與跟隨蜂交流，與標(biāo)準(zhǔn)ABC算法一樣，利用輪盤賭的方式計算跟隨概率，然后選擇是否更新。

(5)判斷蜜源Xi經(jīng)歷的局部搜索次數(shù)k是否大于限制次數(shù)limit：若是，則該引領(lǐng)蜂轉(zhuǎn)變?yōu)閭刹旆洌瑫r按照初始蜜源的生成方式產(chǎn)生新蜜源;若否，保留該蜜源。

(6)轉(zhuǎn)到步驟(2)進(jìn)行下一次循環(huán),直到滿足最大迭代次數(shù)。

1.2 算法測試

本文選擇三種典型的測試函數(shù)[25]檢驗改進(jìn)ABC算法的性能。

Griewank函數(shù)是典型的非線性多模態(tài)函數(shù)，其全局最小值0在(x1,x2,…,xn)=(0,0,…,0)處取得。此函數(shù)存在許多局部極小值點(diǎn)，并隨著問題維數(shù)的增加而增加。此函數(shù)在優(yōu)化過程中，具有廣闊的搜索空間，通常被認(rèn)為是優(yōu)化算法很難處理的復(fù)雜多模態(tài)問題。其函數(shù)形式為

-600≤xn≤600

(4)

式中N表示函數(shù)的維度。

Resenbrock函數(shù)是一個測試最優(yōu)化算法性能的非凸函數(shù)，全局最小值0在(x1,x2,…,xn)=(1,1,…,1)處取得，可用來檢測算法的精度。其函數(shù)形式為

-15≤xn≤15

(5)

Rastrigin函數(shù)存在多個波峰與波谷，易陷入局部極小，全局最小值0在(x1,x2,…xn)=(0,0,…,0)處取得，可用來檢測算法跳出局部極小的能力。其函數(shù)形式為

-15≤xn≤15

(6)

分別利用標(biāo)準(zhǔn)和改進(jìn)的ABC算法對上述三種函數(shù)進(jìn)行尋優(yōu)測試，迭代次數(shù)設(shè)為500，每種算法獨(dú)立重復(fù)試驗30次，試驗結(jié)果見表1和圖2。可看出，與標(biāo)準(zhǔn)ABC算法相比，改進(jìn)后的ABC算法不僅有較快的收斂速度，收斂精度也有較大提高，具有較強(qiáng)的全局尋優(yōu)能力。

表1 兩種算法尋優(yōu)結(jié)果對比

圖2 基于標(biāo)準(zhǔn)的和改進(jìn)的ABC算法的不同函數(shù)進(jìn)化曲線(a)Griewank函數(shù)；(b)Resenbrock函數(shù)；(c)Rastrigin函數(shù)

2 MT和重力數(shù)據(jù)聯(lián)合反演目標(biāo)函數(shù)

根據(jù)正則化反演理論[26]，MT數(shù)據(jù)和重力數(shù)據(jù)正則化聯(lián)合反演的目標(biāo)函數(shù)為

(7)

本文采取雙正則化因子自適應(yīng)調(diào)整方案[27-28]。MT數(shù)據(jù)正演采用有限差分法，重力數(shù)據(jù)正演采用網(wǎng)格累加法。聯(lián)合反演建模采用于鵬等[29]提出的物性參數(shù)隨機(jī)分布的共網(wǎng)格模型建模技術(shù)，該技術(shù)已被應(yīng)用于多種地球物理方法的二維聯(lián)合反演[5-6,30-33]。

3 模型測試

3.1 模型試驗一

為了驗證ABC算法在二維MT和重力數(shù)據(jù)聯(lián)合反演中的應(yīng)用效果，設(shè)計了圖3所示的電阻率和密度模型。模型的背景電阻率為1000Ω·m，密度為2.68g/cm3。模型中包含四個異常體。MT數(shù)據(jù)的模擬頻率范圍為0.001～320Hz，按對數(shù)等間距取38個頻點(diǎn)，測點(diǎn)間距為1km。

圖3 模型一電阻率(a)和密度(b)模型

設(shè)計了三種反演方案：方案一為基于標(biāo)準(zhǔn)ABC算法的聯(lián)合反演；方案二為基于雙種群架構(gòu)的ABC算法單獨(dú)反演；方案三為基于雙種群架構(gòu)的ABC算法聯(lián)合反演。這三種方案的蜜源數(shù)量最大值均設(shè)為300，最大迭代次數(shù)為50。種群中初始蜜源位置按照真實(shí)模型上下擾動一定范圍隨機(jī)生成：異常體的上下界面深度擾動范圍為真實(shí)值±30%，電阻率擾動范圍為真實(shí)值±30%，密度擾動范圍為真實(shí)值±0.02g/cm3。

初始模型及反演結(jié)果見圖4。圖4a、圖4b是隨機(jī)生成的電阻率和密度的初始擾動模型。首先，對比基于改進(jìn)ABC算法的聯(lián)合反演(方案三)和單獨(dú)反演(方案二)結(jié)果，可以看到，雖然重力異常曲線(圖4j)的擬合效果都較好，但是單獨(dú)反演結(jié)果(圖4f)的物性及界面還原效果較差。不同方案的反演結(jié)果與真實(shí)模型的殘差(圖5)也表明，聯(lián)合反演優(yōu)于單獨(dú)反演。

圖4 模型一及不同方案反演結(jié)果(a)初始電阻率擾動模型；(b)初始密度擾動模型；(c)方案一電阻率聯(lián)合反演結(jié)果；(d)方案一密度聯(lián)合反演結(jié)果；(e)方案二電阻率單獨(dú)反演結(jié)果；(f)方案二密度單獨(dú)反演結(jié)果；(g)方案三電阻率聯(lián)合反演結(jié)果；(h)方案三密度聯(lián)合反演結(jié)果；(i)反演目標(biāo)函數(shù)迭代誤差曲線；(j)重力異常擬合曲線

對比基于標(biāo)準(zhǔn)ABC算法的聯(lián)合反演(方案一)與基于改進(jìn)ABC算法的聯(lián)合反演(方案三)結(jié)果可見，在標(biāo)準(zhǔn)ABC算法的反演結(jié)果(圖4c、圖4d)中，界面起伏嚴(yán)重，物性差異大；而采用基于雙種群架構(gòu)的聯(lián)合反演(圖4g、圖4h)，物性界面起伏小且物性分布連續(xù)。由圖5的殘差圖也可以看出，標(biāo)準(zhǔn)ABC算法的反演結(jié)果偏離真實(shí)模型較大。此外，從迭代誤差曲線(圖4i)也可以看出，改進(jìn)ABC算法收斂快、尋優(yōu)精度高。

圖5 不同方案電阻率(左)和密度(右)反演結(jié)果與理論模型的殘差(a)方案一；(b)方案二；(c)方案三

3.2 模型試驗二

為了檢驗本文反演算法對復(fù)雜模型的適應(yīng)性，本文設(shè)計了圖6所示的電阻率界面與密度界面不完全一致的地球物理模型。該模型的背景電阻率為500Ω·m，背景密度為2.65g/cm3。反演過程中，假設(shè)電阻率和密度都是完全共界面進(jìn)行擾動，其余參數(shù)設(shè)置與前述模型試驗相同。初始擾動模型見圖6c和圖6d，基于改進(jìn)ABC算法的聯(lián)合反演結(jié)果見圖6e和圖6f。可以看出，反演的電阻率和密度異常分布與理論模型基本相符。從圖6i所示的迭代誤差曲線及圖6j所示的重力異常擬合曲線也可以看出，反演得到的異常體電阻率和密度分布基本可靠。

圖6 模型二及聯(lián)合反演結(jié)果(a)電阻率模型；(b)密度模型；(c)初始電阻率擾動模型；(d)初始密度擾動模型；(e)電阻率聯(lián)合反演結(jié)果；(f)密度聯(lián)合反演結(jié)果；(g)電阻率殘差；(h)密度殘差；(i)目標(biāo)函數(shù)迭代誤差曲線；(j)重力異常擬合曲線

4 實(shí)測數(shù)據(jù)應(yīng)用

為進(jìn)一步驗證改進(jìn)ABC算法的實(shí)用性，選擇中國下?lián)P子地區(qū)某測線MT和重磁實(shí)測數(shù)據(jù)進(jìn)行聯(lián)合反演。

依據(jù)測線實(shí)測數(shù)據(jù)的OCCAM反演結(jié)果、物性統(tǒng)計及文獻(xiàn)[33]所揭示的構(gòu)造格架建立初始模型，并對待解參數(shù)進(jìn)行隨機(jī)擾動(圖7)。反演過程對Tg印支面以上的界面不做反演，對其以下界面進(jìn)行開放反演。物性約束范圍根據(jù)該測區(qū)的巖石物性資料統(tǒng)計給定，見表2。

圖7 實(shí)際測線電阻率(a)和密度(b)初始擾動模型

表2 工區(qū)物性統(tǒng)計結(jié)果

對比電阻率和密度的初始擾動結(jié)果，可見電阻率和密度的聯(lián)合反演剖面(圖8)具有明顯的成層性。此外，在反演剖面橫向185～200km、縱向-5～-8km范圍內(nèi)，存在物性分布變化劇烈的地層，結(jié)合兩側(cè)的物性分布特征，推測此處可能存在斷裂。反演的密度和電阻率模型的正演響應(yīng)與實(shí)測數(shù)據(jù)的對比分別見圖9和圖10。上述認(rèn)識與文獻(xiàn)[33]基于模擬退火算法的聯(lián)合反演結(jié)果也基本相符，驗證了本文提出的雙種群ABC算法的實(shí)用性。

圖8 基于雙種群ABC算法的實(shí)測數(shù)據(jù)聯(lián)合反演的電阻率(a)和密度(b)剖面

圖9 重力異常擬合曲線

圖10 實(shí)測視電阻率剖面(a)及反演結(jié)果對應(yīng)的視電阻率剖面(b)

5 結(jié)束語

(1)本文通過雙種群架構(gòu)將蜂群一分為二。種群1結(jié)合標(biāo)準(zhǔn)ABC算法和變異交叉操作策略，保證了算法的全局尋優(yōu)能力，又提高了解的搜索效率；種群2則利用最優(yōu)解鄰域的搜索，增強(qiáng)了算法的局部尋優(yōu)能力。算法測試和模型試驗結(jié)果證明雙種群架構(gòu)可以提高ABC算法的尋優(yōu)能力。

(2)MT和重力數(shù)據(jù)聯(lián)合反演模型試驗及實(shí)際數(shù)據(jù)處理結(jié)果表明，MT與重力數(shù)據(jù)的聯(lián)合反演優(yōu)于單獨(dú)MT或重力數(shù)據(jù)反演，改進(jìn)ABC算法可以應(yīng)用于實(shí)際數(shù)據(jù)處理，具有一定的實(shí)用性。

需要指出的是，ABC算法具有很強(qiáng)的開發(fā)潛力，多種群的設(shè)置方案也是不唯一的，如何進(jìn)一步提升該算法的尋優(yōu)能力值得進(jìn)一步研究。