去除爆炸沖擊波信號高頻噪聲的聯(lián)合處理方法

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

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

0 引言

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

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

1.1 HHT(Hilbert-Huang)

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

1.1.1EMD

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

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

2) 求上下包絡線的均值曲線m1(t)

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

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

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

(2)

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

(3)

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

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表示本征模態(tài)分量分解次數(shù)。

余量r(t)同理表示為

(11)

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

(12)

1.2 Bessel低通濾波器

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

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

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

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

圖1 HHT與Bessel聯(lián)合方法流程圖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為裝藥量的質(zhì)量(單位:kg)。

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

式(14)中,H為裝藥的高度(單位:m),w為裝藥量的質(zhì)量(單位: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聯(lián)合的信號處理過程

3.2.1EMD

原始沖擊波信號被分解成9個IMF和1個余量,并按時間尺度從高頻到低頻依次分解,圖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和余量與原始沖擊波信號的互相關(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

對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重構(gòu)結(jié)果分析

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

圖11 處理前后信號曲線對比圖Fig.11 Signal curve comparison before and after processing

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

3.3 濾波方法對比分析

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

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

SNR和RMSE計算公式如下:

(14)

(15)

式中,X表示處理前沖擊波信號,x表示處理后沖擊波信號,N為信號長度。

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

表2 濾波指標對比Tab.2 Comparison of filtering indexes

4 結(jié)論

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

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