利用背景噪声研究地震波速度变化及其在长白山火山监测中的应用

刘国明

(长白山天池火山监测站，安图 133613)

0 引言

火山监测主要依靠火山地震活动探测、地表形变及流体地球化学等方法。应用这些方法，可以探测到岩浆升压和运移过程中的一些信息。地表形变学方法的局限性在于它们都是基于对地表应变及变形的解释，这就使得它们对深部物质的变化不灵敏。通过对火山地震的精确定位，可以掌握火山地震的时空分布特征，从而使我们了解到火山区地下介质属性的变化信息。然而，在没有发生岩浆活动的时间段，即在没有火山地震发生的情况下，使用火山地震学的方法就很难探测到火山区地下介质属性的变化信息。

近年来，利用地震背景噪声恢复格林函数的方法进行高分辨率的面波成像技术被广泛应用于空间结构探测。而在一些地球物理研究领域中，探测地震波速度的连续瞬时变化比探测地下波速结构更为重要，例如，研究与强震相关的波速变化(Brenguier et al.，2008a;Chen et al.，2010;Cheng et al.，2010;Froment et al.，2013)。对于火山监测而言，即使在地下波速结构未知的情况下，掌握地下介质的相对波速变化也是非常重要的。Sens-Schöenfelder等(2006)提出了应用地震背景噪声互相关技术反演火山区地下介质瞬时波速变化的方法，并探测到了印度尼西亚默拉皮火山地震波相对速度的季节性变化;Brenguier等(2008b)报道了留尼汪岛火山喷发前＜0.1% 的相对地震波速变化;Mordret等(2010)报道了新西兰佩鲁阿胡火山喷发2天前0.8%的相对地震波速度下降。应用地震背景噪声研究地震波速度相对变化的方法已经被广泛应用于活动火山监测研究中，但其完整的数据处理过程却未见详细报道。

长白山火山位于中国与朝鲜边境，作为中国最具潜在喷发危险的活火山，近年来很多学者对这座火山进行了研究。吴建平等(2005)和刘国明等(2006，2011)研究了长白山火山的地震活动性，并指出自1999年以来长白山火山活动经历了3个明显的阶段;上官志冠(1997)和上官志冠等(2006)研究了长白山火山区的He同位素比值变化并指出长白山火山正处于岩浆活跃状态;李克等(2009)研究了长白山火山区流动GPS和流动水准监测资料，发现在2002—2005年期间，长白山火山区垂直形变和水平形变都存在明显的异常。

目前长白山火山监测的主要方法有流动GPS、流动水准监测、定点应变和定点倾斜连续观测、流体地球化学定期观测和地震观测等。流动GPS和流动水准测量能精确地反映火山锥体的水平和垂直变化信息，但是每年只进行一期观测，时间分辨率太低，无法监控到岩浆的动态变化过程。定点应变和定点倾斜测量时间分辨率很高，但是目前的定点形变严重受当地雷雨的影响，数据的连续性非常不好。而且所有这些形变观测都是基于对地表变形的解释，无法很好地反映地下岩浆的微观动态过程。流体地球化学方法可以反映深部的岩浆活动信息，但是目前采样不连续，而且分析精度也不高，资料的利用价值不大。通过对火山地震的精确定位，可以了解岩浆活动的过程，但是用这种方法也无法掌握非震时段的岩浆活动情况。

虽然很多学者已经从不同角度应用不同的方法对长白山火山的活动状态进行了深入的研究，但长白山火山岩浆活动在时间尺度上的动态过程至今仍不被人们了解。我们通过对地震背景噪声进行互相关计算从而恢复格林函数的方法，获得了长白山火山区地下介质相对地震波速度变化过程，为研究长白山火山岩浆动态活动过程提供了一个新的方法。

1 单台数据预处理

单台数据预处理包括对每一个台站数据的预处理过程，其目的是去除诸如地震、畸变等可能对正常的背景噪声产生影响的信号，以突出背景噪声信号。通常，先把原始地震数据转换成以天为单位的SAC格式的波形数据，然后再进行数据处理。单台数据预处理的操作流程:去除仪器响应，去除均值，去除趋势变化和带通滤波，时间域归一化和谱白化处理。其中一些数据处理过程，比如时间域归一化和谱白化，对原始波形进行了非线性操作，因而这些数据处理过程的顺序通常是不能更改的。

2 互相关计算和数据叠加

互相关是通过识别时间延迟来判断2个波形相似程度的一种计算，所谓时间延迟是指一个信号相对另一个信号移动的尺度，通过这样的移动使之与另一个信号达到最大程度的相似，这个移动尺度就叫时间延迟。2个信号a和b的互相关函数Ca，b是时间延迟τ的函数，通常被定义为

这里的积分长度是整个记录的长度，u是时间信号函数的振幅。

互相关的作用在于它突出了地震波的走时。一个穿行于2个台站间的波场会在每个台站的地震记录中产生相似的信号，只是在记录到达的时间上会有偏移。因而台站对记录的互相关函数在某个时间延迟下会到达其峰值，这个峰值与台站对间波场的地震波走时一致。

通常，在时间域里对2个台站每天的数据进行互相关计算，然后对这些互相关函数进行叠加(stack)操作，以产生更长的时间序列，其目的是获得更稳定的互相关函数。互相关函数是包含正半部分和负半部分的双边的时间函数，因此出现在互相关函数中的所谓面波的经验格林函数(GFs)也包含了穿行于台站对2个方向间的射线路径信息。

为了提高信噪比(SNR)，采用适当的叠加长度是有必要的。图1显示了CBS－YNB台站对在不同叠加长度下获得的面波波形。图1a是经过0.01～0.2Hz带通滤波后获得的，图1b是经过0.02～0.1Hz带通滤波后获得的。通过比较，图1b中的面波波形更清晰，而图1a中的面波中有少量的短周期信号干扰。图1也表明选择适当的滤波频带对获取更稳定的面波信号是非常重要的。

图1 互相关函数中面波与叠加时间长度的关系Fig.1 Relationship between surface wave and stacking time.

3 测定相对地震波速度变化

在对所有的台站对数据进行了互相关计算和叠加计算后，就可以测定相对地震波速度变化了。台站对每天的互相关函数可以视为这2个台站间的面波格林函数。为了提高计算的信噪比，通常把经过连续的若干天叠加的互相关函数当作当天的经验格林函数。而把经过更长时间(通常为1a或2a)叠加的互相关函数当作参考格林函数。

假设台站对间地壳介质的相对地震波速度变化d v/v在空间上是均匀变化的，那么经验格林函数和参考格林函数间的相对走时偏移dτ/τ是不依赖于流逝时间τ的，并且d v/v=－dτ/t(Brenguier et al.，2008b)。当测定某一个以τ为中心的小窗口的经验格林函数与参考格林函数之间的走时偏移dτ时，如果dτ相对于τ是线性变化的，那么就可以通过测定走时偏移dτ的斜率来估计经验格林函数和参考格林函数之间的相对走时变化dτ/τ，而其相反数就是地震波速度的相对变化d v/v。

在具体的计算过程中，使用“滑动窗互相关法(MWCCA)”得到d v/v(Snieder et al.，2002;Pandolfi et al.，2006)。如图2所示，假设有2条波形曲线，fcur代表经验格林函数，fref代表参考格林函数，这2个波形的相关系数可以代表它们的相似程度。在要计算的波形中取一系列的小窗口，计算该窗口内的走时扰动(Δτi)，根据尾波干涉理论，2个波形互相关函数达到最大值的位置所对应的走时偏移即为走时扰动，相关系数R由式(1)求得:

式(1)中，ti为窗口中心点位置，窗口长度为2tw，ts为与最大相关系数对应的走时偏移值，即走时扰动Δτi。在一系列小窗口内重复上述计算，求得各个小窗口中心点位置的走时扰动Δτi，然后通过线性拟合得到整个经验格林函数相对参考格林函数的走时偏移Δτ/τ，其测量误差为

图2 相对走时偏移测量示意图Fig.2 Schematic diagram showing relative shift time measurement.

4 质量控制

评价互相关计算结果质量好坏的一个主要指标是信噪比。通常信噪比被定义成某信号窗口内最大振幅与噪声窗口内噪声振幅均方根之比。质量控制部分的目的是通过数据处理提高格林函数的信噪比，以得到稳定可信的测量结果。

4.1 经验格林函数数据叠加天数的选取

在进行经验格林函数计算时，选取较长的数据天数可以获得更稳定的计算结果，提高格林函数的信噪比。但是，如果数据叠加天数过长，也会平滑掉可能有用的短期变化信息。因此选取适当的数据叠加天数，对获得有意义的观测结果是非常重要的。

图3a和b显示了不同叠加时间计算得到的信噪比和相关系数。随着叠加时间的增长，测量结果的信噪比和相关系数得到了明显的提高。图3c显示了在不同叠加时间情况下获得的相对地震波速度变化值。11d数据叠加得到的测量结果杂乱无章，31d叠加和61d叠加获得的观测结果比较稳定，且二者差别不大，所以在最后的数据处理过程中我们采用31d数据叠加作为经验格林函数。

从图3c中还可以看到一个有趣的现象，在日本大地震前61d叠加得到的测量结果中有个明显的震前上升过程，地震后又转为急速下降。而在31d叠加的测量结果中，显示的是震前下降，震后上升，这主要是因为61d数据叠加显示的是长趋势变化，相对而言，31d数据叠加显示的是短期变化。2008年汶川地震时，波速变化也出现了类似变化。考虑到汶川地震和日本地震距离长白山火山区都非常远，还不能把这种变化认定为异常现象。

4.2 钟差问题

一些研究在进行走时扰动计算前还进行了所谓的钟差校正(Stehly et al.，2007;刘志坤等，2010)，该钟差表现为每天的经验格林函数与参考格林函数之间的整体走时偏移，即图2a中的拟合直线在通过走时0点时与0点值有较大偏差，而计算结果表明，本研究所用台站钟差较小(＜0.2s)，对波速测量结果基本没有影响，因此本研究未进行钟差校正处理。

5 长白山火山相对地震波速度变化

5.1 数据准备

本研究使用的数据是长白山火山区CBS、YNB、DHT和FST等4个宽频带地震台记录的2007—2012年的连续数字地震观测资料，其中CBS－YNB台站对计算了2001—2012年的地震波速度。对所有的台站对进行了互相关计算，最终获得了6个台站对的相对地震波速度变化结果。图4显示了长白山火山地理位置和本次研究所用到的地震台站分布情况。

5.2 数据处理过程

文中第2至第5部分讨论了应用背景噪声研究相对地震波速度变化的方法和基本流程。本节将总结数据处理流程并重点介绍本次研究中使用的具体处理方法和数据处理流程。基本的数据处理步骤如下:

(1)用RDSEED软件将原始数据(SEED格式)转换成SAC格式。

(2)将SAC数据文件合并成以天为单位的数据文件。

图3 叠加时间长度对计算结果的影响Fig.3 Influence of calculating results from length of stacking time.

图4 长白山火山地理位置和火山监测地震台网Fig.4 The geographical position of Changbaishan volcano and its seismic network.

(3)单台数据预处理及互相关计算。选定台站的数据以10为采样率进行重采样，以减少后续数据处理的计算量。去除仪器响应，去除均值和趋势变化，在0.02～0.1Hz频带范围内进行带通滤波，归一化处理，最后进行谱白化处理。完成单台数据预处理后，对台站对的数据进行互相关计算。

(4)对所有已经获得的互相关数据进行信噪比计算。现在计算信噪比的目的是保证互相关函数的质量，所有信噪比低于某一数值的互相关函数将被删除，具体取舍哪些数据，由研究者根据数据情况自行掌握。

(5)对选定台站对的互相关函数进行31d数据叠加，并产生经验格林函数。对地震活动相对平静的2a时间段的数据进行叠加，产生参考格林函数。本研究中，对2009—2010年的数据进行叠加，作为台站对的参考格林函数。

(6)应用MWCCA法，计算台站对之间地下介质相对地震波速度变化，同时，利用式(3)获得计算误差。

(7)计算经验格林函数的信噪比。

(8)计算结果选取。计算误差反映了相对地震波速度计算过程中的不确定性，根据最后计算结果的实际情况，将计算误差＞0.4% 的结果舍去。

5.3 讨论

图5显示了全部参与计算的台站对2007—2012年地震波速度的计算结果。图5d—f的计算精度明显高于图5a—c，而且图5d—f的长趋势变化曲线有很好的相关性，显示出很好的季节性变化规律。图5a—c中的相对地震波速度计算结果都使用了FST的数据，我们认为导致图5a—c计算结果不好的可能原因是FST的原始数据不理想，因而在接下来的讨论中，我们只讨论图 5d—f。

在研究时段内，发生了3个可能会对研究区地震波速度变化产生影响的地震。第1个是2008年5月12日汶川MW7.9地震，第2个是2010年10月9日发生的位于火山口20km的ML3.7和平林场地震，第3个是2011年3月11日发生的MW9.0日本大地震。因第1和第3个地震距离研究区太远(均＞1 300km)，第2个地震的震级太小，因而在图5 d—f的表现不明显。在2007—2012年时间段内，长白山火山区并未发生足以影响到地震波速度变化的地震事件，因而图5所获得的波速变化信息实际上是长白山火山区地震波速度变化的背景信息。

图5 2007—2012年相对地震波速度变化Fig.5 The relative seismic velocity variation from 2007 to 2012.

2001—2007年6月以前研究区内仅CBS台和YNB台有连续的数字化地震观测资料，其数据保存格式都是EVT格式。2007年6月以后台站进行完“十五”改造，各台站数据都改为由测震服务器存储，导出的数据也变成了SEED格式的数据。“十五”改造前后CBS台和YNB台的观测设备并未发生改变，仅仅是数据存储的格式发生了变化，并不影响前后数据的连续性。为了研究这段时间内长白山火山区地震波速度的变化情况，单独针对这些EVT格式的数据再开发了相关的代码，对这些数据进行了处理，从而获得了CBS－YNB台站对2001—2012年相对地震波速度变化情况(图6)。

从图6可以看出，2001年6月至2007年3月地震波速度的变化幅度达到2%，而其他时间段地震波速度变化幅度最大为1.5%，与此同时，长白山火山区出现岩浆扰动迹象，记录到了大量的火山地震，精密水准测量结果和流体地球化学的研究结论也都证实了在这个时间段内长白山火山经历了一次岩浆活动过程(刘国明等，2011)。

图6 CBS－YNB台站对2001—2012年相对地震波速度变化Fig.6 The relative seismic velocity variation of the CBS－YNB station pair from 2001 to 2012.

一个地区的应力或应变的变化会引起介质中孔隙的扩张或者压缩，可能在介质中产生新的孔隙，也可能会造成孔隙中流体的重新分布，从而导致该地区地震波速度发生变化。同震断层位移会导致区域应力场的变化，从而在应力拉张区域产生介质孔隙的扩张或者在介质中产生新的孔隙，而在应力压缩区域会产生介质孔隙的压缩(Nishimura et al.，2000)。因此，在应力拉张区域可能会产生地震波速度的下降，而在应力压缩区域可能会产生地震波速度的上升。

岩浆活动是一个复杂的动态过程。岩浆上升过程是一个增压的过程，体现在地震波速度变化上应该是一个上升的过程;岩浆下降过程是一个减压的过程，体现在地震波速变化上应该是一个相对下降的过程;同时，岩浆上升或下降的过程也会导致火山区局部应力场的重新调整，都会造成地震波速发生变化。我们推断，此期间地震波速度的异常大幅度变化，可能是岩浆活动造成火山区地下介质属性变化的一种外在表现，是岩浆运移动态过程的一种表现。目前在长白山火山区开展的精密水准测量和流体地球化学测量都是季节性的观测，基本上每年进行1期测量，然后将历年的观测结果进行对比分析，这些测量结果从不同角度证实了这次岩浆扰动事件，但是这些观测每年只进行一期，无法反映岩浆运移的动态过程。而我们此次工作获得了以天为单位的地震波速度变化，反映了岩浆活动的动态过程，从而弥补了其他观测方法的不足。

6 结论

应用背景噪声恢复地震波速度变化进行活动火山监测是近年来刚兴起的一个比较新的研究方法，其数据处理过程总体上分为4个步骤:1)单台数据预处理;2)互相关计算和数据叠加计算;3)相对地震波速度计算;4)质量控制，包括误差计算和计算结果的选取等。

计算了长白山火山区2007—2012年6个台站对的相对地震波速度变化。在0.02～0.1Hz频带范围内，通过背景噪声互相关方法获得了比较理想的面波信号。其中3个台站对的长趋势速度变化具有非常好的相关性，显示了明显的季节性变化规律，认为这是长白山火山区地震波速度变化的背景信息。

2001—2007年CBS－YNB台站对相对地震波速度变化幅度远高于正常背景值，与此同时，长白山火山发生了一次岩浆扰动事件，因而推断此期间地震波速度的异常变化可能与岩浆活动有关，岩浆的升压或降压作用导致了大量火山地震的发生，造成了火山锥体的异常升降变化，导致了温泉逸出气体浓度发生异常变化，同时也使火山区地下介质的物理性质发生了微观的变化，这可能就是2001—2007年CBS－YNB台站对相对地震波速度变化幅度异常的根本原因。

利用背景噪声互相关研究地震波速度变化是一个很有前景的研究领域，因为实现起来比较方便，在预测火山活动方面具有巨大的潜力，已经在活动火山监测方面获得了很多成功的应用。因此，利用背景噪声恢复地震波相对波速变化的研究将会在中国活动火山岩浆演化过程研究和火山活动趋势预测中获得越来越广泛的应用。