基于谱减法及小波分析的海洋磁测噪声抑制方法

摘要

本公开提供了一种基于谱减法及小波分析的海洋磁测噪声抑制方法，包括：步骤S1，建立海浪磁场模型；步骤S2，对带有海浪磁场噪声的磁异常检测信号进行的谱减去噪；步骤S3，确定分解层次，对处理后的信号进行小波分析。利用模型仿真海浪磁场的时域和频域信息不会引入外部因素干扰，在确保模型准确的前提下该方法更具有普适性；同时通过利用谱减法和小波分析结合的方法，能够较大幅度地改善时域波形，抑制海浪噪声，并在信噪比的提升方面表现优异。

说明书

技术领域

本公开涉及磁异常目标探测领域，尤其涉及一种基于谱减法及小波分析的海洋磁测噪声抑制方法。

背景技术

随着近些年我国海洋经济的蓬勃发展，近海海域常有他国潜艇、无人潜航器等目标进行入侵、渗透、侦查和破坏等活动，对我国领土安全构成了严重威胁。含有磁性材料的目标体因其本身的永久磁性和地磁场的磁化作用，会导致物体周围的地磁场发生局部畸变，形成一种准静态磁场，即磁异常信号(Magnetic Anomaly Signatures)。相对于传统的声学信号，海洋磁异常检测具有较快的传播速率、较小的信号衰减以及浅海区域较强的抗干扰能力，常用作水下入侵目标等的探测、定位和追踪，可为我国海域内重要设施提供安全警戒和防御保障，是近海地区可靠的探测手段。

磁异常目标探测可通过机载平台(有人机、无人机)和水下磁场传感器网络实现。由于平台传感器选择的类型不同，输出信号特征会有很大差异。如机载平台测得磁异常信号所在频段为0.04～0.6Hz，而水下磁异常信号一般集中在0.0015～0.5Hz。大量的实验分析表明海水运动噪声与目标信号在量级、频带上非常接近，尤其在飞行器执行低空磁探飞行时，海浪噪声是主要的干扰源，因此对海浪磁场噪声进行抑制和处理很有必要。

七、八十年代美国海军研究了机载平台的水下目标磁异常信号的频带范围，验证了正交基匹配滤波方法去噪的可行性，推导了探测概率与地磁场的关系。以色列学者Ginzburg B等研究了基于自回归模型的白化滤波方法，能有效地抑制地磁背景噪声，提高检测率。OTNES R等人利用参考磁场传感器形成磁场梯度，在海洋背景下噪声抑制水平可达到20 dB。加拿大国防部等相关机构在海底固定站利用水压传感器去除海浪磁场噪声，背景噪声可降低10～20dB。瑞典学者Richardsen SK等人通过海上实验测得无磁异常目标时水下背景磁场数据，分别采用参考滤波法和线性预测法对地磁场噪声的抑制效果进行了对比。前者噪声抑制能力和参考测量点的位置有关，而后者无需传感器之间进行数据交换。有研究人员为解决系统采集时钟偏差问题，利用两个实验点信号进行谱减抑制背景磁场噪声，但由于存在参考噪声测量点选择问题，处理后的目标信号会引入新的干扰噪声，系统探测虚警概率增大。

公开内容

(一)要解决的技术问题

本公开提供了一种基于谱减法及小波分析的海洋磁测噪声抑制方法，以至少部分解决以上所提出的技术问题。

(二)技术方案

根据本公开的一个方面，提供了一种基于谱减法及小波分析的海洋磁测噪声抑制方法，包括：步骤S1，建立海浪磁场模型；步骤S2，对带有海浪磁场噪声的磁异常检测信号进行的谱减去噪；步骤S3，确定分解层次，对处理后的信号进行小波分析。

在本公开一些实施例中，所述步骤S2包括：步骤S201，得到带噪声的磁异常信号；步骤S202，利用海浪磁场模型估算海浪噪声，对带有海浪磁场噪声的磁异常检测信号进行的谱减去噪。

在本公开一些实施例中，所述步骤S2中，带噪磁异常信号为：

y(n)＝s(n)+e(n)0≤n≤N-1

其中，s(n)为纯净的磁异常信号，e(n)为磁异常检测中的各种噪声信号的和，假设Y(k)，S(k)，E(k)分别为y(n)，s(n)，e(n)的傅里叶变换，

Y(k)＝|Y_k|exp(jθ_k)

S(k)＝|S_k|exp(jγ_k)

其中，θ_k为Y(k)的相位，γ_k为S(k)的相位。

在本公开一些实施例中，所述步骤S202中，磁异常检测中的噪声e(n)为海浪感生的磁场信号，利用海浪磁场模型估算海浪噪声，得到E(k)。

在本公开一些实施例中，令x轴为海水水平前进方向，z轴竖直向上，x，y，z满足右手螺旋关系，建立坐标系，海浪运动产生的磁场信号表示为：

H＝h·exp(iωt-jkx)

其中，ω为海浪频率，k为波数，h为一待定系数且满足h＝h_xi+h_yj+h_zk，

求解关于h_z的二元齐次微分方程，得到

其中：

N＝cosIcosθ，

R＝sinI，

其中，T为海浪周期，σ是海水电导率，μ是磁导率，H为海浪高度，d为海洋深度，I为磁倾角，θ为海浪前进方向与地磁北极的夹角。

在本公开一些实施例中，当d/T＞0.5或d/gT²＞0.0792时，忽略海底对海浪磁场的影响，得到深海条件下：

其中，γ＝σω，F是地磁场强度的幅值，i，j，k分别为3个坐标轴的单位矢量。

在本公开一些实施例中，y(n)的功率谱为：

|Y²(k)|＝|S²(k)|+|E²(k)|+S(k)E^*(k)+S^*(k)E(k)

通过对海浪噪声E(k)进行估算，并根据估计得到的噪声功率谱得到

忽略海浪磁场模型相位的影响，运用Y(k)的相位θ_k代替S(k)的相位γ_k，再进行傅里叶反变换有

在本公开一些实施例中，所述步骤S3包括以下子步骤：子步骤S301，确定分解层次，对处理后的信号s进行小波j层分解；子步骤S302，对第一层到第j层的高频系数进行压缩，去除高频噪声；子步骤S303，重构信号，得到去噪后的磁异常信号。

在本公开一些实施例中，所述信号s分解的过程中，利用高通滤波器H1和低通滤波器H0对信号s进行第一重分解，得到高频信号cD1和低频信号cA1，其频带宽度分别占信号s的一半，再利用高通滤波器H1和低通滤波器H0对cA1进行高频和低频滤波分解，同样分解的信号各占cA1一半的频带宽度，依此类推，直到达到设定的分解层数；在重构信号过程中，分解后的信号经过去噪处理后逆序重新组合成处理后的原始信号s。

在本公开一些实施例中，通过j取不同整数值即对应不同的尺度实现对信号的多尺度分解，待分解的信号满足：S＝cA2+cD2+cD1，其中cA2继续分解，记为

其中，d_j为细节系数，c_j为尺度系数，H₀为低通滤波器，H₁为高通滤波器；

待达到分解层数后，将最终得到的各高频信号cD_j，j＝1，2，3...的系数进行压缩，达到去噪的目的；

当对信号进行重构时，有

其中滤波器和分别为H₀和H₁的逆序。

(三)有益效果

从上述技术方案可以看出，本公开基于谱减法及小波分析的海洋磁测噪声抑制方法至少具有以下有益效果其中之一：

(1)相比于采用参考测量点进行谱减法抑制海浪噪声，利用模型仿真海浪磁场的时域和频域信息不会引入外部因素干扰，在确保模型准确的前提下该方法更具有普适性；

(2)在信噪比的提升方面，有限深模型的处理结果信噪比可由-10dB提升至1dB，而无限深海浪磁场模型的处理结果信噪比提升至0.8dB。且在海浪的集中频段0.4-0.8Hz范围内海浪分布明显减少；

(3)通过利用谱减法和小波分析结合的方法，对两种不同深度的海浪磁场模型分别进行分析，都能够较大幅度地改善时域波形，抑制海浪噪声。

附图说明

图1为本公开第一实施例基于谱减法及小波分析的海洋磁测噪声抑制方法的流程图。

图2为本公开实施例磁性目标运动检测模型示意图。

图3为本公开实施例海浪运动产生感应磁场的数学模型示意图。

图4为本公开实施例谱减法的基本原理框图。

图5为本公开实施例小波分析分解与重构原理框图。

图6为本公开实施例实验地点(十字标注原点)以及海面浪高示意图。

图7为本公开实施例磁偶极子产生的z方向磁场信号的示意图。

图8为本公开实施例带有磁异常信号的海洋磁场信息的时域及频域显示的示意图。

图9为本公开实施例海浪运动磁场模型的时域及频域表示的示意图。

图10为本公开实施例经过谱减法及小波分析得到的磁异常信号检测结果。

具体实施方式

本公开提供了一种基于谱减法及小波分析的海洋磁测噪声抑制方法，在改进型的海浪磁场模型基础上，利用谱减法和小波分析相结合的方法对海浪磁场噪声进行抑制，并对比了应用不同海浪磁场模型时的去噪效果。通过实验观测的带噪磁异常目标数据验证了该方法的有效性，信噪比可提升约11dB。该方法计算复杂度低，实时性强且易于实现，可为海洋磁异常检测提供一种更加有效的的噪声抑制手段。

为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。

本公开某些实施例于后方将参照所附附图做更全面性地描述，其中一些但并非全部的实施例将被示出。实际上，本公开的各种实施例可以许多不同形式实现，而不应被解释为限于此数所阐述的实施例；相对地，提供这些实施例使得本公开满足适用的法律要求。

在本公开的第一个示例性实施例中，提供了一种基于谱减法及小波分析的海洋磁测噪声抑制方法。图1为本公开第一实施例基于谱减法及小波分析的海洋磁测噪声抑制方法的流程图。如图1所示，本公开基于谱减法及小波分析的海洋磁测噪声抑制方法包括以下步骤：

步骤S1，建立海浪磁场模型；

步骤S2，对带有海浪磁场噪声的磁异常检测信号进行的谱减去噪；

步骤S3，确定分解层次，对处理后的信号进行小波分析。

以下分别对本实施例基于谱减法及小波分析的海洋磁测噪声抑制方法的各个步骤进行详细描述。

为了验证算法的有效性，可以通过建立磁异常目标模型进行仿真。磁异常目标在地磁场中产生的畸变可以看作具有一定磁距的独立的磁场源作用的结果。在工程实践中，探测距离基本能达到目标体最大尺寸的三倍以上。此时，磁异常体可看作磁偶极子，其矢量表达式如下：

其中，μ₀为自由空间磁导率，m为目标总体磁矩矢量，r为目标到测量点处的距离矢量。假设矢量磁场传感器位于坐标轴零点位置，其三方向指向分别与坐标轴x，y，z方向平行。图2为本公开实施例磁性目标运动检测模型示意图。如图2所示，磁异常目标沿着x轴正方向以速度v匀速直线经过磁场传感器，磁矩大小和方向始终保持不变。

假定目标本身尺寸可忽略，磁场传感器在零点处测得磁异常信号可以表示为：

式中，B_x、B_y、B_z分别为目标产生的磁异常信号在坐标轴三个方向上的投影，m_k，r_k(k＝x，y，z)分别表示磁性目标磁距和目标与磁场传感器距离在x，y，z方向上的三个分量，磁性目标运动轨迹为时间的函数(c₁，c₂为常数)，其中磁异常目标沿着x轴正方向匀速直线经过磁场传感器的速度为v。

在海洋磁异常信号检测过程中时，传感器不可避免的会接收到来自于环境的噪声，如海浪、海流、岩石圈等磁场信号。为了对带噪的观测信号进行处理，需要建立海浪磁场模型对海浪噪声进行估计。图3为本公开实施例海浪运动产生感应磁场的数学模型示意图。令x轴为海水水平前进方向，z轴竖直向上，x，y，z满足右手螺旋关系，建立如图3所示的坐标系。

其中，H为海浪高度，d为海洋深度，I为磁倾角，θ为海浪前进方向与地磁北极的夹角，u为海浪的运动速度矢量，F＝F(cosI cos θi+cosIsin θj-sin Ik)为地磁场的方向矢量(i，j，k分别为3个坐标轴的单位矢量)。空气和海水的介电常数分别为ε₁＝ε₀＝10^-9/36πF/m，ε₂＝ε₀ε_r(ε_r＝78)，ω为海浪频率，k为波数。海浪运动产生的磁场信号可以表示为：

H＝h·exp(iωt-jkx) (3)

h为一待定系数且满足h＝h_xi+h_yj+h_zk，求解关于h_z的二元齐次微分方程，得到

其中：

N＝cos I cosθ，

R＝sin I，

T为海浪周期，σ是海水电导率，μ是磁导率。

当d/T＞0.5或d/gT²＞0.0792时，海底对海浪磁场的影响可忽略，可得到深海条件下：

其中，γ＝σω，

实际观测到的海浪并不是标准的正弦波形式，但可以将其看成不同频率正弦波的叠加。假设各个频率之间彼此不相关，可以得到海浪磁场表达式为：

步骤S2中对带有海浪磁场噪声的磁异常检测信号进行的谱减去噪时，具体包括以下步骤。

步骤S201，得到带噪声的磁异常信号；

假设带噪磁异常信号为：

y(n)＝s(n)+e(n)0≤n≤N-1 (7)

其中，s(n)为纯净的磁异常信号，e(n)为磁异常检测中的各种噪声信号的和，该处主要认为是海浪感生的磁场。假设Y(k)，S(k)，E(k)分别为y(n)，s(n)，e(n)的傅里叶变换，其中

其中，θ_k为Y(k)的相位，γ_k为S(k)的相位；

步骤S202，利用海浪磁场模型估算海浪噪声，对带有海浪磁场噪声的磁异常检测信号进行的谱减去噪。

谱减法实时性强，易于实现，且算法复杂度低。它的基本思想是假定加性高斯噪声与短时平稳信号彼此相互独立，将带噪信号的功率谱与已知的噪声功率谱相减，从而得到较为纯净的目标信号频谱。由于海浪磁场信号与磁异常信号在时域和频域上都有重叠，滤除难度较大。但是由于两者不相关，且在模型仿真中海浪的相位信息可以忽略不计，因此可以考虑采用谱减法进行处理。

谱减法的原理如图4所示。如图4所示，带噪的磁异常信号经过傅里叶变换得到Y(k)，将其幅值部分|Y(k)|进行平方，计算能量，相位部分θ(k)保留，然后用上文中的海浪磁场模型仿真得到的时域信号做傅里叶变换，估算E(k)，得到|E(k)|²，经过与|Y(k)|²相减后得到纯净的磁异常信号|S(k)|²，再通过开方得到纯净磁异常信号的傅氏变换的幅值部分，与保留的θ(k)组合再通过反变化即得到纯净的时域磁异常信号。

由于实际环境中海浪信号与磁异常信号频带相近，且占主要部分。因此假设磁异常检测中的噪声e(n)为海浪感生的磁场信号，利用海浪磁场模型估算海浪噪声。具体地，E(k)可以通过式(4)(5)(6)计算得到。

根据式(7)，y(n)的功率谱为

|Y²(k)|＝|S²(k)|+|E²(k)|+S(k)E^*(k)+S^*(k)E(k)>

由于海浪噪声与磁异常信号不相关，通过对海浪噪声E(k)进行估算，并根据估计得到的噪声功率谱得到

忽略海浪磁场模型相位的影响，运用Y(k)的相位θ_k代替S(k)的相位γ_k，再进行傅里叶反变换有

步骤S3中采用小波分析对谱减去噪后的信号进行处理。频谱相减法中傅里叶变换采用的是固定窗函数，时频分辨率不变，信号的细节特征相对模糊。而小波变换通过伸缩平移运算，将信号在多个尺度上进行子波分解，实现高频处时间细分，低频处频率细分，所得各尺度子波变换系数代表原信号在不同分辨率上的信息。其处理的实质在于减小甚至完全剔除由噪声产生的系数，同时最大限度地保留真实信号的系数，最后由经过处理的小波系数重构原信号，得到真实信号的最优估计。设信号为s(t)，则离散小波变换(记作DWT，即Discrete Wavelet Transform)表示如下：

其中，ψ(t)∈L²(R)为基小波函数，a和b分别为尺度参数和平移参数。其重构公式为

所述步骤S3具体包括以下子步骤：

子步骤S301，确定分解层次，对处理后的信号进行小波j层分解。

子步骤S302，对第一层到第j层的高频系数进行压缩，去除高频噪声；

子步骤S303，重构语音信号，得到去噪后的磁异常信号。

小波分解与重构的原理框图如图5所示。如图5所示，上图中为信号s分解的过程，利用高通滤波器H1和低通滤波器H0对信号进行第一重分解，得到高频信号cD1和低频信号cA1，其频带宽度分别占信号s的一半，再对cA1进行高频和低频滤波分解，同样分解的信号各占cA1一半的频带宽度，依此类推，直到达到设定的分解层数；下图为信号重建的过程，将上图中分解后的信号经过去噪处理后逆序重新组合成处理后的原始信号s。

从分解框图中可以看出，当j取不同整数值即对应不同的尺度，可以实现对信号的多尺度分解。待分解的信号满足：S＝cA2+cD2+cD1，其中cA2可以继续分解，分解公式可以记为

其中，d_j为细节系数，c_j为尺度系数。其中，低通滤波器H₀的作用是实现对信号的逼近，高通滤波器H₁的作用是提取信号的细节信息。待达到分解层数后，将最终得到的各高频信号cDj(j＝1，2，3...)的系数进行压缩，达到去噪的目的。

当对信号进行重构时，有

其中滤波器和分别为H₀和H₁的逆序。该过程即是分解过程的逆序。

为了验证上述算法的有效性，在对我国南海海域的海洋磁异常信号进行了实验采集。实验中所用的三分量海底磁测仪噪声水平约为10pT/sqrt(Hz)@1Hz，采样频率为500Hz。图6为本实施例实验地点(十字标注原点)以及海面浪高示意图，展示了实验地点的海浪情况。

因本次实验中磁测仪水平方向未采用角度计进行校正，因此选用z轴数据进行分析。图7为本公开实施例磁偶极子产生的z方向磁场信号的示意图，该磁异常信号即是根据(2)式模型仿真得到的，根据式(2)计算可以得到纯净的z方向磁异常信号，其中c₁，c₂为0，磁异常目标运动速度为0.1m/s，磁距为(10000000，0，0)，在t＝0时目标位于(-1000，100，200)的位置，产生的磁异常信号如图7所示。将仿真的磁异常信号加入实验观测数据，得到半仿真的磁异常信号。

图8分别展示了在南海区域观测到的磁异常信号的时域和频域信息，从图中并未能观测到明显的磁异常信号，主要表现为海浪信号及地磁信号，在信号频谱中可以很明显的看到海浪信号波峰，中心频率为f＝0.06886Hz，磁异常信号与海浪磁场的信噪比约为-10dB。

考虑到实验进行过程中天气、海面状况及地磁场特性参数，在利用海浪磁场模型估算海浪噪声时，令F＝43333.7nT，θ＝-1°27′6″，I＝23°47′7″。海浪频率在0.04-0.9Hz范围内，海浪高度H＝1m，根据式(6)分别计算海底有限深和无限深条件下海浪磁场的时域及频域仿真结果如图9所示。

由于此次实验采样率为500Hz，选取小波分解层数J＝10。图10为本公开实施例经过谱减法及小波分析得到的磁异常信号检测结果，如图10所示，为通过谱减法及小波分解分别得到有限深及无限深条件下的磁异常信号检测结果。

至此，本公开第一实施例基于谱减法及小波分析的海洋磁测噪声抑制方法介绍完毕。

至此，已经结合附图对本公开实施例进行了详细描述。需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。

并且图中各部件的形状和尺寸不反映真实大小和比例，而仅示意本公开实施例的内容。另外，在权利要求中，不应将位于括号之间的任何参考符号构造成对权利要求的限制。

此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

在此提供的算法和显示不与任何特定计算机、虚拟系统或者其它设备固有相关。各种通用系统也可以与基于在此的示教一起使用。根据上面的描述，构造这类系统所要求的结构是显而易见的。此外，本公开也不针对任何特定编程语言。应当明白，可以利用各种编程语言实现在此描述的本公开的内容，并且上面对特定语言所做的描述是为了披露本公开的最佳实施方式。

本公开可以借助于包括有若干不同元件的硬件以及借助于适当编程的计算机来实现。本公开的各个部件实施例可以以硬件实现，或者以在一个或者多个处理器上运行的软件模块实现，或者以它们的组合实现。本领域的技术人员应当理解，可以在实践中使用微处理器或者数字信号处理器(DSP)来实现根据本公开实施例的相关设备中的一些或者全部部件的一些或者全部功能。本公开还可以实现为用于执行这里所描述的方法的一部分或者全部的设备或者装置程序(例如，计算机程序和计算机程序产品)。这样的实现本公开的程序可以存储在计算机可读介质上，或者可以具有一个或者多个信号的形式。这样的信号可以从因特网网站上下载得到，或者在载体信号上提供，或者以任何其他形式提供。

本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。并且，在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。

类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

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