CN115456028A - 一种基于多模块滤波器集成的震动事件初至的自动识别算法 - Google Patents

Abstract

本发明公开了一种基于多模块滤波器集成的震动事件初至的自动识别算法，主要步骤：首先将连续记录原始信号进行预处理，然后进行切割，针对切割信号逐段进行振动能量有效频率分布判别；然后依次进行嵌入抗混叠低通滤波器、信号降采样、嵌入自适应去趋势项滤波器、嵌入中值滤波、嵌入无边缘延迟滤波等操作，再基于无边缘延迟滤波计算信号包络能量，最后嵌入鉴别滤波器，并将信号进行归一化，进而得到最终的识别特征函数。本发明建立了一套由多重滤波器集成的自动拾取算法，在降噪的同时，构建特征函数，一方面仅需要最基本的信号参数，另一方面有效降低了噪音干扰，提取灾害事件初至时刻，提高了自动拾取灾害事件发生的准确率。

Description

一种基于多模块滤波器集成的震动事件初至的自动识别算法

技术领域

本发明涉及灾害预警预测技术领域，尤其是一种基于多模块滤波器集成的震动事件初至的自动识别算法。

背景技术

目前，灾害事件的预警预测越来越依赖于动力学特征观测。自然界中的天然灾害如地震、火山、滑坡、崩塌落石、泥石流，抑或是人工灾害如矿山岩崩、管道破裂、高压气体泄漏均伴随着能量的震动传播，因此，通过仪器记录其振动时序特征，基于数学物理方法便可分析其灾害源的爆发特征。其中，最为基础的分析便是事件发生时刻的确定，即初至的拾取。理论上，由足够数量监测仪器(宽/窄频地震仪、长周期地震仪，加速度计、声波传感器等)组成观测台阵即可获取事件发生的初至时刻、事件的能量规模及时间的能量源方位角等特征信息。但是在实际的监测中，检测仪器往往数量较少，因此如何通过少量监测仪器便可获取事件的初值时刻一直是各领域的研究热点，同时也存在大量技术难点，如何从繁杂随机的震动记录中寻找出灾害事件所携带的源能量对于灾害的从而形成多尺度灾害的预警意义重大。

利用动力学振动信号揭示不同尺度灾害事件的方法的主要目的在于精确监测灾害事件时空分布特征，其基础在检测出灾害事件的初至，即灾害事件发生的时刻。目前针对于初至时刻的拾取主要是基于动力学时序信号的分析，通常灾害事件发生前期，检测传感器所记录的时序信息基本为背景噪音；灾害事件发生时刻，通常伴随着较为强烈的脉冲信号，进而成为源在介质中传播，直至被不同距离布设的传感器接收；而随着时间延续，信号强度随之渐进衰减，直至震动幅值减小至背景噪音尺度。关于初至拾取整体上是基于信号时域或频域内的包络强度、绝对振幅或信号的功率谱特征。目前较为流行的算法为短时平均/长时平均比算法(STA/LTA)与功率谱阈值方法(PSD)。

短时平均/长时平均比算法虽然可以拾取事件发生初至，但是过于依赖预先设置参数，尽管部分研究通过改变特征函数已修正算法的敏感度。此外，短时与长时时窗长度影响着阈值，低阈值可能导致许多错误初至触发(假阳性)，而高阈值可能导致丢失弱事件(假阴性)。为了改进短时平均/长时平均比算法过于依赖设定参数的不足，功率谱阈值方法应运而生。背景噪音的频谱含量通常在各频带较为均匀，而灾害事件的频谱含量却在某些特定频带内非常大，因此功率谱阈值算法的核心就是寻找频率含量异常大的时刻，拾取灾害事件初至。

针对功率谱阈值算法，虽然其前期参数设置较少，但是算法对于信号信噪比的要求较高。信噪比高的信号，例如强震能量，上述两种方法均有着不错的识别能力，但对于信噪比低，甚至是信号能量湮没与背景噪音的震动能量，如微震、火山微动、管道管涌等灾害事件，上述算法就会有很大程度的错误识别。因此，建立一种既包含降噪模块，又包含特征拾取模块的较为稳定的算法对于振动信号初至拾取，进而对灾害发生时刻的预警有着重要现实意义。

发明内容

针对现有的两种震动事件初至拾取算法存在的依赖参数多，识别低信噪比信号准确度低等问题，本发明提供一种基于多模块滤波器集成的震动事件初至的自动识别算法。

本发明的算法的核心为两部分内容，第一部分主要是进行信号的降噪，即抑制噪音与尖峰失真现象；第二部分主要是基于包络函数与鉴别函数构建特征函数，完成初至拾取。

本发明提供的基于多模块滤波器集成的震动事件初至的自动识别算法，具体步骤如下：

S1、将连续记录原始信号进行预处理，然后进行切割，针对切割信号逐段进行振动能量有效频率分布判别。具体如下：

S11、将连续记录原始信号进行预处理，包含基线调平、去均值及去趋势项步骤；

S12、将预处理后的原始信号进行切割，采用重叠搭接方式进行裁剪；

S13、针对切割信号逐段进行振动能量有效影响时间判别，设置时间的影响持时范围：d_min：灾害的最短持续时间；d_max：灾害的最长持续时间；N_d：持时可分解的时间窗数；

S14、针对切割信号逐段进行振动能量有效频率分布判别，选择低频频段，设置频带范围：f_low：低频截止频率；f_high：高频截止频率；N_f：频率可分解的频带窗数。

S2、嵌入抗混叠低通滤波器，滤波器表示为：

式中：n为滤波器阶数，ω_c为截止角频率；ω_p为通带边缘频率。

S3、信号降采样：在上一步骤低通滤波的基础上，进行降采样，引入降采样的比例系数R，表示为：

fsp＝2f_high

式中，fs为信号的原始采样频率；fsp为降采样后的采样频率；f_high为预先设定的高频截止频率。

S4、嵌入自适应去趋势项滤波器；该滤波器表示为：

B＝FIR(M,ω_c)

M＝round(d_max×fsp)

trend＝filtfit(B,1,x_r)

x_rd＝x_r-trned

式中，B为滤波器，M为截止角频率阶数，ω_c是截止角频率，d_max为预先设定的灾害的最长持续时间；f_low为预先设定的低频截止频率；trend表示去趋势的意思，意思是这个步骤是去趋势滤波；x_r为降采样后的信号；x_rd是趋势项滤波后的信号。

S5、嵌入中值滤波，中值滤波表示为：

y＝medfilt(x_rd,round(fsp×d_min/10))

式中，d_min为预先设定的灾害的最短持续时间。

S6、嵌入无边缘延迟滤波；该滤波器表示为：

y′＝filter(H,1,y)

式中，y’为无边缘延迟滤波后的信号；H为滤波器参数；y为中值滤波后的信号。

滤波器参数H设置如下：

H＝acos(2π×fc×th)

a＝1:length(th):0

k＝1,2,……,N_f

式中，f_c为[f_low，f_high]的中心频率，th是时间序列，a是振幅变换系数，表示从0到1，th的个数，k是频率点数，N_f为预先定义的频率可分解的频带窗数。

S7、计算信号包络能量：基于无边缘延迟滤波，将信号依据不同中心频率信号分为不同频段的子信号，进而计算每一段信号的包络能量E_k，计算公式如下：

k＝1,2,……,N_f

式中，H为信号的希尔伯特变换；y’为无边缘延迟滤波后的信号。

S8、先嵌入鉴别滤波器，然后建立特征函数；

鉴别滤波器的表达式如下：

D(m)＝filter(T,1,E_k)

m＝1,2,……,N_d

k＝1,2,……,N_f

式中，D为鉴别滤波后的信号；T为含有非脉冲信号的惩罚因子β的中间变量；E_k是信号包络能量，m是代表自然数序列,表示从第1个一直取到第N_d个，N_d是持时可分解的时间窗数；

中间变量T表示为：

T＝β×((1:L)-(L+1))

L＝round(D(m)×fsp)

D＝[d_min,N_d,d_max]

式中，β是非脉冲信号的惩罚因子，L是表示D的数据长度。

将信号进行归一化，进而得到最终的识别特征函数，特征函数的表达式如下：

CF＝10×(log₁₀(1+D/V)/log₁₀(2))

式中，CF为初至识别特征函数；V为信号的最大振幅，D为鉴别滤波后的信号。

与现有技术相比，本发明的有益之处在于：

(1)本发明的算法在较宽的振动时窗、信号频率、惩罚因子等参数范围内，均可以拾取稳定的初至时刻，不同的设定仅影响着事件的拾取分辨率，但对于重要事件的拾取结果基本一致。有效避免了特征函数对于阈值参数的依赖。

(2)本发明嵌入了多重滤波器，同时规避了背景噪音、高频仪器噪音及低频环境噪音的影响，有效提高了记录信号的质量，为计算提取特征参数提高了精度。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1、本发明方法基于多模块滤波器集成算法的自动拾取震动事件初至设计计算流程图。

图2、基于多模块滤波器集成算法的自动拾取震动事件示例图。图中，a是原始振动信号时间序列；b是抗混叠低通滤波处理的信号时间序列；c是降采样处理的信号时间序列；d是自适应去趋势项处理的信号时间序列；e是中值滤波处理的信号时间序列；f是无边缘延迟滤波处理的信号时间序列；g是多频段信号能量包络曲线；h是能量包络曲线叠加；i是归一化特征函数。

图3、不同滤波频段的特征函数及初至时间提取(d_min＝1s，d_max＝2s，N_d＝1，N_f＝1，β＝1)。

图4不同振动持时的特征函数及初至时间提取(f_low＝1Hz，f_hight＝2Hz，N_d＝1，N_f＝1，β＝1)。

图5不同频段个数(Nf)的特征函数及初至时间提取(f_low＝1Hz，f_high＝3Hz，d_min＝1s，d_max＝3s，N_d＝1，β＝1)。

图6不同惩罚因子β的特征函数及初至时间提取(f_low＝1Hz，f_high＝3Hz，d_min＝1s，d_max＝3s，N_d＝1，N_f＝1)。

图7不同噪音强度影响下的(SNR＝1,2,3,5)信号特征函数及初至时间提取。

图8地下管道破裂事件振动时序记录及初至时刻。

图9地震事件振动时序记录及初至时。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

如图1-9所示，本发明提供的基于多模块滤波器集成的震动事件初至的自动识别算法，具体步骤如下：

S1、设置时序信号基本参数：具体如下：

S11、将连续记录原始信号进行预处理，包含基线调平、去均值及去趋势项等步骤，如图2a所示。

S12、将预处理后的原始信号进行切割，为尽可能保证切割信号的连续性与能量泄露较少，采用重叠搭接方式(overlapping)进行裁剪，重叠度大于50％。

S13、针对切割信号逐段进行振动能量有效影响时间判别，通常灾害发生引发的能量为脉冲能量，此处设置时间的影响持时范围：d_min：灾害的最短持续时间；d_max：灾害的最长持续时间；N_d：持时可分解的时间窗数；

S14、针对切割信号逐段进行振动能量有效频率分布判别，通常灾害诱发振动能量的初至为压缩体波震相类型振动(P波)，体波不具有频散效应，即任何频段的能量传播速度相同，因此可选择低频频段，以保证信号的简约清晰。此处设置频带范围：f_low：低频截止频率；f_high：高频截止频率；N_f：频率可分解的频带窗数。

以下步骤以一段切割信号为例进行说明。

S2、嵌入抗混叠低通滤波器(Anti-aliasing FIR filter)：为消除高频噪音对于灾害能量产生振动信号的影响，该步骤设计嵌入低频滤波器(如图2b所示)，为保证信号在通带频域内的频响曲线最大限度平坦，而在阻带频域内频响函数逐渐趋近于零，可采用巴特沃斯低通滤波器(Butterworth low pass filter)，以保证信号观测频段内的信息显著。滤波器表示为：

式中：n为滤波器阶数，ω_c为截止角频率，即振幅下降为-3dB时的对应频率；ω_p为通带边缘频率。角频率与线频率的转换为：ω＝2πf.例如：信号的采样频率(fs)为100Hz，信号的最大频率即为50Hz，若要保留10Hz以下的能量，即滤波10Hz以上的能量，ω_c＝2*π*10＝62.8rad/s。

S3、信号降采样：为尽可能简化针对信号分析的计算量，并降低信号复杂度，保留的信号的重要初至信息，在上一步骤低通滤波的基础上，进行降采样(如图2c所示)，引入降采样的比例系数R，表示为：

fsp＝2f_high

式中，fs为信号的原始采样频率；fsp为降采样后的采样频率；f_high为预先设定的高频截止频率。

S4、嵌入自适应去趋势项滤波器(Adaptive de-trending filter)；

为消除低频噪音(气象噪音或仪器缓慢偏移)对于信号监测所带来的影响，该步骤设计嵌入去趋势项滤波器，缘于低频噪音主要会造成信号对于基线的持续偏离(如图2d所示)。该滤波器为线性相位滤波器与零相移滤波函数组成。滤波器表示为：

B＝FIR(M,ω_c)

M＝round(d_max×fsp)

trend＝filtfit(B,1,x_r)

x_rd＝x_r-trned

式中，B为滤波器，M为截止角频率阶数，ω_c是截止角频率，d_max为预先设定的灾害的最长持续时间；f_low为预先设定的低频截止频率；x_r为降采样后的信号；x_rd是趋势项滤波后的信号。

S5、嵌入中值滤波(Median filter)；

进一步消除信号当中由与监测仪器异常或数据传输引发的尖峰噪音，此步骤设计嵌入中值滤波器，旨在进一步获取平滑的信号(如图2e所示)，中值滤波表示为：

y＝medfilt(x_rd,round(fsp×d_min/10))

式中，d_min为预先设定的灾害的最短持续时间。

S6、嵌入无边缘延迟滤波(No edge delay filter)；

灾害作为振动源释放能量，形式多为脉冲信号，为增强脉冲信号的显著性，此步骤嵌入无边缘延迟滤波器，进一步减小背景噪音对于信号的干扰(如图2f所示)；滤波器表示为：

y′＝filter(H,1,y)

式中，y’为无边缘延迟滤波后的信号；H为滤波器参数；y为中值滤波后的信号。

滤波器参数H设置如下：

H＝a cos(2π×fc×th)

a＝1:th:0

k＝1,2,……,N_f

式中，f_c为[f_low，f_high]的中心频率，N_f为预先定义的频率可分解的频带窗数。例如定义信号的分析频段为[1Hz，3Hz]，N_f取值3，那么f_c则为[1.33Hz，2Hz，2,67Hz]。

S7、计算信号包络能量：

基于无边缘延迟滤波，将信号依据不同中心频率信号分为不同频段的子信号，进而计算每一段信号的包络能量E_k(如图2g所示)，计算公式如下：

k＝1,2,……,N_f

式中，H为信号的希尔伯特变换；y’为无边缘延迟滤波后的信号。

S8、先嵌入鉴别滤波器，然后建立特征函数；

(1)为了保证灾害事件引发的振动能量初至的辨识程度，此步骤嵌入鉴别滤波器，通过急剧上升和逐渐下降的速度增强信号衰减，同时再次减轻脉冲噪声的影响，并且引入一个非脉冲信号的惩罚因子β，用于提高滤波器压制背景噪音的作用，鉴别滤波器的表达式如下：

D(m)＝filter(T,1,E_k)

m＝1,2,……,N_d

k＝1,2,……,N_f

式中，D为鉴别滤波后的信号；T为含有非脉冲信号的惩罚因子β的中间变量；E_k是信号包络能量，N_d是持时可分解的时间窗数。

中间变量T表示为：

T＝β×((1:L)-(L+1))

L＝round(D(m)×fsp)

D＝[d_min,N_d,d_max]

式中，β是非脉冲信号的惩罚因子。

(2)为减小鉴别滤波后的信号中初至振幅强度损失，将信号进行归一化，进而得到最终的识别特征函数(如图2i所示)，特征函数的表达式如下：

CF＝10×(log₁₀(1+D/V)/log₁₀(2))

式中，CF为初至识别特征函数；V为信号的最大振幅，D为鉴别滤波后的信号。

本发明方法基于多模块滤波器集成算法的自动拾取震动事件初至设计计算流程见图1所示。

为验证本发明提出的识别算法在参数依赖及压制噪音两方面性能的效果，做如下试验：

针对识别特征函数对于参数设定的依赖研究，开展正交验证试验，可以通过控制信号持时特征(d_min、d_max、N_d参数)、信号频率特征(f_low、f_high、N_f参数)及非脉冲信号的惩罚因子β，三个类别参数进行振动初至时间对比。就信号的持时特征而言，灾害事件引发的振动能量信号的初至一般在10秒以内，因此设置d_min＝1s，d_max＝3～10s，N_d＝1～5。灾害事件信号不具有频散特征，任何频段的能量传播速度相同，因此可选择更为稳定清晰地低频频段，f_low＝1Hz，f_high＝3～10Hz，N_f＝1～5。而非脉冲信号的惩罚因子β可取1～3。参数设置如表1所示。

表1、正交验证试验特征参数

图3显示了在信号持时特征与非脉冲信号的惩罚因子β两类参数变量固定的情况下，改变信号频带宽度，识别特征函数的特征。整体而言，随着频带宽度的增加，识别特征函数的峰值逐渐增多，意味着拾取的振动事件逐渐增多。但是随着频带宽度变化，携带更多振动能量的主要事件均可以被识别。本试验所选取的频带宽度(1-10Hz)包含了大部分天然灾害(火山地震、崩滑灾害)震源及工程灾害(管道破裂、高压气体泄漏)所释放的能量频带。

图4显示了在信号滤波频段与非脉冲信号的惩罚因子β两类参数变量固定的情况下，改变信号持时宽度，识别特征函数的特征。随着持时的宽度改变，特征函数的峰值逐渐发生延迟现象，特别是对于持时宽度在7s以上。而较小持时宽度的特征函数，其形态基本一致。对于灾害源释放的振动，其能量释放是高频脉冲，因此震动时间较短，对于7s甚至更高的持时宽度而言，其对应的低频(<1Hz)振动源，并不是本发明所涉及的以爆发性震源特征的灾害。因此，持时宽度的选择并不会太多影响振动事件初至时刻的拾取。图5显示了在信号不同频段个数(N_f)对于特征函数的影响效应，对于选定的频段宽度，不同的频带个数，几乎对于特征函数的形态没有影响，尤其是在主要事件的初至时刻。但较少的频段个数可以减少算法运算时间。同样的情况也出现在不同持时个数(N_d)的情况。而图6则显示了信号持时特征与频率特征不变情况下，非脉冲信号的惩罚因子β对于特征函数的影响效应。对于较大的β值，仅可以实现强能量事件的初至拾取，而随着β值的逐渐减小，更多的弱能量事件会被会被识别。随着惩罚因子β的变大，也对于脉冲信号的选拔也为严格。

针对识别特征函数在不同噪音强度干扰下的稳定性，进行相同信号，不同噪音等级的验证测试。本试验通过对原始信号混合高斯随机白噪音，分别设定信噪比(SNR)分别为5、3、2、1。其信号时序与特征函数如图7所示。当信噪比为5时，噪音强度对于信号影响较小，信号的振动特征清晰，基于特征函数可以清晰拾取四次主要事件的初至时刻，与信号时序对应。随着信噪比降至3时，信号中的噪音含量增大，但特征函数的形态基本没有改变，四次主要事件的初至时刻依旧明显。当信噪比进一步减小至2时，特征曲线仅有微小改变。当信噪比为1时，即信号的振幅与噪音的振幅相等，这种信号通常质量较低，振动时序特征发生明显改变，信号几乎湮灭与噪音，此时特征曲线的形状较之前发生明显改变。但依旧可以拾取重要事件依然，特别是第一次事件的初至，依旧与高信噪比情况下保持一致。此结果说明本发明算法中的降噪模块有效还原信号形态。

为进一步应用本算法在提取振动事件初至，将其应用于不同震源诱发的类型振动信号的初至拾取，并与传统的短时平均/长时平均比算法(STA/LTA)进行比较。如图8所示，为地下水管破裂时形成微震信号，通过本发明提出的多滤波模块集成算法提取其能量初至时刻为5.03s，同时利用传统的短时平均/长时平均比算法(STA/LTA)常用的3阶阈值算法得到初至时刻为5.16s，而4阶阈值算法得到初至时刻则为5.27s。此信号有着良好的信噪比，因此噪音的干扰较小，两种算法的计算结果基本一致。

如图9所示，为日本北海道Iburiku以此震级为3.1的近震事件，因传地震仪埋设于附近污水处理工厂近地表松散覆盖层，因此地震时序记录，存在较强的干扰，地震事件几乎湮没与工厂噪音中，信号的信噪比较低。通过本发明提出的多滤波模块集成算法提取其能量初至时刻为3.17s，但是传统的短时平均/长时平均比算法(STA/LTA)不同阶参数阈值计算的结果却出现明显差异，3阶阈值算法得到初至时刻为1.29s，而4阶阈值算法得到初至时刻则为6.53s。此处震源距地震台站的距离约为3Km，而大区域的基岩基本为强风化质凝灰岩，初至P波波速在800～1000m/s，理论的初至时刻在3～3.75s，此理论值虽然为一维速度模型计算，但是和本发明提出的多滤波模块集成算法的计算结果接近，因此本发明算法计算较为合理。

综上，本发明建立了一套由多重滤波器集成的自动拾取算法，在降噪的同时，构建特征函数，一方面仅需要最基本的信号参数，另一方面有效降低了噪音干扰，提取灾害事件初至时刻，进一步提高了自动拾取灾害事件发生的准确率。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (6)

1.一种基于多模块滤波器集成的震动事件初至的自动识别算法，其特征在于，步骤如下：

S1、将连续记录原始信号进行预处理，然后进行切割，针对切割信号逐段进行振动能量有效频率分布判别；

S2、嵌入抗混叠低通滤波器；

S3、信号降采样：在上一步骤低通滤波的基础上，进行降采样，引入降采样的比例系数R，表示为：

fsp＝2f_high

式中，fs为信号的原始采样频率；fsp为降采样后的采样频率；f_high为预先设定的高频截止频率；

S4、嵌入自适应去趋势项滤波器；该滤波器表示为：

B＝FIR(M,ω_c)

M＝round(d_max×fsp)

trend＝filtfit(B,1,x_r)

x_rd＝x_r-trned

式中，B为滤波器，M为截止角频率阶数，ω_c是截止角频率，d_max为预先设定的灾害的最长持续时间；f_low为预先设定的低频截止频率；trend是去趋势滤波，x_r为降采样后的信号；x_rd是趋势项滤波后的信号；

S5、嵌入中值滤波，中值滤波表示为：

y＝medfilt(x_rd,round(fsp×d_min/10))

式中，d_min为预先设定的灾害的最短持续时间；

S6、嵌入无边缘延迟滤波；该滤波器表示为：

y′＝filter(H,1,y)

式中，y’为无边缘延迟滤波后的信号；H为滤波器参数；y为中值滤波后的信号；

S7、计算信号包络能量；

S8、先嵌入鉴别滤波器，然后建立特征函数，特征函数的表达式如下：

CF＝10×(log₁₀(1+D/V)/log₁₀(2))

式中，CF为初至识别特征函数；V为信号的最大振幅，D为鉴别滤波后的信号。

2.如权利要求1所述的基于多模块滤波器集成的震动事件初至的自动识别算法，其特征在于，步骤S6中，滤波器参数H设置如下：

H＝acos(2π×fc×th)

a＝1:length(th):0

k＝1,2,……,N_f

式中，f_c为[f_low，f_high]的中心频率，th是时间序列，a是振幅变换系数，表示从0到1，th的个数；k是频率点数，N_f为预先定义的频率可分解的频带窗数。

3.如权利要求1所述的基于多模块滤波器集成的震动事件初至的自动识别算法，其特征在于，步骤S8中，鉴别滤波器的表达式如下：

D(m)＝filter(T,1,E_k)

m＝1,2,……,N_d

k＝1,2,……,N_f

式中，D为鉴别滤波后的信号；T为含有非脉冲信号的惩罚因子β的中间变量；E_k是信号包络能量，m是自然数序列，表示从第1个一直取值到第N_d个，N_d是持时可分解的时间窗数；

中间变量T表示为：

T＝β×((1:L)-(L+1))

L＝round(D(m)×fsp)

D＝[d_min,N_d,d_max]

式中，β是非脉冲信号的惩罚因子，L是表示D的数据长度。

4.如权利要求1所述的基于多模块滤波器集成的震动事件初至的自动识别算法，其特征在于，所述步骤S1具体如下：

S11、将连续记录原始信号进行预处理，包含基线调平、去均值及去趋势项步骤；

S12、将预处理后的原始信号进行切割，采用重叠搭接方式进行裁剪；

S13、针对切割信号逐段进行振动能量有效影响时间判别，设置时间的影响持时范围：d_min：灾害的最短持续时间；d_max：灾害的最长持续时间；N_d：持时可分解的时间窗数；

S14、针对切割信号逐段进行振动能量有效频率分布判别，选择低频频段，设置频带范围：f_low：低频截止频率；f_high：高频截止频率；N_f：频率可分解的频带窗数。

5.如权利要求4所述的基于多模块滤波器集成的震动事件初至的自动识别算法，其特征在于，所述步骤S2中，抗混叠低通滤波器表示为：

式中：n为滤波器阶数，ω_c为截止角频率；ω_p为通带边缘频率。

6.如权利要求5所述的基于多模块滤波器集成的震动事件初至的自动识别算法，其特征在于，所述步骤S7具体是：基于无边缘延迟滤波，将信号依据不同中心频率信号分为不同频段的子信号，进而计算每一段信号的包络能量E_k，计算公式如下：

k＝1,2,……,N_f

式中，H为信号的希尔伯特变换；y’为无边缘延迟滤波后的信号。

Images