面波頻散譜多模式高分辨率成像的多道信號比較法

易佳， 劉伊克 ， 胡昊， 張郁山， 楊宗奇

1 中國地震災害防御中心， 北京 100029 2 中國科學院地質與地球物理研究所， 北京 100029 3 University of Houston, Houston Texas, 77204-5007, USA 4 中國石油集團東方地球物理勘探責任有限公司， 河北 涿州 072750

0 引言

在淺層地震勘探領域，由于面波在地震波場能量中占主導部分，而面波攜帶有豐富的地下介質的橫波速度信息，因此橫波速度被認為是一種十分重要的地震學參數(Aki and Richards, 1980).在地震工程領域，橫波速度是確定場地土類別的重要參數和主要依據，從而用于場地地震反應分析計算，最終服務于場地的地震安全性評價(尤紅兵等, 2009; 張郁山, 2009, 2010; Wang et al., 2019).因此，探明淺層精細的橫波速度結構對于以上兩個領域工作的開展意義重大.

自20世紀50年代面波的頻散性質被發現后，在淺層速度結構的分析中，研究者們便利用面波的頻散特性反演估計淺層橫波速度結構.隨著瞬態面波法的興起和應用，面波頻譜分析法(Stokoe and Nazarian, 1983; Stokoe et al., 1994)、多道面波分析法(Park et al., 1998; Xia et al., 1999, 2012)等通過提取面波的頻散曲線，采用最小二乘局部優化算法或全局優化方法，反演頻散曲線估計一系列的橫波速度隨深度變化的剖面(Beaty et al., 2002; Luo et al., 2007).通過對一維剖面的插值處理，可重建擬二維的橫波速度結構(Boiero and Socco, 2010; Bergamo et al., 2012; Mi et al., 2017).隨著波形反演(Pérez Solano et al., 2014; Liu et al., 2018)和波動方程走時反演(Yi et al., 2019)的飛速發展，基于波動方程的頻散曲線反演(Li et al., 2017; Liu et al., 2019)和二維頻散譜反演方法(Masoni et al., 2014; Zhang and Alkhalifah, 2019)隨即發展起來，該類方法基于波形反演和波動方程走時反演的思想，通過建立頻散曲線或二維的頻散譜相關的目標函數，反傳頻散曲線或頻散譜的殘差直接反演獲得二維的橫波速度結構.

基于面波的頻散特性反演獲取橫波速度結構的方法，均涉及一個關鍵步驟：對面波的頻散能量進行成像，獲得面波的頻散譜.因此，高分辨率的面波頻散譜成像對于準確提取頻散曲線十分關鍵.此外，面波的高階模式相比于基階模式對于地層參數、軟弱夾層等更加敏感，且高階模式具有更深的勘探深度，因此，考慮精確的多模式的面波頻散譜成像能有效提高反演的精度以及探測深度(Luo et al., 2007; Zhang and Alkhalifah, 2019).在主動源面波勘探中，常用的面波頻散譜成像方法有：τ-p變換法(McMechan and Yedlin, 1981)，f-k變換法(Yilmaz, 1987)，相移法(Park et al., 1998)，頻率分解傾斜疊加法(Xia et al., 2007)，高分辨率線性拉東變換法(Luo et al., 2008)，以及矢量波數變換法(楊振濤等, 2019).τ-p變換法將數據按不同的截距和斜率的直線進行疊加，對高階模式成像質量較好，但對于基階模式的成像分辨率較差，且易出現端點效應和假頻問題 (邵廣周和李慶春, 2010).f-k變換法對數據進行二維傅里葉變換，原理簡單，但要求時間和空間采樣等間隔，對采集方式要求較高，不能有空道和壞道，否則將大大影響頻散譜的成像精度.相移法將各道數據沿時間方向做傅里葉變換，再沿空間方向進行積分，不要求檢波器等間隔，但由于每個頻率只能得到一個相位值，單個頻率的頻譜為基階和高階的耦合，因此相移法對于高階模式不能很好地成像.頻率分解傾斜疊加法首先通過頻率掃描，得到偽炮集記錄，再通過傾斜疊加的方式獲得面波的頻散譜，該方法在疊加前對頻率進行了分解，因此，相比于τ-p變換法有更高的成像精度，但該方法在頻率掃描階段存在許多無效計算及頻率偏差的問題(沈超, 2017).高分辨率線性拉東變換法將頻散譜成像設計成一個反演過程，通過迭代可以逐漸提高成像的分辨率(Luo et al., 2008).矢量波數變換法通過對觀測波場的頻譜進行矢量波數變換，在頻率波數域對變換后的波場進行掃描獲得面波的頻散譜.該方法對基階和高階模式均有較高的成像分辨率，但需要已知震源子波的信息，因此在實際資料震源未知的情況下，還需要對震源做近似處理(楊振濤等, 2019).

在天然地震探測領域，雙臺法(Knopoff et al., 1966; Yao et al., 2006; Foster et al., 2014)又被稱為LSC方法(Zheng and Hu, 2017)，被廣泛用來對面波的頻散能量進行成像，提取頻散曲線并用于大尺度的面波層析成像.該方法通過對兩個臺站記錄的來自同一方位的面波進行互相關信號比較，即可獲得面波的頻散譜，但互相關的成像方式對于低頻信號極不敏感，會造成頻散譜在低頻端較低的分辨率，以及整個頻段分辨率的不均勻，導致無法準確拾取頻散曲線.NLSC方法(Zheng and Hu, 2017)利用指數函數克服了低頻端分辨率較低的問題，同時引入了一個可調參數，極大地提高了頻散譜的成像分辨率.主動源面波勘探觀測系統與天然地震觀測系統不同，通常沿著直線布設多個檢波器接收地震信號，因此，在信號比較法的基礎上，考慮多道的MSC方法(Hu et al., 2019)被提出用于多道地震數據的面波頻散譜成像，對于缺道以及噪聲問題，MSC方法也具有較強的穩健性，與此同時，MSC方法也被用于地震數據的面波分離當中，取得了良好的分離效果(Hu and Zheng, 2019).

在研究中我們進一步發現，LSC和NLSC方法利用兩道地震記錄，并不能將面波的頻散特性完整地考慮在內，僅對基階模式成像準確，高階模式的成像存在較大的誤差.MSC方法考慮了完整的面波頻散信息，能獲得準確且分辨率較高的面波頻散譜多模式成像，但該方法需要計算任意兩道的頻散譜并進行疊加，對于多道記錄，存在冗余的計算耗時較長.因此，我們對MSC方法進行了相應的改進，以求高效獲得頻散譜的多模式高分辨率精確成像，從而有助于多模式頻散曲線的精確提取，以及利用多模式頻散譜進行高精度的淺層地震勘探.本文系統性地對LSC、NLSC方法的理論、關鍵技術以及存在的問題進行了闡述分析，并給出了改進的MSC方法的原理.通過對含低速夾層的水平層狀模型及一套實際地震資料的Rayleigh波頻散譜的成像，深入對比研究了幾種基于信號比較理論方法在面波頻散譜成像中的應用效果，從而為不同的面波頻散譜分析研究提供合適的解決方案.

1 基本原理

1.1 傳統線性信號比較法 (LSC)

在天然地震面波頻散性質的分析當中，LSC方法利用互相關的方式比較兩個臺站接收到的地震信號，獲得它們之間的相位差，這個相位差與兩個接收點之間的距離及面波的相速度有關，通過這種方式即可獲得面波的頻散譜.其方法原理如圖1所示.震源激發的地震波，傳播經過同一個方位的兩個臺站Rec 1和Rec 2被接收到，兩者之間的距離為x.Rec 1處的地震記錄表示為R1(t)，Rec 2處的地震記錄表示為R2(t).如果能測得Rec 2處特定頻率的面波到達的延遲時τ，那么就可以反推獲得面波的相速度v.

圖1 傳統線性信號比較法原理示意圖Fig.1 The conventional LSC method

然而，通常兩個臺站之間的距離是已知的，因此，只需給出一定的相速度范圍，移動地震記錄R2(t)，求取與R1(t)的互相關，即可獲得面波的頻散譜(Zheng and Hu, 2017)，記為ELSC(ω,v)：

(1)

1.2 非線性信號比較法 (NLSC)

為克服傳統LSC方法分辨率的限制，并保證整個頻帶范圍內均勻分布的分辨率.Zheng和Hu (2017)提出了NLSC方法：

(2)

其中，ENL表示未歸一化的頻散譜，ε∈[0,+∞)是一個非負實數，用來控制頻散譜的整體分辨率.Zheng和Hu(2017)定義了一個歸一化算子Eπ：

Eπ(ω,v;ε)=

(3)

對ENL進行歸一化處理，即可得到歸一化的高分辨率頻散譜，記為ENLSC：

(4)

歸一化算子Eπ可認為是兩個臺站記錄的地震信號的相位差為π時，由(2)式計算得到.NLSC方法巧妙地運用指數函數的性質，克服了傳統LSC方法生成的頻散譜在整個頻帶范圍分布不均勻的問題，提高了低頻端的分辨率，同時引入參數ε可以有效提高頻散譜的整體成像分辨率.

1.3 多道信號比較法(MSC)

由于探測精度的需求，主動源面波勘探觀測系統與天然地震觀測系統不同，通常沿著直線布設多個檢波器接收地震信號，從而可利用多道信息反演獲得小尺度的更為精細的地下結構.在研究中我們發現，只利用兩道地震記錄的LSC和NLSC方法，并不能將面波的頻散特性完整地考慮在內，只對基階模式有較準確的成像，對于高階模式的成像存在較大的誤差.因此，Hu等(2019)基于勘探地震中多道采集的方式提出了MSC方法：

(5)

其中，i,j表示地震道，n表示總的接收道數，Ei,j表示利用信號比較法計算的任意兩道的頻散譜，E′MSC表示將所有兩道頻散譜疊加的頻散譜，采用歸一化方法，對頻散譜進行歸一化處理，即可得到MSC方法計算的面波頻散譜：

EMSC(ω,v)=

(6)

圖2 改進的多道信號比較法原理示意圖Fig.2 The modified multichannel signal comparison method

(7)

同理，改進的多道非線性信號比較法(MNLSC)的可表示為

(8)

改進的MSC方法追蹤炮集記錄上的面波波組，考慮了完整的面波頻散信息同時不增加冗余的信息，可以有效提高計算效率.

2 數值實驗

2.1 傳統LSC和NLSC方法的頻散譜成像分析

為了研究幾種信號比較法在淺層地震勘探中的成像能力，采用含低速夾層的水平層狀模型(LVL)進數值實驗，模型參數如表1所示.利用垂直震源激發地震波，子波是主頻為20 Hz的Ricker子波，自由邊界條件，采用SPECFEM2D(Komatitsch and Tromp, 1999)模擬含有Rayleigh波的地震波場.線性觀測系統，100個檢波點均勻分布在地表，道間距1 m，采樣間隔0.1 ms，一共記錄時長為1 s.獲得的炮集地震記錄的垂直分量如圖3所示，可以看到地震記錄中Rayleigh波的能量占了主導成分，且可觀測到明顯的面波的頻散特性.

表1 含低速夾層的水平層狀模型參數Table 1 Parameters of the LVL model

圖3 含低速夾層的水平層狀模型的炮集地震記錄的垂直位移分量Fig.3 The vertical displacement component seismograms of the LVL model

抽取其中的兩道地震記錄(如圖3中紅色地震道所示)，利用LSC和NLSC方法進行面波頻散譜的成像計算.此外，還利用Knopoff算法(Schwab and Knopoff, 1970)計算得到該模型的理論頻散曲線疊加在頻散譜上，用以判斷數值方法計算的頻散譜的準確性.所得的Rayleigh波頻散譜能量圖如圖4所示，可以看到LSC方法(圖4a)的頻散譜在給出的頻帶范圍內呈現出非均勻的分辨率，高頻端的分辨率較高，低頻端分辨率較低，在小于10 Hz區域成像的分辨率尤其較低.而NLSC方法(圖4c，ε=0.01時)成像的分辨率得到了有效地提高，尤其在低頻端得到了顯著地提升.當調整分辨率的參數ε=0.001時(圖4d)，成像分辨率達到極高.在頻散譜的準確性上，LSC和NLSC方法均呈現出強而清晰的基階模式Rayleigh波頻散能量，并且與理論頻散曲線的基階模式有較好地吻合.然而，兩種方法的成像結果，在高頻端均出現了扭曲，與理論頻散曲線吻合較差，且高階模式雖有著清晰的成像，卻與理論頻散曲線相差甚遠.因此，只利用兩道地震記錄的信號比較法，無法對面波的多模式進行準確地成像，因此無法利用其進行高精度的淺層結構反演.

圖4 LSC和NLSC方法的Rayleigh波頻散譜能量圖(a) LSC方法； (b) NLSC方法,ε=0.01； (c) NLSC方法,ε=0.001，圖中白色圓點為Knopoff算法計算獲得的理論頻散曲線.Fig.4 The Rayleigh wave dispersion spectrum maps of the LVL model using (a) LSC method an (b) NLSC method with ε=0.01 and (c) NLSC method with ε=0.001. The white dots in the figure are the theoretical dispersion curves calculated by the Knopoff′s algorithm

2.2 改進的MSC方法頻散譜成像分析

分別采用24和100道的地震數據，應用改進的MLSC和MNLSC方法對該地震記錄(圖3)進行頻散譜的成像計算，為了對比展示成像效果，利用相移法也獲得了面波的頻散譜，將三種方法與理論頻散曲線疊加，結果如圖5和圖6所示.可以看到，由于低速夾層的存在，三種方法均呈現出了清晰的多模式面波能量.圖5顯示，當參與計算的道數較少時，頻散譜的成像效果均較差，相移法在低頻端與理論頻散曲線吻合較好，但分辨率較低，而MLSC和MNLSC方法的成像分辨率以及高階模式的成像精度均明顯高于相移法.當參與頻散譜成像的地震道數增加時(圖6)，三種方法的成像質量都獲得了顯著地提升，表明面波頻散譜的分辨率和準確性會隨著道數的增加而提高.圖6a顯示，相移法的成像分辨率在整個頻段上分布極不均勻，低頻端(<10 Hz)的成像分辨率較低，使得低頻端能量峰值與理論頻散曲線出現較大偏差.與相移法相比，MLSC(圖6b)和MNLSC方法(圖6c，6d)的成像質量具有明顯優勢，整個頻段的分辨率均有所提升，尤其在低頻端有顯著的提高，與理論頻散曲線吻合較好.MNLSC方法相比于MLSC方法在成像分辨率上則有更近一步的提升.因此，改進的MLSC和MNLSC方法，利用多道地震信號，在提升計算效率的同時，不僅能有效地對面波的多階模式進行精確地成像，而且能獲得均勻且分辨率較高的頻譜，這在面波頻散譜的成像上具有巨大的潛力.

圖5 24道數據的Rayleigh波頻散譜能量圖(a) 相移法； (b) MLSC方法； (c) MNLSC方法,ε=0.03； (d) MNLSC方法,ε=0.01.Fig.5 The Rayleigh wave dispersion spectrum maps with 24 traces using (a) phase shift method (b) MLSC method (c) MNLSC method with ε=0.03 and (d) MNLSC method with ε=0.01

圖6 100道數據的Rayleigh波頻散譜能量圖(a) 相移法； (b) MLSC方法； (c) MNLSC方法,ε=0.03； (d) MNLSC方法,ε=0.01.Fig.6 The Rayleigh wave dispersion spectrum maps with 100 traces using (a) phase shift method (b) MLSC method (c) MNLSC method with ε=0.03 and (d) MNLSC method with ε=0.01

3 實際地震資料測試

為研究改進的MSC方法在實際地震資料中的適用性，將其應用在一套紅海海岸跨越Qademah斷層采集的含Rayleigh波的實際地震資料中(Hanafy et al., 2015; Li et al., 2017; Zhang and Alkhalifah, 2019).采集參數為：垂直震源，120個檢波器，道間距5 m，采樣間隔 0.05 ms, 一共記錄時長1.2 s.圖7a展示了該資料沿著時間方向進行能量補償(Claerbout, 1985)之后的典型單炮記錄，我們對其進行了包括剔除壞道，開窗切除保留面波成分等預處理工作，預處理之后的炮集記錄如圖7b所示.由于該炮集記錄遠偏移距的地震道未能記錄到完整的面波，因此，我們僅利用預處理之后的前64道進行多道面波頻散譜的計算.

圖7 實際資料的典型單炮地震記錄(a) 能量補償后的單炮記錄； (b) 預處理之后的單炮記錄.Fig.7 The common shot gathers of the field data after (a) energy compensation, and (b) after preprocessing

分別采用相移法、改進的MLSC和MNLSC方法對該實際地震資料進行Rayleigh波的頻散譜成像計算，結果如圖8所示，可以看到MLSC(圖8b)和MNLSC(圖8c和圖8d)頻散譜成像的分辨率，在高頻端和低頻端，均優于相移法(圖8a).此外，對比圖中紅色箭頭處，MLSC和MNLSC相比于相移法對高階模式有更清晰的成像，對比白色箭頭處，MLSC和MNLSC則清晰地指示了模式的分離.相比于相移法，MLSC和MNLSC方法頻散譜的成像細節，能更好地指示地下橫波速度的改變.

圖8 實際地震資料Rayleigh波頻散譜能量圖(a) 相移法； (b) MLSC方法； (c) MNLSC方法,ε=0.05； (d) MNLSC方法,ε=0.02.Fig.8 The Rayleigh wave dispersion spectrum maps of the field data using (a) phase shift method (b) MLSC method (c) MNLSC method with ε=0.05 and (d) MNLSC method with ε=0.02.

4 討論

在數值試驗和實際地震資料的測試中我們可以看到，改進的MSC方法不需要計算每兩道地震信號的頻散譜，提高了計算效率，但同時充分利用接收的地震信息，獲得了準確的多模式高分辨率的頻散譜成像.在抗噪性方面，從方法原理可以看出，改進的MSC方法的成像效果較為依賴參考道的數據質量，若參考道數據質量較差，則會對成像結果有較大的影響，因此，在實際資料的處理中，可在近偏移距中選取數據質量較高的地震道作為參考道進行成像計算.對于MNLSC方法中的參數ε可依據反演需求來設置：當采用直接基于面波頻散譜的反演方法獲得二維橫波速度結構時(Zhang and Alkhalifah, 2019)，可將ε設置為一個稍大的數值，或者直接采用MLSC方法來獲得頻散能量清晰且分布均勻的頻散譜；當采用傳統提取頻散曲線的方法，估計一維橫波速度(Xia et al., 1999; Boiero and Socco, 2010) 或基于波動方程走時反演理論，利用頻散曲線直接反演獲得二維速度結構時(Li et al., 2017)，可將ε設置為一個較小的數值，以獲得高分辨率的頻散譜.因此，改進的MSC方法的另一個優勢是，可依據反演需求獲得合適的面波頻散譜的分辨率.理論上，當ε趨于零時，MNLSC方法的頻散譜分辨率將達到極致，但當ε取值較小時，低頻端的成像質量也會有所損失，如圖9所示.因此，在利用MNLSC方法計算時，合理地對參數ε取值，對于頻散譜的成像質量也十分關鍵.

圖9 LVL模型MNLSC方法 (ε=0.002)Rayleigh波頻散譜能量圖Fig.9 The Rayleigh wave dispersion spectrum maps of the LVL model using MNLSC method with ε=0.002

5 結論

天然地震探測中的兩道信號比較法(LSC方法和NLSC方法)在面波頻散譜的計算中存在一定的局限性：低頻端較低的分辨率，整個頻段非均勻的分辨率，高階模式不準確的成像，這些問題導致無法準確拾取頻散曲線，以及無法有效利用面波的高階模式反演獲得較為準確的橫波速度結構.本文根據主動源面波勘探中多道采集的方式，提出的改進的多道信號比較法(MSC方法)通過追蹤面波波組，考慮了完整的面波頻散信息的同時不增加冗余的計算，可以獲得高精度、高分辨率的面波多模式頻散譜的成像，該方法原理簡單，易于實現，并且理論上不需要觀測系統等間隔排列，降低了實際數據采集時對場地的要求.因此，在淺地表地震勘探以及地震工程的橫波速度結構探測中，多道信號比較法將具有廣闊的應用前景.

致謝感謝兩位評審專家提出的寶貴意見，感謝沙特阿卜杜拉國王科技大學Schuster教授及其課題組成員在采集地震數據過程中的辛苦奉獻，感謝Zhengdong Zhang博士提供的實際地震數據以及相關討論.