面波頻散能量譜計(jì)算方法

于 涵，劉 財(cái)，王 典，鹿 琪

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

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

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

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

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

1 計(jì)算方法理論基礎(chǔ)

1.1 τ-p變換法

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

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

(1)

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

(2)

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

(3)

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

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

1.2 高分辨率線性Radon變換法

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

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

(4)

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

H=LU。

(5)

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

U=LTH。

(6)

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

(7)

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

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

(8)

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

(9)

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

(10)

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

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

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

1.3 頻率分解法

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

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

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

(11)

其中：

(12)

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

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

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

1.4 相移法

相移法本質(zhì)是在頻率域內(nèi)對相位信息進(jìn)行傾斜疊加，是目前主動(dòng)源面波頻散成像應(yīng)用最廣泛的一項(xiàng)技術(shù)。

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

(13)

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

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

(14)

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

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

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

(15)

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

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

(16)

其中：

式中，vf為不同頻率對應(yīng)的相速度。

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

(17)

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

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

1.5 f-k變換法

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

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

(18)

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

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

1.6 SFK變換法

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

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

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

(19)

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

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

(20)

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

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

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

(21)

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

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

(22)

(23)

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

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

(24)

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

(25)

尋找群速度最大振幅：

(26)

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

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

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

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

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

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

(27)

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

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

(28)

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

(29)

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

(30)

式中，θ為方位角。

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

(31)

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

(32)

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

(33)

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

(34)

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

(35)

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

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

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

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

2 頻散能量譜效果對比

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

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

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

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

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

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

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

表2 頻散能量譜計(jì)算耗時(shí)

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

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

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

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

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

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

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

3 結(jié)論

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

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

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