郝鳳龍，徐更光，黃學(xué)義

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

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

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

本文通過分析當(dāng)前炸藥NQR信號處理方法優(yōu)缺點，提出基于經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)與小波變換聯(lián)合的NQR信號去噪方法。該方法不受傅里葉變換及小波函數(shù)選擇限制，自適應(yīng)能力好，并通過對實測炸藥NQR信號去噪處理，驗證該方法的可行性及有效性。

1 經(jīng)驗?zāi)B(tài)分解去噪原理

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

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

(2) 在極值點間利用三次樣條函數(shù)進行插值，獲得x(t)上包絡(luò)線u0(t)與下包絡(luò)線v0(t)，并計算兩包絡(luò)線均值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)，重復(fù)以上步驟，直到滿足IMF條件為止，求得第一階IMF分量c1(t)；

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

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

(3)

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

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

(4)

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

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

2 小波閾值去噪原理

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

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

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

(6)

(7)

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

3 經(jīng)驗?zāi)B(tài)分解與小波變換聯(lián)合去噪原理

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

(8)

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

(9)

其中：N為信號長度；IMFk(ti)為信號分解所得第k個分量?；谝陨戏治觯盘柲芰糠纸琰c可確定為

(10)

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

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

4 炸藥NQR信號去噪分析

4.1 NQR信號模型

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

(11)

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

4.2 黑索金NQR信號去噪實驗

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

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

圖2 原始黑索金NQR信號

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

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

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

圖4 固有模態(tài)函數(shù)分量及殘余量

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

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

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

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

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

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

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

5 結(jié) 論

本文基于NQR信號非線性與非平穩(wěn)性特點，提出基于經(jīng)驗?zāi)B(tài)分解(EMD)與小波變換聯(lián)合的炸藥NQR信號去噪方法。通過對探測的黑索金NQR信號進行經(jīng)驗?zāi)B(tài)分解，分析各階IMF分量特征信息，對噪聲起主導(dǎo)作用的高頻IMF分量進行小波軟閾值去噪，將去噪后高頻IMF分量與低頻IMF分量重構(gòu)信號，既可保留有用信號特征信息又能有效抑制噪聲、提高信噪比，很大程度上克服小波閾值去噪及直接EMD去噪的缺陷。該方法在炸藥NQR信號去噪中的良好性能，可為NQR信號處理奠定理論、技術(shù)基礎(chǔ)。

[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] 李康寧，俞碩，李興，等. 核四極共振技術(shù)在黑火藥探測中的研究[J]. 核科學(xué)與工程, 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] 李志強，金余桓. 基于核電四極矩共振技術(shù)的爆炸物檢測系統(tǒng)的數(shù)據(jù)處理及信號識別算法[J]. 核電子學(xué)與探測技術(shù), 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] 趙振維，婁揚，金燕波，等．基于自適應(yīng)濾波技術(shù)的NQR信號處理[J]. 電波科學(xué)學(xué)報, 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]. 原子能科學(xué)技術(shù),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] 王宏強，尚春陽，高瑞鵬，等. 基于小波系數(shù)變換的小波閾值去噪算法改進[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.