瑞雷面波頻散曲線的粒子群蟻群混合優化反演

王一鳴 宋先海*② 張學強

(①中國地質大學(武漢)地球物理與空間信息學院，湖北武漢 430074； ②中國地質大學(武漢)湖北省地球內部多尺度成像重點實驗室，湖北武漢 430074)

0 引言

瑞雷面波勘探具有能耗低、效率高、經濟便捷、在淺地表分辨率高等優點，在近地表勘探、災害地質調查、工程無損檢測等方面獲得了廣泛的應用[1-3]。

非均勻半空間層狀介質的橫波速度和地層厚度對瑞雷面波的頻散曲線的結構和形態影響最大[4-6]，對瑞雷面波頻散曲線反演可獲得一維橫波速度曲線，能對地下地質結構進行分層。因此頻散曲線反演是瑞雷波勘探中最重要的步驟之一[7-8]。

瑞雷面波頻散曲線反演是典型的高度非線性、多參數和多極值的地球物理優化問題。目前反演方法主要分為局部線性優化算法和非線性全局優化算法。局部線性算法包括最小二乘法、最光滑模型反演(OCCAM)算法等。該類方法計算速度快，但其主要缺點是反演結果的可靠性嚴重依賴初始模型，只有當初始模型接近真實模型時才能獲得較好的反演結果，否則很容易陷入局部極值，得出錯誤的地質解釋； 同時該類方法還需要計算偏導數信息，雅克比矩陣的求取精度也將直接影響反演結果的可靠性。非線性全局優化算法有蟻群優化(Ant Colony Optimization，ACO)算法、遺傳算法、粒子群優化(Particle Swarm Optimization，PSO)算法和模擬退火算法等，該類算法對初始模型要求較低，在全局搜索最優解，避免了局部線性算法存在的問題[9-11]。

粒子群優化法是一種非線性算法，它通過簡單的策略不斷地更新位置和速度引導算法搜索到最優解。蔡偉等[12]將該算法應用于瑞雷波頻散曲線反演，結果表明PSO反演具有全局尋優能力強、快速、穩定且適用性強的特點，但是隨著搜索次數的增加和反演參數的增加，群體易陷入平衡狀態，粒子的進化過程隨著時間的推移將會停滯不前。蟻群優化法同樣是一種非線性算法，模擬了螞蟻覓食的行為，通過不斷釋放信息素縮短覓食路徑。王天琦等[13]將該算法應用于瑞雷波頻散曲線反演，結果表明ACO法在求解單峰函數時收斂速度較快、局部尋優能力強、解的精度高，但面對多峰函數卻經常發生收斂早熟，無法求出全局最優值。

本文綜合考慮了PSO算法和ACO算法各自的優、缺點，在Shelokar等[14]的研究基礎上提出了利用粒子群與蟻群相結合的優化(PSACO)算法進行瑞雷波頻散曲線反演。該反演策略充分發揮PSO算法較強的全局尋優能力，在可行解空間進行大范圍全局搜索，有效避免了算法陷入局部極小的可能性； 又充分利用ACO算法不斷更新PSO算法中粒子的位置，粒子的進化不會停滯，迭代后期具有較強的局部尋優能力，能提高解的精度，加快后期收斂速度，從而有效克服了PSO和ACO算法的各自缺點。理論模型的測試結果表明，與單獨的PSO算法和ACO算法相比，PSACO算法更優越、更有效； 實測資料反演結果驗證了PSACO算法的實用性。

1 基本原理及反演流程

1.1 PSO算法

在PSO反演中，種群的候選解稱為粒子，它與相鄰粒子經驗共享，同時共存和進化。每個粒子在問題的搜索空間中飛行時(每個粒子為指定空間的一個多維向量，代表一個可行性解)，通過使用自己的飛行經驗(早期飛行中找到的最佳位置的記憶)和鄰近粒子的飛行經驗修改其速度以找到更好的解決方案(最優解)。粒子的位置和速度更新可表示為

(1)

(2)

(3)

式中：M為最大迭代次數；ws、we分別為初始和終止慣性權重系數。

1.2 ACO算法

ACO算法模擬螞蟻的覓食行為。當螞蟻從食物源到巢穴的往返過程中，它們會在地上釋放信息素，并形成信息素軌跡； 螞蟻能感知信息素的氣味和濃度，當它們在選擇路徑時，會選擇信息素濃度強(較短)的路徑。該算法可以解決復雜的組合優化問題，例如旅行商(TSP)問題和二次分配問題。以經典的TSP問題為例，螞蟻算子的移動搜索策略遵循

(4)

τij(t+1)=(1-ρ)τij(t)+Δτij(t)

(5)

式中：ρ為信息素揮發系數，經驗取值一般約為0.9； Δτij(t)為本次循環路徑(i，j)上信息素的增量，表示每只螞蟻算子在本次循環留在路徑上信息素含量之和。

1.3 PSACO算法的原理

在一次總迭代過程中交替運行PSO算法和ACO算法實現PSACO算法。ACO算法是一種局部搜索算法，螞蟻利用信息素引導機制更新粒子在早期發現的位置。在PSACO中，一種簡單的信息素引導的ACO機制應用于局部最優值的搜索。ACO算法中生成的螞蟻個數等于PSO算法中粒子的個數。每個螞蟻根據下式圍繞全局最佳位置生成對應的粒子

(6)

1.4 混合優化算法瑞雷面波頻散曲線反演步驟

瑞雷面波頻散曲線反演實質是一個對目標函數求全局最優解的過程。衡量反演結果優劣最直接的因素是反演的目標函數大小。本文瑞雷面波頻散曲線反演的目標函數(適應度)定義為

(7)

混合算法中的粒子個數應與蟻群個數保持相同，一般解向量的維度越低，設置的種群數較少，反之則更多。針對不同維度的問題一般設置種群的數量為2L(2L-1)，其中L表示地層層數。通過該經驗公式可更好地平衡反演精度與反演用時。

PSACO算法實現流程如下：

(1)初始化。設置粒子個數I、最大迭代次數M等參數；

(2)拾取待擬合頻散曲線信息并確定搜索邊界，在搜索范圍內隨機產生I個粒子；

(3)計算每個解的目標函數，并令最小目標函數值對應的粒子為全局最優解，單個粒子為局部最優解；

(4)利用粒子群優化算法更新粒子并計算適應度，若優于局部最優解和全局最優解，則進行粒子替換；

(5)利用蟻群算法在全局最優解附近產生I個螞蟻(可能解)，并計算每個解的適應度；

(6)比較螞蟻與粒子的適應度，若螞蟻解更優，則用螞蟻替代粒子；

(7)若全局最優解的適應度大于精度要求或沒有達到最大迭代次數，則返回步驟(4)；

(8)輸出全局最優解。

2 理論模型試算

應用PSACO算法對瑞雷面波勘探中常遇到的地質結構進行反演，分析不同模型的反演結果，評估算法的應用效果。

首先，設置每個模型的層數、每層的縱橫波速度(VP、VS)、密度以及層厚度(H)； 然后，用快速矢量算法正演模擬瑞雷波頻散曲線(搜索步長為2，頻帶范圍為0～100Hz)； 最后，用PSACO算法對頻散數據進行反演。通過分析目標函數的大小、各層橫波速度和厚度隨搜索次數的變化過程等，探究PSACO算法的實用性和穩定性。

2.1 三層地質模型

為研究算法的適用性，建立三個典型的三層地質模型：速度遞增型、含低速軟夾層型、含高速硬夾層型[15]。三個模型都采用相同的算法參數：初始慣性權重系數ws=0.7，終止慣性權重系數we=0.4，加速系數c1=2、c2=2，初始標準差為σ1=1，控制標準差變化的參數d=0.9。

2.1.1 模型1

模型1 為三層速度遞增型，相關參數以及搜索范圍如表1所示。利用PSACO算法對模型1的頻散曲線進行反演，共計反演20次，每次反演迭代100次(每個粒子都進行100次位置更新)，去除兩次與平均值相差最大的反演結果后，剩余結果取平均值(圖1)。

表1 模型1參數及搜索范圍

圖1a為目標函數隨迭代次數變化曲線，可以看出，前15次目標函數下降較快，到第70次時幾乎下降為0。由圖 1b～圖 1f可知： 當搜索到100次時，反演的每個參數與真值基本一致； 收斂曲線在60次后僅有極小幅度的波動。

由圖2可知，反演的頻散曲線和理論頻散線擬合很好，反演的橫波速度和層厚度與真實模型高度吻合。表 2為反演20次的結果統計。

由表2可知，模型1經過多次反演，PSACO算法所得結果的相對誤差都接近于0，標準差同樣很小，驗證了本文算法精度高、穩定性好。

圖1 模型1平均反演結果隨迭代次數的變化曲線(a)目標函數； (b)VS1； (c)VS2； (d)VS3； (e)H1； (f)H2

圖2 模型1的反演結果(a)反演與觀測頻散曲線對比； (b)橫波速度剖面

表2 模型1多次反演結果統計

2.1.2 模型2

用PSACO算法對三層含低速軟夾層地質模型模擬的頻散曲線進行反演，其模型的相關參數以及搜索范圍如表3所示。共計反演20次，每次反演迭代100次，去除兩次與平均值相差最大的反演結果，剩余結果取平均如圖3所示。

表3 模型2參數及搜索范圍

圖3 模型2平均反演結果隨迭代次數的變化曲線(a)目標函數； (b)VS1； (c)VS2； (d)VS3； (e)H1； (f)H2

圖3a為目標函數隨迭代次數的變化曲線，可見：前15次目標函數下降較快，到第80次時幾乎下降到零； 與模型1相比，搜索難度有所增加。由圖3b～圖3f可以看出，VS1、VS3收斂結果與真值非常接近，但VS2、H1、H2收斂精度有所下降，穩定性也有所降低，但仍與真值接近。

模型2的最終反演結果如圖4所示，反演的橫波速度剖面和理論模型幾乎重合。PSACO算法20次反演的統計結果如表4所示。

由表4可知，經過多次反演，模型2的PSACO算法各參數反演的相對誤差以及標準差均接近于0，說明PSACO算法對三層含低速軟夾層模型的反演結果精度高、穩定性好。

圖4 模型2的反演結果(a)反演與測量頻散曲線； (b)橫波速度剖面

表4 模型2多次反演結果統計

2.1.3 模型3

模型3為三層含高速硬夾層型，相關參數和參數的搜素范圍如表5所示。利用PSACO算法對頻散曲線反演共計20次，每次反演迭代100次，去除兩次與平均值相差最大的反演結果后，剩余結果取平均，如圖5所示。

表5 模型3參數及搜索范圍

由圖5a目標函數隨搜索次數變化曲線可以看出，前15次目標函數下降較快，到50次左右時趨近于0。由圖5b～圖5f可見，每個參數反演都能在搜索100次時與真值重合，收斂曲線在后期未發生明顯波動。

模型3反演的最終結果如圖6所示，反演的橫波速度剖面和理論模型幾乎重合。統計PSACO算法模型3反演20次的結果，如表6所示。

由表6可知，經過多次反演，PSACO算法的各參數反演均值與理論值十分接近，相對誤差較低，標準差控制在合理范圍內。

圖5 模型3平均反演結果隨迭代次數的變化曲線(a)目標函數； (b)VS1； (c)VS2； (d)VS3； (e)H1； (f)H2

圖6 模型3的反演結果(a)反演與測量頻散曲線； (b)橫波速度剖面

表6 模型3多次反演結果統計

三個典型模型數據反演結果均表明，PSACO算法精度高、穩定性好。

2.2 PSACO與單獨的PSO、ACO算法對比

四層含低速軟夾層地質模型(表7)參數較多，適合用來對比不同算法的穩定性、準確性。為有效對比，PSACO、ACO、PSO算法采用的反演參數保持一致：初始慣性權重系數ws=0.7，終止慣性權重系數we=0.4； 加速系數c1=2、c2=2； 初始標準差為1，控制σ變化的參數d=0.9； 影響因子α=0.4，β=0.4； 揮發系數ρ=0.9。

表7 四層低速軟夾層型模型參數及搜索范圍

圖7a為三種方法目標函數隨迭代次數的變化曲線對比，可見： PSACO算法能在目標函數快速下降的同時也能保證算法不陷入局部最優； PSO算法雖然在整個迭代過程中保持下降的趨勢，但是后期的搜索能力不強，解的精度不高； ACO算法在前期收斂速度快，但很明顯在后期出現收斂停滯，解的誤差偏大。由圖7b～圖7h可知，與PSO和ACO算法相比，PSACO算法搜索到100次時的結果更接近真值，在搜索的過程中穩定性更高，表明PSACO算法具有更高的穩定性和更強適用性。

圖8為四層模型三種方法反演結果的對比，正演頻散曲線和反演剖面基本相同，都接近真值。表8為四層模型三種方法反演結果的統計，可見，PSACO算法的7個反演參數的相對誤差和標準差都低于PSO和ACO算法。相對誤差低，說明該算法的全局尋優能力強； 標準差小，說明該算法的穩定性更高。

圖7 四層模型三種方法平均反演結果隨迭代次數的變化曲線對比(a)目標函數； (b)VS1； (c)VS2； (d)VS3；(e)VS4； (f)H1； (g)H2； (h)H3

圖8 四層模型三種方法的反演結果對比(a)反演與測量頻散曲線； (b)橫波速度剖面

表8 四層模型不同優化算法反演結果的統計

3 實例分析

通過對意大利北部某垃圾填埋場的實例反演試算，檢驗PSACO算法對實測數據的反演能力，驗證該算法的實用性。

3.1 實例背景介紹與參數取值

實例數據來自意大利東北邊的垃圾填埋場[16]，該區域的基巖層是灰巖，上層由未膠結的堆積物構成，約18m厚。由該地點的鉆孔資料可知每層的厚度及其分層情況。

同前文對理論模型反演的方法類似，泊松比和密度視為已知(通過已掌握的實測資料進行估算[15])，橫波速度和厚度為待反演參數，初始慣性權重系數ws=0.7，終止慣性權重系數we=0.4； 加速系數c1=2、c2=2； 初始標準差為1，控制σ變化的參數d=0.9。模型參數及參數搜索范圍如表9所示。

表9 實際數據的模型參數和反演搜索范圍

3.2 反演效果分析

由PSACO、ACO、PSO反演算法目標函數隨迭代次數變化曲線(圖9a)可見：PSACO算法能在目標函數快速下降的同時也能保證算法不陷入局部最優； PSO算法在后期的搜索能力不強，解的精度不高； ACO算法在前期收斂速度快，但很在后期出現收斂停滯，解的誤差偏大。由于地層厚度和介質速度是未知的，無法直觀對比三種算法反演參數VS、H與實測值的近似程度，由圖9b～圖9j曲線的穩定性可知，PSACO算法穩定性最高。

圖9 實際數據三種方法平均反演結果隨迭代次數的變化曲線對比(a)目標函數； (b)VS1； (c)VS2； (d)VS3； (e)VS4；(f)VS5； (g)H1； (h)H2； (i)H3； (j)H4

由頻散曲線對比(圖10a)可知，反演的頻散曲線與實測頻散點擬合度高。通過將不同算法反演的橫波速度剖面圖(圖10b)和鉆孔資料(圖10c)對比可以發現，PSACO算法與鉆孔資料的匹配度更高，反演的橫波速度剖面更接近實際地層。

圖10 實測資料三種算法反演結果及與鉆孔剖面對比(a)反演頻散曲線與實測頻散曲線對比； (b)反演的橫波速度剖面； (c)地層柱狀圖(鉆孔資料)

4 結束語

本文通過理論模型試算、不同算法之間的對比、實測資料分析驗證了PSACO算法應用于瑞雷面波頻散曲線非線性反演的可行性、適用性、穩定性。PSACO算法具有優秀的全局尋優能力，在迭代后期局部尋優能力也較強； 該算法的成功率、精度、穩定性在總體上都優于PSO算法和ACO算法； PSACO算法應用至實測資料的瑞雷面波頻散曲線反演結果較其他算法更接近真值。