劉 宇,李新娥,崔春生,李 順

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

0 引言

爆炸沖擊波是典型的非線(xiàn)性、非平穩(wěn)信號(hào),上升沿陡峭且上升時(shí)間短,由于測(cè)試環(huán)境惡劣,測(cè)試信號(hào)混有大量高頻干擾信號(hào)。沖擊波的能量主要集中在低頻段,干擾噪聲分布在高頻段,其能量幅值較小,頻帶較寬。因此為了得到爆炸沖擊波的本征特點(diǎn),有必要對(duì)采集到的爆炸沖擊波信號(hào)進(jìn)行去噪。目前,國(guó)內(nèi)外爆炸沖擊波信號(hào)的數(shù)據(jù)處理方法主要有Bessel[1-3]、小波閾值[1-3]、Butterworth[1-4]、經(jīng)驗(yàn)?zāi)B(tài)分解[4-7](empirical mode decomposition,EMD)、希爾伯特黃變換[8-9](Hilbert-Huang,HHT)、EMD與小波閾值聯(lián)合[10]。上述信號(hào)處理方法均各有利弊,其中,Bessel和Butterworth為經(jīng)典濾波器,設(shè)置階數(shù)和截止頻率,進(jìn)行低通濾波,去除高頻信號(hào),設(shè)計(jì)簡(jiǎn)單。Bessel針對(duì)整個(gè)信號(hào),存在去除混疊的有用信號(hào)的缺點(diǎn)。Butterworth濾波結(jié)果會(huì)出現(xiàn)較大的相位滯后,不符合沖擊波上升時(shí)間短的特點(diǎn)。EMD根據(jù)時(shí)間尺度特征將原始信號(hào)分解成為一系列IMF分量,能夠反映出非平穩(wěn)信號(hào)的局部特征。但也存在模態(tài)混疊、包絡(luò)線(xiàn)擬合偏差、端點(diǎn)效應(yīng)等缺點(diǎn)。直接剔除含有高頻噪聲的IMF,易丟失有用信號(hào),造成失真。小波閾值、EMD與小波閾值聯(lián)合都存在小波基、分解層數(shù)、閾值函數(shù)難以確定的問(wèn)題。針對(duì)上述情況,提出了HHT與Bessel聯(lián)合的信號(hào)處理方法,無(wú)需預(yù)先設(shè)定任何基函數(shù),避免了經(jīng)典低通濾波器針對(duì)整個(gè)信號(hào),易去除有用信號(hào)的缺點(diǎn),結(jié)合了Bessel設(shè)計(jì)簡(jiǎn)單和HHT自適應(yīng)性好且適用于分析非線(xiàn)性非平穩(wěn)信號(hào)的優(yōu)點(diǎn),該方法在剔除噪聲的同時(shí),還能較好地保留沖擊波信號(hào)的本征特點(diǎn)。

1 HHT和Bessel基礎(chǔ)理論

1.1 HHT(Hilbert-Huang)

HHT包含EMD和HT(Hilbert transform,希爾伯特變換)兩部分,適用于分析非線(xiàn)性非平穩(wěn)信號(hào)。

1.1.1EMD

EMD是時(shí)頻域的處理方法,將原始信號(hào)分解成為一系列IMF分量,IMF分量是具有時(shí)變頻率的振蕩函數(shù),能夠反映出非平穩(wěn)信號(hào)的局部特征。利用EMD對(duì)原始沖擊波信號(hào)X(t)進(jìn)行分解,步驟如下:

1) 求X(t)極大值和極小值,擬合上包絡(luò)線(xiàn)u1(t)、下包絡(luò)線(xiàn)v1(t)。

2) 求上下包絡(luò)線(xiàn)的均值曲線(xiàn)m1(t)

3) 求中間曲線(xiàn)h1(t),h1(t)滿(mǎn)足下面兩個(gè)條件即為本征模態(tài)函數(shù)IMF1,進(jìn)行下一步,若不滿(mǎn)足,重復(fù)步驟1)—3)。

①所有極值點(diǎn)和零點(diǎn)的數(shù)量相同或相差1;②上下包絡(luò)線(xiàn)的局部均值為0。

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

(2)

4) 此時(shí)均值曲線(xiàn)m1(t)為IMF1的殘余量,對(duì)m1(t)進(jìn)行步驟1)—3)分解步驟,得到IMF2,以此類(lèi)推,直到分解為單調(diào)信號(hào)為止,分解得到所有IMF和余量,最終的剩余分量為r(t)。

(3)

式(3)中,n為分量個(gè)數(shù)。

1.1.2HT(希爾伯特變換)

HT是分析和處理信號(hào)的一個(gè)重要的理論工具,可以提供90°的相位變換而不影響頻率分量的幅度。IMF和余量的瞬時(shí)頻率有著明確的物理意義。因此,EMD后, 對(duì)IMF和余量作希爾伯特變換(HT),求得IMF和余量的瞬時(shí)頻率和瞬時(shí)幅值。

希爾伯特變換將s(t)與h(t)作卷積。

將IMFi(t)作希爾伯特變換,得式(6),繼而求得每個(gè)IMF的瞬時(shí)頻率wi(t)和瞬時(shí)幅值ai(t),如式(8)—式(9)。

所有IMF求和的時(shí)間-頻率-能量三維頻譜表示為

式(10)中,i表示本征模態(tài)分量分解次數(shù)。

余量r(t)同理表示為

(11)

每一個(gè)IMF和余量重構(gòu),將原始信號(hào)表示為

(12)

1.2 Bessel低通濾波器

Bessel濾波器具有平時(shí)延特性,較平坦的幅度特性,可以減小固有的非線(xiàn)性相位失真。在信號(hào)阻帶衰減方面存在劣勢(shì),一般設(shè)計(jì)為高階的濾波器來(lái)達(dá)到阻帶衰減水平。

2 HHT與Bessel聯(lián)合方法

HHT與Bessel聯(lián)合的信號(hào)處理方法,結(jié)合了Bessel設(shè)計(jì)簡(jiǎn)單和HHT自適應(yīng)性好且適用于分析非線(xiàn)性非平穩(wěn)信號(hào)的優(yōu)點(diǎn)。

利用 EMD將信號(hào)分解成一系列IMF和一個(gè)余量。首先對(duì)各IMF和余量與原始信號(hào)的互相關(guān)系數(shù)進(jìn)行分析,初步判斷含高頻噪聲的IMF,并采用Hilbert頻譜分析最終確定含高頻噪聲的IMF分量;然后采用Bessel僅對(duì)含高頻噪聲的IMF處理;最后將處理后的IMF、純凈的IMF和余量進(jìn)行重構(gòu),從而得到最終處理后的信號(hào)。流程圖如圖1所示。

圖1 HHT與Bessel聯(lián)合方法流程圖Fig.1 Joint flow chart of HHT and Bessel

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

3.1 試驗(yàn)方案

自由場(chǎng)爆炸沖擊波測(cè)試過(guò)程中,針對(duì)測(cè)試點(diǎn)距離爆點(diǎn)近,采集的沖擊波信號(hào)信噪比差的問(wèn)題,有必要對(duì)距離爆心較近的測(cè)點(diǎn)處采集到的信號(hào)進(jìn)行去噪處理。60 kg球形TNT,比例距離為 1.02 m/kg1/3,裝藥的高度距地面 1.5 m,三個(gè)測(cè)點(diǎn)距爆心距離相等。測(cè)點(diǎn)分布如圖2所示。

圖2 測(cè)點(diǎn)分布示意圖Fig.2 Distribution diagram of measuring points

比例距離R的計(jì)算公式為

式(13)中,R為比例距離(單位:m/kg1/3),r為測(cè)試點(diǎn)離爆心的距離(單位:m),w為裝藥量的質(zhì)量(單位:kg)。

炸藥在無(wú)邊界的空中爆炸,裝藥的高度H符合

式(14)中,H為裝藥的高度(單位:m),w為裝藥量的質(zhì)量(單位:kg)。

采集的沖擊波信號(hào)如圖3所示,測(cè)點(diǎn)1、測(cè)點(diǎn)2、測(cè)點(diǎn)3的超壓峰值分別為2.431、2.310、2.430 MPa。以測(cè)點(diǎn)1為例,詳細(xì)分析信號(hào)處理過(guò)程,測(cè)點(diǎn)2和測(cè)點(diǎn)3的分析過(guò)程不再贅述。

圖3 原始沖擊波信號(hào)曲線(xiàn)Fig.3 Original shock wave signal curve

3.2 HHT與Bessel聯(lián)合的信號(hào)處理過(guò)程

3.2.1EMD

原始沖擊波信號(hào)被分解成9個(gè)IMF和1個(gè)余量,并按時(shí)間尺度從高頻到低頻依次分解,圖4表示IMF1—IMF5,圖5表示IMF6—IMF9和余量,重構(gòu)誤差為5.144 0×10-17。

圖4 IMF1—IMF5Fig.4 IMF1—IMF5

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

3.2.2互相關(guān)系數(shù)分析

IMF1—IMF9和余量與原始沖擊波信號(hào)的互相關(guān)系數(shù)如表1所示,IMF1—IMF3互相關(guān)系數(shù)小于0.1。初步判定IMF1—IMF3是含噪IMF。

表1 互相關(guān)系數(shù)表Tab.1 Table of the cross-correlation coefficient

3.2.3確定含噪IMF

對(duì)IMF1—IMF3分別作Hilbert頻譜分析,如圖6—圖8所示。通過(guò)三維時(shí)間-頻率-能量譜分析頻率在某段具體時(shí)間的能量譜密度,確定真實(shí)含噪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以不同的頻率顯示原始信號(hào)特征。IMF1主要頻率范圍在0～400 kHz, IMF2主要頻率范圍在0～200 kHz,IMF3主要頻率范圍在0～70 kHz。確定真實(shí)含高頻噪聲的IMF為IMF1—IMF2。

3.2.4對(duì)含噪IMF進(jìn)行Bessel處理

對(duì)含噪IMF進(jìn)行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重構(gòu)結(jié)果分析

將處理后的IMF、純凈的IMF和余量進(jìn)行重構(gòu),得到重構(gòu)信號(hào)。與原始沖擊波信號(hào)進(jìn)行對(duì)比,處理前后信號(hào)曲線(xiàn)對(duì)比如圖11所示,重構(gòu)信號(hào)與原始沖擊波信號(hào)趨勢(shì)較為貼近。對(duì)重構(gòu)信號(hào)進(jìn)行Hilbert譜分析,如圖12所示,可見(jiàn)去除了高頻噪聲干擾,保留了沖擊波信號(hào)的本征特點(diǎn)。

圖11 處理前后信號(hào)曲線(xiàn)對(duì)比圖Fig.11 Signal curve comparison before and after processing

圖12 重構(gòu)信號(hào)的Hilbert譜Fig.12 Hilbert spectrum of reconstructed signal

3.3 濾波方法對(duì)比分析

對(duì)EMD、Bessel、小波閾值、EMD與小波閾值聯(lián)合方法以及提出的HHT與Bessel聯(lián)合方法進(jìn)行對(duì)比分析。

SNR(信噪比)、RMSE(均方根誤差)和超壓峰值是信號(hào)處理的評(píng)價(jià)指標(biāo),SNR反映信號(hào)和噪聲的比值, RMSE對(duì)信號(hào)處理前后的特大或特小誤差尤為敏感。SNR越大,RMSE越小,超壓峰值越貼近原始數(shù)據(jù)超壓峰值,去噪精度越高。

SNR和RMSE計(jì)算公式如下:

(14)

(15)

式中,X表示處理前沖擊波信號(hào),x表示處理后沖擊波信號(hào),N為信號(hào)長(zhǎng)度。

通過(guò)以上濾波方法分別對(duì)測(cè)點(diǎn)1、測(cè)點(diǎn)2、測(cè)點(diǎn)3的原始信號(hào)進(jìn)行處理,結(jié)果如圖13—圖15所示,濾波指標(biāo)如表2所示。

圖13 測(cè)點(diǎn)1濾波方法對(duì)比Fig.13 Comparison of filtering methods of measuring point 1

圖15 測(cè)點(diǎn)3濾波方法對(duì)比Fig.15 Comparison of filtering methods of measuring point 3

由表2分析可知:由于Bessel針對(duì)整個(gè)信號(hào)設(shè)置階數(shù)和截止頻率,易去除有用信號(hào),超壓峰值相比其他幾種方法較低;小波閾值自適應(yīng)差,預(yù)先設(shè)置不同的小波基、分解層數(shù)、閾值函數(shù)等參數(shù),濾波效果不同,需要人為選擇搭配參數(shù),得到最優(yōu)效果;EMD分解后,判斷含高頻噪聲的IMF,直接去除,將其余純凈IMF重構(gòu),未篩選其中有用信息,去噪精度低;EMD與小波閾值聯(lián)合提高了精度,EMD分解后,判斷含高頻噪聲的IMF分量,利用小波閾值對(duì)含高頻噪聲的IMF分量處理,最后重構(gòu),得到處理后的信號(hào),但同時(shí)存在小波基等參數(shù)的選擇不當(dāng),造成去噪效果較差的結(jié)果。相較于其他濾波方法,HHT與Bessel聯(lián)合采用EMD分解,通過(guò)互相關(guān)系數(shù)和Hilbert譜分析確定含噪IMF,采用Bessel低通數(shù)字濾波器對(duì)含噪IMF進(jìn)行處理,最后重構(gòu),得到處理后的信號(hào),在以上3個(gè)測(cè)點(diǎn)的信號(hào)處理中都占據(jù)優(yōu)勢(shì), SNR更大,RMSE更小,處理后超壓峰值接近原始信號(hào)的超壓峰值。

表2 濾波指標(biāo)對(duì)比Tab.2 Comparison of filtering indexes

4 結(jié)論

HHT與Bessel聯(lián)合的沖擊波信號(hào)處理方法,結(jié)合了Bessel設(shè)計(jì)簡(jiǎn)單和HHT自適應(yīng)性好且適用于分析非線(xiàn)性非平穩(wěn)信號(hào)的優(yōu)點(diǎn)。采用EMD,通過(guò)互相關(guān)系數(shù)和Hilbert譜分析確定含噪IMF,并利用Bessel低通數(shù)字濾波器對(duì)含噪IMF進(jìn)行處理,最后重構(gòu),得到處理后的沖擊波信號(hào)。

相比于目前已有的EMD、Bessel、小波閾值、EMD與小波閾值聯(lián)合,結(jié)合SNR、RMSE、超壓峰值各項(xiàng)濾波指標(biāo)綜合分析,該方法SNR更大,RMSE更小,超壓峰值接近原始信號(hào)的超壓峰值。在剔除噪聲的同時(shí),還能較好地保留信號(hào)中的本征特點(diǎn),提高了爆炸沖擊波信號(hào)分析的準(zhǔn)確性,去噪精度更高。