利用祁連山主動源資料研究2016年門源6.4級地震前后波速變化

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

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

0 引言

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

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

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

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

1 資料選取

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

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

2 原理及方法

2.1 原理

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

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

圖2 所選臺站走時剖面

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

當波速為v的介質發生d v的到時變化時，其走時延遲d t與流逝時間t間有如下關系

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

2.2 數據處理方法

2.2.1 利用互相關技術獲取激發時刻

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

2.2.2 數據預處理

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

2.2.3 互相關時延檢測

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

3 地震前后走時變化測量結果

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

圖3 互相關時延檢測示意

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

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

圖4 震源區內所選臺站的走時變化及時間窗選取

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

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

4 結論與討論

4.1 誤差分析

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

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

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

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

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

4.2 結論與討論

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

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

祁連山主動源氣槍激發觀測數據的處理和應用研究才剛剛起步，本文得到的結果較為初步，但研究結果將有助于觀測數據結果的快速產出和地震相關波速變化信息的分析，以期盡快為地方防震減災提供服務。

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