利用自適應混沌遺傳粒子群算法反演瑞雷面波頻散曲線

楊 博 熊章強 張大洲* 楊振濤

(①中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室，湖南長沙 410083； ②有色資源與地質災害探查湖南省重點實驗室，湖南長沙 410083； ③中南大學地球科學與信息物理學院，湖南長沙 410083； ④南方科技大學地球與空間科學系，廣東深圳 518055)

0 引言

瑞雷面波在速度垂向不均勻介質中呈頻散特性[1]，因此近年來在地震工程、工程地質勘查和無損檢測等領域得到了廣泛的應用[2-3]。瑞雷面波勘探的流程主要包括數據采集、頻散曲線提取[4-5]、頻散曲線反演及解釋等四個環節。頻散曲線的反演方法主要有局部線性化和非線性全局優化[6]兩大類。局部線性化反演方法中最具代表性的是最小二乘法[7]，該方法對初始模型依賴較大，因此要求數據處理人員具有較豐富經驗。常用的非線性全局優化方法包括遺傳算法[8]、BP神經網絡法[9]、模擬退火法[10]、粒子群算法(Particle Swarm Optimization，PSO)[11]和模式識別法[12]等，這類方法對初始模型的依賴較小，求解精度高，但針對復雜地層數據反演時易出現收斂速度慢、陷入局部極小值等問題，一定程度上影響了其實際應用效果。

粒子群算法作為一種全局優化算法，因收斂速度快、求解精度高等特點被廣泛應用于瑞雷面波頻散曲線的反演，但傳統的粒子群算法存在易陷入局部最優、過于早熟等問題[13]，因此對其進行了諸多改進。如慣性權重改進[14]、學習因子改進、粒子位置和速度改變[15]、加入混沌局部搜索[16]和混合遺傳算法[17]等。

本文基于上述研究成果，設計了一種能同步提高全局和局部搜索能力的自適應混沌遺傳粒子群算法(Adaptive Chaos Genetic Particle Swarm Optimization Algorithm，ACGPSO)，該算法采用自適應慣性權重和混合遺傳算法的交叉和變異操作以增強基本粒子群算法的全局搜索能力，約束全部粒子的行進速度，對最優粒子進行單維全分量的混沌局部搜索以增強其局部搜索能力。

瑞雷面波頻散包括基階和高階頻散曲線，目前瑞雷面波的反演主要還是采用基階模式的頻散曲線，這使得反演多解性問題非常突出。若在反演過程中加入高階模式頻散曲線，就可顯著降低反演結果的多解性。另一方面，當地層中含有低速軟夾層或高速硬夾層時，高階模式的面波能量相對于基階模式的面波在中—高頻段漸漸占據主導優勢[18-19]，對于這類地層采用瑞雷波勘探時，必須將基階與多階瑞雷波頻散曲線進行聯合反演[20-21]才能獲取準確的地下橫波速度結構[22]。因此，在瑞雷面波頻散曲線反演算法設計時須考慮進行基階與高階頻散曲線聯合反演，以進一步提高反演精度。

為了檢驗ACGPSO算法的準確性，利用兩種測試函數對該算法進行測試，并將試算結果與遺傳算法(GA)和PSO的結果進行對比。在此基礎上對理論模型的無噪和含噪數據的基階和多階頻散曲線進行反演，驗證其有效性、穩定性和抗噪性； 最后將該方法應用于實際數據反演，進一步驗證其適用性。

1 粒子群算法的基本原理及其改進

1.1 粒子群算法的基本原理

粒子群算法也稱鳥群覓食算法，是由Kennedy等[23]開發的一種進化算法。在用戶給定的d維搜索范圍內，所有粒子通過跟蹤其本身找到的最優解(pi)和整個種群當前找到的最優解(g)更新自己，完成迭代尋優。在對比適應度值尋找到這兩個解后，通過下式更新自己的速度(vi)和位置(xi)

(1)

1.2 粒子群算法的改進

本文對粒子群算法的改進主要體現在四個方面： ①采用自適應慣性權重； ②設置粒子節速度； ③引入GA算法的交叉和變異操作； ④對與歷史最優解連續3次重復相同的g進行單維全分量的混沌局部搜索。采用此四項改進即可獲得一種能同步提高全局和局部搜索能力的自適應混沌遺傳粒子群算法(ACGPSO)。

1.2.1 自適應慣性權重

在PSO算法的可調參數中，慣性權重w是最重要的參數，較大的w有利于提高算法的全局搜索能力，而較小的w會增強算法的局部搜索能力?；镜腜SO算法采用的是固定權重法，即采用一個固定的w。在面對復雜問題時，由于難以預知一個適宜的w使得算法收斂到正確的解，因此固定權重法在處理實際問題時通常并不可取。為了平衡PSO算法的全局搜索能力和局部搜索能力，本文采用非線性的動態慣性權重系數，其表達式為[24]

(2)

式中：wmax、wmin分別表示w的最大和最小值，本文取0.9和0.6；f為粒子當前目標函數值；favg、fmin分別為當前所有粒子的平均目標函數值和最小目標函數值。分析式(2)可知： 當某粒子的目標函數值小于平均目標函數值時，對應慣性權重較小，從而保護了該粒子； 反之，當某粒子目標函數值大于平均目標函數值時，對應的慣性權重較大，使該粒子向較好搜索區域靠攏。因此，該自適應慣性權重能根據當前粒子整體狀況動態地給出不同粒子的慣性權重，全局性地改善算法的搜索能力。

1.2.2 設置節速度

粒子在解空間中搜索的速度須得到控制，太大導致局部搜索能力不夠，太小則收斂過慢。為使粒子具有合理的搜索速度，粒子的節速度設定為

(3)

1.2.3 交叉和變異操作

PSO算法收斂速度快、過于早熟的缺陷緣于進化后期基因多樣性的缺失。針對此現象本文引入遺傳算法中的交叉和變異操作，使其能較好地彌補該缺陷。若某次迭代滿足交叉概率Pc(本文取0.9)，則選取指定數量的粒子放入雜交池內，池中粒子隨機兩兩雜交，產生同樣數目的子代，并用子代粒子替代父代粒子。針對交叉得到的變異概率小于預設閾值(本文取0.05)的子代，進行特定基因位置的突變。子代位置(xchild)和速度(vchild)分別為

(4)

(5)

1.2.4 單維全分量的混沌局部搜索

PSO算法本身具有記憶性，能隨時將搜尋到的最優解記錄在g中。如果當前種群最優解g已與歷史最優解連續3次重復相同，就對第3次的g執行單維全分量混沌局部搜索，即從解向量的第1分量到第d分量，分別只對其中一個分量增加一個混沌擾動量，從而形成d個混沌解。計算這d個混沌解的目標函數值，并將其與g的目標函數值比較，選出最小目標函數值對應的解替換g，并將其作為本次迭代搜尋到的最優解?；煦缃獾木唧w公式為

(6)

1.3 ACGPSO算法測試

為了驗證ACGPSO算法處理復雜問題的性能，選擇Rastrigin函數和Griewank函數對其進行測試。Rastrigin函數由于有非常多的局部最小值點和局部最大值點，很容易使算法陷入局部最優，而得不到全局最優解，因此對模擬退火等進化算法具有很強的欺騙性。Griewank函數是典型的非線性多模態函數，具有廣泛的搜索空間，通常被認為是優化算法很難處理的復雜多模態問題。利用上述兩種函數對ACGPSO算法進行測試可驗證算法的準確性。兩者的二維表達式分別為[22]

R(x,y)=20+x2+y2-10[cos(2πx)+

cos(2πy)]x,y∈[-5.12,5.12]

x,y∈[-100,100]

兩函數均在(x,y)=(0,0)處取得全局最小值0。圖1為兩函數的曲面，從中可直觀地看到它們有很多局部最優解。

分別運用GA、PSO和ACGPSO算法對Rastri-gin函數和Griewank函數尋優求解，設定迭代次數為 1000，種群數為40。三種算法均求取10次平均值作為該算法的最終值，最終計算結果(全文均保留小數點后兩位)如表1和圖2所示。

由表1可知： GA和PSO兩種算法在當前設置的搜索范圍內未搜索到理論解，尋優能力過弱； 而ACGPSO算法能100%收斂到測試函數的理論解，尋優結果明顯優于另兩種算法。就求解速度而言，ACGPSO算法略慢于PSO算法的主要原因是前者增加了交叉和變異操作及混沌局部搜索，但其速度快于GA算法。從圖2可見： GA和PSO兩種算法迭代曲線下降緩慢，收斂不到理論解附近； 而ACGPSO算法在1～100代目標函數值持續下降，迅速收斂，基本不存在GA和PSO算法出現的多代停滯進化的現象，100代后尋找到理論解。

通過上述對比，證明綜合改進后的ACGPSO算法克服了PSO算法疏于開發、易陷入局部最優解的缺陷，增強了其全局和局部搜索能力，提升了對全局最優解的收斂速度，且顯著提高了最終解的精度。因此，在面對多參數、多極值的復雜問題時，ACG-PSO算法比常規PSO或GA算法具有更強尋優能力和更快收斂速度，這為瑞雷面波頻散曲線的反演提供了很好的求解方案。

圖1 Rastrigin函數(a)和Griewank函數(b)的曲面表1 三種算法尋優結果對比

優化算法Rastrigin函數Griewank函數最優解最劣解平均解10次求解用時/s最優解最劣解平均解10次求解用時/sGA5.17×10-21.034.71×10-11.901.06×10-21.43×10-14.48×10-22.14PSO6.47×10-47.46×10-33.06×10-31.212.08×10-32.41×10-21.21×10-21.30ACGPSO0(100%)001.300(100%)001.34理論解000000

圖2 三種算法對Rastrigin函數的尋優迭代對比

2 理論模型試算

在完成Rastrigin和Griewank函數測試后，將ACGPSO算法應用于瑞雷面波理論模型頻散曲線反演。采用實際工程勘查中常見的層狀模型作為測試理論模型。模型A和B分別為四層含高速硬夾層、低速軟夾層的地質模型，其參數及反演搜索范圍如表2所示。由于瑞雷面波頻散曲線對于層厚和橫波速度較敏感[7]，因此選取這兩個參數作為反演目標。為了減少初始模型對反演的影響，利用智能優化算法反演時通常設置較大的橫波速度和層厚搜索范圍。本次測試設定的搜索上、下限與理論參數值相差60%。為了充分驗證ACGPSO算法的有效性和穩定性，反演均取10次的均值作為最終結果。反演迭代次數為80，粒子數是40。

2.1 基階頻散曲線反演測試

在瑞雷面波勘探中，通常從能量相對較強的基階頻散曲線中提取較準確的頻散曲線，即僅利用基階頻散曲線作為反演對象。模型A提取5～100Hz頻段基階頻散曲線，模型B提取5～70Hz頻段基階頻散曲線進行反演。

反演采用如下目標函數

(7)

式中：M為某一反演模型；N為實測頻散點個數；Vobs、Vcal分別為實測的與反演模型M的基階頻散曲線。采用該目標函數的ACGPSO算法反演結果如圖3和表3所示。

表2 模型參數及反演搜索范圍

圖3 模型A、B不含噪聲數據的反演結果(a)模型A頻散曲線； (b)模型A反演的速度剖面； (c)模型B頻散曲線； (d)模型B反演的速度剖面

從模型A和模型B計算得到的無噪聲頻散曲線的反演結果(圖3)可見：在預設的較大搜索范圍內(圖3b和圖3d中短劃線)，ACGPSO算法能對理論模型的實測數據(圖3a和圖3c中圓點)進行有效反演，反演模型的頻散曲線(圖3a和圖3c中虛線)與實測數據幾乎完全擬合； 各理論模型的參數(圖3b和圖3d中實線)也被很好地重建(圖3b和圖3d中虛線)，平均相對誤差分別僅為1.03%和0.70%，平均標準差分別為2.21和2.93。 這些充分表明，ACGPSO算法能對理論模型的頻散曲線進行有效而準確的反演。

實際應用中瑞雷面波所提取的頻散曲線中不可避免地含有一定的噪聲，因此反演方法的抗噪能力是非常重要的。為了檢驗ACGPSO算法對含噪聲頻散曲線反演的能力，對模型A和模型B的頻散曲線(圖3a和圖3c)加入10%的白噪聲，反演結果見圖4和表3所示。

從圖4可看出：反演所得模型的頻散曲線(圖4a和圖4c中虛線)能較好地擬合加噪后的頻散曲線(圖4a和圖4c中實線)； 同時反演所得到的模型參數(圖4b和圖4d中虛線)與理論模型參數(圖4b和圖4d中實線)非常接近，平均相對誤差分別為3.25%和2.23%，平均標準差分別為4.01和3.85。雖然加入噪聲后的平均相對誤差和平均標準差均有所增大，但反演結果依然具有可靠性和穩定性。因此，ACGPSO算法具有一定的抗噪性。以模型A為例，無論是對無噪還是對含噪數據反演，ACGPSO算法均能快速地在20代內收斂到最優解附近，隨后的迭代曲線下降緩慢，直至趨于平緩，表明目標函數已收斂到最優解(圖5)。這說明ACGPSO算法反演頻散曲線具有很好的收斂性。

2.2 多階頻散曲線聯合反演

對于類似模型A和模型B的含異常夾層的地層， 在利用瑞雷面波法勘探時所獲頻散曲線常呈現呈現“之”字形回折現象[3，25]，該頻散曲線隨著頻率增大，高階模式的面波能量相對于基階模式的面波漸漸占據主導優勢。另外，通過對一些模型的測試可知，不同階次的頻散曲線對同一層的敏感度會有所不同，同階次的頻散曲線在不同層和不同頻率上的敏感度也有所差異，這種情況下僅反演基階頻散曲線很難獲取所有層的準確信息。因此在反演過程中加入高階模式的頻散曲線進行聯合反演是非常有必要的。ACGPSO算法進行多階頻散曲線聯合反演時的目標函數為[21]

表3 不含噪與含噪數據反演結果統計表

圖4 模型A、B含噪聲數據的反演結果(a)模型A頻散曲線； (b)模型A反演的速度剖面； (c)模型B頻散曲線； (d)模型B反演的速度剖面

圖5 不含噪和含噪時模型A的目標函數變化曲線

(8)

多階頻散曲線聯合反演所得頻散曲線(圖6a和圖6c中虛線)能與各階實測頻散曲線(圖6a和圖6c中點線)很好地擬合； 反演模型參數(圖6b和圖6d中虛線)與理論模型參數(圖6b和圖6d中實線)的平均相對誤差分別僅為0.22%和0.92%。 平均標準差分別為1.42和2.79； 與僅用基階反演的結果相比，平均相對誤差分別下降了3.03%和1.31%，平均標準差分別下降了2.59和1.06。特別是模型A中反演最差的高速層的相對誤差和標準差分別由6.13%和15.16降至0.09%和1.62； 模型B中反演最差的第三層的相對誤差和標準差分別由2.92%和16.74降至2.43%和11.17。由此可見，運用多階頻散曲線聯合反演能有效降低解的多解性，提高解的穩定性和精度； 同時也證明ACGPSO算法不僅能應用于基階頻散曲線反演，也能應用于多階頻散曲線聯合反演。

圖6 模型A和模型B多階頻散曲線聯合反演結果(a)模型A頻散曲線； (b)模型A反演的速度剖面； (c)模型B頻散曲線； (d)模型B反演的速度剖面表4 多階頻散曲線聯合反演結果統計表

反演參數真實值模型A模型B反演均值相對誤差/%標準差反演均值相對誤差/%標準差h1/m21.990.580.012.000.180.12h2/m22.000.170.182.010.720.09h3/m22.000.060.251.990.250.14VS1/(m·s-1)150149.570.290.26198.480.763.12VS2/(m·s-1)250250.220.091.62152.901.943.60VS3/(m·s-1)200199.670.176.12156.072.4311.17VS4/(m·s-1)400400.800.201.49400.680.171.29

2.3 與PSO算法反演的對比

為了證明ACGPSO算法在反演應用中優于PSO算法，現用PSO算法對模型A含噪數據進行多階頻散曲線的聯合反演，反演的粒子數、迭代數和反演次數均與ACGPSO算法相同。將PSO算法與ACGPSO算法的最佳反演結果進行對比(圖7)。從圖7a可見：PSO算法收斂速度低于ACGPSO算法，最終的目標函數值也大于后者。從圖7b可見：ACGPSO算法反演得到的當前最佳模型在迭代前期能迅速向理論模型靠近，在后期能保持相對的穩定； 而PSO算法反演得到的當前最佳模型不僅向理論模型靠近緩慢，還存在反復跳躍的不穩定性現象，因此所得最終反演結果比ACGPSO算法差。這充分說明在頻散曲線反演中，ACGPSO算法比PSO算法具有更快的收斂速度、更強的穩定性和更高的反演精度。

圖7 ACGPSO算法與PSO算法對比(a)目標函數變化曲線； (b)當前最優模型與理論模型的相對誤差變化曲線

3 實測數據反演

對理論模型的基階和多階頻散曲線的反演，驗證了ACGPSO算法在頻散曲線反演中的可行性。為了進一步驗證ACGPSO算法的實用性，現對美國懷俄明M地區現場實測地震數據進行反演和分析。本次現場數據采集采用了48個8Hz垂直分量檢波器，道間距為0.9m，最小炮檢距為0.9m，選用錘擊震源。采集所得地震記錄如圖8a所示； 圖8b為從炮集記錄提取的頻散能量圖，通過采集能量最大點提取所需頻散點。從圖8b中可見頻散點分為兩段：第1段為10～27Hz的基階頻散曲線； 第2段為35～55Hz的高階模式頻散曲線。為了充分利用高階頻散信息，對該實測數據做兩步法反演： ①僅對第1段的頻散曲線進行單純的基階頻散曲線反演； ②將第2段高階頻散點與第①步反演所得模型的理論頻散曲線進行對比，確定其模態的歸屬，再對兩段頻散點做多階頻散曲線聯合反演，將此結果作為最終反演結果。反演的初始模型數據取自測井資料，將該區劃分為10個地層并設定搜索范圍(表5)。

表5 懷俄明M地區反演搜索范圍、密度和泊松比

對最終反演結果(圖9)進行對比分析。圖9a為頻散曲線對比圖，其中虛線為僅用基階面波頻散曲線反演的模型頻散曲線； 第2段實測頻散點歸屬于三階的頻散曲線； 實線為采用多階頻散曲線反演得到的模型頻散曲線，該頻散曲線與實測數據提取的兩段頻散點都能較好地對應。在橫波速度剖面(圖9b)中，多階(實線)比僅用基階(虛線)頻散曲線反演得到的橫波速度更接近于測井曲線(粉紅點線)。與Xia等[26]的局部線性優化結果(綠星點線)相比，在深度為8m范圍內兩者基本一致，但在較大深度上更能體現測井速度的變化，特別是在12～15m段內，該反演模型準確地刻畫出測井低速段。因此，ACGPSO算法中多階聯合反演相比于單純基階反演和局部線性化方法，能獲取地下介質的更精準地質信息，這也充分證明ACGPSO算法在實際資料處理中具有較好的適用性。

圖8 懷俄明地區實際地震數據與頻散能量圖(a)實測地震記錄； (b)頻散能量圖

圖9 懷俄明M地區頻散曲線反演結果(a)實測數據(圓點)、基階(虛線)和多階聯合(實線)反演模型頻散曲線； (b)橫波速度剖面與測井資料

4 結論

本文從四個方面改進了PSO算法的全局和局部搜索能力，提出了自適應混沌遺傳粒子群算法(ACGPSO)； 采用兩種測試函數檢測了新算法的優化能力，并將其運用于瑞雷波頻散曲線理論模型和實測數據的反演與分析中，取得了以下認識。

(1)ACGPSO算法對兩種測試函數均能100%收斂到其理論解，比常規PSO算法和GA算法具有更高的尋優求解能力和收斂速度。

(2)ACGPSO算法對無噪和含噪的頻散曲線反演均能取得很好效果，充分表明了該算法的有效性、穩定性和抗噪性， 同時也證明基于ACGPSO算法的反演比PSO算法具有更快的收斂速度，更強的穩定性和更高的反演精度。

(3)兩步法反演可較好地判斷高階頻散曲線對應的模態，以此進行的多模式聯合反演更能反映地下地層真實情況，更具可信度及適用性。

感謝夏江海教授團隊和中國地質大學(武漢)地球內部成像與探測實驗室提供的實測數據。