周 曉,喻 顏,阮楊帆,岳子沖,萬幸源,劉宏岳

(1.武漢理工大學(xué) 機(jī)電工程學(xué)院,湖北 武漢430070;2.福建省建筑設(shè)計研究院,福建 福州 350001)

地球表面無時無刻都存在一種人類無法感知的微弱振動。這些微弱振動產(chǎn)生的信號叫微動信號[1]。微動勘探技術(shù)就是通過微震臺陣觀測,以平穩(wěn)隨機(jī)過程理論為依據(jù),從微動信號中提取面波(Rayleigh波)頻散曲線,通過對頻散曲線的反演,從而推斷地下結(jié)構(gòu)的一種方法[2]。微動勘探技術(shù)以其無需人工源、抗干擾能力強(qiáng)、分辨率高、便捷環(huán)保等優(yōu)勢,在城市建設(shè)和地下空間開發(fā)等領(lǐng)域發(fā)揮重要作用。

在工程微動勘探技術(shù)中提取面波頻散信息常用的方法有空間自相關(guān)法(spatial auto-correlation,SPAC)[3]和頻率-波數(shù)法(frequency wavenumber,F-K)[4]。地震干涉法[5](seismic interferometry,SI)是近年來地球物理研究的一個熱點,其理論提出較早,但是并不統(tǒng)一,后續(xù)學(xué)者有從波動互易、模式均分、時間反轉(zhuǎn)和穩(wěn)相近似等不同角度對地震干涉法進(jìn)行了證明和解釋[6]。目前有學(xué)者已經(jīng)證明SI法與SPAC法其物理本質(zhì)是一致的[7-8],SPAC法是在時間域中的描述,而SI法是在頻域中的描述。

這3種頻散曲線的提取方法各有優(yōu)劣。SPAC法在數(shù)據(jù)分析上較為簡單,但是對于臺陣布置要求嚴(yán)格不能設(shè)置為線性排列系統(tǒng),常用多重圓形臺陣,生產(chǎn)作業(yè)中受到地理條件的制約較大[9]。F-K法對于臺陣擺放要求低,但F-K法提取的頻散曲線反映的是臺陣下方區(qū)域的平均效應(yīng),在同樣的測量精度下,需要使用更多的檢波器。SI法計算的是兩個檢波器路徑間的平均效應(yīng),在同等數(shù)量檢波器下,SI法可以得到地表下方地質(zhì)結(jié)構(gòu)的更多頻散信息,覆蓋的勘探區(qū)域更加全面。對于地下結(jié)構(gòu)的二維剖面探測,SI法只用沿測線將檢波器布置成直線陣列,而SPAC法和F-K法往往需要沿測線布置多個檢測臺陣或者采用沿測線平移臺陣多次測量的方式[10]。因此,SI法可以使用更少的檢波器和更短的時間完成勘探任務(wù),極大增加了工作效率,應(yīng)用前景更加廣泛。

筆者研究頻率域下兩兩臺站微動信號間求解面波相速度的地震干涉算法,針對傳統(tǒng)的貝塞爾函數(shù)擬合一般只使用單調(diào)遞減部分的相干系數(shù)[11],提出一種基于漢克爾變換分段的方法,能夠有效地把相干系數(shù)劃分為不同的單調(diào)區(qū)間段,分段擬合貝塞爾函數(shù)求解面波相速度,解決相速度頻散曲線高頻范圍缺失、探測盲區(qū)深度較大的問題。

1 方法原理

1.1 地震干涉算法原理

干涉在物理學(xué)中指的是兩列及其以上的波空間中相遇發(fā)生疊加形成新的波。地震干涉法是基于兩個檢波器記錄的微動信號,其核心思想是通過互相關(guān)理論計算出兩個信號間的經(jīng)驗格林函數(shù)。當(dāng)噪聲源隨機(jī)分布時,表征噪聲源與波場關(guān)系的經(jīng)驗格林函數(shù)符合真實的格林函數(shù)。

對于單一經(jīng)驗格林函數(shù),可以通過地震干涉結(jié)果頻譜實部擬合零階第一類貝塞爾函數(shù),來求解Rayleigh波相速度頻散信息[12]。

根據(jù)Yokoi等[7]在證明SPAC法與SI法的一致性中推導(dǎo)的理論公式,地震干涉法兩兩之間的互相干函數(shù)滿足:

γA,B=J0(krA,B)

(1)

式中:γA,B為A、B臺站垂直方向微動信號的互相干函數(shù);r為A、B臺站之間的距離;k為波數(shù);J0為第一類零階貝塞爾函數(shù)。

在頻率域中,兩道信號地震干涉法的互相干計算公式為:

(2)

式中:uA、uB為A、B兩點信號的傅里葉變換;*為復(fù)數(shù)共軛;|·|為復(fù)數(shù)的模。

結(jié)合式(1)、式(2),有:

(3)

根據(jù)式(3)相干系數(shù)擬合貝塞爾函數(shù),可求出頻率為f時,對應(yīng)的波數(shù)k,相速度滿足c=2πf/k,即可求出兩道微動信號的瑞雷波相速度頻散曲線。上述公式推導(dǎo)適用于以基階瑞雷波能量為主,利用垂直分量觀測數(shù)據(jù),提取面波頻散曲線的情況。對于水平分量觀測數(shù)據(jù)則不適用。

1.2 多段漢克爾變換分段原理

漢克爾變換[13]是一種積分變換,又稱作傅里葉-貝塞爾變換,是由數(shù)學(xué)家Hankel提出的。對給定核函數(shù)f(r),以v階第一類貝塞爾函數(shù)Jv(kr)作無窮級數(shù)展開,級數(shù)的各項k變化,各項Jv(kr)前的系數(shù)Fv構(gòu)成了變換函數(shù),就是一種v階漢克爾變換。k為自變量,把連續(xù)函數(shù)Fv(k)在(0,∞)上表示為:

(4)

地震干涉法在估算瑞雷波相速度頻散曲線時,傳統(tǒng)的貝塞爾函數(shù)擬合方法只使用了單調(diào)遞減部分相干系數(shù),而震蕩衰減部分的相干系數(shù)沒有被利用,對應(yīng)的有效頻率范圍集中在低頻區(qū)間,高頻數(shù)據(jù)缺失,反映在深度上即淺層數(shù)據(jù)缺失,探測盲區(qū)大。筆者利用漢克爾變換,把地震干涉法提取的相干系數(shù)γ(f)作為核函數(shù),其零階第一類貝塞爾函數(shù)的漢克爾變換為:

(5)

當(dāng)k>ks(ks為控制參數(shù),通常取ks=0.15),并且F(k)取得最大值時,對應(yīng)的k=ξ,在滿足這個條件下,擬合的貝塞爾函數(shù)效果最佳,對應(yīng)擬合的貝塞爾函數(shù)表達(dá)式為:

(6)

圖1 漢克爾變換擬合貝塞爾曲線

通過相干系數(shù)的極點劃分單調(diào)區(qū)間段,容易受到局部特殊點的影響,造成分段錯誤。而漢克爾變換利用了所有頻點進(jìn)行積分運算,不受特殊點的影響,其結(jié)果更具有魯棒性。

圖2 多段漢克爾變換擬合貝塞爾曲線

將不同區(qū)間段的相干系數(shù)分別與圖3所示的標(biāo)準(zhǔn)零階第一類貝塞爾函數(shù)對應(yīng)段擬合,求解相速度。

圖3 標(biāo)準(zhǔn)零階第一類貝塞爾函數(shù)

2 微動數(shù)據(jù)采集及其處理

2.1 微動數(shù)據(jù)采集

圖4為本次實驗用到的13個檢波器,由實驗室自主研發(fā)。檢波器收集3個分量的微動信號,存儲在SD(secure digital)卡中,采集完成后通過WiFi將數(shù)據(jù)導(dǎo)出。

圖4 微動智能勘探儀檢波器

本次實驗設(shè)計了直線型陣列,一共有13個測點,相鄰測點間隔5 m,布置在馬路旁。微動觀測陣列布置如圖5所示。

圖5 觀測陣列

實驗選擇在工作日的下午兩點至四點進(jìn)行,記錄采集時長60 min,采樣頻率250 Hz,測得的13個檢波器垂直分量原始波形如圖6所示。

圖6 微動波形

2.2 相干系數(shù)計算

通過式(2)計算10 m間距兩兩檢波器微動數(shù)據(jù)之間的相干系數(shù),結(jié)果如圖7所示。從曲線形態(tài)上分析,11組數(shù)據(jù)具有良好的一致性,特別是在1.5～22 Hz區(qū)間段,曲線形狀接近,并且近似滿足貝塞爾函數(shù)的形狀,包含單調(diào)遞減部分和震蕩衰減部分。從數(shù)據(jù)有效性的角度分析,1.5 Hz以下和40 Hz以上的部分相干系數(shù)震蕩幅度大,數(shù)據(jù)的可靠性較低。后續(xù)進(jìn)行相干系數(shù)分段,貝塞爾函數(shù)擬合時,需要把這些不可靠的數(shù)據(jù)剔除。計算選取的有效數(shù)據(jù)段是2～40 Hz。

圖7 相干系數(shù)

2.3 相干系數(shù)分段

圖8 S1-S3擬合貝塞爾曲線

對10 m距離的11道相干系數(shù)進(jìn)行漢克爾變換,求得各段的最優(yōu)解如表1所示。

表1 漢克爾變化的最優(yōu)解

利用式(6)擬合相干系數(shù)曲線,得到擬合的貝塞爾函數(shù),基于擬合的貝塞爾函數(shù)確定的分段點信息如表2所示。

表2 基于多段漢克爾變換的分段區(qū)間

2.4 頻散曲線提取

利用表2所示信息,把相干系數(shù)分成4段,使用分段后的數(shù)據(jù)分別與標(biāo)準(zhǔn)貝塞爾函數(shù)對應(yīng)單調(diào)段擬合,求解相速度,結(jié)果如圖9所示。

圖9 多段漢克爾變換后的相速度頻散曲線

圖10為傳統(tǒng)的地震干涉算法處理得到的相速度頻散曲線。從圖10可知,使用漢克爾變換分段擬合的方法,將頻率范圍從8 Hz左右擴(kuò)展到32 Hz左右,并且擴(kuò)展頻散點的相速度在200 m/s范圍以內(nèi)。

圖10 單調(diào)遞減段的相速度頻散曲線

根據(jù)半波長深度定理,可以簡單估算面波穿透深度。使用漢克爾變換分段擬合計算的頻散曲線的半波長深度曲線如圖11所示,探測深度范圍約為2.5～110 m。傳統(tǒng)方法的半波長深度曲線如圖12所示,高頻信息缺失,探測的深度范圍約為10～110 m。

圖11 多段漢克爾變換后的半波長深度曲線

圖12 單調(diào)遞減段的半波長深度曲線

傳統(tǒng)方法的半波長曲線揭示的最小探測深度約為10 m,本文的基于漢克爾變換分段擬合方法得到的半波長曲線揭示的最小探測深度約為2.5 m。本文提出的方法探測的最淺深度更小,對應(yīng)的地下淺層面波結(jié)構(gòu)細(xì)節(jié)更多。

3 勘探結(jié)果分析

使用相速度頻散曲線,經(jīng)插值運算繪制面波相速度等值線圖,可以作為地質(zhì)解釋的基本依據(jù)。圖13為使用單調(diào)遞減段擬合求得的相速度繪制的面波相速度等值線圖。圖14為基于漢克爾變換分段擬合求得的相速度繪制的面波相速度等值線圖。圖中顏色代表面波相速度,速度大小如色標(biāo)所示。

圖13 單調(diào)遞減段擬合的剖面圖

圖14 基于漢克爾變換擬合的剖面圖

使用傳統(tǒng)的單調(diào)遞減段擬合的地震干涉法和基于多段漢克爾變換的分段擬合的地震干涉法,得到的面波相速度等值線圖,在大于12 m深度的區(qū)域基本一致,但是傳統(tǒng)方法在淺部的探測盲區(qū)深度約為12 m,改進(jìn)的方法探測盲區(qū)約為2.5 m。

4 結(jié)論

筆者針對SI算法求取面波相速度高頻數(shù)據(jù)缺失,探測盲區(qū)大的問題,提出了基于漢克爾變換的分段擬合方法。將相干系數(shù)根據(jù)零點初步分段,然后在每一段使用漢克爾變換,確定擬合的分段區(qū)間,最后與零階貝塞爾函數(shù)擬合求解相速度。實驗表明,本方法可以有效擴(kuò)展高頻的頻散數(shù)據(jù),減小臺陣的探測盲區(qū)深度,探測能力可以達(dá)到臺站距離的1/4。