楊潤(rùn)海，譚俊卿，向 涯，姜金鐘，王 彬

(1.云南省地震局，云南 昆明 650224；2.云南大學(xué)，云南 昆明 650091；3.中國(guó)地震局地震研究所 中國(guó)地震局地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430071)

0 引言

地震波是照亮地球內(nèi)部的一盞明燈，隨著地震觀測(cè)和數(shù)字信號(hào)處理技術(shù)的發(fā)展，地震波不僅可以提供地球內(nèi)部的結(jié)構(gòu)信息，其波速變化還能反映地球內(nèi)部的介質(zhì)變化信息(陳颙，朱日祥，2005；陳颙等，2007a，b；丘學(xué)林等，2007；Chenetal，2008，2010)。而其中，基于背景噪聲格林函數(shù)的波速變化和基于精確定位的相似地震研究是了解介質(zhì)性質(zhì)變化的重要方法之一(王偉濤，王寶善，2011；Xu，Song，2009)。

為了通過(guò)地震波了解地殼介質(zhì)性質(zhì)變化，在云南、甘肅省和新疆地區(qū)建立了大型陸地水體氣槍震源觀測(cè)基地(陳颙等，2017；王偉濤等，2017)，并積累了多年的觀測(cè)數(shù)據(jù)。如何精確測(cè)量地震波速變化就成為了解地殼介質(zhì)性質(zhì)變化及其過(guò)程的關(guān)鍵問(wèn)題之一(王彬等，2012；王寶善等，2011，2016；Chenetal，2014；Liuetal，2019；Suetal，2019；Wang，Chen，2019)。

地震波速變化測(cè)量本質(zhì)上是要解決數(shù)字信號(hào)處理的時(shí)延估計(jì)問(wèn)題，關(guān)鍵是要提高時(shí)延估計(jì)的精度(劉自鳳等，2015；葉泵等，2017；欒奕等，2016；向涯等，2017；林建民等，2010；Knapp，Carter，2003)。對(duì)于確定性的周期或非周期信號(hào)、或平穩(wěn)隨機(jī)噪聲，增加取樣長(zhǎng)度有利于提高時(shí)延估計(jì)的精度；但是對(duì)于非平穩(wěn)隨機(jī)信號(hào)，增加取樣長(zhǎng)度將帶來(lái)新的誤差。盡管大容量氣槍震源具有較高的重復(fù)性，而實(shí)際地震記錄信號(hào)由于存在隨機(jī)干擾和能量衰減，具有一定的非平穩(wěn)性。同時(shí)氣槍震源能量有限，傳播距離較遠(yuǎn)時(shí)信號(hào)信噪比低，而傳播距離較近的臺(tái)站震相波列發(fā)育較短，有效信號(hào)較短，這些因素決定了進(jìn)行相關(guān)計(jì)算時(shí)只能進(jìn)行較短取樣(蔣生淼等，2017；李孝賓等，2016，2017；徐逸鶴等，2016)。

鑒于互相關(guān)的時(shí)延檢測(cè)方法受相關(guān)窗口長(zhǎng)度和信號(hào)信噪比影響較大，本文將改進(jìn)的線性調(diào)頻Z變換(Modified Chirp Z Transform，簡(jiǎn)稱MCZT)譜細(xì)化和相關(guān)峰精確插值(Fine Interpolation of Correlation Peak，簡(jiǎn)稱FICP)引入到陸地大容量氣槍主動(dòng)源觀測(cè)中的波速干涉測(cè)量中，介紹了一種能夠獲得高精度走時(shí)差的短時(shí)窗時(shí)延估計(jì)方法。

1 氣槍源地震波速干涉測(cè)量方法

陸地大型水體大容量氣槍的能量釋放由壓力脈沖和氣泡脈沖組成，震源持續(xù)時(shí)間約為2 s，不是典型的脈沖源，同時(shí)存在水位和水深等影響因素。為了消除或減弱這些影響因素，通常對(duì)接收臺(tái)站記錄與源附近的臺(tái)站(參考臺(tái))記錄做反褶積運(yùn)算(王寶善等，2012；翟秋實(shí)等，2016；陳佳等，2017)，得到接收臺(tái)近似于脈沖震源的記錄，即格林函數(shù)。

目前陸地大型水體大容量氣槍的波速測(cè)量，多數(shù)采用處理流程和方法如下(王寶善等，2012；Wangetal，2019)：

(1)參考格林函數(shù)選?。簩⒔邮张_(tái)站的格林函數(shù)疊加得到參考格林函數(shù)或者選取某一次激發(fā)的臺(tái)站格林函數(shù)作為參考函數(shù)。

(2)當(dāng)前格林函數(shù)選?。簩⒚恳淮渭ぐl(fā)的格林函數(shù)或者相鄰幾次激發(fā)的格林函數(shù)疊加形成信噪比更高的格林函數(shù)，作為該段時(shí)間(內(nèi))的格林函數(shù)。

(3)震相選?。焊鶕?jù)震源距和傳播特征選取具有明確物理意義的震相，用于干涉測(cè)量，通常選取P波和S波震相。

(4)干涉測(cè)量：將當(dāng)前格林函數(shù)與參考格林函數(shù)的震相信號(hào)做互相關(guān)進(jìn)行時(shí)延估計(jì)，得到時(shí)延量(即走時(shí)差)，然后將走時(shí)差換算成波速變化率。也可用伸縮法(行鴻彥，唐娟，2008)直接得到波速變化率。

(5)波速變化率：根據(jù)震相走時(shí)將走時(shí)差轉(zhuǎn)換為波速變化率。

地震波速變化測(cè)量常用方法有互相關(guān)法、伸縮法和尾波干涉法(王寶善等，2016；Alexandre，2014)。這些方法都是基于互相關(guān)運(yùn)算，相關(guān)運(yùn)算可以在時(shí)域或頻域進(jìn)行，本文主要討論在頻域計(jì)算中如何提高計(jì)算精度。而在頻域進(jìn)行互相關(guān)計(jì)算時(shí)，需要較精確的頻譜估計(jì)。

標(biāo)準(zhǔn)快速傅里葉變換(Fast Fourier Transform，簡(jiǎn)稱FFT)的計(jì)算速度快，但對(duì)于實(shí)數(shù)序列的FFT計(jì)算存在2個(gè)顯著的缺陷：一是在計(jì)算結(jié)果中，序列左半邊代表信號(hào)的頻譜，右半邊是左半邊的共扼對(duì)稱，只有一半是有用信息，浪費(fèi)了另一半；二是FFT計(jì)算出的信號(hào)的頻譜相當(dāng)于對(duì)頻譜的離散抽樣，存在柵欄效應(yīng)。譜分辨率受取樣長(zhǎng)度tp制約，Δf=1/tp。在滿足信號(hào)平穩(wěn)性的前提下，為了減小柵欄效應(yīng)的影響，可以增加取樣長(zhǎng)度，但這也增加了計(jì)算量。另外，標(biāo)準(zhǔn)FFT計(jì)算的相關(guān)函數(shù)，是從N點(diǎn)序列到N點(diǎn)序列的計(jì)算，相關(guān)函數(shù)序列的點(diǎn)的時(shí)間間隔與采樣周期一致，會(huì)降低相關(guān)波形的分辨率(間隔增加)，難以得到精確的時(shí)延。針對(duì)此問(wèn)題，傳統(tǒng)的方法是對(duì)相關(guān)峰進(jìn)行插值計(jì)算找峰值點(diǎn)，因插值結(jié)果不能完全反映相關(guān)峰的形狀，因此，插值計(jì)算會(huì)帶來(lái)新的誤差，且增加了運(yùn)算量(楊亦春等，2002)。綜上認(rèn)為，精確的互譜估計(jì)和互譜插值對(duì)時(shí)延的精確估計(jì)很重要。

在大容量氣槍源觀測(cè)中，地震儀的采樣率一般為100 Hz，實(shí)際時(shí)延量遠(yuǎn)小于一個(gè)采樣間隔，互相關(guān)得到的時(shí)延量只有一個(gè)采樣間隔的分辨率。為了提高分辨率，對(duì)互相關(guān)函數(shù)的相關(guān)峰通過(guò)余弦插值等方法進(jìn)行插值，可以得到更精確的時(shí)延估計(jì)值，取得較好的結(jié)果(Wangetal，2019)。但余弦插值計(jì)算過(guò)程與信號(hào)的信息無(wú)關(guān)，即與互譜無(wú)關(guān)，插值結(jié)果只能近似反映相關(guān)峰的形狀。因此要達(dá)到計(jì)算精度更高而計(jì)算量小的目的，有必要尋求同時(shí)提高互譜和相關(guān)峰計(jì)算精度的方法。為此，本文探討利用MCZT譜細(xì)化和FCIP在陸地大容量氣槍主動(dòng)源觀測(cè)中的波速干涉測(cè)量的應(yīng)用，并與余弦插值方法進(jìn)行了對(duì)比。

2 MCZT譜細(xì)化

頻譜的細(xì)化計(jì)算常見(jiàn)的方法有ZOOM算法、頻率抽取法、降采樣法、ZFFT(降采樣)、cFFT(級(jí)聯(lián)FF)等(楊亦春等，2002；侯朝煥等，1990)，這幾種方法的共同特點(diǎn)是計(jì)算量大。MCZT算法可在不增加計(jì)算量的前提下提高頻譜計(jì)算精度(馬少春等，2014；楊亦春等，2003；唐娟，行鴻彥，2007；杜娟，程擂，2010；王莉等，2011；徐世友，2003；Benestyetal，2004)。

對(duì)長(zhǎng)度為N的時(shí)域序列x1(n)和x2(n)，其相關(guān)函數(shù)為：

R(m)=E[x1(n)·x2(n+m)]

(1)

式中：R表示相關(guān)函數(shù)；m表示時(shí)延；E表示數(shù)學(xué)期望；n表示時(shí)間。

相關(guān)函數(shù)的主峰對(duì)應(yīng)的時(shí)間就是信號(hào)x1(n)和x2(n)之間的時(shí)延量，R(n)的頻譜R(k)與x1(n)和x2(n)的頻譜X1(k)和X2(k)的關(guān)系為：

R(k)=X1(k)·X*2(k)

(2)

故頻譜計(jì)算精度直接影響相關(guān)函數(shù)的精度。

MCZT定義為：

(k= 0，1，…，N-1)

(3)

定義MCZT的逆變換為IMCZT：

(4)

由式(3)(4)可看出，x′(n)不等于x(n)，IMCZT與MCZT的計(jì)算結(jié)果在數(shù)值上是共軛關(guān)系，MCZT是CZT(線性調(diào)頻Z變換)的一個(gè)特例。頻譜的間隔由N1決定，Δf=fs/N，fs為信號(hào)采樣率，當(dāng)N1=fs時(shí)，頻譜的分辨率就達(dá)到1 Hz，N1的選擇不受N的限制，可以采用FFT實(shí)現(xiàn)快速計(jì)算。采用MCZT計(jì)算出的信號(hào)頻譜，是從0 Hz開(kāi)始到fmax=fs*(N/N1)Hz，即是截取了一段低頻譜，為保證不影響相關(guān)精度，應(yīng)該覆蓋信號(hào)的主要頻率范圍(Sorensen，Burrus，1993)，該方法適用于對(duì)低頻信號(hào)分析。

主動(dòng)源臺(tái)陣接收到的主要信號(hào)頻率為2～7 Hz，相對(duì)于100 Hz的采樣率，MCZT不僅能計(jì)算出了被FFT丟失的譜，而且譜峰的高度也更加精確。另外由于MCZT只計(jì)算了有用信息所在的頻段，因此還具有抑制寬帶噪聲的作用。

3 FICP算法

FICP最基本的思想是對(duì)2個(gè)信號(hào)的互譜補(bǔ)零，從而提高相關(guān)函數(shù)的分辨率。從采樣定理可知，采樣信號(hào)的頻譜相當(dāng)于原始信號(hào)頻譜的周期重復(fù)，且重復(fù)間隔等于采樣頻率。離散傅里葉變換(Discrete Fourier Transfrom，簡(jiǎn)稱DFT)計(jì)算出的信號(hào)頻譜，其右半部分是左半部分的共扼對(duì)稱，總長(zhǎng)度對(duì)應(yīng)采樣率數(shù)值。提高采樣率并增加序列點(diǎn)數(shù)(tp不增加)不會(huì)改變頻譜的形態(tài)，只會(huì)拉開(kāi)左右兩段譜之間的間隔。因此，如果把計(jì)算出的頻譜中間加零，拉開(kāi)左右距離，逆變換計(jì)算出的時(shí)域波形的采樣率就會(huì)提高(侯朝煥等，1990)。

對(duì)于2個(gè)信號(hào)x1(n)，x2(n)(n=0，1，…，N-1)，加窗函數(shù)W(n)(n=0，1，…，N-1)后，用MCZT分別計(jì)算其細(xì)化頻譜，根據(jù)相關(guān)定理計(jì)算得到2個(gè)信號(hào)的互譜：

R(k)=X1(k)·X2(k) (k=0，1，…，N-1)

(5a)

式中，R(k)只是譜的前一部分，通過(guò)譜的共軛對(duì)稱性，對(duì)R1(x)補(bǔ)零擴(kuò)展為N1點(diǎn)R(k)：

(5b)

通常信號(hào)的時(shí)延處在有限范圍內(nèi)，相關(guān)函數(shù)的主峰處在零值附近，因此只需要計(jì)算相關(guān)波形的左邊(后段)有限點(diǎn)和右邊(前段)有限點(diǎn)。實(shí)現(xiàn)這一算法的方法可以采用有限點(diǎn)時(shí)域(或頻域)序列計(jì)算在有限點(diǎn)頻域序列(或時(shí)域)的算法(Sorensen，Burrus，1993)該方法雖然考慮到了減少計(jì)算量的問(wèn)題，但是不能對(duì)變換結(jié)果進(jìn)行插值運(yùn)算。在計(jì)算時(shí)延時(shí)，只需要計(jì)算相關(guān)峰附近的點(diǎn)，為此可以對(duì)逆變換計(jì)算進(jìn)行簡(jiǎn)化。

按照式(5b)將互譜拉開(kāi)，設(shè)定互譜的長(zhǎng)度為N2， 對(duì)于n=0，1，…，N-1的時(shí)域范圍，相關(guān)函數(shù)前段為：

(6)

將式(6)兩項(xiàng)分別計(jì)算可以發(fā)現(xiàn)，這兩項(xiàng)分別為IMCZT計(jì)算過(guò)程：

(7)

式中：

(8)

B(n)=IMCZT[R3(k)]

(9)

式中：

(10)

同理，對(duì)于n=N2-N，N2-N+1，…，N2-1的時(shí)域范圍， 取n=N2-N+n1， 對(duì)應(yīng)的n1=0，1，…，N-1相關(guān)函數(shù)為：

(11)

式中：

(n=N2-N+n1；n1=0，1，…，N-1；k=0，1； …，N-1)

(12)

式中：

(13)

F(n)=G(n)·IMCZT[M(k)]

(14)

(15)

(16)

由此可組合出相關(guān)波的峰為：

智慧交通公共數(shù)據(jù)與服務(wù)支撐平臺(tái)，違章違停數(shù)據(jù)已有30多萬(wàn)條，公交車位置數(shù)據(jù)有3206萬(wàn)條，虛擬卡口過(guò)車數(shù)據(jù)約有1.6億條，路邊停車次數(shù)有50多萬(wàn)次。隨著各項(xiàng)數(shù)據(jù)的積累，未來(lái)的吉首將運(yùn)用大數(shù)據(jù)分析平臺(tái)提升吉首市的現(xiàn)代化水平，建立健全的大數(shù)據(jù)輔助科學(xué)決策和社會(huì)治理機(jī)制，推進(jìn)政府管理和社會(huì)治理模式的創(chuàng)新，實(shí)現(xiàn)政府決策科學(xué)化、社會(huì)治理精準(zhǔn)化和公共服務(wù)高效化。

(17)

通過(guò)設(shè)置N2，可以使相關(guān)函數(shù)的峰值分辨率提高N2/N1倍。

一方面，提高互譜的頻譜分辨率可以使相關(guān)波形更光滑，相關(guān)峰頂更圓，且頻譜的細(xì)化不受取樣長(zhǎng)度的影響；另一方面，提高插值的倍數(shù)，可以使相關(guān)峰頂分辨率提高?？梢?jiàn)采用FICP算法可以提高時(shí)延估計(jì)的精度，特別是對(duì)于低采樣的信號(hào)，其優(yōu)點(diǎn)更明顯。

4 計(jì)算分析

4.1 MCZT譜細(xì)化的計(jì)算

為考察MCZT和FFT在較短取樣長(zhǎng)度下譜的細(xì)化表達(dá)能力，將賓川氣槍源的氣泡振蕩頻率在2～7 Hz疊加合成一個(gè)仿真信號(hào)，由相對(duì)振幅分別為0.2，0.4，0.6，0.3，頻率分別為2.2，3.1，4.4，5.4 Hz，采樣率為100 Hz的4個(gè)正弦函數(shù)疊加而成，時(shí)長(zhǎng)為4 s(圖1)。

圖1 仿真信號(hào)Fig.1 The simulated signals

取1 s長(zhǎng)度(圖1中綠色部分)和2 s長(zhǎng)度分別做FFT和MCZT，得到的振幅譜如圖2所示。從圖2可以看出，在100 Hz采樣率，1 s取樣長(zhǎng)度下，雖然有一定的誤差，但MCZT能分辨出4個(gè)頻率成份，F(xiàn)FT只能分辨出2個(gè)頻率成份；在2 s取樣長(zhǎng)度，MCZT和FFT都能分辨出4個(gè)頻率成份，但FFT的頻率分辨率明顯比MCZT的低，誤差比MCZT大。

圖2 仿真信號(hào)不同長(zhǎng)度信號(hào)的MCZT細(xì)化振幅譜和FFT振幅譜Fig.2 MCZT amplitude spectrum and FFT amplitude spectrum of simulated signals with different lengths

4.2 基于MCZT和FICP的時(shí)延估計(jì)

為了評(píng)價(jià)基于MCZT和FICP的時(shí)延估計(jì)方法，先用模擬信號(hào)進(jìn)行仿真實(shí)驗(yàn)，仿真信號(hào)由距氣槍源5.5 km的烏稍臺(tái)(臺(tái)站代碼53272)的5 000槍記錄的格林函數(shù)疊加而成參考格林函數(shù)(圖3a)，走時(shí)在1～2 s之間的震相為P波震相，3～4 s之間的震相為S波震相。

根據(jù)伸縮法原理利用樣條函數(shù)進(jìn)行插值，產(chǎn)生1 000個(gè)具有時(shí)延在±0.5個(gè)采樣間隔之間3個(gè)周期正弦變化時(shí)延(圖3b)的氣槍模擬記錄(圖3c)，所有仿真記錄都將最大振幅歸一化為1。

圖3 參考格林函數(shù)(a)、時(shí)延變化(b)和具有時(shí)延變化的仿真信號(hào)(c)Fig.3 Reference Green function(a)，time delay change(b)，and simulated signals with time delay variation(c)

圖4 1 000道仿真信號(hào)(a)及其在191～192采樣點(diǎn)間的放大圖(b)Fig.4 Simulated signals at track 1 000(a)and it’s amplification between sampling sites of 191～192(b)

從圖4可以看出，最大變化量為5‰的速度變化，這在波形變化上是很細(xì)微的，不足一個(gè)采樣間隔，仿真信號(hào)僅有振幅上的微小差異，延時(shí)計(jì)算需要精細(xì)的譜估計(jì)。

Srec=max[xcorr(Dr，Dt)]/norm(Dt-Dr)

(18)

式中：xcorr(Dr，Dt)表示互相關(guān)；norm(Dt-Dr)表示殘差的絕對(duì)值。波速變化的恢復(fù)度值越大，波速變化的恢復(fù)效果越好。

對(duì)模擬信號(hào)添加高斯白噪聲，信噪比為-20～40 dB，得到不同信噪比下的仿真信號(hào)，當(dāng)信噪比≤5 dB時(shí)，使用曲波變換去噪。計(jì)算波速變化率時(shí)，參考格林函數(shù)為所有當(dāng)前格林函數(shù)的線性疊加，分別對(duì)每一條當(dāng)前格林函數(shù)進(jìn)行計(jì)算，格林函數(shù)震相采用P波震相，信號(hào)長(zhǎng)度約0.8 s共80個(gè)采樣點(diǎn)，得到波速變化率(圖5)。

圖6為信噪比10 dB時(shí)3種方法計(jì)算得到的波速變化率，圖7為計(jì)算與理論波速變化率的相關(guān)系數(shù)、殘差和波速變化恢復(fù)精度。從這2個(gè)圖可知，伸縮法(Stretch)、余弦插值法(CosInter)和本文方法(MCZT)都能較好地恢復(fù)波速變化率的形態(tài)，但伸縮法與理論波速變化率(Theory)殘差大，余弦插值法居中，余弦插值法和本文方法無(wú)論在形態(tài)或者量值上都較好地恢復(fù)了波速的變化率，但本文方法要優(yōu)于余弦插值法。

圖5 信噪比10 dB時(shí)的第1和84道格林函數(shù)Fig.5 Green functions at track 1 and 84 when the signal to noise ratio is 10db

數(shù)值模擬表明，本文方法相比余弦插值法較好地克服了信號(hào)窗較短譜估計(jì)精度低的問(wèn)題，因此提高了互相關(guān)估計(jì)走時(shí)差的精度。

圖6 信噪比10 dB時(shí)不同方法計(jì)算得到的波速變化率與理論波速變化率Fig.6 The wave velocity change rate of different methods and the theoretical wave velocity change rate when signal to noise ratio is 10db

圖7 在不同信噪比下測(cè)得的相關(guān)系數(shù)(a)、殘差(b)和波速變化恢復(fù)精度(c)Fig.7 The correlation coefficient(a)，residual(b) and recovery accuracy(c) of evaluation parameters measured under different signal to noise ratio of three methods

4.3 恢復(fù)度與窗長(zhǎng)的關(guān)系比較

分別計(jì)算伸縮法、余弦插值法和本文方法在同一信噪比下波速變化的恢復(fù)度與時(shí)窗長(zhǎng)度的關(guān)系，如圖8所示。

在無(wú)噪聲情況下(圖8a)，0.8～1.5 s對(duì)應(yīng)的窗長(zhǎng)范圍內(nèi)，波速變化率恢復(fù)度的自然對(duì)數(shù)為：伸縮法最低，本文方法最高，余弦插值法略低于本文方法。本文方法在窗長(zhǎng)大于1 s后，恢復(fù)度基本維持水平狀態(tài)，余弦插值法和伸縮法隨窗長(zhǎng)的增加，恢復(fù)度會(huì)呈緩慢增長(zhǎng)。綜合對(duì)比可知，本文方法可以在較短的窗長(zhǎng)下獲得更高恢復(fù)度及更高精度的波速變化率。從30 dB噪聲情況下計(jì)算的恢復(fù)度與窗長(zhǎng)關(guān)系(圖8b)來(lái)看，同樣是伸縮法最低，本文方法和余弦插值法較高，且本文方法優(yōu)于余弦插值法。在有噪聲情況下，波速的恢復(fù)程度與無(wú)噪聲情況具有相似的特征，隨著窗長(zhǎng)增加，恢復(fù)程度都緩慢升高。在計(jì)算窗長(zhǎng)范圍時(shí)，本文方法始終優(yōu)于余弦插值法、明顯高于伸縮法，說(shuō)明了本文方法在短時(shí)窗互相關(guān)時(shí)延估計(jì)中的優(yōu)勢(shì)。

圖8 無(wú)噪聲(a)和信噪比為30 db(b)情況下恢復(fù)度的自然對(duì)數(shù)與窗長(zhǎng)的關(guān)系Fig.8 Recovery degree and window length of no noise(a)and the signal to noise ratio is 30 db(b)

4.4 實(shí)際氣槍記錄計(jì)算

選取源檢距50 km的53284臺(tái)，共8 181條接收記錄。因源檢距較大，信號(hào)信噪比低，經(jīng)相位加權(quán)疊加和曲波變換去噪后，得到8 181條去噪記錄，分別用3種方法對(duì)11.0～11.30 s時(shí)窗信號(hào)進(jìn)行波速變化率計(jì)算，結(jié)果如圖9所示。從圖中看到，在信噪比對(duì)信噪較高的記錄，3種方法計(jì)算結(jié)果變化形態(tài)基本一致，伸縮法得到的波速變化率范圍大于余弦插值法，本文方法得到的波速變化率范圍較小。伸縮法和余弦插值法對(duì)信噪比過(guò)低的記錄所得到計(jì)算結(jié)果誤差較大，明顯超出合理范圍，而本文方法得到的波速變化率相對(duì)連續(xù)、跳躍不大、抗噪性更好，所計(jì)算的波速變化更穩(wěn)定。

圖9 53284臺(tái)波速變化率Fig.9 Wave velocity change rate of the station 53286

對(duì)得到的去噪信號(hào)剔除信噪比較低的信號(hào)后，得到6 816條記錄，最后得到信噪比較高、一致性較好的6 816條去噪記錄。對(duì)時(shí)窗范圍為10.0～11.3 s 信號(hào)，用3種方法分別計(jì)算其波速變化，以三者的結(jié)果平均值作為參考(標(biāo)準(zhǔn))波速變化，然后以10.60 s為開(kāi)始、步長(zhǎng)為0.02 s，與數(shù)值模擬方法一樣，分別計(jì)算不同窗長(zhǎng)與恢復(fù)度的關(guān)系(圖10)。計(jì)算結(jié)果表明：隨窗長(zhǎng)的增大，3種方法得到的波速變化的恢復(fù)度均升高，本文方法略優(yōu)于余弦插值法，伸縮法最低。

圖10 53284臺(tái)波速恢復(fù)度與計(jì)算信號(hào)窗長(zhǎng)的關(guān)系Fig.10 The relationship between wave velocity recovery and the calculated signal window length at the station 53284

5 結(jié)論

本文根據(jù)地下介質(zhì)波速變化測(cè)量精度的需求和賓川氣槍主動(dòng)源數(shù)據(jù)特點(diǎn)，介紹了一種能夠在短時(shí)窗互相關(guān)計(jì)算中獲得高精度時(shí)延估計(jì)的方法，通過(guò)實(shí)驗(yàn)對(duì)比分析，得到以下結(jié)論：

(1)基于MCZT和FICP的時(shí)延估計(jì)方法可以得到較高精度的時(shí)延估計(jì)值，通過(guò)模擬實(shí)驗(yàn)表明，該方法較余弦插值法有一定的提升，明顯優(yōu)于伸縮法。

(2)本文方法在短時(shí)窗互相關(guān)時(shí)差計(jì)算中優(yōu)勢(shì)較大，且在走時(shí)差計(jì)算過(guò)程中，通過(guò)使用MCZT，只計(jì)算有用信號(hào)所在頻段，實(shí)現(xiàn)對(duì)噪聲干擾的抑制，提高了地震信號(hào)時(shí)延估計(jì)的精度。

(3)利用FICP相關(guān)峰精確插值方法，克服了傳統(tǒng)方法插值結(jié)果不能完全反映相關(guān)峰形狀、插值計(jì)算而帶來(lái)信號(hào)誤差、增加運(yùn)算量等問(wèn)題，提高了相關(guān)峰頂?shù)姆直媛剩貏e是對(duì)低采樣信號(hào)優(yōu)勢(shì)明顯。

(4)實(shí)際數(shù)據(jù)計(jì)算表明，本文方法略優(yōu)于余弦插值方法，也具有一定的抗噪能力。