关于地震自动触发检测研究,吴治涛等(2010 年)曾利用STA/LTA算法进行过拾取微地震事件的研究,主要研究不同时窗下P波检测的敏感度;王继等(2006)曾用AIC的方法对初至震相进行过检测,主要分析了AIC方法和多台最小二乘互相关法的优势;不少科学家利用地震科学台阵(ChinaArray)数据进行研究(王未来等,2009;杜海林等,2009; 鲁来玉等,2009),还有些科学家利用流动台阵进行科学研究(刘启元等,1993;刘启元等,1997;李顺成等,2005;杜海林,2007);另外一些科学家研究了永久性小孔径地震台阵,并利用小孔径地震台阵进行研究(朱元清等,2002;郝春月等,2003;于海英等,2004;郝春月等,2006);但很少有人利用精确震相到时检测与到时修正法进行实际的地震数据库检测,并且是中国台阵项目记录的海量地震波形数据.中国台阵项目,即ChinaArray计划,是利用大孔径宽频带地震台阵,联合其他地球物理方法,对地壳和地幔结构特征进行研究,研究地震活动性与地球物理场演化规律. 在15~20年里分为7个阶段,完成中国国内的流体计划.针对此项目,我们已经记录了大量原始地震波形数据.对于海量的原始数据,用人工检索地震信号非常麻烦、耗时、费力.所以本文为了寻找更简洁更快速的方法从海量数据库中检测所有地震,实现初步数据处理,作者利用STA/LTA方法进行事件检测触发,用AIC初动到时修正方进行到时修正,结果成功地完成了对海量数据库进行地震检测的工作.
1 STA/LTA基本原理STA/LTA地震事件检测方法比较经典,具有原理简单,检测准确的优点.
该方法是取一个短时窗,长度为lsw,按下式计算任意时刻的短时平均值STA为(short term average)
再取一长时窗,长度为llw,按下式计算任意时刻的长时平均值LTA为(long term average)
假定llw 远大于lsw (通常10倍以上),LTA的变化要比STA的变化平稳很多,这样就区分了STA与LTA.
当发生一个地震事件,平稳的背景噪声中突然产生大振幅,则STA发生明显的变化,LTA也会变化,但不会那么快.跟事件发生前比,则STA/LTA的比值也会发生强烈变化,设定一个值,当这个变化超过此值,则认为检测到了该事件,这就是STA/LTA地震事件检测法的基本原理.
用STA/LTA的比值变化来检测地震事件还有一种解释,考虑到LTA变化不明显,可以理解为背景噪声,则STA/LTA也可以看作为信噪比,这样STA/LTA比值触发,也可以说是信噪比触发,即当信噪比达到或者超过一定水平后,并持续一段时间,可认为是地震事件触发.
我们要做的就是对事件触发的处理和事件发生过程中数据的处理,以及事件结束的判断.需要用到LTA延迟计算和LTA锁定来保证背景噪声的相对稳定.
2 AIC原理日本统计学家赤池弘次创立了AIC(Akaike information criterion)信息准则,该准则是衡量统计模型拟合好坏的一个标准,又称为赤池信息量准则.AIC是建立在熵的概念基础上,可以用来评估模型拟合数据的好坏程度和该模型的复杂程度.
多数情况下,AIC可表示为
式中:K为参数的数量,模型的误差服从独立正态分布是假设条件.设n为观察数,RSS为剩余平方和,则AIC为
要达到高度拟合,必须增加自由参数的数量,但是也要防止出现过度拟合的情况发生.所以AIC值最小的模型通常被选为应用模型.AIC准则的目的是寻找一个模型,该模型能够很好地解释数据,但是包含尽量少的自由参数.
对于地震波形,Maeda于1985年直接计算了地震序列的AIC函数,对地震波形x(i)(i=1,…,L)来说,定义AIC检测器为
为了能够自动检测事件,作者编写了事件自动触发程序,图 1显示了程序的流程图.该程序的第一步即进入初始化状态,在此状态,该程序要处理若干个数据包,数据包个数相当于长窗长度(单位为秒).初始化会赋予STA和LTA合理的初值.
![]() | 图 1 检测流程图 Fig. 1 Monitoring flow chart |
第二步,该程序进入检测事件状态,当STA/LTA>Td时(Td为设定的阈值),该程序进入验证状态,在此状态,如果STA/LTA>Tv(Tv为触发阈值)在整个窗内都成立,则表示了成功地验证,即事件触发成功.
设置好验证状态的参数可以抑制干扰信号引起的误触发,干扰信号的持续时间一般较短,导致短窗的长度明显小于验证窗长,如此,干扰信号造成的STA波动会在验证窗中过早结束,从而造成失败验证,减少了误触发.
事件触发成功后,该程序进入事件追踪状态.在此状态下,如果STA/LTA < Te(Te为结束阈值)在指定长度的结束窗内均成立,则表示事件结束,该程序重新进入检测事件状态.为防止持续时间较长的噪声源和大风天气造成的噪声等引起的误触发,设定了事件超长自动结束功能.如果事件启用了事件超长结束功能,则该程序回到初始化状态,赋予STA、LTA新的初始值,如此,持续时间长的干扰便不会再次引起误触发.如果干扰信号变小,则把LTA的值设小,从而提高了地震事件的检测能力.
短窗的长度应该大于输入信号的卓越周期.近震选择0.5~1.5 s,远震选择1.5~5 s比较合理.如果短窗长度短,则可以增加检测灵敏度,但太短则可能随信号周期的变化而影响检测效果.如果短窗长度长,则可以使检测更稳定,但过长则可能使检测灵敏度下降.卓越周期可根据地震记录统计得到,随软硬程度的不同,地基土有不同的卓越周期,可划分为四个等级:一级为稳定基岩,卓越周期为0.1~0.2 s,平均0.15 s.二级为一般土层,卓越周期为0.21~0.4 s,平均0.27 s.三级为松软土层,卓越周期在二级和四级之间.四级为异常松软的土层,卓越周期为0.3~0.7 s,平均0.5 s. 长窗长度应远远大于短窗长度,最少应在10倍以上,尤其以30倍以上为佳.另外,长窗长度应大于背景噪声的变化周期,否则,它作为背景噪声参考值的意义将变小.
图 2表示的是单台触发显示,表示的是2012年9月24日1时2分,ML 2.8级地震事件的波形.每张图中第一道波形为原始波形图,第二道为2~20 Hz滤波后的波形图,第三道为STA/LTA值图.图中第一道竖线是事件触发点,第二道竖线是检测到的事件结束点.第三道波形中的横线是STA/LTA的触发阀值.
![]() | 图 2 发生于2012年9月24日1时2分,震级为ML 2.8的地震事件检测图,第一道波形为原始波形图,第二道为2~20Hz滤波后的波形图,第三道为STA/LTA值图.图中第一道竖线是事件触发点,第二道竖线是检测到的事件结束点.第三道波形中的横线是STA/LTA的触发阀值. Fig. 2 The detecting diagram of the event occurred on Sep. 24, in 2012 whose magnitude is Ml 2.8. The first waveform is original waveform of the event, the second waveform is the one after filter between 2~20 Hz, the third waveform is the STA/LTA value. The first vertical line is the touch off point of the event, the second vertical line is the end point of the event in the upper two waveforms. The horizontal line in the third picture represent the touch off value for STA/LTA. |
科学台震中很多台站的数据有大量不同频率的背景噪声,使用特定频段的滤波处理,有的台站数据信噪比提升比较明显,有的则不是很有效.而同一台站在不同时间段的背景噪声频率也不一样,所以选择适合大多数台站的频率段(2~20 Hz)来滤波(图 3).也可以对同一台站在不同时间段进行不同频率的滤波,但处理起来比较繁琐.
![]() | 图 3 自动检测试例,每幅图中有两行波形,第一道为原始数据波形,第二道为滤波后数据波形,红色竖线为自动检测到的地震初动到时(a)原始波形信噪比较高的情况 (b) 原始波形含有仪器噪声(c) 原始波形信噪比低的情况 Fig. 3 Examples of the auto detecting, there are two waveforms in each picture, the first is the original waveform, the second one had been filtered, red vertical line represent the arrival time of the onset detected automatically (a) the situation of the original waveform with high SNR (b) the situation of equipment noises contained in the original waveform (c) the waveform with a low SNR |
如果触发检测参数设置偏小,则误触发率很高.为降低误触发率,测试中对不同种类地震设置了不同的参数,触发效果比较理想.具体参数如下:
(1)检查波形记录的每个时间段内是否有足够的触发台站,如果触发台站个数少于4个,则无法定位.设触发筛选窗的移动步长为30 s,窗内触发台站的个数为1,少于此值,一个事件的检测就结束了.
(2)检查各台触发的时间是否合理,任意两台P波的到时差不应该大于或小于台间距除以P波速度的值,过大或过小,则设为误触发.把两台的到时差与台间距除以P波速度的值比较,如果到时差在一定的误差范围内,则给两台设定权重为1.
(3)根据STA/LTA方法检测到初动时间后,截取合适的数据长度进行AIC计算,AIC算法则给出修正后的初动时间.
5 地震事件触发测试结果根据以上的原理、 步骤与参数设定,我们对科学台阵的数据进行了处理.首先选取18个地震事件进行了分析与测试(表 1).以2012年9月28日00时34分的ML 3.2级地震为例,科学台阵184个台站中有57个台明显记录到该事件.在利用上述方法对该事件进行自动触发测试中,有56个台检测到该事件,1台漏检.检测成功率为98%.在触发成功的56个台中,触发时间位置明显错误的有2个台,触发时间正确率为96%(图 4a).又如2012年9月24日1时2分的ML 2.8级地震,科学台阵184个台中有33个台明显记录到该事件.在对该事件的自动触发测试中,33个台检测到该事件,没有台漏检,检测成功率为100%.在触发成功的33个台中,触发时间位置明显错误的有1个台,触发时间正确率为97%(图 4b,图 4c).对于2012年10月8日15时5分发生的ML 3.7级地震,检测成功率为100%,触发时间正确率为98%.由于篇幅限制,在这里不再赘述.
|
|
表 1 用于触发测试的地震目录 Table 1 The catalog of earthquakes using in the test |
![]() | 图 4 使用AIC方法对STA/LTA方法检测事件初动的修正.每幅图上的标题是台站编码和事件的发生时间,上图表示STA/LTA方法检测事件初动,下图表示AIC方法对初动的修正, 红色竖线表示检出的初动 (a)-(d)分别为53233、53054、53155和53148台记录的2012年9月24日事件 Fig. 4 The correction of the onset by the using of AIC method. The title of every picture is the station code and the occurrence time of the event, the upper picture is the onset by the using of STA/LTA, the lower part is the correction of the onset by the using of AIC method, the red vertical line is the onset detected automatically (a)-(d) the event occurred on SEP 24 in 2012 recorded by 53233、53054、53155 and 53148 station |
选中测试的地震事件中包含10个2级以下的小震.小震具有震级低,振幅小,信噪比低等特点.信噪比低的事件很容易被漏检,为了解决该问题,我们把触发参数设置得比较灵敏,所以对于频率和持续时间与小震相似的信号在测试中都能被检测到.在小震检测过程中发现,部分地震台站的背景噪声干扰较多,造成很多误触发.过多的误触发会影响多台联合触发检测效果,产生多台联合检测的误触发和漏触发.为了解决该问题,我们统计分析了台阵记录的多个小震的功率谱,得出了小震的主要频段为2~4 Hz,从而把滤波器带宽从2~20 Hz改成2~4 Hz,滤除了尽量多的高频干扰信号,最后通过测试,得到了适合小震的STA/LTA窗长.对于震级为1.4~1.7的小震触发测试表明,检测成功率为98%,触发时间正确率为96%.
6 结论6.1 利用AIC方法比利用STA/LTA方法检测的初动到时准确率大大提高.触发测试表明,在本文处理的2级以上的地震事件中,检测成功率平均为99%,触发时间正确率平均为96.5%.而对于2级以下的小震触发测试表明,检测成功率为98%,触发时间正确率为96%.
6.2 该研究大大缩短了工作人员的时间与精力,使工作人员从繁琐的检验地震图的工作中解放出来,利用更多的时间进行科研工作.该研究大大提高了工作效率,从海量原始地震图中自动抽取包括小震的地震事件,对地震研究和地球物理研究的进一步工作,提供了有力的工具.
致 谢 感谢赵仲和研究员对本研究工作的指导!| [1] | Du H L. 2007. Analysis of the energy radiation sources of the 2004 Sumatra-Andaman Earthquake using time-domain array techniques[M. S.theis](in Chinese). Beijing:Institute of Geophysics, China Earthquake Administration. |
| [2] | Du H L, Xu L S, Chen Y T. 2009. Rupture process of the 2008 great Wenchuan earthquake from the analysis of the Alaska-array data[J]. Chinese J. Geophys. (in Chinese), 52(2):372-378. |
| [3] | Hao C Y, Zheng Z. 2006. Application of Correlation and Coherency Methods to Naqu Seismic Array Design[J]. Earthquake Research in China, 22(1):34-42. |
| [4] | Hao C Y, Zheng Z, Zhou G W. 2003. The calculation and analysis of correlation for Lanzhou seismic array site survey and the estimation of array response[J]. Acta Seismologica Sinica, 25(6):608-614. |
| [5] | Haykin S. 1996. Adaptive Filter Theory[M]. Upper Saddle River, New Jersey:Prentice Hall, 1-989. |
| [6] | Li S C, Liu Q Y, Chen J H, et al. 2005. Passive broadband seismic array observation across Tianshan[J]. Progress in Geophysics (in Chinese), 20(4):955-960. |
| [7] | Liu Q Y, Li S C, Shen Y, Chen J H. 1997. Broadband seismic array study of the crust and upper mantle velocity structure beneath yanhuai basin and its neighbouring region[J]. ACTA GEOPHYSICA SINICA, 40(6):763-773. |
| [8] | Liu Q Y, Li Z M, Shen L R, et al. 1993. GDS-1000 General portable digital seismic observation system[J]. ACTA GEOPHYSICA SINICA, 36(5):600-608. |
| [9] | Lu L Y, He Z Q, Ding Z F, et al. 2009. Investigation of ambient noise source in North China array. Chinese J. Geophys. (in Chinese), 52(10):2566-2572, DOI:10.3969/j.issn.0001-5733.2009.10.015. |
| [10] | Maeda N. 1985. A method for reading and checking phase times in autoprocessing system of seismic wave data[J]. Zisin, 38:365-379. |
| [11] | Wang J, Chen J H, Liu Q Y, et al. 2006. Automatic onset phase picking for portable seismic array observation[J]. ACTA SEISMOLOGICA SINICA, 28(1):42-51. |
| [12] | Wang W L, Wu J P, Fang L H. 2009. Crust and upper mantle S-wave velocity structure beneath Tanghai-Shangdu seismic array profile[J].Chinese J. Geophys. (in Chinese), 52(1):81-89. |
| [13] | Wu Z T, Li S X. 2010. Comparison of STA/LTA P-pickers for micro seismic monitoring[J]. Progress in Geophysics, 25(5):1577-1582. |
| [14] | Yu H Y, Zhu Y Q, Qin H W, Wang X P. 2004. Research on Result and the Calibration for the Location of Seismic Array in Shanghai[J]. Progress in Geophysics (in Chinese), 19(1):045-051. |
| [15] | Zhu Y Q, Tong Y X, Yu H Y, et al. 2002. Automatic recognition of seismic phases of regional earthquake[J].Northwestern Seismological Journal, 24(1):5-12. |
| [16] | 杜海林.2007. 2004年苏门答腊-安达曼大地震能量辐射源的时间域台阵技术分析[硕士论文].北京:中国地震局地球物理研究所. |
| [17] | 杜海林,许力生,陈运泰. 2009.利用阿拉斯加台阵资料分析2008年汶川大地震的破裂过程[J].地球物理学报,52(2):372-378. |
| [18] | 郝春月,郑重.2006.信号相关性方法在西藏那曲台阵设计中的应用[J].中国地震, 22(1):34-42. |
| [19] | 郝春月,郑重,周公威.2003.兰州台阵勘址测点相关值曲线的计算分析与初选台阵评估[J].地震学报,25(6):608-614. |
| [20] | 李顺成,刘启元,陈九辉,郭飚,赖院根,王继.2005.横跨天山的宽频带流动地震台阵观测[J].地球物理学进展,20(4):955-960. |
| [21] | 刘启元,李顺成,沈扬,陈九辉.1997.延怀盆地及其邻区地壳上地幔速度结构的宽频带地震台阵研究[J].地球物理学报,40(6):763-773. |
| [22] | 刘启元,李志明,沈立人,等.1993.GDS-1000宽频带通用流动数字地震观测系统[J].地球物理学报,36(5):600-608. |
| [23] | 鲁来玉,何正勤,丁志峰,等. 2009.华北科学探测台阵背景噪声特征分析[J].地球物理学报,52(10):2566-2572,doi:10.3969/j.issn.0001-5733.2009.10.015. |
| [24] | 王继,陈九辉,刘启元,等. 2006.流动地震台阵观测初至震相的自动检测[J].地震学报,28(1):42-51. |
| [25] | 王未来,吴建平,房立华. 2009.唐海-商都地震台阵剖面下方的地壳上地幔S波速度结构研究[J].地球物理学报,52(1):81-89. |
| [26] | 吴治涛,李仕雄. 2010. STA/LTA算法拾取微地震事件P波到时对比研究[J].地球物理学进展,25(1):577-1582. |
| [27] | 于海英,朱元清,秦浩文,王小平.2004.上海地震台阵数据处理及其在地震研究中的进展[J].地球物理学进展,19(1):045-051. |
| [28] | 朱元清,佟玉霞,于海英,等.2002.数字化台网的近震震相自动识别[J].西北地震学报,24(1):5-12. |
2015, Vol. 30





