一种基于改进的经验模态分解的地磁测量信号去噪方法与流程

本发明涉及地磁测量信号干扰抑制领域，尤其是一种基于改进的经验模态分解的地磁测量信号去噪方法。

背景技术：

在地磁测量中，由于载体上的电子设备、惯性导航设备以及载体周围环境的随机干扰普遍存在且无法有效屏蔽，需通过滤波等方法去除噪声，提高测量精度。常见的卡尔曼滤波、小波变换等方法无法满足磁场测量中的非线性不确定性需求；经验模态分解可根据信号本身的时间尺度特征对信号进行分解，有效区分信号与噪声干扰。而传统的经验模态分解法在端点非极值情况下，会产生端点效应，得到的本征信号函数会严重失真失去意义，重构得到的信号也无法满足真实表达磁场真实特征，需对信号进行预处理，抑制端点效应。

常见的抑制端点效应方法有信号延拓和序列预测。信号延拓包括有：以信号端点为对称中心的偶延拓、周期延拓以及镜像延拓等。此类延拓方法均基于信号整体具有对称性(偶延拓)或周期性(周期延拓)，无法满足磁场的非线性、非周期性变化。信号的序列预测主要有基于过去信号以及当前干扰信号的自回归预测模型(AR预测模型)和基于神经网络的预测方法。这两种方法在处理非线性信号时，会使得模型阶数很高，大大增加运行时间和预测时间。

当地磁信号存在间断的跳跃性变化时，经验模态分解结果将产生模态混叠效应，无法依据时间特征尺度分离不同的模态分量。形态学滤波可以对测量信号的突变进行平滑处理，抑制模态混叠。鉴于传统形态学滤波在平滑信号的异常突变时，会抑制信号的细节信息，使得滤波结果的信噪比降低，影响后续的数据处理，故需对其进行改进。

技术实现要素：

发明目的：本发明的目的是提供一种能够解决现有技术中存在的缺陷的基于改进的经验模态分解的地磁测量信号去噪方法。

技术方案：为达到此目的，本发明采用以下技术方案：

本发明所述的基于改进的经验模态分解的地磁测量信号去噪方法，包括以下步骤：

S1：获取地磁场测量数据序列{x₁(t)}；

S2：多次进行形态学滤波处理；

S3：对步骤S2得到的信号进行经验模态分解。

进一步，所述步骤S2包括以下步骤：

S2.1：第一次腐蚀运算；

利用形态学滤波器运算中的腐蚀运算对数据进行第一步处理，去除粗大误差，获得初步的信号{x₂(t)}，处理过程如式(1)所示：

{x₂(t)}＝ξ_S[x₁(t)]＝min[x₁(t+m)-S(m)] (1)

式(1)中，x₁(t)为待处理信号，t为信号长度，S(m)为形态学滤波器的结构元素，m为结构元素长度，结构元素由磁场测量信号决定；

S2.2：第一次膨胀运算；

利用形态学滤波器运算中的膨胀运算对数据进行进一步的处理，改善信号失真情况，获得信号{x₃(t)}，处理过程如式(2)所示：

{x₃(t)}＝ζ_S[x₂(t)]＝max[x₂(t-m)+S(m)] (2)

S2.3：第二次腐蚀运算；

利用形态学滤波器运算中的腐蚀运算对数据进行进一步处理，进一步去除粗大误差，获得信号{x₄(t)}，处理过程如式(3)所示：

{x₄(t)}＝ξ_S[x₃(t)]＝min[x₃(t+m)-S(m)] (3)

S2.4：第二次膨胀运算；

利用形态学滤波器运算中的膨胀运算对数据进行进一步的处理，进一步改善信号失真情况，获得信号{x₅(t)}，处理过程如式(4)所示：

{x₅(t)}＝ζ_S[x₄(t)]＝max[x₄(t-m)+S(m)] (4)。

进一步，所述步骤S3包括以下步骤：

S3.1：信号延拓；

选择步骤S2得到的信号的左端长度为t₁的待匹配信号段{η_l(t₁)}，以及右端长度为t₂的待匹配信号段{η_r(t₂)}，其中，{η_l(t₁)}＝{η_l(1),…,η_l(i),…η_l(t₁)}，{η_r(t₂)}＝{η_r(1),…,η_r(j),…η_r(t₂)}；步骤S2得到的信号中，除掉待匹配信号段{η_l(t₁)}和待匹配信号段{η_r(t₂)}之后的信号为剩余信号；

首先，将{η_l(t₁)}向右平行移动，在移动的过程中，使{η_l(t₁)}的极大值点依次与剩余信号中的每个极大值点重合，记相应的匹配信号段分别为θ₁,θ₂…θ_i…θ_n，并通过式(5)计算各个匹配波形与待匹配信号段{η_l(t₁)}的相似距离：

其中，

式(8)中的cov[·]是协方差，σ(·)是方差；

选择P_i最小时对应的匹配信号段作为第一最佳匹配信号段，并利用第一最佳匹配信号段对待匹配信号段{η_l(t₁)}进行向左延拓，以待匹配信号段{η_l(t₁)}的极大值点为起点，将第一最佳匹配信号段极大值点左侧的信号段复制到待匹配信号段{η_l(t₁)}极大值点的左边，使得第一最佳匹配信号段的极大值点与待匹配信号段{η_l(t₁)}的极大值点重合；

然后，将{η_r(t₂)}向左平行移动，在移动的过程中，使{η_r(t₂)}的极大值点依次与剩余信号中的每个极大值点重合，记相应的匹配信号段分别为α₁,α₂…α_j…α_n，并通过式(9)计算各个匹配波形与待匹配信号段{η_r(t₂)}的相似距离：

其中，

式(12)中的cov[·]是协方差，σ(·)是方差；

选择P_j最小时对应的匹配信号段作为第二最佳匹配信号段，并利用第二最佳匹配信号段对待匹配信号段{η_r(t₂)}进行向右延拓，以待匹配信号段{η_r(t₂)}的极大值点为起点，将第二最佳匹配信号段极大值点右侧的信号段复制到待匹配信号段{η_r(t₂)}极大值点的右边，使得第二最佳匹配信号段的极大值点与待匹配信号段{η_r(t₂)}的极大值点重合；

步骤S2得到的信号经过向左延拓和向右延拓处理之后，得到的信号记为{x₆(t₃)}；

S3.2：模态分解；

包括以下步骤：

S3.21：令k＝1，p＝1；

S3.22：寻找信号{x₆(t₃)}中所有极大值点，通过三次样条插值拟合得到极大值包络线m_k+(t₃)，寻找信号{x₆(t₃)}中所有极小值点，通过三次样条插值拟合得到极小值包络线m_k-(t₃)，取均值包络线为

S3.23：将信号{x₆(t₃)}减去m_k(t₃)得到新信号判断新信号是否满足“极值点个数和过零点个数相差不大于1并且在任意点处上下包络均值为0”这个准则：如果满足，则令为第p个本征信号函数，记为IMF_p(t₃)，然后进行步骤S3.24；如果不满足，则令k＝k+1，p＝1，然后返回步骤S3.22；

S3.24：将信号{x₆(t₃)}减去IMF_p(t₃)，得到一个去掉高频信号的新信号r_p(t₃)；

S3.25：判断新信号r_p(t₃)是否还能再进行分解：如果能，则令p＝p+1，k＝1，然后返回步骤S3.22；如果不能，则记录分解的总次数为n，最终得到的不能分解的新信号为r(t₃)，信号{x₆(t₃)}如式(13)所示：

S3.26：寻找IMF₁(t₃),…,IMF_n(t₃)中的突变点：按照式(14)计算IMF₁(t₃),…,IMF_n(t₃)中相邻两个本征信号函数之间的均方差CMSE，如果q＝u时，CMSE(IMF_u(t₃),IMF_u+1(t₃))发生了突变，则将x(t₃)作为重构信号，如式(15)所示；

有益效果：本发明采用改进的经验模态分解方法，实现了对地磁测量噪声的抑制作用；为实现测量信号粗大误差的抑制，采用了多次形态学滤波方法，可大大提高信噪比；针对传统经验模态分解存在的端点效应，提出了新的信号延拓方法，实现了对端点效应的抑制作用，加强了本征信号函数的有效性；针对分解后的信号重构问题，提供了详细的区分信号噪声的规则，加强了重构函数的信噪比，抑制了随机噪声。

附图说明

图1为本发明的方法流程图；

图2为本发明的模态分解步骤的流程图。

具体实施方式

下面结合具体实施方式对本发明的技术方案作进一步的介绍。

本发明公开了一种基于改进的经验模态分解的地磁测量信号去噪方法，如图1所示，包括以下步骤：

S1：获取地磁场测量数据序列{x₁(t)}；

S2：多次进行形态学滤波处理；

S3：对步骤S2得到的信号进行经验模态分解。

其中，步骤S2包括以下步骤：

S2.1：第一次腐蚀运算；

利用形态学滤波器运算中的腐蚀运算对数据进行第一步处理，去除粗大误差，获得初步的信号{x₂(t)}，处理过程如式(1)所示：

{x₂(t)}＝ξ_S[x₁(t)]＝min[x₁(t+m)-S(m)] (1)

式(1)中，x₁(t)为待处理信号，t为信号长度，S(m)为形态学滤波器的结构元素，m为结构元素长度，结构元素由磁场测量信号决定；

S2.2：第一次膨胀运算；

利用形态学滤波器运算中的膨胀运算对数据进行进一步的处理，改善信号失真情况，获得信号{x₃(t)}，处理过程如式(2)所示：

{x₃(t)}＝ζ_S[x₂(t)]＝max[x₂(t-m)+S(m)] (2)

S2.3：第二次腐蚀运算；

利用形态学滤波器运算中的腐蚀运算对数据进行进一步处理，进一步去除粗大误差，获得信号{x₄(t)}，处理过程如式(3)所示：

{x₄(t)}＝ξ_S[x₃(t)]＝min[x₃(t+m)-S(m)] (3)

S2.4：第二次膨胀运算；

利用形态学滤波器运算中的膨胀运算对数据进行进一步的处理，进一步改善信号失真情况，获得信号{x₅(t)}，处理过程如式(4)所示：

{x₅(t)}＝ζ_S[x₄(t)]＝max[x₄(t-m)+S(m)] (4)。

步骤S3包括以下步骤：

S3.1：信号延拓；

选择步骤S2得到的信号的左端长度为t₁的待匹配信号段{η_l(t₁)}，以及右端长度为t₂的待匹配信号段{η_r(t₂)}，其中，{η_l(t₁)}＝{η_l(1),…,η_l(i),…η_l(t₁)}，{η_r(t₂)}＝{η_r(1),…,η_r(j),…η_r(t₂)}；步骤S2得到的信号中，除掉待匹配信号段{η_l(t₁)}和待匹配信号段{η_r(t₂)}之后的信号为剩余信号；

首先，将{η_l(t₁)}向右平行移动，在移动的过程中，使{η_l(t₁)}的极大值点依次与剩余信号中的每个极大值点重合，记相应的匹配信号段分别为θ₁,θ₂…θ_i…θ_n，并通过式(5)计算各个匹配波形与待匹配信号段{η_l(t₁)}的相似距离：

其中，

式(8)中的cov[·]是协方差，σ(·)是方差；

选择P_i最小时对应的匹配信号段作为第一最佳匹配信号段，并利用第一最佳匹配信号段对待匹配信号段{η_l(t₁)}进行向左延拓，以待匹配信号段{η_l(t₁)}的极大值点为起点，将第一最佳匹配信号段极大值点左侧的信号段复制到待匹配信号段{η_l(t₁)}极大值点的左边，使得第一最佳匹配信号段的极大值点与待匹配信号段{η_l(t₁)}的极大值点重合；

然后，将{η_r(t₂)}向左平行移动，在移动的过程中，使{η_r(t₂)}的极大值点依次与剩余信号中的每个极大值点重合，记相应的匹配信号段分别为α₁,α₂…α_j…α_n，并通过式(9)计算各个匹配波形与待匹配信号段{η_r(t₂)}的相似距离：

其中，

式(12)中的cov[·]是协方差，σ(·)是方差；

选择P_j最小时对应的匹配信号段作为第二最佳匹配信号段，并利用第二最佳匹配信号段对待匹配信号段{η_r(t₂)}进行向右延拓，以待匹配信号段{η_r(t₂)}的极大值点为起点，将第二最佳匹配信号段极大值点右侧的信号段复制到待匹配信号段{η_r(t₂)}极大值点的右边，使得第二最佳匹配信号段的极大值点与待匹配信号段{η_r(t₂)}的极大值点重合；

步骤S2得到的信号经过向左延拓和向右延拓处理之后，得到的信号记为{x₆(t₃)}；

S3.2：模态分解；

如图2所示，包括以下步骤：

S3.21：令k＝1，p＝1；

S3.22：寻找信号{x₆(t₃)}中所有极大值点，通过三次样条插值拟合得到极大值包络线m_k+(t₃)，寻找信号{x₆(t₃)}中所有极小值点，通过三次样条插值拟合得到极小值包络线m_k-(t₃)，取均值包络线为

S3.23：将信号{x₆(t₃)}减去m_k(t₃)得到新信号判断新信号是否满足“极值点个数和过零点个数相差不大于1并且在任意点处上下包络均值为0”这个准则：如果满足，则令为第p个本征信号函数，记为IMF_p(t₃)，然后进行步骤S3.24；如果不满足，则令k＝k+1，p＝1，然后返回步骤S3.22；

S3.24：将信号{x₆(t₃)}减去IMF_p(t₃)，得到一个去掉高频信号的新信号r_p(t₃)；

S3.25：判断新信号r_p(t₃)是否还能再进行分解：如果能，则令p＝p+1，k＝1，然后返回步骤S3.22；如果不能，则记录分解的总次数为n，最终得到的不能分解的新信号为r(t₃)，信号{x₆(t₃)}如式(13)所示：

S3.26：寻找IMF₁(t₃),…,IMF_n(t₃)中的突变点：按照式(14)计算IMF₁(t₃),…,IMF_n(t₃)中相邻两个本征信号函数之间的均方差CMSE，如果q＝u时，CMSE(IMF_u(t₃),IMF_u+1(t₃))发生了突变，则将x(t₃)作为重构信号，如式(15)所示；