基于相位差法和多重濾波法結(jié)合的瑞雷面波頻散曲線(xiàn)研究

諶 強(qiáng)，熊章強(qiáng),2，張大洲,2

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083；.中南大學(xué) 有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410083)

1 引 言

瑞雷面波法作為一種工程勘察手段，由于具有抗干擾能力強(qiáng)、垂向分辨率高和野外數(shù)據(jù)采集相對(duì)簡(jiǎn)便等優(yōu)點(diǎn)廣泛應(yīng)用于淺層地質(zhì)勘查、地質(zhì)災(zāi)害調(diào)查和無(wú)損檢測(cè)等方面[1-3]。瑞雷面波勘探中頻散曲線(xiàn)提取是非常關(guān)鍵的一個(gè)環(huán)節(jié)，其提取的正確與否直接關(guān)系到最終反演結(jié)果的可靠程度，因此如何準(zhǔn)確提取頻散曲線(xiàn)對(duì)于瑞雷面波勘探顯得尤為重要。目前在瑞雷面波勘探頻散曲線(xiàn)提取中大多采用多道瞬態(tài)面波分析法(Multi-channel Analaysis of Surface Wave，簡(jiǎn)稱(chēng)MASW)[4]，這種面波分析方法利用多道疊加效應(yīng)能夠得到較為準(zhǔn)確的頻散曲線(xiàn)，因而在垂向上探測(cè)精度較高。但由于多道疊加效應(yīng)，一般在工程上應(yīng)用12道或24道得到的地震記錄只能計(jì)算出一條頻散曲線(xiàn)，這樣使得MASW方法由于其平均效應(yīng)導(dǎo)致在橫向上分辨率相對(duì)較低，因此只能適用于介質(zhì)速度橫向變化比較小的地層探測(cè)。為了提高M(jìn)ASW方法的橫向分辨率，Koichi Hayashi等人[5]提出并借鑒了地震反射波勘探中共反射點(diǎn)多次覆蓋思想的共中心點(diǎn)(Common Mid-Point,CMP)瞬態(tài)多道面波方法，該方法利用密集的多覆蓋數(shù)據(jù)能獲得相對(duì)于常規(guī)瞬態(tài)面波方法更密集的頻散曲線(xiàn)，從而提高橫向分辨率。胡晉等人[6]通過(guò)CMP解析法再利用相位差法提取頻散提高了頻散曲線(xiàn)的精度。但對(duì)于共中心點(diǎn)疊加數(shù)據(jù)，當(dāng)炮檢距較大時(shí)，所獲得的地層的頻散信息也會(huì)同樣存在平均效應(yīng)。

隨著瑞雷面波應(yīng)用范圍的不斷擴(kuò)大和需要解決的地質(zhì)問(wèn)題復(fù)雜程度的增大，對(duì)于探測(cè)如孤石、巖溶及城市塌陷等地質(zhì)異常體時(shí)需要有較高的橫向分辨率[7-9]，在這種情況下，MASW方法往往得不到較好的效果。面波譜分析法(Spectral Analysis of Surface Wave，簡(jiǎn)稱(chēng)SASW)是通過(guò)計(jì)算不同頻率下的兩道信號(hào)相位差，得到該頻率下的相速度[10]。相較于MASW法，利用SASW方法對(duì)于一個(gè)單炮記錄可提取多條頻散曲線(xiàn)，顯然可以提高瑞雷面波探測(cè)的橫向分辨率[11]。利用兩道地震數(shù)據(jù)提取頻散曲線(xiàn)的方法除SASW外，多重濾波法也是一種常用方法，該方法避免了利用兩道信號(hào)相位差計(jì)算相速度時(shí)相位的不穩(wěn)定性以及相位差需要小于2π的要求。Xiao-Ping Li等人[12]詳細(xì)介紹了多重濾波方法的基本原理。周青云等人[13]證明了多重濾波法在分析面波的頻散特性是一種有效而快速的方法，并研究了多重濾波法的可靠性。陳杰等人[14]利用多重濾波法提取了基階和高階瑞雷面波頻散曲線(xiàn)。

在瑞雷面波法勘探中，理想的情況是盡可能利用少的地震道數(shù)據(jù)提取到較為可靠的頻散曲線(xiàn)，這樣既能夠滿(mǎn)足垂向探測(cè)的需求，在橫向上也有較高的分辨率，達(dá)到能準(zhǔn)確探測(cè)地下地質(zhì)結(jié)構(gòu)的目的[15,16]。但由于目前采用的是瞬態(tài)瑞雷面波勘探方法，震源在激發(fā)后會(huì)產(chǎn)生不同頻率、不同模式的面波信號(hào)，這些不同類(lèi)型疊加的信號(hào)在利用兩道數(shù)據(jù)提取頻散曲線(xiàn)時(shí)會(huì)互相影響，導(dǎo)致所得到的頻散曲線(xiàn)產(chǎn)生較大誤差[17]，這也是目前采用多道數(shù)據(jù)來(lái)計(jì)算頻散曲線(xiàn)的一個(gè)重要原因，即通過(guò)多道數(shù)據(jù)累加來(lái)增強(qiáng)頻散曲線(xiàn)的準(zhǔn)確性和穩(wěn)定性。兩道數(shù)據(jù)計(jì)算頻散曲線(xiàn)的相位差法是通過(guò)計(jì)算兩道的相位差得到不同頻率的相速度，而多重濾波法是通過(guò)求得兩道數(shù)據(jù)不同頻率信號(hào)的時(shí)間差來(lái)計(jì)算其相速度。

本文通過(guò)對(duì)原始地震數(shù)據(jù)進(jìn)行模式分離，然后對(duì)分離后的兩道數(shù)據(jù)分別利用相位差法和多重濾波法計(jì)算出相速度后進(jìn)行疊加，提高了頻散曲線(xiàn)的可靠性，并通過(guò)理論模型和實(shí)測(cè)數(shù)據(jù)說(shuō)明了該方法的有效性。

2 基本原理

2.1 相位差法

假設(shè)地面x1,x2兩點(diǎn)處地震記錄分別為u1(t),u2(t)，其中x2-x1≤λR，λR為瑞雷波波長(zhǎng)(m)，則x1,x2兩點(diǎn)處瑞雷面波的垂直位移方程可分別簡(jiǎn)寫(xiě)為：

(1)

(2)

式中，A1和A2分別為兩點(diǎn)的幅值(m);t為時(shí)間(s);w為角頻率(rad/s);VR為相速度(m/s)。則可得它們之間的相關(guān)系數(shù)可寫(xiě)為：

(3)

式中，τ為u1(t)與u2(t)之間的延遲時(shí)間(s);T為延遲最大時(shí)間(s);對(duì)求出的相關(guān)系數(shù)進(jìn)行傅里葉變換，可得：

(4)

(5)

式中，Δx為兩點(diǎn)距離差(m);通過(guò)公式(5)可計(jì)算出不同頻率下瑞雷面波的相速度VR。

2.2 多重濾波法

多重濾波法將兩道地震數(shù)據(jù)通過(guò)窄帶濾波，利用兩道信號(hào)之間某一頻率下的能量極值對(duì)應(yīng)的時(shí)間差來(lái)計(jì)算相速度。具體在計(jì)算時(shí)，先分別對(duì)兩道地震信號(hào)S(t1)和S(t2)進(jìn)行傅里葉變換至頻率域：

(6)

式中,f為頻率(Hz);t為時(shí)間(s)。接著設(shè)計(jì)帶通濾波器H(wn,w)如公式(7)所示：

(7)

式中,fn為中心頻率(Hz);f1n為上限頻率(Hz)，在本文中取0.75fn;f2n為下限頻率(Hz)，取1.25fn;α為衰減因子，取50.3。在獲得頻率域信號(hào)u(f)和濾波器H(fn,f)后將兩者相乘再做傅里葉反變換到時(shí)間域，如式(8)所示：

(8)

再對(duì)式(8)進(jìn)行希爾伯特變換得到式(9)：

(9)

最后對(duì)其求幅值：

(10)

式中,real表示取復(fù)數(shù)的實(shí)部。對(duì)于兩道數(shù)據(jù)通過(guò)窄帶濾波得到單頻信號(hào)的幅值A(chǔ)n后，分別可以得到兩道記錄能量最大值A(chǔ)1和A2所對(duì)應(yīng)的時(shí)間t1和t2，進(jìn)而計(jì)算出時(shí)間差Δt。利用公式V=Δx/Δt便可以求出對(duì)應(yīng)頻率的速度值，式中Δx為兩道之間的距離(m)。

兩道數(shù)據(jù)分別利用相位差法和多重濾波方法計(jì)算出不同頻率對(duì)應(yīng)的相速度值。為了能夠?qū)崿F(xiàn)將兩種不同方法計(jì)算得到的結(jié)果進(jìn)行疊加，同時(shí)將不同道間距及不同震源位置多次激發(fā)地震波所計(jì)算得到的結(jié)果進(jìn)行累加，本文通過(guò)將給定的一個(gè)速度范圍值再根據(jù)速度間隔進(jìn)行離散化，形成一個(gè)F-V的二維數(shù)組，對(duì)于不同頻率計(jì)算得到的速度值找到對(duì)應(yīng)的數(shù)組并在其中賦值為1，從而形成一個(gè)頻散譜。對(duì)于不同方法、不同道間距及不同炮點(diǎn)所計(jì)算的頻散譜進(jìn)行疊加，最后在頻率方向上做歸一化，這樣就形成了最終的頻散譜，然后從頻散譜中的不同頻率值拾取能量最大位置所對(duì)應(yīng)的速度值就能得到頻散曲線(xiàn)。

3 理論模型測(cè)試

3.1 兩層層狀介質(zhì)模型

為了對(duì)比分析相位差法和多重濾波法相結(jié)合計(jì)算頻散曲線(xiàn)的有效性，特設(shè)計(jì)一個(gè)兩層層狀介質(zhì)模型，模型大小為100 m×50 m。第一層厚度H1=10 m,縱波速度為VP=800 m/s，橫波速度VS=200 m/s，密度ρ=2.0 g/cm3；第二層縱波速度VP=1 200 m/s，橫波速度VS=400 m/s，密度ρ=2.0 g/cm3。采用高精度交錯(cuò)網(wǎng)格有限差分法對(duì)該模型進(jìn)行全波場(chǎng)正演模擬[18,19]，震源為主頻20 Hz的雷克子波，位于地表(0,0)處。模擬得到的地震記錄如圖1所示，地震記錄共90道(繪圖時(shí)以2 m道間距繪制)，道間距為1 m，偏移距10 m，采樣間隔為0.20 ms，采樣點(diǎn)數(shù)為5 001。

從圖1中可以清晰地看出,模擬得到的地震記錄中有反射波、直達(dá)波、折射波和面波等，其中瑞雷面波能量最強(qiáng)，同相軸最明顯。為驗(yàn)證將相位差法和多重濾波方法結(jié)合的有效性，分別對(duì)原始地震記錄以及F-K模式分離后的地震記錄進(jìn)行處理。首先取原始地震記錄的第44道和第46道用相位差法和多重濾波方法分別計(jì)算其頻散譜，如圖2所示。圖2(a)為利用相位差法計(jì)算得到的頻散譜，從圖中可以看出通過(guò)計(jì)算得到的相速度值與理論值部分能夠?qū)?yīng)，但在20～60 Hz段誤差相對(duì)較大，最大誤差達(dá)10 %。存在誤差的主要原因應(yīng)是各模態(tài)面波信號(hào)之間的相互影響所致。張大洲等人[20,21]通過(guò)模式分離方法也對(duì)此進(jìn)行過(guò)同樣的分析。圖2(b)為利用多重濾波法計(jì)算得到的頻散譜，與通過(guò)相位差法計(jì)算的結(jié)果相比較，其誤差相對(duì)較小，最大誤差為8 %。從圖2(a)和圖2(b)的結(jié)果對(duì)比來(lái)看，利用多重濾波方法計(jì)算頻散曲線(xiàn)時(shí)各模式之間的影響相對(duì)較小。將兩種方法計(jì)算得到的頻散譜進(jìn)行疊加得到如圖2(c)所示結(jié)果，從圖中可以看出，將兩種方法計(jì)算得到的頻散譜進(jìn)行疊加后其相速度相同的位置所對(duì)應(yīng)的能量就會(huì)有所加強(qiáng)，從而使得計(jì)算結(jié)果更加穩(wěn)定。

圖1 兩層模型原始地震記錄Fig.1 Original seismic record of two-layer model

圖2 兩層模型模擬數(shù)據(jù)第44道和第46道不同方法計(jì)算頻散曲線(xiàn)對(duì)比Fig.2 Comparison of the dispersion curves of the 44th and 46th channels of the two-layer model calculated by different methods

為避免反射波、直達(dá)波和折射波等干擾波以及瑞雷面波各模式間對(duì)計(jì)算結(jié)果的影響，因此在計(jì)算之前需要F-K進(jìn)行模式分離。由于采集的數(shù)據(jù)為等間距線(xiàn)性排列，因此可以在頻率-波數(shù)(F-K)域進(jìn)行模式分離。圖3(a)為將圖1所示的地震記錄轉(zhuǎn)換到F-K的結(jié)果，從圖中可看出不同模式的瑞雷波相互分離。在F-K域中將其他模式切除，僅保留基階模態(tài)，如圖3(b)所示。將F-K域中的基階模態(tài)數(shù)據(jù)進(jìn)行傅里葉反變換，便可得到如圖4所示的時(shí)間-空間域基階瑞雷波數(shù)據(jù)。

圖3 兩層介質(zhì)模型頻率波數(shù)譜Fig.3 Frequency wavenumber spectrum of two-layer model

圖4 F-K模式分離后的兩層模型基階地震記錄Fig.4 Fundamental model seismic records of two-layer model after F-K mode separation

圖5 不同方法計(jì)算兩層模型F-K分離后的第44道和第46道的頻散曲線(xiàn)對(duì)比Fig.5 Comparison of the dispersion curvesof channels 44 and 46 after the F-K separation of the two-layer model calculated by different methods

對(duì)于F-K模態(tài)分離后的基階數(shù)據(jù)，選取第44道和第46道采用相位差法和多重濾波法計(jì)算頻散譜，結(jié)果如圖5所示，圖5(a)～圖5(c)分別依次代表相位差法、多重濾波法和兩種方法計(jì)算結(jié)果疊加的頻散譜與理論值的對(duì)比圖。從圖中可見(jiàn)，利用模式分離后的基階地震波數(shù)據(jù)分別用兩種方法計(jì)算得到的頻散曲線(xiàn)與理論值吻合程度均較好。在8～20 Hz頻率段，相位差法與理論值之間的誤差最大值僅為3%，多重濾波法的誤差最大值僅為2.4%，而兩者疊加后的誤差最大值僅為1.8%；頻率在21～60 Hz時(shí)，三者的頻散結(jié)果基本與理論值重合。通過(guò)將F-K分離后的數(shù)據(jù)采用兩種方法計(jì)算的疊加譜與理論值對(duì)比可知，相對(duì)于單一方法，前者可以提供更多的信息。而通過(guò)提高面波頻散曲線(xiàn)提取的準(zhǔn)確度，也能說(shuō)明將兩種數(shù)據(jù)計(jì)算的頻散信息疊加是可行的。

3.2 斜坡模型

為了進(jìn)一步驗(yàn)證通過(guò)相位差法和多重濾波法相互結(jié)合計(jì)算頻散曲線(xiàn)來(lái)提高瑞雷面波探測(cè)的橫向分辨率的情況，特設(shè)計(jì)如圖6所示的斜坡模型。模型大小為90 m×50 m，左端厚度5 m,右端厚度10 m。第一層縱波速度為VP=800 m/s，橫波速度VS=200 m/s，密度ρ=2.0 g/cm3；第二層縱波速度VP=1 200 m/s，橫波速度VS=400 m/s，密度ρ=2.0 g/cm3。數(shù)值模擬采用高精度交錯(cuò)網(wǎng)格有限差分法，震源為主頻25 Hz的雷克子波，分別位于(4,0)、(9,0)、(14,0)、(76,0)、(81,0)和(86,0)處。檢波器位于地表，從19 m到71 m處；道間距為1 m，共有53道。六個(gè)震源位置模擬得到的地震記錄選取位于(4,0)處的地震記錄如圖7所示。

圖6 斜坡模型和參數(shù)Fig.6 Slope model and parameters

圖7 斜坡模型原始地震記錄(震源位于地表4 m處)Fig.7 Original seismic record of slope model(the seismic focus is located at 4 m on the ground)

對(duì)于上述六個(gè)不同震源位置獲得的地震記錄進(jìn)行F-K模式分離得到基階面波，然后利用多道面波分析方法(MASW)計(jì)算得到的頻散曲線(xiàn)如圖8所示，從圖8中可見(jiàn)，盡管震源位置不同，但計(jì)算得到的頻散曲線(xiàn)基本吻合，這也進(jìn)一步說(shuō)明了模擬結(jié)果的正確性。

從另一方面來(lái)看，多道面波方法雖然提取的頻散曲線(xiàn)的準(zhǔn)確度較高，但橫向分辨率較低。接下來(lái)采用相位差法和多重濾波法結(jié)合的方式，利用兩道數(shù)據(jù)計(jì)算出頻散曲線(xiàn)譜。為了使計(jì)算的頻散譜更加穩(wěn)定，在計(jì)算過(guò)程中采用三道數(shù)據(jù)兩兩計(jì)算后再平均，這樣依次向后滑動(dòng)每三道就可以提取到51條頻散曲線(xiàn)。由于本文主要研究的是頻散曲線(xiàn)的提取，因此沒(méi)有對(duì)頻散曲線(xiàn)進(jìn)行反演，避免由于反演帶來(lái)相關(guān)的誤差而影響對(duì)結(jié)果的分析。最后將頻散曲線(xiàn)按照半波長(zhǎng)法進(jìn)行轉(zhuǎn)換得到如圖9所示的相速度擬二維剖面圖。為了進(jìn)一步說(shuō)明該方法的有效性，同樣采用多道(滾動(dòng)式)的方法對(duì)該模型數(shù)據(jù)進(jìn)行處理，結(jié)果如圖10所示。從圖9和圖10中均可以看出相速度在橫向上的變化與傾斜界面的傾向是一致的，并能分辨出界面在橫向上的變化。另外，從圖10中還可以看出，采用多道滾動(dòng)式處理的界面更加平滑，尤其在淺層內(nèi)更為明顯。但這種方法由于不斷移動(dòng)排列會(huì)大大延長(zhǎng)工作時(shí)間，耗費(fèi)更多的人力和物力。因此相對(duì)于采用常規(guī)MASW方法僅得一條頻散曲線(xiàn)而言，利用相位差法和多重濾波法相結(jié)合計(jì)算兩道的方式，不僅可以計(jì)算出高精度的頻散曲線(xiàn)，而且只用一個(gè)排列便可以提取多條頻散曲線(xiàn)，大大節(jié)省了人力和物力，具有良好的實(shí)用性。

圖9 疊加處理的斜坡模型面波速度剖面Fig.9 Surface wave velocity profile of slope model with superposition processing

圖10 多道滾動(dòng)式分析方法斜坡模型面波速度剖面Fig.10 Surface wave velocity profile of slope modelunder multi-channel rolling analysis method

4 混凝土球模型測(cè)試

為測(cè)試本文所提出方法的實(shí)用性，對(duì)某實(shí)驗(yàn)場(chǎng)地埋設(shè)的混凝土球模型采集面波數(shù)據(jù)進(jìn)行處理分析。該模型為直徑1 m的混凝土球，埋深3.5 m，如圖11所示。混凝土周?chē)鸀槊軐?shí)的土層，面波數(shù)據(jù)采集現(xiàn)場(chǎng)布置的排列共有24道，道間距為0.5 m，混凝土球位于排列中間第12個(gè)檢波器的下方。震源方式為錘擊，分別在兩端6 m和8 m偏移距處用18磅大錘激發(fā)，圖12為在6 m偏移距下的實(shí)際地震記錄，從圖中可以看出面波能量較強(qiáng)。

對(duì)于采集得到的地震數(shù)據(jù)通過(guò)模式分離將干擾波去除，然后采用相位差法和多重濾波法計(jì)算其頻散譜并進(jìn)行疊加。為驗(yàn)證該方法在不同偏移距下也能夠進(jìn)行疊加處理，分別對(duì)同一排列兩側(cè)的6 m和8 m處的偏移距提取同一位置的頻散曲線(xiàn)，并對(duì)其進(jìn)行疊加，結(jié)果如圖13所示。從圖中可見(jiàn)，在同一位置處盡管偏移距不同，但是提取的頻散曲線(xiàn)基本一致。相比于疊加處理結(jié)果，兩處偏移距下的最大誤差僅4.3 %，證明了疊加的可行性，并且發(fā)現(xiàn)進(jìn)行疊加處理后的頻散曲線(xiàn)更加穩(wěn)定。

圖13 不同偏移距的頻散曲線(xiàn)Fig.13 Dispersion curves with different offset distances

圖14為計(jì)算得到的該排列的23條頻散曲線(xiàn)空間分布圖，頻散曲線(xiàn)間隔為0.5 m。對(duì)于計(jì)算得到的頻散曲線(xiàn)同樣利用半波長(zhǎng)法轉(zhuǎn)換得到如圖15所示的擬二維相速度剖面，從圖中可以清晰地看出在水平位置5.8 m、深度3.0 m處相速度值明顯高于周?chē)?圖中用黑色虛線(xiàn)圈定位置)，這一異常與混凝土球在土層中表現(xiàn)為相對(duì)高速特征相吻合。異常在水平位置與埋設(shè)的混凝土球?qū)嶋H位置一致，僅在深度上略有差異，這主要是由于半波長(zhǎng)法的近似性所導(dǎo)致?？傊摲椒ㄋ崛〉念l散曲線(xiàn)較好，能夠比較準(zhǔn)確地反映地下地質(zhì)異常體特征，在水平方向有較高的分辨率。另一方面，利用該方法進(jìn)行探測(cè)，一條測(cè)線(xiàn)就可以得到一個(gè)二維斷面，其探測(cè)效率大大提高。

圖14 23條頻散曲線(xiàn)的空間分布Fig.14 Spatial distribution of the 23 dispersion curves

圖15 混凝土球模型擬二維相速度剖面Fig.15 Proposed two-dimensional phase velocity profile for concrete ball model

5 結(jié) 論

通過(guò)采用相位差法和多重濾波法相結(jié)合的方式針對(duì)不同地質(zhì)模型以及實(shí)驗(yàn)場(chǎng)地的面波數(shù)據(jù)提取頻散曲線(xiàn)，可得出以下結(jié)論：

1)利用相位差法和多重濾波法都可以提取頻散曲線(xiàn)，而前者由于相位變化容易受到不同模式的干擾，在利用兩道數(shù)據(jù)計(jì)算頻散曲線(xiàn)時(shí)，多重濾波法相對(duì)于相位差法的抗干擾能力更強(qiáng)。

2)瑞雷面波探測(cè)中可采用多震源激發(fā)瑞雷面波，采用相位差法和多重濾波法相結(jié)合的方式對(duì)不同震源激發(fā)的數(shù)據(jù)所提取的頻散曲線(xiàn)進(jìn)行疊加，可提高頻散曲線(xiàn)的準(zhǔn)確性和可靠性。

3)利用相位差法和多重濾波法相結(jié)合的方式，對(duì)于單個(gè)排列的數(shù)據(jù)可提取到多條頻散曲線(xiàn)，既能提高探測(cè)的橫向分辨率，也能提高其探測(cè)效率。