融合CEEMDAN-GMM和小波阈值方法的爆炸冲击波降噪算法

杨志飞，崔春生，杜桂云，刘双峰，赵海霞

(1 中北大学省部共建动态测试技术国家重点实验室，山西 太原 030051；2 中北大学电气与控制工程学院，山西 太原 030051；3 中北大学环境与安全工程学院，山西 太原 030051)

0 引言

爆炸冲击波信号在军事、航天等领域具有重要意义[1-3]。在采集爆炸冲击波超压信号的过程中,由于高温环境、磁场干扰、测试设备误差等多种因素的影响,实测冲击波信号中往往混有大量环境噪声和系统自身的高频噪声。为了提高信号处理的准确性,人们对爆炸冲击波信号的降噪方法进行了广泛研究[4-18]。

张衍芳等[6]针对不同类型的典型实测冲击波信号进行了滤波处理,发现傅里叶变换等传统滤波方法无法完全消除信号中的毛刺。刘嘉慧[7]对实测数据曲线进行了贝塞尔滤波与巴特沃斯滤波处理,结果表明,截止频率为40 kHz的6阶贝塞尔滤波器在冲击波信号中去噪效果较好。由于传统滤波方法在处理非平稳信号时容易造成信号失真现象,如傅里叶变换只能在频域内进行分析,若信号在时域上某处发生突变,则无法分辨信号的尖峰是突变还是噪声导致,限制了其在爆炸冲击波信号降噪中的应用[8]。近年来,基于经验模态分解(EMD)的信号处理方法在信号降噪领域得到了广泛关注。自适应白噪声完整集合经验模态分解(CEEMDAN)是经验模态分解(EMD)的一种改进方法,通过引入白噪声信号,增强了模态分解的鲁棒性,能够更准确地分解非线性、非平稳信号,但EMD在去噪过程中分解出的固有模态函数分量之间出现模态混叠现象[15],通过经验模态分解去除高频模态分量的方式实现去噪,然而这种“一刀切”式的降噪方式会在降噪的同时消除信号中的有用成分,造成信号失真,对去噪效果会产生影响。梁晶等[16]采用CEEMDAN-小波阈值联合降噪算法对含噪模型和实测数据进行处理,发现该方法能够有效去除冲击波信号中的高频噪声,提高信号质量。张冉等[9]采用CEEMDAN-SG联合降噪算法解决了其他传统算法在信号重构后局部信号与原始信号差异过大的问题。张健等[10]采用CEEMDAN-PCA-CEEMDAN算法,利用PCA对CEEMDAN算法分解后的IMF分量正交化,以抑制模态混叠现象。

尽管上述方法在一定程度上提高了爆炸冲击波信号的降噪效果,但在信号失真、模态混叠等方面仍存在一定问题。因此,文中提出了一种基于自适应白噪声完整集合经验模态分解(CEEMDAN)-高斯混合模型(GMM)和小波阈值方法的爆炸冲击波降噪算法。首先,通过CEEMDAN对原始爆炸冲击波信号进行分解,得到多个本征模态分量(IMF)。然后,采用GMM聚类算法对IMF进行分类,并利用赤池信息准则(AIC)和轮廓系数(SC)确定最佳分类类别。剔除高频噪声类别的IMF分量后,对受噪声干扰的IMF分量采用小波阈值方法进行滤波去除噪声。最后,将滤波后的IMF分量与未受噪声干扰的IMF分量一起重构,从而获得降噪后的爆炸冲击波信号。

与文献[13,15]相比,最大的不同便是文中在获取了模态分解之后采用了GMM聚类算法对其进行分类并对此进行小波阈值滤波、重构等。

1 信号降噪算法实现

1.1 CEEMDAN-GMM

CEEMDAN-GMM算法是GMM算法[11]针对CEEMDAN降噪时容易造成信号失真所改进的算法,虽然CEEMDAN方法在信号降噪过程中减轻了EMD方法的模态混叠问题和避免了残余噪声影响,但该方法容易造成部分有效信号信息缺失,而CEEMDAN-GMM降噪算法充分利用了CEEMDAN在处理非平稳信号方面的优势,结合GMM聚类方法在数据分析领域的强大能力,而且能充分发挥小波阈值去噪法在不同尺度上分别处理信号中的噪声,从而更好地保留信号中有用信息的能力,实现对爆炸冲击波信号的高效降噪。

CEEMDAN-GMM算法流程如图1所示。其主要步骤为:

图1 CEEMDAN-GMM算法流程图Fig.1 CEEMDAN-GMM algorithm flowchart

1)在获取的原始爆炸冲击波信号x(t)中加入噪声序列,得到新信号x′(t),采用CEEMDAN算法分解新信号x′(t)获得m个不同的分量,每个IMF分量包含n个样本点IMFi(t);

2)利用高斯混合模型(GMM)[12]对IMF进行分类。首先,初始化高斯混合模型的参数,包括每个高斯分布的均值μk、协方差σk和权重πk。

3)使用EM算法[13]迭代更新GMM的参数,包括E-step和M-step两步。其中,E-step用于计算每个样本属于各高斯分布的后验概率γ(Zik);M-step根据后验概率更新高斯分布的均值μk、协方差σk和权重πk。其数学表达式分别为:

(1)

式中:i=1,2,…,n;k=1,2,…,m;N(xi|μk,σk)为高斯分布的概率密度函数。

(2)

(3)

(4)

4)反复执行步骤3),当满足收敛条件时停止迭代,完成聚类过程。根据后验概率将不同IMF分量划分至不同类别。

对于每个IMF分量IMFi(t),计算其属于每个高斯分布k的后验概率γ(Zik);

找出具有最大后验概率的高斯分布k*,将IMF分量IMFi(t)划分至对应的类别k*。

k*=arg maxkγ(zik)

(5)

式中:arg max表示使后验概率γ(Zik)取最大值的高斯分布。通过将每个IMF分量划分至具有最大后验概率的高斯分布,可以将不同IMF分量划分至不同类别。

5)计算不同IMF分量的方差贡献率校核聚类结果,剔除噪声类别的IMF分量,将剩余IMF分量重构获得降噪后的爆炸冲击波信号。

1.2 降噪效果评价指标

文中以信噪比RSNR和均方根误差MRMSE[15]作为评价指标来衡量信号降噪效果。信噪比用于衡量信号与噪声的比值,均方根误差用于衡量过滤后的信号与原始信号的紧密程度。其计算公式分别为:

(6)

(7)

式中:RSNR为信噪比;MRMSE为均方根误差;X为原始信号;Y为降噪后的信号;n为信号的序列长度。

一般情况下,降噪后信号RSNR越大,MRMSE越小,降噪效果越好。但降噪效果不仅要考虑评价指标,还要考虑降噪算法对信号超压峰值、正压时间等信号特征的影响较小。

1.3 相关算法信息

1.3.1 赤池信息AIC

为了在降噪过程中更好地区分信号与噪声,从而提高降噪算法的性能,文中引入赤池信息准则AIC和轮廓系数SSC作为评价指标来确定最佳聚类类别,从而有效地剔除噪声类别的IMF分量,保留信号中的有用信息,从而达到较好的降噪效果。

AIC[18]是评估统计模型的复杂度和衡量统计模型“拟合”资料之优良性的一种标准。

在一般情况下,AIC量化值NAIC可以表示为:

NAIC=2k-2lnL

(8)

式中:k是参数的数量,L是似然函数。NAIC鼓励数据拟合的优良性,但是尽量避免出现过度拟合的情况。因此优先考虑的模型应是NAIC值最小的那一个。

1.3.2 轮廓系数SSC

轮廓系数SSC[17]使用数据集中对象之间的相似性度量来评估聚类的质量,是簇的密集与分散程度的评价指标。其计算公式为:

(9)

式中:a为该数据帧与簇内其他数据帧的平均距离;b为该数据帧与距离它最近的另一个簇中样本的平均距离。

轮廓系数SSC的取值范围为[-1,1],其值越大,聚类结果越合理。如果轮廓系数SSC=-1,则表明该数据帧应该被分到其他的类;如果SSC接近0,则表明该数据帧在两个类的交叉点上。对全部样本的轮廓系数求平均值,可以得到聚类结果整体轮廓系数STCS,即:

(10)

STSC的取值范围为[-1,1],STSC值越大,聚类效果越好。

1.3.3 方差贡献率η

为保证GMM算法聚类分析结果的准确性,采用方差贡献率η[14]进行校核,方差贡献率计算公式为:

(11)

由式(11),根据方差贡献率η来判断类别数据受高频噪声干扰的高低,方差贡献率η越低说明此类别受到的噪声越大,反之为有效数据。

2 实测数据处理及分析

2.1 实验数据与预处理

实测爆炸冲击波信号来自60 kg级TNT爆炸冲击波试验,实测爆炸冲击波信号曲线如图2所示。采用CEEMDAN算法分解原始爆炸冲击波信号,分解后的曲线如图3所示。其中,高斯白噪声标准差与原始信号标准差之比为0.2,加入噪声的次数为100,最大迭代次数为800,获得17个不同频率的IMF分量和1个残余分量。根据图3可知,噪声对不同分量的影响逐渐降低,其中,IMF1～IMF8分量具备噪声分量的特征,IMF9～IMF10分量含有部分信号特征,但是也受到了噪声干扰,IMF11之后的分量具有和原始信号相似的波形。

图2 实测爆炸冲击波信号Fig.2 Measured explosion shock wave signal

图3 实测爆炸冲击波信号IMF分量Fig.3 IMF component of measured explosion shock wave signal

不同聚类类别的NAIC变化如图4所示,随着聚类类别个数不断增加,NAIC逐渐减小。由图4可知,NAIC曲线的转折聚类类别为5,在聚类类别为5之后,NAIC曲线下降曲线逐渐平缓,NAIC可以对模型复杂性和拟合优度进行权衡,以防止模型过于复杂而导致的过拟合。不同聚类类别的SSC变化如图5所示,SSC衡量了类别间的分离度和类别内部的紧密度。SSC值范围为-1～1,值越接近1,表示类别间的分离度越高,类别内部的紧密度越好。通过比较不同类别数量下的SSC值,可以找到一个最佳的类别数量,使得类别间的分离度和类别内部的紧密度都达到最优。同样地,当聚类类别为5时,SSC值最接近1。因此为了获得更好的聚类效果,将聚类类别设置为5。在CEEMDAN-GMM算法中,使用赤池信息准则(AIC)和轮廓系数(SSC)确定最佳聚类类别的意义在于:为高斯混合模型(GMM)提供一个合适的类别数量,以便在降噪过程中更好地区分信号与噪声。这两个评价指标有助于评估聚类结果的质量,并找到最佳的聚类划分。

图4 不同聚类类别AIC变化Fig.4 Changes in AIC of different cluster categories

图5 不同聚类类别SC变化Fig.5 Changes in SC of different cluster categories

根据GMM算法聚类分析结果,IMF1～IMF6分量划分为类别一;IMF7～IMF10分量划分为类别二;IMF11～IMF12分量划分为类别三;IMF13、IMF16、IMF17划分为类别四;IMF14、IMF15、IMF18属于类别五。类别一重构信号如图6所示。与原始爆炸冲击波信号图2相比,类别一IMF分量重构波形相似度较差且具有噪声信息特征,其他类别IMF分量重构波形相似度较好且包含较多原始信号信息。因此,可将类别一IMF分量重构波形初步定义为高频噪声序列,其他类别IMF分量重构波形定义为有效信息序列。

图6 类别一重构信号Fig.6 Reconstructed signal of category 1

根据式(11)计算不同类别的方差贡献率如表1。

表1 不同类别的方差贡献率Table 1 Variance contribution rates of different categories

由表1可知,类别五的方差贡献率最高,其值为75.89%,类别一的方差贡献率最低,仅有0.46%,说明方差贡献率的分析结果和GMM算法聚类的分析结论相同。因此,可以认定类别一位高频噪声序列,类别二～类别四是受噪声干扰的序列,需要降噪处理,类别五为有效信息序列。类别二～类别四经小波阈值去噪重构后的波形如图7所示。

图7 类别二～类别四经小波阈值去噪重构后的波形Fig.7 Waveforms reconstructed by wavelet threshold denoising for categories No.2 to No.4

将类别一高频噪声去除,重构类别二～类别四小波阈值降噪信号以及类别五IMF分量得到降噪后的纯净爆炸冲击波信号,同时与CEEMDAN降噪算法、CEEMDAN-小波阈值降噪算法相比的时域和频域图形如图8所示。

图8 爆炸冲击波实测信号降噪方法对比图Fig.8 Comparison of noise reduction methods for measured signals of explosion shock waves

2.2 降噪效果评价结果

根据式(6)～式(7)可以计算方法处理后各模型的RSNR指标和MRMSE指标如表2所示。

表2 降噪效果评价指标Table 2 Noise reduction effect evaluation indicators

由表2可知:三种方法均具有良好的去噪效果,但是,CEEMDAN-GMM和小波阈值算法的信号降噪效果最好,分别比CEEMDAN算法和CEEMDAN-小波阈值算法分别提高了10.19 dB和2.99 dB。由表3可知,该方法对毁伤评估指标的影响最小,相比于其他两种方法分别降低了0.42×10-4和0.14×10-4,降噪后更好地保留了原始信号信息,与原始信号具有较高的相似度。

表3 不同降噪方法对毁伤评估参数的影响Table 3 Impact of different noise reduction methods on damage assessment parameters

3 结论

文中所提出的融合CEEMDAN-GMM和小波阈值方法的爆炸冲击波降噪算法充分利用了CEEMDAN在非平稳信号处理中的优势,并结合了GMM聚类方法在数据挖掘领域的数据分析能力。与其他降噪方法相比,文中提出的CEEMDAN-GMM与小波阈值方法相结合的爆炸冲击波信号降噪算法在去除高频噪声方面具有较好的效果。通过与实测数据对比,CEEMDAN-GMM与小波阈值算法SNR值和RMSE值分别比CEEMDAN算法和CEEMDAN-小波阈值算法提高了10.19 dB和2.99 dB,降低了0.42×10-4和0.14×10-4。同时该算法能够在降低信号失真程度的同时保留信号中的关键特征。尚需完善的地方有:

1)文中采用的实测爆炸冲击波数据集数量有限,可能无法充分反映所有可能的噪声情况和信号特性。因此,需要扩大数据集,包括更多类型和来源的爆炸冲击波信号,以提高算法的通用性和稳健性。

2)文中固定了GMM聚类的参数和小波阈值方法的阈值参数,然而,这些参数可能会对降噪效果产生显著影响。需通过参数优选和策略调整来提高降噪性能。