自適應權重蜻蜓算法及其在瑞雷波頻散曲線反演中的應用

高 旭 于 靜 李學良 胡天躍 何 川* 岳永強

(①北京大學地球與空間科學學院，北京 100871； ②中石化石油工程地球物理有限公司，北京 100020； ③中國科學院地質與地球物理研究所油氣資源研究室，北京 100029； ④保定北奧石油物探特種車輛制造有限公司，河北保定 072550)

0 引言

近地表地層模型的建立對地球物理研究和巖土工程具有重要意義[1]。瑞雷波的頻散曲線對橫波速度具有較強敏感性，因此常利用它反演地層的橫波速度[2]。與常規勘探方法相比，瑞雷波勘探技術具有分辨率高、應用范圍廣、檢測設備簡單且速度快等優點[3]，在淺層勘探中得到了廣泛應用。

對瑞雷波頻散曲線進行反演是瑞雷波勘探的關鍵環節[4-5]，反演算法可分為線性和非線性兩大類。由于反演本身具有非線性和多極值特點[6]，因此諸如最小二乘法[7]和Occam算法[8]等傳統線性算法在求解時易陷入局部極小值。遺傳算法[9-10]、模擬退火法[11-12]、BP神經網絡法[13]、粒子群算法[14-15]等非線性算法雖不需滿足觀測值與未知量之間是線性關系的假設條件[16-17]，且具有全局搜索的優勢，但在求解復雜問題時會出現參數不易調整等弊端，限制了這類方法的應用范圍[18]。

蜻蜓算法(Dragonfly algorithm，DA)是Mirjalili[19]于2016年提出的一種新型非線性算法，即以定義多個行為向量及其權重的方式模仿蜻蜓的飛行規律與協作模式，實現對搜索空間的“探索”與“開發”。然而在原算法的每步迭代過程中，所有蜻蜓個體的同類權重被賦予相同數值，這種無差別賦值方式可能在實際求解過程中降低算法尋找全局極小值的效率。為此，本文提出基于自適應權重的蜻蜓算法(Adaptive weight dragonfly algorithm，AWDA)，以蜻蜓個體的適應度為依據，賦予與其相匹配的自適應權重，及時調整蜻蜓個體的搜索策略，達到提高求解效率的目的。

通常，基階模式波的頻散曲線更易于拾取，是求取地層橫波速度的主要反演對象。當地下含有高速硬夾層或低速軟夾層時，高階模式波的頻散能量在高頻區域占據主導地位。研究[20-22]表明，高階頻散曲線蘊含豐富地質信息，若對基階、高階頻散曲線進行聯合反演可顯著提高結果的精度。

因此，本文首先利用測試函數，對AWDA算法、改進的蜻蜓算法[23](Improved dragonfly algorithm，IDA)、改進的粒子群算法[24](Improved particle swarm optimization algorithm，IEPSO)和改進的自適應遺傳算法[25](New improved adaptive genetic algorithm，NIAGA)分別進行對比測試，驗證了本文所提算法具有更強尋找全局極小值能力。隨后使用這一新型算法，通過理論及實際面波記錄，分別對基階和高階瑞雷波頻散曲線進行反演，新算法可顯著提高地層橫波速度反演結果的準確性和穩定性，取得了較為理想的應用效果。

1 基于自適應權重的蜻蜓算法

1.1 原理描述

蜻蜓以群組方式進行遷徙和覓食兩種運動，這兩種運動規律能分別與非線性算法中的“探索”、“開發”階段一一對應：遷徙中的蜻蜓群體通過遠距離飛行可遍布整個搜索空間，這是“探索”的主要目標； 處于覓食狀態下的蜻蜓群體通過短距離移動可對其周圍區域進行詳細搜尋，這是“開發”的主要特征。在蜻蜓算法中[19]，蜻蜓個體的運動軌跡通過避撞、結隊、聚集、覓食和避敵這五種行為的共同作用，并賦予其相應的權重，從而實現對搜索空間的“探索”和“開發”。

蜻蜓在飛行過程中受某種行為的影響程度由相應行為的權重決定，故對其進行合理的賦值尤為關鍵。針對原蜻蜓算法存在的同類權重無差別賦值問題，本文提出AWDA算法，即根據迭代過程中個體適應度之間的差異，自動調整原聚集、避撞和結隊等三種權重，其具體運算流程(圖1)如下。

圖1 AWDA算法流程

(1)初始化參數。假設蜻蜓的數量為n，搜索空間的維度為D，則蜻蜓的位置X及飛行速度V為

(1)

式中：i=1，2，…，n，表示蜻蜓編號；d為D中任一維度。

根據具體問題，設定聚集權重c、避撞權重s、結隊權重a、覓食權重f、避敵權重e、慣性權重w及鄰域半徑rc的初始值。

(2)計算適應度。將所有蜻蜓個體的位置代入目標函數，計算其適應度。

(3)更新蜻蜓個體的位置。當兩只蜻蜓的距離小于鄰域半徑時，就認為這對蜻蜓是相鄰的，與同一只蜻蜓相鄰的所有蜻蜓就構成了一個子群體。

若某蜻蜓周圍沒有鄰近個體時則將其命名為孤立個體，且通過Lévy飛行完成位置更新，其表達式為

(2)

Lévy飛行使個體進行大量的短距離局部搜索以及少量的長距離全局搜索，增強了個體飛行過程中的隨機性和“探索”性[26]。

若該蜻蜓存在鄰近個體，AWDA算法是通過各行為向量更新該蜻蜓的位置和飛行速度。因此，根據步驟(1)可定義避撞向量Si、結隊向量Ai、聚集向量Ci、覓食向量Fi和避敵向量Ei的表達式為

(3)

式中：Xj、Vj分別為第j個與當前個體相鄰的蜻蜓的位置和飛行速度；m是相鄰蜻蜓的數量；X+和X-分別表征食物和天敵的位置，將已記錄的適應度最大的個體作為食物，適應度最小的個體作為天敵。

隨后據下式更新蜻蜓的飛行速度和位置

(4)

(5)

式中：ci、si和ai分別是原算法第i個蜻蜓個體的聚集權重、避撞權重和結隊權重；b1和b2是常數，本文分別取為1.5和0.5，可根據具體問題對其進行相應的調整； 函數rankfit1(i)代表將同一群體內所有蜻蜓的適應度從低到高排序時，第i個蜻蜓在其中的排名； 函數rankfit2(i)代表將當前迭代中所有個體的適應度從低到高排序時，第i個蜻蜓個體在其中的排名。

在AWDA算法中，適應度大的個體會分配到較小聚集權重，抑制該個體向子群體中其他個體移動； 而適應度小的個體則會分配到較大聚集權重，促使其飛向適應度更大的個體。此外，適應度較大的蜻蜓個體將保持小間距、低速度的飛行狀態，對周圍做細致搜尋，“開發”是其主要任務；而那些適應度較小的個體，擴大間距、提高飛行速度才能增加獲取食物的可能性，“探索”對它們來說更為重要。如此，ADWA算法中的蜻蜓個體在迭代前期不但能對食物進行更高效的“開發”，中后期也會通過持續的“探索”以保存跳出局部極小值的能力。

(4)終止迭代。當最優蜻蜓個體的適應度滿足給定條件或達到最大迭代次數時，算法終止迭代并輸出結果，否則返回步驟(2)繼續迭代(圖1)。

1.2 算法測試

本文采用四個常用的測試函數[27]對AWDA、IDA、IEPSO及NIAGA等算法進行對比測試，以檢驗本文提出算法的有效性。測試函數的表達式及搜索區間如表1所示，空間維度D均為2。其中： Shpere和Schwefel 2.22為單極值函數，該類函數僅含有一個極小值，主要測試算法的收斂性和“開發”能力； Griewank和Ackley為多極值函數，它們包含了一個全局極小值和多個局部極小值，檢測算法的全局“探索”能力。四個測試函數的二維示意圖如圖2所示。所有算法的個體總數設定為50，最大迭代次數均為200。IDA、IEPSO和NIAGA三種算法的其他參數則分別與參考文獻[23-25]中的最佳參數組合設置保持一致。

圖2 四種測試函數的示意圖

表1 測試函數

以各算法對Griewank函數的一次運算結果為例，統計了所有個體在迭代過程中的位置分布，據該結果能更清晰直觀地分析算法對空間的“探索”與“開發”情形。IEPSO、NIAGA、IDA及AWDA四種算法求得的最小誤差分別為0.0038、0.0364、0.0071和0.0003。如圖3所示，IEPSO算法在迭代前期無法對空間進行充分“探索”，有相當多的區域僅分布零星個體甚至沒有個體，這可能會使算法最終無法尋找到全局極小值。NIAGA算法大概從第40次迭代開始，個體在若干局部極小值處呈現較強的聚集、重疊現象并一直持續到迭代結束，即表現出 “早熟收斂”特征。由于該算法的選擇、復制機制大幅度削弱了個體的多樣性，導致算法在中后期出現嚴重的個體趨同現象，這意味著計算機在大部分時間內做著無謂的重復運算。此外，由于變異的觸發概率較低，故通過變異跳出局部極小值從而繼續搜尋全局極小值是非常困難的。IDA算法和AWDA算法憑借避撞機制在迭代過程中能很好地抑制蜻蜓之間的重疊，因而在迭代前期能對整個區間進行廣泛的“探索”。在迭代后期，由于對食物周圍的“探索”不足，導致IDA算法最終收斂于局部極小值(x1=-3.13，x2=4.43)。但AWDA算法能較好地解決上述問題，在迭代后期仍能保持一定的“探索”能力，最終收斂于全局極小值。

圖3 四種算法對Griewank函數尋優的個體歷史位置分布圖

考慮到求解的穩定性因素，各算法對測試函數均進行50次獨立運算，并統計其誤差均值和標準差(表2)。

表2 各算法的尋優結果

對比可見，NIAGA算法求取的誤差均值在上述四種算法結果中是最大的，IEPSO算法的求解能力優于IDA算法。相較于其他三種算法，AWDA算法的誤差均值和標準差更小，顯示出該算法具有更出色的搜尋全局極小值的能力。

2 理論地質模型測試

頻散曲線的正演是瑞雷波勘探技術中的重要部分。地質模型一般由縱波速度vP、橫波速度vS、層厚H及密度ρ進行表征。將上述參數及頻率代入頻散方程中，求解該方程即可計算該頻率的瑞雷波相速度，以此獲得頻散曲線。反演的目的是尋找與已知頻散曲線匹配最好的地質模型。由于瑞雷波頻散曲線的反演問題是多極值的，為了降低反演難度，又因為橫波速度和層厚與瑞雷波頻散曲線的關系最為緊密，本文僅對這兩個參數進行計算，而影響較小的縱波速度和密度可根據先驗信息確定[28]。

2.1 基階頻散曲線反演結果

模型A為四層橫波速度遞增模型[18]，各地層參數和反演搜索范圍如表3所示。利用交錯網格有限差分方法[29]生成其理論面波記錄(圖4a)，采用相移法[30]獲得該地震記錄的瑞雷波頻散能量譜(圖4b)，可見頻散能量的極大值能準確對應由地層模型正演得到的理論頻散曲線(圖4b中黑點)，這說明使用相移法提取頻散曲線是可靠的。

圖4 模型A的理論面波炮集記錄(a)及相移法獲得的頻散能量譜(b)黑點是根據Knopoff方法計算的理論頻散點。圖7同

表3 模型A地質參數及反演搜索范圍

對個體做定量評估，即構建目標函數是頻散曲線反演的關鍵部分。目前最常用的是均方差函數，其基階頻散曲線反演的個體目標函數為

(7)

式中：cobs是實際提取的瑞雷波相速度；ccal是正演算出的瑞雷波相速度；M為頻散曲線的頻率點總數。

前文述及，IEPSO算法具有相對較強的“尋解”能力，因此本文分別采用該算法和AWDA算法反演頻散曲線，以對比驗證AWDA算法的反演效果。以上述兩種算法分別進行50次獨立計算，結果顯示該兩種算法反演得到的頻散曲線(圖5)與理論值均有較高擬合度，但通過AWDA算法所得地層橫波速度剖面(圖6b)的效果明顯優于IEPSO算法結果(圖6a)，各獨立運算結果也更集中地分布于理論值周圍。這從側面反映出瑞雷波頻散曲線的反演具有很強的多極值特性，而AWDA算法能有效地克服此不足，尋找到適應度更高的結果。從表4的反演結果統計可知，AWDA算法計算的地層參數的平均相對誤差為2.75%，顯著低于IEPSO算法結果的9.03%。

圖5 模型A的IEPSO(a)和AWDA(b)頻散曲線反演結果

表4 模型A反演結果統計

圖6 模型A的IEPSO(a)和AWDA(b)橫波速度反演結果

2.2 多階頻散曲線反演

層狀介質中瑞雷波存在多種傳播模式[31-32]，即在同一頻率下存在多個相速度值，速度從低到高依次為基階、一階高階、二階高階……等模式。在某些地質條件下，高階頻散能量在中、高頻段占據主導地位。若將高階頻散曲線納入反演中，就能增強對反演結果的約束，獲得更準確的橫波速度數據。

多階頻散曲線反演中個體的目標函數為

(8)

式中：L為頻散曲線總數；Mk為第k條頻散曲線的頻率點總數。

模型B為含軟夾層的五層模型，人工鋪設路面的地層結構即與此類似。地層參數和反演搜索范圍如表5所示； 所得理論面波記錄及生成的頻散能量譜如圖7所示，可見基階頻散能量主要分布于低頻區域，而高階頻散能量在50Hz以上更活躍(提取的基階、一階高階和二階高階頻散曲線如黑點連線)。

圖7 模型B的理論面波炮集記錄(a)及相移法獲得的頻散能量譜(b)

表5 模型B地質參數及反演搜索范圍

本文首先分別計算了IEPSO算法對基階、多階頻散曲線的反演結果(圖8)，可見IEPSO算法反演基階頻散曲線的結果與理論值的擬合度較好，但仍能觀察到20～60Hz頻段的反演結果與理論值存在一定程度偏差，雖呈現出反映低速夾層的“之”字形結構[33]，但并不如理論值那樣明顯。此外，多階頻散曲線反演結果中40～60Hz頻段的二階高階頻散曲線與理論值存在明顯偏差。

圖8 模型B的IEPSO算法的基階(a)和多階(b)頻散曲線反演結果

橫波速度反演結果如圖9所示，可見IEPSO算法的兩類反演結果與理論值的擬合度均較低。

圖9 模型B的IEPSO算法對基階(a)和多階(b)頻散曲線反演得到的橫波速度

隨后應用AWDA算法分別求得基階(圖10a)及多階(圖10b)頻散曲線反演結果，可見兩種情形的頻散曲線與理論值均有較高擬合度；從橫波速度反演結果(圖11)不難發現，多階頻散曲線的反演結果與理論值的擬合情況明顯更好。對上述反演結果的統計(表6)表明，相比僅反演基階頻散曲線的結果，反演多階頻散曲線結果中的第一層和第三層橫波速度的標準差分別由4.15和8.43降至0.91和2.24。

表6 模型A反演結果統計

圖10 模型B的AWDA算法的基階(a)和多階(b)頻散曲線反演結果

圖11 模型B的AWDA算法對基階(a)和多階(b)頻散曲線反演得到的橫波速度

以上結果說明AWDA算法的求解能力明顯高于IEPSO算法，且反演多階頻散曲線可進一步提高結果的準確度和穩定性。

3 實際測試

依據實際地震數據，進一步驗證AWDA算法反演瑞雷波頻散曲線的效果。數據來自歐洲研究機構發起的InterPacific項目[34-35]，測量地點位于意大利北部城鎮Mirandola。為增加勘探深度，同時使用主動源和被動源地震記錄提取頻散曲線。從檢波器陣列及鉆孔位置(圖12a)可見： 主動源陣列由48個4.5Hz的垂向檢波器以1m間隔線性擺放，震源采用8kg重錘敲擊地面方式，最小炮檢距為15m，采樣間隔為0.25ms，采集時長為2s； 被動源陣列采用兩個半徑分別為5m和15m的共中心圓環，采樣間隔為5ms，采集時長為60min。得到主動源(圖12b)和被動源(圖12c)地震記錄。

圖12 Mirandola地區實際面波勘探的檢波器陣列(a)及主動源(b)、被動源(c)地震記錄

圖13a是主動源地震記錄通過相移法生成的頻散能量譜，易見頻散能量可分為兩部分，即7～20Hz的基階頻散能量和22～45Hz的高階頻散能量。為準確利用該高階頻散能量，對實際數據做兩步法反演[36]：①對上述基階頻散曲線進行反演，得到橫波速度模型； ②對步驟①中反演得到的橫波速度模型做正演，計算其多階理論頻散曲線，隨后與上述高階頻散曲線進行比照，確定該高階頻散曲線的模態，再對基階和高階頻散曲線進行聯合反演，作為多階頻散曲線的反演結果。利用Geopsy軟件中的MSPAC模塊[37]提取被動源地震數據的頻散曲線，并得到相應的頻散能量譜(圖13b)。主、被動源方法提取的多階頻散曲線如圖13c所示，可見包含豐富的低頻信息，能顯著提高瑞雷波勘探深度。

圖13 主動源(a)、被動源(b)的頻散能量譜以及聯合方法提取的頻散曲線(c)

在實際面波勘探中，根據所提取的頻散曲線需確定最大勘探深度。通常最大勘探深度為基階頻散曲線最大波長的1/3～1/2[38]，本文選擇42m為此輪反演最大勘探深度。由于橫波速度和層厚與瑞雷波頻散曲線的關系最緊密[2]，故通過反演頻散曲線僅計算最大勘探深度以上地層的橫波速度和層厚，而將影響較小的縱波速度和密度設為合理的固定值。

根據橫波測井數據將勘探空間劃分為4層，各層橫波速度和層厚的搜索范圍、縱波速度及密度的設定如表7所示。若缺少橫波測井資料，可根據現場實際地質情況，結合提取的頻散曲線所包含的速度信息，設置一組搜索模型并分別進行反演，以獲得擬合度較高結果。之后，可對擬合度較高的反演結果的搜索模型進行更細致合理的調整并再做反演運算，以便得到更能反映地下真實情況的結果。

表7 反演搜索范圍及模型參數設置

本文使用IEPSO算法對提取的基階頻散曲線進行反演，并與AWDA算法的結果進行對比。算法的參數設置與理論地層模型測試一致。首先計算兩種算法對基階頻散曲線的反演結果，根據該反演結果正演計算多階頻散曲線，經過比照確定所提取的高階頻散曲線為二階高階頻散曲線。隨后利用AWDA算法對基階、二階高階頻散曲線進行聯合反演。對于頻散曲線，IEPSO算法計算的結果(圖14a)在10Hz以下與基準值的擬合度相較好，但在10Hz以上出現明顯偏差，最大誤差達2.12%。AWDA算法反演基階頻散曲線的結果(圖14b)與基準值的吻合情況優于前者，最大誤差不超過0.73%。AWDA算法反演多階頻散曲線的結果(圖14c)同樣能很好地擬合提取的頻散曲線，最大誤差僅為0.55%。

圖14 實際資料頻散曲線反演結果

從IEPSO算法反演所得橫波速度剖面(圖15a)可見，各獨立運算結果在淺層橫波速度及第三層厚度的確定與橫波測井數據存在一定程度偏差和波動。再看AWDA算法基階頻散曲線反演結果(圖15b)，淺層橫波速度的穩定性得到提升，且更接近于橫波測井資料。第三層厚度值雖仍存在一定波動，但與IEPSO算法結果相比更準確。從圖15c可見，利用AWDA算法對多階頻散曲線進行反演可進一步提高地層速度結構的穩定性，且對第三層厚度穩定性的提升尤為明顯。

圖15 不同方法針對實際資料反演的橫波速度

總之，實際資料測試結果再次驗證，相比本文提及的其他非線性算法，使用AWDA算法反演瑞雷波頻散曲線能顯著提高結果的準確性。

4 結論與討論

本文根據蜻蜓個體的適應度對行為權重進行調控，構建了AWDA算法。即通過各行為機制的相互配合，不斷調整各行為所占比重，從而靈活地在迭代過程中更改“搜索”策略。將該算法用于瑞雷波頻散曲線反演，與IEPSO算法相比，AWDA算法計算的頻散曲線與基準值的誤差更小，反演得到的地層速度結構與橫波測井數據更接近，證明了該算法的可行性及有效性。

此外，對多階瑞雷波頻散曲線的反演能獲得更準確、穩定的地層速度結構。

當然，AWDA算法也存在一些缺陷，如計算耗時相對較長，對自適應行為權重賦值范圍的設定多基于經驗。對這些問題的持續關注和深入研究無疑會大幅度提高AWDA算法的求解能力和計算效率。