利用祁連山主動(dòng)源資料研究2016年門源6.4級(jí)地震前后波速變化

鄒銳 郭曉 張?jiān)?秦滿忠 顏文華

中國地震局蘭州地震研究所，蘭州市東崗西路450號(hào) 730000

0 引言

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

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

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

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

1 資料選取

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

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

2 原理及方法

2.1 原理

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

圖1 選取臺(tái)站、主動(dòng)源位置及門源地震震源區(qū)位置

圖2 所選臺(tái)站走時(shí)剖面

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

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

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

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

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

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

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

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

2.2.3 互相關(guān)時(shí)延檢測(cè)

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

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

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

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

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

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

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

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

圖5 遠(yuǎn)離震中所選臺(tái)站的走時(shí)變化及時(shí)間窗選取

4 結(jié)論與討論

4.1 誤差分析

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

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

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

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

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

4.2 結(jié)論與討論

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

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

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

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