基于改进EWT和LogitBoost集成分类器的地震事件分类识别算法

孟 娟,张家声,李亚南

(1.防灾科技学院 电子科学与控制工程学院,河北 三河 065201;2.中国地质科学院 地球物理地球化学勘查研究所,河北 廊坊 065000)

0 引言

根据产生机理的不同,地震事件一般可分为天然地震和非天然地震两大类,人工爆破就是最常见的非天然地震类型之一。随着地震监测能力的不断提升,以及矿藏勘探和开发等人类生产活动干扰的急剧增加,地震台站每天都会监测到大量天然与人工地震数据,如果不能及时、准确地对这些记录进行分类,不仅会造成地震目录污染,降低地震观测资料质量,也会影响地震监测预警的准确度。因此,快速、准确地从震源事件波形中识别出天然地震与人工爆破事件一直是地震学研究及相关应用领域的一个热点问题。

围绕天然地震和人工爆破事件的分类识别,国内外众多学者展开了广泛而深入的研究,通过提取地震信号波形、震相、频谱等方面的特征,基于统计模式识别[1-3]、机器学习[4-6]、深度学习[7-8]等方法进行事件类型判别,取得了一定成果。其中,统计模式识别严重依赖于特征提取,深度学习需要大量的训练数据,且对硬件要求高[9],而机器学习对数据和硬件的依赖较小,且识别准确率高[10],因此利用机器学习算法进行地震事件分类成为趋势。尤其是反向传播(Back Propagation,BP)神经网络[11]、支持向量机(Support Vector Machine,SVM)算法[6,12-13]被广泛应用于地震事件分类,成为主流分类识别方法。但BP神经网络容易陷入局部最优,且隐层数、隐层节点等参数无法预先确定,识别准确率受限;而SVM性能优劣取决于核函数,若缺乏有效手段选取合适核函数将影响分类准确率。

在地震事件性质分类过程中,如何提取有效的地震信号属性特征是分类准确与否的关键。国内外研究学者对地震事件提取了多种判别依据,如波形时域特征、频谱特征等。这些方法各有特色,也有较好效果,但由于地震信号本身的复杂性,使得单独采用某种特征进行识别的普适性不够,识别效果受限,因此研究人员开始综合提取时频域或变换域特征进行识别,取得了较好的识别效果。如王婷婷等[14]利用小波包对信号分解,使用不同分解系数中的P波和S波能量进行地震与爆破的识别;Beccar-Varela等[15]通过小波变换进行信号分解,基于分解后的时频能量特征进行地震事件识别;杨千里等[16]利用广义S变换获取地震与爆炸时频谱图二阶矩进行识别。以上方法证明提取不同变换域特征能有效进行地震识别。毕明霞等[17]利用经验模态分解(Empirical Mode Decomposition,EMD)对地震信号进行分解后，提取地震记录特征进行SVM识别,识别率高。但小波和小波包变换的小波基选择缺乏自适应性,广义S变换依赖于窗函数调节参数,EMD分解存在模态混叠、伪模态等问题,于是兼具EMD和小波变换优势的经验小波变换(Empirical Wavelet Transform,EWT)被提出[18]。EWT能自适应分解信号,非常适合处理非线性、非平稳的地震信号[19-20],但传统EWT需要根据信号频谱幅值进行频带划分,易受噪声影响,因此有必要根据地震信号的特点进行频谱分割优化。Wang等[21]基于尺度空间表示提取傅里叶谱的慢变分量,再根据尺度参数和中心频率进行频谱分割,但该方法中尺度参数依经验获取,这样会影响频谱分割的准确性。

为提高识别的准确率和稳定性,集成学习被应用于地震事件分类识别,将识别效果受限的基分类器按照一定规则组合为强分类器,以提高识别效果。赵刚等[22]提取波形时域、小波变换特征等,以BP作为弱分类器,通过AdaBoost方法组合成强分类器,得到较高的识别率。蔡杏辉等[23]提取波形复杂度、频谱比、自相关系数、波形复杂度和自相关系数比值为特征,用决策树作为基分类器,用Bagging集成方法集合成一个强分类器,准确率可达85%以上。Bagging只是对分类器进行简单整合,AdaBoost能对错误分类的训练数据进行处理,准确率更高,但AdaBoost易受样本中的噪声影响,从而产生过拟合。LogitBoost算法是在AdaBoost的基础上改进而来的,对样本噪声不太敏感,泛化能力更强。

为准确提取地震信号特征,提高地震事件分类识别准确率,本文提出基于改进EWT和LogitBoost集成分类器的天然地震和人工爆破事件分类识别算法。首先,对地震记录进行筛选,剔除噪声干扰严重的记录,并对筛选后的地震记录计算P波与S波最大振幅比;在此基础上,利用S变换时间频率分析精度高的优势对信号进行S变换,获取能反映信号频率与能量特征的S谱能量曲线,再基于S谱能量曲线对传统EWT进行改进,将地震信号自适应分解为若干按频率和能量分布的本征模函数(Intrinsic Mode Function,IMF)后去除噪声,然后提取各IMF的香农熵、对数能量熵,及去噪后重构记录的主频等特征;最后,基于LogitBoost机制的决策树集成分类器进行事件识别,以有效提高地震事件识别的准确率。

1 基于S谱能量曲线的改进EWT

1.1 S谱能量曲线

作为小波变换和短时傅里叶变换的发展,S变换采用带有频率变量的高斯窗函数截取信号,实现信号的局部分析。由于其窗函数可随频率自适应改变,具有多分辨率的特点,能较好地适应非平稳信号频率不规律变化的特点,在地震信号处理领域应用较广。对某地震信号h(t)进行S变换:

(1)

式中:f为频率;τ为平移因子,用于控制高斯窗在时间轴t上的位置;j为虚数单位。

对信号进行S变换得到信号时频分布S谱,可同时从时域和频域分析其特性和能量分布。对频点i定义S谱频点能量:

(2)

式中:fi为S变换后频点i的频率,S(fi,τ)为该频点对应的S变换值即S谱值。

获取所有频点的能量谱值LEi,能够得到依频率分布的S谱能量曲线。该S谱能量曲线能清晰描绘不同频率处信号的能量分布情况及不同频率成分的能量强弱。

1.2 基于S谱能量曲线的改进EWT

设地震信号s(t)的傅里叶频谱为F(w)w∈[0,π],S谱为S(f,τ),按式(2)计算S谱能量曲线后对地震信号进行改进EWT分解。

S谱能量曲线能真实体现原信号低频到高频部分各频率成分的能量大小分布,基于S谱能量曲线的改进EWT可将原地震信号准确分解为一系列依频率和能量分布的IMF。分解后的IMF瞬时频谱和幅值具有较高的分辨率,不仅能刻画原信号的时频特征,也能更好地描述信号的局部突变或渐变特征,从而更好地获取原信号特征。

2 数据预处理与特征提取

本文数据来源于中国地震数据共享中心。其中,人工爆破记录是2013—2019年发生于广西平果、武鸣、贵港、乐业等地的矿藏勘探、开采等人工爆破微震事件波形记录,震级为ML1.3～3.0,共973条;天然地震记录是发生于2019—2020年每年8—10月、震级ML1.3～3.0的天然微震波形记录,共867条。由于某些地震事件的台站波形记录中噪声干扰较多,因此需要对地震记录进行预处理,对被噪声淹没的波形进行筛选剔除,以更精确地提取地震信号特征,消除噪声影响。

2.1 波形筛选

对地震事件进行分类识别前需要对记录进行筛选并拾取P波到时。目前应用较广泛的初至拾取方法有短时平均值/长时平均值(Short Term Averaging/Long Term Averaging,STA/LTA)、赤池信息准则(Akaike Information Criterion,AIC)、无监督聚类算法等，本文利用快速简单的STA/LTA算法进行波形记录筛选,具体流程如下:

(1) 截取包含P 波和S波在内的全波段波形记录;

(2) 设定STA,LTA和阈值,取STA=0.2 s,LTA=1 s,阈值为1.5;

(3) 基于STA/LTA进行P波拾取,取初至时间点的前10 s至后100 s的记录波形。

按以上步骤进行波形筛选,得到天然地震记录648条,人工爆破记录575条。

2.2 改进EWT分解

对筛选后的记录进行计算，得到其S谱能量曲线,取能量曲线上的前4个极大值点,以其为频谱分割初始边界。根据ε邻域法重新确定新边界，进行改进EWT分解,最终将每条记录分解为5个IMF分量。由于随机噪声在地震记录中广泛存在,且一般以高频为主,因此将最后一个IMF视为高频噪声并丢弃,尽可能降低噪声影响。去噪后的每条地震记录对应4个IMF,这4个IMF主要为有效波,既保留原信号时频特征,也能刻画信号局部变化细节,有利于提取信号特征属性。

2.3 特征提取

为了消除各维数据之间的数量级差别,对所有输入数据进行归一化处理,将数据归一化到[-1,1]区间。信号P波与S波的最大振幅比是天然地震和人工爆破波形信号中的一个重要特征[23-24],经STA/LTA识别P波初至位置A后,按式(3)计算P波与S波最大振幅比ps:

ps=max(A-50:A+50)/max(A+50:A+1 900)

(3)

熵用于衡量信号的复杂程度,利用香农熵可以分析信号时频分布的聚集性,利用对数能量熵可反映信号的暂态或突变特性,因此香农熵和对数能量熵成为微震和爆破信号识别的重要参数[25]。对分解出的4个IMF,按式(4)分别提取每个IMF的香农熵Esh(i)和对数能量熵Elg(i)。

(4)

式中:si,j为EWT分解后分量i(即IMFi)的第j个值;M为每个分量的长度。

此外,由于爆破信号的频谱特征与天然地震明显不同,其能量较集中,且频率更低[26]。基于这4个IMF进行原信号重构,对重构后的信号进行希尔伯特变换。希尔伯特变换是一种适合非线性、非平稳信号的时频分析方法,常用于求解信号瞬时频率。按式(5)求取主频Fm:

(5)

式中:f为频率;s(f)为希尔伯特变换幅度谱。

如上所述,从筛选后的每条记录中提取10个特征值,所有记录特征共组成1 223×10的波形特征向量矩阵。

3 基于LogitBoost决策树集成分类器的天然地震与人工爆破事件分类识别

3.1 基于LogitBoost机制的决策树集成分类器构建

Boosting算法是一种经典集成学习算法[27],被广泛应用于分类识别。该算法用初始权重从初始训练集训练出一个基分类器,根据基分类器的训练误差对训练样本分布进行调整,使得先前基分类器分类错误的训练样本在下一轮的训练中被选中的概率增加。如此重复进行,直到完成T轮训练,将这T个基分类器组合成一个强分类器,使强分类器在样本集上具有最小的预测损失,从而获得比单个学习器更好的性能。

Boosting算法中最经典的是AdaBoost,它采用的损失函数为ELOSS,损失随分类错误呈指数变化,容易受样本中噪声影响而产生过拟合。针对这一不足,Friedman提出LogitBoost[28]。LogitBoost采用二项式对数似然损失函数LLOSS,损失随分类错误呈线性变化,对样本噪声较不敏感,泛化能力更强。大量研究表明,利用LogitBoost进行分类识别能取得较好的效果[29-30]。LogitBoost对基分类器要求不高,只需要比随机猜测略好。本文选用决策树作为基分类器,决策树是一种非参的有监督学习方法,可从有特征和标签的训练样本数据中学习出决策规则,自动构造决策树进行分类。基于LogitBoost的决策树集成分类器构建过程如下:

步骤(1):选取训练样本集d=(xi,yi),i=1,2,…,N,其中N为样本个数；xi为判别变量,xi=(xi,1,xi,2,…,xi,m),m=10;yi为类别标签,设1为正类(天然地震),-1为负类(人工爆破)。

步骤(4):样本权重的更新。

(6)

(7)

(8)

3.2 分类器识别流程

为剔除干扰严重的记录,有效提取地震信号特征,在识别前先对地震记录进行筛选,基于STA/LTA进行P波初至拾取,截取初至点前10 s至后100 s共110 s记录,计算P波与S波最大振幅比;其次,进行基于S谱能量曲线的改进EWT,将地震信号分解为5个IMF,视最后一个IMF为噪声,取前4个IMF的香农熵和对数能量熵,并计算由前4个IMF重构的地震记录主频,每条记录共提取10个特征值;最后,对所有记录特征组成的波形特征向量矩阵,以决策树为基分类器,基于LogitBoost机制合成强分类器进行事件识别。算法流程如图1所示。

图1 本文分类识别算法流程图Fig.1 Flow chart of the proposed classification and recognition algorithm

3.3 分类器识别效果评价

为验证算法的稳定性及泛化能力,采用交叉验证和独立样本测试方法进行检验。其中交叉验证采用非重叠的随机选取划分策略,训练集和测试集独立,选取10%～100%的样本为训练集,其余为测试集(100%样本时,全部样本数据既做训练集也做测试集);独立样本测试是用新数据集进行测试。为评估算法性能,根据表1 所示的分类识别结果混淆矩阵,计算敏感度SE、特异性SP值和准确率ACC。

表1 地震事件分类识别结果Table 1 Classification and recognition results of seismic events

(9)

敏感度SE也称召回率,代表所有正例被正确分类的比例,能衡量分类器对正类的识别能力;特异性SP表示所有负例被正确分类的比例,能衡量分类器对负例的识别能力;准确率ACC是对分类器整体正确率的评价。SE越高,越容易鉴定出负类,即越灵敏;SP越高,则越不容易误分类,即筛选能力强。准确率、召回率、特异性是分类器的常用性能度量指标,一般认为具有高准确率、高召回率和高特异性的分类器是一个优秀的分类器。

4 实验及分析

4.1 改进EWT信号分解

为验证改进EWT的信号分解效果,分别对天然地震和人工爆破记录进行分解。图2和图3分别为经筛选(初至点前10 s至后100 s)并归一化后的天然地震和人工爆破记录,采样率均为100 Hz。天然地震是2020年8月16日福建安溪(25.391°N,117.742°E)发生的ML2.2微地震在QZH 台站的记录波形。人工爆破记录是2014年8月2日广西武鸣地区(23.006°N,108.21°E)发生的M2.6人工爆破事件在DHX台站的记录波形。图中可见,人工爆破持续时间相对较短,天然地震持续时间更长。基于此波形记录可计算P波与S波最大振幅比ps。

图2 天然地震记录Fig.2 Natural seismic record

图3 人工爆破记录Fig.3 Artificial blasting record

以图2所示天然微震波形记录为例进行改进EWT分解,所得S谱能量曲线绘于图4。该曲线能较清晰地描绘各频率点处原信号的能量大小,可见天然地震能量主要集中在低频端。

图4 S谱能量曲线Fig.4 S-transform spectrum energy curve

基于该曲线进行频带分割得到子频带,如图5所示；传统EWT频谱分割结果如图6所示。可见,传统EWT以傅里叶频谱的极值点为频谱分割边界,频谱分割杂乱无章,从而无法有效区分有效波和噪声;而本文改进EWT的子频带能根据频率和能量大小将原信号频谱分割开,便于信号分解。

图5 改进EWT频带分割结果Fig.5 Intersection of frequency spectrum using improved EWT

图6 传统EWT频带分割结果Fig.6 Intersection of frequency spectrum using traditional EWT

频谱分割后,基于改进EWT和传统EWT所得的IMF分别如图7、图8所示。改进EWT能较好地根据信号的频率和能量进行模态分解,如图7中IMF1～IMF4主要为有效波,而IMF5主要为高频随机噪声,将其视为噪声去除有利于实现有效波与随机噪声的有效分离。而传统EWT分解信号散乱,无法实现有效波与噪声的有效分离。用改进EWT对图3中的人工爆破信号进行分解,所得IMF如图9所示。同理,前4个IMF为有效信号,最后一个IMF为噪声。

图7 改进EWT分解IMF结果(天然地震)Fig.7 Decomposed IMFs by improved EWT (natural event)

图8 传统EWT分解IMF结果(天然地震)Fig.8 Decomposed IMFs by traditional EWT (natural event)

图9 改进EWT分解IMF结果(人工爆破)Fig.9 Decomposed IMFs based on improved EWT (manual blasting)

对IMF1～IMF4提取香农熵和对数能量熵,将IMF5视为噪声舍弃。再对IMF1～IMF4进行重构,计算重构信号主频,获取每条记录的特征属性值。

4.2 基于LogitBoost的决策树集成分类器性能

将1 223条记录提取的特征组成1 223×10的特征矩阵,采用基于LogitBoost机制的决策树集成分类器进行分类识别,迭代次数100次,算法基于Weka开源平台实现,环境为Matlab2016Ra,Win10系统。

基分类器个数是集成分类器的重要参数之一,理论上,样本数据量增大,基分类器个数应适当增加。如图10所示,以训练样本比例70%为例,以决策树为基分类器,将不同数量基分类器集成为强分类器,得到的分类准确率均在93.1%以上,说明本文算法有较好的分类效果。随着基分类器数量增加,分类准确率有一定提高,但分类器数量越多,识别时间也将增加。综合考虑算法时间与性能,本文以10个基分类器组成强分类器。

图10 不同数量基分类器对应的识别准确率Fig.10 Recognition accuracy corresponding to different amounts of base classifiers

采用交叉验证检验,随机抽取一定比例(10%～100%)的数据用于训练,其余数据用于测试,且同一训练样本比例随机抽取、运行100次,计算SE、SP和ACC值,并取均值,得到的分类识别结果列于表2。

表2 不同训练样本比例分类结果(单位：%)Table 2 Classification results with different proportions of training samples (Unit:%)

由表2可见,本文集成分类器方法中训练集的划分具有较好的鲁棒性,即训练样本数量对分类结果影响较小,训练样本数量取不同比例时,分类准确率ACC较高,且差异性较小,同时SE、SP也较高,说明分类器具有较好的稳定性。

4.3 与其他分类器对比

为了进一步分析分类器的性能,分别用BP神经网络、SVM、AdaBoost等分类算法进行分类识别,并与本文分类器的识别效果进行对比。使用BP神经网络分类时,为增强网络分类性能,采用双隐层,每个隐层4个节点,网络迭代100次,学习率为0.1,训练误差设置为10-5;SVM分类时,选用RBF核函数,选择有最高分类准确率的惩罚参数C和核函数参数g为最佳参数;AdaBoost也采用决策树为基分类器,基分类器数为10个。同上,选择不同的比例样本,分别用BP神经网络、SVM、AdaBoost、本文算法各运行100次,取均值后得到的分类准确率如图11所示。

图11 不同分类算法在不同样本比例下的分类准确率Fig.11 Classification accuracy of different algorithms under different proportions of samples

由图11可见,本文方法与其他3种方法相比有1%以上的分类准确率提升。SVM 基于结构风险最小化进行网络训练,BP神经网络基于经验风险最小化进行网络训练,SVM训练结果略优于BP神经网络。AdaBoost集成分类器和本文方法均采用组合策略,能减小单个分类器误差,具有较高识别率;二者对训练集的划分均具有较好稳定性,即训练样本数量多少对分类结果影响较小。BP神经网络和SVM对训练集数量比较敏感,当训练样本较少时分类准确率较低,随着训练样本增加,准确度提升。因此,当样本数量较少时更适合采用集成分类器,而Logitboost由于采用二项式对数似然损失函数,抗噪性较强,比AdaBoost集成分类器具有更好的性能,分类正确率更高一些。

为进一步验证本文分类器的实用性及稳定性,对另外两组独立的样本数量不同的数据集进行分类识别。其中一组:人工爆破为2014年11—12月广西田阳、2015年2—3月广西柳江的ML1.7～2.7矿藏开采记录,共452条,采样率100 Hz;天然地震为2020年8—10月ML1.8～2.8的天然微震记录,共481条,采样率100 Hz。用本文方法识别该数据集,准确识别人工爆破记录422条,识别准确率为93.4%;准确识别天然地震记录455条,识别准确率为94.6%,整体识别准确率94%。另外一组:天然地震为2021年5—8月187条ML3.0以下天然微震记录,人工爆破为2018年4—5月广西临桂202条采矿爆破记录。对该数据集进行分类预测,结果显示对天然地震的识别准确率为94.5%,人工爆破的识别准确率为93.4%,整体识别准确率94.2%。

实验表明,本文所提分类器的整体分类准确率较高,达94%,且训练样本数量对分类结果的影响较小,说明分类器模型具有较好的稳定性,能解决样本不足的问题。

5 结论

本文利用S变换优良的时频特性,提出了基于S谱能量曲线的改进EWT,以有效获取地震信号特征;研究了基于LogitBoost机制的集成分类器分类识别算法,并基于真实的天然微震和人工爆破记录对不同算法性能进行测试,得到如下结论:

(1) S谱能量曲线能真实完整地反映信号能量和频率特征,从而根据原信号的频率和能量完成信号分解,更好地提取信号特征。

(2) 根据原信号P波与S波最大振幅比、分解后IMF的香农熵和对数能量熵,及去噪后重构信号主频等特征,能有效区分天然和人工爆破的不同,可作为识别天然和人工爆破的特征参数,应用于实际的地震事件识别系统。

(3) 基于LogitBoost的决策树集成分类器准确率要优于AdaBoost、BP神经网络和SVM方法,且拥有1%以上的性能提升,说明本文算法有效可靠,具有一定的实用性。

(4) 通过集成学习方法可提高分类效果,基于LogitBoost机制的集成分类器方法具有较好的鲁棒性,即训练样本数量多少对分类结果影响较小,能有效解决样本不足的问题,适用性较强。

文中对分类器相关参数的研究不够细致,下一步将进行分类器相关参数优化,如基分类器种类、迭代次数等参数对分类结果的影响,从而获取最优参数,提高分类器分类效果。