CN109828318B - 一种基于变分模态分解的磁共振测深信号噪声滤除方法 - Google Patents

Abstract

本发明涉及一种磁共振测深(Magnetic Resonance Sounding,MRS)信号噪声滤除领域，具体是一种基于变分模态分解的磁共振测深信号噪声滤除方法主要用于处理磁共振测深信号中的工频谐波噪声及随机白噪声。提出“三VMD”分解的方式以更好地实现含噪MRS信号中噪声的有效去除。对磁共振测深找水仪采集到的MRS信号经带通滤波后进行傅里叶变换，确定其中含有的工频谐波干扰及单频干扰的频率及个数，通过第一、二、三次VMD分解分别去除含噪MRS信号中高斯白噪声、大部分工频和距信号最近的工频，最终提取得到目标MRS信号。本发明不仅解决了应用传统模态分解方法后惯有的模态混叠问题，与传统MRS信号去噪方法相比，具有信噪比高、自适应性强等优点。

Description

一种基于变分模态分解的磁共振测深信号噪声滤除方法

技术领域

本发明涉及一种磁共振测深(Magnetic Resonance Sounding,MRS)信号噪声滤除领域，具体是一种基于变分模态分解的磁共振测深信号噪声滤除方法。

背景技术

地面核磁共振又称磁共振测深(Magnetic Resonance Sounding,MRS)，由于其直接对地下水中的氢质子敏感，因而是目前世界上唯一能够定量确定地下水存储状态的方法。该方法通过人工产生的交变电磁场激发地下水中的氢质子，受到激发的氢质子产生能级跃迁。当撤销人工场后，采集系统接收氢质子释能回迁时感应出的MRS信号，实现对地下水的探测。全波核磁共振信号的一般表达式为：

该式包含MRS信号的四个关键参数：E₀、

f_L、

分别表示与地层含水量相关的信号初始振幅、与含水层孔隙度相关的弛豫时间(自由水的弛豫时间范围为30ms-1000ms)、不同区域具有不同数值的拉莫尔频率、与含水层导电性相关的初始相位。

目前常用的磁共振测深信号的降噪滤波方法有经验模态分解、集成经验模态分解和小波分析等，这些方法运行速度慢、抑噪性能差、提取精度低、实时性较弱。如在分解过程中会产生模态混叠现象，即分解得到的一个或多个模态分量中包含差异极大的特征时间尺度,信号和噪声混叠在一个或多个模态分量中,很难达到有效的降噪滤波效果。

MRS响应信号非常微弱，一般为纳伏级。该方法在实际应用中会受到复杂环境噪声的影响，导致MRS信号质量恶化，甚至信号完全被噪声所淹没。因而，磁共振探测的噪声抑制问题一直备受国内外专家学者的关注。对于工频谐波与随机噪声的滤除：Legchenko和Valla在论文《Removal of power-line harmonics from proton magnetic resonancemeasurements》(《Journal of Applied Geophysics》,2003年第53卷2期:103-120页.)提出采用区块对消法等用于消除MRS信号中工频谐波干扰。Dalgaard等人在论文《Adaptivenoise cancelling of multichannel magnetic resonance sounding signals》(《Geophysical Journal International》,2012年第191卷1期:88-100页.)提出采用自适应噪声对消算法进行多通道MRS探测信号的处理，实现了主通道与参考通道中相关噪声的抵消。田宝凤等人在论文《基于参考线圈和变步长自适应的磁共振信号噪声压制方法》(《地球物理学报》,2012年第55卷7期:2462-2472页.)提出了基于变步长LMS(Least MeanSquare,LMS)算法的自适应噪声对消算法进行MRS信号中工频谐波和部分随机噪声的滤除。Ghanati等人在论文《Joint application of a statistical optimization processempirical mode decomposition to MRS noise cancelation》(《Journal of AppliedGeophysics》,2014年第111卷5期:110-120页.)提出了基于EMD以及CEEMD方法进行有效MRS信号衰减趋势的提取。Larsen等人在论文《Noise cancelling of MRS signals combiningmodel-based removal of powerline harmonics and multichannel Wiener filtering》(《Geophysical Journal International》,2014年,第196卷2期:828-836页.)通过工频谐波建模结合多通道维纳滤波法实现了噪声的有效抑制。Müller-Petke和Costabel在论文《Comparison and optimal parameter setting of reference-based harmonic noisecancellation in time and frequency domain for surface-NMR》(《Near SurfaceGeophysics》,2014年第12卷2期:199-210页.)将基于远端参考的自适应噪声对消技术分别应用于时、频域方法上，结果证明频域方法具有更优的性能。Annette Hein等人在论文《Symmetry based frequency domain processing to remove harmonic noise fromsurface nuclear magnetic resonance measurements》(《Geophysical JournalInternational》,2017年第208卷4期:724–736页.)根据磁共振信号的频域对称性去除工频谐波噪声，将含噪数据的信噪比提升了10％。尽管上述算法在一定条件下取得了一定效果，但是由于噪声源的无规律性和不稳定性以及空间分布性，限制了算法的应用。

专利CN1967656公开了“一种基于VMD音乐格式的手机铃声合成方法”，属于声音信号处理领域；专利CN106019271A公开了“一种基于变分模态分解的多人穿墙时变呼吸信号检测方法”，属于雷达探测技术领域；专利CN107515424A公开了“一种基于VMD与小波包的微震信号降噪滤波方法”，属于地震信号处理技术领域；专利CN107832525A公开了“一种信息熵优化VMD提取轴承故障特征频率的方法及其应用”，属于故障检测技术领域。可见匹配追踪算法已被成功应用到了信号处理的各个领域，但尚未见其应用于MRS信号的噪声滤除中。

发明内容

针对MRS信号中所包含噪声的复杂性以及传统EMD算法容易产生模态混叠问题，本发明提供了一种基于变分模态分解的磁共振测深信号噪声滤除方法，该方法不仅能够解决模态混叠问题，还能实现工频谐波或某一单频干扰的快速、有效滤除。

本发明是这样实现的，

一种基于变分模态分解的磁共振测深信号噪声滤除方法，该方法如下步骤：

步骤1：对磁共振测深探水仪采集到的一组观测MRS信号X(t)，利用带通滤波的方式对其进行预处理，得到目标频带范围内含噪MRS信号x(t)；

步骤2：对预处理后的含噪MRS信号x(t)进行傅里叶变换，得到其频谱，确定信号中含有的工频谐波干扰或某一单频干扰f₁、f₂……f_N；

步骤3：根据VMD分解中根据步骤2中确定的工频谐波干扰数量N，对x(t)采用三VMD分解的方式，获得有效的信噪分离结果。

进一步地，步骤3包括：

步骤31：对x(t)进行第一次VMD算法，令K₁＝3，去除含噪MRS信号中的随机白噪声，并对分解后得到的模态，进行频谱分析，根据每个模态的频谱来判断各模态中的成分；

步骤32：根据步骤31选取同时含有MRS信号和工频的模态作为第二次VMD算法的输入，根据该模态中工频谐波干扰和单频干扰的数量N，取K₂＝N，并对分解后得到的模态，进行频谱分析，根据每个模态的频谱来判断各模态中的成分；

步骤33：根据步骤32选取只含有MRS信号和距离其最近工频的模态作为第三次VMD算法的输入，取K₃＝2，经分解后得到最终去噪后的目标MRS信号s(t)。

进一步地，对含噪MRS信号进行VMD分解，具体步骤为：

步骤3a：选取变分模态分量个数K值与惩罚因子α的值；

步骤3b：初始化

其中

指初始循环中第k个模态函数；

指的初始循环中

对应的根据功率谱中心估计的中心频率，

指的是初始循环中的

对应的拉格朗日算子；

步骤3c：令迭代次数n＝n+1，执行整个循环；

步骤3d：执行内层循环，根据式(1)、(2)更新u_k和ω_k；

其中，

为含噪MRS信号x(t)的傅里叶变换；

是第n+1次循环中第k个模态函数，

是第n+1次循环中u_k对应的根据功率谱中心估计的中心频率；

步骤3e：令k＝k+1，重复步骤3d直到k＝K，结束内层循环；

步骤3f：执行外层循环，根据式(3)更新λ；

其中,τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤3g：重复步骤步骤3c至步骤3f，直到满足迭代停止条件如式(4)所示，结束整个循环，输出结果，得到K个变分模态分量；

其中,ε为求解精度。

本发明与现有技术相比，有益效果在于：

本发明选择变分模态分解作为去噪方法，该方法整体框架是求解变分问题，利用交替方向乘子法不断更新各本征模态函数及其中心频率，并将各本征模态函数解调到相应的基频带，最终提取各个本征模态函数及相应的中心频率，实现信号的有效分解。鉴于变分模态分解的上述优势能够解决传统EMD存在的模态混叠问题，因此，本发明提出在磁共振测深信号处理过程中采用变分模态分解进行信号与噪声的分离。

本发明通过三VMD分解的方式，借助其算法自适应性的特点，自适应地选择频带的特性，进而提高了时频分辨率。此方法不仅解决了传统EMD方法存在的模态混叠问题，而且通过三次VMD分别选取不同的模态数和惩罚因子的数值，依次实现了对随机噪声以及工频谐波和其他单频干扰的去除问题。该方法的有效应用，对于提高磁共振探测方法的抗干扰能力具有重要的理论意义，对于拓展仪器系统的应用范围具有重要价值。

附图说明

图1基于变分模态分解的磁共振测深信号噪声滤除方法的流程框图；

图2变分模态分解算法的流程框图；

图3理想MRS信号及频谱，图3A为时间-幅值图，图3B为频率-幅值图；

图4含噪MRS信号及其频谱，图4A为时间-幅值图，图4B为频率-幅值图；

图5第一次VMD后MRS信号及其频谱，其中图5a和图5b分别为第一次VMD后MRS信号及其频谱，图5c和图5d分别为第一模态的时间-幅值图以及频谱，图5e和图5f分别为第二模态的时间-幅值图以及频谱，图5g和图5h分别为第三模态的时间-幅值图以及频谱；

图6第二次VMD后MRS信号及其频谱，其中图6a和图6b分别为第一次VMD后MRS信号及其频谱，图6c和图6d分别为第一模态的时间-幅值图以及频谱，图6e和图6f分别为第二模态的时间-幅值图以及频谱，图6g和图6h分别为第三模态的时间-幅值图以及频谱，图6i和图6j分别为第三模态的时间-幅值图以及频谱；

图7第三次VMD后MRS信号及其频谱，其中图7a和图7b分别为第一次VMD后MRS信号及其频谱，图7c和图7d分别为第一模态的时间-幅值图以及频谱，图7e和图7f分别为第二模态的时间-幅值图以及频谱；

图8VMD前后MRS信号及其频谱，图8A为MRS信号图谱，图8B为频谱；

图9实测MRS信号及其频谱，图9A为MRS信号图谱，图9B为频谱；

图10带通滤波后实测MRS信号及其频谱，图10A为MRS信号图谱，图10B为频谱；

图11VMD处理前后实测MRS信号及其频谱，图11A为MRS信号图谱，图11B为频谱。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，基于变分模态分解的磁共振测深信号噪声滤除方法，包括以下步骤：

步骤1：对磁共振测深(MRS)探水仪采集到的一组观测MRS信号X(t)，利用带通滤波的方式对其进行预处理，得到目标频带范围内含噪MRS信号x(t)；

步骤2：对预处理后的含噪MRS信号x(t)进行傅里叶变换，得到其频谱，确定信号中含有的工频谐波干扰或某一单频干扰f₁、f₂……f_N；

步骤3：对x(t)采用三VMD分解的方式，惩罚因子α依据经验值给出，获得有效的信噪分离结果。

其中步骤3中包括：

步骤31：对x(t)进行第一次VMD算法，令K₁＝3，以有效去除含噪MRS信号中的随机白噪声，并对分解后得到的模态，进行频谱分析；

步骤32：根据步骤31频谱分析结果选取同时含有MRS信号和工频的模态作为第二次VMD算法的输入，根据该模态中工频谐波干扰和单频干扰的数量N，此时取K₂＝N，并对分解后得到的模态，进行频谱分析；

步骤33：根据32选取只含有MRS信号和距离其最近工频的模态作为第三次VMD算法的输入，此时取K₃＝2，经分解后得到最终去噪后的目标MRS信号s(t)。

如图2所示，本发明提供的变分模态分解算法(VMD)进行含噪MRS信号x(t)的去噪，具体步骤为：

步骤3a：选取变分模态分量个数K值与惩罚因子α的值；

步骤3b：初始化

其中

指初始循环中第k个模态函数；

指的初始循环中

对应的根据功率谱中心估计的中心频率，

指的是初始循环中的

对应的拉格朗日算子；

步骤3c：令n＝n+1，执行整个循环；

步骤3d：执行内层循环，根据式(1)、(2)更新u_k和ω_k；

其中，

为含噪MRS信号x(t)的傅里叶变换。

步骤3e：令k＝k+1，重复步骤D，直到k＝K，结束内层循环；

步骤3f：执行外层循环,根据式(3)更新λ；

其中,τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤3g：重复步骤C至步骤F,直到满足迭代停止条件如式(4)所示,结束整个循环,输出结果,得到K个变分模态分量；

其中,ε为求解精度。

实施例1

本实施例是在MATLAB 7.0编程环境下开展的本发明方法的仿真实验。

基于变分模态分解的磁共振测深信号噪声滤除方法的仿真算法，参照图1，包括以下步骤：

步骤(1)：利用式

构造拉莫尔频率为2347Hz，幅值e₀为150nV，弛豫时间

为180ms的理想MRS信号，如图3所示其中，图3A为时间-幅值图，图3B为频率-幅值图。在该信号拉莫尔频率附近添加2250Hz、2300Hz、2350Hz和2400Hz的工频干扰和强度为-5dB的随机噪声经过某种线性组合形成信噪比为-11.32dB的观测MRS信号x(t)(为行向量)，如图4所示，图4A为时间-幅值图，图4B为频率-幅值图；

步骤(2)：按照权利要求2中所述步骤进行第一次VMD算法，令K₁＝3,α＝50000，以有效去除含噪MRS信号中的随机白噪声，并对分解后得到的模态，进行频谱分析，得到三个模态，模态数量根据K的取值而定，本步骤中K取值为3，是将第一模态的低频部分和第三模态的高频部分积聚了大量的白噪声，而第二模态中则保留了MRS信号和四个工频谐波，如图5所示，图5a和图5b分别为第一次VMD后MRS信号及其频谱，图5c和图5d分别为第一模态的时间-幅值图以及频谱，图5e和图5f分别为第二模态的时间-幅值图以及频谱，图5g和图5h分别为第三模态的时间-幅值图以及频谱。

步骤(3)：根据步骤(2)选取第一次VMD分解后得到的第二模态，根据步骤(2)中的分析，该模态里面有工频和MRS信号，需要进一步去除工频获得MRS信号，作为第二次VMD算法的输入，根据该模态中工频谐波干扰和单频干扰的数量N＝4，此时取K₂＝4,α＝50000，并对分解后得到的模态，进行频谱分析，如图6所示，图6a和图6b分别为第一次VMD后MRS信号及其频谱，图6c和图6d分别为第一模态的时间-幅值图以及频谱，图6e和图6f分别为第二模态的时间-幅值图以及频谱，图6g和图6h分别为第三模态的时间-幅值图以及频谱，图6i和图6j分别为第三模态的时间-幅值图以及频谱，可以看到在第二模态、第三模态、第四模态中分别含有2300Hz、2400Hz和2250Hz的工频，而MRS信号和距其最近的2350Hz工频仍然没有得到分离；

步骤(4)：根据步骤(3)选取第二次VMD分解后得到的第一模态含有MRS信号和距离其最近工频的模态作为第三次VMD算法的输入，此时取K₃＝2,α＝100000，经分解后得到最终去噪后的目标MRS信号如图7所示，图7a和图7b分别为第一次VMD后MRS信号及其频谱，图7c和图7d分别为第一模态的时间-幅值图以及频谱，图7e和图7f分别为第二模态的时间-幅值图以及频谱，可以看出经过第三次VMD分解后，MRS信号与距其最近的2350Hz工频得到有效分离。图8给出了经VMD处理前后MRS信号的时、频域对比图，图8A为MRS信号时域图,图8B为频谱。

为了验证本发明方法的实用性，将去噪后MRS信号s(t)进行了信噪比(SNR)估计。经计算，其SNR＝10.12dB，较分离前的SNR提高了21.44dB；此外，对s(t)进行了包络提取和数据拟合，获得关键参数初始振幅和弛豫时间分别为E₀＝154.31nV，T₂ ^*＝192.2ms，相对误差分别为2.87％、6.77％，均满足应用要求。

实施例2

本实施例是以长春市文化广场实地采集的MRS信号作为本发明方法的处理对象。

基于变分模态分解的磁共振测深信号噪声滤除方法，参照图1，包括以下步骤：

步骤1：对磁共振测深(MRS)探水仪采集到的一组观测MRS信号X(t)，如图9所示，图9A为MRS信号图谱，图9B为频谱，计算其信噪比为-8.35dB，利用带通滤波的方式对其进行预处理，得到目标频带范围内含噪MRS信号x(t)(为行向量)，计算其信噪比为SNR＝11.41dB；

步骤2：对带通滤波后的观测信号x(t)进行傅里叶变换，得到其频谱，如图10所示，图10A为MRS信号图谱，图10B为频谱，可以看到该信号在f₁＝2300Hz、f₂＝2350Hz、f₃＝2400Hz处均有较强的工频谐波干扰，此外，在f₄＝2265Hz、f₅＝2368Hz处受到了强烈的单频干扰；

步骤3：对x(t)进行第一次VMD算法，令K₁＝3,α＝50000，以有效去除含噪MRS信号中的随机白噪声，并对分解后得到的模态，进行频谱分析；

步骤4：根据步骤3选取第一次VMD分解后得到的第二模态作为第二次VMD算法的输入，根据该模态中工频谐波干扰和单频干扰的数量N＝5，此时取K₂＝5,α＝50000，并对分解后得到的模态，进行频谱分析；

步骤5：根据步骤3选取第二次VMD分解后得到的第一模态作为第三次VMD算法的输入，此时取K₃＝2,α＝500000，经分解后得到最终去噪后信噪比为SNR＝24.02dB的目标MRS信号s(t)，如图11所示，图11A为MRS信号图谱，图11B为频谱。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于变分模态分解的磁共振测深信号噪声滤除方法，其特征在于，该方法如下步骤：

步骤1：对磁共振测深探水仪采集到的一组观测MRS信号X(t)，利用带通滤波的方式对其进行预处理，得到目标频带范围内含噪MRS信号x(t)；

步骤2：对预处理后的含噪MRS信号x(t)进行傅里叶变换，得到其频谱，确定信号中含有的工频谐波干扰和某一单频干扰f₁、f₂……f_N；

步骤3：根据VMD分解中根据步骤2中确定的工频谐波干扰和单频干扰的数量N，对x(t)采用三次 VMD分解的方式，获得有效的信噪分离结果；

步骤3包括：

步骤31：对x(t)进行第一次VMD算法，令K₁＝3，去除含噪MRS信号中的随机白噪声，并对分解后得到的模态，进行频谱分析，根据每个模态的频谱来判断各模态中的成分；

步骤32：根据步骤31选取同时含有MRS信号和工频的模态作为第二次VMD算法的输入，根据该模态中工频谐波干扰和单频干扰的数量N，取K₂＝N，并对分解后得到的模态，进行频谱分析，根据每个模态的频谱来判断各模态中的成分；

步骤33：根据步骤32选取只含有MRS信号和距离其最近工频的模态作为第三次VMD算法的输入，取K₃＝2，经分解后得到最终去噪后的目标MRS信号s(t)。

2.按照权利要求1所述的方法，其特征在于，对含噪MRS信号进行VMD分解，具体步骤为：

步骤3a：选取变分模态分量个数K值与惩罚因子α的值；

步骤3b：初始化

其中

指初始循环中第k个模态函数；

指的初始循环中

对应的根据功率谱中心估计的中心频率，

指的是初始循环中的

对应的拉格朗日算子；

步骤3c：令迭代次数n＝n+1，执行整个循环；

步骤3d：执行内层循环，根据式(1)、(2)更新u_k和ω_k；

其中，

为含噪MRS信号x(t)的傅里叶变换；

是第n+1次循环中第k个模态函数，

是第n+1次循环中u_k对应的根据功率谱中心估计的中心频率；

步骤3e：令k＝k+1，重复步骤3d直到k＝K，结束内层循环；

步骤3f：执行外层循环，根据式(3)更新λ；

其中,τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤3g：重复步骤3c至步骤3f，直到满足迭代停止条件如式(4)所示，结束整个循环，输出结果，得到K个变分模态分量；

其中,ε为求解精度。

Images