蚱蜢算法在瑞雷波頻散曲線反演中的應(yīng)用

于東凱 宋先海*② 張學(xué)強(qiáng) 趙素濤 蔡 偉

(①中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北武漢 430074；②中國(guó)地質(zhì)大學(xué)(武漢)湖北省地球內(nèi)部多尺度成像重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074)

0 引言

瑞雷波由縱波與橫波的干涉形成，沿著自由表面?zhèn)鞑ァＳ?guó)學(xué)者瑞雷首先發(fā)現(xiàn)并證明均勻半空間中瑞雷面波的存在[1]。瑞雷波勘探具有無(wú)損、高效、經(jīng)濟(jì)等特點(diǎn)，被視為未來(lái)最有希望的技術(shù)之一[2]，現(xiàn)已被廣泛用于淺地表、巖石圈橫波速度結(jié)構(gòu)信息[3-5]及地層品質(zhì)因子[6]的計(jì)算、路基路面檢測(cè)[7]、基巖面刻畫(huà)[8]等。與常規(guī)勘探方法相比，瑞雷波勘探具有分辨率高、應(yīng)用范圍廣、檢測(cè)設(shè)備簡(jiǎn)單、檢測(cè)速度快等優(yōu)點(diǎn)[9]。近年來(lái)，諸多學(xué)者進(jìn)行了相關(guān)研究，如瑞雷波反演[10-11]、瑞雷波頻散研究[12-14]、瑞雷波數(shù)值模擬研究[15-17]、地形對(duì)瑞雷波傳播的影響[18]等。

瑞雷波勘探可分為三個(gè)過(guò)程：面波數(shù)據(jù)采集、頻散曲線提取和頻散曲線反演[19-20]。拾取瑞雷波頻散曲線后，利用瑞雷波頻散特性可求得近地表橫波速度，是許多工程的關(guān)鍵參數(shù)之一[21-22]。多數(shù)局部線性化優(yōu)化方法被用來(lái)反演頻散曲線，并在各種商業(yè)軟件中占主導(dǎo)地位，如最速下降法、L-M算法等。然而，與大多數(shù)其他地球物理問(wèn)題一樣，瑞雷波頻散曲線反演是高度非線性、多參數(shù)、多極值的，局部線性化方法易收斂到局部極小值，反演結(jié)果在很大程度上取決于初始模型的選擇以及偏導(dǎo)數(shù)計(jì)算的準(zhǔn)確性。全局最優(yōu)化算法由于沒(méi)有上述限制，被越來(lái)越多地應(yīng)用于頻散曲線反演，如Yamanaka等[23]借助遺傳算法反演頻散曲線研究地下結(jié)構(gòu)信息； Beaty等[24]利用模擬退火算法反演多模式瑞雷波頻散曲線； Shirazi等[25]通過(guò)人工神經(jīng)網(wǎng)絡(luò)算法反演頻散曲線以獲得路面結(jié)構(gòu)； Song等[26-27]先后將模式識(shí)別算法、粒子群算法應(yīng)用于頻散曲線反演。傳統(tǒng)算法，如遺傳算法，在解決反演問(wèn)題中具有較高的應(yīng)用前景，然而往往存在計(jì)算量偏大、效率不高并且容易早熟收斂等缺陷，在反演中受到限制。

隨著智能算法的發(fā)展，新型算法不斷被應(yīng)用于地球物理反演[28-31]。本文針對(duì)傳統(tǒng)算法中出現(xiàn)的問(wèn)題，引入一種新型智能算法——蚱蜢算法，對(duì)瑞雷波頻散曲線進(jìn)行反演以獲取近地表橫波速度剖面。對(duì)不同地質(zhì)模型的理論頻散曲線進(jìn)行反演試驗(yàn)，并進(jìn)行多次、多角度、深層次的探索，最后利用實(shí)測(cè)數(shù)據(jù)檢驗(yàn)其實(shí)用性。

1 蚱蜢算法基本原理及驗(yàn)證

1.1 基本原理

受蚱蜢不同時(shí)期覓食行為的啟發(fā)，澳洲學(xué)者Saremi等[32]提出了一種新的仿生學(xué)算法——蚱蜢算法(Grasshopper Optimization Algorithm，GOA)，算法模擬了蚱蜢在幼蟲(chóng)和成蟲(chóng)兩個(gè)階段搜索食物時(shí)所表現(xiàn)出的不同行為。蚱蜢是自然界中一種常見(jiàn)的昆蟲(chóng)，它們成群結(jié)隊(duì)地出現(xiàn)，在幼蟲(chóng)時(shí)依靠強(qiáng)有力的雙腿進(jìn)行跳躍和移動(dòng)，雖移動(dòng)緩慢，但能夠較為詳盡地對(duì)所經(jīng)區(qū)域是否存在食物進(jìn)行判斷；當(dāng)其進(jìn)入成蟲(chóng)期后長(zhǎng)出翅膀，可在空中飛翔，移動(dòng)變得十分迅速，且移動(dòng)范圍迅速擴(kuò)大，搜索范圍得到大幅度的擴(kuò)大，但同時(shí)對(duì)食物的搜索精確度會(huì)得到一定程度的削弱。

蚱蜢算法將蚱蜢個(gè)體模擬為搜索算子，將成年蚱蜢大范圍且迅速移動(dòng)和幼蟲(chóng)期局部移動(dòng)模擬為算法尋優(yōu)過(guò)程中隨迭代次數(shù)的增加所進(jìn)行的全局探查與局部搜索兩個(gè)階段。在全局探查階段，算子在模型內(nèi)進(jìn)行大范圍的移動(dòng)，獲取全局信息，并最終固定為某局部區(qū)域；在局部搜索階段，算子的搜索范圍縮減至某局部區(qū)域，在該區(qū)域內(nèi)不斷迭代優(yōu)化所得解，不斷提高解的精度。算法主要包括兩個(gè)要素：自身所處位置以及與其他算子的位置關(guān)系。其中，自身位置體現(xiàn)了算子當(dāng)前的優(yōu)劣，蚱蜢與其他算子的位置關(guān)系可表示為：搜索算子的移動(dòng)受群體中其他算子的影響，當(dāng)算子相隔較近(距離小于2.079個(gè)單位)時(shí)會(huì)出現(xiàn)彼此排斥現(xiàn)象，此時(shí)稱蚱蜢處于排斥域范圍內(nèi)；當(dāng)算子相隔大于該范圍時(shí)又會(huì)出現(xiàn)吸引現(xiàn)象，稱其位于吸引域內(nèi)；當(dāng)距離恰好相差2.079個(gè)單位時(shí)，彼此之間不會(huì)發(fā)生吸引或排斥現(xiàn)象，蚱蜢位于舒適域。算法通過(guò)不斷調(diào)整自身位置及與其他算子的相對(duì)位置實(shí)現(xiàn)優(yōu)化。需要說(shuō)明的是，這兩個(gè)階段之間沒(méi)有明顯的界限，隨著搜素區(qū)間逐步縮小，算法由全局階段逐漸轉(zhuǎn)至局部階段。蚱蜢算法的優(yōu)化原理描述如下。

蚱蜢算子xi受蚱蜢算子xj的影響因子定義為

(1)

假設(shè)i個(gè)蚱蜢算子當(dāng)前位置為xi,k， 移動(dòng)時(shí)會(huì)受到所有其他算子的影響，因此移動(dòng)后的位置為

]}+xi,k

(2)

式中：i和j表示不同蚱蜢算子；NL表示算子總數(shù)；di,j表示第i個(gè)與第j個(gè)算子的空間距離；h和g均為評(píng)價(jià)算子受其他算子影響的參數(shù)，前者表示吸引力強(qiáng)度，后者是衡量空間距離的尺度，二者與影響因子F(di,j)緊密相關(guān)，不同組合會(huì)導(dǎo)致算子不同的行為習(xí)慣，算法中g(shù)通常取1.5，h取0.5；u(k)和l(k)分別表示算子的第k個(gè)分量在區(qū)間內(nèi)的上、下限；xi,k表示當(dāng)前蚱蜢算子所處位置；c1和c2為自適應(yīng)的遞減變量，c1與粒子群算法中慣性權(quán)重相似，作用是減少蚱蜢算子圍繞目標(biāo)值的移動(dòng)，c2作用是收縮算子群體的舒適域、排斥域和吸引域，通過(guò)對(duì)自身的調(diào)節(jié)可在有限迭代次數(shù)中達(dá)到全局探查以及局部搜索的相對(duì)平衡。

總體來(lái)說(shuō)，式(2)的前半部分實(shí)現(xiàn)了其他蚱蜢算子的約束作用，后半部分則對(duì)當(dāng)前算子移動(dòng)至全局最優(yōu)解的行為趨勢(shì)起到一定的控制作用。c1和c2雖在算法中所起作用不盡相同，但均為隨迭代進(jìn)行自適應(yīng)衰減的變量，可統(tǒng)一表示如下

(3)

式中：Iiter代表當(dāng)前迭代次數(shù)；NI為最大迭代次數(shù)；cmax為自適應(yīng)參數(shù)最大值，本文cmax取值為1；cmin為最小值，文中取1.0×10-5。可根據(jù)解決問(wèn)題的不同設(shè)定變量c1或c2的起始值。

本算法主要受三個(gè)參數(shù)控制：種群數(shù)量NL、迭代次數(shù)NI以及自適應(yīng)參數(shù)c1和c2。種群數(shù)量與迭代次數(shù)對(duì)算法的影響相似：隨著種群數(shù)量取值的增大，算法在變量空間中搜索會(huì)變得詳盡，但同時(shí)也可能出現(xiàn)過(guò)多的冗余搜索，影響算法的效率；若種群數(shù)量取值過(guò)小，則會(huì)出現(xiàn)搜索不充分，無(wú)法找到全局最優(yōu)解。因此，取值要根據(jù)所解決的具體問(wèn)題來(lái)設(shè)定。自適應(yīng)參數(shù)的取值與算法的尋優(yōu)功能緊密相關(guān)，具體來(lái)說(shuō)，為避免所得解陷入局部最優(yōu)，設(shè)置如下機(jī)制：蚱蜢群體對(duì)某算子會(huì)產(chǎn)生排斥或吸引現(xiàn)象，其中對(duì)算子的排斥作用能有效避免局部最優(yōu)解的產(chǎn)生，阻止算法陷入局部最優(yōu)以及停滯不前。同時(shí)，由于自適應(yīng)參數(shù)的引入，算法會(huì)縮小與迭代次數(shù)成比例的排斥面積，并不會(huì)對(duì)算法的收斂性造成過(guò)大影響，因此蚱蜢算子在開(kāi)始階段能有效避免局部最優(yōu)解，在最后的局部搜索階段也能改進(jìn)優(yōu)化所得解。

蚱蜢優(yōu)化算法的基本流程如下。

(1)參數(shù)初始化。在模型搜索空間內(nèi)隨機(jī)生成蚱蜢算子的初始位置xi(i=1,2,…,NL)，每個(gè)算子的位置可表示為xi=(xi,1,xi，2，…,xi，k,…,xi，D)，這里D表示算子中包含的變量個(gè)數(shù)。

(2)根據(jù)蚱蜢算子所處位置，計(jì)算各個(gè)算子的目標(biāo)函數(shù)值，并找出最優(yōu)解及相應(yīng)的算子位置。

(3)更新變量c1和c2，計(jì)算蚱蜢之間空間距離，并對(duì)所有蚱蜢算子遍歷，根據(jù)式(2)更新空間位置； 同時(shí)，若新位置躍出上下邊界，則給予約束。

(4)當(dāng)達(dá)到最大搜索次數(shù)或解滿足給定條件時(shí)，轉(zhuǎn)至步驟(5)； 否則，搜索次數(shù)增加1，轉(zhuǎn)至步驟(2)，進(jìn)行下一次搜索。

(5)輸出全局最優(yōu)解及其對(duì)應(yīng)的最優(yōu)擬合值。

1.2 算法測(cè)試

不同算法在處理問(wèn)題上的能力存在差異，同一算法處理不同問(wèn)題的能力也有高低。故在全局優(yōu)化算法領(lǐng)域，通常使用一些特定的函數(shù)對(duì)算法進(jìn)行測(cè)試，基于生成解的精度對(duì)算法的優(yōu)化能力進(jìn)行評(píng)價(jià)。本文主要研究將蚱蜢算法應(yīng)用于多極值的頻散曲線反演，多模態(tài)復(fù)合測(cè)試函數(shù)與反演中的目標(biāo)函數(shù)存在一定的相似性，即具有多個(gè)局部最優(yōu)解，因此選取兩個(gè)含有多峰值的函數(shù)進(jìn)行全局尋優(yōu)能力測(cè)試。

(1)Griewank函數(shù)。如圖1a所示(為方便顯示，僅顯示局部結(jié)果)，該函數(shù)存在許多局部極小值，其數(shù)目與問(wèn)題的維數(shù)有關(guān)，全局最小值0在(x1,x2,…,xD)=(0,0,…,0)處，因此該函數(shù)是典型的非線性多模態(tài)函數(shù)，具有廣泛的搜索空間，通常被認(rèn)為是優(yōu)化算法很難處理的復(fù)雜多模態(tài)問(wèn)題，在程序中搜索算子個(gè)數(shù)NL取500。其函數(shù)表達(dá)式為

xi∈[-600,600]

(4)

(2)Rastrigrin函數(shù)。如圖1b所示，該函數(shù)是一個(gè)多峰值函數(shù)，在(x1,x2,…,xD)=(0,0,…,0)處取得全局最小值0， 在{xi∈(-5.12, 5.12),i=1,2,…,D}范圍內(nèi)大約有10D個(gè)局部極小值，此函數(shù)與Griewank函數(shù)類似，也是一種典型的非線性多模態(tài)函數(shù)，峰形呈高低起伏不定跳躍性地出現(xiàn)，所以很難優(yōu)化查找到全局最優(yōu)解。其函數(shù)表達(dá)式為

xi∈[-5.12, 5.12]

(5)

測(cè)試中，種群數(shù)量NL取100，迭代次數(shù)NI取700。表1列出了蚱蜢算法對(duì)兩個(gè)測(cè)試函數(shù)的優(yōu)化結(jié)果，該結(jié)果為5次試驗(yàn)的平均值，可看出算法的所得解均收斂到函數(shù)最小值附近，收斂誤差接近于零，說(shuō)明算法具備較強(qiáng)的全局尋優(yōu)能力。圖2a為Griewank函數(shù)的算法擬合誤差隨迭代次數(shù)的變化，可見(jiàn)迭代逐漸收斂并最終趨近于零，說(shuō)明GOA法通過(guò)多次迭代能夠有效提高解的精度。另外，圖中也顯示出擬合誤差發(fā)生了不確定性變化，一方面是因?yàn)樗惴ㄖ幸肓素澙窚?zhǔn)則，保留當(dāng)前最優(yōu)解，算法在迭代之后所得解的精度未得到提高時(shí)會(huì)出現(xiàn)誤差不變的情況； 另一方面，高排斥率導(dǎo)致蚱蜢算子在全局范圍內(nèi)發(fā)生大范圍的移動(dòng)，使迭代誤差發(fā)生較大的改變。圖2b所示為算法在多峰值、多模態(tài)測(cè)試函數(shù)Rastrigrin函數(shù)中的表現(xiàn)，圖中每個(gè)點(diǎn)均代表一個(gè)蚱蜢算子在某次迭代中的搜索位置。可以看出：搜索過(guò)程中算子已基本覆蓋所有目標(biāo)區(qū)域，可見(jiàn)算法具有較高的搜索能力；圖中紅點(diǎn)表示全局最優(yōu)解，算子在該點(diǎn)附近大量聚集，說(shuō)明算法并未陷入局部最優(yōu)解，而是逐步逼近最優(yōu)解。另外，圖中有“抱團(tuán)”現(xiàn)象，即算子在某些區(qū)域聚集，原因可能是該點(diǎn)為局部最優(yōu)點(diǎn)，蚱蜢算子在該處附近迭代次數(shù)較其他位置稍大。

圖1 測(cè)試函數(shù)圖

圖2 蚱蜢算法中Griewank函數(shù)的迭代誤差曲線(a)和Rastrigrin函數(shù)的算子位置圖(b)

測(cè)試函數(shù)Griewank函數(shù)Rastrigrin函數(shù)理論極值00實(shí)際極值0.07880.0037算法所得解(4.17×10-4,7.01×10-3)(2.30×10-6,-1.35×10-6)

2 蚱蜢算法在面波反演中的實(shí)現(xiàn)

在瑞雷波勘探的應(yīng)用中，依據(jù)面波頻散資料推算工程場(chǎng)地的剪切波速度剖面時(shí)，通常認(rèn)為地下為層狀介質(zhì)。同時(shí)，基階波在面波頻散資料中能量最高、最易觀測(cè)，也應(yīng)用最為廣泛，故本文主要采用基階模式頻散曲線反演獲取地下層狀地質(zhì)結(jié)構(gòu)。夏江海等[6]的研究表明，與瑞雷波頻散曲線特征關(guān)系最密切的是地層的橫波速度(vS)和厚度(H)，這兩個(gè)變量可通過(guò)反演求取，地層密度、泊松比為已知參數(shù)(對(duì)頻散曲線影響微弱，可通過(guò)工區(qū)信息大致估算)。另外，為更符合實(shí)際情況(通常沒(méi)有足夠的先驗(yàn)信息估計(jì)淺地表的橫波速度和厚度準(zhǔn)確范圍)，文中使用較大的搜索范圍，在下文理論模型測(cè)試中，搜索區(qū)域的下限和上限設(shè)定為與真實(shí)值相差50%或更多。如前文所述，蚱蜢算法主要受到三個(gè)參數(shù)控制：NL、NI和c。反演程序中NL取值為10D；NI取值200；cmax和cmin分別取1和1.0×10-5。

頻散曲線反演是一個(gè)求解目標(biāo)函數(shù)最小值的優(yōu)化過(guò)程，算法的好壞取決于能否在給定范圍內(nèi)充分搜索以找到反演參數(shù)的最優(yōu)解。同時(shí)，由于反演中存在局部極值的干擾，算法能否跳出局部極值、找到全局最優(yōu)解也是評(píng)價(jià)其質(zhì)量的標(biāo)準(zhǔn)之一。目標(biāo)函數(shù)是通過(guò)算法反演得到的模型能否精確解釋觀測(cè)資料所確定的，應(yīng)用最為廣泛的均方差函數(shù)是

(6)

(7)

目標(biāo)函數(shù)具有以下幾個(gè)特點(diǎn)：

①目標(biāo)函數(shù)為M個(gè)頻點(diǎn)(即頻散曲線上的采樣點(diǎn))上的相速度共同作用。值得注意的是，采樣數(shù)目過(guò)大并不會(huì)對(duì)反演過(guò)程產(chǎn)生顯著影響，反而會(huì)影響反演效率[26]。由于所研究面波的頻散曲線頻率多在5～100Hz范圍，試驗(yàn)發(fā)現(xiàn)頻點(diǎn)數(shù)取25～30能平衡算法的高效性與采樣的保真度。

針對(duì)目標(biāo)函數(shù)的特點(diǎn)，蚱蜢算法在處理反演問(wèn)題上具有獨(dú)特優(yōu)勢(shì)，體現(xiàn)在算法能夠保證目標(biāo)函數(shù)不困于局部極值，而是進(jìn)行全局范圍的搜索；同時(shí)不斷優(yōu)化所得解，提高反演精度。算法引入排斥域、吸引域和舒適域機(jī)制：當(dāng)某位置上多個(gè)算子之間距離過(guò)小時(shí)(在反演中表現(xiàn)為迭代誤差相近)會(huì)相互排斥，算子在下次迭代中遠(yuǎn)離該位置，在算法初期需深入分析變量空間，避免受困于局部極值；當(dāng)距離過(guò)大時(shí)認(rèn)定進(jìn)入吸引域，算子之間會(huì)產(chǎn)生吸引力，體現(xiàn)在中后期的“開(kāi)發(fā)工作”中，隨著算法不斷運(yùn)行，大量的非可行區(qū)間被排除，搜索范圍會(huì)不斷被壓縮至最優(yōu)解附近，算法對(duì)此區(qū)域不斷開(kāi)發(fā)，對(duì)所得解不斷優(yōu)化以求提高反演精度；當(dāng)距離適當(dāng)時(shí)，算子彼此之間不會(huì)吸引或排斥，此時(shí)位于舒適域。為更好地平衡“開(kāi)發(fā)”(即全局探查)與“改進(jìn)”(即局部搜索)兩個(gè)過(guò)程，本算法配備了一個(gè)適應(yīng)性縮減三個(gè)作用域面積的系數(shù)，隨著迭代的進(jìn)行，算法會(huì)逐漸由全局探查轉(zhuǎn)為針對(duì)最優(yōu)解附近區(qū)域的局部搜索。由此可知，通過(guò)群體獲得的最佳解決方案是蚱蜢開(kāi)發(fā)和改進(jìn)的結(jié)果。因此，在解決一系列反演問(wèn)題時(shí)，蚱蜢算法具有超越當(dāng)前幾種算法的潛力。

2.1 理論模型反演

為檢驗(yàn)蚱蜢算法在瑞雷波頻散曲線反演中效果，使用兩個(gè)理論模型進(jìn)行測(cè)試，二者均為實(shí)際淺層工程勘察中經(jīng)常遇到的地質(zhì)模型。如表2所示，模型A代表橫波速度遞增的四層地質(zhì)模型；模型B代表在兩個(gè)較硬地層之間夾雜低速軟弱地層，該地層剪切波速度低于上覆地層速度，常見(jiàn)的有路面結(jié)構(gòu)等。

對(duì)理論模型進(jìn)行蚱蜢算法反演試驗(yàn)，模型A的數(shù)值模擬及反演結(jié)果見(jiàn)圖3和圖4，模型B的數(shù)值模擬及反演結(jié)果見(jiàn)5和圖6。圖3a、圖5a分別為兩個(gè)模型波場(chǎng)正演得到的地震記錄，以此模擬實(shí)測(cè)資料并作為頻散曲線反演的依據(jù)。圖4a、圖6a為擬合頻散曲線，圖3b、圖5b是通過(guò)高分辨率線性拉東變換提取的模擬地震記錄得到的頻散能量圖。值得注意的是，由于含低速軟弱夾層模型中頻散曲線的特殊性，即基階模式在高頻范圍內(nèi)(該模型中為頻率大于60Hz時(shí))攜帶的能量并不是最高的，因此文中拾取的基階波頻散曲線也有一定的范圍。

表2 兩種地質(zhì)模型參數(shù)及反演搜索范圍

圖3 模型A模擬地震記錄(a)及頻散能量圖(b)

圖4 基于蚱蜢算法的模型A反演結(jié)果

圖5 模型B模擬地震記錄(a)以及頻散能量圖(b)

圖6 基于蚱蜢算法的模型B反演結(jié)果

可以看出，在沒(méi)有足夠先驗(yàn)信息的情況下，頻散曲線吻合度較高，且模型中的參數(shù)均可被精確地反演和重建。就解的精確性而言，兩個(gè)模型中參數(shù)的最大相對(duì)誤差相對(duì)分別為2.50%、2.89%。由此可見(jiàn)，通過(guò)蚱蜢算法獲得的地下橫波速度與真實(shí)情況相差不大，說(shuō)明本算法可以對(duì)理論瑞雷波頻散曲線進(jìn)行很好地反演，對(duì)頻散曲線影響較大的參數(shù)(各層橫波速度以及層厚度)均能被精確地反演和重建。

2.1.1 頻散曲線聯(lián)合反演測(cè)試

瑞雷波是頻散且多階的，某頻率下的最小傳播速度被稱為基階模式，其余按相速度大小依次稱為一階高階模式、二階高階模式……。在瑞雷波頻散曲線的應(yīng)用中，由于基階波能量較強(qiáng)且易觀測(cè)等原因，通常使用基階波頻散曲線進(jìn)行非線性反演。但在某些地形條件下，如文中的模型B，高階模式的頻散曲線在高頻范圍內(nèi)攜帶了更多的能量，此時(shí)若只使用基階波進(jìn)行反演可能存在有效信息不足的問(wèn)題。為進(jìn)一步驗(yàn)證算法的反演能力，本次測(cè)試針對(duì)模型B，利用蚱蜢算法對(duì)基階和高階模式頻散曲線進(jìn)行聯(lián)合反演。聯(lián)合反演的頻散曲線共N條，故目標(biāo)函數(shù)相應(yīng)變?yōu)?/p>

(8)

圖7b所示為橫波速度剖面的對(duì)比。表3為基階波反演與聯(lián)合反演的對(duì)比，可以看出聯(lián)合反演結(jié)果較基階波反演得到進(jìn)一步優(yōu)化，各參數(shù)基本接近真實(shí)值，反演誤差趨于零，此結(jié)果進(jìn)一步證實(shí)蚱蜢算法反演精度較高。

圖7 基于蚱蜢算法的模型B聯(lián)合反演結(jié)果

反演方式vS1m·s-1vS2m·s-1vS3m·s-1vS4m·s-1vS5m·s-1H1mH2mH3mH4m基階波反演195.48147.14204.30304.01398.572.012.013.996.01聯(lián)合反演 199.27149.58201.29301.74399.502.002.014.006.01

2.1.2 含噪測(cè)試

由于野外環(huán)境復(fù)雜，采集到的地震記錄中不可避免含有噪聲，含噪數(shù)據(jù)可能降低算法精度及穩(wěn)定性，并影響反演結(jié)果。對(duì)一種新的反演算法需檢驗(yàn)其抗噪能力。本文對(duì)理論頻散曲線進(jìn)行加噪處理(加入5%、10%及20%噪聲)，即對(duì)各頻率下的相速度以理論值為基準(zhǔn)在一定區(qū)間內(nèi)進(jìn)行隨機(jī)程度的增減，以模擬實(shí)際情況中混雜了環(huán)境噪聲的相速度，含噪頻散曲線處理結(jié)果見(jiàn)圖8。

圖8a中雜亂不規(guī)則的藍(lán)色、藍(lán)綠色以及黑色實(shí)線所示為加入不同程度隨機(jī)噪聲的頻散曲線，圖8b為算法對(duì)含噪數(shù)據(jù)的反演結(jié)果。由圖8可見(jiàn)，橫波速度與真實(shí)模型之間有一定的偏離，三種不同噪聲水平的數(shù)據(jù)反演結(jié)果最大相對(duì)誤差分別增加到3.5%、9.5%和13.0%，由此可見(jiàn)噪聲對(duì)反演解的準(zhǔn)確性有一定的影響，但反演速度剖面均未實(shí)質(zhì)性地偏離真實(shí)模型，因此由算法得到的反演結(jié)果仍然具有較高的可信度，說(shuō)明蚱蜢算法對(duì)隨機(jī)噪聲具有一定的容忍度，用于反演瑞雷波頻散曲線具有較高的穩(wěn)定性。

2.1.3 搜索區(qū)間測(cè)試

前文設(shè)定的搜索范圍是根據(jù)前期測(cè)井等數(shù)據(jù)對(duì)整個(gè)工區(qū)層位進(jìn)行推斷，并以此估算地震記錄下對(duì)應(yīng)各層橫波速度及厚度的大致范圍。倘若地質(zhì)環(huán)境復(fù)雜，地層可能發(fā)生突變，則上述層位信息在當(dāng)前地震記錄下是不適用的，有必要設(shè)定更多數(shù)量的薄層以求準(zhǔn)確模擬近地表地質(zhì)情形。這是因?yàn)閼?yīng)用快速矢量傳遞算法計(jì)算頻散曲線時(shí)需要輸入各層橫波速度和厚度，比如n層模型需輸入(2n-1)個(gè)參數(shù)，故在反演之前有必要對(duì)模型層數(shù)進(jìn)行設(shè)定。本文為進(jìn)一步檢測(cè)算法的反演能力，以模型A為例，將模型層數(shù)設(shè)定為7，各層層厚均設(shè)定為2m，且橫波速度搜索范圍均設(shè)定為100～800m/s，頻散曲線反演是在給定的范圍內(nèi)找出最優(yōu)模型，因此要預(yù)先設(shè)定反演中橫波速度的搜索范圍，反演結(jié)果如圖9所示。由圖可見(jiàn)，在先驗(yàn)信息有限、搜索空間更為廣泛的前提下，反演過(guò)程耗時(shí)增加，反演結(jié)果仍十分接近真實(shí)值，并未發(fā)生顯著變化，說(shuō)明算法的全局搜索能力較強(qiáng)，反演結(jié)果具備較高的可信度。

圖8 含不同噪聲水平的模型A理論數(shù)據(jù)反演結(jié)果

圖9 擴(kuò)大搜索空間時(shí)基于蚱蜢算法的模型A反演結(jié)果

2.1.4 頻散曲線缺失頻段測(cè)試

在實(shí)地勘測(cè)中獲得的瑞雷波頻散曲線具有先天的不完整性，存在部分頻段受干擾或缺失等情況。實(shí)際上，在松散的沉積物上施加低頻源(如大錘)所產(chǎn)生的頻譜往往會(huì)缺乏高頻率數(shù)據(jù)。針對(duì)丟失頻段的瑞雷波頻散曲線反演時(shí)，首先應(yīng)考慮目標(biāo)層的敏感頻率范圍。本文針對(duì)模型A首先分析基階頻散曲線對(duì)各層橫波速度的敏感度。從圖10a可以看出，基階波對(duì)各層橫波速度變化敏感的頻段隨深度變化而變化，表層的敏感頻段為高頻，由于穿透深度的影響，深部層敏感頻段的頻率變得更低更窄，即面波的高頻成分對(duì)淺部地層的橫波速度較敏感，而低頻部分對(duì)深度層和半空間的橫波速度更敏感?；诖?，截取0～50Hz范圍的頻散曲線進(jìn)行基階波頻散曲線反演，通過(guò)對(duì)比橫波速度剖面分析蚱蜢算法反演能力。從圖10b可見(jiàn)，在缺失部分頻段的情況下，多次試驗(yàn)結(jié)果雖然出現(xiàn)一定程度的波動(dòng)，誤差較之前也增大，但均未實(shí)質(zhì)性地偏離真實(shí)值，且多次反演的均值與真實(shí)剖面吻合度較高，說(shuō)明算法的反演能力值得信賴。

2.2 反演結(jié)果分析

通過(guò)理論模型運(yùn)算及四種測(cè)試全方位地分析了蚱蜢算法的反演能力，其中穩(wěn)定性、收斂性和有效性是評(píng)價(jià)反演算法的重要指標(biāo)。為避免單次反演的偶然性，對(duì)上述測(cè)試均進(jìn)行20次獨(dú)立反演運(yùn)算，并取均值討論。受篇幅限制，只列出速度遞增模型A的反演結(jié)果統(tǒng)計(jì)表(表4)?？梢钥闯觯簩?duì)于無(wú)噪數(shù)據(jù)，模型反演結(jié)果與真實(shí)值相符，相對(duì)誤差和均方根誤差均趨近于零，說(shuō)明結(jié)果準(zhǔn)確可靠； 對(duì)于含噪測(cè)試中含噪數(shù)據(jù)(10%隨機(jī)噪聲)，模型的7個(gè)參數(shù)中相對(duì)誤差最大不超過(guò)9.5%，均方差最大不超過(guò)14.32，這說(shuō)明噪聲會(huì)對(duì)反演精度造成一定的影響，但所得解與真實(shí)解相差并不大，精度能滿足實(shí)際需求。原因在于蚱蜢算法引入舒適域、排斥域和吸引域三個(gè)作用域，算子會(huì)受到群體中其余所有算子的共同作用，全局尋優(yōu)能力更強(qiáng)，不易出現(xiàn)早熟收斂，能夠有效避免其他算法中由于隨機(jī)化個(gè)體單一、受當(dāng)前最優(yōu)解影響過(guò)大導(dǎo)致的全局尋優(yōu)能力不足等問(wèn)題。

圖10 基于蚱蜢算法的模型A反演結(jié)果(缺失頻段)

參數(shù)真實(shí)值不含噪聲含10%隨機(jī)噪聲反演均值相對(duì)誤差/%標(biāo)準(zhǔn)差反演均值相對(duì)誤差/%標(biāo)準(zhǔn)差vS1/(m·s-1)200.00203.841.921.57208.114.0614.32vS2/(m·s-1)250.00252.400.963.02257.402.965.60vS3/(m·s-1)300.00303.241.082.14290.243.2510.50vS4/(m·s-1)400.00396.810.803.26394.011.507.86H1/(m)2.001.952.501.202.042.001.02H2/(m)2.002.031.500.372.199.500.79H3/(m)4.004.061.500.414.082.000.85

另外，由搜索區(qū)間測(cè)試中的收斂曲線可以看出，無(wú)論搜索區(qū)間的上、下邊界是否靠近真實(shí)值，算法的最優(yōu)擬合值在初期的迭代過(guò)程中均迅速下降，說(shuō)明算法在此階段迅速排除搜索空間內(nèi)大量非可行區(qū)間，類比于蚱蜢在成蟲(chóng)期的覓食中能夠迅速排除非食物源區(qū)域。在隨后迭代過(guò)程中誤差下降趨勢(shì)變緩，且在后期趨近于零，說(shuō)明算法的反演過(guò)程在此之前已基本完成，各個(gè)參數(shù)已被精準(zhǔn)地推算。值得注意的是，任何一種群智能算法都會(huì)出現(xiàn)多次迭代中解的精度沒(méi)有得到提高的情況，這是因?yàn)樵谀硯状位蚨啻蔚^(guò)程中，種群中沒(méi)有找到更優(yōu)個(gè)體，根據(jù)貪婪準(zhǔn)則保留當(dāng)前最優(yōu)解，解的精度沒(méi)有得到提高，故而誤差不變。

同時(shí)，頻散曲線反演是一個(gè)非線性、多極值問(wèn)題，不同的模型參數(shù)組合也可能獲得相似的擬合值，因此還需要分析每次反演所得解中模型參數(shù)的分布概率，以驗(yàn)證反演結(jié)果的有效性。圖11所示為模型A中各個(gè)參數(shù)在真實(shí)值附近的分布概率。圖11a所示為第一層橫波速度vS1在20次反演結(jié)果中的分布概率，真實(shí)值為200m/s，反演所得解在真實(shí)值±10m/s范圍內(nèi)的概率接近0.80。由圖11可見(jiàn)，其余模型參數(shù)在其各自真實(shí)值附近(速度±10m/s，厚度±0.1m)的分布概率分別為0.75、0.60、0.70、0.65、0.70和0.65。可見(jiàn)反演所得模型參數(shù)均大致分布在真實(shí)值附近，進(jìn)一步說(shuō)明蚱蜢算法是一種可靠、準(zhǔn)確、收斂、有效的反演算法，能夠有效地應(yīng)用于瑞雷波頻散曲線的反演。

圖11 模型A反演參數(shù)的分布概率直方圖

3 與粒子群算法的對(duì)比

為驗(yàn)證蚱蜢算法在處理地球物理反演問(wèn)題中的適用性，將該算法與粒子群算法(Particle Swarm Optimization，PSO)進(jìn)行比較。粒子群算法是由Kennedy等[34]提出，后得到眾多學(xué)者的改進(jìn)和優(yōu)化[35]，已成為當(dāng)今最流行、最成熟的自然啟發(fā)優(yōu)化算法之一，并在地震勘探中已獲得廣泛應(yīng)用[36-38]。前人已成功將其應(yīng)用于頻散曲線反演[27]，本文就基于新型算法與粒子群算法的頻散曲線反演結(jié)果進(jìn)行對(duì)比，證明前者具有較高的可信度。

這兩種算法具備如下幾個(gè)相同點(diǎn)：均為全局優(yōu)化算法，算法在全局解空間內(nèi)進(jìn)行搜索，并將搜索重點(diǎn)集中在性能高的部分；都是隨機(jī)搜索算法，且具備記憶性能，最優(yōu)解都被保留下來(lái)；具有隱含并行性搜索特性，搜索過(guò)程均從解空間內(nèi)的一個(gè)集合開(kāi)始。但蚱蜢算法較粒子群算法而言，在處理地球物理反演問(wèn)題上具有不易陷入局部極值的優(yōu)勢(shì)。后者根據(jù)粒子本身歷史最優(yōu)位置和當(dāng)前群體內(nèi)的最優(yōu)算子的位置決定下一步的移動(dòng)方向，因此可能會(huì)出現(xiàn)如下現(xiàn)象：由于目標(biāo)函數(shù)存在大量的局部極小值，若最優(yōu)算子進(jìn)入某個(gè)極值位置，則算法接下來(lái)會(huì)在該位置附近多次移動(dòng)，反演效率不高。而蚱蜢算法中算子會(huì)受到群體內(nèi)其他所有算子的共同作用，體現(xiàn)在程序中算子每一步的移動(dòng)方向不僅受當(dāng)前最優(yōu)算子的影響，也會(huì)受到群體內(nèi)其余算子的共同作用，一定程度上降低了上述現(xiàn)象產(chǎn)生的可能性。同時(shí)，當(dāng)算法運(yùn)行至中后期時(shí)，算法通過(guò)縮小與迭代次數(shù)成比例的排斥面積，會(huì)減少排斥現(xiàn)象的產(chǎn)生，故不會(huì)對(duì)算法的收斂性造成過(guò)大的影響。

以模型B為例，分別采用上述兩種算法對(duì)模型B不含噪聲的理論數(shù)據(jù)進(jìn)行10次獨(dú)立反演測(cè)試。為排除其他因素的干擾，采用相同的群體規(guī)模、搜索空間和迭代次數(shù)，M取值30。測(cè)試分析一次反演的結(jié)果而不是10次反演迭代誤差的均值，目的是更準(zhǔn)確直觀地分析算法的性能。反演結(jié)果如圖12所示，由圖可以看出，如果迭代截至50～100次(圖12a)，迭代誤差的整體變化情形與算法優(yōu)化所得反演解的能力基本一致，兩種算法效果不相上下，說(shuō)明這兩種算法均具備較強(qiáng)的收斂能力和收縮搜索區(qū)間性能；但迭代進(jìn)入中后期(圖12b)，即兩種算法運(yùn)行100～200次，蚱蜢算法比粒子群算法收斂速度更快，且收斂精度也有一定提高，說(shuō)明蚱蜢算法在處理地球物理反演問(wèn)題時(shí)具備較高的尋優(yōu)性能。

圖12 模型B不含噪聲時(shí)的GOA與PSO迭代誤差對(duì)比

4 實(shí)測(cè)資料反演

前文充分驗(yàn)證了蚱蜢算法在理論瑞雷波頻散曲線反演中可行且有效，所得解精度高。為進(jìn)一步檢驗(yàn)算法的適用性，利用蚱蜢算法對(duì)美國(guó)懷俄明某地獲得的地震數(shù)據(jù)進(jìn)行了分析。圖13a所示為采集的地震記錄，數(shù)據(jù)采集采用了48個(gè)8Hz垂直分量檢波器，道間距為0.9m，最小炮檢距為0.9m，震源采用錘擊震源。圖13b所示為對(duì)應(yīng)的頻散能量圖，通過(guò)采集能量最大點(diǎn)提取基階波頻散曲線。與反演理論模型類似，固定地層的密度及泊松比，只反演橫波速度vS與厚度H兩個(gè)參量。根據(jù)測(cè)井?dāng)?shù)據(jù)將地層大致劃分為5個(gè)層位，并粗略估計(jì)搜索范圍。反演模型的搜索空間、泊松比以及密度等參數(shù)設(shè)置見(jiàn)表5。

圖14a所示為算法所得解與實(shí)測(cè)數(shù)據(jù)的基階波頻散曲線擬合情況，通過(guò)拾取頻散能量圖中能量最大值(圖13b中圓點(diǎn))得到不同頻率下的相速度，可以看出頻散曲線擬合情況較好。文中采用基階波反演的原因在于能量圖中出現(xiàn)較為嚴(yán)重的“模式接吻”以及“模式跳躍”現(xiàn)象，在拾取高階頻散曲線的過(guò)程中很容易出現(xiàn)模式誤判，從而影響反演結(jié)果，若采用基階波反演精度會(huì)更高。圖14b為兩種算法反演時(shí)迭代擬合誤差隨迭代次數(shù)的變化曲線，為保證反演的精準(zhǔn)性、降低因單次反演中出現(xiàn)的偶然誤差對(duì)結(jié)果造成影響，針對(duì)野外數(shù)據(jù)進(jìn)行了多次反演，舍棄異常的反演解，并對(duì)剩余解進(jìn)行均值化處理，因此迭代誤差曲線較圖12更為平滑。圖14c為橫波速度剖面對(duì)比，可見(jiàn)通過(guò)有限頻段的頻散曲線反演所得速度模型與測(cè)井資料吻合較好，尤其是淺地表部分。相較于PSO算法，GOA算法能夠準(zhǔn)確地分辨厚度為3～5m的軟弱夾層。值得注意的是，較前文測(cè)試，此次實(shí)測(cè)資料反演的信息量有限，與模型試算結(jié)果相比反演效果并不突出，但借助有限的頻散信息所得的反演結(jié)果能基本滿足實(shí)際工作需求。另外，較深地層反演速度剖面與測(cè)井資料出現(xiàn)了一定偏離，這是由于采用錘擊震源，瑞雷波頻散能量在低頻段分辨率低，導(dǎo)致提取誤差增大；另外一個(gè)原因是算法只設(shè)定5個(gè)層位，較深的地層被認(rèn)定為均勻介質(zhì)，即橫波速度沒(méi)有變化。

表5 GOA反演中模型搜索范圍及模型參數(shù)

圖13 懷俄明某地區(qū)地震記錄(a)與頻散能量圖(b)

圖14 懷俄明某地區(qū)地震數(shù)據(jù)蚱蜢算法反演結(jié)果

5 結(jié)論

本文介紹了一種全新的仿生學(xué)算法——蚱蜢算法，并通過(guò)兩個(gè)測(cè)試函數(shù)驗(yàn)證其全局搜索能力，證明了其優(yōu)異的收斂能力以及全局尋優(yōu)能力。將蚱蜢算法應(yīng)用于瑞雷波頻散曲線反演，通過(guò)設(shè)置較大的搜索范圍以模擬實(shí)際工作中缺少足夠先驗(yàn)信息的情況，并進(jìn)行多角度、深層次的探索試驗(yàn)，深入分析了算法的反演能力。最后通過(guò)實(shí)測(cè)數(shù)據(jù)檢測(cè)了本算法的適用性。理論數(shù)據(jù)和實(shí)測(cè)資料反演的結(jié)果表明：

(1)蚱蜢算法具有較高的全局尋優(yōu)能力以及抗局部最優(yōu)解能力，所得解與函數(shù)最小解十分接近，說(shuō)明這是一種有效的全局優(yōu)化算法；

(2)蚱蜢算法可應(yīng)用于瑞雷波頻散曲線反演，與傳統(tǒng)算法相比，反演過(guò)程中目標(biāo)函數(shù)曲線收斂更加穩(wěn)定，最終迭代誤差更接近于零。證明基于蚱蜢算法的瑞雷波頻散曲線反演具有較高的精度，是一種簡(jiǎn)單、穩(wěn)健的反演算法。

同時(shí)，蚱蜢算法也存在一定的缺陷，比如算法的收斂過(guò)程相對(duì)較長(zhǎng)。在今后的研究中，可與其他算法如遺傳算法、粒子群算法或人工神經(jīng)網(wǎng)絡(luò)算法相結(jié)合，探討如何進(jìn)一步提高算法的收斂效率。