鄒銳 郭曉 張元生 秦滿忠 顏文華

中國地震局蘭州地震研究所，蘭州市東崗西路450號 730000

0 引言

隨著觀測技術(shù)的發(fā)展和地震臺站的密集化，地震學家對天然地震波的傳播規(guī)律以及地球內(nèi)部結(jié)構(gòu)的研究也越來越多。地震通常被認為是由地下應(yīng)力的積累和釋放所造成的，因此，地下介質(zhì)應(yīng)力狀態(tài)及其變化是地震發(fā)生的關(guān)鍵因素（陳蒙，2015），而介質(zhì)應(yīng)力狀態(tài)的變化則會引起地震波波速的變化（Birch，1960、1961；Scholz，1968；Nur et al，1969）。迄今為止，地震波是所知的唯一能夠穿透地球內(nèi)部的振動（陳颙等，2005），因此，通過測量地震前后的波速變化，認識地下介質(zhì)應(yīng)力狀態(tài)的演化，對地震預(yù)測研究具有重要意義。

目前，國內(nèi)外有關(guān)學者分別利用重復(fù)地震、噪聲和人工震源等研究地下介質(zhì)的動態(tài)變化（Schaff et al，2004；Brenguier et al，2008；Siliver et al，2007；Niu et al，2008；陳颙等，2006、2007a；林建民等，2006；Wang et al，2008；王偉濤等，2009）。但重復(fù)地震發(fā)生的時間和地點不可控制，加之天然地震定位精度有限，因此，無法進行長期的動態(tài)監(jiān)測；而噪聲源能量較弱，需要長時間疊加才能獲得可靠的測量結(jié)果。上述因素都限制了利用天然源測量地下介質(zhì)時空變化的精度及分辨率（王寶善等，2011）。而人工震源具有震源位置已知、激發(fā)時間可控、觀測系統(tǒng)分布靈活方便、可以進行密集觀測等優(yōu)點，在區(qū)域尺度研究中有望彌補天然源在精度上的不足（楊微等，2013）。利用人工震源進行主動地震探測已成為高精度探測地下介質(zhì)波速變化的方法之一，而利用大容量氣槍震源作為主動震源對區(qū)域尺度進行地下介質(zhì)探測和監(jiān)測則更具有高度的可重復(fù)性、能量強、探測距離遠、綠色環(huán)保等優(yōu)勢。

通常，地下介質(zhì)的應(yīng)力變化所引起的地震波波速變化很小，故需要不斷提高地震波波速測量的精度。測量波速的幾種常用方法中，干涉法是測量波速精度最高的方法，而基于互相關(guān)時延檢測的干涉法已被廣泛用于地震波波速測量中，并取得一些顯著的進展。Niu等（2008）在美國San Andreas斷層上的SAFOD井內(nèi)采用主動源的方法，利用互相關(guān)延遲檢測技術(shù)，在井深1m處成功觀測到了由大氣壓和固體潮引起的波速變化，同時還觀測到在距井3km以外的1次M3.0地震震前，波速明顯地增加（8%）。Wang等（2008）在云南昆明進行的主動源探測中，利用基于互相關(guān)時延檢測技術(shù)的尾波干涉方法，也觀測到了介質(zhì)波速變化與大氣壓變化間有良好的相關(guān)性。

地震波走時與波速間密切相關(guān)，精確測量走時變化可獲取地震孕育過程中的有關(guān)信息。本文利用祁連山主動源2015年7月～2016年5月的氣槍激發(fā)數(shù)據(jù)，選取2016年青海門源6.4級地震震中附近以及距門源地震震中較遠的相關(guān)臺站接收到的氣槍激發(fā)信號，采用互相關(guān)時延檢測技術(shù)，對青海門源6.4級地震前后的走時變化進行了測量，并對地震前后波速的變化進行了初步分析。

1 資料選取

祁連山主動源激發(fā)系統(tǒng)自2015年7月運行以來，為盡可能減少噪聲對激發(fā)信號的干擾，選擇每周一晚至周二凌晨進行連續(xù)激發(fā)，每次激發(fā)間隔時間約為12min，每周可以獲得40～50次的激發(fā)數(shù)據(jù)。根據(jù)秦滿忠等（2017）的研究，祁連山主動源每周激發(fā)數(shù)據(jù)經(jīng)過疊加之后可被200km以外的臺站觀測到。其中，2015年9月30日～11月6日進行了24h的連續(xù)激發(fā)實驗。截至2016年7月1日，累計有效激發(fā)約6400次，數(shù)據(jù)基本連續(xù)，為后續(xù)的研究提供了可靠的數(shù)據(jù)保障。

2016年1月21日，青海門源發(fā)生了6.4級地震，此次地震震中恰好位于祁連山主動源監(jiān)測區(qū)域的邊緣，距祁連山主動源激發(fā)場地178km，最近的主動源臺站（ZDY38號臺）距震中25km。白超英（1999）研究認為，1個6級地震的震源區(qū)范圍一般為數(shù)十千米，為此選取震源區(qū)內(nèi)ZDY38、ZDY37、喜馬拉雅二期流動62430臺以及距地震震中較遠（距氣槍震源較近）的臺站ZDY28、ZDY30、ZDY31進行處理、對比分析（圖1、2）。

2 原理及方法

2.1 原理

移動窗互相關(guān)是互相關(guān)延時檢測中最常用的方法，其原理是對比2列高度相似的波形，由參考波形干涉另一列波形，給定1個時間窗，通過移動該時間窗求取互相關(guān)函數(shù)最大值，以獲得走時延遲進而求得波速變化。定義記錄點的初始波形為u（t），擾動后的波形為～u（t），選取時間窗的長度為T，擾動前后的時間窗T內(nèi)的互相關(guān)函數(shù)R可表示為

圖1 選取臺站、主動源位置及門源地震震源區(qū)位置

圖2 所選臺站走時剖面

其中，ts為時間變量；t為時間窗T的中心位置，當ts＝t、互相關(guān)函數(shù)取最大值時，則稱t為時間窗T內(nèi)的走時延遲。

當波速為v的介質(zhì)發(fā)生d v的到時變化時，其走時延遲d t與流逝時間t間有如下關(guān)系

假設(shè)介質(zhì)中速度變化均勻，如果介質(zhì)波速降低（或增大），則走時延遲將隨流逝時間而線性增加（或減?。⊿nieder，2006），線性擬合走時延遲與流逝時間之間所得直線斜率的負數(shù)即為介質(zhì)波速的相對變化。

2.2 數(shù)據(jù)處理方法

2.2.1 利用互相關(guān)技術(shù)獲取激發(fā)時刻

對2015年7月～2016年5月間的數(shù)據(jù)進行篩選，手動去除信噪比較差和記錄錯誤的個別波形。由于波形存在一些非介質(zhì)變化所引起的信號影響（如直流量、趨勢項等），需要同時對數(shù)據(jù)進行去均值、去趨勢等處理。然后，利用互相關(guān)技術(shù)獲取激發(fā)時刻，以祁連山主動源發(fā)射場地旁邊的ZDY22號臺（距氣槍激發(fā)源約80m）為參考臺，截取某個晚上參考臺接收到的激發(fā)信號為原始模板，分別對參考臺每天的數(shù)據(jù)進行互相關(guān)，得到激發(fā)時刻（圖3）。

2.2.2 數(shù)據(jù)預(yù)處理

以激發(fā)時刻為零時，對每個臺的數(shù)據(jù)截取激發(fā)時刻后200s的數(shù)據(jù)，采樣率為100Hz。為了消除觸發(fā)誤差、水位等因素對氣槍震源的影響，將每個臺截取的數(shù)據(jù)分別對參考臺接收到的激發(fā)信號進行反褶積。由于單次激發(fā)的氣槍信號會湮沒在噪聲中，因此，需要通過疊加來提高信噪比，以獲得有效的氣槍激發(fā)信號，故將所選臺站的所有數(shù)據(jù)進行疊加作為參考模板。對震源區(qū)內(nèi)臺站每2周的數(shù)據(jù)進行疊加，可得到時間間隔為14天的格林函數(shù)；對于距主動源發(fā)射場較近的臺站，則可得到時間間隔為7天的格林函數(shù)。在2015年10～11月不間斷激發(fā)實驗期間，數(shù)據(jù)疊加的時間間隔為3天。最后，帶通濾波至3～7Hz，以減少噪聲的影響。

2.2.3 互相關(guān)時延檢測

通過對每2周得到的格林函數(shù)與相應(yīng)臺站的參考模板進行移動窗互相關(guān)，在得到最大相關(guān)系數(shù)時可獲取走時延遲。原始數(shù)據(jù)的采樣率為100Hz，這意味著獲得的最小走時延遲為0.01s。但有時相關(guān)函數(shù)的最高峰值不一定在采樣點上，而是偏離采樣點一定的距離。為此，我們選擇余弦插值的方法來恢復(fù)相關(guān)函數(shù)的峰值位置，從而得到亞采樣點的時間精度（王偉濤等，2009）。圖3（d）、3（e）分別為插值前、后走時延遲圖。由圖3（d）、3（e）可見，使用余弦插值后可得到更高精度的走時延遲。根據(jù)式（2），對走時延遲求取平均值與方差，即得到相對平均波速的相對變化量與計算誤差。

3 地震前后走時變化測量結(jié)果

ZDY38臺站距門源地震震中25km，選取P波最大振幅33.06～33.94s、S波最大振幅56.46～57.60s、S波58.26～59.37s等3段時間窗，通過互相關(guān)時延檢測計算其走時變化（圖4）。由圖4可見，2015年8月初至9月底走時變化出現(xiàn)持續(xù)近2個月的低值異常，隨后又逐漸回升。而在2015年11月4日，走時變化又出現(xiàn)第2次下降，幅度最大達－18ms，低值異常持續(xù)至2016年1月21日門源6.4級地震發(fā)生，震后走時變化逐漸回升，至2016年3月25日，走時變化逐漸恢復(fù)正常。

圖3 互相關(guān)時延檢測示意

對于距門源地震震中41km的ZDY37臺站，計算其P波最大振幅26.02～26.95s、S波最大振幅47.10～48.63s、S波49.03～51.87s等3組時間窗內(nèi)的走時變化（圖4）。由圖4可見，ZDY37臺站與ZDY38臺站走時變化趨勢較為一致，2015年8月初至9月底出現(xiàn)第1次低值異常，之后恢復(fù)正常。2015年11月4日走時變化再次出現(xiàn)下降，最低值出現(xiàn)在2016年12月28日，走時變化幅度最大幅度為－16ms，隨后走時變化開始上升，在上升過程中發(fā)生了門源6.4級地震，至2016年3月14日恢復(fù)正常。

對于距門源地震震中53.8km的62430臺站，計算其P波最大振幅前后38.15～39.04s、S波最大振幅前后63.63～64.64s、S波65.77～67.15s等3組時間窗內(nèi)的走時變化（圖5）。由圖5可見，2015年7月28日走時變化出現(xiàn)第1次低值，并持續(xù)至9月底，2015年9月30日恢復(fù)正常。而在2015年11月4日，走時變化第2次出現(xiàn)低值，這與ZDY38、ZDY37兩個臺有所不同，62430臺的走時變化在門源地震震前一直持續(xù)低值，最大變化幅度達－13ms，門源地震發(fā)生后，走時變化于2015年1月25日出現(xiàn)最低值，隨后走時變化開始逐漸上升，直到2015年3月7日恢復(fù)正常。

圖4 震源區(qū)內(nèi)所選臺站的走時變化及時間窗選取

楊微等（2010）研究綿竹MS5.6地震前后波速的變化后認為，地震發(fā)生前，斷裂帶附近區(qū)域應(yīng)力逐漸積累，波速增加，而地震發(fā)生后應(yīng)力的釋放導(dǎo)致波速相對下降。2016年門源6.4級地震前，遠離震中的ZDY28、ZDY30、ZDY31等臺站走時變化趨勢較小，沒有出現(xiàn)明顯的低值異常過程（圖5）。地震前約6個月，震中附近的ZDY38、ZDY37、62430等臺站的相對走時出現(xiàn)下降變化（走時減少），至震前約3個月低值異?；謴?fù)，之后再次出現(xiàn)走時下降變化，這可能與孕震區(qū)內(nèi)地下介質(zhì)應(yīng)力持續(xù)積累而導(dǎo)致穿過該區(qū)域的地震波速度增加有關(guān)。

圖5 遠離震中所選臺站的走時變化及時間窗選取

4 結(jié)論與討論

4.1 誤差分析

走時延遲的誤差下限στ可由Cramer-Rao Lower Bound法則來計算（劉自鳳等，2015），即

其中，f0為信號主頻；T為時間窗長度；B為信號的頻寬比；ρ為波形的相關(guān)系數(shù)；SNR為信噪比。由式（3）可見，波形的相似度、信噪比和信號的頻寬與主頻之比等是影響走時延遲測量精度的主要因素（張金川等，2014）。

祁連山主動源的主頻基本為5Hz，信號的能量主要集中在4.5～6.5Hz，故選取f0≈5Hz，B≈0.4，在進行互相關(guān)時延檢測時，選取窗長T＝0.4s，選取ρ≈1，SNR≈15，則式（3）可簡化為

式（4）可得，走時延遲的理論誤差下限約為2.37×10－4s。而在進行互相關(guān)時求得走時延遲的均方差約為6.15×10－4s，可見實測誤差與理論誤差非常接近，但比ZDY38、ZDY37、62430等臺觀測到的走時變化值小1個數(shù)量級。

除計算誤差以外，測量的精度還受觀測系統(tǒng)和環(huán)境因素等的影響（楊微等，2010）。祁連山主動源臺站采用GPS連續(xù)授時方式，其授時精度為10－7μs，比計算得到的直達S波走時變化值的精度高3個數(shù)量級。固體潮和氣壓的影響是一個復(fù)雜的問題，數(shù)據(jù)的疊加對這些影響有一定的消弱作用，故疊加處理后的波形數(shù)據(jù)受固體潮和氣壓的影響可能不大。

4.2 結(jié)論與討論

本文利用祁連山主動源2015年7月～2016年5月的氣槍激發(fā)數(shù)據(jù)，選取2016年青海門源6.4級地震震中附近、距震中較遠的相關(guān)臺站接收到的氣槍激發(fā)信號，采用互相關(guān)時延檢測技術(shù)對青海門源6.4級地震前后的走時變化進行了測量，結(jié)果表明，遠離震中的ZDY28、ZDY30、ZDY31等臺站的走時變化趨勢較小，沒有出現(xiàn)明顯的低值異常過程；在地震前約6個月，震中附近3個臺站的相對走時出現(xiàn)下降變化（走時減少），至震前約3個月低值異?；謴?fù)，之后再次出現(xiàn)走時下降變化，門源6.4級地震發(fā)生于走時變化恢復(fù)過程中。S波走時變化最大下降幅度達－18ms，震后走時變化逐漸恢復(fù)正常，且3個臺站的變化趨勢較為一致。

地震發(fā)生前，斷裂帶附近區(qū)域應(yīng)力逐漸積累，波速增加，而地震發(fā)生后應(yīng)力釋放而導(dǎo)致波速相對下降（楊微等，2010）。走時縮短意味著速度增加，這可能與區(qū)域應(yīng)力積累間存在一定的關(guān)系。遠離震中的ZDY28、ZDY30、ZDY31等臺站走時變化一直較小，震中附近ZDY38、ZDY37、62430等3個臺站的震前走時變化持續(xù)低值，可能與2016年門源6.4級地震孕震過程中由于應(yīng)力的積累而導(dǎo)致的穿過該區(qū)域的地震波速度的增加有關(guān)。ZDY38臺站震前S波平均波速相對變化為0.38%，ZDY37號臺站為0.27%，62430臺站為0.15%，由此可見，距昌馬-俄博斷裂和門源6.4級地震震中最近的ZDY38臺站波速相對變化最明顯。

祁連山主動源氣槍激發(fā)觀測數(shù)據(jù)的處理和應(yīng)用研究才剛剛起步，本文得到的結(jié)果較為初步，但研究結(jié)果將有助于觀測數(shù)據(jù)結(jié)果的快速產(chǎn)出和地震相關(guān)波速變化信息的分析，以期盡快為地方防震減災(zāi)提供服務(wù)。

致謝：互相關(guān)時延計算程序由中國地震局地球物理研究所地震觀測與地球物理成像重點實驗室王寶善研究員提供，在此表示衷心感謝。