CN109885903B - 一种基于模型的地面核磁共振信号尖峰噪声去除方法 - Google Patents

Abstract

本发明为一种基于模型的地面核磁共振信号噪声滤除方法，主要针对核磁共振信号中的尖峰噪声进行滤除。首先，利用核磁共振探水仪采集MRS信号，采用基于谐波建模方法去除工频谐波。其次，通过两个二阶级联带通滤波器建立尖峰噪声模型，同时，采用NEO算法定位出尖峰噪声的位置并截取相应的数据段。然后，基于最小二乘法进行参数提取，获得尖峰噪声。最后，将去除工频谐波的MRS数据与建模获得的尖峰噪声相减，得到去除尖峰噪声的MRS信号。本发明解决了MRS信号中尖峰噪声的去除问题，与传统MRS尖峰信号去噪方法相比，具有实用性强、明显提升工作效率，改善信噪比等优点。

Description

一种基于模型的地面核磁共振信号尖峰噪声去除方法

技术领域

本发明属于地面核磁共振(Surface Nuclear Magnetic Resonance，SNMR)，又称磁共振测深(Magnetic Resonance Sounding,MRS)信号的噪声滤除领域，具体是为一种基于模型的方案从地面核磁共振信号中去除尖峰噪声的方法。

背景技术

地面核磁共振(Surface Nuclear Magnetic Resonance，SNMR)，又称为磁共振测深(Magnetic Resonance Sounding,MRS)技术是近年发展起来的一种非侵入式的、直接、定量探测地下水的地球物理方法，其基本原理是通过探测地下水中氢质子共振跃迁产生的核磁共振信号来实现地下水的探测。核磁共振信号的表达式为：

式(1)中，E₀表示地面核磁共振信号的初始振幅，其与地层的含水量相关；

表示核磁共振探测信号的平均弛豫时间，与含水层的孔隙度相关，其范围一般在30ms-1000ms范围内；ω₀表示磁共振信号的拉莫尔角频率，与实验环境的地磁场强度成比例；

表示磁共振信号的初始相位，与地下介质的导电性相关。

相比于传统的探测地下水的方式，地面核磁共振探测的优点体现在分辨率高、效率高、信息量丰富，能够对地下含水层深度和厚度、含水量的大小、地下介质孔隙度等信息给出定量解释(Marian Hertrich.Imaging of groundwater with nuclear magneticresonance.Progress in Nuclear Magnetic Resonance Spectroscopy,2008,53(4):227-248)。

MRS方法虽然具有上述优势，但是其应用的主要瓶颈在于信号微弱，通常为纳伏级，在实际探测过程中，由于工作环境为野外天然地磁场，存在大量的自然和人为噪声，会对测量结果产生不可忽视的影响，造成测量数据信噪比低，进而影响反演结果的精度，最终导致含水层或者是含水量的判断误差。尖峰噪声作为人为噪声中的一类典型干扰，主要是由天然或人工放电引起的脉冲干扰，具有偶然发生、持续时间较短、幅值相对较大等特点，其存在将导致特征参数的提取存在较大误差，严重影响反演解释的准确性。

国内外的专家和学者在核磁共振信号尖峰噪声去除方面也开展了大量研究工作。Strehl和Rommel等提出根据不同噪声的特点不同，采用不同的策略进行噪声滤除的方法，其中采用小波变换手段对尖峰噪声进行处理，但小波基和分解级数的选择会对尖峰噪声去除效果造成影响(Strehl,S.,Rommel,I.,Hertrich,M.,Yaramanci,U..New strategiesfor filtering and fitting of MRS signals.Proceedings 3rd International MRSWorkshop,Madrid,Spain,2006,65-68.)。王中兴在吉林大学学报(工学版)发表的一篇文章中提出了奇异噪声的数学建模及其滤波算法(王中兴,荣亮亮,林君.地面核磁共振找水信号中的奇异干扰抑制.吉林大学学报(工学版),2009,39(5):1282-1287)，其采用差阈值替代叠加方法对奇异噪声进行削弱与抑制，提高了对信号参数的提取精度；万玲提出基于能量运算的磁共振信号尖峰噪声抑制方法(万玲,张扬,林君,蒋川东,林婷婷.基于能量运算的磁共振信号尖峰噪声抑制方法[J].地球物理学报,2016,59(06):2290-2301.)，通过推导能量域信号表达式并结合基于中位数的绝对偏差法确定阈值，实现了对尖峰噪声的有效识别与剔除；但这两种方法均需要多次重复采集测量数据才能得以实现。

专利号103150707A公开了一种消除磁共振图像中尖峰噪声的方法，能够修复信号强度较高的低频区域的尖峰噪声，专利号CN104835123A公开了一种基于先验模型的光片显微成像条纹噪声去除，均属于图像处理应用领域；专利号CN101523722A公开了一种尖峰噪声消除电路，采用硬件手段直接消除尖峰噪声，属于数字系统及IIC总线技术领域；专利号CN104218917A公开了一种用于消除APD雪崩信号输出端的尖峰噪声的自差分滤波装置，属于红外探测领域；专利号CN106452470A公开了一种消除信道尖峰脉冲噪声的无线通信系统及除噪方法，采用单匝环形天线接收尖峰脉冲噪声去抵消多匝环形天线信号中的尖峰噪声，属于无线通信技术领域。目前尚未见基于模型的方法去除尖峰噪声用于地面核磁共振地下水探测领域。

发明内容

本发明所要解决的技术问题在于提供一种基于模型的地面核磁共振信号尖峰噪声去除方法，解决由于尖峰噪声存在难以提取有效地面核磁共振信号的问题。

本发明是这样实现的，

一种基于模型的地面核磁共振信号尖峰噪声去除方法，该方法包括：

步骤(1)：利用核磁共振测深探水仪采集到一组观测MRS信号X₁(t)；

步骤(2)：对步骤(1)中获得的X₁(t)采用谐波建模方法去除工频谐波，得到X₂(t)；

步骤(3)：利用两个二阶级联带通滤波器与单位脉冲函数卷积建立尖峰噪声模型Y₁(t)，Y₁(t)＝Aδ(t-t₀)*g_rec(t)，A是尖峰噪声的幅值，t₀是尖峰信号的初始时间，g_rec(t)是滤波器传递函数的脉冲响应；

步骤(4)：对步骤(2)中得到的X₂(t)，通过NEO算法进一步强化尖峰，通过设置阈值找出尖峰所在位置x₁、x₂...；

步骤(5)：基于步骤(3)的模型Y₁(t)和步骤(4)定位的位置，截取尖峰数据段，采用最小二乘法提取特征参数A，t₀，从而获得尖峰噪声Y₂(t)；

步骤(6)：采用步骤(2)获得的X₂(t)减去步骤(5)获得的Y₂(t)，实现尖峰噪声的去除，获得输出MRS信号。

进一步地，所述的步骤(3)中建立尖峰噪声模型Y₁(t)＝Aδ(t-t₀)*g_rec(t)中求取g_rec(t)的具体步骤为：

构造四阶带通滤波器的传递函数G(s)＝G_a(s)G_b(s)，其是由两个不同的二阶带通滤波器级联组成的：

其中，s表示复数频率，ω_a为滤波器a的中心频率，σ_a为滤波器a的带宽；ω_b为滤波器b的中心频率，σ_b为滤波器b的带宽；

采用拉普拉斯求逆变换，计算其脉冲响应为：

进一步地，所述的步骤(5)中的最小二乘法提取特征参数的包括：

限定特征参数t₀，利用最小二乘公式，拟合出特征参数A，通过对比均方根误差：

选择具有最小均方根误差对应的特征参数A作为尖峰噪声模型的幅值参数，而此时的t₀为提取出的尖峰噪声初始时间参数，设采样率为f_s，将建立的模型按f_s采样，采样点数为N，记为向量y_m，将含噪数据段中的数据记为y_r，k＝1,2,...,N表示数据点数。

进一步地，步骤(5)所述的利用最小二乘法提取特征参数的具体步骤为：

51)数据预处理，将含噪数据段中的数据记为M，表示为[m₁...m_S]^T，将建立的模型以同样的采样率处理成S个数据点，记为向量H，表示为[h₁...h_S]^T；

52)提取特征参数A，设t_l是设定的扫描点，t₀为提取出的尖峰噪声起始时间参数；

根据式(7)求取在t₀＝t_l下的A值，记为A_l，存入向量A中；

A＝(M^TM)^-1M^TH (7)

53)利用式(8)计算t₀＝t_l下的均方根误差RMSE，存入向量RMSE中；

54)利用等间隔扫描形式求取特征参数t₀，设置扫描间隔为1/f_s，扫描次数取N/2，在每个t₀＝t_l下进行重复上述51)-53)步骤，获得向量A和RMSE，找出RMSE向量中最小值所对应的A值和t₀值。

进一步地，步骤(2)所述的谐波建模方法去除工频谐波的具体步骤为：

21)建立谐波噪声模型V_har，V_har是以工频f₀为基频的个谐波的累加结果：

其中，A_n和

分别是第n次谐波的幅度和初始相位；N是谐波个数；

22)通过自适应扫描方式确定基频f₀；

23)求解线性系数，将表达式(2)中的第n项线性化为式(4)：

其中，

由此估计A_n和

的问题转换成求解线性系数α_n和β_n，表达式(4)整理成线性方程组的形式：

其中，V_p(p＝1,2,…P)是t_p时刻的接收数据，表达式(5)整理成矩阵形式为：

Ax＝b (6)

其中x＝[α₁...α_n,β₁...β_n]^T,b＝[v₁,v₂...v_p]^T，表达式(6)是标准的线性方程组，采用最小二乘方法求解。

24)经过上述步骤求出V_har，用V_R减去V_har完成了消噪过程，其中V_R，表示测量数据X₁(t)。

进一步地，步骤(4)所述的利用NEO算法进一步强化尖峰噪声的具体步骤为：

41)采用NEO算法强化尖峰幅度：

其中x(k)表示测量数据X₁(t)离散化后的结果，经过NEO进行运算之后，再通过低通滤波器，就会使得尖峰信号突出；

42)设定阈值；

43)对尖峰噪声进行定位。

7.按照权利要求6所述的方法，其特征在于，设定阈值包括：

421)根据MRS持续时间T设定截取数据段时间t₁，取t₁为尖峰信号持续时间的1.5-2倍；

422)根据采样率和t₁设置两截取数据段之间的时间间隔t₂，保证t₂的整数n倍与t₁之和等于持续时间T，T＝t₁+n×t₂；

423)通过t₁、t₂两参数将经过低通滤波器后的NEO信号分为n+1段，将每段信号作为一行，将所有行组合成一数据矩阵，用ME表示；

424)利用median函数求取ME矩阵中每行的中位数，并将数据存入向量y中，作为拟合曲线的Y轴；

425)在横轴上取n+1个点，以T/(n+1)为时间间隔，存入向量x中，作为拟合曲线的X轴；

426)利用matlab自带的Curve Fitting Tool，选择向量x为其中的X data，向量y为Y data，类型为多项式Polynomial，进行曲线拟合，取3阶或4阶，得到多项式系数向量P；

427)利用系数向量P和阶数n建立该信号的趋势曲线，得到拟合函数，其表达式为

n是阶数；

428)拟合函数乘上一个阈值因子作为阈值函数。

本发明与现有技术相比，有益效果在于：本发明基于地面核磁共振探测信号中尖峰噪声具有信号幅值较大，持续时间较短，且由电网噪声产生的尖峰符合特定的传递函数特征，因而提出基于噪声建模的方式进行核磁共振探测信号中尖峰噪声的去除。该方法只需单次采集数据即可进行，避免了大量重复性测量带来的工作效率低下等问题，运算速度快，且明显改善探测信号的信噪比。

附图说明

图1为本发明实施例提供的基于模型的地面核磁共振信号尖峰噪声去除方法流程图

图2为本发明实施例提供的含噪声的SMRS信号X₁(t)；

图3为本发明实施例提供的去除工频谐波的信号X₂(t)；

图4为本发明实施例提供的NEO算法强化后提取MRS信号包络与低通滤波后结果；

图5为本发明实施例提供的阈值曲线；

图6为本发明实施例提供的参数辨识；

图7为本发明实施例提供的实现尖峰噪声的去除图，图7a为为去除工频谐波后的含噪MRS信号，图7b为为去除尖峰噪声后的MRS信号。

具体实施方式

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

参见图1所示，一种基于模型的地面核磁共振信号尖峰噪声去除方法，包括如下的步骤

步骤(1)：利用核磁共振测深(MRS)探水仪采集到一组X₁(t)观测MRS信号，如图2所示，核磁信号仪器以采样率f_s采集接收线圈中的MRS信号，该信号可以分解成四部分：

X₁(t)＝NMR(t)+h(t)+sp(t)+w(t) (10)

其中NMR(t)表示MRS信号，h(t)是电力线谐波噪声，sp(t)代表尖峰噪声，w(t)代表随机噪声。

步骤(2)：对步骤(1)中获得X₁(t)采用基于谐波建模方法去除工频谐波，得到X₂(t)，如图3，此时：

X₂(t)＝NMR(t)+sp(t)+w(t) (11)

谐波建模方法去除工频谐波的具体步骤为：

21)建立谐波噪声模型V_har，V_har是以工频f₀为基频的个谐波的累加结果：

其中，A_n和

分别是第n次谐波的幅度和初始相位；N是谐波个数；根据接收信号的带宽(1kHz～3kHz)，取N＝100。

22)通过自适应扫描方式确定基频f₀；首先在49.9～50.1Hz的范围内以0.03Hz为步长进行粗扫，利用表达式(3)获得2阶范数最小的3个频率点作为下一次扫描的范围。

||V_R-V_har(f₁)||₂→min (3)

其中，V_R为测量数据即X₁(t)，V_har(f_l)为第l个扫描值。在第2组扫描中，扫描范围为49.99～50.05Hz，扫描步长为0.0075Hz，再次获得范数最小的3个频点。当第3组扫描时，扫描范围缩小为50.0125～50.0275Hz，扫描步长为0.001875Hz。此时2阶范数的最小值对应的频点为50.018125Hz，即经过3组扫描即可达到精度要求。

23)求解线性系数，将表达式(2)中的第n项线性化为式(4)：

其中，

由此估计A_n和

的问题转换成求解线性系数α_n和β_n，表达式(4)整理成线性方程组的形式：

其中，V_p(p＝1,2,…P)是t_p时刻的接收数据，表达式(5)整理成矩阵形式为：

Ax＝b (6)

其中x＝[α₁...α_n,β₁...β_n]^T,b＝[v₁,v₂...v_p]^T，表达式(6)是标准的线性方程组，采用最小二乘方法求解。

24)经过上述步骤求出V_har，用V_R减去V_har即完成了消噪过程。

步骤(3)：由于电网的尖峰近似于地面NMR仪器的脉冲激励，故利用两个二阶级联带通滤波器与单位脉冲函数卷积来建立尖峰噪声模型Y₁(t)，则：

sp(t)＝Aδ(t-t₀)*g_rec(t) (12)

其中，A是尖峰噪声的幅值，t₀是尖峰信号的初始时间，g_rec(t)是滤波器传递函数。

其中，步骤3中，建立尖峰噪声模型Y₁(t)＝Aδ(t-t₀)*g_rec(t)中求取g_rec(t)的具体步骤为：

构造四阶带通滤波器的传递函数G(s)＝G_a(s)G_b(s)，其是由两个不同的二阶带通滤波器级联组成的：

s表示复数频率，ω_a为滤波器a的中心频率，σ_a为滤波器a的带宽；ω_b为滤波器b的中心频率，σ_b为滤波器b的带宽；

采用拉普拉斯求逆变换，计算其脉冲响应为：

步骤(4)：对步骤(2)中得到的X₂(t)，通过NEO算法进一步强化尖峰，如图4所示，41)首先利用NEO算法强化尖峰幅度：

其中x(k)表示测量数据X₁(t)离散化后的结果，t＝k/fs，k＝1,2,…,N-1，经过NEO进行运算之后，再通过适当的低通滤波器，就会使得尖峰信号突出；

42)然后设置阈值：包括：

411)根据MRS持续时间T设定截取数据段时间t₁，一般取t₁为尖峰信号持续时间的1.5-2倍；

412)根据采样率和t₁设置两截取数据段之间的时间间隔t₂，应保证t₂的整数n倍与t₁之和等于持续时间T，即T＝t₁+n×t₂。

413)通过t₁、t₂两参数将经过低通滤波器后的NEO信号分为n+1段，将每段信号作为一行，将所有行组合成一数据矩阵，用ME表示。

414)利用median函数求取ME矩阵中每行的中位数，并将数据存入向量y中，作为拟合曲线的Y轴。

415)在横轴上取n+1个点，以T/(n+1)为时间间隔，存入向量x中，作为拟合曲线的X轴。

416)利用matlab自带的Curve Fitting Tool，选择向量x为其中的X data，向量y为Y data，类型为多项式(Polynomial)，进行曲线拟合，一般取3阶或4阶即可，得到多项式系数向量P。

417)利用系数向量P和阶数n建立该信号的趋势曲线，其表达式为

其中n是阶数。

43)最后定位尖峰所在的位置:得到拟合函数envelop_line之后，乘上一个阈值因子(经验值1)作为阈值函数(如图5所示)，选取超过此阈值的部分作为尖峰噪声的定位点x₁、x₂...，利用这些定位点，在X₂(t)对应位置的前后截取27.5ms的数据，以便后续对尖峰噪声处理。

步骤(5)：基于步骤(3)的模型Y₁(t)和步骤(4)定位的位置截取的尖峰数据段，采用最小二乘法提取特征参数A，t₀，从而获得尖峰噪声Y₂(t)，如图6所示。对截取出来的数据进行最小二乘法处理。

提取特征参数A。可以求取在t₀＝t_l下的A值，记为A_l，存入向量A中。

A＝(M^TM)^-1M^TH (7)

计算t₀＝t_l下的均方根误差(RMSE)，存入向量RMSE_l中。

评价得到最小均方根误差对应的A,t₀为特征参数。具体为：

51)数据预处理，将含噪数据段中的数据记为M，表示为[m₁...m_S]^T，将建立的模型以同样的采样率处理成S个数据点，记为向量H，表示为[h₁...h_S]^T；

52)提取特征参数A，设t_l是设定的扫描点，t₀为提取出的尖峰噪声起始时间参数；

根据式(7)求取在t₀＝t_l下的A值，记为A_l，存入向量A中；

A＝(M^TM)^-1M^TH (7)

53)利用式(8)计算t₀＝t_l下的均方根误差RMSE，存入向量RMSE中；

54)利用等间隔扫描形式求取特征参数t₀，设置扫描间隔为1/f_s，扫描次数取N/2，在每个t₀＝t_l下进行重复上述51)-53)步骤，获得向量A和RMSE，找出RMSE向量中最小值所对应的A值和t₀值。

步骤(6)：采用步骤(2)获得的X₂(t)减去步骤(5)获得的Y₂(t)，实现尖峰噪声的去除，如图7所示。由于采集信号所采用的滤波器是两个二阶级联带通滤波器，故通过用两个二阶级联带通滤波器与脉冲函数的拟合来构造出尖峰噪声的模型，代入特征参数实现对该尖峰的模拟。利用之前定位的位置x₁、x₂...，将采集到的含噪信号与建模得到的尖峰噪声相减来实现尖峰噪声的去除；将截取的信号与原信号拼接到一起便实现了相关算法。

Y(t)＝X₂(t)-Y₂(t) (9)

获得输出MRS信号。

实施例

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

步骤(1)：构建一个仿真的地面核磁共振信号用X₁(t)表示，如图2所示，该信号可以分解成四部分：

X₁(t)＝NMR(t)+h(t)+sp(t)+w(t) (10)

其中NMR(t)表示纯净的NMR信号：

在式(14)中，拉莫尔频率：f_l＝2326Hz，幅值：E₀＝200nV，弛豫时间：T₂＝200nV；h(t)是电力线谐波噪声：

sp(t)代表尖峰噪声，具体参数设置如步骤(3)。

w(t)代表随机噪声，此处添加了5dB的随机白噪声。

步骤(2)：对步骤(1)中获得X₁(t)采用基于谐波建模方法去除工频谐波，谐波建模的识别结果为：中心频率：50.01Hz，得到X₂(t)，如图3，此时：

X₂(t)＝NMR(t)+sp(t)+w(t) (11)

步骤(3)：由于电网的尖峰近似于地面NMR仪器的脉冲激励，故利用两个二阶级联带通滤波器与单位脉冲函数卷积来建立尖峰噪声模型Y₁(t)，则：

sp(t)＝Aδ(t-t₀)*g_rec(t) (12)

其中：A＝9μv、t0＝0.0917s，g_rec(t)是滤波器传递函数，g_rec(t)的参数为σ_a＝150Hz、ω_α＝2148Hz、σ_b＝150Hz、ω_b＝2152Hz；

步骤(4)：对步骤(2)中得到的X₂(t)，首先利用NEO算法进一步强化尖峰：

再通过设计的低通滤波器得到Y_NEO。低通滤波器参数如下：通带截频wp₁＝1800×2/f_s，阻带截频wr₁＝2000×2/f_s，如图4所示。

为得到阈值曲线，根据MRS持续时间t＝250ms设定截取数据段时间t₁＝0.0125s，根据采样率和t₁设置两个数据段之间的时间间隔t₂＝0.00475s,通过t₁、t₂两参数将经过低通滤波器后的NEO信号分为n+1＝51段并存入51*500阶矩阵ME，利用函数求取ME矩阵的每行中位数存入向量y中，在横轴上取51个点，以0.005s为时间间隔，存入向量x。利用matlab自带的Curve Fitting Tool曲线拟合，得到多项式系数向量P,利用系数向量P和阶数n建立该信号的趋势曲线，其表达式：

如图5，定位尖峰所在的位置:得到拟合函数envelop_line之后，乘上一个阈值因子(经验值1)作为阈值函数，选取超过此阈值的部分作为尖峰噪声的定位点x₁、x₂...，利用这些定位点，在X₂(t)对应位置的前后截取27.5ms的数据，以便后续对尖峰噪声处理。

步骤(5)：基于步骤(3)的模型Y₁(t)和步骤(4)定位的位置截取的尖峰数据段，采用最小二乘法提取特征参数A,t₀，从而获得尖峰噪声Y₂(t)，如图6所示。对截取出来的数据段进行最小二乘法处理。

提取特征参数A。可以求取在t₀＝t_l下的A值，记为A_l，存入向量A中。

A＝(M^TM)^-1M^TH (7)

计算t₀＝t_l下的均方根误差(RMSE)，存入向量RMSE_l中。

评价得到最小均方根误差对应的A,t₀为特征参数。参数辨识结果：T₀＝0.0024s,A＝2.1479×10^-5，从而获得尖峰噪声Y₂(t)，其中两参数为t₀0.0024+0.0893＝0.0917s,A＝2.1479×10^-5。如图6所示。

步骤(6)：采用步骤(2)获得的X₂(t)减去步骤(5)获得的Y₂(t)，实现尖峰噪声的去除，如图7所示，其中图7a为为去除工频谐波后的含噪MRS信号，图7b为为去除尖峰噪声后的MRS信号。由于采集信号所采用的滤波器是两个不同二阶带通滤波器级联，故通过用上述滤波器与脉冲函数的拟合来构造出尖峰噪声的模型，代入特征参数实现对该尖峰的模拟，利用之前定位的位置x₁、x₂...，将采集到的含噪信号与建模得到的尖峰噪声相减来实现尖峰噪声的去除；将截取的信号与原信号拼接到一起便实现了相关算法。

Y(t)＝X₂(t)-Y₂(t) (9)

获得输出MRS信号。

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

Claims (3)

1.一种基于模型的地面核磁共振信号尖峰噪声去除方法，其特征在于，该方法包括：

步骤(1)：利用核磁共振测深探水仪采集到一组观测MRS信号X₁(t)；

步骤(2)：对步骤(1)中获得的X₁(t)采用谐波建模方法去除工频谐波，得到X₂(t)；

步骤(3)：利用两个二阶级联带通滤波器与单位脉冲函数卷积建立尖峰噪声模型Y₁(t)，Y₁(t)＝Aδ(t-t₀)*g_rec(t)，A是尖峰噪声的幅值，t₀是尖峰信号的初始时间，g_rec(t)是滤波器传递函数的脉冲响应；

步骤(4)：对步骤(2)中得到的X₂(t)，通过NEO算法进一步强化尖峰，通过设置阈值找出尖峰所在位置x₁、x₂；

步骤(5)：基于步骤(3)的模型Y₁(t)和步骤(4)定位的位置，截取尖峰数据段，采用最小二乘法提取特征参数A，t₀，从而获得尖峰噪声Y₂(t)；

步骤(6)：采用步骤(2)获得的X₂(t)减去步骤(5)获得的Y₂(t)，实现尖峰噪声的去除，获得输出MRS信号；

所述的步骤(3)中建立尖峰噪声模型Y₁(t)＝Aδ(t-t₀)*g_rec(t)中求取g_rec(t)的具体步骤为：

构造四阶带通滤波器的传递函数G(s)＝G_a(s)G_b(s)，其是由两个不同的二阶带通滤波器级联组成的：

其中，s表示复数频率，ω_a为滤波器a的中心频率，σ_a为滤波器a的带宽；ω_b为滤波器b的中心频率，σ_b为滤波器b的带宽；

采用拉普拉斯求逆变换，计算其脉冲响应为：

所述的步骤(5)中的最小二乘法提取特征参数的包括：

限定特征参数t₀，利用最小二乘公式，拟合出特征参数A，通过对比均方根误差：

选择具有最小均方根误差对应的特征参数A作为尖峰噪声模型的幅值参数，而此时的t₀为提取出的尖峰噪声初始时间参数，设采样率为f_s，将建立的模型按f_s采样，采样点数为N，记为向量y_m，将含噪数据段中的数据记为y_r，k＝1,2,...,N表示数据点数；

步骤(5)所述的采用最小二乘法提取特征参数A，t₀的具体步骤为：

51)数据预处理，将含噪数据段中的数据记为M，表示为[m₁...m_S]^T，将建立的模型以同样的采样率处理成S个数据点，记为向量H，表示为[h₁...h_S]^T；

52)提取特征参数A，设t_l是设定的扫描点，t₀为提取出的尖峰噪声起始时间参数；

根据式(7)求取在t₀＝t_l下的A值，记为A_l，存入向量A中；

A＝(M^TM)^-1M^TH (7)

53)利用式(8)计算t₀＝t_l下的均方根误差RMSE，存入向量RMSE中；

54)利用等间隔扫描形式求取特征参数t₀，设置扫描间隔为1/f_s，扫描次数取N/2，在每个t₀＝t_l下进行重复上述51)-53)步骤，获得向量A和RMSE，找出RMSE向量中最小值所对应的A值和t₀值；

步骤(2)所述的谐波建模方法去除工频谐波的具体步骤为：

21)建立谐波噪声模型V_har，V_har是以工频f₀为基频的个谐波的累加结果：

其中，A_n和

分别是第n次谐波的幅度和初始相位；N是谐波个数；

22)通过自适应扫描方式确定基频f₀；

23)求解线性系数，将表达式(2)中的第n项线性化为式(4)：

其中，

由此估计A_n和

的问题转换成求解线性系数α_n和β_n，表达式(4)整理成线性方程组的形式：

其中，V_p(p＝1,2,…P)是t_p时刻的接收数据，表达式(5)整理成矩阵形式为：

Ax＝b (6)

其中x＝[α₁...α_n,β₁...β_n]^T,b＝[v₁,v₂...v_p]^T，表达式(6)是标准的线性方程组，采用最小二乘方法求解；

24)经过上述步骤求出V_har，用V_R减去V_har完成了消噪过程，其中V_R，表示测量数据X₁(t)。

2.按照权利要求1所述的方法，其特征在于，步骤(4)所述的通过NEO算法进一步强化尖峰的具体步骤为：

41)采用NEO算法强化尖峰幅度：

其中x(k)表示测量数据X₁(t)离散化后的结果，经过NEO进行运算之后，再通过低通滤波器，就会使得尖峰信号突出；

42)设定阈值；

43)对尖峰噪声进行定位。

3.按照权利要求2所述的方法，其特征在于，设定阈值包括：

421)根据MRS持续时间T设定截取数据段时间t₁，取t₁为尖峰信号持续时间的1.5-2倍；

422)根据采样率和t₁设置两截取数据段之间的时间间隔t₂，保证t₂的整数n倍与t₁之和等于持续时间T，T＝t₁+n×t₂；

423)通过t₁、t₂两参数将经过低通滤波器后的NEO信号分为n+1段，将每段信号作为一行，将所有行组合成一数据矩阵，用ME表示；

424)利用median函数求取ME矩阵中每行的中位数，并将数据存入向量y中，作为拟合曲线的Y轴；

425)在横轴上取n+1个点，以T/(n+1)为时间间隔，存入向量x中，作为拟合曲线的X轴；

426)利用matlab自带的Curve Fitting Tool，选择向量x为其中的X data，向量y为Ydata，类型为多项式Polynomial，进行曲线拟合，取3阶或4阶，得到多项式系数向量P；

427)利用系数向量P和阶数n建立该信号的趋势曲线，得到拟合函数，其表达式为

n是阶数；

428)拟合函数乘上一个阈值因子作为阈值函数。

Images