本发明属于地面核磁共振(surfacenuclearmagneticresonance，snmr)，又称磁共振测深(magneticresonancesounding,mrs)信号的噪声滤除领域，具体是为一种基于模型的方案从地面核磁共振信号中去除尖峰噪声的方法。
背景技术：
地面核磁共振(surfacenuclearmagneticresonance，snmr)，又称为磁共振测深(magneticresonancesounding,mrs)技术是近年发展起来的一种非侵入式的、直接、定量探测地下水的地球物理方法，其基本原理是通过探测地下水中氢质子共振跃迁产生的核磁共振信号来实现地下水的探测。核磁共振信号的表达式为：
相比于传统的探测地下水的方式，地面核磁共振探测的优点体现在分辨率高、效率高、信息量丰富，能够对地下含水层深度和厚度、含水量的大小、地下介质孔隙度等信息给出定量解释(marianhertrich.imagingofgroundwaterwithnuclearmagneticresonance.progressinnuclearmagneticresonancespectroscopy,2008,53(4):227-248)。
专利号103150707a公开了一种消除磁共振图像中尖峰噪声的方法，能够修复信号强度较高的低频区域的尖峰噪声，专利号cn104835123a公开了一种基于先验模型的光片显微成像条纹噪声去除，均属于图像处理应用领域；专利号cn101523722a公开了一种尖峰噪声消除电路，采用硬件手段直接消除尖峰噪声，属于数字系统及iic总线技术领域；专利号cn104218917a公开了一种用于消除apd雪崩信号输出端的尖峰噪声的自差分滤波装置，属于红外探测领域；专利号cn106452470a公开了一种消除信道尖峰脉冲噪声的无线通信系统及除噪方法，采用单匝环形天线接收尖峰脉冲噪声去抵消多匝环形天线信号中的尖峰噪声，属于无线通信技术领域。目前尚未见基于模型的方法去除尖峰噪声用于地面核磁共振地下水探测领域。
技术实现要素：
本发明所要解决的技术问题在于提供一种基于模型的地面核磁共振信号尖峰噪声去除方法，解决由于尖峰噪声存在难以提取有效地面核磁共振信号的问题。
本发明是这样实现的，
一种基于模型的地面核磁共振信号尖峰噪声去除方法，该方法包括：
步骤(1)：利用核磁共振测深探水仪采集到一组观测mrs信号x1(t)；
步骤(2)：对步骤(1)中获得的x1(t)采用谐波建模方法去除工频谐波，得到x2(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，通过对比均方根误差：
进一步地，步骤(5)所述的利用最小二乘法提取特征参数的具体步骤为：
51)数据预处理，将含噪数据段中的数据记为m，表示为[m1...ms]t，将建立的模型以同样的采样率处理成s个数据点，记为向量h，表示为[h1...hs]t；
根据式(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所述的方法，其特征在于，设定阈值包括：
423)通过t1、t2两参数将经过低通滤波器后的neo信号分为n+1段，将每段信号作为一行，将所有行组合成一数据矩阵，用me表示；
424)利用median函数求取me矩阵中每行的中位数，并将数据存入向量y中，作为拟合曲线的y轴；
426)利用matlab自带的curvefittingtool，选择向量x为其中的xdata，向量y为ydata，类型为多项式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)
谐波建模方法去除工频谐波的具体步骤为：
其中，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组扫描即可达到精度要求。
24)经过上述步骤求出vhar，用vr减去vhar即完成了消噪过程。
步骤(3)：由于电网的尖峰近似于地面nmr仪器的脉冲激励，故利用两个二阶级联带通滤波器与单位脉冲函数卷积来建立尖峰噪声模型y1(t)，则：
sp(t)＝aδ(t-t0)*grec(t)(12)
其中，步骤3中，建立尖峰噪声模型y1(t)＝aδ(t-t0)*grec(t)中求取grec(t)的具体步骤为：
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)然后设置阈值：包括：
413)通过t1、t2两参数将经过低通滤波器后的neo信号分为n+1段，将每段信号作为一行，将所有行组合成一数据矩阵，用me表示。
414)利用median函数求取me矩阵中每行的中位数，并将数据存入向量y中，作为拟合曲线的y轴。
416)利用matlab自带的curvefittingtool，选择向量x为其中的xdata，向量y为ydata，类型为多项式(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中。
计算t0＝tl下的均方根误差(rmse)，存入向量rmsel中。
评价得到最小均方根误差对应的a,t0为特征参数。具体为：
y(t)＝x2(t)-y2(t)(9)
获得输出mrs信号。
实施例
本实施例是在matlab2015b编程环境下开展的本发明方法的仿真实验。
步骤(1)：构建一个仿真的地面核磁共振信号用x1(t)表示，如图2所示，该信号可以分解成四部分：
其中nmr(t)表示纯净的nmr信号：
sp(t)代表尖峰噪声，具体参数设置如步骤(3)。
w(t)代表随机噪声，此处添加了5db的随机白噪声。
步骤(2)：对步骤(1)中获得x1(t)采用基于谐波建模方法去除工频谐波，谐波建模的识别结果为：中心频率：50.01hz，得到x2(t)，如图3，此时：
其中：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所示。
如图5，定位尖峰所在的位置:得到拟合函数envelop_line之后，乘上一个阈值因子(经验值1)作为阈值函数，选取超过此阈值的部分作为尖峰噪声的定位点x1、x2...，利用这些定位点，在x2(t)对应位置的前后截取27.5ms的数据，以便后续对尖峰噪声处理。
步骤(5)：基于步骤(3)的模型y1(t)和步骤(4)定位的位置截取的尖峰数据段，采用最小二乘法提取特征参数a,t0，从而获得尖峰噪声y2(t)，如图6所示。对截取出来的数据段进行最小二乘法处理。
评价得到最小均方根误差对应的a,t0为特征参数。参数辨识结果：t0＝0.0024s,a＝2.1479×10-5，从而获得尖峰噪声y2(t)，其中两参数为t00.0024+0.0893＝0.0917s,a＝2.1479×10-5。如图6所示。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。