基于动态卡尔曼滤波的极低频海洋电磁信号探测方法

1.本发明涉及信号检测技术领域，特别涉及基于动态卡尔曼滤波的极低频海洋电磁信号探测方法。


背景技术：

2.为抵御海水的腐蚀，舰艇通常都会加装阴极保护措施，由此产生的保护电流被螺旋桨调制而形成以螺旋桨的旋转频率为基频的时谐交变电流，继而形成具有显著谐波特征的目标极低频电磁场。作为一种舰艇目标的固有物理属性，极低频电磁场在海水中具有良好的穿透性、传播衰减较慢，并且携带有明显的目标运动信号特征，可作为海洋目标远距离探测、识别的重要依据。
3.虽然极低频电磁信号传播衰减较慢，但是在实际海洋环境下包含各种类型的电磁噪声干扰，致使所获实测数据的信噪比严重降低，给极低频信号的远距离探测带来了很大的困难。作为对极低频电磁信号探测影响较大的主要噪声源，海浪感应噪声和低频随机噪声都具有时变和随机特点，其产生机理和信号域特征各不相同。
4.鉴于海浪感应噪声和低频随机噪声对极低频电磁信号探测带来的影响，需要从信号处理角度来抑制噪声以提升数据的信噪比，为目标极低频电磁信号的远距离探测提供行之有效的方法。已有海洋电磁噪声的抑制技术大多基于噪声和信号在频域的特征差异来实现对噪声的过滤，在很多问题中取得了较为理想的信号探测效果。然而，现有方法的前提是基于噪声和信号在频域特征的可分性，很难直接用于对极低频信号的探测问题。海浪感应噪声的频谱为宽带谱，与低频随机噪声的频谱类似，都覆盖了目标信号所在的极低频段，通过频域滤波技术很难将目标信号和噪声有效分离。因此，有必要针对极低频电磁信号和海洋环境噪声的特点提出有效的去噪算法，提升在低信噪比条件下对极低频电磁信号的探测能力。


技术实现要素：

5.本发明实施例提供了基于动态卡尔曼滤波的极低频海洋电磁信号探测方法，用以解决现有去噪技术基于噪声和信号在频域特征的可分性，很难直接用于对极低频信号的探测问题。
6.一方面，本发明实施例提供了基于动态卡尔曼滤波的极低频海洋电磁信号探测方法，包括：
7.获取含噪极低频海洋电磁信号时间序列；
8.将含噪极低频海洋电磁信号时间序列输入动态卡尔曼滤波器，获得滤波后的电磁信号时间序列；
9.对滤波后的电磁信号时间序列进行短时傅里叶变换，获得相应的去噪时频谱；
10.提取去噪时频谱中极低频海洋电磁信号的特征参数；
11.其中，动态卡尔曼滤波器对含噪极低频海洋电磁信号时间序列的滤波处理包括：
12.利用动态卡尔曼滤波器对含噪极低频海洋电磁信号时间序列进行系统状态和协方差更新；
13.根据更新后的系统状态和协方差确定当前状态的最优估计值；
14.更新与最优估计值对应的协方差。
15.本发明中的基于动态卡尔曼滤波的极低频海洋电磁信号探测方法，具有以下优点：
16.1、本发明联合信号的时域和时频域分析方法，对海洋电磁噪声的特性进行了研究，并挖掘了海洋电磁噪声和极低频电磁信号的可分离属性，基于此提出了基于卡尔曼模型的电磁噪声动态滤波算法，为在时频特征混叠情况下的微弱信号增强探测问题提供了新思路。
17.2、本发明所提出的滤波算法具有较强的鲁棒性，能够根据环境的变化动态调整和优化滤波模型参数，与传统的时频域滤波技术相比，对信号和噪声先验知识的依赖程度较低，适用于具有时变和随机海洋电磁干扰情况的远距离目标探测场景。
18.3、本发明在实测海洋电磁数据上进行了验证，所提方法能够在负信噪比条件下成功提取出感兴趣的目标极低频信号，同时设计的迭代滤波算法计算效率较高，能够实现对海洋电磁数据的实时处理，在实际探测问题中具有一定的应用价值和潜力。
附图说明
19.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
20.图1为本发明实施例提供的基于动态卡尔曼滤波的极低频海洋电磁信号探测方法的流程图；
21.图2为本发明实施例提供的实测海洋含噪电磁信号时间序列；
22.图3为本发明实施例提供的实测海洋含噪电磁信号频谱图；
23.图4为本发明实施例提供的实测海洋含噪电磁信号时频图；
24.图5为本发明实施例提供的时域冲击信号的示意图；
25.图6为本发明实施例提供的与时域冲击信号对应的冲击时频图；
26.图7为本发明实施例提供的实验一中海洋含噪电磁信号的时间序列；
27.图8为本发明实施例提供的实验一滤波后的电磁信号的时间序列；
28.图9为本发明实施例提供的实验二中海洋含噪电磁信号时频图；
29.图10为本发明实施例提供的实验二滤波后的电磁信号时频图。
具体实施方式
30.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
31.图1为本发明实施例提供的基于动态卡尔曼滤波的极低频海洋电磁信号探测方法的流程图。本发明实施例提供了基于动态卡尔曼滤波的极低频海洋电磁信号探测方法，该方法包括：
32.s100，获取含噪极低频海洋电磁信号时间序列。
33.示例性地，可将获取到的含噪极低频海洋电磁信号时间序列记为x(m)，如图2所示。
34.s120，将含噪极低频海洋电磁信号时间序列输入动态卡尔曼滤波器，获得滤波后的电磁信号时间序列。
35.示例性地，对于包含微弱极低频电磁信号的时间序列x(m)，将其按采样顺序依次输入动态卡尔曼滤波器中，经过迭代运算得到滤波后的电磁信号时间序列
36.s140，对滤波后的电磁信号时间序列进行短时傅里叶变换，获得相应的去噪时频谱。
37.示例性地，对滤波后的电磁信号时间序列做短时傅里叶变换(short time fourier transform，stft)，得到相应的去噪时频谱
[0038][0039]
其中，w(m)是窗函数，n为其样本偏移量，m是含噪极低频海洋电磁信号时间序列x(m)的长度。根据去噪时频谱绘制的去噪时频图如图4所示。
[0040]
s160，提取去噪时频谱中极低频海洋电磁信号的特征参数。
[0041]
示例性地，基于去噪时频谱借助图像可视化分析方法提取经增强后的极低频电磁信号的频率以及其在所观测时间序列中出现的起始时间ts和截止时间te，以作为对海洋极低频电磁信号的探测结果，继而为目标信息的准确感知提供依据。
[0042]
其中，s120具体包括：
[0043]
s121，利用动态卡尔曼滤波器对含噪极低频海洋电磁信号时间序列进行系统状态和协方差更新。
[0044]
示例性地，首先建立动态卡尔曼滤波器的状态方程x(k)和观测方程z(k)：
[0045]
x(k)＝ax(k-1) bu(k) w(k)
[0046]
z(k)＝hx(k) v(k)
[0047]
其中，x(k)、u(k)和w(k)是k时刻系统的状态、控制量和过程噪声，z(k)和v(k)是k时刻系统的观测值和测量噪声，a和b是系统参数，h是测量系统的参数。在本发明的实施例中，w(k)和v(k)被假设成高斯白噪声，对应的协方差q和r不随系统状态变化而变化。
[0048]
在滤波处理前，还需要设置观测和测量过程中的噪声协方差参数q和r的值，并初始化系统状态x(1|1)和相应的协方差p(1|1)。为减小滤波后起始过程的起伏波动，初始系统状态x(1|1)选取为电磁信号时间序列的均值。同时为保证迭代算法的收敛性，滤波器初始状态的协方差p(1|1)取不为0的数。
[0049]
接着对时刻k进行迭代。令k＝2，利用建立的动态卡尔曼滤波器，基于系统上一时
刻的状态，即x(1|1)预测当前时刻的状态x(k|k-1)和相应的协方差p(k|k-1)，获得对系统当前状态的预测结果：
[0050]
x(k|k-1)＝ax(k-1|k-1) bu(k)
[0051]
其中，x(k|k-1)是利用上一时刻状态预测的结果，x(k-1|k-1)是上一时刻状态最优估计值。
[0052]
然后更新x(k|k-1)的协方差p(k|k-1)的值：
[0053]
p(k|k-1)＝ap(k-1|k-1)a
t
q
[0054]
其中，p(k-1|k-1)是x(k-1|k-1)对应的协方差，a
t
表示a的转置矩阵。
[0055]
s122，根据更新后的系统状态和协方差确定当前状态的最优估计值。
[0056]
示例性地，结合当前时刻的状态预测结果和电磁传感器获得的实际测量结果z(k)，计算当前时刻下的最优化估算值x(k|k)：
[0057]
x(k|k)＝x(k|k-1) kg(k)[z(k)-hx(k|k-1)]
[0058]
其中，当前时刻的卡尔曼增益项kg(k)按下式计算：
[0059]
kg(k)＝p(k|k-1)h
t
/[hp(k|k-1)h
t
r]
[0060]
s123，更新与最优估计值对应的协方差。
[0061]
示例性地，更新当前时刻下最优估算值x(k|k)的协方差p(k|k)，为下一轮迭代运算做准备：
[0062]
p(k|k)＝[i-kg(k)h]p(k|k-1)
[0063]
其中，i为单位矩阵。
[0064]
在完成以上更新处理后，令k＝k 1，循环执行以上步骤的操作，直到k大于含噪极低频海洋电磁信号时间序列x(m)的长度m时为止，动态卡尔曼滤波器的迭代过程结束。
[0065]
在一种可能的实施例中，在s160之后，方法还包括：
[0066]
s180，对含噪极低频海洋电磁信号时间序列进行快速傅里叶变换和短时傅里叶变换，获得相应的含噪频谱含噪时频谱。
[0067]
示例性地，对于图2所示的电磁信号时间序列x(m)，做快速傅里叶变换(fast fouriertransform，fft)得到其含噪频谱x(k)：
[0068][0069]
根据含噪频谱计算结果x(k)绘制的含噪频谱图如图3所示。
[0070]
s200，结合含噪频谱和含噪时频谱，确定海浪感应噪声的出现时刻以及低频随机噪声的能量频段范围。
[0071]
示例性地，结合以上含噪频谱图和含噪时频图，对观测数据段内海浪感应噪声的时频关联特性进行分析，并估计出各海浪感应噪声干扰的出现时刻ti,i＝1,2,
…
,n。同时，结合含噪频谱和含噪时频谱，对海洋低频随机噪声的时频分布特性进行分析，并估计其能量存在的主要频段范围[f
min
,f
max
]。
[0072]
s220，建立时域冲击信号。
[0073]
示例性地，仿真建立如图5所示的时域冲击信号h(t)＝5δ(t-4)-5δ(t-6)。
[0074]
s240，对时域冲击信号进行短时傅里叶变换，获得相应的冲击时频谱。
[0075]
示例性地，对时域冲击信号做短时傅里叶变换，得到冲击时频谱h(t,ω)：
[0076]
h(t,ω)＝∫h(τ)g(τ-t)e-jωτ
dτ
[0077]
其中，g(τ)为窗函数，τ为积分变量，t和ω分别为移位变量和频率变量。根据冲击时频谱计算结果h(t,ω)绘制冲击时频图如图6所示。
[0078]
s260，根据冲击时频谱，确定海浪感应噪声的幅度。
[0079]
示例性地，按照所建立的时域冲击信号和冲击时频谱，估计含噪极低频海洋电磁信号时间序列x(m)中出现的各海浪感应噪声的幅度ai,i＝1,2,
…
,n。
[0080]
本发明实施例还提供了极低频海洋电磁信号探测系统，包括：
[0081]
序列获取模块，用于获取含噪极低频海洋电磁信号时间序列；
[0082]
滤波模块，用于将含噪极低频海洋电磁信号时间序列输入动态卡尔曼滤波器，获得滤波后的电磁信号时间序列；
[0083]
短时傅里叶变换模块，用于对滤波后的电磁信号时间序列进行短时傅里叶变换，获得相应的去噪时频谱；
[0084]
参数提取模块，用于提取去噪时频谱中极低频海洋电磁信号的特征参数；
[0085]
其中，滤波模块包括：
[0086]
系统更新子模块，用于利用动态卡尔曼滤波器对含噪极低频海洋电磁信号时间序列进行系统状态和协方差更新；
[0087]
最优估计子模块，用于根据更新后的系统状态和协方差确定当前状态的最优估计值；
[0088]
协方差更新子模块，用于更新与最优估计值对应的协方差。
[0089]
在一种可能的实施例中，还包括：信号转换模块，用于含噪极低频海洋电磁信号时间序列进行快速傅里叶变换和短时傅里叶变换，获得相应的含噪频谱和含噪时频谱；噪声确定模块，用于结合含噪频谱和含噪时频谱，确定海浪感应噪声的出现时刻以及低频随机噪声的能量频段范围；信号建立模块，用于建立时域冲击信号；冲击变换模块，用于对时域冲击信号进行短时傅里叶变换，获得相应的冲击时频谱；幅度确定模块，用于根据冲击时频谱，确定海浪感应噪声的幅度。
[0090]
本发明实施例还提供了一种计算机可读存储介质，该计算机可读存储介质中存储有多条计算机指令，多条计算机指令用于使计算机执行上述的方法。
[0091]
本发明的效果通过以下仿真实验进一步说明：
[0092]
(1)实验仿真条件：
[0093]
本实验采用的数据是典型的海洋电磁信号时间序列，该数据是2018年8月在山东省青岛市冯家河码头由感应式电磁传感器采集，数据总时长1h、采样率为200hz，包含了海浪感应噪声和低频随机噪声，以及感兴趣的目标极低频谐波信号。极低频谐波信号所在基频段被低频随机噪声干扰严重，同时伴随海浪出现的感应噪声也会将极低频谐波信号淹没，致使所获电磁观测数据的整体信噪比偏低，对于极低频信号的探测具有一定的挑战性。本实验在主频为3.20ghz的core
tm i7-8700 cpu、ram内存为16gb的64位windows 10操作系统上，基于matlab r2013b软件开展。
[0094]
(2)噪声抑制评价标准：
[0095]
(2a)时域评价标准：峰峰值下降比η
[0096]
峰峰值下降比的定义式如下：
[0097][0098]
其中，为滤波前噪声时域峰峰值，为滤波后噪声时域峰峰值。η的值越大，滤波后噪声峰峰值较滤波前下降越多，表明在时域对噪声的抑制效果越好。
[0099]
(2b)频域评价标准：噪声频带上限δ
[0100]
噪声频带上限的定义式如下：
[0101][0102]
其中，h(f)为含噪数据频谱，f
up
为频率上限变量。上式通过寻找满足括号内不等式的频率上限最小值，将其定义为噪声频带上限δ。噪声频带上限δ越大，表示低频随机海洋噪声的频谱能量分布范围越广，对于极低频信号频段的干扰也越严重。
[0103]
(3)实验内容：
[0104]
实验一
[0105]
在系统参数a和b取值为1和0、测量系统参数h取值为1，观测过程和测量过程噪声的协方差q和r分别取值为1
×
10-6
和4
×
10-4
，滤波器初始状态x(1|1)设置为0.0406nt、初始状态的协方差p(1|1)设置为1的条件下，利用本发明设计的动态卡尔曼滤波器对极低频海洋电磁信号时间序列进行滤波，实验结果如图7-8所示，其中图7是含噪电磁信号时间序列，图8是滤波后的电磁信号时间序列。
[0106]
从图7-8可见，本发明通过结合历史状态来动态校正测量结果，可以取得较好的噪声时域能量抑制效果。
[0107]
分别计算含噪电磁数据和滤波后电磁数据在时域的峰峰值以及峰峰值下降比，计算结果如表1所示。
[0108]
表1含噪和滤波后电磁数据时域峰峰值及峰峰值下降比
[0109]
指标含噪数据峰峰值滤波数据峰峰值峰峰值下降量峰峰值下降比数值53.0022nt38.0096nt14.9926nt0.2829
[0110]
从表1可见，滤波后电磁数据的峰峰值较含噪电磁数据下降了14.9926nt，峰峰值下降比达0.2829。本发明能够通过动态噪声过滤算法，对海洋噪声在时域的能量进行有效抑制，为海浪感应噪声干扰下的极低频信号探测提供支撑。
[0111]
实验二
[0112]
对含噪电磁数据做短时傅里叶变换得到其含噪时频图，并与经过本发明方法滤波后的电磁数据的去噪时频图结果进行对比，实验结果如图9-10所示，其中图9是含噪电磁数据的含噪时频图，图10是经本发明方法滤波后电磁数据的去噪时频图。
[0113]
从图9-10可见，本发明可以有效抑制低频随机噪声，并顺利探测出被噪声淹没的3hz微弱极低频基频信号，实现了在低信噪比条件下对目标信号的准确感知。
[0114]
对本发明方法滤波后数据和含噪电磁数据的噪声频带上限进行对比，数值指标结果如表2所示。
[0115]
表2本发明方法滤波后数据和含噪电磁数据的噪声频带上限指标对比
[0116]
待评价数据本发明方法滤波后数据含噪电磁数据噪声频带上限δ2hz5hz
[0117]
从表2可见，本发明方法滤波后数据的噪声频带上限为2hz，较原始含噪电磁数据下降了3hz。本发明能够通过动态噪声过滤方法，对海洋噪声在频域的能量进行有效抑制，为低频随机噪声干扰下的极低频信号探测提供支撑。
[0118]
综上，本发明通过将历史状态信息引入滤波器之中，可以有效抑制包括海浪感应噪声和低频随机噪声在内的典型噪声干扰，极大地提升了在极低频段海洋电磁数据的整体信噪比，并基于时频域特征关联分析对目标极低频信号基频实现了增强探测，在海洋遥感信息感知领域具有广阔的应用前景和价值。
[0119]
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0120]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。