瑞雷波多階模式頻散曲線稀疏正則化反演方法研究

崔巖,王彥飛

1 山東科技大學地球科學與工程學院，青島 266590 2 河北地質(zhì)大學，京津冀城市群地下空間智能探測與裝備重點實驗室，石家莊 050031 3 中國科學院地質(zhì)與地球物理研究所,中國科學院油氣資源研究重點實驗室，北京 100029 4 中國科學院大學,北京 100049 5 中國科學院地球科學研究院，北京 100029

0 引言

近年來，與環(huán)境、工程、能源密切相關的近地表地球物理探測技術及方法成為研究熱點，橫波速度是近地表最關鍵的地震參數(shù)之一，目前主要用面波來提取橫波速度.面波速度建模有三個步驟：(1)數(shù)據(jù)的采集和處理；(2)頻散曲線拾取；(3)頻散曲線反演.頻散曲線反演由頻散曲線獲取地下橫波速度，是其中的關鍵步驟.目前人們最常使用的面波是瑞雷波，它有很強的刻畫淺地表/淺水介質(zhì)(沉積物)的能力(夏江海，2015；Wang et al.,2020).瑞雷波在層狀介質(zhì)中發(fā)生頻散，并且會發(fā)育多階模態(tài)，即一個頻率對應多個不同的相速度值(牛濱華和孫春巖，2002).對某給定頻率，它所對應的最小相速度值稱為該頻率的基階模式相速度(一階模式相速度)，其余相速度值可統(tǒng)稱為該頻率的高階模式相速度.高階模式相速度按照從低到高的順序，又分別稱為一階高階模式相速度(二階模式相速度)、二階高階模式相速度(三階模式相速度)等.在對瑞雷波頻散曲線進行反演時，既可以使用單階模式，也可以使用多階模式.一般情況下，地層剛度隨深度逐層增加，此時瑞雷波能量大部分集中在基階模式，使用基階模式就能夠取得好的效果.但由于實際地層結(jié)構(gòu)是復雜多變的，瑞雷波能量并不總是集中在基階模式，這就使得僅用基階模式頻散曲線反演出的橫波速度偏離真實值.我們可以從數(shù)學角度把這看作一個加權(quán)問題，各階模式的加權(quán)系數(shù)是其能量占總能量的百分比，它決定了各階模式頻散曲線在反演橫波速度模型時影響的大小.前人通過分析指出，在某些特殊地質(zhì)結(jié)構(gòu)中(例如含低速夾層)，高階模式較基階模式在某些頻段可能更加發(fā)育，僅用基階模式頻散曲線反演出的橫波速度存在失真(徐佩芬等，2020).而多階模式頻散曲線聯(lián)合反演不僅能使反演過程穩(wěn)定、提高橫波速度的反演精度和計算結(jié)果的可靠性，還能夠提高模型縱向分辨率(宋先海，2008；潘冬明等，2010；夏江海等，2015).在實際工作中，頻散曲線提取的誤差會隨著階數(shù)的增加而增大，為了保證結(jié)果的可靠性，本文研究的瑞雷波多階模式頻散曲線反演使用基階、一階高階和二階高階模式頻散曲線.

瑞雷波多階模式頻散曲線反演是一個典型的地球物理反問題，任何反問題的研究都涉及正演模擬、反演建模和數(shù)值實現(xiàn)三個方面(王彥飛等，2011).利用頻散曲線反演地下速度結(jié)構(gòu)，頻散曲線的計算是首先要解決的正問題，對此前人已進行了大量研究，研究基于水平層狀介質(zhì)模型，主要算法有Haskell算法(Haskell,1953)、Schwab-Knopoff算法(Schwab and Knopoff,1970)、廣義反射-透射系數(shù)算法(Chen,1993,1999)等.其中廣義反射-透射系數(shù)算法在高頻更為穩(wěn)定、精確，在求取頻散曲線的同時還能給出瑞雷波本征位移、應力及理論地震圖的求解公式.在反演模型的建立上，現(xiàn)有文獻大多僅考慮數(shù)據(jù)的擬合，將頻散曲線理論計算值與實測值之間誤差的L2范數(shù)作為目標函數(shù)(蔡偉等，2017).但從連續(xù)介質(zhì)物理學的概念出發(fā)，地球物理模型最好采用“逐塊”光滑的間斷函數(shù)集來描述，允許空間存在間斷面，間斷面之間的介質(zhì)參數(shù)是連續(xù)且光滑的(楊文采，1997)，這種物理上的間斷可以通過某些模型參數(shù)在數(shù)學上的稀疏性進行描述.L1范數(shù)正則化將待反演參數(shù)的L1范數(shù)或其稀疏表示后的L1范數(shù)最小作為約束條件，目前已應用于反射系數(shù)反演(Wang et al.，2016)、波阻抗反演(Liu et al.，2018)、疊前彈性參數(shù)反演(劉曉晶等，2016)和基階模式頻散曲線反演(Cui and Wang,2021)等，研究表明L1范數(shù)正則化可以提高分辨率，增強對非高斯噪聲的魯棒性(王彥飛等，2021).在反問題的數(shù)值實現(xiàn)上，現(xiàn)有方法可以分為線性和非線性兩大類(Wang et al.，2011b)，在瑞雷波頻散曲線反演中使用的線性方法主要有最小二乘法(Gabriels et al.，1987)、奇異值分解和Levenberg-Marquardt法(Xia et al.，1999)、Occam算法(Lai et al.，2002)等；非線性方法有遺傳算法(Yamanaka and Ishida，1996)、小波分析(Tillmann，2005)、粒子群算法(Song et al.，2012)，洗牌蛙跳算法(Sun et al.，2017)、蜂群算法(于東凱等，2018)等.線性方法結(jié)構(gòu)簡單、計算效率高、所需內(nèi)存少，但對初始模型依賴性較強，容易陷入局部極值；非線性方法對初始模型的要求不高，具有較強的全局尋優(yōu)能力，但收斂速度較慢，運算量大，目前線性方法的應用更為廣泛.

本文針對目前瑞雷波多階模式頻散曲線反演中僅考慮數(shù)據(jù)的擬合，缺乏對模型約束的問題，研究了瑞雷波多階模式頻散曲線的稀疏正則化反演方法.首先討論瑞雷波多階模式頻散曲線正演的原理及數(shù)值計算，原理采用經(jīng)典的廣義反射-透射系數(shù)法，數(shù)值計算上采用具有三階收斂速度的快速求根方法.然后是反演模型的建立及反問題的數(shù)值實現(xiàn)，通過在目標函數(shù)中引入模型的L1范數(shù)施加稀疏正則化約束，使反演結(jié)果更加符合地質(zhì)實際，并針對該模型提出一種隱式迭代正則化算法.最后將新的反演方案應用于理論模型，并與傳統(tǒng)方法進行對比.

1 正問題

1.1 瑞雷波頻散曲線正演原理

所有反演方法都涉及理論數(shù)據(jù)的模擬，利用頻散曲線反演地下速度結(jié)構(gòu)，頻散曲線的計算是首先要解決的正問題，本文選用廣義反射-透射系數(shù)法.簡正振型(Seismic Normal Modes)是無源彈性動力學方程在特定邊界條件下的非零解，從這個基本原理出發(fā)，計算任意給定頻率和任意固體介質(zhì)參數(shù)下水平層狀介質(zhì)中的簡正振型，從而得到頻散曲線.

圖1 廣義反射-透射系數(shù)法示意圖Fig.1 Schematic diagram of generalized reflection-transmission coefficient method

(1)

當j=1,2,…,N時

(2a)

當j=0時(自由界面處)

(2b)

(3)

由此可得

(4)

線性齊次代數(shù)方程組存在非零解的充要條件是系數(shù)行列式為零，因此

上式左端是瑞雷波頻散函數(shù)(久期函數(shù))，它是一個關于模型參數(shù)(拉梅常數(shù)、密度、層厚)、頻率和瑞雷波相速度的函數(shù).給定模型參數(shù)，選取某一頻率，只有有限個相速度能夠滿足方程(5)，它們就是頻散函數(shù)的根.這一組根由小到大分別對應于該頻率成分的基階模式相速度(一階模式相速度)、一階高階模式相速度(二階模式相速度)、二階高階模式相速度(三階模式相速度)，以此類推，所有頻率成分和它們所對應的某模式相速度之間的關系就是該模式的頻散曲線.

1.2 瑞雷波多階模式頻散曲線的數(shù)值計算

頻散函數(shù)的快速求解直接影響反演過程的穩(wěn)定性和結(jié)果的準確性.傳統(tǒng)的二等分方法從某個初始值(通常是第一層S波速度的0.8倍)開始，并沿相速度方向以一定步長執(zhí)行根搜索.當滿足給定的誤差精度時，認為獲得了該頻率下的相速度.然后以相同的方式滑動頻率以獲取整個頻散曲線.該方法不具有超線性收斂性，計算效率低.當復雜的地層模型引起“模式接吻”或“模式跳變”時，可能會發(fā)生嚴重的缺陷，例如丟根和模式誤判(夏江海，2015).因此有必要在正演模擬中引入先進的尋根方法，以提高計算的效率和準確性，故本文考慮快速根搜索方法，該方法具有三階收斂速度(Wang and Xiao，2001).

對一個地質(zhì)模型，給定頻率ω，瑞雷波相速度VR要滿足方程(5).為了簡化公式，將方程(5)的左側(cè)定義為非線性函數(shù)φ(x)，其中x指相速度VR.方程(5)轉(zhuǎn)換為尋找非線性方程φ(x)=0的根.眾所周知，φ(x)相對于x是無限可微的，在項(x-xn)3處截斷泰勒展開式，可以得到

(6)

其中ξn的取值范圍在x和xn之間.

通過式(6)可以獲得迭代公式

xn+1=xn-

(7)

上述頻散曲線求根的迭代過程可以通過圖2所示的頻散曲線計算流程實現(xiàn).

圖2 瑞雷波多階模式頻散曲線計算流程圖Fig.2 Flowchart of Rayleigh wave multiple mode dispersion curve extraction

通過上述算法得出Rayleigh波多階模式頻散曲線，該曲線作為理論模擬數(shù)據(jù)用于反演地下橫波速度.

2 反問題

2.1 瑞雷波多階模式頻散曲線反演建模

通過上面的討論可以看出，瑞雷波頻散曲線與模型參數(shù)有關.正演是根據(jù)模型計算頻散曲線，反演則是根據(jù)頻散曲線反推模型.要討論瑞雷波多階模式頻散曲線聯(lián)合反演的建模問題，首先需要對單階模式頻散曲線反演建模.

2.1.1 單階模式

根據(jù)式(5)，對于瑞雷波單階模式的頻散曲線，頻率成分ωi和它所對應的相速度VRi之間的關系可表示為

Fi(ωi,VRi,λ,μ,ρ,h)=0i=1,2,…,M

(8)

其中，M為頻散曲線數(shù)據(jù)點的個數(shù)，ωi為頻率，VRi為頻率ωi對應的瑞雷波相速度，λ=[λ(1),λ(2),…,λ(N+1)]、μ=[μ(1),μ(2),…,μ(N+1)]、ρ=[ρ(1),ρ(2),…,ρ(N+1)]、h=[h(1),h(2),…,h(N)]分別為每一層的拉梅常數(shù)、密度、層厚組成的向量，它們共同描述了該水平層狀模型.

Fi(ωi,VRi,VP,VS,ρ,h)=0i=1,2,…,M

(9)

由上式可知，橫波速度、縱波速度、密度和層厚都會對瑞雷波頻散曲線產(chǎn)生影響.但如果同時反演橫波速度、縱波速度、密度和層厚，反演過程將極不穩(wěn)定，存在嚴重的多解性.為此，我們希望減少方程中的未知量，專注于主要參數(shù).

研究表明，橫波速度是瑞雷波頻散曲線最主要的影響因素.因此在反演過程中可以假定密度和縱波速度已知(夏江海，2015)，將層厚固定為合理的薄層(常數(shù)).這樣只有橫波速度是未知的，式(9)可以簡化為

Fi(ωi,VRi,VS)=0,i=1,2,…,M

(10)

式中，VRi為該頻率對應的瑞雷波某階模式相速度，VS為橫波速度向量.

2.1.2 多階模式

本文研究基階模式、一階高階模式和二階高階模式頻散曲線的聯(lián)合反演.由于高階模式是在頻率大于一定數(shù)值(截止頻率)后才逐漸出現(xiàn)的，因此式(10)改寫為

Fi0(ωi0,VR0,i0,VS)=0 (i0=1,2,…,M0),

Fi1(ωi1,VR1,i1,VS)=0 (i1=1,2,…,M1),

Fi2(ωi2,VR2,i2,VS)=0 (i2=1,2,…,M2),

(11)

2.2 非線性反演問題的正則化

2.2.1 Tikhonov正則化

廣義Tikhonov正則化是解下面的最優(yōu)化問題

(12)

方程右側(cè)第一項表示數(shù)據(jù)的擬合，第二項Ω(·)是穩(wěn)定因子，α>0是正則化因子，m是模型參數(shù)，d是觀測數(shù)據(jù)，F(xiàn)(m)是由模型到數(shù)據(jù)的正演算子.

2.2.2 稀疏約束正則化

選擇不同的穩(wěn)定因子Ω(·)會導致不同的正則化方案，考慮解的稀疏約束，這里的“稀疏”一詞用于描述具有相對較少非零值的時間序列，在其他量為常數(shù)的情況下，稀疏度越大，反演結(jié)果中的層數(shù)越少.由于反演的橫波速度模型是一維的，因此考慮具有一維稀疏性的L1范數(shù)正則化.L1范數(shù)正則化將待反演參數(shù)的L1范數(shù)或其稀疏表示后的L1范數(shù)最小作為約束條件，可以提高分辨率，增強對非高斯噪聲的魯棒性(王彥飛等,2011；Wang et al.，2011a).此處選取離散模型的l1范數(shù)作為Ω(·)，即

=J1(m)+J2(m).

(13)

依據(jù)式(13)，我們構(gòu)造瑞雷波基階、一階高階和二階高階模式頻散曲線聯(lián)合反演的目標函數(shù)

(14)

注：(1)瑞雷波多階模式頻散曲線反演在方法上與基階模式頻散曲線反演區(qū)別不大，只需依據(jù)數(shù)據(jù)精度給予高階模式數(shù)據(jù)相同或不同的權(quán)重，然后將其作為數(shù)據(jù)集加入反演過程(Xia et al.，2003).在數(shù)值實驗中，數(shù)據(jù)通常由正演模擬產(chǎn)生，不同模式頻散曲線數(shù)據(jù)的精度一致，故選擇相同權(quán)重具有合理性；而在實際資料中，由于數(shù)據(jù)的采集、處理和頻散曲線提取都會對數(shù)據(jù)精度產(chǎn)生影響，不同模式頻散曲線數(shù)據(jù)的精度可能存在較大差異，故需要根據(jù)具體情況合理設置權(quán)重.(2)正則化因子影響反演結(jié)果的穩(wěn)定性，根據(jù)反問題的正則化理論，正則化因子應當是噪聲水平δ的等價無窮小(王彥飛，2007).在數(shù)值實驗中，噪聲水平δ是已知的，當信噪比遠遠大于1時，可以直接通過δ2-ε(ε是一個小的正數(shù))計算正則化因子；在實際資料中，噪聲水平通常不確定，但可依據(jù)上述參數(shù)選擇理論對正則化因子進行幾何級數(shù)選取.

2.2.3 隱式迭代正則化算法

編程最簡單，最容易實現(xiàn)的梯度方法是最速下降法

mk+1=mk+ωksk,(15)

在第k步，梯度可由下式計算

(16)

mk+1=mk+Δm,(17)

其中Δm是更新項.

Fl(m)的線性化格式可以寫作

Fl(mk+1)=Fl(mk+Δm)=Fl(mk)+HlkΔm

+O‖(Δm)2‖,(18)

(19)

J2可近似為

(20)

(21)

其中γ(·)是向量函數(shù)，表示為

(22)

因此目標函數(shù)Jα(mk+1)的梯度可以近似為

+αγ(Δmp),(23)

從而得到隱式迭代公式

(24)

即

(25)

給定Δmp的初始猜測值，下一步迭代將滿足

(26)

這個線性系統(tǒng)可以通過共軛梯度法求解，且‖Δmp+1-Δmp‖≤δ時迭代終止，其中δ>0是一個很小的正數(shù).

獲得最佳解Δmp后，可以通過以下方式更新外部迭代點

mk+1=mk+Δmp.

(27)

由于極小化問題的凸性，采用上述共軛梯度算法可以收斂到極小化問題的解.上述迭代求解橫波速度的方法定義為瑞雷波多階模式頻散曲線稀疏正則化反演算法(SIRWDC，Sparse inversion of Rayleigh wave dispersion curve)，算法描述如下：

第一步，給定最大外部迭代步數(shù)kmax，最大內(nèi)部迭代步數(shù)pmax，Δmp的初始猜測值，終止誤差δ，p=0，k=0；

第二步，使用共軛梯度法根據(jù)式(26)進行迭代，直到‖Δmp+1-Δmp‖≤δ或p=pmax，得到Δmp；

第三步，根據(jù)式(27)更新模型；

第四步，如果‖mk+1-mk‖≤δ，迭代終止；否則令k:=k+1，回到第二步.

3 數(shù)值實驗

為了測試新算法的有效性和可靠性，設計半空間上覆四套地層的理論模型，模型參數(shù)如表1所示.基于廣義反射-透射系數(shù)法對上述模型的瑞雷波多階模式頻散曲線(基階、一階高階和二階高階)進行正演模擬，總頻帶范圍取實際生產(chǎn)中常用的1～100 Hz，在數(shù)值計算時采用快速根搜索方法，結(jié)果如圖3所示.

表1 地層模型參數(shù)Table 1 Layer model′s parameters

圖3 正演模擬得到的瑞雷波多階模式(基階、一階高階和二階高階)頻散曲線Fig.3 Forward simulation of Rayleigh wave multiple mode dispersion curves (fundamental mode,first-order mode,second-order mode)

3.1 理想數(shù)據(jù)

將正演模擬得到的多階模式頻散曲線作為觀測數(shù)據(jù)，分別采用SIRWDC方法和傳統(tǒng)方法反演橫波速度模型，3次迭代后的結(jié)果如圖4所示.從圖中可以看出，相較于傳統(tǒng)方法，SIRWDC方法對速度變化的刻畫更為精細，反演得到的速度模型更逼近真實.但由于反問題存在多解性，在缺乏先驗信息(例如地質(zhì)認識)的情況下，反演結(jié)果不會隨著迭代次數(shù)的增加而得到顯著改善.

圖4 SIRWDC與傳統(tǒng)方法反演結(jié)果的對比(理想數(shù)據(jù))Fig.4 Comparison of SIRWDC with traditional algorithm

3.2 含噪聲的數(shù)據(jù)

理想情況下，頻散分析應該完整地識別和提取每個模式的頻散曲線.而實際中，利用頻散能量成像提取頻散曲線會遇到各種挑戰(zhàn)，因此我們對理想頻散曲線數(shù)據(jù)進行刪減和加噪，使其更加接近真實情況.令提取的基階頻散曲線的頻率范圍在5～50 Hz，一階高階頻散曲線的頻率范圍在40～100 Hz，二階高階頻散曲線的頻率范圍在70～100 Hz，并且在頻率范圍10～30 Hz處添加噪聲，噪聲用來模擬基階模式數(shù)據(jù)被高階模式波或體波污染的情況，結(jié)果如圖5所示.

圖5 含噪聲的瑞雷波多階模式(基階、一階高階和二階高階)頻散曲線Fig.5 Multiple mode dispersion curves in noisy case (fundamental mode,first-order mode and second-order mode)

圖6 SIRWDC與傳統(tǒng)方法反演結(jié)果的對比(含噪聲數(shù)據(jù))Fig.6 Comparison of SIRWDC with the traditional method (noisy case)

4 結(jié)論

本文針對目前瑞雷波多階模式頻散曲線反演中僅考慮數(shù)據(jù)的擬合，缺乏對模型的約束，不能很好地刻畫地層間斷面的問題，研究了瑞雷波多階模式頻散曲線稀疏正則化反演方法，從正演模擬、反演建模和數(shù)值實現(xiàn)三個方面展開工作，提出了一套完整的瑞雷波多階模式頻散曲線反演方案：

(1)頻散曲線正演基于經(jīng)典的廣義反射-透射系數(shù)法，采用一種快速求根方法進行瑞雷波相速度數(shù)值計算，與文獻中常用的二分類方法相比，能夠在很短的時間內(nèi)達到最優(yōu)的收斂效果.

(2)反演建模時考慮了地球物理模型的“逐塊”光滑特性，由于這種物理上的間斷可以通過模型參數(shù)在數(shù)學上的稀疏性進行描述，通過在目標函數(shù)中引入模型的L1范數(shù)對模型進行稀疏性刻畫，使反演出的橫波速度模型更加符合地質(zhì)實際，并且能夠增強反演算法對非高斯噪聲的魯棒性.

(3)在反問題的數(shù)值實現(xiàn)上，針對上述稀疏正化模型提出一種隱式迭代正則化算法，其迭代算子具有非膨脹特性，可以收斂到極小化問題的解，是一種結(jié)構(gòu)簡單、計算效率較高的反演計算方法.

數(shù)值實驗結(jié)果表明，新的反演方案具有計算效率高，模型“逐塊”光滑的特性刻畫好、對非高斯噪聲魯棒性強的特點.在下一步的研究工作中，我們將針對實際資料進一步完善方案(包括地震記錄中噪聲影響的分析及壓制，不同模態(tài)頻散曲線的分離和提取，以及各階模態(tài)加權(quán)系數(shù)的確定方法)，使其能夠更好地應用于實際資料.

致謝感謝兩位審稿人提出的非常具體的修改意見.