MSNoise：利用地震背景噪声监测地震波速变化的Python程序包

MSNoise：利用地震背景噪声监测地震波速变化的Python程序包

Thomas LecocqCorentin CaudronFlorent Brenguier

0引言

最近，从成像到监测，地震背景噪声在不同领域都证明了它的效能。两个传感器之间的脉冲响应[或格林函数(GF)]可以通过所记录地震噪声的互相关加以重构(Campillo and Paul，2003)。从全球到区域范围或局部地点的地球内部成像来看，这个方法都给出了非常出色的结果。近来，这个方法又被拓展用于研究这些格林函数随时间的变化。其走时延迟随时间的变化可能源于介质速度的改变，也可能源于一个/多个散射体或噪声源的位置的显著移动。一些利用地震背景噪声的研究揭示出，火山体内部结构的小扰动可用地震波性质的变化来探测(Sens-Schönfelder and Wegler，2006；Brenguier，Shapiro，etal，2008；Duputeletal，2009；Mordretetal，2010；Brenguieretal，2011；Anggonoetal，2012)。与使用主动源或地震尾波方法不同的是：这个技术使用地震台站24小时不间断的记录对介质进行连续采样。这个方法已经证明了它能够发现在断层带(Wegler and Sens-Schönfelder，2007；Brenguier，Campilloetal，2008)和月球环境(Sens-Schönfelder and Wegler，2011)下物理随时间变化的证据，或者能检测出仪器问题(Stehlyetal，2007；Sens-Schönfelder，2008)。

在诸如地震分析代码(Seismic Analysis Code，SAC)(Goldsteinetal，2003)或地震学计算程序(Computer Programs in Seismology，CPS)3.3版(Herrmann，2002)中，都有一些计算地震噪声互相关的程序。据我们所知，至今还没有发布一个集成的解决方案，能够对原始波形处理得到走时变化，包括自动检测归档数据中的变化，并预先设定处理计划(每小时，每天)只进行必要的计算，同时它可作为一个工具来研究归档的波形库。自动化检测很重要，因为，例如来自实时遥测台站的数据流包含断点，这些断点随后会被填补，并且数据流提供用于分析的重要数据。这种工具必须能够处理任何常用的地震格式，并且尽可能快速和优化。为了可插性和可扩展，它应该与数据档案和具有高级辅助功能的数据库进行交互。最后，它必须产出可输出的数据，可以是波形格式、表格式的文本文件，也可以是高质量的图形。这就是地震噪声用于监测(即MSNoise，Monitoring using Seismic Noise)程序的目的。

我们不愿意将MSNoise系统作为另一个黑箱提供给用户，因此它作为一个集成的解决方案，其源代码是公开的，有注释和说明文档。处理流程分成若干步骤，只要遵循某一模块的输入和输出格式，用户就可以很容易地用其他代码代替原有模块中的程序。

图1　MSNoise工作流程示意图。一次性的安装部分在新任务定义步骤之前进入到日常工作流程

在这篇文章中，我们从信息技术和科学算法的两个角度给出了软件的实现步骤，并配有图形输出的例子来进行展示和描述。最后我们使用富尔奈斯火山(法国留尼旺岛)的数据验证了MSNoise的可靠性；该火山在UnderVolc研究计划(Brenguieretal，2012)期间分别于2010年10月和12月发生了两次喷发，并在两次喷发前观测到了前兆性的地震波速变化。UnderVolc是一项隶属法国国家科研署(ANR)的研究项目，名字的缩写代表“理解火山过程(UNDERstanding VOLCanic Processes)：致力于喷发预测和风险降低的研究，应用于法国留尼旺岛的富尔奈斯火山地区(2009—2013)”。

MSNoise用Python语言编写，是完全跨平台的。由于是用于非盈利性的教育和研究，软件本身是开源免费的。MSNoise在欧盟公共许可证(EUPLv1.1)下被许可使用。未经与原著者签署协议，禁止与MSNoise相关的商业服务，诸如提供安装、维护方面的咨询或支持。MSNoise可在http：//www.msnoise.org/网页(最后访问时间2014年3月)上获得。

1软件功能

MSNoise的工作流程非常简单明了(图1)，由一些简单的步骤组成，在配套的网站上有更多的描述。

1.1　安装和配置

MSNoise需要Python和其他几个程序包，以便在任何操作系统(OS)下运行。波形归档文件必须是已知的或自定义的格式，例如已知的格式可以是SeisComP数据结构(SDS)或者是统一数据缓冲区(BUD)结构。MSNoise使用数据库来储存配置参数、波形元数据和任务。数据库可以是MySQL或者sqlite，但由于性能的原因，我们推荐使用MySQL数据库。安装脚本随程序包一起提供，用于初始化数据库连接，存储本地用户名和密码，使软件配置器以及其他处理进程能够与数据库交互。软件配置器的图形用户接口[GUI，使用Enthought工具套件(Enthought，2008)构建]显示三个不同的面板：MSNoise的一般配置、台站配置和滤波器配置。在线文档中描述了所有的配置细节。

图2　UnderVolc项目的数据可用性图：上图中灰色部分表示包含数据的日期，下图汇总了每天的可用台站数N。这幅图是利用脚本plot_data_availability.py自动绘制的

1.2　自动任务定义

1.2.1数据发现

为了每晚自动运行前一天获得的数据，MSNoise需要检查数据档案，查找新的或修改过的文件。这些文件可能是前一天获得的，也可能是从以前的离线台站获得的包含有用信息的数据，比如一个月以前的数据。搜索时间在config(配置)中定义。MSNoise使用带有-mtime参数的find(查找)命令(Windows系统下为gnufind)来定位新的或修改过的文件。一旦有文件被定位并读取[使用obspy(Beyreutheretal，2010；Me-

giesetal，2011)]，这些文件将被插入(若是新的)或者更新(若是修改过的)到data availability(数据可用性)表中。之后这张表可被转化为一张数据可用性图(图2)。脚本首次运行调用时，必须带着初始化参数init，目的是避免在find命令中使用mtime参数，这样才能将所有被发现的文件都插入到data availability(数据可用性)表中。

1.2.2任务定义

MSNoise一旦识别出所有新的或修改过的文件，就会确定需要处理哪些天哪些台站对。对于N个台站，要计算M组台站对(公式1a)。如果需要波形自相关，M则由公式(1b)来定义：

(1a)

(1b)

1.3　输出、扩展性和可插性

MSNoise最初的目的是用来提供随时间变化的δv/v图。这些数据可以用几种图形来表示，例如每对台站对应一张图、所有台站对平均后的图，等等。MSNoise与各种例图一起提供，这样用户可以根据需求取舍。MSNoise使用pandas模块进行数据分析(McKinney，2012)，matplotlib模块(Hunter，2007)进行静态图形输出，pandas提供了强大的辅助函数，其中包括对时间序列数据进行移动平均、重采样和去除包含非数值的序列的趋势等。

MSNoise在设计时考虑了可插性和扩展性。所有与数据库、波形数据档案或互相关函数文件的通信都要通过辅助函数完成。任何处理步骤都不直接与数据库或数据文件相连。建议功能扩展开发者查看函数文档来建立其扩展模块。查看程序包提供的输出例图是一个良好的开端。我们要强调的是，尽管目前证明了MSNoise能够按预期来工作(见下面的“检验”一节)，但人们也许会对修改或使用其他计算子程序感兴趣。这可以通过更换工作流程中的某些步骤很容易地做到。

2计算

我们连续测量地球内部地震波速变化的方法依赖于三个步骤。(1)第一步是计算每对台站间不同时期地震噪声时间序列的互相关函数(CCF)；(2)测量每个互相关函数与确定的参考互相关函数之间不同到达波(直达波或尾波)之间的走时延迟；(3)对不同台站对间不同相关时间滞后的走时延迟进行平均，并用研究区域内相对速度均匀一致变化的简单模型来解释这些走时延迟(δv/v=常数)。

2.1　互相关函数计算

对于每个互相关的任务，对波形预处理后再进行分析。任务是独立的，因此几个计算任务可以并行执行，这通常依赖于处理器-内存-磁盘三者速度的配合。

2.1.1波形预处理

对每个台站，所有可能包含有关这一天数据的文件都被打开。数据道被合并和拆分，以获得最为连续的波形段。之后不同的波形段被去平均、尖灭处理，并被重新合并为一天长度的波形。如果长度小于一天，则补零。如果长度大于一天，则波形被截断为一天。随后，对每天的波形进行低通、高通滤波，然后进行抽样或降采样。抽样较快，但仅允许以整数因子进行抽样，而降采样支持任何因子，可使用不同采样率台站的配置。抽样/降采样是可配置的，建议用户对两种方式均作尝试。MSNoise中使用的重采样方法来自音频处理领域，可以提供优质的结果(de Castro Lopo，2013)。低通和降采样的值必须同时适当地设置。

2.1.2处理

所有波形一旦被加载到内存里，就对所有波形片段(其默认时长为30分钟，可配置)根据不同台站对、不同分量和不同的滤波器进行循环计算。

R=Ncos(Az)+Esin(Az)

T=Nsin(Az)-Ecos(Az)

(2)

径向(R)和切向(T)分量由旋转东(E)和北(N)向分量计算得到。旋转角度由两个传感器间的方位角(Az)确定。

为了减弱地震事件(地方震、区域地震或全球地震)产生的寄生信号，在第一步处理中，我们对波形应用了三倍于RMS的Windsorizing去离群点处理(Tukey，1962)。第二步是对两个给定的频率之间的信号幅度进行白噪声化处理。

互相关计算在频率域内进行。令x(t)和y(t)为两个时间序列，X(f)和Y(f)是它们的傅里叶变换。相关运算被定义为：

(3)

式中，X*(f)是X(f)的复共轭。如果x(t)=y(t)，则上述运算称为自相关。互相关函数c(t)是C(f)的傅里叶反变换。如果配置中设定(非默认值)，每30分钟波形片段进行互相关获得的c(t)就可以被保存在计算机上的产出文件夹里。同样，如配置(为默认值)，由非零、有限、非空的48个30分钟的c(t)叠加成的一天的结果也会被保留下来。

2.2　叠加策略和输出

MSNoise能够使用由绝对或相对时间范围确定的参考函数(REF)。例如，绝对范围可能是从2010年1月1日到2011年12月31日，而相对范围可能是最近的200天。为了确定这个参考函数，我们建议画出已计算得到的所有互相关函数，并检查它们的一致性。一旦识别出一个足够稳定的时间范围，它便可被确定为参考函数，用于后续分析。

通过给定参考函数，可以计算每一天的互相关函数与参考函数的相关系数；这已是一个很好的指示器，用于说明某个传感器下方、周围或者本身是否有事情发生(图3)。如图3所示，正负时间滞后的信号有时会不同(Mordretetal，2010)。速度的变化会引起相关性的平滑变化，而地震噪声源或一些散射体位置的变化应引起相关性显著的改变(Wegler and Sens-Schönfelder，2007)。

由每一天的互相关函数与参考函数的比较计算出的速度变化，会因互相关函数中的噪声成份而出现强烈起伏。将越来越多天的互相关叠加起来，再与参考函数作互相关，有助于确定使互相关系数达到0.96的合理的天数。使用过长动窗进行叠加的缺点是分辨率随即降低。Liu等(2010)已证明窗口大小的选择会突显或隐藏数据中的特征，因此我们使MSNoise能够输出多种长度动窗叠加，例如2天、5天、10天和30天的叠加。一旦用户确定了参考函数和一个(或若干)合理的动窗平均结果，在计算相对走时变化之前，两者都需要输出。同样要注意对参考函数的确定。如果将所有归档数据的平均作为参考，但是实际上它包含了强烈的速度变化，那么参考函数的质量将会大大降低(Duputeletal，2009；Mordretetal，2010)。强烈推荐在数据总时长度内调查比较稳定的时间范围。

只有全新的或更新过的那些天的数据才需要输出。如果前一天内任何互相关(CC)的任务已被标为“Done”(完成)，那么将计算叠加，并且会向数据库中插入新的走时变化(DTT)任务。注意，只对有新的或修订的互相关的日期，或已确定出新的参考函数的日期，才会计算后续的δt/t。如果定义了动窗的参考函数，那么每次都要计算参考函数的叠加，并且也要完全重新计算走时变化。对于一固定的参考函数，仅需要计算前一天的δt/t；这使得整个处理过程更快。

正如Hadziioannou等(2009)所述，地震波速变化测量的准确性依赖于时间域互相关函数的稳定重构。这样，作为第一步，必须确保有足够长的噪声时间序列进行互相关，目的是使每对台站的单个和参考互相关函数的相干性达到一定的阈值。一个明显的折中是时间序列越长，相干性就越高，而时间分辨率就越低。

2.3　相对走时变化

2.3.1时间延迟测量

可用于估计时间延迟的技术有两种：动窗互谱分析(MWCS；最早由Ratdomopurbo and Poupinet，1995提出)和被动干涉法(Sens-Schönfelder and Wegler，2006)。尽管有两项研究都强调了第二种技术的稳定性超过第一种(Duputeletal，2009；Hadziioannouetal，2009)，但是动窗互谱分析技术的优势在于频率域内的操作可以很明确地确定相关函数中相干信号的带宽(Clarkeetal，2011)。而且，Clarke等(2011)最近研制的方法已经对动窗互谱分析作了一些改进，尤其是对统计参数的估计，它可以评估变化的程度，在何种程度上对扰动的观测是可靠的。在MSNoise中，按照Clarke等(2011)的方法计算时间延迟，为保证程序包的一致性，用Python语言进行了重写。在文中的“检验”一节，我们比较了MSNoise和使用Daniel Clarke测量的FORTRAN程序(私人通讯，2011)获得的结果。

图3　UnderVolc计划中任意两个台站(YA.UV02和YA.UV05)不同时期的互相关函数(CCF)(上图)。水平虚线代表正负时间滞后的范围，这个范围确定了与参考函数(此处是所有可用数据的叠加)比较的两个区域。通过比较不同动窗叠加的互相关函数和参考函数来确定相关系数(下面4幅图)。富尔奈斯火山的喷发用红色垂直条带表示

将当前的互相关函数与参考函数进行了比较。两个时间序列被划分成若干存在重叠部分的窗。在经傅里叶变换转换到频率域之前，对每个片段做去平均和余弦尖灭(在两端进行1%的尖灭)处理。Fcur(ν)和Fref(ν)是傅里叶变换后埃尔米特对称的频谱的前一半，长度被补零至等于下一个2的整数幂的长度。互谱X(ν)被定义为：

(4)

式中，星号表示复共轭。之后X(ν)通过与一个汉宁窗卷积进行平滑。两个时间序列的相似程度可用频率域中两个能量密度间的互相干性来评估：

(5)

(6)

对每一处理窗口，两个信号间的时移是研究频带内采样加权线性回归后获得的斜率m。不像Poupinet等(1984)，这里所说的权重由Clarke等(2011)提出，包含了互谱的振幅和互相干性。误差通过使用权重(也就是相干性)和拟合差的平方进行估计：

(7)

(8)

这个过程的输出是一张表格，对于每个动窗都包含中心时间滞后、测量的延迟及其误差，以及该段的平均相干性。表格以文本格式[逗号分隔值(csv)]保存至硬盘上。对于任何给定的天，MSNoise计算M组(公式1)台站对的延迟时间和滞后时间，并把结果保存到文本文件中。

2.3.2波速变化计算

如果假设相对速度变化δv/v在空间上是均匀的，并且参考函数与当前相关函数之间存在相对时移δt/t，那么有(Ratdomopurbo and Poupinet，1995)：

(9)

以下各节中，我们将主要讨论δt/t，因为这是MSNoise保存下来的原始信息，并且必要时可以很容易地将其转换成δv/v。

上一步保存的动窗互谱分析延迟时间表可作为延迟矩阵来可视化(图4)，每天一张图显示了台站对之间微小的延迟差异是关于时间滞后的函数。MSNoise不仅能够计算每对台站的δt/t，也可以计算滤波后所有台站对加权平均的δt/t。当所有台站对被载入后，矩阵的底部会额外增加一行。这一行在MSNoise中叫做“ALL”(全部)行，包含了每一列的延迟时间在滤波后的加权平均值。对于每一列(滞后时间)，延迟时间必须满足一定的规则才能参与加权平均计算。在“检验”一节给出了对于UnderVolc火山数据中数据取舍规则的例子。

图4　UnderVolc计划期间2010年10月12日所有可用台站对的延迟、误差、相位相干性、数据选择和延迟时间变化矩阵的实例。矩阵的每一行相当于一对台站，而前4个矩阵的每一列是动窗互谱分析处理过程中动窗的中心时间滞后。每个矩阵都有自己的色标。矩阵底部的行是对所在列进行滤波后再加权平均的结果(详见正文)，它在MSNoise中被称为“ALL”行。“ALL”行仅在时间滞后-left_maxlag：-left_minlag和right_minlag：right_maxlag(虚线)之间才被计算。由选中的数据通过加权最小二乘回归(强制或没有强制穿过原点)确定的每个δt/t值被表示为画有误差棒的离散点。红色竖线是由“ALL”行确定的δt/t(不强制过原点：0.074%±0.009%，强制：0.074%±0.008%)，而绿线是所有δt/t的加权平均(没有强制：0.081%±0.065%，强制：0.074%±0.069%)。强制和没有强制两种方式获得的δt/t值非常接近，但是后者的误差比前者大一个数量级

图5　图4的实例矩阵中δt/t值以离散点的方式呈现。本图中，每一行或每一列都代表唯一的台站，因此每个单元格就代表一对台站。单元的颜色依赖于每对台站计算出的δt/t(斜率)。这里，对角线为空，是因为没有计算自相关。根据互易定理，每个台站对的互相关函数仅需计算一次

对于计算时移的最小和最大滞后时间的选择是十分重要的。由于沿较长路径传播的散射波会积累较长的时间延迟，所以在火山环境中开展的大部分研究使用了互相关函数的晚至部分(例如，Duputeletal，2009，5~20s或Mordretetal，2010，10~30s)。最近Anggono等(2012)的研究是个例外，他们在研究中使用了±10s到±1s的滞后时间。由于尾波部分(>10s)的波形相似度较差，他们既没使用动窗互谱分析方法，也没使用拉伸方法，而是比较了每天互相关函数和参考互相关函数的最大值。我们建议用户预先测试不同的窗口，并且避免在零时滞附近操作。而且，互相关函数仅在火山区，特别是较短滞后时间处(即由Duputeletal，2009描述的弹道波)对称。这或许说明了源分布不均匀，或者在某一特定位置存在次级噪声源的散射体(Paul，2005；Stehlyetal，2006)，亦或反映出有时钟问题(Sens-Schönfelder，2008)。Mordret等(2010)平均了互相关函数的因果和非因果部分，得到重构的波场有更好的相干性和稳定性。然而，同样也可以决定仅对重构互相关函数中因果的或非因果的部分进行单独处理。由于背景噪声源会随时间变化(其频谱、位置和振幅均有变化)，因此地震波速的变化会存在偏差。噪声源性质上的这些不同主要影响重构直达波。因而，为了估计这种偏差，可以初步只考虑测量互相关函数直达波或早至波的走时时移，从而测量地震波速的变化。

接下来对矩阵的M+1行通过满足要求的点(图4中数据选择矩阵中黑色部分)进行加权最小二乘回归(WLS)计算。加权最小二乘回归计算进行两次，第一次允许有非零常值截距(公式10a)，第二次要求强制回归经过原点([0，0])(公式10b)。计算每个参数的误差，在公式(10a)中，a和m对应的误差分别是ea和em；在公式(10b)中，m0对应的误差是em0。存在较大的非零常数a反映出仪器可能存在零点漂移。单天的结果可以画成台站对序号所对应的δt/t的散点图(图4右边两幅图)，或是表示为一个矩阵图(图5)。前者还展示出了由“ALL”行确定的斜率(红色竖线)和由所有台站对确定的斜率的加权平均值(绿色竖线)。在这个例子中，两个结果高度相似，但是，即使从直观上看，有些台站对也似乎表现异常。矩阵(图5)有助于识别出可以放在一起分析的台站类群。以UnderVolc项目的情况为例，我们可以定义两个台站类群，例如“火山口台站”和“大间距台站”，

δt=a+mt

(10a)

δt0=m0t

(10b)

式中，下标0用于强调回归直线强制通过[0，0]。

对于每个台站对和“ALL”行，这6个数值(斜率m和m0，常数a和它们的误差)都被存至硬盘，以便绘图程序或输出程序很容易地读取。

3检验

为了检验我们的处理流程的正确性，我们用相同的数据集进行了测试，该数据集在UnderVolc项目期间被成功用于识别2010年10月和12月两次喷发前富尔奈斯火山下的前兆速度变化。

3.1　数据和参数

数据集由布设在富尔奈斯火山区的21个宽频带三分量地震台站组成，台站间距从<1km到>10km不等(图6)。大部分台站于2009年11月初架设，并在2011年6月底撤除(图2)。事实上，大部分台址由富尔奈斯火山观测台网接管，并使用了新的台站代码。本次测试中，共处理了11 393个Z分量文件。为了计算互相关函数，对100Hz的波形进行了0.01至8.0Hz的带通滤波，然后抽样到20Hz(因子为5)。对于不同的台站对，仅计算了ZZ分量(未作自相关)。每天的数据都被分割成每30分钟一段，产出的互相关函数与按天叠加的互相关函数被一起存入硬盘。按2天、5天、10天和30天动窗进行叠加计算。动窗互谱分析的频带是0.2~0.85Hz，以10s动窗50%重叠计算。对于δt/t计算，都使用了因果和非因果部分，滞后时间选在±5s和±50s之间。其余的选择参数为：最小相干性是0.5，最大误差是0.1s，最大的dt是0.5s。表1和表2给出了在4v CPU的虚拟机上运行的时间和数据量的示例。

3.2　结果

应用相同的计算延迟时间的方法，分别利用Clarke等(2011)的FORTRAN程序和MSNoise计算了UnderVolc项目期间火山的δv/v变化，两者的结果非常相似(图7)。大部分时间两者的差异较小，明显低于0.02%，而火山喷发时要高于0.05%。总体趋势是存在的，并明显表现为每年速度有0.2%的微小变化。可以去除Brenguier，Shapiro等(2008)意义上的这一长期成分(LTV)，以突出更短时期的变化。深入地分析这两种程序，我们认为微小的差异可能来自余弦平滑函数的构造方式或是相位计算程序。我们选择使用了numpy，scipy(Jonesetal，2001；Oliphant，2006)和statsmodels(Seabold and Perktold，2010)程序包中久经考验的数值计算子程序，而没有重写Clarke的FORTRAN程序中的相关函数。

图6　UnderVolc计划的台站分布图。白色三角形是台站，各台站对间的路径以白线表示，并覆盖在富尔奈斯火山的数字高程模型(DEM)上。没有数字高程模型信息的台站分布图是使用脚本plot_station_map.py自动生成的

步骤运行时间安装10分钟配置5分钟扫描数据档案30分钟识别102626个新任务30分钟处理互相关任务(1个滤波器)21小时叠加获得210对台站的参考函数(REF)33分钟按4种动窗对所有台站对叠加6小时计算动窗互谱(MWCS)8小时计算δt/t1小时

表2　数据档案和MSNoise输出的数据量

2009年11月、2010年1月、2011年10月和12月的喷发之前δv/v都减小(图8中向上的部分)。在2009年10月喷发前似乎没有前兆信号。我们可以看到，2009年11月、2010年10月和12月喷发的前兆时间大约为15天。对于后两次喷发，在真正喷发前1天的δv/v增大很明显。这可以解释为深部速度的增大，也可解释为喷发前几小时内强烈的地震活动或是数分钟内出现的喷发振颤。最后两种情况导致了噪声源区的变化，并干扰了这个系统。值得注意的是，移动叠加函数在喷发前5天与参考函数的相关程度变弱(图3)。正如Liu等(2010)所指出的，改变叠加持续时间的长度可能会在δv/v曲线演化中突出或去除信号。图3中以2天、5天、10天或30天叠加的互相关函数与同样的参考函数比较时，均出现了相关程度变弱的情况，但表现不同。例如，2009年11月喷发前δv/v的前兆变化在不超过10天的叠加持续时间是可见的，却在30天叠加时完全被去除(图8)。短周期(等于或少于几天)的变化可在1天、2天和5天的叠加中看到(图8)，可能具有不同的来源，例如与“局部”原因有关，如气压、降雨或地震活动的变化，或者还有像噪声源位置变化这样的“区域”原因。虽然都有可能，但几天的短周期使“局部”原因比“区域”原因更为可能。

4结论和展望

我们给出了第一个使用背景噪声计算和监测相对波速变化的完整软件包。MSNoise是一个集成的解决方案，当执行预定任务时，可以自动扫描数据档案，确定哪些任务需要处理。整个软件包均由Python语言编写，是开源的，并可通过网页http：//www.msnoise.org(最后访问时间2014年3月)免费获得。其处理步骤与在UnderVolc计划期间(Brenguieretal，2012)成功识别富尔奈斯火山(法国留尼旺岛)下前兆速度变化的另一程序(Clarkeetal，2011)进行了比较测试。对于相同的数据集，两种处理程序的结果高度相似。MSNoise是可插入和扩展的；对于研究者来说，这实际上是一个十分有用的工具，因为他们只需关注实现他们自己的处理程序，而这得益于该系统稳健的架构。

图7　(a)UnderVolc计划期间富尔奈斯火山地区的相对地震波速变化。绿色曲线为使用D.Clarke的FORTRAN程序计算的结果，红色曲线为使用我们的程序计算的结果，红色垂直条带表示喷发。(b)红色曲线为两种结果之间的绝对差异，用百分比表示，蓝色曲线为项目中每天可用的台站对数

图8　使用MSNoise对5种不同动窗叠加(1，2，5，10和30天)计算的δv/v的结果去趋势后进行比较，且在δt/t估计程序中没有强制斜率估计穿过原点(公式10)。富尔奈斯火山的喷发表示为红色垂直条带。整个台网的δv/v，即“ALL”行(红色)，与所有单对台站的加权平均(绿色)匹配。UnderVolc计划期间，5次喷发中的4次都以δv/v的变化为前兆(见正文)

我们的方法基于下述假设，即重构的互相关函数尾波部分的每个波包在其传播路径上都经历了均匀一致的相对波速变化(δv/v为常数)。这是一个强假设，可能不是在任何情况下都有效。其结果是测量得到的相对地震波速变化可能与研究介质中“真实的”δv/v有很大不同。主要的局限性产生于尾波波包可能经过非均匀的速度变化区域这一事实。在这种情况下，测量的δv/v就是真实δv/v的平均结果；例如，鉴于δv/v的值在空间范围内是常数这个假设，人们可能在δv/v增大时错误地解释波速的变化δv。

最后，尽管我们给出的对MSNoise的有效检验是在火山环境下进行的，但其应用不应当局限于这种特殊情况，我们期望看到它在不同环境中得到应用，并期望得到来自其他研究小组的评论意见和贡献。例如，这些贡献可能是一种新的波形相关算法、新的δt/t计算技术、衰减计算的实现(如在Prietoetal，2011中)、延迟矩阵在定位δv/v变化中的运用，或者在不同时间滞后上的解相关计算(如在Laroseetal，2010；Obermannetal，2013中)。其他方面的贡献可能包括分析更高频数据(高于1kHz，如在Olivieretal，2012中)，或者将互相关函数输入到频率-时间分析(FTAN)(如Levshinetal，1989中描述的)，用于计算频散和随后的一维面波速度变化的反演。我们正计划实现定义多种参考函数以及“台站分类”的可能性；前者就像定义多种叠加方式一样，后者则类似于定义“ALL”行(见“波速变化计算”一节)。

参考文献

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.，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.，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.，P.Kowalski，T.Staudacher，V.Ferrazzini，F.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.，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.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.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

(中国地震台网中心安艳茹、张莹莹译；马延路、邹立晔校；吕春来复校)

安艳茹(1986—)，女，中国地震台网中心工程师，主要从事地震监测和地震背景噪声方面的研究，Tel：59959348，E-mail：anyanru8@126.com。

译 者 简 介

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