MSNoise：利用地震背景噪聲監(jiān)測(cè)地震波速變化的Python程序包

Thomas LecocqCorentin CaudronFlorent Brenguier

0引言

最近，從成像到監(jiān)測(cè)，地震背景噪聲在不同領(lǐng)域都證明了它的效能。兩個(gè)傳感器之間的脈沖響應(yīng)[或格林函數(shù)(GF)]可以通過所記錄地震噪聲的互相關(guān)加以重構(gòu)(Campillo and Paul，2003)。從全球到區(qū)域范圍或局部地點(diǎn)的地球內(nèi)部成像來看，這個(gè)方法都給出了非常出色的結(jié)果。近來，這個(gè)方法又被拓展用于研究這些格林函數(shù)隨時(shí)間的變化。其走時(shí)延遲隨時(shí)間的變化可能源于介質(zhì)速度的改變，也可能源于一個(gè)/多個(gè)散射體或噪聲源的位置的顯著移動(dòng)。一些利用地震背景噪聲的研究揭示出，火山體內(nèi)部結(jié)構(gòu)的小擾動(dòng)可用地震波性質(zhì)的變化來探測(cè)(Sens-Sch?nfelder and Wegler，2006；Brenguier，Shapiro，etal，2008；Duputeletal，2009；Mordretetal，2010；Brenguieretal，2011；Anggonoetal，2012)。與使用主動(dòng)源或地震尾波方法不同的是：這個(gè)技術(shù)使用地震臺(tái)站24小時(shí)不間斷的記錄對(duì)介質(zhì)進(jìn)行連續(xù)采樣。這個(gè)方法已經(jīng)證明了它能夠發(fā)現(xiàn)在斷層帶(Wegler and Sens-Sch?nfelder，2007；Brenguier，Campilloetal，2008)和月球環(huán)境(Sens-Sch?nfelder and Wegler，2011)下物理隨時(shí)間變化的證據(jù)，或者能檢測(cè)出儀器問題(Stehlyetal，2007；Sens-Sch?nfelder，2008)。

在諸如地震分析代碼(Seismic Analysis Code，SAC)(Goldsteinetal，2003)或地震學(xué)計(jì)算程序(Computer Programs in Seismology，CPS)3.3版(Herrmann，2002)中，都有一些計(jì)算地震噪聲互相關(guān)的程序。據(jù)我們所知，至今還沒有發(fā)布一個(gè)集成的解決方案，能夠?qū)υ疾ㄐ翁幚淼玫阶邥r(shí)變化，包括自動(dòng)檢測(cè)歸檔數(shù)據(jù)中的變化，并預(yù)先設(shè)定處理計(jì)劃(每小時(shí)，每天)只進(jìn)行必要的計(jì)算，同時(shí)它可作為一個(gè)工具來研究歸檔的波形庫(kù)。自動(dòng)化檢測(cè)很重要，因?yàn)椋鐏碜詫?shí)時(shí)遙測(cè)臺(tái)站的數(shù)據(jù)流包含斷點(diǎn)，這些斷點(diǎn)隨后會(huì)被填補(bǔ)，并且數(shù)據(jù)流提供用于分析的重要數(shù)據(jù)。這種工具必須能夠處理任何常用的地震格式，并且盡可能快速和優(yōu)化。為了可插性和可擴(kuò)展，它應(yīng)該與數(shù)據(jù)檔案和具有高級(jí)輔助功能的數(shù)據(jù)庫(kù)進(jìn)行交互。最后，它必須產(chǎn)出可輸出的數(shù)據(jù)，可以是波形格式、表格式的文本文件，也可以是高質(zhì)量的圖形。這就是地震噪聲用于監(jiān)測(cè)(即MSNoise，Monitoring using Seismic Noise)程序的目的。

我們不愿意將MSNoise系統(tǒng)作為另一個(gè)黑箱提供給用戶，因此它作為一個(gè)集成的解決方案，其源代碼是公開的，有注釋和說明文檔。處理流程分成若干步驟，只要遵循某一模塊的輸入和輸出格式，用戶就可以很容易地用其他代碼代替原有模塊中的程序。

圖1　MSNoise工作流程示意圖。一次性的安裝部分在新任務(wù)定義步驟之前進(jìn)入到日常工作流程

在這篇文章中，我們從信息技術(shù)和科學(xué)算法的兩個(gè)角度給出了軟件的實(shí)現(xiàn)步驟，并配有圖形輸出的例子來進(jìn)行展示和描述。最后我們使用富爾奈斯火山(法國(guó)留尼旺島)的數(shù)據(jù)驗(yàn)證了MSNoise的可靠性；該火山在UnderVolc研究計(jì)劃(Brenguieretal，2012)期間分別于2010年10月和12月發(fā)生了兩次噴發(fā)，并在兩次噴發(fā)前觀測(cè)到了前兆性的地震波速變化。UnderVolc是一項(xiàng)隸屬法國(guó)國(guó)家科研署(ANR)的研究項(xiàng)目，名字的縮寫代表“理解火山過程(UNDERstanding VOLCanic Processes)：致力于噴發(fā)預(yù)測(cè)和風(fēng)險(xiǎn)降低的研究，應(yīng)用于法國(guó)留尼旺島的富爾奈斯火山地區(qū)(2009—2013)”。

MSNoise用Python語言編寫，是完全跨平臺(tái)的。由于是用于非盈利性的教育和研究，軟件本身是開源免費(fèi)的。MSNoise在歐盟公共許可證(EUPLv1.1)下被許可使用。未經(jīng)與原著者簽署協(xié)議，禁止與MSNoise相關(guān)的商業(yè)服務(wù)，諸如提供安裝、維護(hù)方面的咨詢或支持。MSNoise可在http：//www.msnoise.org/網(wǎng)頁(最后訪問時(shí)間2014年3月)上獲得。

1軟件功能

MSNoise的工作流程非常簡(jiǎn)單明了(圖1)，由一些簡(jiǎn)單的步驟組成，在配套的網(wǎng)站上有更多的描述。

1.1　安裝和配置

MSNoise需要Python和其他幾個(gè)程序包，以便在任何操作系統(tǒng)(OS)下運(yùn)行。波形歸檔文件必須是已知的或自定義的格式，例如已知的格式可以是SeisComP數(shù)據(jù)結(jié)構(gòu)(SDS)或者是統(tǒng)一數(shù)據(jù)緩沖區(qū)(BUD)結(jié)構(gòu)。MSNoise使用數(shù)據(jù)庫(kù)來儲(chǔ)存配置參數(shù)、波形元數(shù)據(jù)和任務(wù)。數(shù)據(jù)庫(kù)可以是MySQL或者sqlite，但由于性能的原因，我們推薦使用MySQL數(shù)據(jù)庫(kù)。安裝腳本隨程序包一起提供，用于初始化數(shù)據(jù)庫(kù)連接，存儲(chǔ)本地用戶名和密碼，使軟件配置器以及其他處理進(jìn)程能夠與數(shù)據(jù)庫(kù)交互。軟件配置器的圖形用戶接口[GUI，使用Enthought工具套件(Enthought，2008)構(gòu)建]顯示三個(gè)不同的面板：MSNoise的一般配置、臺(tái)站配置和濾波器配置。在線文檔中描述了所有的配置細(xì)節(jié)。

圖2　UnderVolc項(xiàng)目的數(shù)據(jù)可用性圖：上圖中灰色部分表示包含數(shù)據(jù)的日期，下圖匯總了每天的可用臺(tái)站數(shù)N。這幅圖是利用腳本plot_data_availability.py自動(dòng)繪制的

1.2　自動(dòng)任務(wù)定義

1.2.1數(shù)據(jù)發(fā)現(xiàn)

為了每晚自動(dòng)運(yùn)行前一天獲得的數(shù)據(jù)，MSNoise需要檢查數(shù)據(jù)檔案，查找新的或修改過的文件。這些文件可能是前一天獲得的，也可能是從以前的離線臺(tái)站獲得的包含有用信息的數(shù)據(jù)，比如一個(gè)月以前的數(shù)據(jù)。搜索時(shí)間在config(配置)中定義。MSNoise使用帶有-mtime參數(shù)的find(查找)命令(Windows系統(tǒng)下為gnufind)來定位新的或修改過的文件。一旦有文件被定位并讀取[使用obspy(Beyreutheretal，2010；Me-

giesetal，2011)]，這些文件將被插入(若是新的)或者更新(若是修改過的)到data availability(數(shù)據(jù)可用性)表中。之后這張表可被轉(zhuǎn)化為一張數(shù)據(jù)可用性圖(圖2)。腳本首次運(yùn)行調(diào)用時(shí)，必須帶著初始化參數(shù)init，目的是避免在find命令中使用mtime參數(shù)，這樣才能將所有被發(fā)現(xiàn)的文件都插入到data availability(數(shù)據(jù)可用性)表中。

1.2.2任務(wù)定義

MSNoise一旦識(shí)別出所有新的或修改過的文件，就會(huì)確定需要處理哪些天哪些臺(tái)站對(duì)。對(duì)于N個(gè)臺(tái)站，要計(jì)算M組臺(tái)站對(duì)(公式1a)。如果需要波形自相關(guān)，M則由公式(1b)來定義：

(1a)

(1b)

1.3　輸出、擴(kuò)展性和可插性

MSNoise最初的目的是用來提供隨時(shí)間變化的δv/v圖。這些數(shù)據(jù)可以用幾種圖形來表示，例如每對(duì)臺(tái)站對(duì)應(yīng)一張圖、所有臺(tái)站對(duì)平均后的圖，等等。MSNoise與各種例圖一起提供，這樣用戶可以根據(jù)需求取舍。MSNoise使用pandas模塊進(jìn)行數(shù)據(jù)分析(McKinney，2012)，matplotlib模塊(Hunter，2007)進(jìn)行靜態(tài)圖形輸出，pandas提供了強(qiáng)大的輔助函數(shù)，其中包括對(duì)時(shí)間序列數(shù)據(jù)進(jìn)行移動(dòng)平均、重采樣和去除包含非數(shù)值的序列的趨勢(shì)等。

MSNoise在設(shè)計(jì)時(shí)考慮了可插性和擴(kuò)展性。所有與數(shù)據(jù)庫(kù)、波形數(shù)據(jù)檔案或互相關(guān)函數(shù)文件的通信都要通過輔助函數(shù)完成。任何處理步驟都不直接與數(shù)據(jù)庫(kù)或數(shù)據(jù)文件相連。建議功能擴(kuò)展開發(fā)者查看函數(shù)文檔來建立其擴(kuò)展模塊。查看程序包提供的輸出例圖是一個(gè)良好的開端。我們要強(qiáng)調(diào)的是，盡管目前證明了MSNoise能夠按預(yù)期來工作(見下面的“檢驗(yàn)”一節(jié))，但人們也許會(huì)對(duì)修改或使用其他計(jì)算子程序感興趣。這可以通過更換工作流程中的某些步驟很容易地做到。

2計(jì)算

我們連續(xù)測(cè)量地球內(nèi)部地震波速變化的方法依賴于三個(gè)步驟。(1)第一步是計(jì)算每對(duì)臺(tái)站間不同時(shí)期地震噪聲時(shí)間序列的互相關(guān)函數(shù)(CCF)；(2)測(cè)量每個(gè)互相關(guān)函數(shù)與確定的參考互相關(guān)函數(shù)之間不同到達(dá)波(直達(dá)波或尾波)之間的走時(shí)延遲；(3)對(duì)不同臺(tái)站對(duì)間不同相關(guān)時(shí)間滯后的走時(shí)延遲進(jìn)行平均，并用研究區(qū)域內(nèi)相對(duì)速度均勻一致變化的簡(jiǎn)單模型來解釋這些走時(shí)延遲(δv/v=常數(shù))。

2.1　互相關(guān)函數(shù)計(jì)算

對(duì)于每個(gè)互相關(guān)的任務(wù)，對(duì)波形預(yù)處理后再進(jìn)行分析。任務(wù)是獨(dú)立的，因此幾個(gè)計(jì)算任務(wù)可以并行執(zhí)行，這通常依賴于處理器-內(nèi)存-磁盤三者速度的配合。

2.1.1波形預(yù)處理

對(duì)每個(gè)臺(tái)站，所有可能包含有關(guān)這一天數(shù)據(jù)的文件都被打開。數(shù)據(jù)道被合并和拆分，以獲得最為連續(xù)的波形段。之后不同的波形段被去平均、尖滅處理，并被重新合并為一天長(zhǎng)度的波形。如果長(zhǎng)度小于一天，則補(bǔ)零。如果長(zhǎng)度大于一天，則波形被截?cái)酁橐惶?。隨后，對(duì)每天的波形進(jìn)行低通、高通濾波，然后進(jìn)行抽樣或降采樣。抽樣較快，但僅允許以整數(shù)因子進(jìn)行抽樣，而降采樣支持任何因子，可使用不同采樣率臺(tái)站的配置。抽樣/降采樣是可配置的，建議用戶對(duì)兩種方式均作嘗試。MSNoise中使用的重采樣方法來自音頻處理領(lǐng)域，可以提供優(yōu)質(zhì)的結(jié)果(de Castro Lopo，2013)。低通和降采樣的值必須同時(shí)適當(dāng)?shù)卦O(shè)置。

2.1.2處理

所有波形一旦被加載到內(nèi)存里，就對(duì)所有波形片段(其默認(rèn)時(shí)長(zhǎng)為30分鐘，可配置)根據(jù)不同臺(tái)站對(duì)、不同分量和不同的濾波器進(jìn)行循環(huán)計(jì)算。

R=Ncos(Az)+Esin(Az)

T=Nsin(Az)-Ecos(Az)

(2)

徑向(R)和切向(T)分量由旋轉(zhuǎn)東(E)和北(N)向分量計(jì)算得到。旋轉(zhuǎn)角度由兩個(gè)傳感器間的方位角(Az)確定。

為了減弱地震事件(地方震、區(qū)域地震或全球地震)產(chǎn)生的寄生信號(hào)，在第一步處理中，我們對(duì)波形應(yīng)用了三倍于RMS的Windsorizing去離群點(diǎn)處理(Tukey，1962)。第二步是對(duì)兩個(gè)給定的頻率之間的信號(hào)幅度進(jìn)行白噪聲化處理。

互相關(guān)計(jì)算在頻率域內(nèi)進(jìn)行。令x(t)和y(t)為兩個(gè)時(shí)間序列，X(f)和Y(f)是它們的傅里葉變換。相關(guān)運(yùn)算被定義為：

(3)

式中，X*(f)是X(f)的復(fù)共軛。如果x(t)=y(t)，則上述運(yùn)算稱為自相關(guān)?；ハ嚓P(guān)函數(shù)c(t)是C(f)的傅里葉反變換。如果配置中設(shè)定(非默認(rèn)值)，每30分鐘波形片段進(jìn)行互相關(guān)獲得的c(t)就可以被保存在計(jì)算機(jī)上的產(chǎn)出文件夾里。同樣，如配置(為默認(rèn)值)，由非零、有限、非空的48個(gè)30分鐘的c(t)疊加成的一天的結(jié)果也會(huì)被保留下來。

2.2　疊加策略和輸出

MSNoise能夠使用由絕對(duì)或相對(duì)時(shí)間范圍確定的參考函數(shù)(REF)。例如，絕對(duì)范圍可能是從2010年1月1日到2011年12月31日，而相對(duì)范圍可能是最近的200天。為了確定這個(gè)參考函數(shù)，我們建議畫出已計(jì)算得到的所有互相關(guān)函數(shù)，并檢查它們的一致性。一旦識(shí)別出一個(gè)足夠穩(wěn)定的時(shí)間范圍，它便可被確定為參考函數(shù)，用于后續(xù)分析。

通過給定參考函數(shù)，可以計(jì)算每一天的互相關(guān)函數(shù)與參考函數(shù)的相關(guān)系數(shù)；這已是一個(gè)很好的指示器，用于說明某個(gè)傳感器下方、周圍或者本身是否有事情發(fā)生(圖3)。如圖3所示，正負(fù)時(shí)間滯后的信號(hào)有時(shí)會(huì)不同(Mordretetal，2010)。速度的變化會(huì)引起相關(guān)性的平滑變化，而地震噪聲源或一些散射體位置的變化應(yīng)引起相關(guān)性顯著的改變(Wegler and Sens-Sch?nfelder，2007)。

由每一天的互相關(guān)函數(shù)與參考函數(shù)的比較計(jì)算出的速度變化，會(huì)因互相關(guān)函數(shù)中的噪聲成份而出現(xiàn)強(qiáng)烈起伏。將越來越多天的互相關(guān)疊加起來，再與參考函數(shù)作互相關(guān)，有助于確定使互相關(guān)系數(shù)達(dá)到0.96的合理的天數(shù)。使用過長(zhǎng)動(dòng)窗進(jìn)行疊加的缺點(diǎn)是分辨率隨即降低。Liu等(2010)已證明窗口大小的選擇會(huì)突顯或隱藏?cái)?shù)據(jù)中的特征，因此我們使MSNoise能夠輸出多種長(zhǎng)度動(dòng)窗疊加，例如2天、5天、10天和30天的疊加。一旦用戶確定了參考函數(shù)和一個(gè)(或若干)合理的動(dòng)窗平均結(jié)果，在計(jì)算相對(duì)走時(shí)變化之前，兩者都需要輸出。同樣要注意對(duì)參考函數(shù)的確定。如果將所有歸檔數(shù)據(jù)的平均作為參考，但是實(shí)際上它包含了強(qiáng)烈的速度變化，那么參考函數(shù)的質(zhì)量將會(huì)大大降低(Duputeletal，2009；Mordretetal，2010)。強(qiáng)烈推薦在數(shù)據(jù)總時(shí)長(zhǎng)度內(nèi)調(diào)查比較穩(wěn)定的時(shí)間范圍。

只有全新的或更新過的那些天的數(shù)據(jù)才需要輸出。如果前一天內(nèi)任何互相關(guān)(CC)的任務(wù)已被標(biāo)為“Done”(完成)，那么將計(jì)算疊加，并且會(huì)向數(shù)據(jù)庫(kù)中插入新的走時(shí)變化(DTT)任務(wù)。注意，只對(duì)有新的或修訂的互相關(guān)的日期，或已確定出新的參考函數(shù)的日期，才會(huì)計(jì)算后續(xù)的δt/t。如果定義了動(dòng)窗的參考函數(shù)，那么每次都要計(jì)算參考函數(shù)的疊加，并且也要完全重新計(jì)算走時(shí)變化。對(duì)于一固定的參考函數(shù)，僅需要計(jì)算前一天的δt/t；這使得整個(gè)處理過程更快。

正如Hadziioannou等(2009)所述，地震波速變化測(cè)量的準(zhǔn)確性依賴于時(shí)間域互相關(guān)函數(shù)的穩(wěn)定重構(gòu)。這樣，作為第一步，必須確保有足夠長(zhǎng)的噪聲時(shí)間序列進(jìn)行互相關(guān)，目的是使每對(duì)臺(tái)站的單個(gè)和參考互相關(guān)函數(shù)的相干性達(dá)到一定的閾值。一個(gè)明顯的折中是時(shí)間序列越長(zhǎng)，相干性就越高，而時(shí)間分辨率就越低。

2.3　相對(duì)走時(shí)變化

2.3.1時(shí)間延遲測(cè)量

可用于估計(jì)時(shí)間延遲的技術(shù)有兩種：動(dòng)窗互譜分析(MWCS；最早由Ratdomopurbo and Poupinet，1995提出)和被動(dòng)干涉法(Sens-Sch?nfelder and Wegler，2006)。盡管有兩項(xiàng)研究都強(qiáng)調(diào)了第二種技術(shù)的穩(wěn)定性超過第一種(Duputeletal，2009；Hadziioannouetal，2009)，但是動(dòng)窗互譜分析技術(shù)的優(yōu)勢(shì)在于頻率域內(nèi)的操作可以很明確地確定相關(guān)函數(shù)中相干信號(hào)的帶寬(Clarkeetal，2011)。而且，Clarke等(2011)最近研制的方法已經(jīng)對(duì)動(dòng)窗互譜分析作了一些改進(jìn)，尤其是對(duì)統(tǒng)計(jì)參數(shù)的估計(jì)，它可以評(píng)估變化的程度，在何種程度上對(duì)擾動(dòng)的觀測(cè)是可靠的。在MSNoise中，按照Clarke等(2011)的方法計(jì)算時(shí)間延遲，為保證程序包的一致性，用Python語言進(jìn)行了重寫。在文中的“檢驗(yàn)”一節(jié)，我們比較了MSNoise和使用Daniel Clarke測(cè)量的FORTRAN程序(私人通訊，2011)獲得的結(jié)果。

圖3　UnderVolc計(jì)劃中任意兩個(gè)臺(tái)站(YA.UV02和YA.UV05)不同時(shí)期的互相關(guān)函數(shù)(CCF)(上圖)。水平虛線代表正負(fù)時(shí)間滯后的范圍，這個(gè)范圍確定了與參考函數(shù)(此處是所有可用數(shù)據(jù)的疊加)比較的兩個(gè)區(qū)域。通過比較不同動(dòng)窗疊加的互相關(guān)函數(shù)和參考函數(shù)來確定相關(guān)系數(shù)(下面4幅圖)。富爾奈斯火山的噴發(fā)用紅色垂直條帶表示

將當(dāng)前的互相關(guān)函數(shù)與參考函數(shù)進(jìn)行了比較。兩個(gè)時(shí)間序列被劃分成若干存在重疊部分的窗。在經(jīng)傅里葉變換轉(zhuǎn)換到頻率域之前，對(duì)每個(gè)片段做去平均和余弦尖滅(在兩端進(jìn)行1%的尖滅)處理。Fcur(ν)和Fref(ν)是傅里葉變換后埃爾米特對(duì)稱的頻譜的前一半，長(zhǎng)度被補(bǔ)零至等于下一個(gè)2的整數(shù)冪的長(zhǎng)度?；プVX(ν)被定義為：

(4)

式中，星號(hào)表示復(fù)共軛。之后X(ν)通過與一個(gè)漢寧窗卷積進(jìn)行平滑。兩個(gè)時(shí)間序列的相似程度可用頻率域中兩個(gè)能量密度間的互相干性來評(píng)估：

(5)

(6)

對(duì)每一處理窗口，兩個(gè)信號(hào)間的時(shí)移是研究頻帶內(nèi)采樣加權(quán)線性回歸后獲得的斜率m。不像Poupinet等(1984)，這里所說的權(quán)重由Clarke等(2011)提出，包含了互譜的振幅和互相干性。誤差通過使用權(quán)重(也就是相干性)和擬合差的平方進(jìn)行估計(jì)：

(7)

(8)

這個(gè)過程的輸出是一張表格，對(duì)于每個(gè)動(dòng)窗都包含中心時(shí)間滯后、測(cè)量的延遲及其誤差，以及該段的平均相干性。表格以文本格式[逗號(hào)分隔值(csv)]保存至硬盤上。對(duì)于任何給定的天，MSNoise計(jì)算M組(公式1)臺(tái)站對(duì)的延遲時(shí)間和滯后時(shí)間，并把結(jié)果保存到文本文件中。

2.3.2波速變化計(jì)算

如果假設(shè)相對(duì)速度變化δv/v在空間上是均勻的，并且參考函數(shù)與當(dāng)前相關(guān)函數(shù)之間存在相對(duì)時(shí)移δt/t，那么有(Ratdomopurbo and Poupinet，1995)：

(9)

以下各節(jié)中，我們將主要討論δt/t，因?yàn)檫@是MSNoise保存下來的原始信息，并且必要時(shí)可以很容易地將其轉(zhuǎn)換成δv/v。

上一步保存的動(dòng)窗互譜分析延遲時(shí)間表可作為延遲矩陣來可視化(圖4)，每天一張圖顯示了臺(tái)站對(duì)之間微小的延遲差異是關(guān)于時(shí)間滯后的函數(shù)。MSNoise不僅能夠計(jì)算每對(duì)臺(tái)站的δt/t，也可以計(jì)算濾波后所有臺(tái)站對(duì)加權(quán)平均的δt/t。當(dāng)所有臺(tái)站對(duì)被載入后，矩陣的底部會(huì)額外增加一行。這一行在MSNoise中叫做“ALL”(全部)行，包含了每一列的延遲時(shí)間在濾波后的加權(quán)平均值。對(duì)于每一列(滯后時(shí)間)，延遲時(shí)間必須滿足一定的規(guī)則才能參與加權(quán)平均計(jì)算。在“檢驗(yàn)”一節(jié)給出了對(duì)于UnderVolc火山數(shù)據(jù)中數(shù)據(jù)取舍規(guī)則的例子。

圖4　UnderVolc計(jì)劃期間2010年10月12日所有可用臺(tái)站對(duì)的延遲、誤差、相位相干性、數(shù)據(jù)選擇和延遲時(shí)間變化矩陣的實(shí)例。矩陣的每一行相當(dāng)于一對(duì)臺(tái)站，而前4個(gè)矩陣的每一列是動(dòng)窗互譜分析處理過程中動(dòng)窗的中心時(shí)間滯后。每個(gè)矩陣都有自己的色標(biāo)。矩陣底部的行是對(duì)所在列進(jìn)行濾波后再加權(quán)平均的結(jié)果(詳見正文)，它在MSNoise中被稱為“ALL”行?！癆LL”行僅在時(shí)間滯后-left_maxlag：-left_minlag和right_minlag：right_maxlag(虛線)之間才被計(jì)算。由選中的數(shù)據(jù)通過加權(quán)最小二乘回歸(強(qiáng)制或沒有強(qiáng)制穿過原點(diǎn))確定的每個(gè)δt/t值被表示為畫有誤差棒的離散點(diǎn)。紅色豎線是由“ALL”行確定的δt/t(不強(qiáng)制過原點(diǎn)：0.074%±0.009%，強(qiáng)制：0.074%±0.008%)，而綠線是所有δt/t的加權(quán)平均(沒有強(qiáng)制：0.081%±0.065%，強(qiáng)制：0.074%±0.069%)。強(qiáng)制和沒有強(qiáng)制兩種方式獲得的δt/t值非常接近，但是后者的誤差比前者大一個(gè)數(shù)量級(jí)

圖5　圖4的實(shí)例矩陣中δt/t值以離散點(diǎn)的方式呈現(xiàn)。本圖中，每一行或每一列都代表唯一的臺(tái)站，因此每個(gè)單元格就代表一對(duì)臺(tái)站。單元的顏色依賴于每對(duì)臺(tái)站計(jì)算出的δt/t(斜率)。這里，對(duì)角線為空，是因?yàn)闆]有計(jì)算自相關(guān)。根據(jù)互易定理，每個(gè)臺(tái)站對(duì)的互相關(guān)函數(shù)僅需計(jì)算一次

對(duì)于計(jì)算時(shí)移的最小和最大滯后時(shí)間的選擇是十分重要的。由于沿較長(zhǎng)路徑傳播的散射波會(huì)積累較長(zhǎng)的時(shí)間延遲，所以在火山環(huán)境中開展的大部分研究使用了互相關(guān)函數(shù)的晚至部分(例如，Duputeletal，2009，5~20s或Mordretetal，2010，10~30s)。最近Anggono等(2012)的研究是個(gè)例外，他們?cè)谘芯恐惺褂昧恕?0s到±1s的滯后時(shí)間。由于尾波部分(>10s)的波形相似度較差，他們既沒使用動(dòng)窗互譜分析方法，也沒使用拉伸方法，而是比較了每天互相關(guān)函數(shù)和參考互相關(guān)函數(shù)的最大值。我們建議用戶預(yù)先測(cè)試不同的窗口，并且避免在零時(shí)滯附近操作。而且，互相關(guān)函數(shù)僅在火山區(qū)，特別是較短滯后時(shí)間處(即由Duputeletal，2009描述的彈道波)對(duì)稱。這或許說明了源分布不均勻，或者在某一特定位置存在次級(jí)噪聲源的散射體(Paul，2005；Stehlyetal，2006)，亦或反映出有時(shí)鐘問題(Sens-Sch?nfelder，2008)。Mordret等(2010)平均了互相關(guān)函數(shù)的因果和非因果部分，得到重構(gòu)的波場(chǎng)有更好的相干性和穩(wěn)定性。然而，同樣也可以決定僅對(duì)重構(gòu)互相關(guān)函數(shù)中因果的或非因果的部分進(jìn)行單獨(dú)處理。由于背景噪聲源會(huì)隨時(shí)間變化(其頻譜、位置和振幅均有變化)，因此地震波速的變化會(huì)存在偏差。噪聲源性質(zhì)上的這些不同主要影響重構(gòu)直達(dá)波。因而，為了估計(jì)這種偏差，可以初步只考慮測(cè)量互相關(guān)函數(shù)直達(dá)波或早至波的走時(shí)時(shí)移，從而測(cè)量地震波速的變化。

接下來對(duì)矩陣的M+1行通過滿足要求的點(diǎn)(圖4中數(shù)據(jù)選擇矩陣中黑色部分)進(jìn)行加權(quán)最小二乘回歸(WLS)計(jì)算。加權(quán)最小二乘回歸計(jì)算進(jìn)行兩次，第一次允許有非零常值截距(公式10a)，第二次要求強(qiáng)制回歸經(jīng)過原點(diǎn)([0，0])(公式10b)。計(jì)算每個(gè)參數(shù)的誤差，在公式(10a)中，a和m對(duì)應(yīng)的誤差分別是ea和em；在公式(10b)中，m0對(duì)應(yīng)的誤差是em0。存在較大的非零常數(shù)a反映出儀器可能存在零點(diǎn)漂移。單天的結(jié)果可以畫成臺(tái)站對(duì)序號(hào)所對(duì)應(yīng)的δt/t的散點(diǎn)圖(圖4右邊兩幅圖)，或是表示為一個(gè)矩陣圖(圖5)。前者還展示出了由“ALL”行確定的斜率(紅色豎線)和由所有臺(tái)站對(duì)確定的斜率的加權(quán)平均值(綠色豎線)。在這個(gè)例子中，兩個(gè)結(jié)果高度相似，但是，即使從直觀上看，有些臺(tái)站對(duì)也似乎表現(xiàn)異常。矩陣(圖5)有助于識(shí)別出可以放在一起分析的臺(tái)站類群。以UnderVolc項(xiàng)目的情況為例，我們可以定義兩個(gè)臺(tái)站類群，例如“火山口臺(tái)站”和“大間距臺(tái)站”，

δt=a+mt

(10a)

δt0=m0t

(10b)

式中，下標(biāo)0用于強(qiáng)調(diào)回歸直線強(qiáng)制通過[0，0]。

對(duì)于每個(gè)臺(tái)站對(duì)和“ALL”行，這6個(gè)數(shù)值(斜率m和m0，常數(shù)a和它們的誤差)都被存至硬盤，以便繪圖程序或輸出程序很容易地讀取。

3檢驗(yàn)

為了檢驗(yàn)我們的處理流程的正確性，我們用相同的數(shù)據(jù)集進(jìn)行了測(cè)試，該數(shù)據(jù)集在UnderVolc項(xiàng)目期間被成功用于識(shí)別2010年10月和12月兩次噴發(fā)前富爾奈斯火山下的前兆速度變化。

3.1　數(shù)據(jù)和參數(shù)

數(shù)據(jù)集由布設(shè)在富爾奈斯火山區(qū)的21個(gè)寬頻帶三分量地震臺(tái)站組成，臺(tái)站間距從<1km到>10km不等(圖6)。大部分臺(tái)站于2009年11月初架設(shè)，并在2011年6月底撤除(圖2)。事實(shí)上，大部分臺(tái)址由富爾奈斯火山觀測(cè)臺(tái)網(wǎng)接管，并使用了新的臺(tái)站代碼。本次測(cè)試中，共處理了11 393個(gè)Z分量文件。為了計(jì)算互相關(guān)函數(shù)，對(duì)100Hz的波形進(jìn)行了0.01至8.0Hz的帶通濾波，然后抽樣到20Hz(因子為5)。對(duì)于不同的臺(tái)站對(duì)，僅計(jì)算了ZZ分量(未作自相關(guān))。每天的數(shù)據(jù)都被分割成每30分鐘一段，產(chǎn)出的互相關(guān)函數(shù)與按天疊加的互相關(guān)函數(shù)被一起存入硬盤。按2天、5天、10天和30天動(dòng)窗進(jìn)行疊加計(jì)算。動(dòng)窗互譜分析的頻帶是0.2~0.85Hz，以10s動(dòng)窗50%重疊計(jì)算。對(duì)于δt/t計(jì)算，都使用了因果和非因果部分，滯后時(shí)間選在±5s和±50s之間。其余的選擇參數(shù)為：最小相干性是0.5，最大誤差是0.1s，最大的dt是0.5s。表1和表2給出了在4v CPU的虛擬機(jī)上運(yùn)行的時(shí)間和數(shù)據(jù)量的示例。

3.2　結(jié)果

應(yīng)用相同的計(jì)算延遲時(shí)間的方法，分別利用Clarke等(2011)的FORTRAN程序和MSNoise計(jì)算了UnderVolc項(xiàng)目期間火山的δv/v變化，兩者的結(jié)果非常相似(圖7)。大部分時(shí)間兩者的差異較小，明顯低于0.02%，而火山噴發(fā)時(shí)要高于0.05%?？傮w趨勢(shì)是存在的，并明顯表現(xiàn)為每年速度有0.2%的微小變化?？梢匀コ鼴renguier，Shapiro等(2008)意義上的這一長(zhǎng)期成分(LTV)，以突出更短時(shí)期的變化。深入地分析這兩種程序，我們認(rèn)為微小的差異可能來自余弦平滑函數(shù)的構(gòu)造方式或是相位計(jì)算程序。我們選擇使用了numpy，scipy(Jonesetal，2001；Oliphant，2006)和statsmodels(Seabold and Perktold，2010)程序包中久經(jīng)考驗(yàn)的數(shù)值計(jì)算子程序，而沒有重寫Clarke的FORTRAN程序中的相關(guān)函數(shù)。

圖6　UnderVolc計(jì)劃的臺(tái)站分布圖。白色三角形是臺(tái)站，各臺(tái)站對(duì)間的路徑以白線表示，并覆蓋在富爾奈斯火山的數(shù)字高程模型(DEM)上。沒有數(shù)字高程模型信息的臺(tái)站分布圖是使用腳本plot_station_map.py自動(dòng)生成的

步驟運(yùn)行時(shí)間安裝10分鐘配置5分鐘掃描數(shù)據(jù)檔案30分鐘識(shí)別102626個(gè)新任務(wù)30分鐘處理互相關(guān)任務(wù)(1個(gè)濾波器)21小時(shí)疊加獲得210對(duì)臺(tái)站的參考函數(shù)(REF)33分鐘按4種動(dòng)窗對(duì)所有臺(tái)站對(duì)疊加6小時(shí)計(jì)算動(dòng)窗互譜(MWCS)8小時(shí)計(jì)算δt/t1小時(shí)

表2　數(shù)據(jù)檔案和MSNoise輸出的數(shù)據(jù)量

2009年11月、2010年1月、2011年10月和12月的噴發(fā)之前δv/v都減小(圖8中向上的部分)。在2009年10月噴發(fā)前似乎沒有前兆信號(hào)。我們可以看到，2009年11月、2010年10月和12月噴發(fā)的前兆時(shí)間大約為15天。對(duì)于后兩次噴發(fā)，在真正噴發(fā)前1天的δv/v增大很明顯。這可以解釋為深部速度的增大，也可解釋為噴發(fā)前幾小時(shí)內(nèi)強(qiáng)烈的地震活動(dòng)或是數(shù)分鐘內(nèi)出現(xiàn)的噴發(fā)振顫。最后兩種情況導(dǎo)致了噪聲源區(qū)的變化，并干擾了這個(gè)系統(tǒng)。值得注意的是，移動(dòng)疊加函數(shù)在噴發(fā)前5天與參考函數(shù)的相關(guān)程度變?nèi)?圖3)。正如Liu等(2010)所指出的，改變疊加持續(xù)時(shí)間的長(zhǎng)度可能會(huì)在δv/v曲線演化中突出或去除信號(hào)。圖3中以2天、5天、10天或30天疊加的互相關(guān)函數(shù)與同樣的參考函數(shù)比較時(shí)，均出現(xiàn)了相關(guān)程度變?nèi)醯那闆r，但表現(xiàn)不同。例如，2009年11月噴發(fā)前δv/v的前兆變化在不超過10天的疊加持續(xù)時(shí)間是可見的，卻在30天疊加時(shí)完全被去除(圖8)。短周期(等于或少于幾天)的變化可在1天、2天和5天的疊加中看到(圖8)，可能具有不同的來源，例如與“局部”原因有關(guān)，如氣壓、降雨或地震活動(dòng)的變化，或者還有像噪聲源位置變化這樣的“區(qū)域”原因。雖然都有可能，但幾天的短周期使“局部”原因比“區(qū)域”原因更為可能。

4結(jié)論和展望

我們給出了第一個(gè)使用背景噪聲計(jì)算和監(jiān)測(cè)相對(duì)波速變化的完整軟件包。MSNoise是一個(gè)集成的解決方案，當(dāng)執(zhí)行預(yù)定任務(wù)時(shí)，可以自動(dòng)掃描數(shù)據(jù)檔案，確定哪些任務(wù)需要處理。整個(gè)軟件包均由Python語言編寫，是開源的，并可通過網(wǎng)頁http：//www.msnoise.org(最后訪問時(shí)間2014年3月)免費(fèi)獲得。其處理步驟與在UnderVolc計(jì)劃期間(Brenguieretal，2012)成功識(shí)別富爾奈斯火山(法國(guó)留尼旺島)下前兆速度變化的另一程序(Clarkeetal，2011)進(jìn)行了比較測(cè)試。對(duì)于相同的數(shù)據(jù)集，兩種處理程序的結(jié)果高度相似。MSNoise是可插入和擴(kuò)展的；對(duì)于研究者來說，這實(shí)際上是一個(gè)十分有用的工具，因?yàn)樗麄冎恍桕P(guān)注實(shí)現(xiàn)他們自己的處理程序，而這得益于該系統(tǒng)穩(wěn)健的架構(gòu)。

圖7　(a)UnderVolc計(jì)劃期間富爾奈斯火山地區(qū)的相對(duì)地震波速變化。綠色曲線為使用D.Clarke的FORTRAN程序計(jì)算的結(jié)果，紅色曲線為使用我們的程序計(jì)算的結(jié)果，紅色垂直條帶表示噴發(fā)。(b)紅色曲線為兩種結(jié)果之間的絕對(duì)差異，用百分比表示，藍(lán)色曲線為項(xiàng)目中每天可用的臺(tái)站對(duì)數(shù)

圖8　使用MSNoise對(duì)5種不同動(dòng)窗疊加(1，2，5，10和30天)計(jì)算的δv/v的結(jié)果去趨勢(shì)后進(jìn)行比較，且在δt/t估計(jì)程序中沒有強(qiáng)制斜率估計(jì)穿過原點(diǎn)(公式10)。富爾奈斯火山的噴發(fā)表示為紅色垂直條帶。整個(gè)臺(tái)網(wǎng)的δv/v，即“ALL”行(紅色)，與所有單對(duì)臺(tái)站的加權(quán)平均(綠色)匹配。UnderVolc計(jì)劃期間，5次噴發(fā)中的4次都以δv/v的變化為前兆(見正文)

我們的方法基于下述假設(shè)，即重構(gòu)的互相關(guān)函數(shù)尾波部分的每個(gè)波包在其傳播路徑上都經(jīng)歷了均勻一致的相對(duì)波速變化(δv/v為常數(shù))。這是一個(gè)強(qiáng)假設(shè)，可能不是在任何情況下都有效。其結(jié)果是測(cè)量得到的相對(duì)地震波速變化可能與研究介質(zhì)中“真實(shí)的”δv/v有很大不同。主要的局限性產(chǎn)生于尾波波包可能經(jīng)過非均勻的速度變化區(qū)域這一事實(shí)。在這種情況下，測(cè)量的δv/v就是真實(shí)δv/v的平均結(jié)果；例如，鑒于δv/v的值在空間范圍內(nèi)是常數(shù)這個(gè)假設(shè)，人們可能在δv/v增大時(shí)錯(cuò)誤地解釋波速的變化δv。

最后，盡管我們給出的對(duì)MSNoise的有效檢驗(yàn)是在火山環(huán)境下進(jìn)行的，但其應(yīng)用不應(yīng)當(dāng)局限于這種特殊情況，我們期望看到它在不同環(huán)境中得到應(yīng)用，并期望得到來自其他研究小組的評(píng)論意見和貢獻(xiàn)。例如，這些貢獻(xiàn)可能是一種新的波形相關(guān)算法、新的δt/t計(jì)算技術(shù)、衰減計(jì)算的實(shí)現(xiàn)(如在Prietoetal，2011中)、延遲矩陣在定位δv/v變化中的運(yùn)用，或者在不同時(shí)間滯后上的解相關(guān)計(jì)算(如在Laroseetal，2010；Obermannetal，2013中)。其他方面的貢獻(xiàn)可能包括分析更高頻數(shù)據(jù)(高于1kHz，如在Olivieretal，2012中)，或者將互相關(guān)函數(shù)輸入到頻率-時(shí)間分析(FTAN)(如Levshinetal，1989中描述的)，用于計(jì)算頻散和隨后的一維面波速度變化的反演。我們正計(jì)劃實(shí)現(xiàn)定義多種參考函數(shù)以及“臺(tái)站分類”的可能性；前者就像定義多種疊加方式一樣，后者則類似于定義“ALL”行(見“波速變化計(jì)算”一節(jié))。

參考文獻(xiàn)

Anggono，T.，T.Nishimura，H.Sato，H.Ueda，and M.Ukawa (2012).Spatio-temporal changes in seismic velocity associated with the 2000 activity of Miyakejima volcano as inferred from cross-correlation analyses of ambient noise，J.Volcanol.Geoth.Res.247/248，93-107，doi：10.1016/j.jvolgeores.2012.08.001.

Beyreuther，M.，R.Barsch，L.Krischer，T.Megies，Y.Behr，and J.Wassermann (2010).ObsPy：A Python toolbox for seismology，Seismol.Res.Lett.81，no.3，530-533，doi：10.1785/gssrl.81.3.530.

Brenguier，F(xiàn).，M.Campillo，C.Hadziioannou，N.M.Shapiro，R.M.Nadeau，and E.Larose (2008).Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations，Science321，no.5895，1478-1481，doi：10.1126/science.1160943.

Brenguier，F(xiàn).，D.Clarke，Y.Aoki，N.M.Shapiro，M.Campillo，and V.Ferrazzini(2011).Monitoring volcanoes using seismic noise correlations，ComptesRendusGeosci.343，nos.8/9，633-638，doi：10.1016/j.crte.2010.12.010.

Brenguier，F(xiàn).，P.Kowalski，T.Staudacher，V.Ferrazzini，F(xiàn).Lauret，P.Boissier，P.Catherine，A.Lemarchand，C.Pequegnat，O.Meric，C.Pardo，A.Peltier，S.Tait，N.M.Shapiro，M.Campillo，and A.Di Muro(2012).First results from the UnderVolc high resolution seismic and GPS network deployed on Piton de la Fournaise Volcano，Seismol.Res.Lett.83，no.1，97-102，doi：10.1785/gssrl.83.1.97.

Brenguier，F(xiàn).，N.Shapiro，M.Campillo，V.Ferrazzini，Z.Duputel，O.Coutant，and A.Nercessian(2008).Towards forecasting volcanic eruptions using seismic noise，NatureGeosci.1，no.2，126.

Campillo，M.，and A.Paul(2003).Long-range correlations in the diffuse seismic coda，Science299，no.5606，547-549，doi：10.1126/science.1078551.

Clarke，D.，L.Zaccarelli，N.M.Shapiro，and F.Brenguier(2011).Assessment of resolution and accuracy of the Moving Window Cross Spectral technique for monitoring crustal temporal variations using ambient seismic noise，Geophys.J.Int.186，no.2，867-882，doi：10.1111/j.1365-246X.2011.05074.x.

de Castro Lopo，E.(2013).SecretRabbitCode(akalibsamplerate)，http：//www.mega.nerd.com/SRC/ (last accessed March 2014).

Duputel，Z.，V.Ferrazzini，F(xiàn).Brenguier，N.Shapiro，M.Campillo，and A.Nercessian (2009).Real time monitoring of relative velocity changes using ambient seismic noise at the Piton de la Fournaise volcano(La Réunion)from January 2006 to June 2007，J.Volcanol.Geoth.Res.184，no.1/2，164-173，doi：10.1016/j.jvolgeores.2008.11.024.

Enthought (2008).TheEnthoughtPythonDistribution，http：//code.enthought.com (last accessed March 2014).

Goldstein，P.，J.A.Snoke，M.Firpo，and L.Minner (2003).SAC2000：Signal processing and analysis tools for seismologists and engineers，inTheIASPEIInternationalHandbookofEarthquakeandEngineeringSeismology，W.H.K.Lee，H.Kanamori，P.C.Jennings，and C.Kisslinger(Editors)，Academic press，London，2000 pp.

Hadziioannou，C.，E.Larose，O.Coutant，P.Roux，and M.Campillo(2009).Stability of monitoring weak changes in multiply scattering media with ambient noise correlation：Laboratory experiments，J.Acoust.Soc.Am.125，no.6，3688-3695.

Herrmann，R.(2002).ComputerProgramsinSeismology3.3，http：//www.eas.slu.edu/eqc/eqccps.html (last accessed March 2014).

Hunter，J.(2007).Matplotlib：A 2D graphics environment，Comput.Sci.Eng.9，no.3，90-95，doi：10.1109/MCSE.2007.55.

Jones，E.，T.Oliphant，P.Peterson，and SciPy deve-lopment team(2001).SciPy：OpensourcescientifictoolsforPython，http：//www.scipy.org/(last accessed March 2014).

Larose，E.，T.Planes，V.Rossetto，and L.Margerin (2010).Locating a small change in a multiple scattering environment，Appl.Phys.Lett.96，no.20，204，101-204，101-3.

Levshin，A.L.，T.Yanovskaya，A.Lander，B.Bukshin，M.P.Barmin，L.Ratnikova，and E.Its (1989).SeismicSurfaceWavesinaLaterallyInhomogeneousEarth，V.I.Keilis-Borok (Editor)，Kluwer Academic Publishers，Norwell，Massachusetts，293 pp.

Liu，Z.，J.Huang，and J.Li (2010).Comparison of four techniques for estimating temporal change of seismic velocity with passive image interfero-metry，Earthq.Sci.23，no.5，511-518，doi：10.1007/s11589-010-0749-z.

McKinney，W.(2012).PythonforDataAnalysis，O’Reilly Media，Cambridge，470 pp.

Megies，T.，M.Beyreuther，R.Barsch，L.Krischer，and J.Wassermann(2011).ObsPy—What can it do for data centers and observatories?Ann.Geophys.54，no.1，47-58，doi：10.4401/ag-4838.

Mordret，A.，A.Jolly，Z.Duputel，and N.Fournier(2010).Monitoring of phreatic eruptions using interferometry on retrieved cross-correlation function from ambient seismic noise：Results from Mt.Ruapehu，New Zealand，J.Volcanol.Geoth.Res.191，nos.1/2，46-59，doi：10.1016/j.jvolgeores.2010.01.010.

Obermann，A.，T.Planes，E.Larose，C.Sens-Sch?nfelder，and M.Campillo(2013).Depth sensitivity of seismic coda waves to velocity perturbations in an elastic heterogeneous medium，Geophys.J.Int.194，no.1，372-382，doi：10.1093/gji/ggt043.

Oliphant，T.E.(2006).GuidetoNumPy，Brigham Young University，Provo，Utah，371 pp.

Olivier，G.，F(xiàn).Brenguier，M.Campillo，R.Lynch，and N.Shapiro(2012).Monitoring of velocity variations in an underground mine using ambient seismic noise，inAbstractBookofAGUFallMeeting2012.

Paul，A.(2005).Empirical synthesis of time-asymmetrical green functions from the correlation of coda waves，J.Geophys.Res.110，no.B8，doi：10.1029/2004JB003521.

Poupinet，G.，W.L.Ellsworth，and J.Frechet(2011).Monitoring velocity variations in the Crust using earthquake doublets：An application to the Calaveras fault，California，J.Geophys.Res.89，no.B7，5719-5731.

Prieto，G.A.，M.Denolle，J.F.Lawrence，and G.C.Beroza(2011).On amplitude information carried by the ambient seismic field，ComptesRendusGeosci.343，nos.8/9，600-614，doi：10.1016/j.crte.2011.03.006.

Ratdomopurbo，A.，andG.Poupinet (1995).Monitoring a temporal change of seismic velocity in a volcano：Application to the 1992 eruption of Mt.Merapi(Indonesia)，Geophys.Res.Lett.22，no.7，775-778.

Seabold，J.，and J.Perktold (2010).Statsmodels：Econometric and statistical modeling with python，in S.van der Walt and J.Millman(Editors)，Proc.ofthe9thPythoninScienceConference，Austin，Texas，http：//conference.scipy.org/proceedings/scipy2010/(last accessed March 2014).

Sens-Sch?nfelder，C.(2008).Synchronizing seismic networks with ambient noise，Geophys.J.Int.174，no.3，966-970.

Sens-Sch?nfelder，C.，and U.Wegler(2006).Passive image interferometry and seasonal variations of seismic velocities at Merapi volcano(Indonesia)，Geophys.Res.Lett.33，L21302，doi：10.1029/2006GL027797.

Sens-Sch?nfelder，C.，and U.Wegler(2011).Passive image interferometry for monitoring crustal changes with ambient seismic noise，ComptesRendusGeosci.343，nos.8/9，639-651，doi：10.1016/j.crte.2011.02.005.

Stehly，L.，M.Campillo，and N.M.Shapiro(2006).A study of the seismic noise from its long-range correlation properties，J.Geophys.Res.111，no.B10，doi：10.1029/2005JB004237.

Stehly，L.，M.Campillo，and N.M.Shapiro(2007).Traveltime measurements from noise correlation：Stability and detection of instrumental time-shifts，Geophys.J.Int.171，no.1，223-230，doi：10.1111/j.1365-246X.2007.03492.x.

Tukey，J.W.(1962).The future of data analysis，Ann.Math.Stat.33，no.1，1-67，doi：10.1214/aoms/1177704711.

Wegler，U.，and C.Sens-Sch?nfelder(2007).Fault zone monitoring with passive image interferome-try，Geophys.J.Int.168，no.3，1029-1033，doi：10.1111/j.1365-246X.2006.03284.x.

譯自：Seismol.Res.Lett.2014.85(3)：715-725

原題：MSNoise，a Python Package for Monitoring Seismic Velocity Changes Using Ambient Seismic Noise

(中國(guó)地震臺(tái)網(wǎng)中心安艷茹、張瑩瑩譯；馬延路、鄒立曄校；呂春來復(fù)校)

安艷茹(1986—)，女，中國(guó)地震臺(tái)網(wǎng)中心工程師，主要從事地震監(jiān)測(cè)和地震背景噪聲方面的研究，Tel：59959348，E-mail：anyanru8@126.com。

譯 者 簡(jiǎn) 介

doi：10.16738/j.cnki.issn.1003-3238.201504005