CN107045149A - 一种基于双奇异值分解的全波核磁共振信号噪声滤除方法 - Google Patents

Abstract

本发明涉及核磁共振测深信号噪声滤除领域，具体是基于双奇异值分解的全波核磁共振信号噪声的滤除方法，包括：利用地面核磁共振地下水探测仪器采集一组观测MRS含噪数据；对数据进行以拉莫尔频率为中心的带阻滤波，实现初步信噪分离；进行第一次奇异值分解，重构出噪声数据；用观测MRS含噪数据减去重构出的噪声数据，得到以信号能量为主、噪声能量相对较小的数据；进行第二次奇异值分解，重构出想要提取的MRS信号；利用非线性拟合方法进行MRS信号的特征参数提取。本发明针对单通道采集的全波MRS数据，可以同时去除工频谐波干扰和随机噪声的影响，通过非线性拟合实现有效MRS信号特征参数的提取。

Description

一种基于双奇异值分解的全波核磁共振信号噪声滤除方法

技术领域

本发明涉及核磁共振测深(Magnetic ResonanceSounding，MRS)信号噪声滤除领域，具体是利用双奇异值分解的方法对全波核磁共振信号噪声的滤除方法。

背景技术

核磁共振测深(magnetic resonance sounding，MRS)技术作为目前唯一一种直接有效探测地下水的地球物理勘探方法，与常规只能通过勘察含水构造和层位来间接找水的地球物理勘探方法相比，能够对地下含水层深度和厚度、含水量的大小、地下介质孔隙度等信息给出定量解释，因而被广泛应用在水资源勘查、滑坡地质灾害检测、矿井突水和隧道涌水等领域。其基本原理是通过探测地下水中氢质子共振跃迁产生的核磁共振信号来实现地下水探测的。全波核磁共振信号的一般表达式为:

式(1)中包含MRS信号的四个关键特征参数，即与地层含水量相关的信号初始振幅E₀、与含水层孔隙大小相关的弛豫时间(范围为30ms-1000ms)、不同区域具有不同数值的拉莫尔角频率ω₀以及与含水层导电性相关的初始相位

MRS方法虽然比其他传统方法具有明显优势，但也存在不足，实际探测到的MRS信号是微弱的纳伏级信号，环境中的尖峰噪声、随机噪声和工频谐波噪声等会被高灵敏度的仪器系统一并采集，进而影响获取信号的质量。尖峰噪声通常具有幅值大于信号，持续时间短等特点，一般采用统计分析、非线性能量算子等方法将其去除，而由电力线干扰引起的工频谐波噪声以及环境中的随机噪声对MRS信号的影响较大，如何实现强干扰下有效信号的提取一直是国内外学者研究的热点问题。目前，针对全波核磁共振信号中噪声去除的理论和方法有很多。2008年，美国推出的GMR仪器和法国推出的NUMIS^poly仪器均采用了自适应噪声对消(Adaptive Noise Cancellation，ANC)算法进行噪声抑制。2012年，田宝凤等人在论文《基于参考线圈和变步长自适应的磁共振信号噪声压制方法》(《地球物理学报》，2012年第55卷7期:2462-2472页.)中直接利用观测数据不断递归更新处理参数，实现不同信噪比和信号强度下MRS信号噪声的去除，但该方法需要远端线圈同时采集噪声数据，对于单通道仪器无法实现，而且只能对消相关噪声，无法去除不相关的工频谐波噪声和随机噪声的影响。2012年，Dalgaard等人在论文《Adaptive noise cancelling of multichannelmagnetic resonance sounding signals》(《Geophysical Journal International》，2012年第191卷1期:88-100页.)中对比了维纳滤波器和ANC算法的滤波性能，结果表明二者滤波效果相当。2013年，Larsen等人在论文《Noise cancelling of MRS signals combiningmodel-based removal of powerline harmonics and multichannel Wiener filtering》(《Geophysical Journal International》，2013年第196卷2期:828-836页.)中提出基于工频谐波建模和多通道维纳滤波相结合的方法进行噪声抑制，同样需要多通道探测方式，且滤波效果取决于噪声源的分布情况。2014年，Costabel和Müller-Petke在论文《Comparisonand optimal parameter setting of reference-based harmonic noise cancellationin time and frequency domain for surface-NMR》(《Near Surface Geophysics》，2014年第12卷2期:199-210页)中通过对比分析时域和频域两种噪声对消算法，证明了频域方法性能更佳。2015年，訾彦勇在论文《基于EMD的磁共振探测信号噪声抑制方法研究》(吉林大学硕士学位论文,2015.)中提出一种新型的信号时频处理方法，可以在没有输入信号任何先验知识的情况下，自适应地将信号分解成若干个固有模态函数(IMF)，实现信号趋势的有效提取，但存在端点效应和模态混叠问题。2015年，田宝凤等在论文《基于独立成分分析的全波核磁共振信号噪声滤除方法研究》(《物理学报》，2015年64卷22期:446-457页)，提出了基于独立成分分析(Independent Component Analysis，ICA)方法进行多源混合噪声下MRS信号中工频谐波干扰或单频干扰的去除。2016年，Ghanati等人在论文《Filtering andparameter estimation of surface-NMR data using singular spectrum analysis》(《Journal of Applied Geophysics》,2016年第130期:118-130页.)提出了基于奇异谱分析(Singular Spectrum Analysis)方法进行地面磁共振数据的滤波和参数估计，滤波性能取决于嵌入窗口长度和奇异值个数两个参数的选择。上述方法均取得了一定程度上MRS信号信噪比的提升，但是由于噪声源众多、空间分布不均匀、复杂多变等因素，使得磁共振测深领域的消噪技术研究依然具有挑战性。

奇异值分解(Singular Value Decomposition，SVD)是一种基于矩阵分解的非线性滤波方法，广泛应用于信号处理、图像压缩以及统计学等领域。该方法通过将包含信号信息的矩阵分解到一系列奇异值以及奇异值矢量对应的时频子空间中，依据信号和噪声对奇异值的贡献不同，从而实现噪声滤除。

专利CN105607125A公开了“一种基于块匹配算法和奇异值分解的地震资料噪声压制方法”，专利CN105319593A公开了“一种基于曲波变换和奇异值分解的联合去噪方法”，均属于地震信号处理领域；专利CN104200441A公开了“一种基于高阶奇异值分解的磁共振图像去噪方法”，属于生物医学图像处理领域；专利CN103810394A公开了“一种旋转设备故障信号奇异值分解降噪的设计方法”，属于机械故障信号处理领域。可见奇异值分解已被成功应用到了信号处理的各个领域，但尚未见其应用于MRS信号的噪声滤除中。

发明内容

本发明所要解决的技术问题在于提供一种基于双奇异值分解的全波核磁共振信号噪声滤除方法，解决磁共振地下水探测工作中由于强工频谐波干扰和随机噪声造成的MRS信号有效提取的问题。

本发明是这样实现的，

一种基于双奇异值分解的全波核磁共振信号噪声滤除方法，包括以下步骤：

步骤(1)：利用地面核磁共振地下水探测仪器采集一组观测MRS含噪数据；

步骤(2)：对步骤(1)的数据进行以拉莫尔频率为中心的带阻滤波，实现初步信噪分离；

步骤(3)：将步骤(2)的结果进行第一次奇异值分解，重构出噪声数据；

步骤(4)：用步骤(1)的观测MRS含噪数据减去步骤(3)的重构出的噪声数据，得到以信号能量为主、噪声能量相对较小的数据；

步骤(5)：对步骤(4)的结果进行第二次奇异值分解，重构出想要提取的MRS信号；

步骤(6)：对步骤(5)后的结果利用非线性拟合方法进行MRS信号的特征参数提取。

进一步地，步骤(3)中的第一次奇异值分解方法的具体步骤为：

步骤3a：确定分解阶数L、延迟步长τ，对步骤(2)得到的数据构建吸引子轨迹矩阵A；

步骤3b：对吸引子轨迹矩阵A进行第一次奇异值分解；

步骤3c：分解后得到K个奇异值(λ₁,λ₂,....λ_P...λ_K)且按照降序排列，前P个较大的奇异值主要代表噪声数据；

步骤3d：选取前P个较大的奇异值，其余奇异值置零；

步骤3e：利用奇异值分解的逆过程，重构出噪声数据。

进一步地，延迟步长τ＝3～5，分解阶数L满足分解阶数L大于非零奇异值的个数,使得吸引子轨迹矩阵A分解充分。

进一步地，步骤(5)中的第二次奇异值分解方法的具体步骤为：

步骤5a：确定分解阶数L1、延迟步长τ1，对步骤(4)得到的数据构建吸引子轨迹矩阵A1；

步骤5b：对吸引子轨迹矩阵A1进行第二次奇异值分解；

步骤5c：第二次分解后，前两个较大的奇异值主要代表想要提取的MRS信号，其余奇异值代表噪声；

步骤5d：选取前两个较大的奇异值，其余奇异值置零；

步骤5e：利用奇异值分解的逆过程，重构得到想要提取的MRS信号。

进一步地，延迟步长τ1＝3～5，分解阶数L1满足分解阶数L1大于非零奇异值的个数，使得吸引子轨迹矩阵A分解充分。

本发明与现有技术相比，有益效果在于：本发明与现有技术相比，其有益效果在于：提出了基于双奇异值分解的全波核磁共振信号噪声滤除方法，针对单通道采集的全波MRS数据，可以同时去除工频谐波干扰和随机噪声的影响，通过非线性拟合实现有效MRS信号特征参数的提取。本发明方法解决了磁共振地下水探测工作中由于强工频谐波干扰和随机噪声造成的MRS信号有效提取的难题，为仪器在城镇、村庄等强干扰环境下开展有效探测提供了技术保障，并且，采用本发明进行MRS信号消噪后获取的初始振幅和弛豫时间等关键特征参数的拟合误差较小，为后续精确反演提供了高质量的数据信息。同时，本发明摆脱了经典消噪方法需要借助多通道探测来提高信噪比的限制，节省大量人力物力资源，提高了工作的时效性。

附图说明

图1：本发明实施例提供的基于双奇异值分解的全波核磁共振信号噪声滤除方法流程框图；

图2：本发明实施例提供的含噪MRS信号带阻滤波后的时频图，图2(a)黑色是理想的MRS信号时域图，参数为e₀＝150nV，f＝2325Hz和灰色是在理想MRS信号中加入-2dB的高斯随机噪声和频率2200Hz～2450Hz之间的6个幅值大小可调的工频谐波的信号时域图，此时数据的信噪比为-10.5271dB，MRS信号完全淹没在噪声中；图2(b)黑色线是理想MRS信号频谱，灰色是含噪MRS信号频谱；图2(c)是含噪MRS信号滤除以拉莫尔频率为中心的有用信号后时域图；图2(d)为2(c)对应的频谱；

图3：本发明实施例提供的双SVD奇异值分布图，图3(a)是第一次SVD的奇异值分布图，前P个较大的奇异值代表噪声数据；图3(b)是第二次SVD的奇异值分布图，前两个较大的奇异值代表想要提取的MRS信号；

图4：本发明实施例提供的双SVD消噪原理及效果图，图4(a)是第一次SVD后重构的噪声数据的时域图；图4(b)为4(a)对应的频谱；图4(c)是原观测MRS含噪数据减去第一次SVD后重构的噪声数据的时域图；图4(d)为4(c)对应的频谱；图4(e)是第二次SVD后重构出的想要提取MRS信号时域图；图4(f)为4(e)对应的频谱；

图5：本发明实施例提供的仿真数据和数据处理结果图，图5(a)深灰色是含噪信号时域图，浅灰色是理想MRS信号时域图，黑色是消噪处理后信号时域图；图5(b)为5(a)对应的频谱；

图6：本发明实施例提供的实测MRS含噪观测数据及带阻滤波后的时频图，图6(a)实测MRS含噪观测数据时域图，图6(b)为6(a)对应频谱；图6(c)是实测MRS含噪数据滤除以拉莫尔频率为中心的有用信号后时域图；图6(d)为6(c)对应的频谱；

图7：本发明实施例提供的双SVD奇异值分布图，图7(a)是第一次SVD的奇异值分布图，前m个较大的奇异值代表噪声数据；图7(b)是第二次SVD的奇异值分布图，前两个较大的奇异值代表想要提取的MRS信号；

图8：本发明实施例提供的双SVD消噪原理及效果图，图8(a)是第一次SVD后重构的噪声数据的时域图；图8(b)为8(a)对应的频谱；图8(c)是原实测观测MRS含噪数据减去第一次SVD后重构的噪声数据的时域图；图8(d)为8(c)对应的频谱；图8(e)是第二次SVD后重构的想要提取的MRS信号时域图；图8(f)为8(e)对应的频谱；

图9：本发明实施例提供的实测数据处理前后时频图，图9(a)灰色是实测数据时域，黑色是双SVD处理后信号时域图；图9(b)为9(a)对应的频谱；

图10本发明实施例提供的选择不同分解阶数L时的奇异值分布情况图；图10(a)中L为100阶时奇异值分布情况图；图10(b)中L为300阶时奇异值分布情况图、图10(c)中L为500阶时奇异值分布情况图、图10(e)中L为800阶时奇异值分布情况图、图10(f)中L为1000阶时奇异值分布情况图，图10(g)中L为1300阶时奇异值分布情况图；

图11本发明实施例提供的不同信噪比下当分解阶数L已确定时，关于延迟步长τ的算法选择图。11(a)为信噪比为-5dB时(下上有一定的浮动，这里为了方便计为整数，以下同理)，延迟步长τ的算法选择图、图11(b)为信噪比为-10dB时，延迟步长τ的算法选择图、图11(c)为信噪比为-15dB时，延迟步长τ的算法选择图、图11(d)为信噪比为-20dB时，延迟步长τ的算法选择图、图11(e)为信噪比为-25dB时，延迟步长τ的算法选择图、图11(f)为信噪比为-30dB时，延迟步长τ的算法选择图。

具体实施方式

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

如图1所示，一种基于双奇异值分解的全波核磁共振信号噪声滤除方法，包括以下步骤：

步骤(1)：利用地面核磁共振地下水探测仪器采集一组观测MRS含噪数据；

步骤(2)：对其进行以拉莫尔频率(f＝2325Hz)为中心的带阻滤波，实现初步信噪分离；

步骤(3)：将步骤(2)的结果进行第一次奇异值分解，重构出噪声数据；

步骤(4)：用步骤(1)的观测MRS含噪数据减去步骤(3)的结果，得到以信号能量为主、噪声能量相对较小的数据；

步骤(5)：对步骤(4)的结果进行第二次奇异值分解，重构出想要提取的MRS信号；

步骤(6)：对步骤(5)后的结果采用非线性拟合方法进行MRS信号参数提取。

如图2所示，步骤(2)中的带阻滤波具体步骤为：

步骤2a：设定带阻滤波器的阻带频率范围2320Hz～2330Hz；

步骤2b：将观测MRS含噪数据(其时频图如图2(a)和图2(b)灰色所示)经过带阻滤波器，滤除以拉莫尔频率为中心的有用信号，得到以强工频谐波和随机噪声为主的数据,其时频图如图2(c)和图2(d)所示。

如图3所示，步骤(3)中的第一次奇异值分解方法的具体步骤为：

步骤3a：确定合适的分解阶数L、延迟步长τ两个关键参数，对步骤(2)得到的数据构建吸引子轨迹矩阵A；

步骤3b：对吸引子轨迹矩阵A进行第一次奇异值分解；

步骤3c：分解后得到K个奇异值(λ₁,λ₂,....λ_P...λ_K)且按照降序排列,此时前P个较大的奇异值主要代表噪声数据，如图3(a)所示；

步骤3d：选取前P个较大的奇异值，其余奇异值置零；

步骤3e：利用奇异值分解的逆过程，重构得到噪声数据，其时频图如图4(a)、4(b)所示；

如图4所示，步骤(4)中用步骤(1)的观测MRS含噪数据减去第一次SVD后重构的噪声数据的时频图如图4(c)、4(d)所示。

步骤(5)中第二次奇异值分解方法的具体步骤为：

步骤5a：确定合适的分解阶数L1、延迟步长τ1两个关键参数，对步骤(4)得到的数据构建吸引子轨迹矩阵A1；

步骤5b：对吸引子轨迹矩阵A1进行第二次奇异值分解；

步骤5c：第二次分解后，前两个较大的奇异值主要代表想要提取的MRS信号，其余奇异值代表噪声，如图3(b)所示；

步骤5d：选取前两个较大的奇异值，其余奇异值置零；

步骤5e：利用奇异值分解的逆过程，重构得到想要提取的MRS信号，如图4(e)、4(f)所示消噪后的时频图；

步骤6：对步骤(5)后的结果利用非线性拟合方法进行MRS信号的特征参数提取。

以第一次奇异值分解的过程为例：构建吸引子轨迹矩阵的过程为：

设仪器系统探测到观测MRS含噪数据序列N＝[x₁ x₂ .... x_n]，构造的吸引子轨迹矩阵下式(2)所示：

其中，A为m×n的矩阵，参数τ称为延迟步长、参数L称为分解阶数(当x_τ+n≥x_n时，则x_τ+n+1＝0)。这两个参数的选择是奇异值分解法用于磁共振测深信号噪声抑制能否取得成功的关键。

步骤3b中对吸引子轨迹矩阵A进行第一次奇异值分解的过程为：

观测含噪MRS数据N是由MRS信号和噪声共同组成，则吸引子轨迹矩阵A也是由MRS信号和噪声共同组成的轨迹矩阵，对矩阵A进行奇异值分解如公式(3)所示：

A＝UΣV^H (3)

Σ表示非零奇异值对角阵：

Σ＝diag(λ₁,λ₂,....λ_K)是对角矩阵，其对角元素λ₁,λ₂,....λ_K为吸引子轨迹矩阵A分解后的奇异值，并按降序排列，即：λ₁≥λ₂≥....≥λ_K。K是吸引子轨迹矩阵A的秩。由于Σ是一对角阵，因此SVD可以将一个秩为K的m×n阶矩阵A表示为K个秩为1的m×n阶子矩阵的和。其中，每个子矩阵分别由2个特征向量u_i和v_i与权值相乘得到，可表示为：

式(5)中，u_i和v_i分别为矩阵U和矩阵V的第i个列向量，且u_i和v_i两两正交，它们分别构成了信号时、频信息；λ_i是矩阵Σ的第i个奇异值；X_i是包含u_i和v_i的子矩阵，即时—频子空间。

提取有效特征重构信号的过程为：

由公式(5)知矩阵A的奇异值λ_i可以反映信号和噪声的能量集中情况，保留矩阵A奇异值分解后与信号相对应的前P个奇异值项，而其它和噪声对应的奇异值置零，即可选取相应的子空间进行逆运算，从而重构出特定成分有效信号，实现噪声的去除。

X′＝UΣ_pV^H (6)

在双SVD算法仿真的过程中，延迟步长τ和分解阶数L的选择至关重要，直接决定算法的消噪性能。对于分解阶数L的选择，选取的原则是使分解阶数L大于非零奇异值的个数，即矩阵充分分解。本实施例中参数选择在不同信噪比下进行了大量仿真实验。图10给出了SNR＝-10.5271dB时分解阶数L选择图。选择不同分解阶数L时的奇异值分布情况如图10所示。当分解阶数L为100、300、500、800、1000时，非零奇异值的个数与分解阶数相等并且最大奇异值代表的能量值在不断增大，当分解阶数N为1300时，此时非零奇异值的个数为1250，矩阵完全分解，符合选取原则。

图11为不同信噪比下当分解阶数L已确定时，关于延迟步长τ的算法选择图。从图11中可看出当SNR＝-5dB时，τ＝2或3平均弛豫时间T₂*的拟合误差在±5％以内，并且从τ＝2开始初始振幅E₀的相对误差均在±0％～±1％内波动。当SNR＝-10dB时，由图11可看出τ＝3～5时初始振幅E₀和平均弛豫时间T₂*的拟合误差均在±5％以内。当SNR＝-15dB--30dB和τ＝3～8时，初始振幅E₀和平均弛豫时间T₂*的拟合误差均在±5％以内，而若τ过大，构建吸引子矩阵时则需要更多的观测数据。因此经过大量仿真实验结合实际不同信噪比情况选择最佳延迟步长为τ＝3～5之间。

实施例1

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

基于双奇异值分解的全波核磁共振信号噪声滤除方法的仿真算法，参照附图1，包括以下步骤：

步骤(1)：利用式(1)构造初始振幅e₀＝150nV，弛豫时间拉莫尔频率f＝2325Hz和相位的理想MRS信号，在信号中加入-2dB的高斯随机噪声和频率2200Hz～2450Hz之间的6个幅值大小可调的工频谐波，此时MRS信号完全淹没在噪声中，构成观测数据的信噪比为-10.5271dB，如图2(a)所示，图2(b)为其对应的频谱。

步骤(2)：对观测数据进行以拉莫尔频率为中心的带阻滤波，实现初步信噪分离，如图2(c)所示为含噪MRS信号滤除以拉莫尔频率为中心的有用信号后时域图，图2(d)为其对应的频谱；从图2(d)中的频谱图可以看出，已滤除以拉莫尔频率为中心的有用信号，剩余主要是工频谐波和随机噪声数据。

步骤(3)：将步骤(2)的结果进行第一次奇异值分解，选择如图3(a)所示的前50个较大奇异值，重构出噪声数据，如图4(a)是第一次SVD后重构的噪声数据的时域图，图4(b)为其对应的频谱；由图4(a)、(b)的时频图可看出，第一次SVD重构出了部分高斯随机噪声和频率2200Hz～2450Hz之间的6个幅值大小可调的工频谐波的主要噪声数据。

步骤(4)：用步骤(1)的观测MRS含噪数据减去步骤(3)的结果，得到以信号能量为主、噪声能量相对小的数据，图4(c)是用步骤(1)的观测MRS含噪数据减去第一次SVD后重构的噪声数据的时域图，图4(d)为其对应的频谱；其中从图4(d)频谱图中可以看出此时得到的数据为以信号的能量为主，同时还有少量随机噪声和残余的工频谐波。

步骤(5)：对步骤(4)的结果进行第二次奇异值分解、选择如图3(b)所示的前2个较大奇异值，重构出想要提取的MRS信号，图4(e)是第二次SVD后重构出的想要提取的MRS信号时域图，图4(f)为其对应的频谱；从图4(e)时域图可以看出在信号中加入-2dB的高斯随机噪声和频率2200Hz～2450Hz之间的6个幅值大小可调的工频谐波均已被滤除，由时频图可看出此时MRS信号完全提取出来。

步骤(6)：为了验证本方法的实用性，将消噪后的MRS信号进行信噪比(SNR)估计。经计算SNR＝20.3334dB，较消噪前信噪比提高30.8605dB；接着利用非线性拟合方法求解消噪后MRS信号的特征参数E₀＝149.93nV、相对误差分别为-0.05％、-2.87％，均控制在±5％以内，满足应用要求。

实施例2

本实例是在长春市吉林大学地质宫开展，该测试中MRS信号采集仪器采用的是吉林大学自主研发的“核磁共振阵列式接收仪(JLMRS-ARRAY)”，信号源设置为频率等于2326Hz，初始振幅0.5V，弛豫时间180ms的全波MRS信号，通过发射线圈发射该模拟信号。发射线圈与接收线圈间保持一定的距离，确保接收线圈能耦合MRS信号和环境噪声。如附图1所示，基于双奇异值分解的全波核磁共振信号噪声滤除方法，包括以下步骤：

步骤(1)：利用地面核磁共振地下水探测仪器采集一组观测MRS含噪数据，SNR＝-11.5324dB,如图6(a)为该观测数据时域图，图6(b)为其对应频谱，从图中可看到其含有大量强工频谐波和随机噪声；

步骤(2)：对观测数据进行以拉莫尔频率为中心的带阻滤波，实现初步信噪分离如图6(c)是滤除主要信号之后的时域图，图6(d)是其对应的频谱；从图6(d)中的频谱图可以看出，已滤除以拉莫尔频率为中心的有用信号，剩余主要是噪声数据。

步骤(3)：将步骤(2)的结果进行第一次奇异值分解、选择如图7(a)所示的前30个较大奇异值，重构出噪声数据，如图8(a)是第一次SVD后重构的噪声数据的时域图，图8(b)为其对应的频谱；由图8(a)、(b)的时频图可看出，第一次SVD重构出了部分随机噪声、工频谐波等主要噪声数据。

步骤(4)：用步骤(1)的观测MRS含噪数据减去步骤(3)的结果，得到以信号能量为主、噪声能量相对小的数据，图8(c)是用步骤(1)的观测MRS含噪数据减去第一次SVD后重构的噪声数据的时域图，图8(d)为其对应的频谱；其中从图8(d)频谱图中可看出此时得到的数据是以想要提取的信号的能量为主，同时还有少量随机噪声和残余的工频谐波等噪声数据。

步骤(5)：对步骤(4)的结果进行第二次SVD、选择如图7(b)所示的前2个较大奇异值，重构出想要提取的MRS信号，图8(e)是第二次SVD后重构MRS信号时域图，图8(f)为其对应的频谱；从图8(e)时域图可以看出在信号中的随机噪声和工频谐波等噪声均已被消除，由时频图可看出此时MRS信号完全提取出来。

步骤(6)：为了验证本方法的实用性，将消噪后的MRS信号进行信噪比(SNR)估计。经计算SNR＝5.6471dB，较消噪前信噪比提高17.1795dB；接着利用非线性拟合方法求解消噪后MRS信号的特征参数相对误差分别为0.897％，满足应用要求。

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

Claims (5)

1.一种基于双奇异值分解的全波核磁共振信号噪声滤除方法，其特征在于，包括以下步骤：

步骤(1)：利用地面核磁共振地下水探测仪器采集一组观测MRS含噪数据；

步骤(2)：对步骤(1)的数据进行以拉莫尔频率为中心的带阻滤波，实现初步信噪分离；

步骤(3)：将步骤(2)的结果进行第一次奇异值分解，重构出噪声数据；

步骤(4)：用步骤(1)的观测MRS含噪数据减去步骤(3)的重构出的噪声数据，得到以信号能量为主、噪声能量相对较小的数据；

步骤(5)：对步骤(4)的结果进行第二次奇异值分解，重构出想要提取的MRS信号；

步骤(6)：对步骤(5)后的结果利用非线性拟合方法进行MRS信号的特征参数提取。

2.按照权利要求1所述的方法，其特征在于，步骤(3)中的第一次奇异值分解方法的具体步骤为：

步骤3a：确定分解阶数L、延迟步长τ，对步骤(2)得到的数据构建吸引子轨迹矩阵A；

步骤3b：对吸引子轨迹矩阵A进行第一次奇异值分解；

步骤3c：分解后得到K个奇异值(λ₁,λ₂,....λ_Ρ...λ_K)且按照降序排列，前P个较大的奇异值主要代表噪声数据；

步骤3d：选取前P个较大的奇异值，其余奇异值置零；

步骤3e：利用奇异值分解的逆过程，重构出噪声数据。

3.按照权利要求2所述的方法，其特征在于，延迟步长τ＝3～5，分解阶数L满足分解阶数L大于非零奇异值的个数,使得吸引子轨迹矩阵A分解充分。

4.按照权利要求1所述的方法，其特征在于，步骤(5)中的第二次奇异值分解方法的具体步骤为：

步骤5a：确定分解阶数L1、延迟步长τ1，对步骤(4)得到的数据构建吸引子轨迹矩阵A1；

步骤5b：对吸引子轨迹矩阵A1进行第二次奇异值分解；

步骤5c：第二次分解后，前两个较大的奇异值主要代表想要提取的MRS信号，其余奇异值代表噪声；

步骤5d：选取前两个较大的奇异值，其余奇异值置零；

步骤5e：利用奇异值分解的逆过程，重构得到想要提取的MRS信号。

5.按照权利要求4所述的方法，其特征在于，延迟步长τ1＝3～5，分解阶数L1满足分解阶数L1大于非零奇异值的个数，使得吸引子轨迹矩阵A分解充分。