雙種群人工蜂群算法及其在MT和重力聯合反演中的應用

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

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

0 引言

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

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

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

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

1 基于雙種群架構的ABC算法

1.1 方法原理

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

1.1.1 變異和交叉操作

在標準ABC算法生成蜜源階段，只對模型分量進行隨機的一維搜索，即

(1)

在標準ABC算法搜索策略更新的蜜源U1a劣于上代蜜源時，按照下式進行搜索以提高尋優效率

(2)

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

1.1.2 最優解鄰域搜索

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

(3)

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

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

1.1.3 基本流程

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

圖1 改進ABC算法流程圖

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

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

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

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

(5)判斷蜜源Xi經歷的局部搜索次數k是否大于限制次數limit：若是，則該引領蜂轉變為偵察蜂，同時按照初始蜜源的生成方式產生新蜜源;若否，保留該蜜源。

(6)轉到步驟(2)進行下一次循環,直到滿足最大迭代次數。

1.2 算法測試

本文選擇三種典型的測試函數[25]檢驗改進ABC算法的性能。

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

-600≤xn≤600

(4)

式中N表示函數的維度。

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

-15≤xn≤15

(5)

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

-15≤xn≤15

(6)

分別利用標準和改進的ABC算法對上述三種函數進行尋優測試，迭代次數設為500，每種算法獨立重復試驗30次，試驗結果見表1和圖2。可看出，與標準ABC算法相比，改進后的ABC算法不僅有較快的收斂速度，收斂精度也有較大提高，具有較強的全局尋優能力。

表1 兩種算法尋優結果對比

圖2 基于標準的和改進的ABC算法的不同函數進化曲線(a)Griewank函數；(b)Resenbrock函數；(c)Rastrigin函數

2 MT和重力數據聯合反演目標函數

根據正則化反演理論[26]，MT數據和重力數據正則化聯合反演的目標函數為

(7)

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

3 模型測試

3.1 模型試驗一

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

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

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

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

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

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

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

3.2 模型試驗二

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

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

4 實測數據應用

為進一步驗證改進ABC算法的實用性，選擇中國下揚子地區某測線MT和重磁實測數據進行聯合反演。

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

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

表2 工區物性統計結果

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

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

圖9 重力異常擬合曲線

圖10 實測視電阻率剖面(a)及反演結果對應的視電阻率剖面(b)

5 結束語

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

(2)MT和重力數據聯合反演模型試驗及實際數據處理結果表明，MT與重力數據的聯合反演優于單獨MT或重力數據反演，改進ABC算法可以應用于實際數據處理，具有一定的實用性。

需要指出的是，ABC算法具有很強的開發潛力，多種群的設置方案也是不唯一的，如何進一步提升該算法的尋優能力值得進一步研究。