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

(1.中北大學(xué)省部共建動(dòng)態(tài)測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,山西 太原 030051;2.中北大學(xué)電氣與控制工程學(xué)院,山西 太原 030051)

0 引言

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

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

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

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

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

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

1.2 奇異值分解(SVD)

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

根據(jù)奇異值分解理論,對(duì)于給定的秩為r的m×n維矩陣H,則H的奇異值分解形式為

式(2)中,U,V分別為m×m,n×n維酉矩陣;Σ為r×r維對(duì)角陣,其對(duì)角線元素為矩陣H的非零奇異值σi,且以降序排列,即σ1≥σ2≥…≥σr。若源信號(hào)Y是由有用信號(hào)和噪聲共同組成,則矩陣H的奇異值σi可以反映信號(hào)和噪聲能量集中的情況。前p個(gè)較大的奇異值將主要反映有用信號(hào),較小的奇異值則主要反映噪聲,把這部分反映噪聲的奇異值置零就可以去除信號(hào)中的噪聲。將重構(gòu)后的矩陣對(duì)應(yīng)的項(xiàng)相加取平均值,就可以還原出降噪后的信號(hào)。

1.3 爆炸沖擊波壓力時(shí)間模型

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

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

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

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

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

(6)

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

正壓持續(xù)時(shí)間為

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

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

2.1 EWT分解效果

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

圖1 沖擊波信號(hào)直接進(jìn)行EWT分解結(jié)果Fig.1 Shockwave signal directly performs EWT decomposition results

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

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

2.3 降噪方法計(jì)算步驟

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

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

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

3 算例分析

3.1 仿真計(jì)算

3.1.1仿真模型的建立

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

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

,

(8)

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

對(duì)信號(hào)進(jìn)行傅里葉變換得到它對(duì)應(yīng)的頻譜圖,仿真信號(hào)及其頻譜圖如圖2、圖3所示,從圖中可以看出構(gòu)造的仿真模型符合沖擊波信號(hào)主要集中在低頻部分,且高頻部分符合高斯白噪聲的分布的特點(diǎn)。

圖2 仿真模型Fig.2 Simulation model

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

3.1.2EWT-SVD分解效果

使用本文提出的方法對(duì)仿真信號(hào)進(jìn)行處理,待處理的沖擊波超壓信號(hào)是由離散點(diǎn)集構(gòu)成的,曲線擬合的目標(biāo)就是尋找這些點(diǎn)所對(duì)應(yīng)的時(shí)間與其具體值之間的函數(shù)關(guān)系P=f(t)。在擬合過程中只對(duì)沖擊波到達(dá)后的衰減過程進(jìn)行擬合,得到的系數(shù)如表1所示。

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

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

圖4 EWT分解結(jié)果Fig.4 EWT decomposition results

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

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

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

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

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

3.2 試驗(yàn)驗(yàn)證

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

作為對(duì)比,對(duì)實(shí)測(cè)信號(hào)直接進(jìn)行EWT分解,其結(jié)果如圖7所示。對(duì)減去擬合數(shù)據(jù)后的實(shí)測(cè)信號(hào)進(jìn)行EWT分解,所得分量如圖8所示。

可以看出各個(gè)分量的飛翼現(xiàn)象得到了明顯抑制。將所得的3個(gè)分量分別進(jìn)行相空間重構(gòu),對(duì)重構(gòu)后的矩陣進(jìn)行奇異值分解,計(jì)算相應(yīng)的奇異值熵增量。以分量2為例,所對(duì)應(yīng)的奇異值熵增量如圖9所示。

可以看出24階之后的奇異值熵增量小于0.025,將24階之后的奇異值置零,然后進(jìn)行分量的重構(gòu)。SVD降噪前的分量頻譜如圖10所示。重構(gòu)后的3個(gè)分量的頻譜如圖11所示。

圖7 直接EWT分解結(jié)果Fig.7 Direct EWT decomposition results

圖8 減去擬合后的信號(hào)EWT分解結(jié)果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

通過對(duì)比,明顯觀察到噪聲得到了極大地濾除。將分量及擬合信號(hào)相加,得到最終降噪后的沖擊波信號(hào)。作為比較對(duì)原始信號(hào)進(jìn)行6階貝塞爾濾波,截止頻率100 kHz。原始信號(hào)和兩種方法濾波后的結(jié)果如圖12所示?；贓WT-SVD聯(lián)合算法的沖擊波降噪方法很好地還原出了沖擊波信號(hào)的超壓峰值、上升時(shí)間以及變化趨勢(shì)且濾除了試驗(yàn)信號(hào)中的噪聲干擾。從6階貝塞爾濾波在超壓峰值處可以觀察到峰值大小的改變,且整個(gè)信號(hào)存在相位的偏移。該方法在提取超壓峰值、抑制相位偏移、還原上升沿時(shí)間方面有著更好的優(yōu)勢(shì)。

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

4 結(jié)論

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