利用背景噪聲研究地震波速度變化及其在長白山火山監(jiān)測中的應(yīng)用

劉國明

(長白山天池火山監(jiān)測站，安圖 133613)

0 引言

火山監(jiān)測主要依靠火山地震活動探測、地表形變及流體地球化學(xué)等方法。應(yīng)用這些方法，可以探測到巖漿升壓和運移過程中的一些信息。地表形變學(xué)方法的局限性在于它們都是基于對地表應(yīng)變及變形的解釋，這就使得它們對深部物質(zhì)的變化不靈敏。通過對火山地震的精確定位，可以掌握火山地震的時空分布特征，從而使我們了解到火山區(qū)地下介質(zhì)屬性的變化信息。然而，在沒有發(fā)生巖漿活動的時間段，即在沒有火山地震發(fā)生的情況下，使用火山地震學(xué)的方法就很難探測到火山區(qū)地下介質(zhì)屬性的變化信息。

近年來，利用地震背景噪聲恢復(fù)格林函數(shù)的方法進行高分辨率的面波成像技術(shù)被廣泛應(yīng)用于空間結(jié)構(gòu)探測。而在一些地球物理研究領(lǐng)域中，探測地震波速度的連續(xù)瞬時變化比探測地下波速結(jié)構(gòu)更為重要，例如，研究與強震相關(guān)的波速變化(Brenguier et al.，2008a;Chen et al.，2010;Cheng et al.，2010;Froment et al.，2013)。對于火山監(jiān)測而言，即使在地下波速結(jié)構(gòu)未知的情況下，掌握地下介質(zhì)的相對波速變化也是非常重要的。Sens-Sch?enfelder等(2006)提出了應(yīng)用地震背景噪聲互相關(guān)技術(shù)反演火山區(qū)地下介質(zhì)瞬時波速變化的方法，并探測到了印度尼西亞默拉皮火山地震波相對速度的季節(jié)性變化;Brenguier等(2008b)報道了留尼汪島火山噴發(fā)前＜0.1% 的相對地震波速變化;Mordret等(2010)報道了新西蘭佩魯阿胡火山噴發(fā)2天前0.8%的相對地震波速度下降。應(yīng)用地震背景噪聲研究地震波速度相對變化的方法已經(jīng)被廣泛應(yīng)用于活動火山監(jiān)測研究中，但其完整的數(shù)據(jù)處理過程卻未見詳細報道。

長白山火山位于中國與朝鮮邊境，作為中國最具潛在噴發(fā)危險的活火山，近年來很多學(xué)者對這座火山進行了研究。吳建平等(2005)和劉國明等(2006，2011)研究了長白山火山的地震活動性，并指出自1999年以來長白山火山活動經(jīng)歷了3個明顯的階段;上官志冠(1997)和上官志冠等(2006)研究了長白山火山區(qū)的He同位素比值變化并指出長白山火山正處于巖漿活躍狀態(tài);李克等(2009)研究了長白山火山區(qū)流動GPS和流動水準監(jiān)測資料，發(fā)現(xiàn)在2002—2005年期間，長白山火山區(qū)垂直形變和水平形變都存在明顯的異常。

目前長白山火山監(jiān)測的主要方法有流動GPS、流動水準監(jiān)測、定點應(yīng)變和定點傾斜連續(xù)觀測、流體地球化學(xué)定期觀測和地震觀測等。流動GPS和流動水準測量能精確地反映火山錐體的水平和垂直變化信息，但是每年只進行一期觀測，時間分辨率太低，無法監(jiān)控到巖漿的動態(tài)變化過程。定點應(yīng)變和定點傾斜測量時間分辨率很高，但是目前的定點形變嚴重受當?shù)乩子甑挠绊懀瑪?shù)據(jù)的連續(xù)性非常不好。而且所有這些形變觀測都是基于對地表變形的解釋，無法很好地反映地下巖漿的微觀動態(tài)過程。流體地球化學(xué)方法可以反映深部的巖漿活動信息，但是目前采樣不連續(xù)，而且分析精度也不高，資料的利用價值不大。通過對火山地震的精確定位，可以了解巖漿活動的過程，但是用這種方法也無法掌握非震時段的巖漿活動情況。

雖然很多學(xué)者已經(jīng)從不同角度應(yīng)用不同的方法對長白山火山的活動狀態(tài)進行了深入的研究，但長白山火山巖漿活動在時間尺度上的動態(tài)過程至今仍不被人們了解。我們通過對地震背景噪聲進行互相關(guān)計算從而恢復(fù)格林函數(shù)的方法，獲得了長白山火山區(qū)地下介質(zhì)相對地震波速度變化過程，為研究長白山火山巖漿動態(tài)活動過程提供了一個新的方法。

1 單臺數(shù)據(jù)預(yù)處理

單臺數(shù)據(jù)預(yù)處理包括對每一個臺站數(shù)據(jù)的預(yù)處理過程，其目的是去除諸如地震、畸變等可能對正常的背景噪聲產(chǎn)生影響的信號，以突出背景噪聲信號。通常，先把原始地震數(shù)據(jù)轉(zhuǎn)換成以天為單位的SAC格式的波形數(shù)據(jù)，然后再進行數(shù)據(jù)處理。單臺數(shù)據(jù)預(yù)處理的操作流程:去除儀器響應(yīng)，去除均值，去除趨勢變化和帶通濾波，時間域歸一化和譜白化處理。其中一些數(shù)據(jù)處理過程，比如時間域歸一化和譜白化，對原始波形進行了非線性操作，因而這些數(shù)據(jù)處理過程的順序通常是不能更改的。

2 互相關(guān)計算和數(shù)據(jù)疊加

互相關(guān)是通過識別時間延遲來判斷2個波形相似程度的一種計算，所謂時間延遲是指一個信號相對另一個信號移動的尺度，通過這樣的移動使之與另一個信號達到最大程度的相似，這個移動尺度就叫時間延遲。2個信號a和b的互相關(guān)函數(shù)Ca，b是時間延遲τ的函數(shù)，通常被定義為

這里的積分長度是整個記錄的長度，u是時間信號函數(shù)的振幅。

互相關(guān)的作用在于它突出了地震波的走時。一個穿行于2個臺站間的波場會在每個臺站的地震記錄中產(chǎn)生相似的信號，只是在記錄到達的時間上會有偏移。因而臺站對記錄的互相關(guān)函數(shù)在某個時間延遲下會到達其峰值，這個峰值與臺站對間波場的地震波走時一致。

通常，在時間域里對2個臺站每天的數(shù)據(jù)進行互相關(guān)計算，然后對這些互相關(guān)函數(shù)進行疊加(stack)操作，以產(chǎn)生更長的時間序列，其目的是獲得更穩(wěn)定的互相關(guān)函數(shù)?；ハ嚓P(guān)函數(shù)是包含正半部分和負半部分的雙邊的時間函數(shù)，因此出現(xiàn)在互相關(guān)函數(shù)中的所謂面波的經(jīng)驗格林函數(shù)(GFs)也包含了穿行于臺站對2個方向間的射線路徑信息。

為了提高信噪比(SNR)，采用適當?shù)寞B加長度是有必要的。圖1顯示了CBS－YNB臺站對在不同疊加長度下獲得的面波波形。圖1a是經(jīng)過0.01～0.2Hz帶通濾波后獲得的，圖1b是經(jīng)過0.02～0.1Hz帶通濾波后獲得的。通過比較，圖1b中的面波波形更清晰，而圖1a中的面波中有少量的短周期信號干擾。圖1也表明選擇適當?shù)臑V波頻帶對獲取更穩(wěn)定的面波信號是非常重要的。

圖1 互相關(guān)函數(shù)中面波與疊加時間長度的關(guān)系Fig.1 Relationship between surface wave and stacking time.

3 測定相對地震波速度變化

在對所有的臺站對數(shù)據(jù)進行了互相關(guān)計算和疊加計算后，就可以測定相對地震波速度變化了。臺站對每天的互相關(guān)函數(shù)可以視為這2個臺站間的面波格林函數(shù)。為了提高計算的信噪比，通常把經(jīng)過連續(xù)的若干天疊加的互相關(guān)函數(shù)當作當天的經(jīng)驗格林函數(shù)。而把經(jīng)過更長時間(通常為1a或2a)疊加的互相關(guān)函數(shù)當作參考格林函數(shù)。

假設(shè)臺站對間地殼介質(zhì)的相對地震波速度變化d v/v在空間上是均勻變化的，那么經(jīng)驗格林函數(shù)和參考格林函數(shù)間的相對走時偏移dτ/τ是不依賴于流逝時間τ的，并且d v/v=－dτ/t(Brenguier et al.，2008b)。當測定某一個以τ為中心的小窗口的經(jīng)驗格林函數(shù)與參考格林函數(shù)之間的走時偏移dτ時，如果dτ相對于τ是線性變化的，那么就可以通過測定走時偏移dτ的斜率來估計經(jīng)驗格林函數(shù)和參考格林函數(shù)之間的相對走時變化dτ/τ，而其相反數(shù)就是地震波速度的相對變化d v/v。

在具體的計算過程中，使用“滑動窗互相關(guān)法(MWCCA)”得到d v/v(Snieder et al.，2002;Pandolfi et al.，2006)。如圖2所示，假設(shè)有2條波形曲線，fcur代表經(jīng)驗格林函數(shù)，fref代表參考格林函數(shù)，這2個波形的相關(guān)系數(shù)可以代表它們的相似程度。在要計算的波形中取一系列的小窗口，計算該窗口內(nèi)的走時擾動(Δτi)，根據(jù)尾波干涉理論，2個波形互相關(guān)函數(shù)達到最大值的位置所對應(yīng)的走時偏移即為走時擾動，相關(guān)系數(shù)R由式(1)求得:

式(1)中，ti為窗口中心點位置，窗口長度為2tw，ts為與最大相關(guān)系數(shù)對應(yīng)的走時偏移值，即走時擾動Δτi。在一系列小窗口內(nèi)重復(fù)上述計算，求得各個小窗口中心點位置的走時擾動Δτi，然后通過線性擬合得到整個經(jīng)驗格林函數(shù)相對參考格林函數(shù)的走時偏移Δτ/τ，其測量誤差為

圖2 相對走時偏移測量示意圖Fig.2 Schematic diagram showing relative shift time measurement.

4 質(zhì)量控制

評價互相關(guān)計算結(jié)果質(zhì)量好壞的一個主要指標是信噪比。通常信噪比被定義成某信號窗口內(nèi)最大振幅與噪聲窗口內(nèi)噪聲振幅均方根之比。質(zhì)量控制部分的目的是通過數(shù)據(jù)處理提高格林函數(shù)的信噪比，以得到穩(wěn)定可信的測量結(jié)果。

4.1 經(jīng)驗格林函數(shù)數(shù)據(jù)疊加天數(shù)的選取

在進行經(jīng)驗格林函數(shù)計算時，選取較長的數(shù)據(jù)天數(shù)可以獲得更穩(wěn)定的計算結(jié)果，提高格林函數(shù)的信噪比。但是，如果數(shù)據(jù)疊加天數(shù)過長，也會平滑掉可能有用的短期變化信息。因此選取適當?shù)臄?shù)據(jù)疊加天數(shù)，對獲得有意義的觀測結(jié)果是非常重要的。

圖3a和b顯示了不同疊加時間計算得到的信噪比和相關(guān)系數(shù)。隨著疊加時間的增長，測量結(jié)果的信噪比和相關(guān)系數(shù)得到了明顯的提高。圖3c顯示了在不同疊加時間情況下獲得的相對地震波速度變化值。11d數(shù)據(jù)疊加得到的測量結(jié)果雜亂無章，31d疊加和61d疊加獲得的觀測結(jié)果比較穩(wěn)定，且二者差別不大，所以在最后的數(shù)據(jù)處理過程中我們采用31d數(shù)據(jù)疊加作為經(jīng)驗格林函數(shù)。

從圖3c中還可以看到一個有趣的現(xiàn)象，在日本大地震前61d疊加得到的測量結(jié)果中有個明顯的震前上升過程，地震后又轉(zhuǎn)為急速下降。而在31d疊加的測量結(jié)果中，顯示的是震前下降，震后上升，這主要是因為61d數(shù)據(jù)疊加顯示的是長趨勢變化，相對而言，31d數(shù)據(jù)疊加顯示的是短期變化。2008年汶川地震時，波速變化也出現(xiàn)了類似變化?？紤]到汶川地震和日本地震距離長白山火山區(qū)都非常遠，還不能把這種變化認定為異常現(xiàn)象。

4.2 鐘差問題

一些研究在進行走時擾動計算前還進行了所謂的鐘差校正(Stehly et al.，2007;劉志坤等，2010)，該鐘差表現(xiàn)為每天的經(jīng)驗格林函數(shù)與參考格林函數(shù)之間的整體走時偏移，即圖2a中的擬合直線在通過走時0點時與0點值有較大偏差，而計算結(jié)果表明，本研究所用臺站鐘差較小(＜0.2s)，對波速測量結(jié)果基本沒有影響，因此本研究未進行鐘差校正處理。

5 長白山火山相對地震波速度變化

5.1 數(shù)據(jù)準備

本研究使用的數(shù)據(jù)是長白山火山區(qū)CBS、YNB、DHT和FST等4個寬頻帶地震臺記錄的2007—2012年的連續(xù)數(shù)字地震觀測資料，其中CBS－YNB臺站對計算了2001—2012年的地震波速度。對所有的臺站對進行了互相關(guān)計算，最終獲得了6個臺站對的相對地震波速度變化結(jié)果。圖4顯示了長白山火山地理位置和本次研究所用到的地震臺站分布情況。

5.2 數(shù)據(jù)處理過程

文中第2至第5部分討論了應(yīng)用背景噪聲研究相對地震波速度變化的方法和基本流程。本節(jié)將總結(jié)數(shù)據(jù)處理流程并重點介紹本次研究中使用的具體處理方法和數(shù)據(jù)處理流程。基本的數(shù)據(jù)處理步驟如下:

(1)用RDSEED軟件將原始數(shù)據(jù)(SEED格式)轉(zhuǎn)換成SAC格式。

(2)將SAC數(shù)據(jù)文件合并成以天為單位的數(shù)據(jù)文件。

圖3 疊加時間長度對計算結(jié)果的影響Fig.3 Influence of calculating results from length of stacking time.

圖4 長白山火山地理位置和火山監(jiān)測地震臺網(wǎng)Fig.4 The geographical position of Changbaishan volcano and its seismic network.

(3)單臺數(shù)據(jù)預(yù)處理及互相關(guān)計算。選定臺站的數(shù)據(jù)以10為采樣率進行重采樣，以減少后續(xù)數(shù)據(jù)處理的計算量。去除儀器響應(yīng)，去除均值和趨勢變化，在0.02～0.1Hz頻帶范圍內(nèi)進行帶通濾波，歸一化處理，最后進行譜白化處理。完成單臺數(shù)據(jù)預(yù)處理后，對臺站對的數(shù)據(jù)進行互相關(guān)計算。

(4)對所有已經(jīng)獲得的互相關(guān)數(shù)據(jù)進行信噪比計算?，F(xiàn)在計算信噪比的目的是保證互相關(guān)函數(shù)的質(zhì)量，所有信噪比低于某一數(shù)值的互相關(guān)函數(shù)將被刪除，具體取舍哪些數(shù)據(jù)，由研究者根據(jù)數(shù)據(jù)情況自行掌握。

(5)對選定臺站對的互相關(guān)函數(shù)進行31d數(shù)據(jù)疊加，并產(chǎn)生經(jīng)驗格林函數(shù)。對地震活動相對平靜的2a時間段的數(shù)據(jù)進行疊加，產(chǎn)生參考格林函數(shù)。本研究中，對2009—2010年的數(shù)據(jù)進行疊加，作為臺站對的參考格林函數(shù)。

(6)應(yīng)用MWCCA法，計算臺站對之間地下介質(zhì)相對地震波速度變化，同時，利用式(3)獲得計算誤差。

(7)計算經(jīng)驗格林函數(shù)的信噪比。

(8)計算結(jié)果選取。計算誤差反映了相對地震波速度計算過程中的不確定性，根據(jù)最后計算結(jié)果的實際情況，將計算誤差＞0.4% 的結(jié)果舍去。

5.3 討論

圖5顯示了全部參與計算的臺站對2007—2012年地震波速度的計算結(jié)果。圖5d—f的計算精度明顯高于圖5a—c，而且圖5d—f的長趨勢變化曲線有很好的相關(guān)性，顯示出很好的季節(jié)性變化規(guī)律。圖5a—c中的相對地震波速度計算結(jié)果都使用了FST的數(shù)據(jù)，我們認為導(dǎo)致圖5a—c計算結(jié)果不好的可能原因是FST的原始數(shù)據(jù)不理想，因而在接下來的討論中，我們只討論圖 5d—f。

在研究時段內(nèi)，發(fā)生了3個可能會對研究區(qū)地震波速度變化產(chǎn)生影響的地震。第1個是2008年5月12日汶川MW7.9地震，第2個是2010年10月9日發(fā)生的位于火山口20km的ML3.7和平林場地震，第3個是2011年3月11日發(fā)生的MW9.0日本大地震。因第1和第3個地震距離研究區(qū)太遠(均＞1 300km)，第2個地震的震級太小，因而在圖5 d—f的表現(xiàn)不明顯。在2007—2012年時間段內(nèi)，長白山火山區(qū)并未發(fā)生足以影響到地震波速度變化的地震事件，因而圖5所獲得的波速變化信息實際上是長白山火山區(qū)地震波速度變化的背景信息。

圖5 2007—2012年相對地震波速度變化Fig.5 The relative seismic velocity variation from 2007 to 2012.

2001—2007年6月以前研究區(qū)內(nèi)僅CBS臺和YNB臺有連續(xù)的數(shù)字化地震觀測資料，其數(shù)據(jù)保存格式都是EVT格式。2007年6月以后臺站進行完“十五”改造，各臺站數(shù)據(jù)都改為由測震服務(wù)器存儲，導(dǎo)出的數(shù)據(jù)也變成了SEED格式的數(shù)據(jù)?！笆濉备脑烨昂驝BS臺和YNB臺的觀測設(shè)備并未發(fā)生改變，僅僅是數(shù)據(jù)存儲的格式發(fā)生了變化，并不影響前后數(shù)據(jù)的連續(xù)性。為了研究這段時間內(nèi)長白山火山區(qū)地震波速度的變化情況，單獨針對這些EVT格式的數(shù)據(jù)再開發(fā)了相關(guān)的代碼，對這些數(shù)據(jù)進行了處理，從而獲得了CBS－YNB臺站對2001—2012年相對地震波速度變化情況(圖6)。

從圖6可以看出，2001年6月至2007年3月地震波速度的變化幅度達到2%，而其他時間段地震波速度變化幅度最大為1.5%，與此同時，長白山火山區(qū)出現(xiàn)巖漿擾動跡象，記錄到了大量的火山地震，精密水準測量結(jié)果和流體地球化學(xué)的研究結(jié)論也都證實了在這個時間段內(nèi)長白山火山經(jīng)歷了一次巖漿活動過程(劉國明等，2011)。

圖6 CBS－YNB臺站對2001—2012年相對地震波速度變化Fig.6 The relative seismic velocity variation of the CBS－YNB station pair from 2001 to 2012.

一個地區(qū)的應(yīng)力或應(yīng)變的變化會引起介質(zhì)中孔隙的擴張或者壓縮，可能在介質(zhì)中產(chǎn)生新的孔隙，也可能會造成孔隙中流體的重新分布，從而導(dǎo)致該地區(qū)地震波速度發(fā)生變化。同震斷層位移會導(dǎo)致區(qū)域應(yīng)力場的變化，從而在應(yīng)力拉張區(qū)域產(chǎn)生介質(zhì)孔隙的擴張或者在介質(zhì)中產(chǎn)生新的孔隙，而在應(yīng)力壓縮區(qū)域會產(chǎn)生介質(zhì)孔隙的壓縮(Nishimura et al.，2000)。因此，在應(yīng)力拉張區(qū)域可能會產(chǎn)生地震波速度的下降，而在應(yīng)力壓縮區(qū)域可能會產(chǎn)生地震波速度的上升。

巖漿活動是一個復(fù)雜的動態(tài)過程。巖漿上升過程是一個增壓的過程，體現(xiàn)在地震波速度變化上應(yīng)該是一個上升的過程;巖漿下降過程是一個減壓的過程，體現(xiàn)在地震波速變化上應(yīng)該是一個相對下降的過程;同時，巖漿上升或下降的過程也會導(dǎo)致火山區(qū)局部應(yīng)力場的重新調(diào)整，都會造成地震波速發(fā)生變化。我們推斷，此期間地震波速度的異常大幅度變化，可能是巖漿活動造成火山區(qū)地下介質(zhì)屬性變化的一種外在表現(xiàn)，是巖漿運移動態(tài)過程的一種表現(xiàn)。目前在長白山火山區(qū)開展的精密水準測量和流體地球化學(xué)測量都是季節(jié)性的觀測，基本上每年進行1期測量，然后將歷年的觀測結(jié)果進行對比分析，這些測量結(jié)果從不同角度證實了這次巖漿擾動事件，但是這些觀測每年只進行一期，無法反映巖漿運移的動態(tài)過程。而我們此次工作獲得了以天為單位的地震波速度變化，反映了巖漿活動的動態(tài)過程，從而彌補了其他觀測方法的不足。

6 結(jié)論

應(yīng)用背景噪聲恢復(fù)地震波速度變化進行活動火山監(jiān)測是近年來剛興起的一個比較新的研究方法，其數(shù)據(jù)處理過程總體上分為4個步驟:1)單臺數(shù)據(jù)預(yù)處理;2)互相關(guān)計算和數(shù)據(jù)疊加計算;3)相對地震波速度計算;4)質(zhì)量控制，包括誤差計算和計算結(jié)果的選取等。

計算了長白山火山區(qū)2007—2012年6個臺站對的相對地震波速度變化。在0.02～0.1Hz頻帶范圍內(nèi)，通過背景噪聲互相關(guān)方法獲得了比較理想的面波信號。其中3個臺站對的長趨勢速度變化具有非常好的相關(guān)性，顯示了明顯的季節(jié)性變化規(guī)律，認為這是長白山火山區(qū)地震波速度變化的背景信息。

2001—2007年CBS－YNB臺站對相對地震波速度變化幅度遠高于正常背景值，與此同時，長白山火山發(fā)生了一次巖漿擾動事件，因而推斷此期間地震波速度的異常變化可能與巖漿活動有關(guān)，巖漿的升壓或降壓作用導(dǎo)致了大量火山地震的發(fā)生，造成了火山錐體的異常升降變化，導(dǎo)致了溫泉逸出氣體濃度發(fā)生異常變化，同時也使火山區(qū)地下介質(zhì)的物理性質(zhì)發(fā)生了微觀的變化，這可能就是2001—2007年CBS－YNB臺站對相對地震波速度變化幅度異常的根本原因。

利用背景噪聲互相關(guān)研究地震波速度變化是一個很有前景的研究領(lǐng)域，因為實現(xiàn)起來比較方便，在預(yù)測火山活動方面具有巨大的潛力，已經(jīng)在活動火山監(jiān)測方面獲得了很多成功的應(yīng)用。因此，利用背景噪聲恢復(fù)地震波相對波速變化的研究將會在中國活動火山巖漿演化過程研究和火山活動趨勢預(yù)測中獲得越來越廣泛的應(yīng)用。