CN109765629B - 一种能够压制同频噪声干扰的地面磁共振信号提取方法 - Google Patents

Abstract

本发明涉及一种能够压制同频噪声干扰的地面磁共振信号提取方法，包括：采用传统分步消噪流程对一组核磁共振响应数据的原始数据进行处理；用一组与地面磁共振信号拉莫尔频率相同的频率的正交参考信号，与消噪后的磁共振信号相乘获得同相通道结果和正交通道结果；将处理后的包含和频成分和差频成分的信号利用低通滤波滤除高频成分；将滤波后的差频成分利用求导等数学变换，得到线性函数；利用最小二乘法对线性函数进行直线拟合，获得直线斜率和纵坐标截距，进一步求取得到获得平均横向弛豫时间、初始振幅和初始相位的数值。解决了传统消噪方法无法处理同频噪声干扰的问题，有效抑制了同频噪声而且对于磁共振信号有很好的保护。

Description

一种能够压制同频噪声干扰的地面磁共振信号提取方法

技术领域

本发明属于涉及地面磁共振(Nuclear Magnetic Resonance，NMR)信号噪声消减领域，具体涉及一种能够压制同频噪声干扰的地面磁共振信号提取方法。

背景技术

地面磁共振探测技术可以对地下水进行直接和定量探测。地下水中氢原子被拉莫尔频率电磁场激发形成宏观磁矩，通过铺设线圈拾取宏观磁矩进动产生的核磁共振地下水探测信号。信号的初始振幅、平均横向弛豫时间、初始相位参数可以反映地下含水量和岩石孔隙度等信息。但核磁共振信号容易受到周围噪声的干扰，主要有工频谐波噪声、尖峰噪声、随机噪声等，此外当拉莫尔频率与工频谐波频率一致或接近时，还有同频噪声的干扰。这些噪声的存在会影响到后期反演解释结果的准确性，因此研究磁共振地下水探测信号的消噪方法十分必要。

针对核磁共振地下水探测数据的噪声，已经有很多方法消除噪声。例如针对尖峰噪声，基于能量运算的方法，采用中位数的绝对偏差法确定阈值来检测尖峰噪声进而剔除尖峰噪声的，而对磁共振探测数据采用统计叠加方法也可以对随机噪声进行抑制。针对工频谐波噪声有自适应抵消技术和基于快速独立分量分析等方法，这些方法虽然能消除核磁共振地下水探测数据中的工频谐波噪声，但是对同频噪声却不能有效剔除干扰。

CN104614778A公开了一种“基于ICA的核磁共振地下水探测信号噪声消除方法”，通过录入三组核磁共振响应数据，确认核磁共振拉莫尔频率附近的工频谐波噪声频率，经过构造与工频谐波干扰频率相同的正弦函数、余弦函数，并将其与核磁共振响应数据一起构成观测信号，采用快速独立分量分析算法对信号进行分离，并保留解混信号中噪声成分少而核磁共振信号理想的成分，并进行重构信号来消除工频谐波噪声。但此方法仅仅针对核磁共振拉莫尔频率附近的工频谐波噪声有效，当拉莫尔频率与工频谐波频率一致或接近时，即存在同频噪声的干扰时，该方法不能保证对同频谐波噪声的识别，这些噪声残余会影响到核磁共振地下水探测信号参数的准确提取，影响到对地下含水量等参数的反演解释结果。

CN108254794A公开的“一种基于建模反恢复技术的磁共振消噪方法”，该方法利用极窄的低通滤波压制噪声，通过后期拉普拉斯和求导变换恢复原始信号。该方法可以有效压制远离信号频率的噪声，同样，当噪声与磁共振信号具有相同频率时，该方法时无能为力的。

发明内容

本发明所要解决的技术问题在于提供一种能够压制同频噪声干扰的地面磁共振信号提取方法，解决了传统消噪方法无法处理同频噪声干扰的问题。

本发明是这样实现的，

一种能够压制同频噪声干扰的地面磁共振信号提取方法，该方法包括：

步骤A，采用传统分步消噪流程对一组核磁共振响应数据的原始数据进行处理；

步骤B，用一组与地面磁共振信号拉莫尔频率相同的频率的正交参考信号，与步骤A消噪后的磁共振信号相乘获得同相通道结果和正交通道结果，两个通道结果都包含和频成分和差频成分的信号；

步骤C，将步骤B处理后的包含和频成分和差频成分的信号利用低通滤波滤除高频成分，保留含有核地面磁共振信号的差频成分；

步骤D，将滤波后的差频成分利用求导等数学变换，得到线性函数；利用最小二乘法对线性函数进行直线拟合，获得直线斜率和纵坐标截距，根据平均横向弛豫时间、初始振幅和初始相位与直线斜率和纵坐标截距关系求取得到获得平均横向弛豫时间、初始振幅和初始相位的数值。

进一步地，所述步骤A利用传统分步消噪流程对原始数据进行处理，包括基于能量算子检测包含尖峰噪声的数据，并剔除该噪声，利用建模法压制工频谐波噪声，再利用叠加平均压制随机噪声。

进一步地，步骤B具体包括：

B1、设置正交参考信号为参考信号r₁(k)和参考信号r₂(k)，参考信号的频率与拉莫尔频率f_L一致；

B2、将核磁共振响应数据与正交的两个参考信号分别相乘，利用三角和差公式获得磁共振信号、同频噪声与参考信号的和频成分和差频成分；

其中和频成分为正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相加的频率成分，差频成分为正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相减的频率成分。

进一步地，步骤C具体包括：

C1、设置低通滤波的截止频率为ω_c，使通过低通滤波后仅通过低频分量，远离拉莫尔频率的噪声与和频成分被滤除，只剩下差频成分；

C2、经过低通滤波后的含有核地面磁共振信号的差频成分经拉普拉斯变换；

C3、经过步骤2的拉普拉斯变换后利用拉普拉斯反变换，得到滤波后的信号u₁(k)和信号u₂(k)，均包含有信号项与同频噪声项。

进一步地，将滤波后的差频成分利用求导等数学变换，得到线性函数包括：

D11、将滤波后获得的差频成分与拉莫尔频率与时间相乘的e指数相乘，得到信号g(k)；

D12、对信号g(k)进行求导得到信号p(k)；

D13、对信号p(k)进行求导得到信号h(k)；

D14、对信号h(k)求导，得到信号q(k)；

D15、对信号q(k)，取底数为e的对数，得到信号f(k)，该信号是关于k的线性函数。

进一步地，获得直线斜率-a，利用纵坐标截距建立方程组，方程组如下：

f₁(0)＝ln(aE_r1w_c)；

f₂(0)＝ln(aE_r1w_c)；

得到E_r1和E_r2，平均横向弛豫时间、初始振幅和初始相位与直线斜率和纵坐标截距关系如下：

其中，E₀为初始振幅、T₂*为平均横向弛豫时间、

为初始相位、E_r1和E_r2分别为同相通道结果和正交通道结果中包含磁共振信号的差频分量。

本发明与现有技术相比，有益效果在于：本发明在传统地面磁共振信号消噪方法的基础上，利用相敏检波技术对包含同频噪声干扰的数据进行处理，通过数学变换，有效提取出了地面磁共振信号，解决了传统消噪方法无法处理同频噪声干扰的问题。

通过对比消噪后核磁共振信号与理论信号初始振幅及消噪后核磁共振信号与理论信号频谱图，该方法有效抑制了同频噪声而且对于磁共振信号有很好的保护。

附图说明

图1为本发明提供的磁共振地下水探测信号消噪方法流程图；

图2为本发明实施例提供的相敏检波技术原理流程图；

具体实施方式

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

图1为一种能够压制同频噪声干扰的地面磁共振信号提取方法的工作流程图。针对在拉莫尔频率处同频噪声的干扰，通过正交参考信号提取在拉莫尔频率处的地面磁共振信号和同频噪声，并通过求导等数学变换及最小二乘法的直线拟合方法去除同频噪声，提取核磁共振地下水探测信号的参数。

图2是相敏检波技术原理流程图。采用正交参考信号来去除地面磁共振数据中的工频谐波噪声和随机噪声。

参见图1所示，一种能够压制同频噪声干扰的地面磁共振信号提取方法，包括以下步骤：

A、录入一组核磁共振响应数据，核磁共振响应数据包含地面磁共振信号和噪声，其中噪声包括随机噪声、工频谐波噪声、同频噪声和尖峰噪声等。利用传统分步消噪流程对原始数据进行处理，包括基于能量算子检测包含尖峰噪声的数据，并剔除该噪声，利用建模法压制工频谐波噪声，再利用叠加平均压制随机噪声；处理后的数据主要包括地面磁共振信号和同频干扰信号，该组核磁共振响应数据，记为x(k)，采样频率为f_s,

x(k)＝s(k)+n_c(k),

其中，

为磁共振信号，E₀为初始振幅，T₂*为平均横向弛豫时间，拉莫尔频率为f_L，

为初始相位，为了简化公式，令

n_c(k)为同频噪声，采样频率为f_s，

B、采用与地面磁共振信号拉莫尔频率相同的频率f_L的正交参考信号，与磁共振信号相乘获得和频成分和差频成分，具体包含以下步骤：

B1设置正交参考信号r₁(k)和r₂(k)，参考信号的频率与拉莫尔频率f_L一致，

其中，V_r是两个参考信号的幅值，一般可以设置为V_r＝1,

B2将核磁共振响应数据分别与参考信号分别相乘，利用三角和差公式获得同相通道结果和正交通道结果，两个通道结果都包括地面磁共振信号、同频噪声与参考信号的和频成分和差频成分：

同相通道结果为：

正交通道结果为：

其中和频成分是正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相加的频率成分，差频成分是正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相减的频率成分；

C、利用低通滤波滤除高频成分，保留含有地面磁共振信号的低频分量，具体步骤如下：

C1、设置低通滤波的截止频率为ω_c，使通过低通滤波后仅通过低频分量，远离拉莫尔频率的噪声与和频成分被滤除，只剩下差频成分；

两个输出通道结果经过低通滤波后的拉普拉斯变换为：

其中，一阶RC低通滤波的传递函数为：

其中，E_r1和N_r1分别为同相通道结果v₁(k)中包含磁共振信号的差频分量和包含同频噪声的差频分量，其中，

同样，E_r2和N_r2分别为正交通道结果v₂(k)中包含磁共振信号的差频分量和包含同频噪声的差频分量，其中，

C2、对C1步骤中的U₁(s)和U₂(s)进行拉普拉斯反变换，得到滤波后的时间域信号u₁(k)和u₂(k)，其中都包含有信号项与同频噪声项：

D、滤波后的时间域信号u₁(k)和信号u₂(k)利用求导等数学变换和最小二乘法进行直线拟合，提取核磁共振地下水探测信号初始振幅、平均横向弛豫时间、初始相位参数。

求导等数学变换包括以下步骤：

D11、将步骤C滤波后的时间域信号u₁(k)和信号u₂(k)与拉莫尔频率与时间相乘的e指数相乘，得到信号g₁(k)和g₂(k)，

D12、信号g₁(k)和信号g₂(k)进行求导得到信号p₁(k)和信号p₂(k)，

p₁(k)＝f_s[g₁(k+1)-g₁(k)]＝E_r1ω_cexp(ω_ck-ak)+N_r1ω_cexp(ω_ck)，

p₂(k)＝f_s[g₂(k+1)-g₂(k)]＝E_r2ω_cexp(ω_ck-ak)+N_r2ω_cexp(ω_ck)；

D13、对求导信号p₁(k)和信号p₂(k)与拉莫尔频率与时间相乘的e指数的倒数相乘，得到信号h₁(k)和信号h₂(k)，

h₁(k)＝p₁(t)exp(-ω_ck)＝E_r1ω_cexp(-ak)+N_r1ω_c，

h₂(k)＝p₂(t)exp(-ω_ck)＝E_r2ω_cexp(-ak)+N_r2ω_c；

D14对信号h₁(k)和h₂(k)求导，得到信号q₁(k)和信号q₂(k)，

q₁(k)＝f_s[h₁(k+1)-h₁(k)]＝-E_r1ω_caexp(-ak)，

q₂(k)＝f_s[h₂(k+1)-h₂(k)]＝-E_r2ω_caexp(-ak)；

D15对信号q₁(k)和信号q₂(k)，取底数为e的对数，得到信号f₁(k)和信号f₂(k)，该信号是关于k的线性函数,

f₁(k)＝ln|q₁(k)|＝ln(aE_r1ω_c)-ak，

f₂(k)＝ln|q₂(k)|＝ln(aE_r2ω_c)-ak；

对上述两个线性函数运用最小二乘法进行直线拟合，先获得直线斜率-a，再利用纵坐标截距建立方程组，方程组如下：

f₁(0)＝ln(aE_r1w_c)；

f₂(0)＝ln(aE_r1w_c)；

由于a和w_c已知，可以获得E_r1和E_r2，进一步地，可以获取平均横向弛豫时间、初始振幅和初始相位关于直线斜率与纵坐标截距关系，如下：

通过以上数学关系计算地面磁共振信号的参数，滤除同频噪声。

应用示例:

以理论磁共振信号加入同频噪声及随机噪声为例进行消噪处理：磁共振信号

其拉莫尔频率为f_L＝2250Hz，初始振幅为E₀＝200nV，平均横向弛豫时间T₂*＝150ms，初始相位

同频噪声为n_c(k)＝200cos(2π×2250×k+3π/12)，白噪声与地面磁共振信号功率相同。采样频率为f_s＝25KHz。在matlab环境下，基于正相敏检波技术及最小二乘数据拟合原理对核磁共振地下水探测数据中的同频噪声进行处理。

具体步骤如下：

A、录入一组核磁共振响应数据，利用传统分步消噪流程对原始数据进行处理，包括基于能量算子检测包含尖峰噪声的数据，并剔除该噪声，利用建模法压制工频谐波噪声，再利用叠加平均压制随机噪声；该组核磁共振响应数据记为x(k)，采样频率为f_s＝25KHz,其中，s(k)为磁共振响应信号，拉莫尔频率为f_L，n_c(k)为同频噪声，而n_r(k)为随机噪声,

x(k)＝s(k)+n_c(k)+n_r(k),

其中n_c(k)为同频噪声，采样频率为f_s，

为了简化公式令

B、正交参考信号采用与核磁共振地下水探测信号拉莫尔频率相同的频率f_L，与磁共振信号相乘获得和频成分和差频成分，具体包含以下步骤：

a、设置正交参考信号r₁(k)和r₂(k)，参考信号的频率与拉莫尔频率f_L一致，

b、将核磁共振响应数据分别与参考信号分别相乘，利用三角和差公式获得磁共振信号、各种噪声与参考信号的和频成分和差频成分：

其中和频成分是正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相加的频率成分，差频成分是正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相减的频率成分；

C、利用低通滤波滤除高频成分，保留含有核磁共振地下水探测信号的低频

分量，具体步骤如下：

a、设置低通滤波的截止频率为ω_c，使通过低通滤波后仅通过低频分量，远离拉莫尔频率的噪声与和频成分被滤除，只剩下差频成分；

b、经过低通滤波后的拉普拉斯变换为：

其中，一阶RC低通滤波的传递函数为：

c、利用拉普拉斯反变换，得到滤波后的信号u₁(k)和u₂(k)，其中包含有信号项与同频噪声项：

D、利用求导等数学变换和最小二乘法进行直线拟合，提取核磁共振地下水探测信号初始振幅、平均横向弛豫时间、初始相位参数。

求导等数学变换包括以下步骤：

a、将滤波后获得的差频成分与拉莫尔频率与时间相乘的e指数相乘，得到信号g₁(k)和g₂(k)，

b、对信号g₁(k)和g₂(k)进行求导得到信号p₁(k)和p₂(k)，

p₁(k)＝f_s[g₁(k+1)-g₁(k)]＝E_r1ω_cexp(ω_ck-ak)+N_r1ω_cexp(ω_ck)，

p₂(k)＝f_s[g₂(k+1)-g₂(k)]＝E_r2ω_cexp(ω_ck-ak)+N_r2ω_cexp(ω_ck)；

c、对求导信号p₁(k)和p₂(k)与拉莫尔频率与时间相乘的e指数的倒数相乘，得到信号h₁(k)和h₂(k)，

h₁(k)＝p₁(t)exp(-ω_ck)＝E_r1ω_cexp(-ak)+N_r1ω_c，

h₂(k)＝p₂(t)exp(-ω_ck)＝E_r2ω_cexp(-ak)+N_r2ω_c；

d、对信号h₁(k)和h₂(k)求导，得到信号q₁(k)和q₂(k)，

q₁(k)＝f_s[h₁(k+1)-h₁(k)]＝-E_r1ω_caexp(-ak)，

q₂(k)＝f_s[h₂(k+1)-h₂(k)]＝-E_r2ω_caexp(-ak)；

e、对信号q₁(k)和q₂(k)，取底数为e的对数，得到信号f₁(k)和f₂(k)，该信号是关于k的线性函数,

f₁(k)＝ln|q₁(k)|＝ln(aE_r1ω_c)-ak,

f₂(k)＝ln|q₂(k)|＝ln(aE_r2ω_c)-ak；

运用最小二乘法进行直线拟合，获得直线斜率和纵坐标截距，获得E_r1、E_r2和a的数值。利用平均横向弛豫时间、初始振幅和初始相位关于直线斜率与纵坐标截距关系，计算如下：

经过该方法滤除同频噪声，获得地面磁共振信号平均横向弛豫时间、初始振幅和初始相位参数数值，其初始振幅和平均横向弛豫时间误差为0.78％和1.15％。

基于相敏检波技术及最小二乘法进行数据拟合的噪声消噪方法可以很好地剔除核地面磁共振数据中的噪声，该方法有效抑制了同频噪声而且对于磁共振信号有很好的保护。

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

Claims (5)

1.一种能够压制同频噪声干扰的地面磁共振信号提取方法，其特征在于，该方法包括：

步骤A，采用传统分步消噪流程对一组核磁共振响应数据的原始数据进行处理；

步骤B，用一组与地面磁共振信号拉莫尔频率相同的频率的正交参考信号，与步骤A消噪后的磁共振信号相乘获得同相通道结果和正交通道结果，两个通道结果都包含和频成分和差频成分的信号；

步骤C，将步骤B处理后的包含和频成分和差频成分的信号利用低通滤波滤除高频成分，保留含有地面核磁共振信号的差频成分；

步骤D，将滤波后的差频成分利用求导等数学变换，得到线性函数；利用最小二乘法对线性函数进行直线拟合，获得直线斜率和纵坐标截距，根据平均横向弛豫时间、初始振幅和初始相位与直线斜率和纵坐标截距关系求取得到平均横向弛豫时间、初始振幅和初始相位的数值；

将滤波后的差频成分利用求导等数学变换，得到线性函数包括：

D11、将滤波后获得的差频成分与拉莫尔频率与时间相乘的e指数相乘，得到信号g(k)；具体包括：将步骤C滤波后的时间域信号u₁(k)和信号u₂(k)与拉莫尔频率与时间相乘的e指数相乘，得到信号g₁(k)和g₂(k)，

D12、对信号g(k)进行求导得到信号p(k)；具体包括：信号g₁(k)和信号g₂(k)进行求导得到信号p₁(k)和信号p₂(k)，

p₁(k)＝f_s[g₁(k+1)-g₁(k)]＝E_r1ω_cexp(ω_ck-ak)+N_r1ω_cexp(ω_ck)，

p₂(k)＝f_s[g₂(k+1)-g₂(k)]＝E_r2ω_cexp(ω_ck-ak)+N_r2ω_cexp(ω_ck)；

D13、对信号p(k)进行求导得到信号h(k)；具体包括：对求导信号p₁(k)和信号p₂(k)与拉莫尔频率与时间相乘的e指数的倒数相乘，得到信号h₁(k)和信号h₂(k)，

h₁(k)＝p₁(t)exp(-ω_ck)＝E_r1ω_cexp(-ak)+N_r1ω_c，

h₂(k)＝p₂(t)exp(-ω_ck)＝E_r2ω_cexp(-ak)+N_r2ω_c；

D14、对信号h(k)求导，得到信号q(k)；具体包括：对信号h₁(k)和h₂(k)求导，得到信号q₁(k)和信号q₂(k)，

q₁(k)＝f_s[h₁(k+1)-h₁(k)]＝-E_r1ω_caexp(-ak)，

q₂(k)＝f_s[h₂(k+1)-h₂(k)]＝-E_r2ω_caexp(-ak)；

D15、对信号q(k)，取底数为e的对数，得到信号f(k)，该信号是关于k的线性函数，具体包括：对信号q₁(k)和信号q₂(k)，取底数为e的对数，得到信号f₁(k)和信号f₂(k)，该信号是关于k的线性函数,

f₁(k)＝ln|q₁(k)|＝ln(aE_r1ω_c)-ak，

f₂(k)＝ln|q₂(k)|＝ln(aE_r2ω_c)-ak；

其中E_r1为同相通道结果中包含磁共振信号的差频分量；N_r1为同相通道结果中包含同频噪声的差频分量；E_r2为正交通道结果中包含磁共振信号的差频分量；N_r2为正交通道结果中包含同频噪声的差频分量。

2.按照权利要求1所述的方法，其特征在于，所述步骤A利用传统分步消噪流程对原始数据进行处理，包括基于能量算子检测包含尖峰噪声的数据，并剔除该噪声，利用建模法压制工频谐波噪声，再利用叠加平均压制随机噪声。

3.按照权利要求1所述的方法，其特征在于，步骤B具体包括：

B1、设置正交参考信号为参考信号r₁(k)和参考信号r₂(k)，参考信号的频率与拉莫尔频率一致；

B2、将核磁共振响应数据与正交的两个参考信号分别相乘，利用三角和差公式获得磁共振信号、同频噪声与参考信号的和频成分和差频成分；

其中和频成分为正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相加的频率成分，差频成分为正交参考信号与核磁共振响应数据相乘后，利用三角和差公式变换得到正交参考信号频率与核磁共振响应数据频率相减的频率成分。

4.按照权利要求1所述的方法，其特征在于，步骤C具体包括：

C1、设置低通滤波的截止频率为ω_c，使通过低通滤波后仅通过低频分量，远离拉莫尔频率的噪声与和频成分被滤除，只剩下差频成分；

C2、经过低通滤波后的含有核地面磁共振信号的差频成分经拉普拉斯变换；

C3、经过步骤C2的拉普拉斯变换后利用拉普拉斯反变换，得到滤波后的信号u₁(k)和信号u₂(k)，均包含有信号项与同频噪声项。

5.按照权利要求1或4所述的方法，其特征在于，获得直线斜率-a，利用纵坐标截距建立方程组，方程组如下：

f₁(0)＝ln(aE_r1ω_c)；

f₂(0)＝ln(aE_r2ω_c)；

得到E_r1和E_r2，平均横向弛豫时间、初始振幅和初始相位与直线斜率和纵坐标截距关系如下：

其中，E₀为初始振幅、T₂*为平均横向弛豫时间、

为初始相位、E_r1和E_r2分别为同相通道结果和正交通道结果中包含磁共振信号的差频分量,V_r是参考信号的幅值。

Images