于 涵，劉 財，王 典，鹿 琪

吉林大學地球探測科學與技術學院，長春 130026

0 引言

淺層地震勘探的目的在于確定較小尺度下地質體的結構特征和地質剖面的物性參數(shù)，對解決淺層地質問題、探測淺層地下空間具有重要意義[1-3]。面波勘探為無損性淺表地震勘探，以面波的頻散特性為基礎。不同周期的面波以不同的波速進行傳播。周期越長的面波其縱向傳播距離越遠，越能反映較深處地下介質的響應特征。通過測量不同周期的面波速度，可以得到面波中不同頻率成分對應速度的關系(即頻散曲線)，這樣的頻率-速度關系受到地下介質的影響[4-5]。采取合理的數(shù)據(jù)處理方法，有效提取面波頻散曲線，是提升面波勘探精度的重要手段[6-9]。通過構建相應的反問題，利用反演方法可以推斷地下不同深度介質的屬性，從而達到工程勘探的目的[10-11]。

19世紀中期，英國物理學家Rayleigh[12]研究發(fā)現(xiàn)一種P波(縱波)與SV波(垂直偏振橫波)在自由表面干涉形成并沿自由表面?zhèn)鞑サ腞ayleigh波，對其產生的數(shù)學物理機理進行了論述。Love[13]從數(shù)學上證明了Love面波的存在。Thomson[14]、Haskell[15]發(fā)現(xiàn)了多層介質中Rayleigh波的頻散特性。受限于當時的計算能力，并未在實驗中得到充分驗證。隨著計算水平不斷的提升，面波的理論研究得到了快速的發(fā)展。Aki[16]通過地震噪聲源得到地下面波的頻散信息，首次利用噪聲源獲取面波信息。von Seggern等[17]從地震爆炸譜比角度，分析Rayleigh波層狀介質傳遞函數(shù)對源譜的影響，定義了Rayleigh波。Campillo等[18]對兩個地震臺站記錄的地震尾波進行了互相關處理，從互相關函數(shù)中識別出了Rayleigh波和Love波的信號，證明隨機波場具有頻散性質。同年，Hayashi 等[19]運用互相關方法從背景噪聲中提取出面波數(shù)據(jù)。Park等[20]在主被動源激發(fā)的地震數(shù)據(jù)中有效提取面波信息，在近地表速度成像方面取得了良好的應用效果。

頻散能量譜主要計算方法有：τ-p變換法(τ為垂直波慢度，p為水平波慢度)、高分辨率線性Radon變換法、頻率分解法、相移法、f-k變換法(f為頻率，k為波數(shù))、SFK(slantf-k)變換法和矢量波數(shù)變換法等。Mcmechan等[21]最早提出τ-p變換計算方法，將離散波列的線性Radon變換方法應用到Rayleigh波地震波場分離中。許多學者針對τ-p變換法的地震數(shù)據(jù)處理效果進行了進一步研究，針對端點效應和假頻效應等問題提出了改進方法，并取得了較好的效果[22-25]。Luo等[26]在線性Radon變換中引入加權系數(shù)，提出高分辨率線性Radon變換法。Xia等[27]提出了利用拉伸函數(shù)與地震數(shù)據(jù)褶積形成頻率掃描數(shù)據(jù)，之后進行傾斜疊加變換的頻率分解法。Park等[28]基于面波的相位信息進行傾斜疊加，加大了對相位信息的探索，抑制了空間假頻，增強了頻散能量譜上的相干信號，提高了頻散能量譜的分辨率，提出相移法求取瞬態(tài)Rayleigh面波頻散能量譜。伍敦仕等[29]和Cheng等[30]提出改進的相移法成像方法，一定程度抑制了空間假頻的干擾，并增強了頻散能量譜上的相干信號，提高了頻散能量譜的分辨率。Gabriels等[31]利用f-k變換提取出Rayleigh波頻散曲線。Serdyukov等[32]利用S變換改進了f-k變換，在時頻域自適應調解時窗函數(shù)，顯著放大面波信息，增強了頻散能量譜的低頻有效信號，提高了頻散能量譜的抗噪性，有效抑制空間假頻現(xiàn)象，稱之為SFK變換。楊振濤等[33]提出了矢量波數(shù)變換方法，在f-k變換法基礎上引入地震矩張量概念，提供了更豐富的高階模態(tài)面波頻散信息。近年來，大量學者結合先進的信號處理方法，利用不同的時頻分析方法擴展波束形成能力，選擇最佳觀測時窗提高識別時空上傳播波的靈活性以及方法的分辨率，進一步提高了面波的成像精度和成像質量，并驗證了其高效性和較好的應用潛力[34-39]。

本文分析了τ-p變換法、高分辨率線性Radon變換法、頻率分解法、相移法、f-k變換法、SFK變換法和矢量波數(shù)變換法等7種頻散能量譜成像方法的技術原理和特點，通過對理論模型合成的Rayleigh波記錄進行處理，比較了不同方法在f-v(速度)域頻散能量譜的成像效果，并對各方法的優(yōu)缺點進行了總結。

1 計算方法理論基礎

1.1 τ-p變換法

τ-p變換法又稱離散線性Radon變換法。用τ-p變換法求取頻散曲線，實質上是一種由兩個線性變換組成平面波的分解方法[17]。

將共炮點面波地震數(shù)據(jù)h(x,t)進行τ-p變換，沿直線方程t=τ+p進行傾斜疊加，其離散形式表示為

(1)

式中：τ為垂直波慢度；p為水平波慢度；x為偏移距；xmax、xmin分別為x的最大值和最小值；t為時間。式(1)沿時間方向進行一維傅里葉變換得

(2)

式中，U(f,p)和H(x,f)分別為u(p,τ)和h(x,t)的一維傅里葉表現(xiàn)形式。根據(jù)v=1/p，得

(3)

實際數(shù)據(jù)處理過程如下：

其中：①對h(x,t)沿不同斜率的直線方程進行積分，得到u(τ,p)；②沿時間方向做一維傅里葉變換，得到U(f,p)；③轉換至頻率速度域。

1.2 高分辨率線性Radon變換法

高分辨率線性Radon變換在Radon變換的基礎上引入加權函數(shù)，用于降低含噪數(shù)據(jù)對頻散能量譜成像效果的限制。

式(2)τ-p反變換的離散形式可以表示為

(4)

式中：pmax、pmin`分別為p的最大值和最小值。式(4)寫為矩陣形式為

H=LU。

(5)

式中：H為數(shù)據(jù)空間向量；U為模型空間向量；L為線性Radon變換算子矩陣。類似地，式(3)也可寫成矩陣形式：

U=LTH。

(6)

通過預加權的共軛梯度算法降低面波數(shù)據(jù)中的噪聲限制，提高分辨率。在方程(5)中引入模型權重矩陣Wu，得

(7)

反演過程中引入阻尼因子λ，得

U=(LTL+λI)-1LTH。

(8)

算法核心是求取目標函數(shù)最小解，求得U。目標函數(shù)M表示為

(9)

對式(9)中的U求導并令導數(shù)為零可得

(10)

式中：I為單位矩陣；Wh為數(shù)據(jù)權重矩陣。可通過共軛梯度法求解式(10)[26]。

實際數(shù)據(jù)處理過程如下：

其中：①對h(x,t)沿時間方向做一維傅里葉變換，得到H(x,f)；②對H(x,f)進行高分辨率線性Radon變換，沿不同斜率的直線方程進行傾斜疊加，得到U(f,p)；③轉換至頻率速度域。

1.3 頻率分解法

頻率分解法是基于τ-p變換原理的改進方法，通過褶積計算，將時間域轉換至頻率域，將原始地震數(shù)據(jù)轉換為連續(xù)頻率下的地震數(shù)據(jù)，并對其進行類似τ-p變換法的傾斜疊加變換，從而得到頻散能量譜。

對h(x,t)進行褶積：

Hr(x,t)=r(t)?h(x,t) 。

(11)

其中：

(12)

式中：Hr(x,t)為瞬時頻率地震記錄；r(t)為掃描樣本函數(shù)；f(t)為時間與頻率一一對應的瞬時函數(shù)。

實際數(shù)據(jù)處理過程如下：

其中：①對h(x,t)進行褶積，生成Hr(x,t)；②對Hr(x,t)進行τ-p變換(具體處理過程如1.1節(jié))，轉換至頻率速度域。

1.4 相移法

相移法本質是在頻率域內對相位信息進行傾斜疊加，是目前主動源面波頻散成像應用最廣泛的一項技術。

對h(x,t)沿時間方向做一維傅里葉變換，得

(13)

將式(13)表示為相位譜與振幅譜的乘積：

H(x,f)=P(x,f)A(x,f)。

(14)

式中：P(x,f)為相位譜，包含了頻散的所有信息；A(x,f)為振幅譜，包含了地震信號的其他信息，如振幅衰減、球面擴散等。

基于頻散特性，面波數(shù)據(jù)可表示為

U(k,f)=cos(2πft+2πkx+φ0)。

(15)

式中，φ0為初始相位。相位譜2πkx+φ0表示各頻率簡諧波在時間原點的相位，與偏移距相關。將相位譜用e-iΦx表示，則式(14)可表示為

H(x,f)=e-iΦx+φ0A(x,f)。

(16)

其中：

式中，vf為不同頻率對應的相速度。

對同一頻率面波的各相位信息進行空間積分，得

(17)

實際數(shù)據(jù)處理過程如下：

其中：①對h(x,t)沿時間方向進行傅里葉變換，得到H(x,f)；②歸一化頻譜，消除振幅譜影響，對各道地震波能量譜沿著空間水平方向疊加求和，得到U(φ,f)；③根據(jù)頻率與相位之間的關系，轉換至頻率速度域。

1.5 f-k變換法

f-k變換法的本質是對面波的頻散特征進行二維傅里葉積分變換。

對h(x,t)沿時間方向做一維傅里葉變換后(式(13))，沿空間方向做一維傅里葉變換，得

(18)

根據(jù)v=f/k，即可得到頻率速度域頻散能量譜U(f,v)。

實際數(shù)據(jù)處理過程如下：

1.6 SFK變換法

近年來時頻分析方法已廣泛應用于主動源和被動源面波頻散成像中，核心在于時頻變換窗函數(shù)的選取，在不同尺度下對面波信號進行分析，從而有效地提取面波頻散信息。時頻分析方法目前有很多種，本文僅詳細介紹基于S變換改進的f-k變換。

對第一道地震數(shù)據(jù)h1(t)進行S變換：

S[h1(t)](τ,f)=

(19)

式中，τ為時間軸上高斯窗的分析點位置。第二道數(shù)據(jù)h2(t)可以用h1(t)表達：

S[h2(t)](τ,f)=

(20)

e-i2πk(f)le-ζ(f)lS[h1(t)]](τ-k′(f)l,f)。

式中:ζ(f)為與頻率相關的衰減參數(shù);k(f)為與頻率相關的波數(shù)，k′(f)為k(f)的頻率導數(shù)；l為兩道數(shù)據(jù)之間的距離。將S變換應用于地震記錄中，形成同一頻率的偽面波地震記錄：

hf(x,τ)=S[h(x,t)](τ,f,x)。

(21)

式中，hf(x,τ)為S變換后的共頻率偽面波記錄。

假設震源在零時間產生一個沖量和偏移量，引入基于S變換的面波分組模型[35]：

(22)

(23)

式中，μf(x,τ)為誤差項，包含其他地震波、噪聲和建模圖像誤差。

固定時間、頻率，在群速度范圍內對偽面波地震記錄進行傾斜切片，式(21)可表示為

(24)

式中，huf(x)為傾斜相位函數(shù)。再沿x方向做傅里葉變換，得

(25)

尋找群速度最大振幅：

(26)

實際數(shù)據(jù)處理過程如下：

其中：①對h(x,t)進行S變換生成hf(x,τ)；②將hf(x,τ)固定時間、頻率，在群速度范圍內進行傾斜切片，得到huf(x)；③沿時間方向進行傅里葉變換并取極大值得到U(f,k)，轉換至頻率速度域。

1.7 矢量波數(shù)變換法

楊振濤等[33]近年來在f-k變換法的基礎上，研究垂直于地面單力點源產生的地震，基于不同階貝塞爾函數(shù)對應的實部與虛部頻散譜不同，以波數(shù)維度對稱性為核心思想，引入地震矩張量概念，對f-k變化進行優(yōu)化。

在地表面觀測到的垂直分量地震波場記錄可表示為

h(x,t)=f(t)*g(x,t)。

(27)

式中，g(x,t)為震源與觀測點間的格林函數(shù)。對式(27)進行傅里葉變換，得

H(x,ω)=F(ω)G(r,ω)。

(28)

式中：G(r,ω)為震源與觀測點間在頻率域的格林函數(shù)；r=|x|；ω為角頻率。自由表面各向同性震源的格林函數(shù)在頻率域的計算公式可以表示為

(29)

式中：g(ω,k)為格林核函數(shù)；J0(kr)是第一類零階貝塞爾函數(shù)。

(30)

式中，θ為方位角。

對觀測波場的頻譜進行矢量波數(shù)變換[33]：

(31)

當空間為水平層狀介質時，矢量波數(shù)變換結果僅依賴波數(shù)矢量的模k(k=|k|)，與波數(shù)矢量的方向無關，式(31)可以表示為

(32)

在頻率域中式(32)可以表示為

(33)

根據(jù)貝塞爾函數(shù)的正交性[40]：

(34)

將式(34)代入式(33)，簡化為

(35)

g(ω,k)與面波頻散特性的久期函數(shù)值成反比，矢量波數(shù)變換依據(jù)此特性進行頻散能量譜計算。

根據(jù)v=ω/k、ω=2πf/v，可以得到頻率速度域頻散能量譜U(f,v)。

實際數(shù)據(jù)處理過程如下：

其中：①對h(x,t)沿時間方向進行傅里葉變換轉換至頻率域，得到H(r,ω)；②對H(r,ω)進行多道地震數(shù)據(jù)加權求和，根據(jù)式(35)獲得g(ω,k)；③對g(ω,k)進行相速度和波數(shù)掃描，獲得近似格林核函數(shù)圖，根據(jù)速度、頻率和波數(shù)之間的關系，轉換至頻率速度域。

2 頻散能量譜效果對比

本文建立了三層速度遞增水平層狀介質模型(表1)，震源為主頻20 Hz、延遲時間0.03 s的高斯子波，將其布置在(10 m, 0 m)處。圖1a為模型合成的Rayleigh波地震記錄，最小偏移距是5 m，道間距為2 m，每道的采樣點數(shù)為6 000，采樣時間間隔為0.2 ms，記錄延續(xù)時間為1.2 s。由圖1a可見，利用該模型模擬的60道Rayleigh波地震記錄(垂直分量)清晰可見，發(fā)散的掃帚狀十分明顯。利用此合成記錄，應用上文所述7種方法變換得到頻散能量譜，并與多模式快速矢量傳遞算法求取的理論頻散曲線[41](頻率為5～100 Hz，間隔1 Hz)進行對比，結果如圖1b—h所示。

表1 三層速度遞增水平層狀介質模型參數(shù)

由圖1可見，文中所述7種方法均可對面波頻譜有效成像，對于基階面波頻散曲線變化趨勢具有相似的特征，識別效果均較為清晰。f-k變換法在頻率20 Hz以內頻散能量譜極值點與解析解誤差較大， 低頻成像精度低，高階模態(tài)頻散能量信息不明顯；由于速度與波數(shù)之間的映射關系，矢量波數(shù)變換法和f-k變換法存在明顯的空間假頻現(xiàn)象，可能由于計算過程中插值不均導致；頻率分解法在15 Hz以內的頻散曲線存在明顯波動，能量較為豐富；對于多模態(tài)下的面波頻散信息，f-k變換法無法有效拾取高階面波頻散曲線，頻率分解法、矢量波數(shù)變換法、相移法、高分辨率線性Radon變換法頻散能量譜分辨率較高，可有效提取出高階面波頻散曲線。

a. 合成Rayleigh波地震記錄。頻散能量譜：b. τ-p變換法；c. 高分辨率線性Radon變換法；d. 頻率分解法；e. 相移法；f. f-k變換法；g. SFK變換法；h. 矢量波數(shù)變換法。圖中白點為頻散曲線的解析解。

對文中7種方法的計算耗時進行比較，結果如表2所示。由表2可見，f-k變換法計算速度最快，頻率分解法、矢量波數(shù)變換法和SFK變換法計算速度較慢。

近年來，各種頻散能量譜計算方法已經在面波工程勘探領域廣泛應用。多種面波勘探方法，如共中心點互相關二維面波法、瞬態(tài)多道面波分析法和面波譜分析法等已被國內外多家生產和科研單位使用，在工程勘探中也有深入探討[42-47]，面波勘探技術已在近地表地質勘察、石油勘探、隧道檢測和大壩穩(wěn)定性等實際應用中取得良好效果[48-53]。本文針對不同面波勘探方法適用的地質環(huán)境、采集方式和適用空間范圍等方面進行整理，結果如表3所示。文中詳細闡述的7種數(shù)據(jù)計算方法，多應用于主動源及被動源面波勘探二維多道采集數(shù)據(jù)處理過程中。近年來，許多學者應用不同頻散能量譜計算方法對實際數(shù)據(jù)進行了處理[54-60]，本文對其應用特點進行總結。

τ-p變換法：計算量小、失真小、易于實現(xiàn)、壓制端點效應好，對多階模態(tài)面波頻散信息成像效果較好；但對基階模態(tài)低頻面波頻散信息的成像效果較差，會出現(xiàn)波場能量延伸和假頻現(xiàn)象，結果易導致反變換后的波場相互疊加干涉。

表2 頻散能量譜計算耗時

表3 面波勘探方法特征對比

高分辨率線性Radon變換法：抗噪性強，可解決缺失數(shù)據(jù)重建問題，有效進行多模式面波頻散能量成像；但對信噪比較低的數(shù)據(jù)難以形成連續(xù)的頻散能量，成像效果不佳。

頻率分解法：抗噪性強，檢波器分布不局限于等距直線排列，即對以任意幾何形狀參數(shù)排列的檢波器采集到的Rayleigh波數(shù)據(jù)均能使用，適用范圍廣；但計算量大，處理數(shù)據(jù)較慢。

相移法：抗噪性強、計算量小、分辨率較高，對基階模態(tài)的面波頻散信息提取效果較好，降低Rayleigh波勘探方法對道間距的要求，使用較少道集即可直接產生分辨率較高的面波頻散能量成像；但對高階模態(tài)面波頻散信息分辨率不高。

f-k變換法：計算量小，計算速度快；但要求采集時間間隔以及空間間隔相等，采集數(shù)據(jù)質量要求嚴格，容易發(fā)生空間假頻問題，頻散能量集中在一個很窄的頻帶范圍，幅值最大原則下提取的頻散曲線不穩(wěn)定。

SFK變換法：抗噪性能強，可有效降低數(shù)據(jù)轉換過程中頻散能量的損失，頻帶能量分布均勻；但計算量大，耗時較長。

矢量波數(shù)變換法：抗噪性強，適用范圍廣，在成像分辨率和高階模態(tài)成像質量上具有明顯的優(yōu)勢，較少的地震道也能得到高質量的高階模態(tài)成像；但計算量大，耗時較長。

3 結論

本文分析了τ-p變換法、高分辨率線性Radon變換法、頻率分解法、相移法、f-k變換法、SFK變換法和矢量波數(shù)變換法等7種頻散能量譜成像方法的技術原理和特點，得到以下認識：

1)7種計算方法均可有效地進行頻散能量譜成像，在基階頻散曲線提取時效果差異較小，τ-p變換法、高分辨率線性Radon變換法、相移法和矢量波數(shù)變換法提取頻散曲線精度較高，f-k變換法精度較低；在高階模態(tài)頻散曲線提取時效果差異較大，τ-p變換法、頻率分解法、高分辨率線性Radon變換法和矢量波數(shù)變化法精度較高，相移法、f-k變換法精度較低。在多模態(tài)面波頻散提取過程中，高分辨率線性Radon變換法和矢量波數(shù)變換法提取精度較高，f-k變換法提取精度最低。

2)τ-p變換法、高分辨率線性Radon變換法、相移法、f-k變換法高效簡易、計算速度較快；頻率分解法、矢量波數(shù)變換法、SFK變換法計算量較大，處理數(shù)據(jù)耗時較長?；诟叻直媛示€性Radon變換法、SFK變換法、相移法、頻率分解法、矢量波數(shù)變換法的面波頻散能量譜計算抗噪性較強，在實際面波數(shù)據(jù)處理中，可更好地拾取面波頻散曲線。