基于經驗模態分解及小波變換的炸藥NQR信號處理

郝鳳龍，徐更光，黃學義

(北京理工大學 爆炸科學與技術國家重點實驗室，北京 100081)

針對恐怖爆炸襲擊，及時有效查出隱藏的爆炸物已成當務之急。核電四極矩共振(Nuclear Quadrupole Resonance, NQR)探測技術以準確率高、誤報率低、無磁污染等優點深受青睞[1-3]。然而，炸藥的NQR頻率較低(0.5～6 MHz)，信號強度極弱，并易受背景噪聲干擾，信噪比較低，探測識別較困難。如何有效提高NQR信號的信噪比，已成為基于NQR技術的炸藥探測系統研制重點。原子核四極矩共振為原子核物理現象，指原子核非球對稱部分因與核外電場梯度相互作用引起能級分裂，在外加射頻場作用下，產生能級躍遷過程[4]。炸藥中普遍含自旋量子數I為1的14N，其原子核具備產生核電四極矩共振的內在因素，可將14N作為炸藥探測的一種特征成分[5]；而不同原子核、同種原子核在不同物質中或同種物質在不同晶型下，NQR頻率君不同，一旦檢測到14N原子核的NQR信號，便可據其特征頻率譜線唯一確定炸藥類型，實現對炸藥的探測與識別。

炸藥產生的NQR信號非常微弱，僅為納伏量級，極易被淹沒在背景噪聲中，極低的信噪比制約了NQR技術在炸藥探測中的應用；因此，須采取有效信號處理方法提高信噪比，進而提高炸藥探測準確率。文獻[6]利用加權傅里葉頻譜分析方法對NQR信號進行處理。該方法充分利用有限信號能量，能降低噪聲能量影響，因而信噪比較高。然而，該方法需信號滿足平穩性假設條件，無法描述信號局部特征。文獻[7]提出自適應濾波算法，利用干擾信號間相關性，有效抑制振鈴拖尾，提取出干擾中的NQR信號，從而提高NQR信號的信噪比。文獻[8]用小波變換方法提高信噪比，利用Mallat算法對NQR信號進行分解及重構，并通過多個閾值函數對小波系數作門限閾值處理，使噪聲得到有效抑制。此方法缺陷為小波基與閾值函數選取較困難，需繁瑣調試才能達到好效果，缺乏自適應能力。文獻[9]提出基于改進閾值函數小波變換的NQR信號處理方法，將軟硬閾值函數結合，動態調節閾值大小，既兼具軟硬閾值法優點又能避免其缺陷，可有效濾除背景噪聲，提高信噪比。

本文通過分析當前炸藥NQR信號處理方法優缺點，提出基于經驗模態分解(Empirical Mode Decomposition, EMD)與小波變換聯合的NQR信號去噪方法。該方法不受傅里葉變換及小波函數選擇限制，自適應能力好，并通過對實測炸藥NQR信號去噪處理，驗證該方法的可行性及有效性。

1 經驗模態分解去噪原理

經驗模態分解方法[10]為針對非線性及非平穩信號的時頻域信號處理方法，為Hilbert-Huang變換 (HHT) 的核心部分。該方法由信號本身時間尺度特征出發，將含噪聲的原始信號分解成一系列保留局部特征信息的數據序列，即固有模態函數 (Intrinsic Mode Function, IMF)，再對該序列進行處理、信號重構，可有效去除信號中混雜噪聲。固有模態函數須滿足兩條件：① 在整個數據序列中極值點數量與過零點數量須相等或最多相差一個；② 在任意時間點由數據序列局部極大值點確定的上包絡線與局部極小值點確定的下包絡線均值為零，即信號關于時間軸局部對稱[11]。EMD方法將信號分解處理過程描述為篩選過程，對給定的時間序列信號x(t)其分解步驟如下：

(1) 確定x(t)所有極大值點、極小值點；

(2) 在極值點間利用三次樣條函數進行插值，獲得x(t)上包絡線u0(t)與下包絡線v0(t)，并計算兩包絡線均值m0(t)，即

(1)

(3) 用x(t)減去平均曲線m0(t)，得差值h1(t)為

h1(t)=x(t)-m0(t)

(2)

判斷差值h1(t)是否滿足IMF的兩條件，若滿足則令h1(t)為x(t)的第一階IMF分量c1(t)，即c1(t)=h1(t)；若不滿足則用h1(t)代替x(t)，重復以上步驟，直到滿足IMF條件為止，求得第一階IMF分量c1(t)；

(4) 求出原信號與第一階IMF分量差值,即

r1(t)=x(t)-c1(t)

(3)

將剩余部分r1(t)作為新信號x(t)處理，重復步驟(1)～(4)，可依次獲得c2(t)，c3(t)，…cn(t)，直到剩余部分rn(t)為單調函數時分解完成。

通過EMD 方法篩選處理，原始信號x(t)可表示為所有IMF分量及殘余量之和，即

(4)

式中：n為IMF分量個數；rn(t)為信號x(t)平均趨勢。

頻域上分量IMF已含原始信號中多個頻率段成分，且據計算次序所得各階IMF分量所含頻率成分由高到低分布，即階數小的IMF代表信號高頻成分，階數大的IMF代表信號低頻成分。基于EMD的去噪方法主要思路為：大多數被噪聲干擾的目標信號主要能量集中在低頻段，頻段越高所含信號能量越少，而噪聲主要集中在高頻段，因此可舍棄部分階數小的IMF，利用剩余IMF重構信號，達到去噪目的。由分析可知，EMD過程實為頻帶篩分過程，完全自適應于被分解信號，無需預先提供分解的基函數，自適應性更好、靈活性更高。該方法也存在不足，如較難確定舍棄的IMF數目，舍棄部分階數小IMF的同時也會丟失信號部分能量。

2 小波閾值去噪原理

小波變換因具有良好的時頻局部化特性及多分辨分析特性得以廣泛應用。Donoho等[12]在小波變換基礎上提出基于小波閾值的信號去噪方法，基本原理即據信號、噪聲在各層小波空間分別具有的不同特性實現信噪分離。信號主要特征由分布在較大尺度上少數幅值較大系數表征；而噪聲主要特征由分布在各層小波空間多個小幅值小波系數表征。將信號據所選小波基函數進行分解，獲得小波域上小波系數。利用閾值函數將各尺度上由噪聲產生的小波分量濾除，重構原始信號完成去噪。

設含噪信號y(t)=s(t)+n(t)。其中s(t)為原始信號；n(t)為高斯白噪聲，服從N(0,σ2)分布。對y(t)作離散小波變換為

wy(j,k)=ws(j,k)+wn(j,k)

(5)

式中：j=0,1,2…,J，J為小波變換最大分解層數；k=0,1,2…,N，N為信號長度；wy(j,k)，ws(j,k)，wn(j,k)為含噪信號、原始信號及噪聲信號在小波空間第j層的小波系數。wy(j,k)小于某一閾值時，wy(j,k)主要由噪聲引起，可將其舍去；wy(j,k)大于閾值時, 小波系數主要由信號引起，需保留。

傳統上對小波系數作門限閾值處理有硬閾值函數法與軟閾值函數法兩種[13],分別定義為

(6)

(7)

式(6)為硬閾值法, 式(7)為軟閾值法。閾值確定規則包括無偏風險規則 (Rigrsure)、固定閾值規則 (Sqtwolog)、啟發式閾值規則 (Heursure) 及極大極小值規則 (Minimaxi)。閾值方法本身存在一定缺陷，其中硬閾值法去噪函數在λ與-λ兩點處存在間斷點，所得小波系數值連續性較差，會導致重構信號振蕩。軟閾值函數雖在小波域連續，但軟閾值對大于閾值的小波系數采取恒定值壓縮，小波系數較大時會給重構信號帶來誤差；因此，閾值估計及閾值函數選取成為閾值去噪方法的關鍵與難點。

3 經驗模態分解與小波變換聯合去噪原理

在分析EMD去噪與小波閾值去噪優缺點基礎上，本文提出將二者聯合的NQR信號去噪方法，以充分發揮其各自優點。經EMD篩選的IMF分量恰好滿足由高頻到低頻的系列分布，低頻段IMF分量由信號主導，高頻段IMF分量由噪聲主導；因此，必存在某個IMFjs分量，其后的IMF分量由信號起主導作用，而之前js個IMF分量由噪聲起主導作用。此處通過連續均方誤差準則確定分界點js。定義重構信號為

(8)

式中：IMFj(t)為EMD分解所得第j個IMF分量；rn(t)為殘余量；n為原始信號分解后所得IMF分量個數。信號的連續均方誤差[14]定義為

(9)

其中：N為信號長度；IMFk(ti)為信號分解所得第k個分量。基于以上分析，信號能量分界點可確定為

(10)

確定信號能量分界點后，對含噪的高頻分量用小波閾值去噪處理，將去噪數據與不含噪聲的低頻IMF分量及殘余量重構原始信號，不損失存于高頻IMF分量中的有用信息，且小波閾值去噪僅作用于高頻IMF分量非直接作用于整個信號，較大程度上能克服小波閾值去噪缺陷。基于EMD與小波閾值聯合的NQR信號去噪處理流程見圖1。

圖1 基于EMD與小波閾值的信號去噪流程

4 炸藥NQR信號去噪分析

4.1 NQR信號模型

NQR信號有兩種類型：① 自由感應衰減(Free Induction Decay, FID)信號，在射頻脈沖結束后立刻呈指數形式衰減；② 自旋回波(Spin-Echo，SE)信號，為FID 信號再重聚結果。實測中NQR信號強度較弱，易受線圈內熱噪聲、發射機射頻及外部電磁波等干擾；通過接收探頭檢測的NQR信號可視為原始FID信號即指數衰減正弦信號、噪聲信號及隨機干擾信號之和，表示成復值時間序列[15]為

(11)

式中：n=0,1,…,N-1；αk，βk分別為第k條正弦曲線振幅及衰減常數，炸藥不同衰減常數亦不同；d為衰減正弦曲線分量個數；ωk(τ)為第k條正弦曲線頻率漂移函數，為炸藥樣品溫度τ的線性函數；ω(n)為加性噪聲，可設為零均值復高斯白噪聲。

4.2 黑索金NQR信號去噪實驗

為驗證本文方法對炸藥NQR信號的去噪效果，由實驗測試獲得黑索金炸藥NQR信號，對實測信號分別進行直接小波閾值去噪、直接EMD去噪及EMD與小波聯合去噪處理，分析比較各種方法的去噪性能。

在室溫條件下將20 g 黑索金粉末樣品密封于螺口玻璃瓶中置于封閉射頻線圈內，用PSL組合脈沖序列，激勵頻率3.41 MHz，采樣點數1 000，接收增益40 dB，累加次數100。探測的時域NQR信號見圖2，可以看出NQR信號極微弱，幾乎被淹沒在背景噪聲中。

圖2 原始黑索金NQR信號

選db3小波對原始NQR信號進行分解，分解層數為3，采用軟閾值去噪方法，閾值據自適應無偏風險規則確定，處理后NQR信號見圖3，NQR信號已從背景噪聲中分離，較光滑，但在整個時域上有較多毛刺。

圖3 直接用小波軟閾值法所得去噪效果

對圖2 NQR信號進行EMD分解，所得分量見圖4，含IMF1～IMF99個固有模態函數與殘余分量RES。由圖4看出，噪聲主要分布在IMF1，IMF2分量中，用直接EMD去噪方法，將IMF1，IMF2分量舍棄，所得NQR信號見圖5，因兩高頻分量中有用信號被刪除, 致信號強度下降；因此該去噪方法較粗糙。

圖4 固有模態函數分量及殘余量

圖5 舍棄IMF1，IMF2分量所得去噪效果

用db3小波分別對IMF1，IMF2分量進行軟閾值去噪處理，所得剩余信號見圖6、圖7。由兩圖看出，兩高頻分量中確含有用信號特征信息。將去噪后IMF1，IMF2分量與低頻分量重構信號，所得NQR信號見圖8，可見信號較光滑，幾乎無毛刺。

對三種去噪方法得到的時域NQR信號進行時頻域轉換，得到的頻域信號如圖9、10、11所示，通過分析對比發現，圖9的噪聲強度最大，圖10的噪聲強度有所下降，但NQR信號強度最小，而圖11的噪聲基本可以忽略，且NQR信號幅值最大，在頻域中較好地保留了有用信號的特征信息，克服了單獨小波軟閾值去噪和EMD去噪的缺陷。

采用信噪比評估各種方法的去噪性能，計算結果見表1。由表1看出，本文方法抑制噪聲更有效，炸藥NQR信號信噪比得以提高，去噪效果較好，從而驗證本文方法的有效性及可行性。

表1 小波、EMD及本文方法去噪性能比較

圖6 IMF1小波軟閾值去噪效果

圖9 小波軟閾值去噪后NQR信號頻譜

5 結 論

本文基于NQR信號非線性與非平穩性特點，提出基于經驗模態分解(EMD)與小波變換聯合的炸藥NQR信號去噪方法。通過對探測的黑索金NQR信號進行經驗模態分解，分析各階IMF分量特征信息，對噪聲起主導作用的高頻IMF分量進行小波軟閾值去噪，將去噪后高頻IMF分量與低頻IMF分量重構信號，既可保留有用信號特征信息又能有效抑制噪聲、提高信噪比，很大程度上克服小波閾值去噪及直接EMD去噪的缺陷。該方法在炸藥NQR信號去噪中的良好性能，可為NQR信號處理奠定理論、技術基礎。

[1] Mozzhukhin G V. Three-frequency composite multipulse nuclear quadrupole resonance technique for explosive detection[J]. Applied Magnetic Resonance, 2012, 43(4): 547-556.

[2] Rati R, Pink R H, Scheicher R H, et al. Nuclear quadrupole interactions in nuclear quadrupole resonance detection of energetic and controlled materials: theoretical study[J]. Applied Magnetic Resonance, 2012, 43(4): 591-617.

[3] 李康寧，俞碩，李興，等. 核四極共振技術在黑火藥探測中的研究[J]. 核科學與工程, 2011, 31 (3): 270-273.

LI Kang-ning, YU Shuo, LI Xing, et al. Research in the detection of black powder based on nuclear quadrupole resonance technology[J]. Chinese Journal of Nuclear Science and Engineering, 2011, 31 (3): 270-273.

[4] Peshkovsky A S, Cattena C J, Cerioni L M, et al. Noise-resilient multi-frequency surface sensor for nuclear quadrupole resonance[J]. Journal of Magnetic Resonance, 2008, 194(2): 222-229.

[5] Mozzhukhin G V, Rameev B Z, Dogan N, et al. Secondary signals in two-frequency nuclear quadrupole resonance on14N nuclei with I=1[J]. Journal of Magnetic Resonance, 2008, 193(1): 49-53.

[6] 李志強，金余桓. 基于核電四極矩共振技術的爆炸物檢測系統的數據處理及信號識別算法[J]. 核電子學與探測技術, 2004,24(6): 587-590.

LI Zhi-qiang, JIN Yu-heng. The data processing of the explosive detection system based on nuclear quadrupole resonance[J]. Nuclear Electronics & Detection Technology, 2004, 24(6): 587-590.

[7] 趙振維，婁揚，金燕波，等．基于自適應濾波技術的NQR信號處理[J]. 電波科學學報, 2008,23(3):429-432.

ZHAO Zhen-wei, LOU Yang, JIN Yan-bo, et al. Signal processing for NQR based on adaptive filtering[J]. Chinese Journal of Radio Science, 2008,23(3):429-432.

[8] Mozzhukhin G V, Molchanov S V. Application of the wavelet transform for detecting signals of nuclear quadrupole resonance[J]. Russian Physics Journal, 2005, 48(1):53-56.

[9] 楊振磊，徐更光，王振華，等. 基于小波變換的炸藥NQR信號處理[J]. 原子能科學技術,2010,44(3):354-357.

YANG Zhen-lei, XU Geng-guang, WANG Zhen-hua, et al. Processing method of14N nuclear quadrupole resonance signal in explosive[J]. Atomic Energy Science and Technology, 2010, 44(3):354-357.

[10] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hibert spectrum for nonlinear and non-stationary time series analysis[J]. Pro R Soc London, Ser A, 1998, 454:903-995.

[11] Huang N E, Wu M C, Long S R, et al. A conFIDence limit for the empirical mode decomposition and the Hibert spectral analysis[J]. Proc R Soc London, Ser A, 2003, 459:2317-2345.

[12] Donoho D L, Johnstone I M, Kerkyacharian G, et al. Wavelet shrinkage:asymptopia[J]. Journal of the Royal Statistical Society, Ser B, 1995, 57(2):301-337.

[13] 王宏強，尚春陽，高瑞鵬，等. 基于小波系數變換的小波閾值去噪算法改進[J]. 振動與沖擊, 2011, 30 (10): 165-168.

WANG Hong-qiang,ShANG Chun-yang,GAO Rui-peng, et al. An improvement of wavelet shrinkage denoising via wavelet coefficient transformation[J]. Journal of Vibration and Shock, 2011, 30 (10): 165-168.

[14] Boudraa A O, Cexus J C.EMD-based signal filtering[J]. IEEE Transactions on Instrumentation and Measurement, 2007, 56(6):2196-2202.

[15] Samuel D S, Andreas J, John A S, et al. Exploiting spin echo decay in the detection of nuclear quadrupole resonance signals[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(4):925-933.