王彤洲,崔春生,刘双峰,郭德月,张 健

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

0 引言

爆炸常常伴随着高温、高压、破片、强光和电磁脉冲等现象,在对爆炸场冲击波参数进行测试时,这些都是测试信号中噪声的来源[1]。爆炸场测试环境复杂,测试信号易受多种因素的噪声影响。爆炸冲击波具有上升沿陡峭、超压峰值高、正压作用时间短、带宽较宽等特点,冲击波信号在达到超压峰值后呈现指数式的衰减。冲击波的测量实质上是对冲击波的超压峰值、正压作用时间和比冲量的测量,是爆炸破坏效应和爆炸威力评估的重要手段。

传统的爆炸冲击波数据处理方法[2]主要有低通滤波法、小波分析法、最小二乘指数拟合。文献[3—4]对四种常见的IIR滤波器与不同窗函数的FIR滤波器进行了对比分析,并对最佳截止频率进行了探究。文献[5]运用小波变换,对适用冲击波信号的最佳分解小波基进行了探讨。

经验模态分解[6](EMD)适用于非线性、非稳定性信号处理,根据被分析信号本身的特点自适应地选择频带,在滤波去噪、机械故障检测等多个领域取得很好的应用;但是其仍存在模态混叠、虚假模态等问题。文献[7]指出间断事件干扰和密集模态相互作用是造成模态混叠的根本原因,并提出了3种改进方法。实测爆炸冲击波信号自身就属于一个幅值较大、呈指数衰减的间断事件,并且存在多种类型的噪声干扰,必然会产生模态混叠现象,而盲目地取舍模态函数可导致重构信号中部分有效信号丢失[8-9]。文献[10]结合EMD自适应性和小波分析的理论,提出经验小波变换(EWT),EWT自适应性好且具有完备的理论基础。针对上述情况,本文提出基于EWT-SVD联合算法的冲击波降噪方法。

1 经验小波变换及奇异值分解概述

1.1 经验小波变换(EWT)

经验小波变换(EWT)[11]借鉴了经验模态分解的自适应性又结合小波分析的理论,解决了EMD存在的模态混叠问题。经验小波变换的目的是把信号f(t)分解成N+1个本征模态函数fk(t)之和。本质上是对信号频谱进行自适应分割,从而建立合适的小波滤波器组,加小波窗后提取紧支撑傅里叶频谱的调幅-调频(AM-FM)成分。

1.2 奇异值分解(SVD)

奇异值分解的信号降噪方法[12-13]是先基于相空间重构理论,将待处理信号Y=[y(1),y(2),…,y(n+m-1)]构造成Hankel矩阵:

根据奇异值分解理论,对于给定的秩为r的m×n维矩阵H,则H的奇异值分解形式为

式(2)中,U,V分别为m×m,n×n维酉矩阵;Σ为r×r维对角阵,其对角线元素为矩阵H的非零奇异值σi,且以降序排列,即σ1≥σ2≥…≥σr。若源信号Y是由有用信号和噪声共同组成,则矩阵H的奇异值σi可以反映信号和噪声能量集中的情况。前p个较大的奇异值将主要反映有用信号,较小的奇异值则主要反映噪声,把这部分反映噪声的奇异值置零就可以去除信号中的噪声。将重构后的矩阵对应的项相加取平均值,就可以还原出降噪后的信号。

1.3 爆炸冲击波压力时间模型

冲击波超压随时间变化规律的数学表达式为Friedlander表达式:

式(3)中,ΔP为峰值压力,τ+为正压作用时间,b为衰减系数。

无限空中爆炸冲击波超压ΔP使用萨多夫斯基公式[14-15]计算:

马赫反射区的超压峰值ΔPM为

ΔPM=ΔP(1+cosα0),

(6)

式(6)中,ΔP是无限空中或近地爆炸对应比例距离处的超压峰值,α0为冲击波入射角。

正压持续时间为

式(7)中,r为测点到爆心的距离(m),ω为装药量(kg)。

2 基于EWT-SVD联合算法的冲击波降噪方法

2.1 EWT分解效果

冲击波信号与振动信号不同,冲击波信号周期性差,且存在较大的趋势项。与EWT方法在振动信号消噪中的应用不同,针对冲击波信号的特点,直接使用EWT对冲击波信号进行处理,分解得到3个分量M1、M2、M3,如图1所示。能明显观察到在信号端点处存在幅值很大的飞翼现象,若在此基础上对信号进一步作消噪处理,将会导致超压峰值出现较大误差、上升沿时间变长,而在冲击波数据处理中,冲击波的超压峰值是其关键的评价指标,因而不能将EWT方法直接运用在冲击波超压数据处理中。因而本文提出使用指数拟合的方法去除信号的趋势项,再在此基础上进一步作数据处理。

图1 冲击波信号直接进行EWT分解结果Fig.1 Shockwave signal directly performs EWT decomposition results

2.2 基于EWT-SVD联合算法的冲击波降噪原理

通过对冲击波信号进行指数拟合,用原始信号减去拟合后对应点的值,避免因冲击波信号上升沿陡峭、有效信号能量远大于噪声的特点而带来的较大的飞翼现象;再对上述信号使用EWT算法进行分解,随后对分解得到的分量采用SVD算法进行降噪重构,最终实现噪声的滤除。

2.3 降噪方法计算步骤

步骤1 对冲击波信号x(t)利用最小二乘指数拟合,拟合函数为Friedlander方程,得到拟合后的函数p(t)。

步骤3 对分解得到的每个MRA分量都进行相空间重构得到其对应的矩阵H,对矩阵H进行奇异值分解,通过对奇异熵增量的变化趋势的判断,保留前k个奇异值,其余的奇异值进行置零,然后进行逆变换,得到去噪后的矩阵H,将矩阵对应位置上的值相加取平均,得到去噪后的MRA分量。

步骤4 将所有去噪后的MRA分量及拟合所得信号p(t)相加,得到降噪后的冲击波信号。

3 算例分析

3.1 仿真计算

3.1.1仿真模型的建立

根据经验公式构建ω=60 kg,r=8 m处的理想冲击波信号P(t),并加入两个不同频率的模态信号从而构造出仿真信号Y(t),其表达式为

Y(t)=P(t)+0.01sin(600 000πt)+
0.01sin(120 000πt)

,

(8)

式(8)中,两个高频信号代表测试冲击波信号中的两种高频模态,并在此基础上加入已知信噪比(SNR)为20 dB的高斯白噪声。

对信号进行傅里叶变换得到它对应的频谱图,仿真信号及其频谱图如图2、图3所示,从图中可以看出构造的仿真模型符合冲击波信号主要集中在低频部分,且高频部分符合高斯白噪声的分布的特点。

图2 仿真模型Fig.2 Simulation model

图3 仿真模型频谱图Fig.3 Simulation model spectrogram

3.1.2EWT-SVD分解效果

使用本文提出的方法对仿真信号进行处理,待处理的冲击波超压信号是由离散点集构成的,曲线拟合的目标就是寻找这些点所对应的时间与其具体值之间的函数关系P=f(t)。在拟合过程中只对冲击波到达后的衰减过程进行拟合,得到的系数如表1所示。

表1 指数拟合所得Friedlander方程参数Tab.1 Parameters of the Friedlander equation obtained by exponential fitting

拟合所得数值与理论公式计算系数差距极小,很好地表示了冲击波信号的变化趋势。对减去拟合后的数据进行EWT处理,其分解结果及频谱图如图4、图5所示。

图4 EWT分解结果Fig.4 EWT decomposition results

图5 EWT分解结果频谱图Fig.5 EWT decomposition results in a spectrogram

该冲击波信号分解出了两个模态,分别对应在仿真信号中添加的两个高频信号。与图1对比可以看出飞翼现象得到明显抑制。对分解得到的模态进行相空间的重构,对于数据长度为N的信号重构矩阵时选取m=N/2为重构矩阵的行数,根据奇异熵增量选取重构奇异值的阶数,得到降噪后的模态,其频谱如图6所示。

图6 降噪后的分量频谱图Fig.6 Component spectrogram after noise reduction

将降噪后的模态与拟合数据进行相加得到最终降噪后的数据,并且有效地分离出加入的两个高频模态。作为比较使用EEMD方法对仿真信号进行分解,得到13个分量,与EEMD方法相比基于EWT-SVD联合算法的冲击波降噪方法很好地分离出加入的模态。对爆炸场冲击波中包含的不同种类的噪声的特性认识及区分具有一定的指导意义。

从去噪的能力方面对本文方法进行评估,已知对信号添加了信噪比为20.073 6 dB的高斯白噪声,对降噪后的信号求取信噪比其值为19.913 4 dB,证明了该方法在降噪方面的有效性。

3.2 试验验证

选取某次装药量60 kg,装药高度1.5 m,测试点距爆心8 m的试验数据,对该信号使用本文方法进行处理,首先对信号进行拟合,得到拟合后的Friedlander方程参数ΔP为0.351 8 MPa,τ+为0.004 s,b为1.113。

作为对比,对实测信号直接进行EWT分解,其结果如图7所示。对减去拟合数据后的实测信号进行EWT分解,所得分量如图8所示。

可以看出各个分量的飞翼现象得到了明显抑制。将所得的3个分量分别进行相空间重构,对重构后的矩阵进行奇异值分解,计算相应的奇异值熵增量。以分量2为例,所对应的奇异值熵增量如图9所示。

可以看出24阶之后的奇异值熵增量小于0.025,将24阶之后的奇异值置零,然后进行分量的重构。SVD降噪前的分量频谱如图10所示。重构后的3个分量的频谱如图11所示。

图7 直接EWT分解结果Fig.7 Direct EWT decomposition results

图8 减去拟合后的信号EWT分解结果Fig.8 Subtract the result of the EWT decomposition of the signal after subtraction

图9 分量2的奇异值熵增量图Fig.9 Singular value entropy delta plot of component

图10 SVD降噪前各分量频谱图Fig.10 Spectrogram of each component before SVD noise reduction

图11 SVD降噪后各分量频谱图Fig.11 Spectrogram of each component after SVD noise reduction

通过对比,明显观察到噪声得到了极大地滤除。将分量及拟合信号相加,得到最终降噪后的冲击波信号。作为比较对原始信号进行6阶贝塞尔滤波,截止频率100 kHz。原始信号和两种方法滤波后的结果如图12所示。基于EWT-SVD联合算法的冲击波降噪方法很好地还原出了冲击波信号的超压峰值、上升时间以及变化趋势且滤除了试验信号中的噪声干扰。从6阶贝塞尔滤波在超压峰值处可以观察到峰值大小的改变,且整个信号存在相位的偏移。该方法在提取超压峰值、抑制相位偏移、还原上升沿时间方面有着更好的优势。

图12 联合降噪后的结果Fig.12 Results after joint noise reduction

4 结论

本文提出了基于EWT-SVD联合算法的冲击波降噪方法。该方法首先对冲击波信号进行指数拟合,用原始信号减去拟合后对应点的值;对上述信号使用EWT算法进行分解,随后对分解得到的分量采用SVD算法进行降噪重构,实现噪声的滤除。仿真及试验验证表明,该方法可有效抑制相位偏移,更好地还原出冲击波信号的超压峰值、上升沿时间,可用于冲击波信号的降噪。该方法不足之处为对相空间重构后的矩阵进行SVD算法计算量过大,尚需改进。