宁文茜 ,王艳华 ,樊黎明,3* ,张晓峻 ,谢志臻

(1.西北工业大学 航海学院,陕西 西安,710077;2.自然资源部 海洋环境探测技术与应用重点实验室,广东 广州,510300;3.西北工业大学 青岛研究院,山东 青岛,266200;4.哈尔滨工程大学 物理与光电学院,黑龙江 哈尔滨,150001)

0 引言

磁性目标产生的磁场叠加在背景磁场上,使得其周围空间的磁场分布发生变化,从而形成磁异常。磁异常探测通过揭示环境磁场中的异常来检测铁磁物体,是一种广泛应用的被动目标探测方法[1]。该技术已在许多领域得到应用,如未爆炸弹药(unexploded ordnance,UXO)检测[2-4]、海底电缆检测[5]、交通监视[6]及人体医学研究[7]等。然而,磁异常会随着距离的增加而迅速衰减。

近年来,学者们提出多种磁异常检测的方法。根据特点,这些方法主要可以分为两类。一类是基于信号分析的磁异常检测方法[8-16]。Sheinker 等[8]提出标准正交基函数(orthonormal basis function,OBF)的匹配滤波,探测深埋在白高斯噪声中的磁异常信号。为了应对功率谱密度为1/fα(0 ＜α ＜2)的背景磁噪声,其设计了白化滤波器,能够提升磁异常的探测性能。该方法可以被认为是自回归(autoregressive,AR)过程,并且可以在高阶过滤器下有效工作。Liu 等[11]分析了空间分布系数的约束,并提出一种改进的基于空间投影系数变换的OBFs 方法用于磁异常检测。Fan 等[13]提出磁异常梯度信号的OBF 探测方法,该方法能够有效降低探测的虚警概率。假设磁噪声的模式变化是由磁异常引起,最小磁熵法[9]可以揭示出该模式的变化,从而达到检测磁异常的目的。该方法的优点是不需要对目标进行先验假设,计算复杂度较低。然而,该方法中噪声概率密度函数通过统计静止状态下的背景噪声获得。针对运动平台的探测需求,Fan 等[14]提出基于经验模态分解(empirical mode decomposition,EMD)和最小熵的磁异常探测方法。该方法通过经验模式分解可快速计算运动状态下的近似噪声概率密度函数,实现对目标的探测。为了提高低信噪比条件下的检测性能,Wan 等[10]提出了一种基于随机共振的磁异常检测方法。然而,随机共振参数的选取影响了该方法的有效性。此外,学者们还提出了混合算法来检测磁异常[15]。

另一类磁异常检测方法是基于机器学习的磁异常检测[17-20]。该类方法的主要思想是将磁异常检测视为分类问题,即确定磁信号中是否存在异常信号。Liu 等[19]提出了一种使用全连接神经网络的磁异常检测方法,其中选择了2 种类型的特征作为输入。Hu 等[17]提出了一种基于卷积神经网络的磁异常检测方法,将特征信息的提取和分类任务结合在一个单一体系中。Fan 等[18]提出了一种使用支持向量机的自适应磁异常检测方法,其中OBF 能量和磁熵作为磁异常特征输入到支持向量机。

通过磁噪声抑制、磁信号增强或磁异常多特征信息融合等方式,上述方法能实现在有效探测范围内识别磁异常信号,判断是否存在磁性目标。当检测到磁性目标后,在进一步定位与识别磁性目标时,需要获取目标的弱磁异常信号。针对低信噪比下弱磁异常信号的获取问题,文中提出基于集合经验模态分解(ensemble empirical mode decomposition,EEMD)和离散小波变换(discrete wavelet transformation,DWT)的弱磁异常提取方法。

1 弱磁异常信号分析

当磁性目标与磁传感器的距离大于目标最大尺度的3 倍时,可以把磁性目标视为一个磁偶极子。磁性目标产生的磁场可以表示为[21]

在实际测量中,测量磁场Bm由地磁场Be和目标产生的磁场Ba组成。由于磁场随距离的快速衰减,当目标远离磁传感器时,Ba远小于Be。当利用标量磁力仪探测远距离的目标时,磁场的标量形式可以表示为

由式(2),目标产生磁异常的标量形式可以表示为

由式(1)和式(3)可知,目标产生的磁异常信号随距离的3 次方衰减。随着探测距离的增大,目标磁异常信号将会变得微弱。因此,目标弱磁异常提取是开展远距离目标定位与识别的基础。

2 弱磁异常提取方法

2.1 方法概述

当探测远距离磁性目标时,由于磁性目标产生的磁场随着距离的增加快速衰减,导致目标的弱磁异常淹没在磁噪声中。通过磁异常探测方法确定目标存在后,需要对弱磁异常信号进行有效的提取,从而为后续的目标定位与识别提供高精度的数据基础。通常认为,磁异常信号的频率位于极低频范围[22]。

为有效抑制背景磁噪声并提取目标产生的弱磁异常信号,提出基于EEMD 和DWT 的弱磁异常提取方法。该方法框架图如图1 所示。首先,利用EEMD 将磁信号分解成一系列固有模态函数(intrinsic mode function,IMF)项和残差项。分解出来的各IMF 分量包含了原信号不同时间尺度的局部特征信号,包括从高到低的不同频段。根据噪声评估结果,IMF 被分为噪声域和信号域。噪声域的IMF 需要进行降噪处理,获得部分磁异常信号。其次,利用DWT 将噪声域IMF 分解为低频带的近似系数(AC)和相对较高频带上的一系列细节系数(DC)[23]。保留低频带的近似系数,利用逆离散小波变换(inverse discrete wavelet transform,IDWT)重新构建噪声域的磁异常。最后,将信号域磁异常信号项和重构的噪声域磁异常信号合成去噪的磁异常信号,从而实现微弱磁异常信号的提取,提高磁异常的信噪比。

图1 基于EEMD-DWT 的弱磁异常提取方法框图Fig.1 Block diagram of weak magnetic anomaly extraction method based on EEDM-DWT

2.2 EEMD

EMD 是由Huang 等[24]提出的一种处理非线性非平稳信号方法。该方法可将复杂的时间序列信号分解为少数几个具有不同特征尺度的时间序列,称为 IMF 和一个趋势项。因此,信号可以表示为

式中:IIMF,i(t)为第i个本征模函数;r(t)为趋势项。

EMD 方法从数据本身特征出发,不依赖于任何外部函数和参数,具有自适应性。然而,由于某些信号具有间歇性,分解时常常出现模态混合现象,即单独一个IMF 分量包含着不同尺度的信号或一种尺度的信号分布在不同的IMF 分量中,导致分解出的IMF 分量缺乏物理意义。

为克服这个问题,Wu 等[25]发展了EEMD 方法。该方法是通过在原信号中添加白噪声序列,再利用EMD 方法将添加了白噪声序列的信号分解成若干个IMF,重复多次且每次添加不同的白噪声序列,最后取多次分解得到的相应IMF 的算术平均值作为最终结果。该方法依据数据自身的时间尺度特征来进行信号分解,可以应用于任何类型信号的分解,在处理非平稳及非线性数据上,具有非常明显的优势,适合于分析非线性、非平稳信号序列,具有很高的信噪比。

通过EEMD 所分解出来的各IMF 分量包含了原信号的不同时间尺度的局部特征信号。根据噪声评估结果,IMF 被分为噪声域和信号域。信号域的IMF 需保留,噪声域中的IMF 需要进行降噪处理。

基于阈值信息,开展IMF 噪声域和信号域的区分。Kopsinis 等[26]提出了一种鲁棒的阈值选择策略,并推导了分解中分形高斯噪声的方差传递模型。根据噪声方差,推导出噪声能量在各阶IMF中的传递模型,然后根据各阶IMF 中的噪声能量确定相应的阈值。最后,根据阈值信息,确定对应的IMF 为信号域或噪声域。

通常,第1 层IMF 主要成分是噪声,噪声的能量可以通过噪声的标准差来估计,即

其他IMF 的噪声能量通过传递模型获得,即

式中,CH=(E1/βH),文献[27]给出了 βH和 ρH的取值说明。根据噪声能量,可通过构造阈值函数获得阈值。常用的阈值计算模型为

式中:Tk为第k层IMF 对应的阈值;Ek为第k层IMF 的噪声能量;N表示信号长度。

根据阈值信息,区分信号域和噪声域IMF 的步骤如下:

2.3 DWT 方法

DWT 具备处理信号低频部分的能力,同时能够较好地分解高频部分,收集信号各个频段的有用信息,进而提高去噪精度。因此,DWT 可将信号x(t)的时域表示分解为母小波的多个时移和缩放。数学上,小波ψa,b(t)由它的母小波ψ(·)通过缩放和平移得到,即

式中:a∈R+为尺度参数;b∈R为平移参数。将参数a和b离散化得到DWT。在DWT 中,根据2 的幂次选择参数:a=2j,b=k2j,其中k,j∈Z。因此,信号的DWT 可写成

通常,小波细节系数表征高频信号部分,小波近似系数表征低频信号部分。由于磁异常信号的频率位于极低频范围。因此,对噪声域部分xn(t)需要进行DWT,获得近似系数和细节系数。通过利用近似系数可重构出低频的弱磁异常信号实现多噪声域IMF 中的高频磁噪声抑制。

2.4 弱磁异常信号提取

低信噪比的弱磁异常信号经过EEMD 和DWT的去噪处理,获得去噪后的磁异常信号为

为了量化该方法的性能,使用信噪比 (signalnoise-ratio,SNR)和均方根误差(root mean squared error,RMSE)进行评价。评价指标的数学表达式为

式中:s(t)为目标的磁异常信号;x(t)为含有噪声的磁异常信号;x～(t)为噪声抑制后的磁异常信号。

3 实验验证

3.1 仿真实验

为了有效评价该方法的性能,设计1 个一维传感器阵列(如图2 所示)来收集磁背景噪声。磁传感器为高灵敏度的光泵磁力仪,其本征噪声约为选择哈尔滨市金沙滩湿地公园采集了真实的磁噪声样本,该区域环境磁活动非常低。地磁场局地倾角和偏角分别为63.3°和-10.1°。地磁场的振幅约为55 200 nT。磁力仪的采样频率设置为10 Hz。

图2 一维磁传感器阵列Fig.2 One-dimensional magnetic sensor array

典型的磁异常信号是由仿真目标产生的。目标以6.4 m/s 的恒定速度平行于X轴运动,从 (191,40,5) m 开始,在 (-192,40,5) m 结束。目标的磁矩在运动过程中保持不变,磁矩的大小为 200 A·m2。通过仿真可以得到典型的磁异常,然后将传感器阵列测量的真实磁噪声添加到模拟信号中,可以获得具有真实噪声和模拟磁异常的合成磁信号,如图3 所示。

图3 合成的磁异常信号Fig.3 Synthetic magnetic anomaly signal

基于EEMD-DWT 方法和基于EEMD 方法用于合成磁异常信号的信号提取。在基于EEMDDWT 方法中,噪声能量传递模型中相关参数设置如下:H=0.5,βH=0.719,ρH=2.010。尝试多组小波函数,确定最优的小波函数为db4。根据磁力仪的采样频率和磁异常的频率范围,确定最优的小波分解层数为5 层。基于EEMD 方法中的参数设置与基于EEMD-DWT 方法中的参数一致,分别利用这2 种方法对合成磁异常信号进行处理,提取磁异常的结果如图4 所示。通过图4 可知,基于EEMD-DWT 方法提取的磁异常信号更接近于目标信号。在基于EEMD 方法中,信噪比为1.72 dB。在基于EEMD-DWT 方法中,信噪比为5.27 dB,比基于EEMD 方法高出3.56 dB。

图4 不同方法的弱磁异常提取结果Fig.4 Results of weak magnetic anomaly signal extraction using different methods

2 种方法获得的RMSE 结果如图5 所示。基于EEMD-DWT 方法的RMSE 值为0.000 9,磁异常最大值处偏差的误差为4.36%。基于EEMD 方法的RMSE 值为0.002 4,磁异常最大值处偏差的误差为29.20%。因此,基于EEMD-DWT 方法能够有效提取弱磁异常信号。

图5 不同方法的RMSE 结果Fig.5 Results of RMSE using different methods

3.2 外场试验

为了验证所提方法在实际磁测中的有效性,在试验中将其应用于实际测量的弱磁异常提取,如图6 所示。在磁噪声采集试验中,磁测量系统由2 个铯光泵磁力仪(CS-L 和Scintrex)、1 个数据处理单元(MMS-4,PICO)和1 台计算机组成。该测量系统以10 Hz 的采样率用于获取磁场信号。磁测量系统放置在河道岸边,河道中行驶的拉沙船、摩托艇等作为磁性目标。

图6 河道岸边磁场测量现场Fig.6 Scene of magnetic field measurement at river bank

拉沙船的体积大,磁性材料多,其产生的磁异常信号较大,背景磁噪声对其磁异常信号影响有限。试验不对拉沙船的磁异常信息进行处理。由于摩托艇的体积较小,磁性材料少,产生的磁异常较弱,磁噪声对其磁异常信号影响大。因此,选择摩托艇产生的磁异常信号进行分析和处理。

利用基于EEMD-DWT 方法对摩托艇产生的磁异常信号进行处理,获得目标的磁异常信号,如图7 所示。试验结果显示,该方法能够有效抑制背景磁噪声的影响,获得目标的磁异常信号。相比于基于EEMD 方法,所提出的基于EEMD-DWT方法在弱磁异常信号提取中具有明显的优势。

图7 实测磁异常信号提取Fig.7 Measured magnetic anomaly signal extraction

4 结束语

为获取远距离目标产生的微弱磁异常,针对弱磁异常信号的极低频特性,文中提出基于EEMDDWT 的弱磁异常提取方法。采用EEMD 将微弱磁异常信号分解为信号域信号和噪声域信号。噪声域信号通过DWT 进行低频信号的提取。利用信号域信号和经DWT 处理噪声域信号的叠加,获得目标的弱磁异常信号。在仿真实验中,基于EEMDDWT 方法提取弱磁异常与真实磁异常的RMSE值为0.000 9,磁异常最大值处的误差为4.36%。外场试验中的结果显示基于EEMD-DWT 方法提取弱磁异常信号与实测数据的变化趋势一致。该方法提取的磁异常信号可为后续目标的定位与识别提供必要的数据基础。