李欣欣，李江，劉軍，沈鴻雁

(1.西安石油大學 地球科學與工程學院，陜西 西安 710065；2.西安石油大學 陜西省油氣成藏地質學重點實驗室，陜西 西安 710065；3.中國石化地球物理重點實驗室，江蘇 南京 211103；4.中煤科工集團西安研究院有限公司，陜西 西安 710077)

0 引言

瑞利波是一類沿地球表面?zhèn)鞑サ牡卣鸩?，其質點在豎直平面內振動，在近地表的淺部其質點振動的軌跡為“逆進”的橢圓[1,2]。瑞利波的能量衰減較穿透地球內部的體波更慢，在垂向地震記錄中有較強的振幅和較高的信噪比，因此利用其進行地下結構探測具有高效便捷和抗干擾能力強的優(yōu)勢[3]。早期的瑞利波探測方法利用兩道檢波器采集信號，然后通過計算兩道之間不同頻率瑞利波的相位差獲得頻散曲線，并反演下方地層的速度結構[4]。這類方法由于僅有兩道數(shù)據(jù)，若記錄中含有較強的反射波、折射波以及噪聲信號時，計算的相位差會包含很大的誤差，導致最終的探測結果不準確。因此研究人員發(fā)展了多道面波分析(multichannel analysis of surface waves)方法[5]，該方法利用多個檢波器采集瑞利波信號并進行處理，相比于兩道面波方法極大地提高了頻散曲線計算的準確性和穩(wěn)定性，且能夠同時獲得基階和高階模式的頻散曲線，為后續(xù)反演處理提供了更多的約束，提高了勘探的質量，因此多道瑞利波方法在煤層采空區(qū)探測[6,7]、工程勘察[8,9]、油氣地震資料靜校正[10,11]等不同領域的地下結構探測中獲得廣泛的應用。近年來快速發(fā)展全波形面波成像方法能夠直接從時間域多道面波(瑞利波或勒夫波)波形記錄反演地下結構，雖然在理論上可以免去頻散曲線的計算，但是在實際應用中成像的質量高度依賴于反演的初始模型，需要借助頻散曲線建立較為合理的初始模型[12-14]，并將頻散曲線作為約束因子引入全波形反演計算，以提高反演過程的穩(wěn)定性與收斂性[15]。因此不論對于常規(guī)的多道面波分析方法還是全波形面波反演方法而言，獲得準確的頻散曲線都是瑞利波勘探中的關鍵所在[16-19]。

針對多道面波頻散曲線計算這一個關鍵問題，目前已有諸如F-K變換[20]、相移法[21]、τ-p變換[22]、傾斜疊加[23]、高分辨率線性拉東變換[24]等各類方法，其中相移法由于對占主導的基階模式瑞利波的頻散能量成像效果較好[25]，在實踐中得到廣泛的應用。然而當采集數(shù)據(jù)中波場較復雜或有效信號信噪比較低時，相移法生成的瑞利波頻散能量會受到影響，導致提取頻散曲線的精度降低，因此本文提出一種改進的相移法對數(shù)據(jù)進行處理。

1 方法原理

相移法通過傅里葉變換生成頻散能量譜，使得多道記錄中不同頻率成分的瑞利波相互分開。頻散能量譜包含相位譜和振幅譜兩部分，瑞利波的走時信息包含在相位譜中，振幅衰減、球面擴散等信息則包含在振幅譜中。通過對相位譜和振幅譜進行處理，就可以獲得瑞利波的頻散曲線，本文對該方法進行改進，現(xiàn)將原理介紹如下。

1.1 相移法

多道地震記錄可視為空間-時間域的二維信號d(x,t)，其中x表示偏移距，t表示記錄時間，對該信號沿時間軸做傅里葉變換可得：

，

(1)

上式計算得到的頻譜D(x,f)可以表示為振幅譜和相位譜的乘積：

D(x,f)=ph(x,f)A(x,f)

，

(2)

式中ph(x,f)表示信號的相位譜，A(x,f)表示信號的振幅譜。若將相位譜利用e-iΦx表示，則上式可寫為：

D(x,f)=e-iΦxA(x,f)

，

(3)

式中Φ=2πf/v(f)，v(f)為隨頻率變化的相速度。若將某一頻率的瑞利波經過φ相位移動并進行振幅歸一化，然后沿著偏移距方向疊加，則可由式(3)得：

，

(4)

由式(4)可知，對于給定頻率f，如果滿足：

，

(5)

則Z(φ,f)將會有一個極大值，那么根據(jù)Z(φ,f)譜上的極大值點就可以確定該極大值點對應的φ值。因此，在計算得到Z(φ,f)后，可利用式(5)將f-φ域的Z(φ,f)能量圖的每一個能量點所對應的φ值換算為v值[25,26]，即可獲得Z(v,f)，再將Z(v,f)繪制在f-v域坐標系中，即可得到f-v域的頻散能量圖，進而提取頻散曲線。

1.2 改進的相移法

實際中采集的數(shù)據(jù)由于包含噪聲和各類型波場，利用相移法計算得到頻率譜后，其中的瑞利波能量受其他信號能量干擾，不易區(qū)分識別，且頻散能量條帶的收斂性和聚焦性較差，使得提取的頻散曲線準確性降低，甚至在部分頻段，瑞利波的頻散能量受到壓制導致無法獲取頻散信息。因此本文提出一種改進的相移法對瑞利波頻散曲線進行提取，由于瑞利波頻散信息包含于相位譜中，而振幅譜表示的是不同成分信號的能量強度，因此我們對相移法生成的頻散能量譜中的振幅分量進行改造，同時保持相位譜不變，就可以達到改善瑞利波頻散能量收斂性和聚焦性的目的，實現(xiàn)對頻散曲線的準確提取。

根據(jù)上文所述我們將相移法計算得到的頻散能量譜記為Z(v,f)，對Z(v,f)譜上每個能量點的振幅數(shù)值進行冪次方運算，且冪值n為大于1的整數(shù)，則有：

Zn(v,f)=[Z(v,f)]n,n>1

，

(6)

其中Zn(v,f)為冪運算后得到新頻散能量譜，由于經過冪運算，Zn(v,f)中每一點的振幅相對于Z(v,f)都獲得了增大，但是并未改變利用式(4)計算Z(v,f)的φ值，因此求冪只針對振幅譜，不會改變Z(v,f)中的相位信息。與此同時，由于經過n(n>1)次冪運算，Z(v,f)中不同點之間的振幅值差異得到擴大，那么在求得的Zn(v,f)中，瑞利波的頻譜振幅得到相對增強，并且頻散能量條帶的收斂性和聚焦性得到改善，進而使頻散曲線的提取更加精確。

我們以圖1為例具體說明求冪運算的作用，首先提取頻散能量譜上某個頻率截面上的振動曲線，如圖1a為某實際數(shù)據(jù)的頻散能量譜上15 Hz頻率處的振幅譜曲線，由圖中可見該頻率處頻譜能量有多個峰值，且不同峰值之間差異較小，其中在約300 m/s處最大振幅為0.459，其右側的次級振幅為0.276，主次振幅值之比為1.66；圖1b為經過n次冪(n=4)計算后的曲線，由圖可見約300 m/s處的曲線峰值得到相對加強，其他峰值相對減弱，使得主峰值更加突出，處理后的最大振幅值0.045，次級振幅為0.0058，主次振幅值之比為7.76。

a—原始多峰振幅曲線；b—經過4次方冪運算處理后的曲線

2 模擬數(shù)據(jù)分析

我們模擬一組地層的理論數(shù)據(jù)，采用高階交錯網(wǎng)格有限差分和PML吸收邊界的方法[27]對理論模型的彈性波場進行數(shù)值模擬計算，模型地層參數(shù)由橫波速度vs、縱波速度vp、密度ρ和層厚度h組成。具體模型參數(shù)見表1，模型大小設置為80 m×50 m，震源為主頻25 Hz的雷克子波，震源與檢波器最小偏移距為5 m，設置51道檢波器接收，道間距1 m，采樣間隔0.5 ms，采樣長度0.5 s。在模擬數(shù)據(jù)中加入隨機噪聲，此外數(shù)據(jù)中也包含反射波、折射波等干擾波場。

表1 理論模型地層參數(shù)

圖2a為加入白噪聲的模擬數(shù)據(jù)結果，利用相移法和改進的相移法分別對這組數(shù)據(jù)進行處理，從炮集記錄中可以看到較強的噪聲信號以及初至波同相軸，圖2b為常規(guī)相移法計算生成的頻散能量，圖中可見清晰的面波頻散能量條帶，在面波能量條帶兩側也顯示出噪聲及其他類型波場的能量團，根據(jù)給定模型計算出的理論相速度頻散曲線(圖2b中的點狀曲線)，我們可以識別出頻帶范圍較窄的部分高階模式面波的頻散能量，由于煤田地震資料中的面波以能量較強的基階模式為主，因此本文主要針對面波基階模式頻散曲線的提取技術進行研究。

圖2c為利用改進的相移法對模擬數(shù)據(jù)進行處理后生成的頻散能量圖，與圖2b對比可以看出基階模式面波能量更加收斂，頻散譜寬度變窄，在實際數(shù)據(jù)處理中，能夠更加準確地從頻散譜的峰值提取頻散曲線；與此同時，噪聲及其他類型波場(初至波、高階面波等)在能量譜中受到了較好的壓制，使得基階面波能量更易識別。

a—模擬地震數(shù)據(jù)；b—常規(guī)相移法生成的頻散能量；c—改進的相移法生成的頻散能量

3 應用實例

以鄂爾多斯盆地南部黃土塬某研究區(qū)的煤田地震勘探資料為例，利用改進相移法對其進行處理。研究區(qū)表層發(fā)育巨厚、干燥、疏松的黃土層，地震激發(fā)和接收條件較差，導致地震資料干擾波發(fā)育、靜校正問題突出、高頻吸收衰減嚴重、資料分辨率及信噪比較低。一般情況下，黃土層的縱波速度為400～1 800 m/s，橫波速度為200～700 m/s，且速度、密度隨深度增加連續(xù)變化，常規(guī)的反射、折射等勘探方法難以對巨厚黃土層的內部結構進行有效探測，為煤田地震勘探資料的后續(xù)處理帶來較大的困難。

本例中的地震數(shù)據(jù)通過布設396道垂向檢波器組成的接收排列采集，檢波器的自然頻率為10 Hz，道間距10 m，每道采樣點數(shù)為2 001個，采樣間隔為1 ms。激發(fā)采用中間放炮的方式，炮點在排列中不同位置處激發(fā)，本研究利用資料中連續(xù)采集的41炮記錄，其中首炮位于184道檢波器位置處，然后炮點向排列增大方向移動，炮間距30 m(由于研究區(qū)地表條件限制，部分相鄰炮間距為20 m或40 m)。

圖3a為采集的第1炮數(shù)據(jù)的160～360道地震記錄，炮點位于184道位置處，從記錄上能夠看到清晰的面波、反射波、直達波以及折射波信號。我們從中抽取炮點右側25道(186～210道)瑞利波記錄較完整的采集道(圖3b)，對其進行處理。

a—第1炮160～360道記錄；b—第1炮186～210道記錄

3.1 頻散測量

圖4a為利用常規(guī)相移法的生成的能量譜，譜中可見兩個較強的能量團，圖4b為利用改進的相移法生成的能量譜，譜中僅有一個明顯的能量條帶，即為瑞利波頻散能量，相比于常規(guī)相移法生成的結果，改進相移法生成的能量譜中其他類型波場及噪聲能量譜均得到了壓制，且瑞利波的能量條帶更加收斂，峰值得到增強，在提取頻散曲線時，能夠提高各頻點對應的相速度的識別精度，我們對該能量條帶上不同頻率處峰值點的橫縱坐標進行提取，就得到了各個頻率對應的面波相速度值。

圖4 常規(guī)相移法(a)和改進的相移法(b)生成的頻散能量譜

3.2 橫波速度反演

利用改進的相移法對41炮數(shù)據(jù)分別進行處理，提取了不同測點位置處的頻散曲線，然后對頻散曲線進行反演，獲得不同測點下方地層的橫波速度結構。圖5為第17炮數(shù)據(jù)的處理結果，包括從該數(shù)據(jù)中提取的頻散曲線(圖5a中藍色虛線)，以及對該頻散曲線進行反演獲得的測點下方的橫波速度結構(圖5b中紅色折線)。本研究采用阻尼最小二乘方法進行反演，反演的初始模型由頻散曲線估計得到，初始模型見圖5b中的藍色折線所示，最終反演模型對應的頻散曲線見圖5a中紅色實線所示，其與觀測頻散曲線吻合良好。

圖5 第1炮記錄的頻散曲線(a)及第1炮測點下方的地層橫波速度結構(b)

3.3 2D剖面

經過頻散測量和反演處理后，得到各測點的橫波速度結構，然后通過插值將所有測點的橫波速度結構聯(lián)合，生成了研究區(qū)內一條1.2 km長的橫波速度剖面(圖6)。剖面圖顯示了測線下方80 m深度范圍內的地層橫波速度變化情況，地層橫波速度值總體在150～550 m/s范圍內，其中0～15 m深度范圍橫波速度低于250 m/s，主要為松散的土壤沉積物，在大約15～50 m深度范圍，橫波速度為250～350 m/s，在大于50～70 m深度范圍，橫波速度為350～500 m/s，主要成分為經過壓實的較堅硬黃土，在大于70 m深度，地層的橫波速度大于500 m/s，推測為較密實的碎石土沉積。

圖6 研究區(qū)2D橫波速度結構剖面

4 結論

1)本文基于改進的相移法對煤田地震中的瑞利波資料進行處理，在求取頻散能量譜時引入求冪運算，對譜上各數(shù)據(jù)點的振幅進行二次處理，使得瑞利波的能量得到增強，其他類型波場及噪聲的能量受到壓制。經過處理，瑞利波的頻散能量條帶的收斂性和聚焦性顯著改善，有效地提高了提取頻散曲線的精度。

2)利用改進相移法對某研究區(qū)煤田地震資料中的瑞利波進行了處理，獲得研究區(qū)的二維橫波速度剖面，建立了覆蓋層模型，結果表明該方法能為復雜波場條件下的瑞利波頻散測量提供技術支持，在黃土覆蓋區(qū)表層速度建模及近地表成像領域具有良好的應用前景。