三分量地震记录的互相关分析

马腾飞

1) 中国北京100081中国地震局地球物理研究所 2) 中国北京100033中国财产再保险有限责任公司



三分量地震记录的互相关分析

1) 中国北京100081中国地震局地球物理研究所 2) 中国北京100033中国财产再保险有限责任公司

地震波场本质上是三维矢量波场， 单分量记录实际上是三维矢量震动在某一方向上的部分投影． 本文基于单分量地震记录互相关公式， 提出了一种新的适合三分量地震波形记录多元互相关运算的简易方法， 并以2008年汶川MS8.0地震余震序列波形为例， 对其进行了效果验证． 结果表明， 相对单分量互相关， 该方法可以得到更为合理的全局最优结果， 解决波形识别匹配工作中不同分量间的差异问题． 该方法还可以利用不同分量间的“同源”信息， 有效压制随机噪声， 并从理论上说明其合理性． 其原理及计算过程均较为简单， 整体运算量较小， 适用于目前“大数据”时代的地震数据处理．

三分量记录 多元互相关 模板匹配方法

引言

波形互相关技术是地震学中一种常用的技术手段， 目前已经在“重复地震”识别(Schaff， Richards， 2004， 2011; Lietal， 2007， 2011; Maetal， 2014)、 余震事件检测(Peng， Zhao， 2009; Wuetal， 2014)、 低频事件观测与识别(Obara， 2002; Shellyetal， 2007)、 地震精定位(Waldhauser， Ellsworth， 2000; Schaffetal， 2004; Schaff， Waldhauser， 2005)等领域得到了广泛应用． 迄今为止， 地震学家已经发展了许多先进的数字记录处理方法来计算两个波列间的互相关系数， 但这些方法绝大部分是针对单分量的时间序列记录． 鉴于地震记录的三维空间属性， 各个分量仅是质点运动在垂直或水平方向的投影， 因此基于各个分量的单分量互相关系数计算可能会丢失地震记录中某些空间相关信息， 不能反映地震波在传播过程中波形和震相变化的全貌， 只得到片面性的结果． 事实上， 只有三分量耦合的空间记录才能真实地反映实际地震波场所包含的全方位信息．

就地震事件波形而言， 由于其不同震相在不同分量上有较大的运动学差异， 因此在采用互相关识别计算时对所用分量均采用统一时间窗口就显得不大“合适”． 例如： 由于地震信号中噪声的存在， 当初至P波到达时， 垂直分量会产生较明显的变化， 但水平分量依旧处于信噪比很低的“噪声模式”； 在S波尾部， 垂向信噪比较低， 水平向却还有较为明显的振动． 在这些互相关运算的时间窗口内， 低信噪比“噪声”的存在无疑会对最终的运算结果产生较大影响， 使事件的识别与检测面临较大困难， 对于震级较小的微震事件来说更是如此． 针对这些问题， 通过采用对不同分量的地震记录选取不同相关运算窗口的方法(如垂直分量时间窗口为P波到达后4 s， 水平分量时间窗口为S波到达后4 s)， 便可在一定程度上提高检测识别的准确率(Peng， Zhao， 2009; Mengetal， 2012; Wuetal， 2014)． 但是， 这种硬性规定的不同分量时间窗口难免会“错杀”一部分不符合这种“标准模式”的地震事件， 从而影响其识别的完整性； 与此同时， 当所获取的各个分量之间的差异性较大时， 如何对所得结果作出合理的解读也是一大难题． 对于台阵记录我们可以通过各种技术手段叠加不同台站的信息来达到压制噪声、 提高信噪比的目的(Leonard， Kennett， 1999; Kennett， 2000)， 但对于单台站地震记录则无法开展． 因此， 如何充分发掘不同分量间的“同源”作用， 得到能全面反映地震波三维属性的相关信息也是本文将要探讨的内容．

1 基本方法与理论

1.1 单分量波形互相关原理

(1)

1.2 三分量多元综合互相关系数

三分量地震记录有3个独立分量(垂直分量V， 切向分量T， 径向分量R)， 不同波列组合后可以形成一个3×n的矩阵. 对于多维矢量矩阵而言， 空间夹角没有意义， 这种情况下则不能用上述向量相关的思路来解决多元的相关问题． 在实际应用中， 通常对不同分量两两相关后计算得出3个独立的相关系数， 该相关系数矩阵可以用来描述两矩阵间的相关关系． 但在具体工作中我们也会遇到诸如不同分量间最大相关位置不一致、 各分量间相关系数差别较大等问题， 这给我们带来了较大的挑战． 事实上， 由于各分量间的的振幅能量、 信噪比(signal noise ratio， 简写为SNR)水平均不相同， 上述问题的出现在匹配识别工作中并不罕见． 如果可以找到一种简单快捷的方法， 能够综合考虑各分量的振幅能量水平， 得到全局最优结果， 无疑对此类工作的开展具有重要意义．

一种可行的办法为矩阵向量化， 即将三分量记录投影展开到一条直线上， 以便我们能继续使用向量相关的计算公式来处理三分量问题． 需要注意的是， 所采用的变换方式必须使各分量之间满足等价互易性， 否则所得结果不唯一．

图1 三分量地震记录投影到直线上的示意图 Fig.1 The diagram of three-component recordings projected onto a line

(2)

其中，

(3)

(4)

(5)

需要注意的是， 式(5)中fi(t)和gi(t)均为归一化前记录到的原始数据， 其中包含各个分量的绝对振幅信息． 由此可以看出： 这种“拼接”处理的实质在于可以将各个分量间不同的振幅及相关信息置于同一参考系下， 从而得出考虑全局后的整体结果； 同时将各分量波形置于更大参考系下也可以有效压制振幅较小、 相位不相关的噪声部分， 增大综合信噪比， 从而提高识别精度． 从形式上看， 式(5)与单分量相关公式很接近， 且其原理技术相对简单， 形式也较为简洁， 适合大规模地震数据资料的处理计算．

2 应用实例

以成都台(CD2)记录到的2008年汶川MS8.0地震余震序列中的两个地震事件波形为例， 详细分析所得三分量整体相关公式在实际中的应用效果． 该地震事件对的震中距为27.385 km， 为典型的近台记录， 目录参数引自中国地震台网中心的《中国地震月报目录》， 两个地震事件均属于微震事件， 震级几乎相等， 波形数据引自中国地震局地球物理研究所“国家数字测震台网数据备份中心”(Zhengetal， 2010)．

由于成都台采用甚宽频地震计， 在较大范围内的频率响应曲线较为平缓， 因此我们采用1—10 Hz四阶巴特沃斯(Butterworth)带通滤波器对去除均值、 线性趋势后的原始波形进行处理， 这也与前人工作的参数选取相一致(Lietal， 2007， 2011; Maetal， 2014)． 以P波到时为互相关计算起点， 所选取的窗口时长为P波与S波走时差的4倍， 滑动时长为±2 s， 这样便可包含全部的尾波序列， 同时避免后续噪声混入影响计算结果(蒋长胜等， 2008)．

采用上述流程， 计算所得的三分量波形互相关系数分别为0.81205(E--W)， 0.86067(N--S)和0.76896(U--D)， 如图2所示． 可以看出， 各分量的波形只在某些不同时段上具有较高的相似性． 由于各分量间的计算结果差异较大， 对于两次地震的判定也成为一大难题．

图3给出了两次地震在成都台(CD2)的三分量整体互相关系数以及各分量间的互相关系数随时间的变化曲线． 取P波到时前1 s起算， 以窗口中心为计算点， 2 s为互相关计算滑动窗口长度， 4倍的S-P走时差为窗口长度， 自左向右， 每次移动1个数据点(即0.01 s)进行互相关运算． 为了避免地震波的谐波特性对此处小窗口计算结果产生较大的影响， 每次进行窗口计算时， 波形只能相对移动±0.01 s (1个数据点)． 由图3可以看出， 南北分量与东西分量在全程运算中均有较高的互相关值(≥0.8)， 而垂向分量在P波和S波到达时间之外振幅较小部分的互相关值有较大的波动， 拖累了整体的相关系数计算， 因此未达到给定的阈值(≥0.8)． 此外， 不同时段的总体相关系数给予不同分量的“权重”不同， 振幅越大的分量权重越高， 因此采用三分量相关可以部分压制不同分量内部低信噪比的“噪声”部分， 提高识别的效率和准确率． 在上述例子中， 运用三分量整体相关公式后可以得出两个地震事件的互相关系数为0.8157， 满足重复地震识别互相关系数≥0.8的阈值条件，可以视为一对“重复地震”事件．

图2 成都(CD2)台记录到的2009-07-08(红色)与2010-09-12(蓝色)两次地震滤波后的归一化三分量波形对比图及其互相关系数， 图中绿色竖线表示P波到达时间

图3 成都(CD2)台不同分量互相关系数随时间的变化曲线

图4 不同分量计算得到的互相关系数与地震事件三分量综合后和水平(或垂直)分量的归一化互相关系数

与上述例子类似， 本文详细统计了成都台记录到的一系列不同地震事件之间不同分量的互相关系数差异， 得到了采用本文方法后所得到的整体互相关系数与各分量互相关系数之间的关系， 如图4所示． 可以看出， 各分量的相关系数总体上与全局相关系数呈线性对应关系， 但仅看某一分量有时会出现较大偏离， 因此在这种条件下有必要根据三分量整体相关系数对事件进行合理判断．

从图4中也可以看出， 即便是对于某一给定的地震事件而言， 3个分量相关计算的结果也会呈现出一定的规律性， 即水平分量的相关系数显著高于垂直分量， 东西分量的相关系数明显高于南北分量， 因此在只能选择一个分量作相关计算时， 垂直分量具有更高的识别可信度． 李宇彤(2012)利用区域台网对海城—岫岩地区“重复地震”识别的研究也得出类似的结论， 这也说明了仅用垂直向的波形数据进行相关运算的合理性(Schaff， Richards， 2004， 2011)．

3 讨论与结论

地震波本身为矢量场， 本质上为不同特性、 不同类型的振动相互叠加干涉的结果， 而单分量记录实际上仅为三维矢量在某一方向上的部分投影， 因此常规的基于各分量的单分量互相关计算可能会丢失信号中部分与空间相关的信息， 不能反映地震波在传播过程中波形和震相变化的全貌， 从而导致结果不一致． 本文基于三分量记录之间的“同源”特性， 提出了一种可以计算三分量记录总体相关系数的简易方法， 能够尽可能地利用数据的内在信息， 压制随机噪声， 并从理论上说明了其合理性． 该方法的原理和计算过程均较为简单， 整体运算量也较小， 无论从经济上还是技术上都适用于未来“大数据”时代的海量资料处理， 值得在实际工作中推广应用．

由文中的实例可以看出， 采用三分量整体相关可以解决互相关系数在不同分量间的差异以及临界识别等问题， 所得结果也不是简单的三分量单独相关运算后的算术平均值， 而是在各分量间(inter-component)相关后综合叠加得到的结果． 这实际上是一种对信号的压噪重构， 可以达到增强有效信号、 压制干扰噪声的目的， 解决了目前“重复地震”以及类似事件识别工作中遇到的问题． 此外， 采用三分量整体相关还可以同步三分量地震波形记录， 避免片面追求各个分量的互相关系数单独最大而造成错误时移(这种现象可能是由地震波的谐波特性所导致)， 因此从理论上来讲， 也可能会存在整体相关系数比3个分量都小的极端情况， 但在实际中由于各种震相混叠、 介质不均性等情况的客观存在， 故难以出现上述情形．

应该看到， 这种综合相关算法似乎对在有效震相外信噪较低的波形记录部分的相关计算效果并不明显， 因此如何压制有效震相外的信号噪声以增加信噪比， 以及充分利用各分量中所包含的地震信息以提高微小地震事件的可探测性也是未来工作的一个可行方向．

吴忠良研究员为本研究进行了分析和指导， 与加州大学圣克鲁兹分校地震学实验室的Emily Brodsky教授、 Lian Xue博士、 Stephen Hernandez博士进行了有益讨论， 中国地震局地球物理研究所“国家数字测震台网数据备份中心” (doi:10.7914/SN/CB)为本研究提供了波形数据， 审稿专家提出了建设性的修改意见， 作者在此一并表示诚挚谢意．

蒋长胜， 吴忠良， 李宇彤. 2008. 首都圈地区“重复地震”及其在区域地震台网定位精度评价中的应用[J]. 地球物理学报， 51(3): 817--827.

Jiang C S， Wu Z L， Li Y T. 2008. Estimating the location accuracy of the Beijing Capital Digital Seismograph Network using repeating events[J].ChineseJournalofGeophysics， 51(3): 817--827 (in Chinese).

李宇彤. 2012. “重复地震”的若干地震学问题[D]. 北京: 中国地震局地球物理研究所: 36--40.

Li Y T. 2012.TheSeismologyof‘RepeatingEarthquakes’[D]. Beijing: Institute of Geophysics， China Earthquake Administration: 36--40 (in Chinese).

Båth M. 1974.SpectralAnalysisinGeophysics[M]. Amsterdam: Elsevier Scientific Publishing Company: 87--94.

Fredricksen H， Kessler I. 1977. Lexicographic compositions and deBruijn sequences[J].JCombTheory：SerA， 22(1): 17--30.

Kennett B L N. 2000. Stacking three-component seismograms[J].GeophysJInt， 141(1): 263--269.

Leonard M， Kennett B L N. 1999. Multi-component autoregressive techniques for the analysis of seismograms[J].PhysEarthPlanetInt， 113(1/2/3/4): 247--263.

Li L， Chen Q F， Cheng X， Niu F L. 2007. Spatial clustering and repeating of seismic events observed along the 1976 Tangshan fault， North China[J].GeophysResLett， 34(23): L23309. doi:10.1029/2007GL031594.

Li L， Chen Q F， Niu F L， Su J. 2011. Deep slip rates along the Longmen Shan fault zone estimated from repeating microearthquakes[J].JGeophysRes， 116(B9): B09310. doi:10.1029/2011JB008406.

Ma X J， Wu Z L， Jiang C S. 2014. ‘Repeating earthquakes’associated with the WFSD-1 drilling site[J].Tectonophy-sics， 619/620: 44--50.

Meng X F， Yu X， Peng Z G， Hong B. 2012. Detecting earthquakes around Salton Sea following the 2010MW7.2 El Mayor-Cucapah earthquake using GPU parallel computing[J].ProcCompSci， 9: 937--946.

Obara K. 2002. Nonvolcanic deep tremor associated with subduction in Southwest Japan[J].Science， 296(5573): 1679--1681. doi:10.1126/science.1070378.

Peng Z G， Zhao P. 2009. Migration of early aftershocks following the 2004 Parkfield earthquake[J].NatGeosci， 2(12): 877--881.

Schaff D P， Bokelmann G H R， Ellsworth W L， Zanzerkia E， Waldhauser F， Beroza G C. 2004. Optimizing correlation techniques for improved earthquake location[J].BullSeismolSocAm， 94(2): 705--721.

Schaff D P， Richards P G. 2004. Repeating seismic events in China[J].Science， 303(5661): 1176--1178.

Schaff D P， Waldhauser F. 2005. Waveform cross-correlation-based differential travel-time measurements at the Northern California Seismic Network[J].BullSeismolSocAm， 95(6): 2446--2461.

Schaff D P， Richards P G. 2011. On finding and using repeating seismic events in and near China[J].JGeophysRes， 116(B3): B03309. doi:10.1029/2010JB007895.

Shelly D R， Beroza G C， Ide S. 2007. Non-volcanic tremor and low-frequency earthquake swarms[J].Nature， 446(7133): 305--307.

Waldhauser F， Ellsworth W L. 2000. A double-difference earthquake location algorithm: Method and application to the northern Hayward fault， California[J].BullSeismolSocAm， 90(6): 1353--1368.

Wu C Q， Meng X F， Peng Z G， Ben-Zion Y. 2014. Lack of spatiotemporal localization of foreshocks before the 1999MW7.1 Düzce， Turkey， earthquake[J].BullSeismolSocAm， 104(1): 560--566.

Zheng X F， Yao Z X， Liang J H， Zheng J. 2010. The role played and opportunities provided by IGP DMC of China National Seismic Network in Wenchuan earthquake disaster relief and researches[J].BullSeismolSocAm， 100(5B): 2866--2872. doi:10.1785/0120090257.

Cross-correlation analysis of three-component seismic recordings

1)InstituteofGeophysics，ChinaEarthquakeAdministration，Beijing100081，China2)ChinaProperty&CasualtyReinsuranceCompanyLtd.,Beijing100033,China

Seismic wave is a three dimensional vector wavefield， the single component recordings are actually the projection of particle motion along certain directions. Based on the single component seismogram cross-correlation formula， this paper presents a novel simple solution which is suitable for the calculation of three-component seismogram cross-correlation， and the effectiveness of this new approach is verifiedviaa practical case from the aftershock sequence of 2008 WenchuanMS8.0 earthquake. Compared with single component seismogram cross-correlation， the new approach can obtain a global optimized result more reasonably and erase the discrepancy between different components in the work of template waveform matching. Also this new formula can take advantage of the congenerous between different components， suppress the ambient seismic noise effectively， and its rationality was demonstrated in theory. This new approach requires rather small computations as its simplicity in principles and procedures， which is suitable for the seismic data processing in the current era of “big data”.

three-component seismogram; multiple cross-correlation; template matching method

国家留学基金委和中国地震局“地震科技青年骨干人才培养项目”(20105007)联合资助.

2015-02-17收到初稿， 2015-05-07决定采用修改稿.

e-mail: matengfei@cpcr.com.cn

10.11939/jass.2016.01.009

P315.63

A

马腾飞. 2016. 三分量地震记录的互相关分析. 地震学报, 38(1): 96--102. doi:10.11939/jass.2016.01.009.

Ma T F. 2016. Cross-correlation analysis of three-component seismic recordings.ActaSeismologicaSinica, 38(1): 96--102. doi:10.11939/jass.2016.01.009.