刘 宇,李新娥,崔春生,李 顺

(1.中北大学省部共建动态测试技术国家重点实验室, 山西 太原 030051;2.中北大学电气与控制工程学院,山西 太原 030051)

0 引言

爆炸冲击波是典型的非线性、非平稳信号,上升沿陡峭且上升时间短,由于测试环境恶劣,测试信号混有大量高频干扰信号。冲击波的能量主要集中在低频段,干扰噪声分布在高频段,其能量幅值较小,频带较宽。因此为了得到爆炸冲击波的本征特点,有必要对采集到的爆炸冲击波信号进行去噪。目前,国内外爆炸冲击波信号的数据处理方法主要有Bessel[1-3]、小波阈值[1-3]、Butterworth[1-4]、经验模态分解[4-7](empirical mode decomposition,EMD)、希尔伯特黄变换[8-9](Hilbert-Huang,HHT)、EMD与小波阈值联合[10]。上述信号处理方法均各有利弊,其中,Bessel和Butterworth为经典滤波器,设置阶数和截止频率,进行低通滤波,去除高频信号,设计简单。Bessel针对整个信号,存在去除混叠的有用信号的缺点。Butterworth滤波结果会出现较大的相位滞后,不符合冲击波上升时间短的特点。EMD根据时间尺度特征将原始信号分解成为一系列IMF分量,能够反映出非平稳信号的局部特征。但也存在模态混叠、包络线拟合偏差、端点效应等缺点。直接剔除含有高频噪声的IMF,易丢失有用信号,造成失真。小波阈值、EMD与小波阈值联合都存在小波基、分解层数、阈值函数难以确定的问题。针对上述情况,提出了HHT与Bessel联合的信号处理方法,无需预先设定任何基函数,避免了经典低通滤波器针对整个信号,易去除有用信号的缺点,结合了Bessel设计简单和HHT自适应性好且适用于分析非线性非平稳信号的优点,该方法在剔除噪声的同时,还能较好地保留冲击波信号的本征特点。

1 HHT和Bessel基础理论

1.1 HHT(Hilbert-Huang)

HHT包含EMD和HT(Hilbert transform,希尔伯特变换)两部分,适用于分析非线性非平稳信号。

1.1.1EMD

EMD是时频域的处理方法,将原始信号分解成为一系列IMF分量,IMF分量是具有时变频率的振荡函数,能够反映出非平稳信号的局部特征。利用EMD对原始冲击波信号X(t)进行分解,步骤如下:

1) 求X(t)极大值和极小值,拟合上包络线u1(t)、下包络线v1(t)。

2) 求上下包络线的均值曲线m1(t)

3) 求中间曲线h1(t),h1(t)满足下面两个条件即为本征模态函数IMF1,进行下一步,若不满足,重复步骤1)—3)。

①所有极值点和零点的数量相同或相差1;②上下包络线的局部均值为0。

h1(t)=X1(t)-m1(t)。

(2)

4) 此时均值曲线m1(t)为IMF1的残余量,对m1(t)进行步骤1)—3)分解步骤,得到IMF2,以此类推,直到分解为单调信号为止,分解得到所有IMF和余量,最终的剩余分量为r(t)。

(3)

式(3)中,n为分量个数。

1.1.2HT(希尔伯特变换)

HT是分析和处理信号的一个重要的理论工具,可以提供90°的相位变换而不影响频率分量的幅度。IMF和余量的瞬时频率有着明确的物理意义。因此,EMD后, 对IMF和余量作希尔伯特变换(HT),求得IMF和余量的瞬时频率和瞬时幅值。

希尔伯特变换将s(t)与h(t)作卷积。

将IMFi(t)作希尔伯特变换,得式(6),继而求得每个IMF的瞬时频率wi(t)和瞬时幅值ai(t),如式(8)—式(9)。

所有IMF求和的时间-频率-能量三维频谱表示为

式(10)中,i表示本征模态分量分解次数。

余量r(t)同理表示为

(11)

每一个IMF和余量重构,将原始信号表示为

(12)

1.2 Bessel低通滤波器

Bessel滤波器具有平时延特性,较平坦的幅度特性,可以减小固有的非线性相位失真。在信号阻带衰减方面存在劣势,一般设计为高阶的滤波器来达到阻带衰减水平。

2 HHT与Bessel联合方法

HHT与Bessel联合的信号处理方法,结合了Bessel设计简单和HHT自适应性好且适用于分析非线性非平稳信号的优点。

利用 EMD将信号分解成一系列IMF和一个余量。首先对各IMF和余量与原始信号的互相关系数进行分析,初步判断含高频噪声的IMF,并采用Hilbert频谱分析最终确定含高频噪声的IMF分量;然后采用Bessel仅对含高频噪声的IMF处理;最后将处理后的IMF、纯净的IMF和余量进行重构,从而得到最终处理后的信号。流程图如图1所示。

图1 HHT与Bessel联合方法流程图Fig.1 Joint flow chart of HHT and Bessel

3 试验验证

3.1 试验方案

自由场爆炸冲击波测试过程中,针对测试点距离爆点近,采集的冲击波信号信噪比差的问题,有必要对距离爆心较近的测点处采集到的信号进行去噪处理。60 kg球形TNT,比例距离为 1.02 m/kg1/3,装药的高度距地面 1.5 m,三个测点距爆心距离相等。测点分布如图2所示。

图2 测点分布示意图Fig.2 Distribution diagram of measuring points

比例距离R的计算公式为

式(13)中,R为比例距离(单位:m/kg1/3),r为测试点离爆心的距离(单位:m),w为装药量的质量(单位:kg)。

炸药在无边界的空中爆炸,装药的高度H符合

式(14)中,H为装药的高度(单位:m),w为装药量的质量(单位:kg)。

采集的冲击波信号如图3所示,测点1、测点2、测点3的超压峰值分别为2.431、2.310、2.430 MPa。以测点1为例,详细分析信号处理过程,测点2和测点3的分析过程不再赘述。

图3 原始冲击波信号曲线Fig.3 Original shock wave signal curve

3.2 HHT与Bessel联合的信号处理过程

3.2.1EMD

原始冲击波信号被分解成9个IMF和1个余量,并按时间尺度从高频到低频依次分解,图4表示IMF1—IMF5,图5表示IMF6—IMF9和余量,重构误差为5.144 0×10-17。

图4 IMF1—IMF5Fig.4 IMF1—IMF5

图5 IMF6—IMF9和余量Fig.5 IMF6—IMF9 and residual

3.2.2互相关系数分析

IMF1—IMF9和余量与原始冲击波信号的互相关系数如表1所示,IMF1—IMF3互相关系数小于0.1。初步判定IMF1—IMF3是含噪IMF。

表1 互相关系数表Tab.1 Table of the cross-correlation coefficient

3.2.3确定含噪IMF

对IMF1—IMF3分别作Hilbert频谱分析,如图6—图8所示。通过三维时间-频率-能量谱分析频率在某段具体时间的能量谱密度,确定真实含噪IMF。

图6 IMF1的Hilbert频谱图Fig.6 Hilbert spectrum of IMF1

图7 IMF2的Hilbert频谱图Fig.7 Hilbert spectrum of IMF2

图8 IMF3的Hilbert频谱图Fig.8 Hilbert spectrum of IMF3

IMF以不同的频率显示原始信号特征。IMF1主要频率范围在0～400 kHz, IMF2主要频率范围在0～200 kHz,IMF3主要频率范围在0～70 kHz。确定真实含高频噪声的IMF为IMF1—IMF2。

3.2.4对含噪IMF进行Bessel处理

对含噪IMF进行Bessel低通滤波处理,处理后 IMF1和IMF2的Hilbert频谱图如图9—图10所示,去除了高频噪声干扰,并提取了IMF1和IMF2中的有用特征信息。

图9 处理后IMF1的Hilbert频谱图Fig.9 Hilbert spectrum of IMF1 after processing

图10 处理后IMF2的Hilbert频谱图Fig.10 Hilbert spectrum of IMF2 after processing

3.2.5重构结果分析

将处理后的IMF、纯净的IMF和余量进行重构,得到重构信号。与原始冲击波信号进行对比,处理前后信号曲线对比如图11所示,重构信号与原始冲击波信号趋势较为贴近。对重构信号进行Hilbert谱分析,如图12所示,可见去除了高频噪声干扰,保留了冲击波信号的本征特点。

图11 处理前后信号曲线对比图Fig.11 Signal curve comparison before and after processing

图12 重构信号的Hilbert谱Fig.12 Hilbert spectrum of reconstructed signal

3.3 滤波方法对比分析

对EMD、Bessel、小波阈值、EMD与小波阈值联合方法以及提出的HHT与Bessel联合方法进行对比分析。

SNR(信噪比)、RMSE(均方根误差)和超压峰值是信号处理的评价指标,SNR反映信号和噪声的比值, RMSE对信号处理前后的特大或特小误差尤为敏感。SNR越大,RMSE越小,超压峰值越贴近原始数据超压峰值,去噪精度越高。

SNR和RMSE计算公式如下:

(14)

(15)

式中,X表示处理前冲击波信号,x表示处理后冲击波信号,N为信号长度。

通过以上滤波方法分别对测点1、测点2、测点3的原始信号进行处理,结果如图13—图15所示,滤波指标如表2所示。

图13 测点1滤波方法对比Fig.13 Comparison of filtering methods of measuring point 1

图15 测点3滤波方法对比Fig.15 Comparison of filtering methods of measuring point 3

由表2分析可知:由于Bessel针对整个信号设置阶数和截止频率,易去除有用信号,超压峰值相比其他几种方法较低;小波阈值自适应差,预先设置不同的小波基、分解层数、阈值函数等参数,滤波效果不同,需要人为选择搭配参数,得到最优效果;EMD分解后,判断含高频噪声的IMF,直接去除,将其余纯净IMF重构,未筛选其中有用信息,去噪精度低;EMD与小波阈值联合提高了精度,EMD分解后,判断含高频噪声的IMF分量,利用小波阈值对含高频噪声的IMF分量处理,最后重构,得到处理后的信号,但同时存在小波基等参数的选择不当,造成去噪效果较差的结果。相较于其他滤波方法,HHT与Bessel联合采用EMD分解,通过互相关系数和Hilbert谱分析确定含噪IMF,采用Bessel低通数字滤波器对含噪IMF进行处理,最后重构,得到处理后的信号,在以上3个测点的信号处理中都占据优势, SNR更大,RMSE更小,处理后超压峰值接近原始信号的超压峰值。

表2 滤波指标对比Tab.2 Comparison of filtering indexes

4 结论

HHT与Bessel联合的冲击波信号处理方法,结合了Bessel设计简单和HHT自适应性好且适用于分析非线性非平稳信号的优点。采用EMD,通过互相关系数和Hilbert谱分析确定含噪IMF,并利用Bessel低通数字滤波器对含噪IMF进行处理,最后重构,得到处理后的冲击波信号。

相比于目前已有的EMD、Bessel、小波阈值、EMD与小波阈值联合,结合SNR、RMSE、超压峰值各项滤波指标综合分析,该方法SNR更大,RMSE更小,超压峰值接近原始信号的超压峰值。在剔除噪声的同时,还能较好地保留信号中的本征特点,提高了爆炸冲击波信号分析的准确性,去噪精度更高。