一种地面核磁共振信号提取方法与流程

本发明涉及地面核磁共振地下水探测信号噪声滤除和参数提取技术领域，具体是利用基于工频谐波建模和自相关的地面核磁共振信号提取方法。

背景技术：

地面核磁共振(Magnetic Resonance Sounding,MRS)作为一种地球物理新技术，近年来被广泛应用到地下水探测中。与传统的地球物理间接找水方法相比，地面核磁共振方法具有高分辨率，高效率，信息丰富和对水文参数具有唯一解释的特点。地面核磁共振方法的基本原理是通过探测地下水中氢质子共振跃迁产生的MRS信号，提取其特征参数，再通过反演解释来获取地下水的水质信息。从一些国内外研究和应用中可以看出，MRS方法不仅在寻找地下水和地下水资源评价等方面取得了良好的应用效果，而且在获取水文地质参数和定量化解释方面也取得进展。

MRS方法虽然比其他传统方法具有明显优势，但也存在一些不足。实际探测到的MRS信号是非常微弱的纳伏级信号，信号中噪声干扰非常大，包括尖峰噪声、随机噪声和工频谐波噪声等，其中由电力线干扰引起的工频谐波噪声对其影响最大。目前，针对滤除核磁共振信号中噪声的理论和方法有很多。訾彦勇在著作《基于EMD的磁共振探测信号噪声抑制方法》，提出一种新型的信号时频处理方法，可以在没有输入信号任何先验知识的情况下，自适应地将信号分解成若干个固有模态函数(Intrinsic Mode Functions,IMF)，实现信号趋势的有效提取，但存在端点效应和模态混叠问题；田宝凤等人在著作《基于参考线圈和变步长自适应的磁共振信号噪声压制方法》中直接利用观测数据不断递归更新处理参数，实现不同信噪比和信号强度下MRS信号噪声的去除，但该方法需要远端线圈同时采集噪声数据做参考，对于单通道仪器无法实现，而且只能对相关噪声处理，不能保证将多个工频谐波噪声全部消除。

CN104459809A公开了“基于独立成分分析的全波核磁共振信号噪声滤除方法”，该算法能够有效地实现全波MRS信号的信噪分离，且数据拟合后初始振幅和弛豫时间的相对误差小于±5.00％，相比于其他经典算法该算法更具有消噪性能的优越性，但该方法只针对去除磁共振信号中的工频谐波噪声。CN104898172A公开了“一种基于互相关的核磁共振全波信号噪声滤除方法”，该发明利用噪声与拉莫尔频率的正弦信号不相关，而MRS幅度衰减正弦信号与拉莫尔频率的正弦信号具有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波形的包络，并重构不含噪声的互相关波形，最后利用解卷积算法提取核磁共振全波数据中的MRS信号。该方法运算数据计算量小，可以同时压制工频及其谐波噪声、随机噪声和尖峰噪声，明显提高核磁共振全波数据信噪比。但是在地磁场不稳定时或信噪比较低时，拉莫尔频率不能准确获得，利用存在频率偏差的信号进行互相关计算后，信号参数提取误差大。CN104636609A中公开了“一种基于经验模态分解与小波分析的信号联合去噪方法”，根据信号的自相关性，对混有高斯白噪声的信号进行EMD分解，由于EMD分解的性质，高斯白噪声已不再是真正的白噪声，但白噪声的统计特性近似存在，即所述混有高斯白噪声的信号的自相关函数在零点取得最大值，幅值随着时间差的变化而变化，但其随着时间的衰减很快。利用这种差异可以选取出噪声起主导作用的IMF分量有效降低噪声对信号的影响。但该方法主要是根据自相关特性对EMD分解后每个模态的属性进行判断，主要用于低信噪比下去除高斯白噪声，所以其应用存在局限性。

技术实现要素：

本发明所要解决的技术问题在于提供一种地面核磁共振信号提取方法，解决了磁共振测深找水工作中由于强工频谐波干扰和随机噪声造成的MRS信号有效提取的问题。

本发明是这样实现的，一种地面核磁共振信号提取方法，该方法包括以下步骤：

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

步骤(2):利用统计方法判断是否存在尖峰噪声，如果存在，去除尖峰噪声并用插值结果代替，如果不存在，则测量数据保持不变；

步骤(3):将去除尖峰噪声的数据，利用谐波建模的方法去除工频谐波噪声；

步骤(4):对步骤(3)的结果进行自相关和叠加处理减小随机噪声；

步骤(5):对步骤(4)后的结果进行MRS信号参数提取。

进一步地，步骤(3)中的建模方法的具体步骤为：

步骤3a：设定工频谐波基频的搜索范围以及步长进行粗扫；

步骤3b：将工频谐波模型转化为Ax＝b的形式，解包含测量数据的线性方程组，得到工频各谐波频点的系数矩阵；

步骤3c：计算工频谐波基频搜索范围内不同扫描值的模型估计值；

步骤3d：计算模型估计值与测量数据差值的2范数；

步骤3e：确定使误差2范数取得最小值的频点，从而获得两侧相邻的频点；

步骤3f：将步骤3e中得到的两侧相邻频点作为下一次扫描范围，选取适当步长再次扫描；

步骤3g：判断是否扫描设定的次数，如果达到则根据误差2范数最小的频点进行工频谐波建模；否则重复步骤3c～3g；

步骤3h:从步骤(2)的结果中减去工频谐波模型，去除工频谐波噪声。

进一步地，步骤3f：步长为(f_m2-f_m1)/M，其中f_m1，f_m2分别为最小值频点的两侧相邻频点，其中M＝3～7。

进一步地，步骤3g扫描设定的扫描次数为3-7次。

进一步地，所述的步骤(4)中对步骤(3)得到的结果进行自相关和叠加处理减小随机噪声的具体步骤为：

步骤4a：将MRS信号代入自相关函数公式，修改积分上下限，求得MRS信号的自相关表达式；

步骤4b：判断采集时间是否大于1s，如果大于1s截取自相关结果的前半段用近似，其中e₀为初始振幅、为弛豫时间、f为拉莫尔频率和τ为自相关时间间隔；如果小于1s自相0关结果用近似，其中t_max为时间最大值；

步骤4c：基于步骤4b的判断后将多次测量数据的自相关结果进行叠加，进一步减小随机噪声。

进一步地，所述的步骤(5)中对步骤(4)得到的结果进行参数提取的具体步骤为：

步骤5a:对步骤(4)中叠加的数据进行希尔伯特变换，再通过低通滤波，转化为两个正交分量；

步骤5b:利用非线性拟合方法求解MRS信号的特征参数A₀和弛豫时间其中再计算初始振幅e₀结果，参数提取结果为初始振幅e₀和弛豫时间

本发明与现有技术相比，有益效果在于：本发明提出了基于工频谐波建模和自相关的地面核磁共振信号提取方法，针对单通道采集的全波MRS数据，不仅可以一次性去除所有谐波干扰，而且还结合了自相关方法压制随机噪声，最后通过非线性拟合实现了有效MRS信号特征参数的提取。本发明方法解决了磁共振测深找水工作中由于强工频谐波干扰和随机噪声造成的MRS信号有效提取的难题，获取的MRS信号关键特征参数的拟合误差较小，同时本发明突破了经典消噪方法需多通道探测等其他条件的限制，节省了大量的人力物力。

附图说明

图1为本发明地面核磁共振信号提取方法流程框图；

图2为工频谐波建模消噪原理及效果图，图2(a)是其时域信号，MRS信号淹没在噪声之中；图2(b)是其功率谱，在谐波的频点上均存在不同幅度的谐波噪声；图2(c)谐波建模后，在测量数据b中减去V_harmonic完成消噪过程，2(d)谐波建模消噪后，MRS信号的轮廓呈现衰减趋势，在功率谱中只存在MRS信号的频点；

图3为自相关消噪效果图；图3(a)中黑色线是理想的MRS信号，参数为e₀＝100nV，f＝2326Hz和灰色线是加入了100nV均方差的高斯白噪声后的测量信号，此时数据的信噪比为0dB，MRS信号基本淹没在噪声中；图3(b)是测量信号和MRS信号的频谱，可见MRS信号的频谱峰值是噪声频谱均值的10倍。3(c)经过自相关函数后，测量信号(灰色线)和MRS信号(黑色线)的结果；3(d)MRS信号自相关后的频谱峰值与噪声频谱；

图4为自相关后非线性拟合效果图，图4(a)是未经过自相关函数的曲线，图4(b)是经过自相关函数的曲线；

图5为仿真数据和数据处理结果图，图5(a)含噪MRS信号；图5(b)从含噪MRS信号中减去工频谐波模型，得到去工频信号图；5(c)自相关和叠加处理得到的消噪后信号图，5(d)MRS信号的实部，5(e)MRS信号的虚部，5(f)MRS信号实部的叠加曲线和虚部见；5(g)MRS信号虚部的叠加曲线；

图6为实测数据处理结果图；图6(a-1)为q＝1.4As的MRS信号；6(a-2)为q＝4.2As的MRS信号；6(a-3)为q＝8.6As的MRS信号；

图6(b-1)为q＝1.4As的信噪比信号；6(b-2)为q＝4.2As的的信噪比信号；6(b-3)为q＝8.6As的的信噪比信号。

具体实施方式

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

如图1所示，一种基于工频谐波建模和自相关的地面核磁共振信号提取方法，包括以下步骤：

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

步骤(2):利用统计方法(如罗曼诺夫斯基准则和能量运算等算法)判断是否存在尖峰噪声，如果存在，去除尖峰噪声并用插值结果代替，如果不存在，则测量数据保持不变；

步骤(3):利用谐波建模的方法去除工频谐波噪声；

步骤(4):对步骤(3)的结果进行自相关和叠加处理减小随机噪声；

步骤(5):对步骤(4)后的结果进行MRS信号参数提取；

如图2所示，步骤(3)中的工频建模方法的具体步骤：

步骤3a:设定工频谐波基频f₀的搜索范围49.9Hz～50.1Hz，f₀以步长为0.03Hz进行粗扫；

步骤3b:工频谐波噪声V_harmonic可表示为式(1)的形式

其中，A_n和分别是第n个谐波的幅度和相位，并有和t是时间。将表达式(1)整理成线性矩阵形式Ax＝b，

其中，b＝[V₁,V₂,…,V_P]^T中的V_p(p＝1,2,…P)是t_p时刻的接收数据；x＝[α₁,…,α_N,β₁,…,β_N]^T是谐波系数，N是谐波个数，根据接收信号的带宽(1～3kHz)取N＝100。求解线性方程组(2)得到谐波系数矩阵x。

步骤3c:分别计算频率搜索范围内不同f₀的模型估计值V_harmonic(f₀)；

步骤3d:计算V_harmonic(f₀)与测量数据b的差值(误差)的2范数：||b-V_harmonic(f₀)||₂；

步骤3e:确定使误差2范数取得最小值时对应的频点f_min，从而获得两侧相邻的频点(f_m1,f_m2)；

步骤3f:确定下一次扫描范围(f_m1,f_m2)，f₀选取步长为(f_m2-f_m1)/M，其中M＝3～7；

步骤g:判断是否扫描5次(可以设置为3-7次，每次搜索的工频基频频率相对变化量小于0.1％即可)，如果达到，则根据误差2范数最小时对应的频点f_min进行工频谐波建模；否则重复步骤3c～3g；

步骤h:从步骤(2)的结果中减去工频谐波模型V_harmonic(f_min)，从而去除工频谐波噪声；

如图2所示的一组接收信号数据，图2(a)是其时域信号，MRS信号淹没在噪声之中；图2(b)是其功率谱，在谐波的频点上均存在不同幅度的谐波噪声。谐波建模后，在测量数据b中减去V_harmonic即完成了消噪过程，如图2(c)和2(d)所示。谐波建模消噪后，MRS信号的轮廓呈现衰减趋势，在功率谱中只存在MRS信号的频点。

本发明提供的自相关和叠加处理的具体步骤为：

步骤4a：MRS信号用V_MRS表示：

其中，e₀表示初始振幅；表示平均弛豫时间，f表示拉莫尔频率(与当地地磁场有关)，表示初始相位。

根据自相关原理，将V_MRS代入自相关函数的表达式中：

其中，τ为时间间隔。由于自相关函数是关于自变量τ的实偶函数，所以τ只取其正半轴。表达式(4)中的积分上下限分别取τ和t_max，t_max是测量时间长度，通过推导得到V_MRS的自相关表达式为：

R(τ)＝RC1(τ)+RC2(τ)+RC3(τ)+RC4(τ)

其中，

由表达式(5)～(7)可知，MRS信号自相关函数包含四个分量RC1～RC4，MRS信号自相关后的前半部分可以只用RC1代替；后半部分可以用RC1+RC2代替；RC3+RC4的结果很小，可以忽略。

步骤4b：判断采集时间是否大于1s，如果大于1s截取自相关结果的前半段用RC1近似；如果小于1s自相关结果用RC1+RC2近似。

步骤4c：基于步骤4b的判断后将多次测量数据b的自相关结果R进行叠加，即其中N是测量次数。N次叠加后信噪比提高倍，进一步减小随机噪声。

图3给出了含噪MRS信号经自相关和叠加处理的结果。图3(a)中黑色线是理想的MRS信号，参数为e₀＝100nV，f＝2326Hz和灰色线是加入了100nV均方差的高斯白噪声后的测量信号，此时数据的信噪比为0dB，MRS信号基本淹没在噪声中。图3(b)是测量信号和MRS信号的频谱，可见MRS信号的频谱峰值是噪声频谱均值的10倍。经过自相关函数后，测量信号(灰色线)和MRS信号(黑色线)的结果见图3(c)所示。此时，灰色线与黑色线类似，也呈现明显的衰减趋势，只出现了少量的随机噪声，信噪比为23.75dB，与自相关之前相比大幅提高。从图3(d)中的频谱分析也可以得出，MRS信号自相关后的频谱峰值约为噪声频谱均值的100倍，与自相关之前提高了一个数量级。因此，自相关函数能够有效压制随机噪声，提高信噪比。

如图4所示，本发明提供的参数提取方法的具体步骤为：

步骤5a:对叠加的数据进行希尔伯特变换，希尔伯特变换表达式为：可将实信号R(t)转化为复信号其包络信号为再通过低通滤波，得到两个正交分量V_RC＝V_real+i·V_imag；

步骤5b:利用非线性拟合方法对步骤a得到的两个正交分量进行参数提取如式(8)，

获得MRS信号的特征参数A₀和其中再计算e₀结果，参数提取结果为初始振幅e₀和弛豫时间

图3中的数据经过Hilbert变换后的两个分量见图4所示，其中图4(a)是未经过自相关函数的曲线，图4(b)是经过自相关函数的曲线。分别利用非线性拟合方法求解得到的两个分量的拟合曲线见图4中黑色实线和黑色虚线所示。未经过自相关函数的数据信噪比低，参数提取结果分别为e₀＝108.3nV，f＝2326.1Hz；而经过自相关函数的数据信噪比提高明显，参数提取结果为e₀＝100.5nV，f＝2326.0Hz，误差分别降低了7.8％、10.0％和0.01％。因此，自相关函数能够提高数据的信噪比，从而提高了参数提取的准确度。

实施例1

基于工频谐波建模和自相关的地面核磁共振信号提取方法，本发明模拟了一组包含16次独立采集的测量数据集，参照图1，包括以下步骤：

步骤(1)：利用式(3)构造含噪MRS信号如图5(a)所示，其参数为e₀＝100nV，f＝2326Hz和每次采集的谐波噪声基频在49.9～50.1Hz随机产生，谐波个数为100，谐波幅度在200nV内随机分布；每次采集的随机噪声为200nV的高斯白噪声；采集数据中未加入尖峰噪声，假设所有的尖峰噪声已通过已有的方法消除。

步骤(2)：利用谐波建模的方法得到工频谐波模型其中f₀＝50.018125Hz，从含噪MRS信号中减去工频谐波模型，得到去工频信号，如图5(b)所示，可见工频谐波消噪后，数据的信噪比大幅提升，平均提高了17.03dB，数据幅度减小到±500nV范围内；

步骤(3):对步骤(2)的结果进行自相关和叠加处理，减小随机噪声，得到消噪后信号，如图5(c)所示，可见自相关处理后，数据呈现明显的衰减趋势，除了信号初期存在脉冲(随机噪声的自相关特性)，其他时刻的数据幅度比自相关之前的数据明显提高而随机噪声的影响较小，因此信噪比有较大程度的提升，为15.64～16.43dB，信噪比平均提高了16.10dB；

步骤(4):对步骤(3)的结果进行MRS信号参数提取，利用希尔伯特变换(参考频率2325Hz)和低通滤波(截至频率200Hz)，得到MRS信号的实部和虚部见图5(d)和5(e)所示。建模消噪处理后，由于还存在大量的随机噪声，数据的实部和虚部曲线区分不明显，只能观察到曲线的整体趋势稍有起伏(1Hz的震荡)。相比自相关处理后，数据的实部和虚部能够清晰地区分开，两者都呈现出1Hz的震荡衰减过程，再对16次采集和处理的数据进行叠加，得到MRS信号实部和虚部的叠加曲线见图5(f)和(g)所示。叠加处理能够进一步降低随机噪声，16次叠加信噪比提高4倍(12.04dB)，利用非线性拟合方法对两组数据进行MRS信号参数提取，通过建模消噪和自相关处理后的参数提取结果为：e₀＝100.45nV，f＝2326.0Hz，相对误差分别为：0.45％，-2.14％和0，均控制在±3％以内，满足应用要求。

实施例2

本实施例以长春市农安县烧锅镇的太平池水库旁(该地地磁场为54720nT，对应的Larmor频率在2330Hz，变化范围小于2Hz)利用自主研制的地面核磁共振地下水探测仪器进行了野外数据采集与数据处理实验。为了从浅到深地对含水层进行探测，在0.2～8.5As的范围内由小到大按对数分布设置了20组发射脉冲矩，每组脉冲矩重复叠加16次，采样率为50kHz，采集时间为1s。将第10、16和20个脉冲矩采集的MRS信号作为本发明方法的处理对象。如图1所示，基于工频谐波建模和自相关的地面核磁共振信号提取方法，包括以下步骤:

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

步骤(2):利用统计方法判断是否存在尖峰噪声，如果存在，去除尖峰噪声并用插值结果代替，如果不存在，则测量数据保持不变；

步骤(3)：利用谐波建模的方法得到模型其中f₀＝50.018125Hz，从步骤(2)的结果中减去工频谐波模型，得到去工频信号，如图6(a-1)为q＝1.4As的MRS信号；6(a-2)为q＝4.2As的MRS信号；6(a-3)为q＝8.6As的MRS信号所示，可见实测数据中存在大量的谐波噪声，噪声水平在2500～3000nV范围内，而MRS信号幅度在80～375nV，因此信噪比较低。经过工频谐波消噪处理后，大多数谐波噪声显著降低，谐波基频在49.975～50.041Hz之间，MRS信号(2330Hz)的峰值比较明显，其中q＝1.4As的信号较大，而q＝4.2As和8.6As的信号较小，信噪比提高了近20dB；

步骤(4):对步骤(3)的结果进行自相关和叠加处理，减小随机噪声，得到消噪后信号；

步骤(5):对步骤(4)后的结果进行MRS信号参数提取，如图6(b-1)，6(b-2)，6(b-3))所示，3个脉冲矩下信噪比增量分别为12.05dB、2.98dB和7.71dB，此时经过非线性拟合后，参数提取结果为：

q＝4.2As时，e₀拟合结果为80.7nV，不准确度为15.9nV，而拟合结果为235.9ms，不准确度为38.2ms；q＝1.4As时，e₀拟合结果为396.8nV，不准确度为17.4nV，而拟合结果为248.6ms，不准确度为10.8ms；q＝8.6As时，e₀拟合结果为175.6nV，不准确度为13.4nV，而拟合结果为263.1ms，不准确度为33.1ms，说明SNR的增加，能够提高非线性拟合结果的准确度。

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