STA/LTA算法在MEMS烈度计中的适用性分析
2021-08-05钟媛胡星星滕云田陈波何朝博沈晓宇
钟媛 胡星星 滕云田 陈波 何朝博 沈晓宇
摘要:为解决用于高密度布设的低成本MEMS烈度计集成软、硬件资源有限,且难以嵌入较为复杂算法的这一问题,基于Matlab的仿真计算,通过讨论在不同特征函数、时窗长度和短窗位置下STA/LTA值的变化趋势、拾取效果和运算时间,以选取能提高算法灵敏性、改善地震事件拾取效果和提高算法计算效率的参数,并将改进的STA/LTA算法应用于实际地震数据处理。结果表明:不同的特征函数对事件拾取率、拾取时间偏差和算法运算时间影响不同;长短时窗长度相差越大,STA/LTA值的变化越明显;时窗越长,算法运算时间越长;短窗置后可以增大STA/LTA值的变化幅度、减少算法计算量,改善算法拾取时间。改进的STA/LTA算法拾取效果更好,计算效率更高,占用内存资源更小,更适用于集成资源有限的MEMS烈度计。
关键词:MEMS烈度计;STA/LTA;拾取;运算时间
中图分类号:P315.61 文献标识码:A 文章编号:1000-0666(2021)02-0208-08
0 引言
地震烈度速报需要布设大量的地震仪器组建高密度监测网络,并要求仪器在断电后能依靠后备电池独立供电继续工作。相比于传统力平衡地震加速度计,MEMS加速度传感器具有体积小、质量轻、成本低、功耗低等特点(张海涛,阎贵平,2003;胡星星等,2013;王浩,丁炜,2013;李昌珑,2013),可大规模布设,适用于需要密集监测的地震预警和烈度速报网络。
王建军等(2009)、蔡莉等(2014),胡星星等(2015)都基于ARM9与Linux系统,采用MEMS加速度传感器研制了数字烈度计,仪器终端结合嵌入式系统和高性能单片机,集成网络传输功能,内嵌地震数据处理算法,實现了地震事件自动识别、地震动参数自动计算和通讯传输一体化,验证了地震数据处理算法嵌入设备的可行性。Peng等(2017,2019)基于Linux操作系统研制了一种可用于高动态范围的MEMS三轴加速度计,内置STA/LTA算法自动检测事件,用于地震预警。
MEMS烈度计以MEMS加速度传感器作为传感单元采集数据,利用嵌入式处理器和Linux操作系统组建控制处理系统,内嵌地震数据处理算法实时处理地震数据。为了降低仪器功耗和网络运行成本,MEMS烈度计应内嵌地震事件识别算法,以保证在识别到地震事件时才进行地震动参数计算和数据传输。
常见的地震事件识别算法方法有能量分析法(Allen,1978;Massan et al,2006)、自回归方法(Akaike,1973;Takanami,Kitagawa,1993;Sleeman,Eck,1999;赵大鹏等,2012)、小波变换法(刘希强等,1998;王喜珍,2004;滕云田等,2006;Xu et al,2019)、偏振分析法(Bai,2000;林建民等,2012)等。其中,长短时窗比方法(STA/LTA)、自回归信息准则法(AR-AIC)被广泛应用于微地震信号自动拾取(吴治涛,李仕雄,2010;刘劲松等,2013;段建华等,2015;Chen et al,2020)。为了综合利用信号与噪声的多个差异特征,马强等(2013)、蒋策等(2018)结合STA/LTA算法和AIC(池赤准则)算法相继提出了多步骤自动拾取方法,改善了算法拾取精度,可以实时处理地震台网传输的地震记录。
STA/LTA算法作为一种简单的能量分析法,计算量少、参数简单、计算效率高、适应性强,无需占用太多资源(刘劲松等,2013;刘晗,张建中,2014),适用于集成资源有限的MEMS烈度计。为了研究STA/LTA算法在MEMS烈度计中的适用性,本文基于MEMS烈度计采集的地震加速度数据对STA/LTA算法进行应用分析,主要从特征函数、时窗长度和短窗位置3个方面,对比了不同参数下STA/LTA值的变化趋势以及STA/LTA算法的拾取效果和运算时间。
1 STA/LTA算法的基本原理
STA/LTA属于能量分析法,是由Stevenson(1976)提出并应用于地震波初至到时的判别,其基本原理为:用STA(Short Time Average)和LTA(Long Time Average)的比值来反应地震序列幅值的变化。LTA反映了采集到的信号的背景噪声的变化趋势,STA反映了地震信号振幅(能量)的变化趋势(牟培杰等,2012)。当地震波初至时,STA早于LTA增大,以至于其比值有一个明显的增加。当某时刻的比值大于预先设定的阈值时,认为该时刻为地震事件的初至时间。传统STA/LTA算法的计算方法如下(Jubran,David,2012;李昌珑,2013;Kwon et al,2018):
STA(i)=1Nsta∑Nsta-1j=0CF(i-j),i=Nsta,Nsta+1,…,n(1)
LTA(i)=1Nlta∑Nlta-1j=0CF(i-j),i=Nlta,Nlta+1,…,n(2)
R(i)=STA(i)LTA(i)i=1,2,3,…,n(3)
式中:CF(i)为地震信号的特征函数;STA(i)、 LTA(i)分别为短时窗和长时窗的平均值;R(i)为长短时窗比值;Nsta、Nlta分别为长短时窗的长度;n为信号序列总长度;i为采样点。
2 STA/LTA算法的参数分析
从STA/LTA算法基本原理可以看出,该算法主要通过识别地震波初至时引起的长短时窗平均值比值的突增来拾取地震事件到时。地震事件初至时刻,STA/LTA值的变化幅度越大,算法的敏感性越强。算法的运算时间越短,资源占用越少。本文通过选用STA/LTA算法不同的特征函数,改变算法的时窗长度以及短窗位置来分析STA/LTA值的变化趋势、STA/LTA算法的运算时间和拾取效果。所有的实验仿真均基于win10系统的MATLAB(2018版)平台进行,计算机CPU为八核心八线程的酷睿i7-9700系列。
2.1 特征函数分析
特征函数(Allen,1978)是包含地震事件信息的重要参数,在地震波到达时应该能放大频率和幅值的变化,并能灵敏地反映这种变化(赵岑,2013)。选取合适的特征函数有利于提高算法灵敏度,改善地震事件到时的拾取精度。李昌珑(2013)用地震信号幅值绝对值作为特征函数,即CF1(i);马强等(2013)用差分放大后的地震信号作为特征函数,即CF2(i);蒋策等(2018)用地震信号幅值平方值作为特征函数,即CF3(i);常见的特征函数还有CF4(i)和CF5(i)等(余建华等,2011)。具体公式如下:
CF1(i)=X(i)(4)
CF2(i)=X(i)2+[X(i)-X(i-1)]2(5)
CF3(i)=X(i)2(6)
CF4(i)=X(i)2-[X(i-1)X(i+1)](7)
CF5=X(i)-X(i-1)(8)
式中:X(i)為地震信号;i为采样点。
设置短时窗长度Nsta为20个采样点(0.1 s),长时窗长度Nlta为100个采样点(0.5 s),计算不同特征函数下的STA/LTA值。本文基于MATLAB设置采样频率为200 Hz、时长为8 s的组合正弦波来模拟地震信号,组合正弦波的频率、幅值各不相同,在4 s处发生模拟地震事件。为了更好地分析,模拟地震信号只截取3.5~4.5 s。
如图1a-1所示,模拟地震信号幅值在4 s处由1变到5,频率为5 Hz。由于特征函数CF2(i)、 CF3(i)、 CF4(i)包含了信号幅值的平方项,其计算的STA/LTA值变化幅度大于特征函数CF1(i)和CF5(i)。设置模拟地震信号幅值在4 s处由1变到5,频率在4 s处由5 Hz变到10 Hz(图1a-2),再将其频率在4 s处由10 Hz变到5 Hz(图1a-3)。从图1b中可以看出,同时改变模拟地震信号的幅值和频率,无论频率改变从大到小还是从小到大,选择CF2(i)、 CF4(i)作为特征函数时,STA/LTA值变化幅度较大。
为了进一步分析特征函数对STA/LTA算法运算时间和拾取效果的影响,本文用不同的特征函数对2019年6月四川长宁地震后采集的10条地震数据进行事件拾取,触发阈值固定设置为10。特征函数CF1(i)~CF4(i)计算结果详见表1,其中偏差表示算法自动拾取到时与人工拾取到时之间的时间差,平均运算时间表示特征函数拾取不同地震数据时的STA/LTA算法的计算时间。当特征函数是CF5(i)时,STA/LTA算法平均运算时间为5.152 2 s,所有地震事件都因为阈值设置过大无法识别,事件拾取率为0。
根据表1以及特征函数CF5(i)的计算情况可得表2。从表2可知,保持STA/LTA算法的触发阈值不变,对相同地震事件进行拾取时,不同特征函数对STA/LTA算法的影响不同。从平均运算时间来看,特征函数CF3(i)最短, CF1(i)和CF2(i)次之, CF4(i)和 CF5(i)最长。从事件拾取率来看,CF2(i)和CF4(i)能100%识别地震事件,CF1(i)和CF3(i)能识别部分地震事件, CF5(i)由于阈值设置过大无法识别地震事件。从拾取时间偏差来看,CF1(i)、 CF2(i)、 CF3(i)拾取效果相当, 且远远优于CF4(i)、 CF5(i)。综合分析得出,5种特征函数中,CF2(i)地震事件拾取率高、运算时间更短、效率更高,更适合STA/LTA算法的嵌入式应用。
2.2 时窗长度分析
STA/LTA算法的关键在于利用长、短时窗内的平均值来反映能量的变化趋势,从而捕捉地震事件的特征。短时窗平均值STA主要反映了地震事件有效信号能量的瞬时变化,时窗短、变化快。长时窗平均值LTA主要反应了背景噪声信号能量的大小,变化相对平稳。分析时窗长度对STA/LTA算法的作用效果,就是分析不同时窗长度下STA/LTA值的变化趋势。以2019年1月30日在四川洪溪村采用MEMS烈度计采集的地震数据(图2a)作为样本,选取CF2(i)为特征函数,分析相同时刻、不同时窗长度下STA/LTA值的变化趋势及运算时间,如图2和表3所示,数据采样率为200 Hz,时长为1 h。
固定短时窗长度为40(0.2 s),分别取长时窗长度为400(2 s)、2 000(10 s)、4 000(20 s),在标注时刻为3 221.735 s时、Nlta=400时计算得到的STA/LTA最大值为9.87,运算时间为1.054 2 s;同时刻,Nlta=2 000和Nlta=4 000对应的STA/LTA值分别为46.88(增大了4.75倍)和87.78(增大了8.89倍),运算时间分别为4.600 3 s和8.893 0 s。由上可见,短时窗固定不变时,长时窗长度的增大,将导致STA/LTA值增大,STA/LTA算法运算时间增加。设置长时窗长度为2 000(10 s)保持不变,分别取短时窗长度为40(0.2 s)、100(0.5 s)、200(1 s),在标注时间为3 221.73 s时、Nsta=40时计算得到的STA/LTA最大值为46.89,运算时间为4.626 9 s;同时刻,Nsta=100对应的STA/LTA值为18.81,缩小了2.49倍,运算时间为4.757 6 s;Nsta=200对应的STA/LTA值为9.44,缩小了4.97倍,运算时间为5.017 6 s。由上可见,长时窗固定不变时,短时窗长度的增大,将导致STA/LTA值减小,STA/LTA算法运算时间增加。
综上所述,长、短时窗的长度对STA/LTA值的变化趋势有很大的影响。同一时刻,当长时窗固定不变时,短时窗越长,STA相对于LTA的变化越小,STA/LTA值越小;当短时窗固定不变时,随着长时窗长度的增大,STA相对于LTA的变化越大,STA/LTA值逐渐增大。
当地震到来时,STA/LTA值变化幅度越大,算法的敏感性越强。时窗选取过长,计算量增多,运算时间增长,资源占用增加;时窗选取过短,短期噪声干扰明显,仪器的抗干扰能力减弱。因此,选取合适的时窗长度尤为重要。
利用MEMS烈度计进行地震监测时,由于采用嵌入式控制器与存储器,资源有限,所用STA/LTA算法要尽可能降低资源占用。实验表明,选择不同的时窗长度,STA/LTA值变化范围为9.44~87.78,运算时间变化范围为1.054 2~8.893 0 s。因此笔者取短时窗长度为100,长时窗长度为2 000,以保证STA/LTA值的变化幅度较大、运算时间较短、资源占用较少。
2.3 短时窗位置分析
传统STA/LTA算法是将短时窗取于长时窗内,即长、短时窗内拥有一部分相同的信号序列,这在一定程度上弱化了短时窗内地震信号的能量特征。有研究者提出延迟STA/LTA算法,即将STA和LTA窗口分离,使STA和LTA窗口之间有一段时间间隔(Ruud,Huserbye,1992;段建华等,2015)。因此,本文将短时窗取于长时窗之后,比较短窗置后的STA/LTA算法与传统STA/LTA算法对STA/LTA值变化趋势的影响。改进的计算公式如下:
STA(i)=1Nsta∑Nsta-1j=0CF(i-j),i=Nsta+Nlta,Nsta+Nlta+1,…,n(9)
LTA(i)=1Nlta∑Nlta-1j=0CF(i-j),i=Nlta,Nlta+1,…,n(10)
R(i)=STA(i)LTA(i-Nsta),i=Nsta+Nlta,Nsta+Nlta+1,…,n(11)
利用上述公式,本文仍采用2019年1月30日在四川洪溪村用MEMS烈度计采集的地震数据作为样本进行分析。选择CF2(i)为特征函数,短时窗长度为100(0.5 s),长时窗长度为2 000(10 s),分析比较了相同地震事件下,两种不同短窗位置下STA/LTA值的变化趋势。从图3b可以看出,在3 222.04 s时,短窗取于长窗之内的STA/LTA最大值为19.18,STA/LTA值运算时间为4.674 6 s。同一时刻,短窗取于长窗之后的STA/LTA值达到452.84,增大了23.61倍,STA/LTA值运算时间为4.626 7s。
将短窗取于长窗之后,STA/LTA值的变化幅度增大,算法运算时间减少,运算效率得到提高。因此可以通过短窗置后来改进STA/LTA算法,使其更适用于MEMS烈度计。
3 算法应用
本文选用MEMS烈度计采集的多组地震数据进行分析,数据采样率为200 Hz,时长为1 h,选择CF2(i)为特征函数,短时窗长度为100(0.5 s),长时窗长度为2 000(10 s)。进一步分析比较短时窗取于长时窗之后的改進的STA/LTA算法与传统STA/LTA算法的拾取效果,结果见表4,其中偏差表示算法拾取到时与人工拾取到时之间的时间差。表4显示,改进的STA/LTA算法拾取时间最大偏差为0.92 s,传统STA/LTA算法拾取时间的最大偏差为1.94 s。对每一组地震数据,改进的STA/LTA算法平均运算时间为4.765 5 s,传统STA/LTA算法平均运算时间为5.111 2 s。由上可见,当选取相同的特征函数和时窗长度,设置相同的触发阈值时,改进的STA/LTA算法拾取效果更好,且该算法提高了计算效率,节约了运算空间,更适用于资源有限的MEMS烈度计的嵌入式应用。
为了进一步检验改进的STA/LTA算法对地震序列的处理能力,本文选取汶川8.0级地震的23个余震(表5)进行实验计算。按照主震-余震型地震序列,本文取其中110条不同震级、不同震中距的余震记录对改进STA/LTA算法与传统STA/LTA算法的拾取效果进行对比分析。选取特征函数CF2(i),短时窗长度为100(0.5 s)、长时窗长度为2 000(10 s),触发阈值设置为4。经过对比发现,传统STA/LTA算法漏识别7条记录,错误识别4条记录,正确识别99条记录,事件识别率为90%;改进STA/LTA算法,漏识别5条记录,错误识别4条记录,正确识别101条记录,识别率为91.8%。相比传统STA/LTA算法,改进的STA/LTA算法拾取到时平均超前0.14 s且运算时间短。这表明短窗位置的改变在一定程度上提高了算法的地震事件识别率,减弱了算法拾取到时的滞后度,缩短了算法的运算时间,减少了仪器资源占用。漏识别率较高反映出地震序列时间间隔短,长窗设置过长,实际检测中可适当减小长窗长度来改善。
4 结论
本文分析了STA/LTA算法在不同特征函数、不同时窗长度和不同短窗位置下的STA/LTA值的变化幅度、算法运算时间以及地震事件识别效果。选择最适合的MEMS烈度计参数,对STA/LTA算法进行改进。改进的STA/LTA算法选择CF2(i)为特征函数,短窗长度为100(0.5 s),长窗长度为 2 000(10 s),短窗置于长窗之后,其中长窗长度还应根据地震序列密集程度进行适当的缩短调整以提高识别率,减少漏检情况。将改进的STA/LTA算法应用于实际地震数据检测,地震到来时,与传统STA/LTA算法相比,改进的STA/LTA算法的STA/LTA值变化幅度更大、运算时间更短,选择合适的时窗长度也使得仪器资源占用更少。短窗置后有效减少了计算量、提高了计算效率、减弱了算法拾取时间的滞后度、改善了事件拾取效果,更适用于资源有限的MEMS烈度计。
感谢中国地震局工程力学研究所为本研究提供数据支持。
参考文献:
蔡莉,占伟伟,王秀,等.2014.基于微机电加速度计的地震烈度计[J].地震地磁观测与研究,35(1/2):240-246.
段建华,程建远,王云宏,等.2015.基于STA/LTA方法的微地震事件自动识别技术[J].煤田地质与勘探,43(1):76-80+85.
胡星星,高孟潭,滕云田,等.2015.地震动参数速报仪的研制[J].地球物理学报,58(3):1024-1034.
胡星星,滕云田,谢凡,等.2013.基于MEMS传感器的速度型地震计技术研究[J].地球物理学进展,28(1):515-522.
蒋策,吴建平,房立华.2018.地震检测与震相自动拾取研究[J].地震学报,40(1):45-57.
李昌珑.2013.基于MEMS加速度传感器的地震烈度监测技术研究与实现[D].北京:中国地震局地震研究所.
林建民,杨微,陈蒙,等.2012.偏振分析在地震信号检测中的应用[J].中国地震,28(2):133-143.
刘晗,张建中.2014.微震信号自动检测的STA/LTA算法及其改进分析[J].地球物理学进展,29(4):1708-1714.
刘劲松,王赟,姚振兴.2013.微地震信号到时自动拾取方法[J].地球物理学报,56(5):1660-1666.
刘希强,周蕙兰,郑治真,等.1998.基于小波变换的信号突变检测、滤波和瞬态谱研究[J].国际地震动态,(1):1-9.
马强,金星,李山有,等.2013.用于地震预警的P波震相到时自动拾取[J].地球物理学报,56(7):2313-2321.
牟培杰,王润秋,李彦鹏,等.2012.微地震事件的自动检测研究[C]//中国地球物理学会第二十八届年会论文集.北京:中国地球物理学会.
滕云田,王喜珍,王晓美,等.2006.用B-样条双正交小波拾取P波到时[J].地震学报,28(3):329-333.
王浩,丁炜.2013.MEMS加速度计与传统地震加速度计的比较研究[J].大地测量与地球动力学,33(S2):93-95.
王建军,吴荣辉,何加勇.2009.基于IPv6和无线网络的地震烈度计开发[J].现代电子技术,32(1):23-25,29.
王喜珍.2004.小波变换在地震数据压缩和震相到时拾取中的应用研究[D].北京:中国地震局地球物理研究所.
吴治涛,李仕雄.2010.STA/LTA算法拾取微地震事件P波到时对比研究[J].地球物理学进展,25(5):1577-1582.
余建华,李丹丹,韩国栋.2011.特征函数响应特性分析及STA/LTA方法的改进[J].国外电子测量技术,30(7):17-19,23,27.
张海涛,阎贵平.2003.MEMS加速度传感器的原理及分析[J].电子工艺技术,(6):260-262,265.
赵岑.2013.P波预警中的震级预测和PGA、PGV估算研究[D].成都:西南交通大学.
赵大鹏,刘希强,李红,等.2012.峰度和AIC方法在區域地震事件和直达P波初动自动识别方面的应用[J].地震研究,35(2):220-225.
Akaike H.1973.Information theory as an extension of the maximum likelihood principle[C]// Proceedings in second international symposium on information theory.Budapest Akademiai Kiado,267-281.
Allen R E.1978.Automatic earthquake recognition and timing form single traces[J].Bulletin of the Seismological Research of America,68(5):1521-1532.
Bai C Y.2000.Automatic phase-detection and identification by full Use of a single three-component broadband seismogram[J].Bulletin of the Seismological Society of America,90(1):187-198.
Chen Y N,Li J,Wang Z J,et al.2020.Quick phase identification for dense seismic array with aid from long term phase records of co-located sparse permanent stations[J].Earthquake Research in China,34(3):328-342.
Jubran A,David E.2012.Adaptive microseismic event detection and automatic time picking[C].Geo Convention 2012:Vision,1-5.
Kwon J,Heo T,Kim J K,et al.2018.A new P-wave detector via moving empirical cumulative distribution function[C].Bulletin of the Seismological Society of America,108(4):2080-2089.
Massan,Ferrettig,Spallarossad,et al.2006.Improving automatic location procedure by waveform similarity analysis:An application in the South Western Alps(Italy)[J].Play Earth Planet Interiors,154:18-29.
Peng C,Chen Y,Chen Q,et al.2017.A new type of tri-axial accelerometers with high dynamic range MEMS for earthquake early warning[J].Computers & Geosciences,100(C):179-187.
Peng C,Jiang P,Chen Q,et al.2019.Performance evaluation of a dense MEMS-based seismic sensor array deployed in the sichuan-Yunnan border region for earthquake early warning[J].Micromachines,10(11):735.
Ruud B O,Huserbye E S.1992.A new three-component detector and automatic single-station bulletin production[J].Bulletin of the Seismological Research of America,82(1):221-237.
Sleeman R,Eck T V.1999.Robust automatic P-phase picking:an on-line implementation in the analysis of broadband seismogram recordings[J].Physics of the Earth and Planetary Interiors,113(1-4):265-275.
Stevenson R.1976.Microearthquakes at Flathead Lake,Montana:A study using automatic earthquake processing[J].Bulletin of the Seismological Research of America,66(5):61-79.
Takanami T,Kitagawa G.1993.Multivariate time-series model to estimate the arrival times of S-waves[J].Pergamon,19(2):295-301.
Xu Z,Wang T,Xu S H,et al.2019.Active source seismic identification and automatic picking of the P-wave first arrival using a convolutional neural network[J].Earthquake Research in China,33(2):288-304.
The Applicability Analysis of STA/LTA Algorithm forMEMS Seismic Intensity Instruments
ZHONG Yuan,HU Xingxing,TENG Yuntian,CHEN Bo,HE Zhaobo,SHEN Xiaoyu
(Institute of Geophysics,China Earthquake Administration,Beijing 100081,China)
Abstract
Though a number of forecasting techniques are used for the detection of seismic events,each of them has its own pros and cons.Since the MEMS seismic intensity instrument which features low cost and follows the principle of high-density deployment has limitation in integrating software and hardware resources,it is difficult to embed some more complex algorithms.Based on the simulation calculation of Matlab,this paper aims to select the parameters that can improve algorithmic sensitivity,picking-effect of seismic events and calculation efficiency according to different factors of STA/LTA algorithm including variation tendency of the STA/LTA value,picking-effect and computing-time in various situations such as different characteristic functions,different window lengths and different short-window locations.And the improvement of STA/LTA algorithm is applied to practical seismic-data processing.The final analysis suggests that diverse characteristic functions have corresponding effects on facets like event picking-rate,picking-time deviation and computing-time of algorithms;the greater difference between the lengths of long-time windows and short-time ones,the more obviously the values of STA/LTA tend to change;the computing-time of algorithm depends on the length of the window (positive correlation);putting short windows after long ones can intensify variation range of STA/LTA values,reduce calculation amount and improve the picking time.In conclusion,the improved STA/LTA algorithm can improve the picking-effect,calculation efficiency and use less memory resources,so it is more suitable for MEMS intensity instruments with limited integrated resources.
Keywords:MEMS seismic intensity instruments;STA/LTA;picking;computing-time
收稿日期:2021-01-19.
基金項目:国家重点研发计划(2019YFC1511001)和中国地震局地球物理研究所基本科研业务费专项(DQJB19A0133)联合资助.
第一作者简介:钟媛(1995-),硕士研究生,主要从事地震观测和数据处理研究工作.E-mail:zhongyuan29@126.com.