基于Monte carlo法的Rayleigh波淺層橫波速度結構反演方法研究

胡昌榮, 黃金強,2, 劉 宏,2, 竇勇碩, 何云川

(1.貴州大學 喀斯特地質資源與環境教育部重點實驗室，貴陽 550025；2.貴州大學 資源與環境工程學院，貴陽 550025)

0 引言

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

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

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

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

1 方法原理

1.1 Monte carlo反演算法

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

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

Vs={Vs1,Vs2,…,Vsn}

h={h1,h2,…,hn}

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

Vp={Vp1,Vp2,…,Vpn}

(1)

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

(2)

其中：c為Rayleigh波相速度;k為波數;n為層數;i為給定波數所對應的Rayleigh波相速度的點數。

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

(3)

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

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

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

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

隨后，在一定空間范圍內隨機生成每一次正演的樣本模型:

Vsj=Vs+Δj

hl=h+Δl

(4)

為獲取相應地層中的Vp值，根據Everett[23]提出的縱橫波速度比與泊松比的關系進行縱波速度值推算：

(5)

通過密度與縱波速度經驗公式對相應地層中的密度值進行計算[24]：

(6)

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

1.2 阻尼最小二乘反演算法

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

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

(7)

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

(8)

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

(9)

表述的矩陣方程式(9)可進一步簡述為：

Δc=AΔx

(10)

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

針對線性最小二乘反演問題，可將目標函數表述為若干個函數的平方和形式:

(11)

式中：ciobs表示實測Rayleigh波相速度值，cical表示理論Rayleigh波相速度值。

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

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

(12)

式中：‖‖2表示l2范數;W表示權重矩陣;λ為阻尼因子。

進一步將目標函數簡化并整理得：

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

(13)

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

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

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

(14)

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

A=UΛVT

(15)

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

(16)

即有：

(17)

進一步整理得：

(18)

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

(19)

(20)

1.3 Monte carlo 檢驗準則

采用均方根誤差定義每一條頻散曲線與實測頻散曲線的誤差，將均方根誤差[26]作為樣本模型的檢驗標準，若誤差達到最小時，則接收該樣本模型是最優模型，反之，則排斥該樣本模型。具體反演算法流程如圖2所示。

(21)

2 模型測試/數值模擬

2.1 三層模型實驗

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

三層模型參數如表1所示，由表1中可知，橫波速度較小，泊松比高達0.47，在工程物探中主要代表未固結的松散碎屑物質(如黏土層、飽和沙土等)。在高泊松比地質結構中，由于傳統的PML邊界條件在Rayleigh波正演數值模擬中將產生數值不穩定?；诖耍捎枚噍S卷積完全匹配層吸收邊界條件處理人工截斷邊界[27]，應力鏡像+二階速度展開處理自由表面邊界[28]，數值模擬方法為交錯網格有限差分正演模擬[29-30]?？v橫向網格大小為2001×601，縱橫向網格間距均為0.1 m，實際模擬區域為200 m×60 m，邊界層厚度為801個網格點，時間采樣間隔為0.05 ms，采樣點數為60 001，總模擬時間為3 s，震源為主頻8 Hz的雷克子波，震源位置位于(0，0)m，由于勘探目的層的深度應大于25 m，本算例共設置61道接收，接收點深度為0 m，第一個檢波器的位置為70 m，道間距為2 m，因此最小偏移距為70 m，最大偏移距為190 m。

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

表 1 三層模型參數表

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

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

表2 初始模型參數值

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

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

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

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

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

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

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

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

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

在同臺PC計算機上經對比內存消耗情況，第I組線性遞增模型中，Monte carlo消耗內存1.09 GB，經典阻尼最小二乘消耗內存1.37 GB。在第II組均勻模型中，Monte carlo消耗內存1.29 GB，經典阻尼最小二乘消耗內存1.35 GB，在第Ⅲ組非等厚薄層遞增模型，Monte carlo消耗內存1.13 GB,經典阻尼最小二乘消耗內存1.28 GB。實驗表明基于Monte carlo反演策略在淺部橫波速度結構反演中具有占用內存小的優勢。

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

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

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

2.2 二維橫波速度剖面

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

表3 三層模型參數表

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

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

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

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

表4初始模型參數值

表 5 第25炮反演結果統計表

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

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

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

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

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

2.3 Monte carlo實際應用案例

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

圖12 現場測試示意圖Fig.12 Field test diagram

沿測線方向并根據多炮面波數據采集，共采集17炮數據，圖13所示為該測試點沿測線采集的第6炮及第12炮的炮集記錄。從圖13的炮集記錄中可見Rayleigh波信噪比最高，能量最強，占據主導地位。

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

進一步，根據第6炮炮集記錄及第12炮炮集記錄，利用相移法提取Rayleigh波實測頻散曲線，根據實測頻散曲線的形態特征初步估計初始模型。初始模型參數估計如表6所示。

表6 初始模型參數Tab.6 Initial model parameter values

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

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

根據一系列反演的橫波速度剖面并排列，可得到該測試點所測試的17 m～49 m范圍內的二維橫波速度剖面，如圖15所示。

圖15 測試點的二維橫波速度剖面Fig.15 2D shear wave velocity profile obtained at the test site

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

4 結論

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

1)基于經典阻尼最小二乘反演策略在反演過程中采用固定約束，未能較好地評價地下介質速度結構，且在反演過程中需利用導數信息，所以對初始模型依賴性較高，而基于Monte carlo反演策略采用函數約束關系，能較好地評價地下介質速度結構，且不需要利用導數信息，故對初始模型依賴性低，在反演過程中也較為穩定，并具有較高的反演精度。

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

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

4)在多條Rayleigh波頻散曲線自動反演中，基于Monte carlo反演的二維橫波速度剖面與真實模型也較為吻合，能較好地刻畫地質模型的精細結構。

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