CN110967774A - 一种基于传感器阵列的磁异常检测方法 - Google Patents

Abstract

本发明提供了一种基于传感器阵列的磁异常检测方法于，包括：步骤S1：数据预处理，利用经验模态分解去除测量传感器与参考传感器的差量磁场中的地磁背景趋势项；步骤S2：信号调制，对预处理后的差量磁场进行调制，实现数据融合；步骤S3：对调制后的差量磁场进行时频分析与量化，实现目标磁异常的检测。

Description

一种基于传感器阵列的磁异常检测方法

技术领域

本发明涉及入侵目标磁法检测，用于针对固定区域或重要据点的实时监测与安全防卫应用中，尤指一种基于传感器阵列的磁异常检测方法。

背景技术

随着陆地资源开发难度与成本的日益增加，各国对海洋资源的开发愈加重视，对公共海域或争议岛屿的主权争夺日趋激烈。我国领海主权及合法海洋权益正面临敌对势力蓄意入侵与日益猖獗破坏活动等造成的严重威胁。传统的航空磁探测技术主要用于重点海域与可疑目标巡航探测，而对于固定海域的实时安全防御则需在敏感水域布设传感器网络实现。由于地磁等背景磁场噪声影响，目标的磁异常常常淹没于噪声中，对目标精准探测带来极大的困难与挑战，研究基于静止平台磁探测技术对水下入侵目标的实时监测与识别技术具有重要价值。传统磁异常检测方法主要分为两大类：(1)基于目标磁异常特征建模法，最典型的代表是正交基匹配滤波检测算法，将磁性目标磁场建模获取磁异常的特征基函数，通过匹配滤波的形式达到检测目标磁异常并抑制噪声目的，；(2)基于地磁背景噪声特征建模法，最典型的代表为最小熵滤波检测算法，通过统计分析背景磁噪声，利用背景熵值高而信号熵值低特性，达到检测目标磁异常的目的。

正交基匹配滤波算法原理如图1所示，其中h(n)为白化滤波器， g_i(n)，i＝1，2，3为正交基函数，

表示卷积运算。信号x(n)经白化滤波后，利用滑动时间窗口与正交基函数匹配，计算各窗口的匹配系数的平方相加作为该时刻的探测指数。根据奈曼-皮尔逊(Neyman-Pearson，N-P)准则设定不同信噪比的磁异常信号的探测指数门限值，以达到最优探测效果。

最小熵滤波检测算法是依据背景噪声与磁异常的统计特性差异来检测目标信号的。算法可分三部分：1.对长期采集噪声样本的分布频率进行统计；2.用统频率代替分别概率，计算样本熵值；3.根据检测性能需求，根据N-P准则确定检测阈值，判断不同滑动时间窗内是否出现磁异常。

用序列x_i，i＝1，2，…，M表示噪声样本，假设其服从正态分布。使用样本期望值μ及样本标准差σ可以将其概率密度函数表达如式(1)。其中样本序列的均值与方差如式(2)。

则每个样本点x_i的概率可表达如式(3)。

磁场信号经传感器采集，经过模/数转换为数字信号。假设量化步长为δ_x，则对于每一个样本点的概率可用式(4)近似计算。

p(x_i)＝f(x_i)δ_x (4)

利用熵滤波器根据式(5)计算长度为L的滑动窗口内的熵值和。

最后，根据N-P准则确定探测阈值，当熵滤波器输出值低于阈值时，触发磁异常检测。

发明内容

(一)要解决的技术问题

常用的正交基匹配与最小熵滤波检测算法分别存在如下不足：(1)基于正交基匹配滤波的检测算法，需要满足磁性目标的磁偶极子目标模型、磁偶极子目标做匀速直线运动、运动过程磁偶极子目标磁矩保持不变，这些假设条件降低了检测算法的适应能力。(2)基于最小熵滤波的检测算法，通过磁异常存在时熵变小特性检测目标，对于非平稳的背景磁场其熵并非像高斯白噪声一样稳定，很难保证在实际地磁背景下对于弱磁异常信号的检测性能。因此，本发明所要解决的技术问题是如何提高检测算法的适应能力、保证在实际地磁背景下对于弱磁异常信号的检测性能。

(二)技术方案

本发明提供了一种基于传感器阵列的磁异常检测方法，

步骤如下：

步骤S1：数据预处理，利用经验模态分解去除测量传感器与参考传感器的差量磁场中的地磁背景趋势项；

在步骤S1中，对于共有N个传感器的阵列，选择居中位置传感器为参考传感器，其他位置的传感器为测量传感器，计算测量传感器与参考传感器的差量磁场，通过地磁场的空间相关性为传感器阵列抑制背景磁场噪声。

其中，采用经验模态分解(Empirical Mode Decomposition，EMD)将差量磁场中不同尺度的波动和趋势逐级分解，产生一系列具有不同尺度的序列，获得稳定的差量磁场。

步骤S2：信号调制，对预处理后的差量磁场进行调制，实现数据融合；

对预处理后的差量磁场进行调制包括：

步骤S2-1：用y_r1(n)，y_r2(n)，…，y_rr(n)，…，y_rN(n)表示预处理后的差量磁场信号序列，其中y_rr(n)＝0；

步骤S2-2：按传感器编号，对N个差量磁场序列进行调制，调制后的信号可表示为：

…，y_r1(n-1)，y_r2(n-1)，y_r3(n-1)，…，y_rr(n-1)，…y_rN(n-1)，y_r1(n)，y_r2(n)，y_r3(n)，…， y_rr(n)，…y_rN(n)，y_r1(n+1)，y_r2(n+1)，y_r3(n+1)，…，y_rr(n+1)，…y_rN(n+1)，……，其中n代表序列的第n个采样点。

其中，调制后的差量磁场为以采样率为中心频率的交变信号。

步骤S3：对调制后的差量磁场进行时频分析与量化，实现目标磁异常的检测。

在步骤S3中，对调制后的差量磁场进行短时傅里叶变换，获得调制后的差量磁场的瞬时频谱特性，对瞬时频谱特性进行实时分析。

其中，短时傅里叶变换的表达式为：

；OUT(m)为磁异常检测方法的量化输出，

OUT(m)＝max(STFT(m,n))，n＝p，…，q (10)

式中，p、q分别为调制后的差量磁场频率的下边界与上边界。

(三)有益效果

从上述技术方案可以看出，本发明具有以下有益效果：

(1)通过对传感器阵列数据进行预处理，成功消除了非平稳背景地磁场的影响。

(2)通过信号调制方法，实现了阵列磁异常信号融合。

(3)与传统的时域检测算法不同，新算法在时频域进行，增强检测弱磁异常的性能。

(4)检测算法不需要预先估计目标的特征及其运动状态，也不需要考虑背景磁场变化的干扰，因此具有良好的环境适应性。

附图说明

图1为正交基匹配滤波算法原理图。

图2为基于传感器阵列的调制磁异常检测方法原理图。

图3为仿真阵列测量磁异常信号。

图4为磁异常信号调制过程。

图5为调制信号时频分析示意图。

图6为本发明实施例的传感器阵列分布与目标运动轨迹示意图。

图7是本发明实施例的示意图，其中，(a)为本发明实施例的不同传感器测得的磁异常示意图；(b)为本发明实施例的差量磁异常示意图。

图8是本发明实施例的另一示意图，其中，(a)、(b)、(c)为本发明实施例的不同信噪比时，检测算法输出示意图，(d)为本发明实施例的不同信噪比时与对应检测概率的误警率示意图。

图9是本发明实施例的又一示意图，其中，(a)为本发明实施例的测量磁场处理前的低通滤波的测量信号示意图；(b)为本发明实施例的测量磁场处理前的计算差量场信号示意图；(c)为本发明实施例的经EMD处理后的差量信号示意图；(d)为本发明实施例的经EMD处理后差量信号的幅度概率分布示意图。

图10是本发明实施例的又一示意图，其中，(a)为本发明实施例的差量场信号EMD处理前的功率谱密度示意图；(b)为本发明实施例的差量场信号EMD处理后的功率谱密度示意图。

图11是本发明实施例的再一示意图，其中，(a)为本发明实施例的传感器I的测量磁场与目标磁异常信号示意图；(b)为本发明实施例的传感器III的测量磁场与目标磁异常信号示意图。

图12为本发明实施例的SNR＝-10.5dB时，目标磁异常的检测结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向。因此，使用的方向用语是用来说明并非用来限制本发明的保护范围。

本发明实施例提供了一种基于传感器阵列的磁异常检测方法，如图2 所示，步骤如下：

步骤S1：数据预处理，利用经验模态分解去除测量传感器与参考传感器的差量磁场中的地磁背景趋势变化项；

磁场传感器检测到磁异常信号的同时也会包含背景噪声，为提高磁异常信号的检测能力和可靠性，本发明通过地磁场的空间相关性为传感器阵列提供了抑制背景磁场噪声的良好途径。

对于共有N个传感器的阵列，选择居中位置传感器作为参考传感器，计算其他测量传感器与参考传感器的差量磁场：假设参考传感器的信号 b_r(n)，记参考传感器外的第i个测量传感器的信号为b_i(n)，i＝1，2，…，r-1，r+1，…，N，计算第i个测量传感器与参考传感器的差量磁场信号表达如式(6)所示，n表示样本序号。

b_ri(n)＝b_i(n)-b_r(n) (6)

当参考传感器与测量传感器距离较近时，各传感器测量磁场中的磁源位于磁层或电离层的变化磁场近乎完全一致，理论上利用计算差量场可完全消除变化地磁场的影响。然而，实际磁场测量环境下由于工作环境(如工作温度等)变化及磁场各传感器自身性能的差异或者发热等情况，阵列中各传感器的输出磁场并非完全一致，计算的差量磁场中出现趋势项，由于差量磁场中趋势项的成因复杂且非平稳，很难通过传统建模方式进行消除。本发明采用经验模态分解(Empirical Mode Decomposition，EMD)将信号中不同尺度的波动和趋势逐级分解，产生一系列具有不同尺度的序列，将趋势项消除。

其中，经验模态分解(Empirical Mode Decomposition，EMD)是黄锷博士(NEHuang)在美国国家宇航局与其他人于1998年创造性提出的。 EMD实质上是对非线性非平稳信号进行平稳化处理的一种手段，将信号中不同尺度的波动和趋势逐级分解，产生一系列具有不同尺度的序列，称为固有模态函数(Intrinsic Mode Function,IMF)。EMD可将时域信号x(t)分解为：

其中，c_i(t)为本征模函数，IMF分量包含了原信号的不同时间尺度的局部特征信号；r(t)频率成分最低，作为趋势项。

磁场传感器测得的磁信号x(t)包含地磁背景场、传感器本底噪声及潜在的目标磁异常，其中地磁场随时间的漂移及传感器的温度漂移等形成变化的趋势项，鉴于磁异常场与背景磁场特征的区别，对原始测量信号进行经验模态分解，实现保留磁异常信号并去除地磁趋势变化，提升磁异常信号的检测能力。

步骤S2：信号调制，对预处理后的差量磁场进行调制，实现数据融合；

通过对阵列磁场测量数据进行预处理后，可获得稳定的差量磁场。当磁异常信号不存在时，传感器阵列的输出表现为稳定噪声；当磁异常存在时，考虑目标和传感器阵列间距离不同，各传感器测量的磁异常强度差异显著，甚至极性不同。

通过对预处理信号进行调制处理，实现阵列测量磁异常信号的融合。

信号调制过程步骤如下：

步骤S2-1：用y_r1(n)，y_r2(n)，…，y_rr(n)，…，y_rN(n)表示预处理后的差量磁场信号系列，其中y_rr(n)＝0。

步骤S2-2：按传感器编号，对N个差量磁场序列进行调制，调制后的信号可表示为：

……，y_r1(n-1)，y_r2(n-1)，y_r3(n-1)，…，y_rr(n-1)，…y_rN(n-1)，y_r1(n)，y_r2(n)，y_r3(n)，…， y_rr(n)，…y_rN(n)，y_r1(n+1)，y_r2(n+1)，y_r3(n+1)，…，y_rr(n+1)，…y_rN(n+1)，……，其中n代表序列的第n个采样点。

其中，当磁异常不存在时，噪声信号调制后与调制前无明显特征差异；当磁异常存在时，准静态的磁异常经调制处理后形成以采样率为中心频率的交变信号。如图3所示，以3个差量磁异常为例，对调制过程进行形象化描述。假设测量传感器I、II、III的差量磁异常信号如图3所示，采样率为f_s＝2Hz，采集时间为25秒。按照步骤S2-2对该信号进行调制，调制信号如图4中黑色实线所示，三个测量传感器信号序列经调制后变为一个序列，因此调制后的信号样本点数为3个测量传感器样本点的总和，在相同时间内调制信号的等效采样率为f_e＝3f_s＝6Hz。由图4可以看出相比于准静态的差量磁异常，调制信号中出现交变单频信号，且该信号中心频率为采样率f_s。

步骤S3：对调制后的差量磁场进行时频分析与量化，实现目标磁异常的检测。

短时傅里叶变换作为非平稳信号分析工具，通过选择适当的窗函数在本地截断信号，对截断信号进行短时傅里叶变换用于获得瞬时频谱特性，最终达到实时分析信号频谱特性的目的。对于连续信号z(u)，它的短时傅里叶变换可表达如式(8)所示，其中g(u-t)表示窗函数，*表示共轭。

在实际应用中，STFT(t，f)需要离散化，即STFT(t，f)要被等间隔的时、频采样，其中T＞0且F＞0分别代表了时域和频域的采样周期，m 和n是整数。为了简化表示，将STFT(mT，nF)用STFT(m，n)表示。因此，对于离散信号z(k)的短时傅里叶变换可以表达如式(9)。

其中，z(u)表示连续信号，z(k)表示离散信号，利用公式(9)对图4 中的步骤S2-2得到的差量磁场序列调制信号进行短时傅里叶变换并进行时频分析，结果如图5所示，可见调制信号为中心频率为2Hz的单频信号，与调制原理理论分析结果一致。

时频分析为磁异常的检测提供良好的视觉效果，但为定量评估算法性能，需量化时频分析结果。定义OUT(m)作为检测算法的量化输出，可用式(10)表示。式中p、q分别表示了调制信号频率的下边界与上边界，由调制原理可得，将差量磁场序列调制信号的短时傅里叶变换结果代入公式 (10)，得到目标磁异常检测的量化输出，调制后的磁异常信号频率在采样频率附近。

OUT(m)＝max(STFT(m，n))，n＝p，…，q (10)

传感器阵列分布和目标运动轨迹如图6所示。传感器阵列由7个成线性分布的磁力仪组成。目标参数如表1所示：

表1磁异常目标参数

其中，目标的磁异常信号如图7(a)所示，以传感器4为参考，则计算差量磁异常如图7(b)所示。

根据图7(a)磁异常的幅度，因此通过在磁异常信号中依次添加标准差不同的高斯白噪声获取信噪比分别为-18.1dB、-20.6dB和-22.5dB的含噪磁异常信号。针对不同信噪比的含噪磁异常信号，进行10000次随机实验。图8中展示含磁异常与不含磁异常信号的量化输出及检测概率密度函数。当误警率为0.001时，检测阈值由N-P准则确定，不同信噪比的磁异常检测结果如表2所示。以信噪比为-20.6dB为例，由N-P准则确定阈值为1.330，检测概率为0.9648。

表2不同信噪比下算法检测性能

为测试本发明在真实环境下的检测性能，开展地磁背景噪声中叠加磁异常信号的检测实验，采用使用CS-3光泵磁力仪作为地磁测量传感器，磁力计固有噪声为

采样率为40Hz，由于磁异常信号频率很低，因此对测量信号首先进行低通滤波处理(截至频率2Hz)。为与仿真一致，将低通处理后的信号降采样到10Hz。实验采用3个光泵磁力仪，间隔距离为15m。

将传感器2作为参考传感器，计算差量磁场并进行预处理。测量传感器的测量信号与处理后信号如图9，其中图9(a)为低通滤波的测量信号，图9(b)为计算差量场信号，图9(c)为EMD处理后的差量信号，图9(d)为 EMD处理后差量信号的幅度概率分布。从图9(c)、图9(d)可知差量场中趋势项经EMD处理后完全消除，且处理后的磁场近似为高斯分布。计算图9中(b)、(c)中信号的功率谱密度，结果分别如图10所示，表明预处理后的差量场信号在0-1Hz间近似为白噪声。综上所述，经预处理后的各传感器输出信号均近似为高斯白噪声，验证了仿真实验中添加噪声的合理性。

将图7(b)中Mag-3、Mag-5的差量磁异常添加至图9(b)差量磁场中进行地磁背景磁异常检测实验，经过EMD去除趋势项获得测量磁场如图11，其中信号的SNR约为-10.5dB。目标磁异常检测结果如图12所示，目标信号可以从时频图中清楚的识别出来。

至此，已经结合附图对本实施例一种基于传感器阵列的磁异常检测方法进行了详细描述。依据以上描述，本领域技术人员应当对本发明的磁异常算法有了清晰的认识。

此外，上述对方法的定义并不仅限于实施例中提到的各种具体方式，本领域普通技术人员可对其进行简单地更改或替换，例如：

(1)本发明中采用EMD的方法去除差量磁场的趋势变化，经实验验证可靠、有效，但预处理的方法不局限于EMD方法，还包括诸如如自适应滤波、带通滤波等方法。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于传感器阵列的磁异常检测方法，其特征在于，包括：

步骤S1：数据预处理，利用经验模态分解去除测量传感器与参考传感器的差量磁场中的地磁背景趋势项；

步骤S2：信号调制，对预处理后的差量磁场进行调制，实现数据融合；

步骤S3：对调制后的差量磁场进行时频分析与量化，实现目标磁异常的检测。

2.根据权利要求1所述的一种基于传感器阵列的磁异常检测方法，其特征在于，在步骤S1中，对于共有N个传感器的阵列，选择居中位置传感器为参考传感器，其他位置的传感器为测量传感器，计算测量传感器与参考传感器的差量磁场，通过地磁场的空间相关性为传感器阵列抑制背景磁场噪声。

3.根据权利要求1或2所述的一种基于传感器阵列的磁异常检测方法，其特征在于，采用经验模态分解将差量磁场中不同尺度的波动和趋势逐级分解，产生一系列具有不同尺度的序列，获得稳定的差量磁场。

4.根据权利要求1或2所述的一种基于传感器阵列的磁异常检测方法，其特征在于，对预处理后的差量磁场进行调制包括：

步骤S2-1：用y_r1(n)，y_r2(n)，…，y_rr(n)，…，y_rN(n)表示预处理后的差量磁场信号系列，其中y_rr(n)＝0；

步骤S2-2：按传感器编号，对N个差量磁场序列进行调制，调制后的信号可表示为：

……，y_r1(n-1)，y_r2(n-1)，y_r3(n-1)，…，y_rr(n-1)，…y_rN(n-1)，y_r1(n)，y_r2(n)，y_r3(n)，…，y_rr(n)，…y_rN(n)，y_r1(n+1)，y_r2(n+1)，y_r3(n+1)，…，y_rr(n+1)，…y_rN(n+1)，……，其中n代表序列的第n个采样点。

5.根据权利要求1或4所述的一种基于传感器阵列的磁异常检测方法，其特征在于，调制后的差量磁场为以采样率为中心频率的交变信号。

6.根据权利要求1所述的一种基于传感器阵列的磁异常检测方法，其特征在于，在步骤S3中，对调制后的差量磁场进行短时傅里叶变换，获得调制后的差量磁场的瞬时频谱特性，对瞬时频谱特性进行实时分析。

7.根据权利要求6所述的一种基于传感器阵列的磁异常检测方法，其特征在于，

短时傅里叶变换的表达式为：

；

OUT(m)为磁异常检测方法的量化输出，

OUT(m)＝max(STFT(m，n))，n＝p，…，q (10)

式中，p、q分别为调制后的差量磁场频率的下边界与上边界。

Images