基于經(jīng)驗小波變換和奇異值分解的沖擊波降噪方法

王彤洲,崔春生,劉雙峰,郭德月,張 健

(1.中北大學省部共建動態(tài)測試技術國家重點實驗室,山西 太原 030051;2.中北大學電氣與控制工程學院,山西 太原 030051)

0 引言

爆炸常常伴隨著高溫、高壓、破片、強光和電磁脈沖等現(xiàn)象,在對爆炸場沖擊波參數(shù)進行測試時,這些都是測試信號中噪聲的來源[1]。爆炸場測試環(huán)境復雜,測試信號易受多種因素的噪聲影響。爆炸沖擊波具有上升沿陡峭、超壓峰值高、正壓作用時間短、帶寬較寬等特點,沖擊波信號在達到超壓峰值后呈現(xiàn)指數(shù)式的衰減。沖擊波的測量實質(zhì)上是對沖擊波的超壓峰值、正壓作用時間和比沖量的測量,是爆炸破壞效應和爆炸威力評估的重要手段。

傳統(tǒng)的爆炸沖擊波數(shù)據(jù)處理方法[2]主要有低通濾波法、小波分析法、最小二乘指數(shù)擬合。文獻[3—4]對四種常見的IIR濾波器與不同窗函數(shù)的FIR濾波器進行了對比分析,并對最佳截止頻率進行了探究。文獻[5]運用小波變換,對適用沖擊波信號的最佳分解小波基進行了探討。

經(jīng)驗模態(tài)分解[6](EMD)適用于非線性、非穩(wěn)定性信號處理,根據(jù)被分析信號本身的特點自適應地選擇頻帶,在濾波去噪、機械故障檢測等多個領域取得很好的應用;但是其仍存在模態(tài)混疊、虛假模態(tài)等問題。文獻[7]指出間斷事件干擾和密集模態(tài)相互作用是造成模態(tài)混疊的根本原因,并提出了3種改進方法。實測爆炸沖擊波信號自身就屬于一個幅值較大、呈指數(shù)衰減的間斷事件,并且存在多種類型的噪聲干擾,必然會產(chǎn)生模態(tài)混疊現(xiàn)象,而盲目地取舍模態(tài)函數(shù)可導致重構信號中部分有效信號丟失[8-9]。文獻[10]結合EMD自適應性和小波分析的理論,提出經(jīng)驗小波變換(EWT),EWT自適應性好且具有完備的理論基礎。針對上述情況,本文提出基于EWT-SVD聯(lián)合算法的沖擊波降噪方法。

1 經(jīng)驗小波變換及奇異值分解概述

1.1 經(jīng)驗小波變換(EWT)

經(jīng)驗小波變換(EWT)[11]借鑒了經(jīng)驗模態(tài)分解的自適應性又結合小波分析的理論,解決了EMD存在的模態(tài)混疊問題。經(jīng)驗小波變換的目的是把信號f(t)分解成N+1個本征模態(tài)函數(shù)fk(t)之和。本質(zhì)上是對信號頻譜進行自適應分割,從而建立合適的小波濾波器組,加小波窗后提取緊支撐傅里葉頻譜的調(diào)幅-調(diào)頻(AM-FM)成分。

1.2 奇異值分解(SVD)

奇異值分解的信號降噪方法[12-13]是先基于相空間重構理論,將待處理信號Y=[y(1),y(2),…,y(n+m-1)]構造成Hankel矩陣:

根據(jù)奇異值分解理論,對于給定的秩為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 爆炸沖擊波壓力時間模型

沖擊波超壓隨時間變化規(guī)律的數(shù)學表達式為Friedlander表達式:

式(3)中,ΔP為峰值壓力,τ+為正壓作用時間,b為衰減系數(shù)。

無限空中爆炸沖擊波超壓ΔP使用薩多夫斯基公式[14-15]計算:

馬赫反射區(qū)的超壓峰值ΔPM為

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

(6)

式(6)中,ΔP是無限空中或近地爆炸對應比例距離處的超壓峰值,α0為沖擊波入射角。

正壓持續(xù)時間為

式(7)中,r為測點到爆心的距離(m),ω為裝藥量(kg)。

2 基于EWT-SVD聯(lián)合算法的沖擊波降噪方法

2.1 EWT分解效果

沖擊波信號與振動信號不同,沖擊波信號周期性差,且存在較大的趨勢項。與EWT方法在振動信號消噪中的應用不同,針對沖擊波信號的特點,直接使用EWT對沖擊波信號進行處理,分解得到3個分量M1、M2、M3,如圖1所示。能明顯觀察到在信號端點處存在幅值很大的飛翼現(xiàn)象,若在此基礎上對信號進一步作消噪處理,將會導致超壓峰值出現(xiàn)較大誤差、上升沿時間變長,而在沖擊波數(shù)據(jù)處理中,沖擊波的超壓峰值是其關鍵的評價指標,因而不能將EWT方法直接運用在沖擊波超壓數(shù)據(jù)處理中。因而本文提出使用指數(shù)擬合的方法去除信號的趨勢項,再在此基礎上進一步作數(shù)據(jù)處理。

圖1 沖擊波信號直接進行EWT分解結果Fig.1 Shockwave signal directly performs EWT decomposition results

2.2 基于EWT-SVD聯(lián)合算法的沖擊波降噪原理

通過對沖擊波信號進行指數(shù)擬合,用原始信號減去擬合后對應點的值,避免因沖擊波信號上升沿陡峭、有效信號能量遠大于噪聲的特點而帶來的較大的飛翼現(xiàn)象;再對上述信號使用EWT算法進行分解,隨后對分解得到的分量采用SVD算法進行降噪重構,最終實現(xiàn)噪聲的濾除。

2.3 降噪方法計算步驟

步驟1 對沖擊波信號x(t)利用最小二乘指數(shù)擬合,擬合函數(shù)為Friedlander方程,得到擬合后的函數(shù)p(t)。

步驟3 對分解得到的每個MRA分量都進行相空間重構得到其對應的矩陣H,對矩陣H進行奇異值分解,通過對奇異熵增量的變化趨勢的判斷,保留前k個奇異值,其余的奇異值進行置零,然后進行逆變換,得到去噪后的矩陣H,將矩陣對應位置上的值相加取平均,得到去噪后的MRA分量。

步驟4 將所有去噪后的MRA分量及擬合所得信號p(t)相加,得到降噪后的沖擊波信號。

3 算例分析

3.1 仿真計算

3.1.1仿真模型的建立

根據(jù)經(jīng)驗公式構建ω=60 kg,r=8 m處的理想沖擊波信號P(t),并加入兩個不同頻率的模態(tài)信號從而構造出仿真信號Y(t),其表達式為

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

,

(8)

式(8)中,兩個高頻信號代表測試沖擊波信號中的兩種高頻模態(tài),并在此基礎上加入已知信噪比(SNR)為20 dB的高斯白噪聲。

對信號進行傅里葉變換得到它對應的頻譜圖,仿真信號及其頻譜圖如圖2、圖3所示,從圖中可以看出構造的仿真模型符合沖擊波信號主要集中在低頻部分,且高頻部分符合高斯白噪聲的分布的特點。

圖2 仿真模型Fig.2 Simulation model

圖3 仿真模型頻譜圖Fig.3 Simulation model spectrogram

3.1.2EWT-SVD分解效果

使用本文提出的方法對仿真信號進行處理,待處理的沖擊波超壓信號是由離散點集構成的,曲線擬合的目標就是尋找這些點所對應的時間與其具體值之間的函數(shù)關系P=f(t)。在擬合過程中只對沖擊波到達后的衰減過程進行擬合,得到的系數(shù)如表1所示。

表1 指數(shù)擬合所得Friedlander方程參數(shù)Tab.1 Parameters of the Friedlander equation obtained by exponential fitting

擬合所得數(shù)值與理論公式計算系數(shù)差距極小,很好地表示了沖擊波信號的變化趨勢。對減去擬合后的數(shù)據(jù)進行EWT處理,其分解結果及頻譜圖如圖4、圖5所示。

圖4 EWT分解結果Fig.4 EWT decomposition results

圖5 EWT分解結果頻譜圖Fig.5 EWT decomposition results in a spectrogram

該沖擊波信號分解出了兩個模態(tài),分別對應在仿真信號中添加的兩個高頻信號。與圖1對比可以看出飛翼現(xiàn)象得到明顯抑制。對分解得到的模態(tài)進行相空間的重構,對于數(shù)據(jù)長度為N的信號重構矩陣時選取m=N/2為重構矩陣的行數(shù),根據(jù)奇異熵增量選取重構奇異值的階數(shù),得到降噪后的模態(tài),其頻譜如圖6所示。

圖6 降噪后的分量頻譜圖Fig.6 Component spectrogram after noise reduction

將降噪后的模態(tài)與擬合數(shù)據(jù)進行相加得到最終降噪后的數(shù)據(jù),并且有效地分離出加入的兩個高頻模態(tài)。作為比較使用EEMD方法對仿真信號進行分解,得到13個分量,與EEMD方法相比基于EWT-SVD聯(lián)合算法的沖擊波降噪方法很好地分離出加入的模態(tài)。對爆炸場沖擊波中包含的不同種類的噪聲的特性認識及區(qū)分具有一定的指導意義。

從去噪的能力方面對本文方法進行評估,已知對信號添加了信噪比為20.073 6 dB的高斯白噪聲,對降噪后的信號求取信噪比其值為19.913 4 dB,證明了該方法在降噪方面的有效性。

3.2 試驗驗證

選取某次裝藥量60 kg,裝藥高度1.5 m,測試點距爆心8 m的試驗數(shù)據(jù),對該信號使用本文方法進行處理,首先對信號進行擬合,得到擬合后的Friedlander方程參數(shù)ΔP為0.351 8 MPa,τ+為0.004 s,b為1.113。

作為對比,對實測信號直接進行EWT分解,其結果如圖7所示。對減去擬合數(shù)據(jù)后的實測信號進行EWT分解,所得分量如圖8所示。

可以看出各個分量的飛翼現(xiàn)象得到了明顯抑制。將所得的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聯(lián)合算法的沖擊波降噪方法很好地還原出了沖擊波信號的超壓峰值、上升時間以及變化趨勢且濾除了試驗信號中的噪聲干擾。從6階貝塞爾濾波在超壓峰值處可以觀察到峰值大小的改變,且整個信號存在相位的偏移。該方法在提取超壓峰值、抑制相位偏移、還原上升沿時間方面有著更好的優(yōu)勢。

圖12 聯(lián)合降噪后的結果Fig.12 Results after joint noise reduction

4 結論

本文提出了基于EWT-SVD聯(lián)合算法的沖擊波降噪方法。該方法首先對沖擊波信號進行指數(shù)擬合,用原始信號減去擬合后對應點的值;對上述信號使用EWT算法進行分解,隨后對分解得到的分量采用SVD算法進行降噪重構,實現(xiàn)噪聲的濾除。仿真及試驗驗證表明,該方法可有效抑制相位偏移,更好地還原出沖擊波信號的超壓峰值、上升沿時間,可用于沖擊波信號的降噪。該方法不足之處為對相空間重構后的矩陣進行SVD算法計算量過大,尚需改進。