基于Monte carlo法的Rayleigh波淺層橫波速度結(jié)構(gòu)反演方法研究

胡昌榮, 黃金強(qiáng),2, 劉 宏,2, 竇勇碩, 何云川

(1.貴州大學(xué) 喀斯特地質(zhì)資源與環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，貴陽 550025；2.貴州大學(xué) 資源與環(huán)境工程學(xué)院，貴陽 550025)

0 引言

近地表橫波速度值，是公路橋隧建筑等抗震設(shè)計(jì)、土壤級(jí)別劃分及地層軟硬程度判定的關(guān)鍵參數(shù)之一，在淺層巖土工程勘察、設(shè)計(jì)與施工中應(yīng)用十分廣泛[1-3]。目前，橫波速度反演主要有：基于多波多分量數(shù)據(jù)的彈性波層析建模方法、基于頻散特征的面波方法以及近幾年提出的面波全波形反演方法。其中面波全波形反演具有反演質(zhì)量高、不受硬表層及低速半空間影響等優(yōu)勢(shì)，但卻存在計(jì)算成本高、反演效果受限于數(shù)據(jù)品質(zhì)及方法并未完全實(shí)用化等缺點(diǎn)；多波多分量彈性波法因震相拾取困難，大多應(yīng)用于孔間地層速度反演，對(duì)施工要求較高；而基于頻散特征的面波法不僅采集方式簡(jiǎn)單易行，操作方便，還具有數(shù)據(jù)處理方法較為穩(wěn)定可靠的優(yōu)勢(shì)。然而，基于非線性全局優(yōu)化反演算法的面波反演策略，仍然存在諸如計(jì)算效率低、收斂迭代慢和內(nèi)存占用高等問題，為此，進(jìn)一步發(fā)展一種快速且穩(wěn)定的反演方法就顯得極為重要。

Rayleigh波勘探技術(shù)具有非損傷性、低成本及勘探高效性等優(yōu)點(diǎn)，備受淺地表地球物理及工程地質(zhì)學(xué)界重視，是當(dāng)前最具潛力的勘探技術(shù)之一[4]。面波譜分析方法(SASW)最早由Nazarian等[5]提出，但該方法在實(shí)際淺地表工程應(yīng)用中存在提取頻散曲線的抗干擾能力差等缺陷。隨后，美國堪薩斯大學(xué)Park等[6]進(jìn)一步提出面波多道分析技術(shù)(MASW)，該方法有效地彌補(bǔ)了SASW方法的不足。由于基于主動(dòng)源的面波多道分析技術(shù)在提取頻散曲時(shí)具有精度低等問題，Park等[7]將主動(dòng)源及被動(dòng)源優(yōu)勢(shì)相聯(lián)合，進(jìn)一步發(fā)展和系統(tǒng)研究了MASW技術(shù)。在MASW技術(shù)中，高頻Rayleigh波勘探技術(shù)大致可分為4個(gè)具體步驟：①面波數(shù)據(jù)采集；②數(shù)據(jù)預(yù)處理與頻散曲線提?。虎鄢跏寄Ｐ凸烙?jì)并進(jìn)行頻散曲線的正演計(jì)算；④頻散曲線反演橫波速度。由于Rayleigh波在非均勻介質(zhì)中傳播具有頻散特征，通過Xia等[8]研究結(jié)果表明，Rayleigh波頻散方程對(duì)橫波速度及層厚最為敏感，因此利用Rayleigh波頻散特征來反演淺地表橫波速度與層厚具有明顯優(yōu)勢(shì)。但是在一些特殊模型中，頻散方程對(duì)模型參數(shù)的敏感性較低，反演還存在一定的技術(shù)問題：①在低速半空間中Rayleigh波頻散方程在某些頻率求解下存在復(fù)數(shù)根，導(dǎo)致在某些頻率下無法正常獲取Rayleigh波相速度，因此，Pan等[9]提出了用復(fù)數(shù)根實(shí)部作為Rayleigh相速度，有效解決了該問題，但低速半空間的反演問題是一個(gè)復(fù)雜問題，局部線性化方法反演低速半空間橫波速度的適用性將大大降低；②若表層的橫波速度高于下伏地層時(shí)，高頻處的Rayleigh相速度將趨近于除半空間以外其他層的最小橫波速度值，而不是趨近表層橫波速度相關(guān)的某個(gè)值，Pan等[10]提出采用“替代”的方法來解決Rayleigh波波長(zhǎng)不足以穿透最低橫波速度的層厚時(shí)引起地計(jì)算錯(cuò)誤，但是基于Monte Carlo反演策略在反演高速表層模型時(shí)結(jié)果卻并不令人滿意；③地下介質(zhì)表現(xiàn)為橫向不均勻時(shí)，在求解水平層狀介質(zhì)中Rayleigh波的本征值問題之上開展的Monte carlo局部線性化反演算法，在反演橫向不均勻介質(zhì)時(shí)，反演效果不佳。

MASW技術(shù)將采集得到的t-x域多道Rayleigh波信號(hào)通過一定方法變換到f-v域，從而提取Rayleigh波頻散能量(如：f-k變換法[11]、高線性拉東變換[12]、τ-p變換法[13]、相移法[14]等)。Rayleigh波頻散曲線是淺地表地下特征的綜合反映，在利用Rayleigh波頻散特性反演橫波速度中，分別為非線性全局優(yōu)化算法及局部線性化反演算法兩大類?；诰植烤€性化反演算法在求解非線性問題時(shí)，每一次迭代過程都將朝向目標(biāo)函數(shù)取極小值，為快速求解全局極小值，前人們提出了各種最優(yōu)化算法，Xia[8]采用改進(jìn)的阻尼最小二乘反演策略提高了Rayleigh波頻散曲線反演的計(jì)算效率，但局部線性化反演算法強(qiáng)烈依賴初始模型，若初始模型選取不當(dāng)易陷入局部極小值。因此，學(xué)者們進(jìn)一步發(fā)展了不依賴于初始模型的非線性全局優(yōu)化反演算法，提出了徹底搜索法、Monte carlo反演算法、改進(jìn)的Monte carlo反演算法等。但是在非線性全局優(yōu)化反演算法中，徹底搜索法對(duì)模型參數(shù)所組合到的所有模型進(jìn)行空間搜索，存在計(jì)算機(jī)的海量計(jì)算[15]。而基于概率論觀點(diǎn)的Monte carlo反演算法[16]，以隨機(jī)方式進(jìn)行模型的空間搜索，有效避免了計(jì)算機(jī)的海量計(jì)算，在非線性全局優(yōu)化反演算法中得到了廣泛應(yīng)用[17-18]，但在Monte carlo反演算法中，若沒有足夠的先驗(yàn)信息引導(dǎo)隨機(jī)搜索，將存在計(jì)算機(jī)的“盲目”及完全隨機(jī)搜索，極大增加了計(jì)算量。為此，隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，一些改進(jìn)的Monte carlo反演算法也得到了廣泛地應(yīng)用(如遺傳算法、退火算法、BP神經(jīng)網(wǎng)絡(luò)等)。但此類方法仍然存在計(jì)算效率低、迭代慢、計(jì)算機(jī)內(nèi)存占用高等問題。如：石耀霖等[19]用非線性遺傳算法反演地球內(nèi)部構(gòu)造時(shí)，該算法存在不能保證找到最優(yōu)解的缺點(diǎn)，同時(shí)易出現(xiàn)“早熟收斂”或收斂速度極慢等問題?；凇皢l(fā)”式隨機(jī)搜索的模擬退火算法對(duì)地層橫波各向異性進(jìn)行反演時(shí)，存在計(jì)算速度慢等問題[20]?；贐P神經(jīng)網(wǎng)絡(luò)的Rayleigh波智能優(yōu)化反演橫波速度時(shí)，存在網(wǎng)絡(luò)訓(xùn)練耗時(shí)較長(zhǎng)以及針對(duì)不同地震模型需要重新網(wǎng)絡(luò)訓(xùn)練等繁瑣問題，其次若訓(xùn)練樣本過大也將影響收斂速度等缺點(diǎn)[21]。因此，在求解水平層狀介質(zhì)中Rayleigh波的本征值問題之上，采用Monte carlo局部線性化反演算法，該算法不僅具有計(jì)算速度快、計(jì)算機(jī)內(nèi)存占用小等優(yōu)勢(shì)，還有效避免了Monte carlo反演算法中的“盲目”及完全隨機(jī)搜索。

綜上所述，筆者在MASW技術(shù)中，借鑒Press的思想，在求解水平層狀介質(zhì)中Rayleigh波的本征值問題之上，采用Monte carlo局部線性化反演算法獲取淺地表橫波速度值[22]，并進(jìn)一步在典型數(shù)值算例中，對(duì)比分析了基于Monte Carlo反演策略與經(jīng)典阻尼最小二乘反演策略，驗(yàn)證了基于Monte Carlo反演策略在求取淺層橫波速度方面的優(yōu)勢(shì)。

1 方法原理

1.1 Monte carlo反演算法

改進(jìn)的Monte carlo反演算法在獲取淺地表橫波速度中存在固有的計(jì)算效率低、迭代慢、計(jì)算機(jī)內(nèi)存占用高等問題(如遺傳算法、BP人工神經(jīng)網(wǎng)絡(luò)、退火算法等)。因此借鑒Press的思想，在求解水平層狀介質(zhì)中Rayleigh波的本征值問題之上，采用Monte carlo局部線性化反演算法獲取淺地表橫波速度值可避免此類問題。在Monte carlo反演算法中為避免“盲目”及完全隨機(jī)搜索，采用空間約束方法，對(duì)解的范圍進(jìn)行空間約束，這極大提高M(jìn)onte carlo反演算法的計(jì)算效率，具體推導(dǎo)過程如下：

在求解水平層狀介質(zhì)中Rayleigh波的本征值問題上，可將水平層狀介質(zhì)模型中的每層橫波速度值、層厚值、密度值及縱波速度值分別定義為Vs、h、ρ、Vp的集合。

Vs={Vs1,Vs2,…,Vsn}

h={h1,h2,…,hn}

ρ={ρ1,ρ2,…,ρn}

Vp={Vp1,Vp2,…,Vpn}

(1)

通過改進(jìn)的Haskell-Thomson傳遞矩陣算法正演頻散曲線[22](詳見附錄A)可計(jì)算Rayleigh波理論頻散曲線。

(2)

其中：c為Rayleigh波相速度;k為波數(shù);n為層數(shù);i為給定波數(shù)所對(duì)應(yīng)的Rayleigh波相速度的點(diǎn)數(shù)。

因此在Haskell-Thomson傳遞矩陣正演算法的基礎(chǔ)上，已知實(shí)測(cè)頻散曲線反推地質(zhì)模型。首先，根據(jù)實(shí)測(cè)頻散曲線的形態(tài)特征邏輯推斷反演的初始模型Vs、h、Vp、ρ(圖1)。其次，定義當(dāng)前初始模型Vs、h的空間范圍：

(3)

圖1 實(shí)測(cè)頻散曲線的參數(shù)估計(jì)Fig.1 Parameter estimation of measured dispersion curve(a)實(shí)測(cè)頻散曲線;(b)模型參數(shù)

Δj=Ω(-γ·Vs,γ·Vs)

Δl=Ω(-δ·h,δ·h)

式中：Δj、Δl分別表示橫波速度及層厚度模型的采樣空間，由于在采樣空間中能隨機(jī)抽取到無數(shù)個(gè)值，因此，j表示在采樣空間中某一層橫波速度能取到的值的個(gè)數(shù)，l表示在采樣空間中與橫波速度所對(duì)應(yīng)的地層中層厚度能取到的值的個(gè)數(shù)，γ、δ為常數(shù)，將對(duì)橫波速度及層厚度采樣空間進(jìn)行約束，即：對(duì)橫波速度及層厚度的解的空間范圍進(jìn)行約束。

隨后，在一定空間范圍內(nèi)隨機(jī)生成每一次正演的樣本模型:

Vsj=Vs+Δj

hl=h+Δl

(4)

為獲取相應(yīng)地層中的Vp值，根據(jù)Everett[23]提出的縱橫波速度比與泊松比的關(guān)系進(jìn)行縱波速度值推算：

(5)

通過密度與縱波速度經(jīng)驗(yàn)公式對(duì)相應(yīng)地層中的密度值進(jìn)行計(jì)算[24]：

(6)

最后，將隨機(jī)生成的每一次迭代樣本模型{Vsj,Vpj,ρj,hl}，利用改進(jìn)的Haskell-Thomson傳遞矩陣算法進(jìn)行正演，以此得到每一個(gè)模型的頻散曲線。

1.2 阻尼最小二乘反演算法

Rayleigh波頻散方程與模型參數(shù)Vs、Vp、ρ、h有關(guān)，但頻散方程對(duì)橫波速度Vs的敏感性極高，因此阻尼最小二乘為減少待反演的模型參數(shù)，在反演時(shí)將其余參數(shù)Vp、ρ、h作為已知變量，并在反演過程中保持不變，只將Vs作為反演過程中的獨(dú)立變

量，這極大提高了反演的計(jì)算效率。為此進(jìn)一步將頻散方程式(2)簡(jiǎn)述為式(7)[8]。

(7)

由于大多數(shù)地球物理問題并非線性問題，為使非線性問題線性化，根據(jù)泰勒級(jí)數(shù)定理，將式(7)中ci在初始值ci,0處展開并取一階近似可得：

(8)

進(jìn)一步采用矩陣方程形式將式(8)表述為式(9)。

(9)

表述的矩陣方程式(9)可進(jìn)一步簡(jiǎn)述為：

Δc=AΔx

(10)

式中：Δc=ci-ci,0,i=1、2、…、M，表示M維長(zhǎng)度的Rayleigh波相速度殘差向量;A表示一階偏導(dǎo)數(shù)組成的M×n維雅格比矩陣;Δx=Vsq-Vsq,0,q=1、2、…、n，表示n維長(zhǎng)度的橫波速度殘差向量。

針對(duì)線性最小二乘反演問題，可將目標(biāo)函數(shù)表述為若干個(gè)函數(shù)的平方和形式:

(11)

式中：ciobs表示實(shí)測(cè)Rayleigh波相速度值，cical表示理論Rayleigh波相速度值。

為使最小二乘的目標(biāo)函數(shù)最小化，采用L-M方法進(jìn)行改進(jìn)[8,15,25]，根據(jù)L-M方法的原理，將式(10)帶入式(11)，可將目標(biāo)函數(shù)表述為矩陣的形式：

φ(Vs)=‖AΔx-Δc‖2W‖AΔx-

(12)

式中：‖‖2表示l2范數(shù);W表示權(quán)重矩陣;λ為阻尼因子。

進(jìn)一步將目標(biāo)函數(shù)簡(jiǎn)化并整理得：

(ATA+λI)Δx=ATΔc

(13)

式中：AT表示雅格比矩陣A的轉(zhuǎn)置;I為單位矩陣。

通過移項(xiàng)并整理得到方程(12)的解：

Δx=(ATA+λI)-1ATΔc

(14)

式(14)中，通過引入阻尼因子λ，可避免“病態(tài)”系數(shù)矩陣導(dǎo)致的求解過程不穩(wěn)定，同時(shí)也控制著Δx的大小，由于λ的每一次變化，將引發(fā)計(jì)算的繁瑣，因此采用奇異值分解方法計(jì)算廣義逆矩陣，將雅格比矩陣A奇異值分解表述為：

A=UΛVT

(15)

將式(15)帶入式(13)可得式(16)。

(16)

即有：

(17)

進(jìn)一步整理得：

(18)

從式(18)中可見，在對(duì)角矩陣中，λ的每一次變化，不會(huì)引發(fā)繁瑣的計(jì)算。

(19)

(20)

1.3 Monte carlo 檢驗(yàn)準(zhǔn)則

采用均方根誤差定義每一條頻散曲線與實(shí)測(cè)頻散曲線的誤差，將均方根誤差[26]作為樣本模型的檢驗(yàn)標(biāo)準(zhǔn)，若誤差達(dá)到最小時(shí)，則接收該樣本模型是最優(yōu)模型，反之，則排斥該樣本模型。具體反演算法流程如圖2所示。

(21)

2 模型測(cè)試/數(shù)值模擬

2.1 三層模型實(shí)驗(yàn)

筆者首先利用三層層狀模型，通過Rayleigh波有限差分正演模擬合成單炮面波正演記錄，隨后通過相移法[14]提取頻散曲線(詳見附錄B)，最后通過設(shè)置兩組不同的初始模型，分別對(duì)比了基于Monte carlo與基于經(jīng)典阻尼最小二乘兩種Rayleigh波反演策略在橫波速度建模中的優(yōu)缺點(diǎn)。

三層模型參數(shù)如表1所示，由表1中可知，橫波速度較小，泊松比高達(dá)0.47，在工程物探中主要代表未固結(jié)的松散碎屑物質(zhì)(如黏土層、飽和沙土等)。在高泊松比地質(zhì)結(jié)構(gòu)中，由于傳統(tǒng)的PML邊界條件在Rayleigh波正演數(shù)值模擬中將產(chǎn)生數(shù)值不穩(wěn)定?；诖?，采用多軸卷積完全匹配層吸收邊界條件處理人工截?cái)噙吔鏪27]，應(yīng)力鏡像+二階速度展開處理自由表面邊界[28]，數(shù)值模擬方法為交錯(cuò)網(wǎng)格有限差分正演模擬[29-30]。縱橫向網(wǎng)格大小為2001×601，縱橫向網(wǎng)格間距均為0.1 m，實(shí)際模擬區(qū)域?yàn)?00 m×60 m，邊界層厚度為801個(gè)網(wǎng)格點(diǎn)，時(shí)間采樣間隔為0.05 ms，采樣點(diǎn)數(shù)為60 001，總模擬時(shí)間為3 s，震源為主頻8 Hz的雷克子波，震源位置位于(0，0)m，由于勘探目的層的深度應(yīng)大于25 m，本算例共設(shè)置61道接收，接收點(diǎn)深度為0 m，第一個(gè)檢波器的位置為70 m，道間距為2 m，因此最小偏移距為70 m，最大偏移距為190 m。

圖2 Monte Carlo反演流程圖Fig.2 Flowchart of Monte Carlo inversion

表 1 三層模型參數(shù)表

利用上述參數(shù)，通過數(shù)值模擬得到如圖3所示的垂向位移分量的地震記錄。從圖3中可知：①地震記錄中沒有明顯的虛假邊界反射；②在時(shí)間較大時(shí)，正演模擬仍然保持穩(wěn)定，這是由于本文采用了多軸卷積完全匹配層吸收邊界條件；③從記錄中可觀察到縱波走時(shí)快、能量低，面波走時(shí)慢、能量強(qiáng)，呈現(xiàn)“掃帚狀”的頻散特征，表明在自由表面采用應(yīng)力鏡像法的處理方式是合理的。進(jìn)一步，根據(jù)相移法頻散曲線提取流程，以圖3(a)所示的地震記錄作為輸入，計(jì)算得到如圖3(b)所示的頻散能量圖，從圖3中可知：①在高頻部分，由于波長(zhǎng)較短，主要受到淺部地層的影響，因此，頻散能量圖中的高頻能量峰值主要指示淺層面波速度，數(shù)值約為96 m/s，趨近淺層橫波速度的0.94倍；②在低頻部分，由于波長(zhǎng)較長(zhǎng)，面波速度將同時(shí)受到來自淺中深地層的影響，故圖3(b)中的能量峰值線能夠反映深層面波速度，數(shù)值約為310 m/s，趨近于半空間橫波速度的0.90倍；③如果地層情況較為復(fù)雜，則頻散曲線是地下介質(zhì)的綜合反映。

圖3 Rayleigh波合成記錄及頻散能量圖Fig.3 Synthetic record and dispersive image of Rayleigh waves(a)垂向位移分量的地震記錄；(b)頻散能量圖

表2 初始模型參數(shù)值

隨后，從圖3(b)中提取實(shí)測(cè)頻散曲線作為反演的觀測(cè)數(shù)據(jù)，并設(shè)計(jì)三組初始模型，前Ⅰ、Ⅱ組模型為等厚薄層模型，第Ⅲ組為非等厚薄層模型，前Ⅰ、Ⅱ組的層數(shù)均為9層，前8層厚度均為4 m，第9層為半空間，由于縱波速度及密度對(duì)Rayleigh波頻散曲線影響較小，因此，縱波速度與密度設(shè)置為常數(shù)，分別為517 m/s和1 500 g/cm3。三組模型的初始橫波速度設(shè)置也較為簡(jiǎn)單，并未采用復(fù)雜的拐點(diǎn)分層方案，而是將第I組模型的橫波速度設(shè)定為線性遞增速度，其中模型第1層的橫波速度由頻散能量圖中高頻能量峰值所指示的橫波速度所決定，其數(shù)值為106 m/s，每層速度的增量為24 m/s，第II組設(shè)為均勻模型，橫波速度值為253 m/s，第Ⅲ組模型的第一層為106 m/s，每層速度為非線性遞增速度。模型參數(shù)見表2。通過對(duì)比三組初始模型與真實(shí)模型可知：兩者差異較大，可用于對(duì)比檢驗(yàn)Monte carlo反演算法與阻尼最小二乘反演算法對(duì)初始模型的依賴性與反演的穩(wěn)定性。

為使反演結(jié)果達(dá)最優(yōu)，第I組模型中Monte carlo反演的空間約束參數(shù)、分別為0.15、0.112 5，而第II組中、分別為0.125、0.1，第Ⅲ組中γ、δ均為0.1。經(jīng)統(tǒng)計(jì)，第I組模型中，Monte carlo反演迭代1 000次，耗時(shí)104 s，經(jīng)典阻尼最小二乘反演迭代20次，耗時(shí)122 s。而第II組模型中，Monte carlo反演迭代1 000次，耗時(shí)109 s，經(jīng)典阻尼最小二乘迭代25次，耗時(shí)304 s，第Ⅲ組模型中，Monte carlo反演迭代1 000次，耗時(shí)28.8 s，經(jīng)典阻尼最小二乘反演迭代10次，耗時(shí)110.7 s。實(shí)驗(yàn)結(jié)果表明：Monte carlo反演策略相比經(jīng)典阻尼最小二乘具有較快的收斂速度。

經(jīng)典阻尼最小二乘反演策略及Monte carlo反演策略從本質(zhì)上來看是基于實(shí)測(cè)頻散數(shù)據(jù)的擬合，為分析經(jīng)典阻尼最小二乘及Monte carlo的反演質(zhì)量，在三組初始模型中分別計(jì)算了經(jīng)典阻尼最小及Monte carlo反演的理論頻散，并進(jìn)行了理論頻散曲線與實(shí)測(cè)頻散曲線的均方根誤差分析。

圖4表示在第I組等厚薄層線性遞增模型中兩種方法的反演結(jié)果及頻散曲線的對(duì)比。為更清晰刻畫在等厚薄層線性遞增模型中的反演細(xì)節(jié)，圖4(a)給出了經(jīng)典阻尼最小二乘反演的橫波速度剖面與真實(shí)模型的對(duì)比，從圖4(a)中可見經(jīng)典阻尼最小二乘存在如下缺點(diǎn)：①經(jīng)典阻尼最小二乘反演法在保證反演精度的條件下，為提高反演解的唯一性，采用固定約束條件，在初始模型為等厚薄層中，層參數(shù)在每次迭代過程保持不變，減少了待反演的參數(shù)，只有橫波速度作為獨(dú)立變量；②在層速度突變界面，反演的橫波速度剖面與真實(shí)模型存在較大差異，未能較好地評(píng)價(jià)地下介質(zhì)速度結(jié)構(gòu)；③對(duì)初始模型依賴性較高。圖4(c)給出了Monte carlo反演的橫波速度剖面與真實(shí)模型的對(duì)比，從圖4(c)中可知，Monte carlo存在如下優(yōu)點(diǎn)：①M(fèi)onte carlo反演過程中采用函數(shù)約束關(guān)系尋求反演的全局最優(yōu)解，即橫波速度與層厚作為獨(dú)立變量;②在層速度突變界面，反演的橫波速度剖面與真實(shí)模型具有較小的差異，能更好地評(píng)價(jià)地下介質(zhì)的速度結(jié)構(gòu);③對(duì)初始模型依賴性較低。

為進(jìn)一步驗(yàn)證兩種方法在等厚薄層線性遞增模型中的反演質(zhì)量。圖4(b)展示了經(jīng)典阻尼最小二乘反演的理論頻散曲線與實(shí)測(cè)頻散曲線的均方根誤差，在Rayleigh波相速度100 m/s ～125 m/s區(qū)域內(nèi)存在較輕微的誤差(黃色框區(qū)域內(nèi))。圖4(d)展示了Monte carlo反演的理論頻散曲線與實(shí)測(cè)頻散曲線的均方根誤差。從圖4(b)及圖4(d)中可見，因在低頻處拾取實(shí)測(cè)頻散曲線存在一定誤差，故兩種反演算法在低頻處均無法有效擬合。最后經(jīng)對(duì)比，經(jīng)典阻尼最小二乘反演的理論頻散與實(shí)測(cè)頻散的均方根誤差為1.6%，Monte carlo反演的理論頻散與實(shí)測(cè)頻散的均方根誤差為0.8%。從擬合的殘差大小表明：在線性遞增模型中，相比與經(jīng)典阻尼最小二乘，Monte carlo反演具有較高的反演質(zhì)量。

圖5表示在第II組均勻模型中經(jīng)典阻尼最小二乘與Monte carlo 的反演結(jié)果及頻散曲線對(duì)比，為更清晰刻畫在均勻模型中的反演細(xì)節(jié)，圖5(a)給出了經(jīng)典阻尼最小二乘反演的橫波速度剖面與真實(shí)模型之間的對(duì)比，在5(a)中可見經(jīng)典阻尼最小二乘存在如下缺點(diǎn)：①經(jīng)典阻尼最小二乘為提高解的唯一性，在均勻模型中采用固定約束關(guān)系，在迭代反演中層厚度保持不變，減少了待反演的參數(shù)，只有橫波速度作為獨(dú)立變量；②在層速度突變界面，反演的橫波速度剖面與真實(shí)模型存在較大差異，未能較好地評(píng)價(jià)地下介質(zhì)速度結(jié)構(gòu)；③對(duì)初始模型依賴性較高。圖5(c)給出了Monte carlo反演結(jié)果與真實(shí)模型之間的對(duì)比，從5(c)中可知,Monte carlo存在如下優(yōu)點(diǎn)：①M(fèi)onte carlo在反演過程中采用函數(shù)約束關(guān)系尋求反演的全局最優(yōu)解，即橫波速度與層厚作為獨(dú)立變量；②在層速度突變界面，反演的橫波速度剖面與真實(shí)模型具有較小的差異，能更好的評(píng)價(jià)地下介質(zhì)的速度結(jié)構(gòu)；③對(duì)初始模型依賴性低。

為更進(jìn)一步驗(yàn)證兩種方法在均勻模型中的反演質(zhì)量。圖5(b)描繪了經(jīng)典阻尼最小二乘反演的理論頻散曲線與實(shí)測(cè)頻散曲線的均方根誤差。從圖5(b)可知，在均勻模型中的Rayleigh波不具備頻散特性，且在Rayleigh波相速度100 m/s～135 m/s區(qū)域內(nèi)存在較大的誤差(黃色框區(qū)域內(nèi))。圖5(d)描繪了Monte carlo反演的理論頻散曲線與實(shí)測(cè)頻散曲線的均方根誤差，紅色實(shí)線表示Monte carlo反演的理論頻散曲線，黃色實(shí)線表示初始模型頻散曲線，黑色實(shí)線表示實(shí)測(cè)頻散曲線。從圖5(b)及5(d)中可見，因在低頻處拾取實(shí)測(cè)頻散曲線存在一定誤差，故兩種反演算法在低頻處均無法有效擬合。最后經(jīng)對(duì)比，經(jīng)典阻尼最小二乘反演的理論頻散與實(shí)測(cè)頻散的均方根誤差為2.8%。Monte carlo反演的理論頻散與實(shí)測(cè)頻散的均方根誤差為0.9%，從擬合的殘差大小表明：在均勻模型中，相比與經(jīng)典阻尼最小二乘，Monte carlo頻散曲線反演具有較高的反演質(zhì)量。

值得注意的是，第Ⅲ組初始模型并非線性遞增模型，每層的層厚度也非等厚薄層，與第Ⅰ組初始模型有明顯區(qū)別。圖6表示在第Ⅲ組非等厚薄層遞增模型中經(jīng)典阻尼最小二乘與Monte carlo 的反演結(jié)果及頻散曲線對(duì)比。圖6(a)為經(jīng)典阻尼最小二乘反演的橫波速度剖面與真實(shí)模型之間的對(duì)比，其線條的表示關(guān)系同上。在6(a)可知，經(jīng)典阻尼最小二乘同樣存在如下缺點(diǎn)：①經(jīng)典阻尼最小二乘為提高解的唯一性，采用固定約束關(guān)系，在迭代反演中層厚度保持不變，減少了待反演的參數(shù)，只有橫波速度作為獨(dú)立變量，因此反演結(jié)果與真實(shí)模型存在較大的差異；②在層速度突變界面，由于層厚度值在反演過程中保持不變，未能較好地評(píng)價(jià)地下介質(zhì)速度結(jié)構(gòu)；③由于采用固定約束關(guān)系且反演過程需要利用導(dǎo)數(shù)信息，故初始模型的依賴性較高。圖6(c)為Monte carlo反演結(jié)果與真實(shí)模型之間的對(duì)比結(jié)果。從6(c)可知，Monte carlo存在如下優(yōu)點(diǎn)：①M(fèi)onte carlo在反演過程中采用函數(shù)約束關(guān)系尋求反演的全局最優(yōu)解，即橫波速度與層厚作為獨(dú)立變量；②在層速度突變界面，由于采用函數(shù)約束關(guān)系，能更好地評(píng)價(jià)地下介質(zhì)的速度結(jié)構(gòu)；③由于采用函數(shù)約束關(guān)系，故對(duì)初始模型的依賴性較低。

為更進(jìn)一步驗(yàn)證兩種方法在非等厚薄層遞增模型中的反演質(zhì)量。圖6(b)描繪了經(jīng)典阻尼最小二乘反演的理論頻散曲線與實(shí)測(cè)頻散曲線的均方根誤差。從圖6(b)可知，經(jīng)典阻尼最小二乘在Rayleigh波相速度100 m/s～110 m/s區(qū)域內(nèi)，存在較大的誤差(黃色框內(nèi))，因此導(dǎo)致總體誤差較大。圖6(d)描繪了Monte carlo反演的理論頻散曲線與實(shí)測(cè)頻散曲線的均方根誤差，因受低頻處拾取的實(shí)測(cè)頻散曲線在200 m/s～300 m/s區(qū)域的錯(cuò)誤影響，存在較大的誤差。但是從圖6(b)及6(d)總體可見，因在低頻處拾取實(shí)測(cè)頻散曲線存在一定誤差，故兩種反演算法在低頻處均無法有效擬合。最后經(jīng)對(duì)比，經(jīng)典阻尼最小二乘反演的理論頻散與實(shí)測(cè)頻散的均方根誤差為0.9%。Monte carlo反演的理論頻散與實(shí)測(cè)頻散的均方根誤差為0.7%，從擬合的殘差大小表明：在非等厚薄層遞增模型中，相比與經(jīng)典阻尼最小二乘，Monte carlo頻散曲線反演具有較高的反演質(zhì)量。

在同臺(tái)PC計(jì)算機(jī)上經(jīng)對(duì)比內(nèi)存消耗情況，第I組線性遞增模型中，Monte carlo消耗內(nèi)存1.09 GB，經(jīng)典阻尼最小二乘消耗內(nèi)存1.37 GB。在第II組均勻模型中，Monte carlo消耗內(nèi)存1.29 GB，經(jīng)典阻尼最小二乘消耗內(nèi)存1.35 GB，在第Ⅲ組非等厚薄層遞增模型，Monte carlo消耗內(nèi)存1.13 GB,經(jīng)典阻尼最小二乘消耗內(nèi)存1.28 GB。實(shí)驗(yàn)表明基于Monte carlo反演策略在淺部橫波速度結(jié)構(gòu)反演中具有占用內(nèi)存小的優(yōu)勢(shì)。

最后為驗(yàn)證Monte carlo反演策略的抗噪性，將實(shí)測(cè)頻散曲線分別加噪0.1、0.2，利用第Ⅲ組非等厚薄層遞增模型進(jìn)行反演，反演結(jié)果如圖7所示。

圖7 Monte carl抗噪性分析Fig.7 Monte Carlo noise resistance analysis(a)Monte carlo反演的理論頻散曲線與加噪0.1后的實(shí)測(cè)頻散曲線對(duì)比；(b)加噪0.1后Monte carlo 反演的橫波速度剖面與真實(shí)模型的對(duì)比；(c)Monte carlo反演的理論頻散曲線與加噪0.2后的實(shí)測(cè)頻散曲線對(duì)比;(d)加噪0.2后Monte carlo 反演的橫波速度剖面與真實(shí)模型的對(duì)比

值得注意的是，由于實(shí)測(cè)頻散曲線的加噪，所以反演的理論頻散曲線與加噪后的實(shí)測(cè)頻散曲線均方根誤差ε>1%，但是從圖7(b)可知，在實(shí)測(cè)頻散曲線加噪0.1后，Monte carlo 反演的橫波速度剖面與真實(shí)值十分接近，由圖7(d)可見，Monte carlo 在實(shí)測(cè)頻散曲線加噪0.2后，反演的橫波速度剖面也與真實(shí)值十分接近，從抗噪實(shí)驗(yàn)結(jié)果表明，Monte carlo反演策略具有較好的抗噪性。

2.2 二維橫波速度剖面

為精細(xì)刻畫層狀介質(zhì)模型結(jié)構(gòu)。首先，在三層層狀介質(zhì)模型中，模型參數(shù)見表3，通過Rayleigh波有限差分正演模擬[29-30]合成Rayleigh波多炮正演記錄，為采集多炮Rayleigh波合成記錄，設(shè)定觀測(cè)系統(tǒng)為：共采集70炮數(shù)據(jù)，每炮移動(dòng)4 m，第一炮震源位置為(0，0)m處，由于勘探目的層的深度應(yīng)大于40 m，本算例共設(shè)置84道接收，接收點(diǎn)深度為0 m，第一個(gè)檢波器的位置為58 m，道間距為2 m，因此最小偏移距為58 m，最大偏移距為224 m， 排列長(zhǎng)度166 m，縱橫向網(wǎng)格大小為501×451，縱橫向網(wǎng)格間距都為1 m，因此模型尺寸為500 m×450 m，邊界層厚度為11個(gè)網(wǎng)格點(diǎn)，采樣時(shí)間間隔為0.000 1 ms，采樣點(diǎn)數(shù)為2001，因此總時(shí)間為2 s，震源采用主頻為9 Hz的雷克子波。自由邊界采取應(yīng)力鏡像法[28]，圖8為多炮面波采集示意圖。將采集的多炮Rayleigh波合成記錄，通過相移法提取多條Rayleigh波頻散曲線。最后基于Monte carlo反演策略，對(duì)多炮Rayleigh合成記錄中提取出的多條頻散曲線進(jìn)行自動(dòng)反演。

表3 三層模型參數(shù)表

圖8 多炮面波采集示意圖Fig.8 Schematic diagram of multi-shot surface wave collection

將采集到的多炮合成記錄根據(jù)相移法提取頻散曲線，以圖9(a)所示的炮記錄作為輸入，計(jì)算得到如圖9(b)所示的頻散能量圖，沿能量峰值線拾取實(shí)測(cè)頻散曲線得到圖9(c)。每一炮地震記錄均可產(chǎn)生頻散能量圖，可沿能量峰值線拾取每一炮地震記錄的實(shí)測(cè)頻散曲線。

圖9 第一炮Rayleigh波合成記錄及實(shí)測(cè)頻散數(shù)據(jù)Fig.9 Rayleigh-wave synthesis record and measured dispersion data of the first shot(a)垂直位移分量的地震記錄；(b)頻散能量圖；(c)實(shí)測(cè)頻散曲線

根據(jù)拾取出的70條實(shí)測(cè)頻散曲線，反演其相速度，可得到各70炮炮集排列中點(diǎn)處所對(duì)應(yīng)的橫波速度隨深度的變化曲線，根據(jù)一系列變化曲線可得到二維橫波速度剖面。為驗(yàn)證Monte carlo在70炮中的反演質(zhì)量，現(xiàn)從70炮合成記錄中抽取第25炮合成記錄，根據(jù)相移法拾取實(shí)測(cè)頻散曲線，最后利用表4所示的初始模型參數(shù)值進(jìn)行反演，迭代次數(shù)共1 000 次，為使Monte Carlo反演法能在眾多局部最優(yōu)解中尋求全局最優(yōu)解，因此設(shè)置橫波速度與層厚的空間約束參數(shù)、分為0.37、0.38。

表4初始模型參數(shù)值

表 5 第25炮反演結(jié)果統(tǒng)計(jì)表

在空間約束參數(shù)、所確定的隨機(jī)樣本空間中，每一個(gè)隨機(jī)樣本空間對(duì)應(yīng)每一個(gè)模型參數(shù)值，將每一個(gè)模型參數(shù)值正演一次，并與實(shí)測(cè)頻散曲線計(jì)算均方根誤差，反演結(jié)果如圖10所示，圖10(a)為Monte carlo反演的理論頻散與實(shí)測(cè)頻散的均方根誤差，最小均方根誤差為0.6%。圖10(b)表示Monte carlo反演的橫波速度剖面與真實(shí)模型的對(duì)比結(jié)果，最小均方根誤差為0.6%。表5為均方根誤差達(dá)0.6%時(shí)的反演結(jié)果。

圖10 排列中點(diǎn)為141 m處的反演結(jié)果Fig.10 Inversion results of the permutation midpoint at 141 m(a)Monte carlo反演的理論頻散與實(shí)測(cè)頻散的對(duì)比；(b)Monte carlo反演的橫波速度剖面與真實(shí)模型的對(duì)比

在第25炮頻散曲線反演中耗時(shí)共計(jì)79.846 s，誤差0.6%，表明在迭代1 000次中，Monte carlo反演法收斂速度快，誤差較小，具有較好的反演質(zhì)量。

提取70炮的實(shí)測(cè)頻散曲線，根據(jù)表4初始模型參數(shù)值利用Monte carlo自動(dòng)反演，可得到二維橫波速度剖面。圖11表示Monte carlo在提取的70條頻散曲線中自動(dòng)反演的結(jié)果并與真實(shí)模型之間的對(duì)比，從圖11反演對(duì)比結(jié)果上來看，Monte carlo反演剖面與真實(shí)模型基本吻合，表明基于Monte carlo反演策略在淺地表橫波速度反演中具有良好的反演質(zhì)量，能較好刻畫模型的精細(xì)結(jié)構(gòu)。

圖11 Monte carlo反演的橫波速度剖面與真實(shí)模型對(duì)比結(jié)果Fig.11 Shear wave velocity profile of Monte Carlo inversion and comparison with the real model(a)Monte carlo反演結(jié)果；(b)真實(shí)模型

2.3 Monte carlo實(shí)際應(yīng)用案例

測(cè)試點(diǎn)位于貴州大學(xué)西校區(qū)圖書館前，如圖12中紅色實(shí)線所示。為在測(cè)試點(diǎn)構(gòu)建地下二維垂直橫波速度場(chǎng)，采用工程地震儀WZG-24A進(jìn)行Rayleigh波數(shù)據(jù)采集，炮激發(fā)源采用重5 kg的大錘敲擊0.2 m的鋁板作為沖擊激發(fā)。共采集17炮數(shù)據(jù)，第一炮震源位置為(0,0)m，采用12道檢波器線性排列接收信號(hào)，檢波器深度為0.1 m，最小偏移距為6 m，最大偏移距為28 m，道間距為2 m，因此排列長(zhǎng)度為22 m，沿測(cè)線方向以2 m的間隔不斷移動(dòng)震源及檢波器位置。根據(jù)MASW原理，在該觀測(cè)系統(tǒng)下可測(cè)定17 m～49 m范圍內(nèi)的地下介質(zhì)情況。

圖12 現(xiàn)場(chǎng)測(cè)試示意圖Fig.12 Field test diagram

沿測(cè)線方向并根據(jù)多炮面波數(shù)據(jù)采集，共采集17炮數(shù)據(jù)，圖13所示為該測(cè)試點(diǎn)沿測(cè)線采集的第6炮及第12炮的炮集記錄。從圖13的炮集記錄中可見Rayleigh波信噪比最高，能量最強(qiáng)，占據(jù)主導(dǎo)地位。

圖13 測(cè)試點(diǎn)獲取的面波記錄Fig.13 Typical surface wave records obtained at the test site(a)第六炮炮集記錄;(b)第12炮炮集記錄

進(jìn)一步，根據(jù)第6炮炮集記錄及第12炮炮集記錄，利用相移法提取Rayleigh波實(shí)測(cè)頻散曲線，根據(jù)實(shí)測(cè)頻散曲線的形態(tài)特征初步估計(jì)初始模型。初始模型參數(shù)估計(jì)如表6所示。

表6 初始模型參數(shù)Tab.6 Initial model parameter values

最后根據(jù)初始模型參數(shù)值，采取Monte carlo反演策略反演其相速度，可得到該排列中點(diǎn)所對(duì)應(yīng)的橫波速度剖面(圖14)。圖14(a)為從第6炮炮數(shù)據(jù)中提取頻散曲線并反演其相速度，所得到的排列中點(diǎn)為27 m的橫波速度剖面，圖14(b)為從第12炮炮集記錄中提取頻散曲線并反演其相速度，所得到的排列中點(diǎn)為39 m的橫波速度剖面。

圖14 Monte carlo 反演的橫波速度剖面Fig.14 Shear wave velocity profile by Monte Carlo inversion(a)排列中點(diǎn)為27 m的反演結(jié)果;(b)排列中點(diǎn)為39 m的反演結(jié)果

根據(jù)一系列反演的橫波速度剖面并排列，可得到該測(cè)試點(diǎn)所測(cè)試的17 m～49 m范圍內(nèi)的二維橫波速度剖面，如圖15所示。

圖15 測(cè)試點(diǎn)的二維橫波速度剖面Fig.15 2D shear wave velocity profile obtained at the test site

圖15中，兩條黃色虛線分別代表圖14(a)中27 m處反演的橫波速度剖面，圖14(b)中39 m處反演的橫波速度剖面。從圖15中總體可見，該地下介質(zhì)在深度為0 m～1.5 m處，橫波速度值較小，從現(xiàn)場(chǎng)的推斷上來看上層覆蓋物應(yīng)為第四紀(jì)松散的土質(zhì)層。而深度為1.5 m～2.8 m，橫波速度值較大，應(yīng)推斷該處為巖石層。深度為2.8 m～9 m處的橫波速度值總體較小，5.5 m～6 m處局部橫波速度又較大，故推斷該深度下應(yīng)為黏土層，局部夾厚度約為0.5 m的巖石層。

4 結(jié)論

筆者在求解水平層狀介質(zhì)中Rayleigh波的本征值問題之上發(fā)展了Monte carlo局部線性化反演算法，在典型數(shù)值模型中，通過三組初始模型參數(shù)值分別用基于經(jīng)典阻尼最小二乘與基于Monte carlo兩種反演策略進(jìn)行頻散曲線反演，并用反演的理論頻散與實(shí)測(cè)頻散進(jìn)行均方根誤差對(duì)比分析，可得到如下認(rèn)識(shí)。

1)基于經(jīng)典阻尼最小二乘反演策略在反演過程中采用固定約束，未能較好地評(píng)價(jià)地下介質(zhì)速度結(jié)構(gòu)，且在反演過程中需利用導(dǎo)數(shù)信息，所以對(duì)初始模型依賴性較高，而基于Monte carlo反演策略采用函數(shù)約束關(guān)系，能較好地評(píng)價(jià)地下介質(zhì)速度結(jié)構(gòu)，且不需要利用導(dǎo)數(shù)信息，故對(duì)初始模型依賴性低，在反演過程中也較為穩(wěn)定，并具有較高的反演精度。

2)Monte carlo反演策略在實(shí)測(cè)頻散曲線加噪0.1、0.2后仍然具有較高的反演精度，表明Monte carlo 反演策略具有較好的抗噪性。

3)在同臺(tái)PC計(jì)算機(jī)條件下，相比于經(jīng)典的阻尼最小二乘，Monte carlo反演策略具有較快的收斂速度和反演質(zhì)量。

4)在多條Rayleigh波頻散曲線自動(dòng)反演中，基于Monte carlo反演的二維橫波速度剖面與真實(shí)模型也較為吻合，能較好地刻畫地質(zhì)模型的精細(xì)結(jié)構(gòu)。

雖然基于Monte carlo頻散曲線反演具有較高的計(jì)算效率，但是在低速半空間、高速表層、橫向不均勻介質(zhì)模型等復(fù)雜模型中該算法存在適用性較低等缺點(diǎn)。因此在今后的研究中，將Monte Carlo反演算法應(yīng)用于低速半空間、高速表層、橫向不均勻介質(zhì)等復(fù)雜模型的橫波速度建模中是下一步的研究重點(diǎn)，特別是針對(duì)Monte Carlo反演算法存在分辨率低，在未知層厚的情況下，適用性變差等缺陷，有待進(jìn)一步研究。