基于经验模态分解及小波变换的炸药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.