基于多段漢克爾變換擬合的地震干涉算法研究

周 曉,喻 顏,阮楊帆,岳子沖,萬(wàn)幸源,劉宏岳

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

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

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

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

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

1 方法原理

1.1 地震干涉算法原理

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

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

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

γA,B=J0(krA,B)

(1)

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

在頻率域中,兩道信號(hào)地震干涉法的互相干計(jì)算公式為:

(2)

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

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

(3)

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

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

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

(4)

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

(5)

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

(6)

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

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

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

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

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

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

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

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

圖4 微動(dòng)智能勘探儀檢波器

本次實(shí)驗(yàn)設(shè)計(jì)了直線型陣列,一共有13個(gè)測(cè)點(diǎn),相鄰測(cè)點(diǎn)間隔5 m,布置在馬路旁。微動(dòng)觀測(cè)陣列布置如圖5所示。

圖5 觀測(cè)陣列

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

圖6 微動(dòng)波形

2.2 相干系數(shù)計(jì)算

通過(guò)式(2)計(jì)算10 m間距兩兩檢波器微動(dòng)數(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í),需要把這些不可靠的數(shù)據(jù)剔除。計(jì)算選取的有效數(shù)據(jù)段是2～40 Hz。

圖7 相干系數(shù)

2.3 相干系數(shù)分段

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

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

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

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

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

2.4 頻散曲線提取

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

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

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

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

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

圖11 多段漢克爾變換后的半波長(zhǎng)深度曲線

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

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

3 勘探結(jié)果分析

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

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

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

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

4 結(jié)論

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