楊潤海，譚俊卿，向 涯，姜金鐘，王 彬

(1.云南省地震局，云南 昆明 650224；2.云南大學，云南 昆明 650091；3.中國地震局地震研究所 中國地震局地震大地測量重點實驗室，湖北 武漢 430071)

0 引言

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

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

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

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

1 氣槍源地震波速干涉測量方法

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

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

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

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

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

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

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

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

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

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

2 MCZT譜細化

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

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

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

(1)

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

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

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

(2)

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

MCZT定義為：

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

(3)

定義MCZT的逆變換為IMCZT：

(4)

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

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

3 FICP算法

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

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

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

(5a)

式中，R(k)只是譜的前一部分，通過譜的共軛對稱性，對R1(x)補零擴展為N1點R(k)：

(5b)

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

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

(6)

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

(7)

式中：

(8)

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

(9)

式中：

(10)

同理，對于n=N2-N，N2-N+1，…，N2-1的時域范圍， 取n=N2-N+n1， 對應的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ù)與服務支撐平臺，違章違停數(shù)據(jù)已有30多萬條，公交車位置數(shù)據(jù)有3206萬條，虛擬卡口過車數(shù)據(jù)約有1.6億條，路邊停車次數(shù)有50多萬次。隨著各項數(shù)據(jù)的積累，未來的吉首將運用大數(shù)據(jù)分析平臺提升吉首市的現(xiàn)代化水平，建立健全的大數(shù)據(jù)輔助科學決策和社會治理機制，推進政府管理和社會治理模式的創(chuàng)新，實現(xiàn)政府決策科學化、社會治理精準化和公共服務高效化。

(17)

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

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

4 計算分析

4.1 MCZT譜細化的計算

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

圖1 仿真信號Fig.1 The simulated signals

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

圖2 仿真信號不同長度信號的MCZT細化振幅譜和FFT振幅譜Fig.2 MCZT amplitude spectrum and FFT amplitude spectrum of simulated signals with different lengths

4.2 基于MCZT和FICP的時延估計

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

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

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

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

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

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

(18)

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

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

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

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

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

圖6 信噪比10 dB時不同方法計算得到的波速變化率與理論波速變化率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 在不同信噪比下測得的相關(guān)系數(shù)(a)、殘差(b)和波速變化恢復精度(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 恢復度與窗長的關(guān)系比較

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

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

圖8 無噪聲(a)和信噪比為30 db(b)情況下恢復度的自然對數(shù)與窗長的關(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 實際氣槍記錄計算

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

圖9 53284臺波速變化率Fig.9 Wave velocity change rate of the station 53286

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

圖10 53284臺波速恢復度與計算信號窗長的關(guān)系Fig.10 The relationship between wave velocity recovery and the calculated signal window length at the station 53284

5 結(jié)論

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

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

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

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

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