基于短时互相关算法对航空磁异常信号的检测

磁异常是指在地磁场背景下，铁磁性物体被磁化后产生了的磁场与地磁场重合在一起，使周围磁场发生一定程度的扰动，从而与原磁场存在变化的现象。使用高精度磁探设备，对产生磁异常的铁磁性目标进行检测、识别和定位的技术称为磁异常检测(Magnetic Anomaly Detection,MAD)。随着人类科技的进步与发展，人类通过对地磁扰动的测量，实现了地下矿藏的探测[1]、地下炸弹的定位[2]、并能够对地震期间地磁场的观测[3]，以期实现地震预报。以上这些领域的进步发展，离不开地磁场的实时、无源、全地域、无辐射等特性。

在地磁场长期的磁化作用下，潜艇会在其周围产生磁场。艇艇体主要由铁磁性物质建造而成，即使地磁场不存在，铁磁性物质天然具有磁性，并且铁磁性物质的磁滞效应而形成的历史磁化产生的磁场，称为固定磁场[4]；潜艇长期处于地磁场背景下，由于磁化过程对其的不断作用，产生了的磁场称为感应磁场，感应磁场的特性与潜艇的地理位置、运动态势有着密切的联系。潜艇在建造过程中所采用的防腐蚀材料以及壳体材料在海水中产生了的稳恒电流，会在潜艇周围激发出的磁场，称为腐蚀相关磁场[5](Corrosion-related magnetic field,CRM)。潜艇艇体内部设备在通信过程中产生的电磁波和其电流产生的辐射激发出的电磁场，也是潜艇磁场的重要来源。潜艇在地磁场下对其周围磁场产生扰动而形成的磁异常现象，会破坏潜艇的隐身性能。通过对地磁场异常扰动的探测，并对扰动信号进行识别，就能够实现对潜艇目标的探测、识别与定位。

磁异信号检测算法有正交基(Orthogonal Basis Function,OBF)分解[6]的方法及其延续算法[7-10]。随机共振检测[11]、主成分检测方法[12]、熵滤波器检测[13]、高阶过零检测[14]等。OBF分解对于目标磁异常信号的检测需要知道水下目标位置的先验信息；随机共振算法对于磁异常信号的检测在信噪比SNR=-3.5 dB，虚警率Far=1.5%情况下,检测概率仅有60%。熵滤波器检测能力低于随机共振检测。HOC检测器概率信噪比SNR=-2 dB，虚警率Far=0.1%情况下，检测概率为50%。

本研究基于潜艇磁异常信号时域特征，提出建立一个合适的参考信号，并基于短时互相关的方法来对磁探信号进行检测。当选取的参考信号与磁异常信号时域特征相近时，互相关算法就能够克服噪声干扰，使得潜艇磁异常信号在互相关输出图上凸显出来。仿真显示，短时互相关算法在信噪幅度比Rsn=3 dB的情况下，仍能检测出潜艇磁异常信号。

1 模型建立

1.1 互相关模型

互相关的表示两组信号的相似程度。当选取出合适的参考信号x(t)后，通过检测x(t)与待检测采样信号y(t)两者的相似程度，就能够对信号进行有效地检测。其中，对于y(t)来说，其信号检测为双择假设：

H0为真时，定义为采样的信号中不包含磁异常信号，只存在噪声；H1为真，定义为采样信号中含有磁异常信号。式中，n为噪声，s为磁异常信号。

互相关运算模型如下：

式中：y为选取的参考信号；x是采样信号，主体是噪声，有可能存在磁异常信号

如图1所示，互相关运算在输出序列尾部会有明显的波动，从而影响到互相关检测器门限的选取。所以需要对经典的互相关运算方法进行改进，使互相关的方法能够应用到磁异常信号检测中。

1.2 短时互相关

为了避免尾部波动对磁异常信号检测的影响，尝试对采样信号加窗，并通过移窗，与参考信号x作互相关运算，并对输出的尾部若干数据点进行截断,从而去除了互相关输出的尾部异常波动，算法示意图如图2。

当信号窗口为200个数据点时，短时互相关对信号的检测能力较好。

其数学模型为：

其中：m为窗口的中心位置，窗口长度为2N+1；x是采样信号；y为选取的参考信号。

图3～图5是对短时互相关的尾部无异常波动的仿真验真信号曲线。仿真得到，在Rsn=-10 dB的情况下，磁异常信号出现在29～31 s处，且基本淹没在噪声之中，其采样信号图如图3所示。

将29～31 s窗口中的采样信号y与参考信号x作互相关运算，其输出曲线如图5所示。

图5中，对采样信号29～31 s的互相关运算输出，于30 s附近出现最值，其幅度为4.67。能够有效监测处信号，并且没有尾部异常波动。

2 参考信号的选取

由于相关检测的方法是用参考信号和待检测信号的相似程度作为检测，但是由于磁异常信号受目标磁矩、磁倾角、磁偏角、飞机与水下目标的航向、航速的影响。变化较为剧烈，波形变化较大，所以对参考信号的选取相当重要。

对信号进行观测后发现，磁异常信号可以看作由若干个正弦脉冲信号的组合。提取5 000次仿真随机情况下的水下目标高空磁异常信号的宽度，其统计直方图如图6所示，均值为2.05 s，方差为0.373。

故把周期T=4、幅度A=1的半周期脉冲正弦信号作为磁异常相关检测的参考信号(图7)，且该参考信号幅度为1。

3 仿真验证

3.1 磁异常信号的获取

用磁偶极子模型仿真目标周围磁场数据，通过等效磁源对潜艇磁场的分布进行模拟的方法叫做等效磁源法，主要的等效磁体模型有单个磁偶极子模型、旋转椭球体模型和磁偶极子阵列模型等方法。张朝阳和肖昌汉在文献[16]中指出，在2.5倍物体长度以上的空间，完全可以吧磁性物体是做磁偶极子处理，满足一般工程应用的需求。潜艇一般与反潜机距离百米以上，磁偶极子模型作为等效磁源对潜艇磁场进行仿真，是完全可行的。

在航空磁探反潜过程中[17]，潜艇作为铁磁性目标，其磁场H可以表示为

式中：M是磁偶极子磁矩，在笛卡尔坐标系下可以表示为[Mx、My、Mz]类型的三分量形式；r0是潜艇的等效磁偶极子中心至预测空间P点的单位向量；

是潜艇的等效磁偶极子中心至预测空间P点的距离；P点坐标为(x,y,z)，o(1/r5)是高阶量，其强度以距离的5次方进行衰减，当r足够大，高阶量o(1/r5)可以在计算过程中忽略不计，即潜艇的等效磁偶极子模型磁场强度可以表示为

通过计算，可以得到

式中：H为等效磁偶极子磁场强度；Hx,Hy,Hz是磁场强度的三分量；H=ixHx+iyHy+izHz，ix、iy、iz为坐标系中的单位向量。

根据式(5)可得，B=μ0H可得，潜艇等效磁偶极子模型的磁感应强度三分量值为Bx,By,Bz可以通过式(7)计算而得，

某地地磁场HE,其三分量分别为(HEX,HEY,HEZ)，其单位方向向量为iE=(ixHEX+iyHEY+ izHEZ)/||HE||2。其中||HE||2表示为地磁场矢量的2-范数，如图8所示。H为磁偶极子模型下的潜艇的磁场矢量，其三分量分别为(Hx,Hy,Hz)。X地磁场的潜艇磁场叠加后形成的干扰磁场的总磁场，表达如下

设HE与H的夹角为γ，则

因HE>>H，所以式(9)可以展开为：

目前航空磁探反潜的探测数据是尾端光泵磁力计所测的标量磁场数据，所以磁偶极子模型下对潜艇磁场的仿真中，将三分量磁场投影到地磁场的坐标系上，得到其磁异常数值X，并表示如下：

即测得的标量磁异常信号X为潜艇信号H在地磁场HE方向的投影。

为了不失一般性。在某地磁场环境下(BEX= 28 900.3 nT, BEY=3 896 nT, BEZ=44 247.5 nT)，磁探平台高度100 m过顶飞过目标，采样率为100 Hz，飞机航速100 m/s，目标出现时间为t=30 s，仿真时间为0～60 s，共6 001个数据。分别仿真目标不同航向下所探测到的磁异常曲线如图9所示。

3.2 检测能力仿真

为了验证算法的有效性，用仿真数据进行仿真分析。在此，定义信噪幅度比Rsn，并表示为对数形式：

式中：Smax为磁异常信号的峰值；std(n)表示为噪声的标准差。

仿真得到了目标航向0°至360°(间隔10°)的磁异常信号，共计36组磁异常信号曲线。每组磁异常信号含6 001个数据点，参考信号为sin(pi/2×t)，0<t≤2的正弦脉冲。脉冲信号共计200个数据点。对磁异常信号进行加窗，逐窗与参考信号(Reference Signal)作短时互相关，如图10。

当目标航向为0°、90°、180°、270°，参考信号为正脉冲信号，Rsn为-2 dB时，对目标进行短时互相关运算，通过仿真发现当窗口位置m在0～2 800处，即0～28 s附近时，互相关系统的输出为随机数值，可以认为是噪声和参考信号的互相关；当窗口位置在2 800～3 200处，即28～32 s时，互相关输出有较大起伏，同时该时间点也是目标出现的时间。可以认为，短时互相关系统在高斯白噪声下对磁异常信号有较好的检测性能；当窗口位置m在3 200～5 000处，窗口不再包含磁异常信号，互相关输出又恢复成随机分布，即参考信号只与噪声作互相关运算。

短时互相关算法通过移窗的运算顺序，依次对采样信号不同时间段的信号进行互相关运算，且不同次序的窗口有重叠的时间段，这就使得信号的检测有一定的规则可循。即当伪彩图连续3～4 s出现峰值区域时，即可判定采样信号x(t)中出现潜艇磁场信号，即式(1)中的双择假设结果为H1。这种判决规则也是短时互相关算法优于互相关算法的体现。

图11为互相关检测器的ROC性能曲线图，当Rsn=3 dB时，检测虚警率为20%时，有90%的检测概率；当Rsn=4 dB，检测虚警率为10%时，有98%的检测概率，检测虚警率为3%时，有90%的检测概率；当Rsn=5 dB，检测虚警率为1%时，有99%的检测概率。互相关检测算法计算速度快，且无需目标先验信息，具有工程价值。

4 结论

短时互相关对于高斯白噪声下的磁异信号有较好的检测效果，能够在Rsn=3 dB的高斯白噪声干扰下检测出信号的存在与否。由于采取移窗运算的方法，使得存在一种新的判决思路，即短时互相关结果图存在连续2s高亮区域，即可判定存在磁异常信号；短时互相关检测无需知晓目标的先验信息，且算法容易实现，具有工程价值。下一阶段的工作，需要克服低信噪比环境，实现信号的高准确率检测。

[1]RIGGS L S,MOONEY J E,LAWRENCE D E J I T O G,et al.Identification of metallic mine-like objects using low frequency magnetic fields[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(1):56-66.

[2]周帅,黄大年,王泰涵.利用磁法数据的解析信号探测未爆炸弹(UXO)[J].地球物理学进展,2016(4):1767-1770.

[3]HAYAKAWA M,KAWATE R,MOLCHANOV O A,et al.Results of ultra-low-frequency magnetic field measurements during the Guam Earthquake of 8 August 1993[J].Geophysical Research Letters,1996,23(3):241-244.

[4]高俊吉,刘大明,肖昌汉,等.舰船固定磁场分量分离的一种新方法[J].海军工程大学学报，2007,19(2):86-9.

[5]陈聪,李定国,蒋治国,等.深海中潜艇腐蚀相关磁场全空间分布特征分析[J].哈尔滨工程大学学报,2014(06):684-689.

[6]GINZBURG B,FRUMKIS L,KAPLAN B Z,et al.Processing of magnetic scalar gradiometer signals using orthonormalized functions[J].Sensors,2002,102(1):67-75.

[7]SHEINKER A,FRUMKIS L,GINZBURG B,et al.Magnetic Anomaly Detection Using a Three-Axis Magnetometer[J].IEEE Transactions on Magnetics,2009,45(1):160-7.

[8]张坚,林春生,黄凡.OBF分解与BP网络在船舶磁场信号检测中的应用[J].船电技术,2011,31(7):13-6.

[9]张坚,林春生,邓鹏,等.基于小波域OBF分解的磁异常信号检测算法[J].弹箭与制导学报,2011,31(6):187-189.

[10]王珺琳,陈博.基于运动目标的OBF方法检测性能分析[J].中国电子科学研究院学报,2017,12(1):47-51.

[11]WAN C,PAN M,QI Z,et al.Magnetic anomaly detection based on stochastic resonance[J].Sensors & Actuators A Physical,2018,278:11-17.

[12]SHEINKER A,MOLDWIN M B,TECHNOLOGY.Magnetic anomaly detection (MAD) of ferromagnetic pipelines using principal component analysis (PCA)[J].Measurement Science，2016,27(4):045104.

[13]SHEINKER A,SALOMONSKI N,GINZBURG B,et al.Magnetic anomaly detection using entropy filter[J].Measurement Science,2008,19(4):045205.

[14]SHEINKER A,GINZBURG B,SALOMONSKI N,et al.Magnetic Anomaly Detection Using High-Order Crossing Method[J].Transactions on Geoscience,2012,50(4):1095-103.

[15]翁行泰，曹梅芬，吴文福，等.磁异探潜中潜艇的数学模型[J].上海交通大学学报,1995,03):27-32.

[16]张朝阳,肖昌汉,高俊吉,等.磁性物体磁偶极子模型适用性的试验研究[J].应用基础与工程科学学报,2010,18(5):862-868.

[17]嵇斗,王向军,柳懿,等.地磁场中潜艇运动感应静电场建模[J].四川兵工学报,2014(7):4-5.