CN109885903A - 一种基于模型的地面核磁共振信号尖峰噪声去除方法 - Google Patents
一种基于模型的地面核磁共振信号尖峰噪声去除方法 Download PDFInfo
- Publication number
- CN109885903A CN109885903A CN201910083580.8A CN201910083580A CN109885903A CN 109885903 A CN109885903 A CN 109885903A CN 201910083580 A CN201910083580 A CN 201910083580A CN 109885903 A CN109885903 A CN 109885903A
- Authority
- CN
- China
- Prior art keywords
- data
- spike noise
- signal
- vector
- noise
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Magnetic Resonance Imaging Apparatus (AREA)
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)中,E0表示地面核磁共振信号的初始振幅,其与地层的含水量相关;表示核磁共振探测信号的平均弛豫时间,与含水层的孔隙度相关,其范围一般在30ms-1000ms范围内;ω0表示磁共振信号的拉莫尔角频率,与实验环境的地磁场强度成比例;表示磁共振信号的初始相位,与地下介质的导电性相关。
相比于传统的探测地下水的方式,地面核磁共振探测的优点体现在分辨率高、效率高、信息量丰富,能够对地下含水层深度和厚度、含水量的大小、地下介质孔隙度等信息给出定量解释(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信号X1(t);
步骤(2):对步骤(1)中获得的X1(t)采用谐波建模方法去除工频谐波,得到X2(t);
步骤(3):利用两个二阶级联带通滤波器与单位脉冲函数卷积建立尖峰噪声模型Y1(t),Y1(t)=Aδ(t-t0)*grec(t),A是尖峰噪声的幅值,t0是尖峰信号的初始时间,grec(t)是滤波器传递函数的脉冲响应;
步骤(4):对步骤(2)中得到的X2(t),通过NEO算法进一步强化尖峰,通过设置阈值找出尖峰所在位置x1、x2...;
步骤(5):基于步骤(3)的模型Y1(t)和步骤(4)定位的位置,截取尖峰数据段,采用最小二乘法提取特征参数A,t0,从而获得尖峰噪声Y2(t);
步骤(6):采用步骤(2)获得的X2(t)减去步骤(5)获得的Y2(t),实现尖峰噪声的去除,获得输出MRS信号。
进一步地,所述的步骤(3)中建立尖峰噪声模型Y1(t)=Aδ(t-t0)*grec(t)中求取grec(t)的具体步骤为:
构造四阶带通滤波器的传递函数G(s)=Ga(s)Gb(s),其是由两个不同的二阶带通滤波器级联组成的:
其中,s表示复数频率,ωa为滤波器a的中心频率,σa为滤波器a的带宽;ωb为滤波器b的中心频率,σb为滤波器b的带宽;
采用拉普拉斯求逆变换,计算其脉冲响应为:
进一步地,所述的步骤(5)中的最小二乘法提取特征参数的包括:
限定特征参数t0,利用最小二乘公式,拟合出特征参数A,通过对比均方根误差:
选择具有最小均方根误差对应的特征参数A作为尖峰噪声模型的幅值参数,而此时的t0为提取出的尖峰噪声初始时间参数,设采样率为fs,将建立的模型按fs采样,采样点数为N,记为向量ym,将含噪数据段中的数据记为yr,k=1,2,...,N表示数据点数。
进一步地,步骤(5)所述的利用最小二乘法提取特征参数的具体步骤为:
51)数据预处理,将含噪数据段中的数据记为M,表示为[m1...mS]T,将建立的模型以同样的采样率处理成S个数据点,记为向量H,表示为[h1...hS]T;
52)提取特征参数A,设tl是设定的扫描点,t0为提取出的尖峰噪声起始时间参数;
根据式(7)求取在t0=tl下的A值,记为Al,存入向量A中;
A=(MTM)-1MTH (7)
53)利用式(8)计算t0=tl下的均方根误差RMSE,存入向量RMSE中;
54)利用等间隔扫描形式求取特征参数t0,设置扫描间隔为1/fs,扫描次数取N/2,在每个t0=tl下进行重复上述51)-53)步骤,获得向量A和RMSE,找出RMSE向量中最小值所对应的A值和t0值。
进一步地,步骤(2)所述的谐波建模方法去除工频谐波的具体步骤为:
21)建立谐波噪声模型Vhar,Vhar是以工频f0为基频的个谐波的累加结果:
其中,An和分别是第n次谐波的幅度和初始相位;N是谐波个数;
22)通过自适应扫描方式确定基频f0;
23)求解线性系数,将表达式(2)中的第n项线性化为式(4):
其中,由此估计An和的问题转换成求解线性系数αn和βn,表达式(4)整理成线性方程组的形式:
其中,Vp(p=1,2,…P)是tp时刻的接收数据,表达式(5)整理成矩阵形式为:
Ax=b (6)
其中x=[α1...αn,β1...βn]T,b=[v1,v2...vp]T,表达式(6)是标准的线性方程组,采用最小二乘方法求解。
24)经过上述步骤求出Vhar,用VR减去Vhar完成了消噪过程,其中VR,表示测量数据X1(t)。
进一步地,步骤(4)所述的利用NEO算法进一步强化尖峰噪声的具体步骤为:
41)采用NEO算法强化尖峰幅度:
其中x(k)表示测量数据X1(t)离散化后的结果,经过NEO进行运算之后,再通过低通滤波器,就会使得尖峰信号突出;
42)设定阈值;
43)对尖峰噪声进行定位。
7.按照权利要求6所述的方法,其特征在于,设定阈值包括:
421)根据MRS持续时间T设定截取数据段时间t1,取t1为尖峰信号持续时间的1.5-2倍;
422)根据采样率和t1设置两截取数据段之间的时间间隔t2,保证t2的整数n倍与t1之和等于持续时间T,T=t1+n×t2;
423)通过t1、t2两参数将经过低通滤波器后的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信号X1(t);
图3为本发明实施例提供的去除工频谐波的信号X2(t);
图4为本发明实施例提供的NEO算法强化后提取MRS信号包络与低通滤波后结果;
图5为本发明实施例提供的阈值曲线;
图6为本发明实施例提供的参数辨识;
图7为本发明实施例提供的实现尖峰噪声的去除图,图7a为为去除工频谐波后的含噪MRS信号,图7b为为去除尖峰噪声后的MRS信号。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
参见图1所示,一种基于模型的地面核磁共振信号尖峰噪声去除方法,包括如下的步骤
步骤(1):利用核磁共振测深(MRS)探水仪采集到一组X1(t)观测MRS信号,如图2所示,核磁信号仪器以采样率fs采集接收线圈中的MRS信号,该信号可以分解成四部分:
X1(t)=NMR(t)+h(t)+sp(t)+w(t) (10)
其中NMR(t)表示MRS信号,h(t)是电力线谐波噪声,sp(t)代表尖峰噪声,w(t)代表随机噪声。
步骤(2):对步骤(1)中获得X1(t)采用基于谐波建模方法去除工频谐波,得到X2(t),如图3,此时:
X2(t)=NMR(t)+sp(t)+w(t) (11)
谐波建模方法去除工频谐波的具体步骤为:
21)建立谐波噪声模型Vhar,Vhar是以工频f0为基频的个谐波的累加结果:
其中,An和分别是第n次谐波的幅度和初始相位;N是谐波个数;根据接收信号的带宽(1kHz~3kHz),取N=100。
22)通过自适应扫描方式确定基频f0;首先在49.9~50.1Hz的范围内以0.03Hz为步长进行粗扫,利用表达式(3)获得2阶范数最小的3个频率点作为下一次扫描的范围。
||VR-Vhar(f1)||2→min (3)
其中,VR为测量数据即X1(t),Vhar(fl)为第l个扫描值。在第2组扫描中,扫描范围为49.99~50.05Hz,扫描步长为0.0075Hz,再次获得范数最小的3个频点。当第3组扫描时,扫描范围缩小为50.0125~50.0275Hz,扫描步长为0.001875Hz。此时2阶范数的最小值对应的频点为50.018125Hz,即经过3组扫描即可达到精度要求。
23)求解线性系数,将表达式(2)中的第n项线性化为式(4):
其中,由此估计An和的问题转换成求解线性系数αn和βn,表达式(4)整理成线性方程组的形式:
其中,Vp(p=1,2,…P)是tp时刻的接收数据,表达式(5)整理成矩阵形式为:
Ax=b (6)
其中x=[α1...αn,β1...βn]T,b=[v1,v2...vp]T,表达式(6)是标准的线性方程组,采用最小二乘方法求解。
24)经过上述步骤求出Vhar,用VR减去Vhar即完成了消噪过程。
步骤(3):由于电网的尖峰近似于地面NMR仪器的脉冲激励,故利用两个二阶级联带通滤波器与单位脉冲函数卷积来建立尖峰噪声模型Y1(t),则:
sp(t)=Aδ(t-t0)*grec(t) (12)
其中,A是尖峰噪声的幅值,t0是尖峰信号的初始时间,grec(t)是滤波器传递函数。
其中,步骤3中,建立尖峰噪声模型Y1(t)=Aδ(t-t0)*grec(t)中求取grec(t)的具体步骤为:
构造四阶带通滤波器的传递函数G(s)=Ga(s)Gb(s),其是由两个不同的二阶带通滤波器级联组成的:
s表示复数频率,ωa为滤波器a的中心频率,σa为滤波器a的带宽;ωb为滤波器b的中心频率,σb为滤波器b的带宽;
采用拉普拉斯求逆变换,计算其脉冲响应为:
步骤(4):对步骤(2)中得到的X2(t),通过NEO算法进一步强化尖峰,如图4所示,41)首先利用NEO算法强化尖峰幅度:
其中x(k)表示测量数据X1(t)离散化后的结果,t=k/fs,k=1,2,….,N-1,经过NEO进行运算之后,再通过适当的低通滤波器,就会使得尖峰信号突出;
42)然后设置阈值:包括:
411)根据MRS持续时间T设定截取数据段时间t1,一般取t1为尖峰信号持续时间的1.5-2倍;
412)根据采样率和t1设置两截取数据段之间的时间间隔t2,应保证t2的整数n倍与t1之和等于持续时间T,即T=t1+n×t2。
413)通过t1、t2两参数将经过低通滤波器后的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所示),选取超过此阈值的部分作为尖峰噪声的定位点x1、x2...,利用这些定位点,在X2(t)对应位置的前后截取27.5ms的数据,以便后续对尖峰噪声处理。
步骤(5):基于步骤(3)的模型Y1(t)和步骤(4)定位的位置截取的尖峰数据段,采用最小二乘法提取特征参数A,t0,从而获得尖峰噪声Y2(t),如图6所示。对截取出来的数据进行最小二乘法处理。
提取特征参数A。可以求取在t0=tl下的A值,记为Al,存入向量A中。
A=(MTM)-1MTH (7)
计算t0=tl下的均方根误差(RMSE),存入向量RMSEl中。
评价得到最小均方根误差对应的A,t0为特征参数。具体为:
51)数据预处理,将含噪数据段中的数据记为M,表示为[m1...mS]T,将建立的模型以同样的采样率处理成S个数据点,记为向量H,表示为[h1...hS]T;
52)提取特征参数A,设tl是设定的扫描点,t0为提取出的尖峰噪声起始时间参数;
根据式(7)求取在t0=tl下的A值,记为Al,存入向量A中;
A=(MTM)-1MTH (7)
53)利用式(8)计算t0=tl下的均方根误差RMSE,存入向量RMSE中;
54)利用等间隔扫描形式求取特征参数t0,设置扫描间隔为1/fs,扫描次数取N/2,在每个t0=tl下进行重复上述51)-53)步骤,获得向量A和RMSE,找出RMSE向量中最小值所对应的A值和t0值。
步骤(6):采用步骤(2)获得的X2(t)减去步骤(5)获得的Y2(t),实现尖峰噪声的去除,如图7所示。由于采集信号所采用的滤波器是两个二阶级联带通滤波器,故通过用两个二阶级联带通滤波器与脉冲函数的拟合来构造出尖峰噪声的模型,代入特征参数实现对该尖峰的模拟。利用之前定位的位置x1、x2...,将采集到的含噪信号与建模得到的尖峰噪声相减来实现尖峰噪声的去除;将截取的信号与原信号拼接到一起便实现了相关算法。
Y(t)=X2(t)-Y2(t) (9)
获得输出MRS信号。
实施例
本实施例是在MATLAB2015b编程环境下开展的本发明方法的仿真实验。
步骤(1):构建一个仿真的地面核磁共振信号用X1(t)表示,如图2所示,该信号可以分解成四部分:
X1(t)=NMR(t)+h(t)+sp(t)+w(t) (10)
其中NMR(t)表示纯净的NMR信号:
在式(14)中,拉莫尔频率:fl=2326Hz,幅值:E0=200nV,弛豫时间:T2=200nV;h(t)是电力线谐波噪声:
sp(t)代表尖峰噪声,具体参数设置如步骤(3)。
w(t)代表随机噪声,此处添加了5dB的随机白噪声。
步骤(2):对步骤(1)中获得X1(t)采用基于谐波建模方法去除工频谐波,谐波建模的识别结果为:中心频率:50.01Hz,得到X2(t),如图3,此时:
X2(t)=NMR(t)+sp(t)+w(t) (11)
步骤(3):由于电网的尖峰近似于地面NMR仪器的脉冲激励,故利用两个二阶级联带通滤波器与单位脉冲函数卷积来建立尖峰噪声模型Y1(t),则:
sp(t)=Aδ(t-t0)*grec(t) (12)
其中:A=9μv、t0=0.0917s,grec(t)是滤波器传递函数,grec(t)的参数为σa=150Hz、ωα=2148Hz、σb=150Hz、ωb=2152Hz;
步骤(4):对步骤(2)中得到的X2(t),首先利用NEO算法进一步强化尖峰:
再通过设计的低通滤波器得到YNEO。低通滤波器参数如下:通带截频wp1=1800×2/fs,阻带截频wr1=2000×2/fs,如图4所示。
为得到阈值曲线,根据MRS持续时间t=250ms设定截取数据段时间t1=0.0125s,根据采样率和t1设置两个数据段之间的时间间隔t2=0.00475s,通过t1、t2两参数将经过低通滤波器后的NEO信号分为n+1=51段并存入51*500阶矩阵ME,利用函数求取ME矩阵的每行中位数存入向量y中,在横轴上取51个点,以0.005s为时间间隔,存入向量x。利用matlab自带的Curve Fitting Tool曲线拟合,得到多项式系数向量P,利用系数向量P和阶数n建立该信号的趋势曲线,其表达式:
如图5,定位尖峰所在的位置:得到拟合函数envelop_line之后,乘上一个阈值因子(经验值1)作为阈值函数,选取超过此阈值的部分作为尖峰噪声的定位点x1、x2...,利用这些定位点,在X2(t)对应位置的前后截取27.5ms的数据,以便后续对尖峰噪声处理。
步骤(5):基于步骤(3)的模型Y1(t)和步骤(4)定位的位置截取的尖峰数据段,采用最小二乘法提取特征参数A,t0,从而获得尖峰噪声Y2(t),如图6所示。对截取出来的数据段进行最小二乘法处理。
提取特征参数A。可以求取在t0=tl下的A值,记为Al,存入向量A中。
A=(MTM)-1MTH (7)
计算t0=tl下的均方根误差(RMSE),存入向量RMSEl中。
评价得到最小均方根误差对应的A,t0为特征参数。参数辨识结果:T0=0.0024s,A=2.1479×10-5,从而获得尖峰噪声Y2(t),其中两参数为t00.0024+0.0893=0.0917s,A=2.1479×10-5。如图6所示。
步骤(6):采用步骤(2)获得的X2(t)减去步骤(5)获得的Y2(t),实现尖峰噪声的去除,如图7所示,其中图7a为为去除工频谐波后的含噪MRS信号,图7b为为去除尖峰噪声后的MRS信号。由于采集信号所采用的滤波器是两个不同二阶带通滤波器级联,故通过用上述滤波器与脉冲函数的拟合来构造出尖峰噪声的模型,代入特征参数实现对该尖峰的模拟,利用之前定位的位置x1、x2...,将采集到的含噪信号与建模得到的尖峰噪声相减来实现尖峰噪声的去除;将截取的信号与原信号拼接到一起便实现了相关算法。
Y(t)=X2(t)-Y2(t) (9)
获得输出MRS信号。
以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。
Claims (7)
1.一种基于模型的地面核磁共振信号尖峰噪声去除方法,其特征在于,该方法包括:
步骤(1):利用核磁共振测深探水仪采集到一组观测MRS信号X1(t);
步骤(2):对步骤(1)中获得的X1(t)采用谐波建模方法去除工频谐波,得到X2(t);
步骤(3):利用两个二阶级联带通滤波器与单位脉冲函数卷积建立尖峰噪声模型Y1(t),Y1(t)=Aδ(t-t0)*grec(t),A是尖峰噪声的幅值,t0是尖峰信号的初始时间,grec(t)是滤波器传递函数的脉冲响应;
步骤(4):对步骤(2)中得到的X2(t),通过NEO算法进一步强化尖峰,通过设置阈值找出尖峰所在位置x1、x2...;
步骤(5):基于步骤(3)的模型Y1(t)和步骤(4)定位的位置,截取尖峰数据段,采用最小二乘法提取特征参数A,t0,从而获得尖峰噪声Y2(t);
步骤(6):采用步骤(2)获得的X2(t)减去步骤(5)获得的Y2(t),实现尖峰噪声的去除,获得输出MRS信号。
2.按照权利要求1所述的方法,其特征在于,所述的步骤(3)中建立尖峰噪声模型Y1(t)=Aδ(t-t0)*grec(t)中求取grec(t)的具体步骤为:
构造四阶带通滤波器的传递函数G(s)=Ga(s)Gb(s),其是由两个不同的二阶带通滤波器级联组成的:
其中,s表示复数频率,ωa为滤波器a的中心频率,σa为滤波器a的带宽;ωb为滤波器b的中心频率,σb为滤波器b的带宽;
采用拉普拉斯求逆变换,计算其脉冲响应为:
3.按照权利要求1所述的方法,其特征在于,所述的步骤(5)中的最小二乘法提取特征参数的包括:
限定特征参数t0,利用最小二乘公式,拟合出特征参数A,通过对比均方根误差:
选择具有最小均方根误差对应的特征参数A作为尖峰噪声模型的幅值参数,而此时的t0为提取出的尖峰噪声初始时间参数,设采样率为fs,将建立的模型按fs采样,采样点数为N,记为向量ym,将含噪数据段中的数据记为yr,k=1,2,...,N表示数据点数。
4.按照权利要求3所述的方法,其特征在于,步骤(5)所述的利用最小二乘法提取特征参数的具体步骤为:
51)数据预处理,将含噪数据段中的数据记为M,表示为[m1...mS]T,将建立的模型以同样的采样率处理成S个数据点,记为向量H,表示为[h1...hS]T;
52)提取特征参数A,设tl是设定的扫描点,t0为提取出的尖峰噪声起始时间参数;
根据式(7)求取在t0=tl下的A值,记为Al,存入向量A中;
A=(MTM)-1MTH (7)
53)利用式(8)计算t0=tl下的均方根误差RMSE,存入向量RMSE中;
54)利用等间隔扫描形式求取特征参数t0,设置扫描间隔为1/fs,扫描次数取N/2,在每个t0=tl下进行重复上述51)-53)步骤,获得向量A和RMSE,找出RMSE向量中最小值所对应的A值和t0值。
5.按照权利要求1所述的方法,其特征在于,步骤(2)所述的谐波建模方法去除工频谐波的具体步骤为:
21)建立谐波噪声模型Vhar,Vhar是以工频f0为基频的个谐波的累加结果:
其中,An和分别是第n次谐波的幅度和初始相位;N是谐波个数;
22)通过自适应扫描方式确定基频f0;
23)求解线性系数,将表达式(2)中的第n项线性化为式(4):
其中,由此估计An和的问题转换成求解线性系数αn和βn,表达式(4)整理成线性方程组的形式:
其中,Vp(p=1,2,…P)是tp时刻的接收数据,表达式(5)整理成矩阵形式为:
Ax=b (6)
其中x=[α1...αn,β1...βn]T,b=[v1,v2...vp]T,表达式(6)是标准的线性方程组,采用最小二乘方法求解。
24)经过上述步骤求出Vhar,用VR减去Vhar完成了消噪过程,其中VR,表示测量数据X1(t)。
6.按照权利要求1所述的方法,其特征在于,步骤(4)所述的利用NEO算法进一步强化尖峰噪声的具体步骤为:
41)采用NEO算法强化尖峰幅度:
其中x(k)表示测量数据X1(t)离散化后的结果,经过NEO进行运算之后,再通过低通滤波器,就会使得尖峰信号突出;
42)设定阈值;
43)对尖峰噪声进行定位。
7.按照权利要求6所述的方法,其特征在于,设定阈值包括:
421)根据MRS持续时间T设定截取数据段时间t1,取t1为尖峰信号持续时间的1.5-2倍;
422)根据采样率和t1设置两截取数据段之间的时间间隔t2,保证t2的整数n倍与t1之和等于持续时间T,T=t1+n×t2;
423)通过t1、t2两参数将经过低通滤波器后的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)拟合函数乘上一个阈值因子作为阈值函数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910083580.8A CN109885903B (zh) | 2019-01-29 | 2019-01-29 | 一种基于模型的地面核磁共振信号尖峰噪声去除方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910083580.8A CN109885903B (zh) | 2019-01-29 | 2019-01-29 | 一种基于模型的地面核磁共振信号尖峰噪声去除方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109885903A true CN109885903A (zh) | 2019-06-14 |
CN109885903B CN109885903B (zh) | 2022-07-08 |
Family
ID=66927239
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910083580.8A Active CN109885903B (zh) | 2019-01-29 | 2019-01-29 | 一种基于模型的地面核磁共振信号尖峰噪声去除方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109885903B (zh) |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110542926A (zh) * | 2019-09-02 | 2019-12-06 | 吉林大学 | 一种地震数据尖峰噪声簇的自主检测和压制方法 |
CN110542925A (zh) * | 2019-09-02 | 2019-12-06 | 吉林大学 | 一种基于峰值包络线的地震数据尖峰噪声识别及压制方法 |
CN112816510A (zh) * | 2021-03-03 | 2021-05-18 | 明峰医疗***股份有限公司 | Ct扫描设备的探测信号处理方法、***及计算机可读存储介质 |
CN112834844A (zh) * | 2020-12-31 | 2021-05-25 | 瑞斯康微电子(深圳)有限公司 | 一种消除避雷器漏电流信号中尖峰奇异信号的方法 |
CN113361436A (zh) * | 2021-06-16 | 2021-09-07 | 中国农业大学 | 采用一阶导数与对抗网络的信号自动识别方法 |
CN114252899A (zh) * | 2022-03-02 | 2022-03-29 | 四川新先达测控技术有限公司 | 一种核信号的级联冲激卷积成形方法和装置 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5939884A (en) * | 1996-04-24 | 1999-08-17 | Siemens Aktiengesellschaft | Method for identifying spikes in MR signals |
CN104898172A (zh) * | 2015-05-19 | 2015-09-09 | 吉林大学 | 一种基于互相关的核磁共振全波信号噪声滤除方法 |
CN106164694A (zh) * | 2014-03-31 | 2016-11-23 | 皇家飞利浦有限公司 | 具有rf噪声检测线圈的磁共振成像 |
CN106646637A (zh) * | 2016-12-27 | 2017-05-10 | 吉林大学 | 一种去除核磁信号中尖峰噪声的方法 |
US9702953B1 (en) * | 2012-08-10 | 2017-07-11 | University Of New Brunswick | Method of sampling in pure phase encode magnetic resonance imaging |
CN107957566A (zh) * | 2017-11-17 | 2018-04-24 | 吉林大学 | 基于频率选择奇异谱分析的磁共振测深信号提取方法 |
-
2019
- 2019-01-29 CN CN201910083580.8A patent/CN109885903B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5939884A (en) * | 1996-04-24 | 1999-08-17 | Siemens Aktiengesellschaft | Method for identifying spikes in MR signals |
US9702953B1 (en) * | 2012-08-10 | 2017-07-11 | University Of New Brunswick | Method of sampling in pure phase encode magnetic resonance imaging |
CN106164694A (zh) * | 2014-03-31 | 2016-11-23 | 皇家飞利浦有限公司 | 具有rf噪声检测线圈的磁共振成像 |
CN104898172A (zh) * | 2015-05-19 | 2015-09-09 | 吉林大学 | 一种基于互相关的核磁共振全波信号噪声滤除方法 |
CN106646637A (zh) * | 2016-12-27 | 2017-05-10 | 吉林大学 | 一种去除核磁信号中尖峰噪声的方法 |
CN107957566A (zh) * | 2017-11-17 | 2018-04-24 | 吉林大学 | 基于频率选择奇异谱分析的磁共振测深信号提取方法 |
Non-Patent Citations (2)
Title |
---|
RAFIK SOLTANI 等: "Advancing the surface-NMR technology for hydrogeophysical exploration", 《2009 INTERNATIONAL CONFERENCE ON INNOVATIONS IN INFORMATION TECHNOLOGY (IIT)》 * |
李同 等: "基于标准差中位数的磁共振信号尖峰去除方法", 《吉林大学学报(信息科学版)》 * |
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110542926A (zh) * | 2019-09-02 | 2019-12-06 | 吉林大学 | 一种地震数据尖峰噪声簇的自主检测和压制方法 |
CN110542925A (zh) * | 2019-09-02 | 2019-12-06 | 吉林大学 | 一种基于峰值包络线的地震数据尖峰噪声识别及压制方法 |
CN110542926B (zh) * | 2019-09-02 | 2020-07-28 | 吉林大学 | 一种地震数据尖峰噪声簇的自主检测和压制方法 |
CN112834844A (zh) * | 2020-12-31 | 2021-05-25 | 瑞斯康微电子(深圳)有限公司 | 一种消除避雷器漏电流信号中尖峰奇异信号的方法 |
CN112834844B (zh) * | 2020-12-31 | 2023-10-20 | 瑞斯康微电子(深圳)有限公司 | 一种消除避雷器漏电流信号中尖峰奇异信号的方法 |
CN112816510A (zh) * | 2021-03-03 | 2021-05-18 | 明峰医疗***股份有限公司 | Ct扫描设备的探测信号处理方法、***及计算机可读存储介质 |
CN112816510B (zh) * | 2021-03-03 | 2022-06-21 | 明峰医疗***股份有限公司 | Ct扫描设备的探测信号处理方法、***及计算机可读存储介质 |
CN113361436A (zh) * | 2021-06-16 | 2021-09-07 | 中国农业大学 | 采用一阶导数与对抗网络的信号自动识别方法 |
CN114252899A (zh) * | 2022-03-02 | 2022-03-29 | 四川新先达测控技术有限公司 | 一种核信号的级联冲激卷积成形方法和装置 |
CN114252899B (zh) * | 2022-03-02 | 2022-05-20 | 四川新先达测控技术有限公司 | 一种核信号的级联冲激卷积成形方法和装置 |
Also Published As
Publication number | Publication date |
---|---|
CN109885903B (zh) | 2022-07-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109885903A (zh) | 一种基于模型的地面核磁共振信号尖峰噪声去除方法 | |
CN107045149B (zh) | 一种基于双奇异值分解的全波核磁共振信号噪声滤除方法 | |
CN104459809B (zh) | 一种基于独立成分分析的全波核磁共振信号噪声滤除方法 | |
CN106772646A (zh) | 一种地面核磁共振信号提取方法 | |
CN106771905B (zh) | 一种适用于高频电流局部放电检测的脉冲提取方法 | |
CN107219555B (zh) | 基于主成分分析的并行震源地震勘探资料强工频噪声压制方法 | |
CN106680874B (zh) | 基于波形形态特征稀疏化建模的谐波噪声压制方法 | |
CN104614778A (zh) | 基于ica的核磁共振地下水探测信号噪声消除方法 | |
CN106646637A (zh) | 一种去除核磁信号中尖峰噪声的方法 | |
CN104237820B (zh) | 一种单扫描获取磁共振二维j‑分解谱的方法 | |
CN104849757B (zh) | 消除地震信号中随机噪声***及方法 | |
CN106814402B (zh) | 瞬变电磁信号叠前去噪方法 | |
CN107179550B (zh) | 一种数据驱动的地震信号零相位反褶积方法 | |
CN103293556A (zh) | 一种地磁异动监测***及方法 | |
CN108345039A (zh) | 一种消除地面核磁共振数据中邻频谐波干扰的方法 | |
CN104133248B (zh) | 一种高保真声波干扰压制方法 | |
CN103837892B (zh) | 微地震自适应独立分量分析弱信号提取方法 | |
CN106950597B (zh) | 基于三边滤波的混合震源数据分离方法 | |
CN103344988B (zh) | 基于k-l分解的可控震源信号相位检测方法 | |
CN103995293B (zh) | 磁共振测深信号检测方法 | |
CN108983321B (zh) | 一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法 | |
CN103941280B (zh) | 基于冲激响应不变法的数字核脉冲高斯成形方法 | |
CN109782363A (zh) | 一种基于时域建模与频域对称的磁共振信号消噪方法 | |
Lin et al. | Removal of a series of spikes from magnetic resonance sounding signal by combining empirical mode decomposition and wavelet thresholding | |
CN115826068A (zh) | 一种基于自适应高斯滤波的mrs信号包络提取方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |