一种基于互相关的核磁共振全波信号噪声滤除方法

摘要

本发明涉及一种基于互相关的核磁共振全波信号噪声滤除方法，是利用噪声与拉莫尔频率的正弦信号不相关，而FID幅度衰减正弦信号与拉莫尔频率的正弦信号具有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波形的包络，并重构不含噪声的互相关波形，最后利用解卷积算法提取核磁共振全波数据中的FID信号。该方法运算数据计算量小，可以同时压制工频及其谐波噪声、随机噪声和尖峰噪声，明显提高核磁共振全波数据信噪比，有利于扩展核磁共振找水仪的应用范围和探测深度，而且该发明不需要额外的参考消噪装置，节约了成本，一次采集的全波数据通过互相关即可达到理想的滤波效果，节省测量时间。

说明书

技术领域

本发明涉及一种核磁共振测深(Magnetic Resonance Sounding,MRS) 信号噪声滤除方法，具体是基于互相关噪声研制原理的核磁共振全波信号噪声的 滤除方法。

背景技术

核磁共振地下水探测利用人工激发的电磁场使地下水中的氢核形成宏观磁 矩，这一宏观磁矩在地磁场中产生旋进运动，使用线圈接收宏观磁矩进动产生的 电磁信号，该电磁信号的幅度随时间按指数形式衰减，成为自由感应衰减(Free  Induction Decay，FID)信号，通过对FID进行分析可以探测地下是否存在地 下水。核磁共振找水是一种快速、直接的地下水探测方法。FID信号的表达式为：

上式中，ω₀表示拉莫尔角频率，由当地的地磁场强度决定，一般为1000- 3000Hz；FID信号的三个特征参数E₀(初始振幅)、T₂^*(平均弛豫时间)和(初始相位)分别与含水层的分布、厚度、平均含水率、渗透率和导电性等信息 有着密切的关系，准确测量FID信号是进行核磁共振地下水探测的关键，测量 FID信号时需要在野外勘探过程中将线圈接收信号的时域波形记录下来，该时域 波形成为核磁共振全波信号。进行地下水探测时，核磁共振信号非常微弱，一般 为nV级，而周围电磁场噪声却非常强烈，如工频及其谐波噪声、随机噪声和尖 峰噪声，即使在信号检测电路中使用了带通滤波器，核磁共振全波信号的信噪比 依然很低，一般低于0dB，这严重限制了核磁共振地下水探测的应用。

目前，核磁共振全波信号一般采用数据叠加方法，该方法能抑制大部分的随 机噪声，但对工频噪声和尖峰噪声抑制效果差，且需要多次采集核磁共振信号， 测量效率低。CN103823244A公开了一种磁共振三分量消噪装置，该发明基于 参考通道噪声与主通道含噪信号中噪声的相关性来实现信噪分离，但由于某些噪 声的无规律性和不稳定性以及混合机制的不确定性，限制了算法的应用，而且增 加参考线圈和参考通道使仪器系统庞大复杂。CN104459809A公开了一种基于 独立成分分析的全波核磁共振信号噪声滤除方法，首先利用核磁共振测深探水仪 采集FID信号，通过频谱分析获得采集信号中含有的工频谐波干扰频率，采用 数字正交法构造输入通道信号解决欠定盲源分离问题；然后，将构造的输入通道 信号与采集的全波信号一并作为输入信号进行独立成分分析，得到分离的FID 信号；最后采用频谱校正法解决分离后MRS信号幅值不定问题，进而提取去噪 后的FID信号。该发明主要针对核磁共振全波信号中工频谐波干扰或某一单频 干扰，不能滤除随机噪声和尖峰噪声，该方法算法过程复杂，计算量大，而且输 出结果不稳定。

发明内容

本发明的目的就在于针对上述现有技术的不足，提供一种基于互相关的核磁 共振全波信号噪声滤除方法。

本发明的目的是通过以下方式实现的：

一种基于互相关的核磁共振全波信号噪声滤除方法，包括以下步骤：

A、读取核磁共振接收机野外采集的全波数据s(k)及其采样率fs；

B、使用快速傅里叶变换求取全波数据的频谱S(ω)；

C、读取频谱S(ω)中FID信号的频率f0所对应的相位θ0；

D、读取频谱S(ω)中，拉莫尔频率f0左右两侧的两个工频谐波频率f1和f2， 及其所对应的幅度a1，a2，相位θ1和θ2；

E、重构两个工频谐波噪声的波形n1(k)和n2(k)；

F、在全波数据s(k)中减去n1(k)和n2(k)，得到x(k)；

G、将x(k)与参考信号cos(2πkf0/fs)进行互相关运算，得到R(m)；

H、利用希尔伯特变换获取R(m)的包络信号，En(m)；

I、将包络信号En(m)分为长度相等的两段，并分别使用双指数函数拟合包 络曲线；

J、使用双指数函数的拟合结果重构包络信号，Env(m)＝a·exp(bmf0/fs)+ c·exp(dmf0/fs)，a、b、c、d为拟合参数；

K、构造不含噪声的互相关波形rf(m)＝Env(m)·cos(2πmf0/fs+θ0)；

使用参考信号cos(2πkf0/fs)对rf(m)进行解卷积运算，得到的解卷积结果ff(k) 即为滤除噪声的核磁共振全波信号。

有益效果：本发明提供的基于互相关核磁共振全波信号噪声滤除方法，是利 用噪声与拉莫尔频率的正弦信号不相关，而FID幅度衰减正弦信号与拉莫尔频 率的正弦信号具有相关性的特点，通过互相关运算滤除噪声，然后拟合互相关波 形的包络，并重构不含噪声的互相关波形，最后利用解卷积算法提取核磁共振全 波数据中的FID信号。该方法运算数据计算量小，可以同时压制工频及其谐波 噪声、随机噪声和尖峰噪声，明显提高核磁共振全波数据信噪比，有利于扩展核 磁共振找水仪的应用范围和探测深度，而且该发明不需要额外的参考消噪装置， 节约了成本，一次采集的全波数据通过互相关即可达到理想的滤波效果，节省测 量时间。

附图说明

附图1是基于互相关的核磁共振全波信号噪声滤除方法的流程图。

附图2是野外采集的不含噪声的核磁共振全波数据的波形图。

附图3是野外采集的含有噪声的核磁共振全波数据的波形图。

附图4是含有噪声的核磁共振全波数据的频谱图。

附图5是核磁共振全波数据与参考信号进行互相关后的波形图。

附图6是互相关波形的包络拟合曲线。

附图7是重构后得到的不含噪声的互相关波形图。

附图8是解卷积得到的滤除噪声的核磁共振全波数据波形图。

具体实施方式

下面结合附图和实施例作进一步详细说明：

如图1所示，一种基于互相关的核磁共振全波信号噪声滤除方法，包括以下 步骤：

A、读取核磁共振接收机野外采集的全波数据，核磁共振找水仪接收机野外 工作时以恒定的采样率f_s采集数据，读取接收机记录的全波数据s(k)，k＝1，2， 3…N；

B、使用快速傅里叶变换求取全波数据的频谱S(ω)，在Matlab中使用fft(s,N) 命令可得到全波的频谱S(ω)，S(ω)包含N个频率所对应的复数，频率间隔为f_s/N；

C、读取频谱S(ω)中FID信号的频率f₀所对应的相位θ₀，S(ω)包含N个频 率所对应的复数代表了该频率的幅度和相位信息，即

${\theta }_0=\arctan \left[\frac{\mathrm{im}\left(S(N{f}_0/f_s)\right)}{\mathrm{re}\left(S(N{f}_0/f_s)\right)}\right]$

上式中，abs表示取模运算符，arctan为反正切运算符，im为取虚部运算符， re为取实部运算符；

D、读取频谱S(ω)中，拉莫尔频率f₀左右两侧的两个工频谐波频率f₁和f₂， 全波数据与参考信号作互相关运算时，如果在频谱有上接近拉莫尔频率的周期信 号，该周期信号会影响互相关运算的结果，这是由于当周期噪声的频率接近信号 频率时，周期噪声与FID信号的相关性增强，互相关噪声抑制效果变差。因此， 需要将频率最接近拉莫尔频率的两个工频谐波滤除掉。读取这两个工频谐波噪声 的幅度a₁，a₂，相位θ₁和θ₂；

E、重构两个工频谐波噪声的波形n₁(k)和n₂(k)，由于这两个工频谐波为频率 已知的周期信号，当其幅度和相位已知时，其波形就可唯一确定，可使用如下公 式重构这两个工频谐波噪声的波形：

$n_1\left(k\right)=a_1\cos (\frac{2\mathrm{\pi k}{f}_1}{\left(f\right)s}+{\theta }_1),n_2\left(k\right)=a_2\cos (\frac{2\mathrm{\pi k}{f}_2}{f_s}+{\theta }_2)$

F、在全波数据s(k)中减去n₁(k)和n₂(k)，得到x(k)

x(k)＝s(k)-n₁(k)-n₂(k)

得到的x(k)不含有频率最接近拉莫尔频率的两个工频谐波噪声，但依然含有 其他频率的工频噪声、随机噪声和尖峰噪声；

G、将x(k)与参考信号cos(2πkf₀/f_s)进行互相关运算，得到R(k)，计算公式为：

$R\left(m\right)=\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}\cos \left(\frac{2\pi f_{0}m}{f_s}\right)x(n-m),m=\mathrm{1,2,3},...,2N-1$

相关检测是一种噪声抑制的有效方法，它基于信号与噪声的统计特性进行检 测，互相关函数是两个时域信号相似性的一种度量，被检测信号与参考信号具有 相关性，而与工频及其谐波噪声、随机噪声和尖峰噪声的相关性较差或没有相关 性，利用这一差异可把被检测信号与噪声区分开了。虽然互相关检测一般应用于 周期的正弦信号，但是核磁共振的FID信号是幅度按指数规律衰减正弦信号，与 频率为拉莫尔频率的参考正弦信号cos(2πkf₀/f_s)具有较强相关性，而与随机噪声 和尖峰噪声没有相关性，与频率不等于拉莫尔频率的周期噪声相关性较差，因此， 互相关检测适用于核磁共振全波数据噪声滤除。两个长度均为N的数据进行互 相关运算后得到的互相关数组长度为2N-1；

H、利用希尔伯特变换获取R(m)的包络信号En(m)，互相关数据R(m)的希尔伯 特变换为：

$R^{\prime }\left(m\right)=\frac{1}{\pi }\underset{n=1}{\overset{2N-1}{\Sigma }}\frac{R(m-n)}{n}$

信号经希尔伯特变换后，在频域各频率分量的幅度保持不变，但相位将出现 90°相移。得到解解析信号为：z(m)＝R(m)+jR′(m)

解析信号z(k)的瞬时振幅即为互相关R(m)波形的包

络En(m)；

I、将En(m)分为两段，每段长度均为N，即En₁(k)＝En(1:N)，En₂(k)＝En (N:2N-1)，分别使用双指数函数拟合En₁(k)和En₂(k)，采用方法为最小二乘拟合， 即

$\min \left\{\underset{k=1}{\overset{N}{\Sigma }}{[{\mathrm{En}}_1\left(k\right)-a_1·\exp (b_{1}k{f}_0/f_s)+c_1·\exp (d_{1}k{f}_0/f_s)]}^2\right\}$

$\min \left\{\underset{k=1}{\overset{N}{\Sigma }}{[{\mathrm{En}}_2\left(k\right)-a_2·\exp (b_{2}k{f}_0/f_s)+c_2·\exp (d_{2}k{f}_0/f_s)]}^2\right\}$

a₁、b₁、c₁、d₁；a₁、b₁、c₁、d₁为拟合得到的参数；

J、使用双指数函数的拟合结果重构包络信号，Env(1:N)＝a₁·exp(b₁kf₀/f_s)+ c₁·exp(d₁kf₀/f_s)，Env(N:2N-1)＝a₂·exp(b₂kf₀/f_s)+c₂·exp(d₂kf₀/f_s)，拟合后的包络信号 将不含有噪声；

K、构造不含噪声的互相关波形rf(m)＝Env(m)·cos(2πmf₀/f_s+θ₀)，Env(m)为拟 合得到的包络，cos(2πmf₀/f_s+θ₀)为FID信号相位相同的正弦信号，其频率为拉莫 尔频率；

使用参考信号cos(2πkf₀/f_s)对rf(k)进行解卷积运算，互相关函数可表示为卷 积运算：

R₁₂＝f₁(t)*f₂(t)，

“*”表示卷积运算，可见将f₂(t)反褶(变量取符号)与f₁(t)做卷积即得到 f₁(t)与f₂的互相关函数，同时，由于参考信号cos(2πkf₀/f_s)的反褶为其本身，R(m) 可以看作x(k)与参考信号cos(2πkf₀/f_s)的卷积运算，重构的互相关波形rf(m)是不 含噪声的核磁全波数据与参考信号的卷积。因此，通过反卷积可以得到不含噪声 的核磁全波数据。在Matlab是使用解卷积命令deconv即可得到ff(k)： ff＝deconv(rf,cos(2πf₀t))，得到的解卷积结果ff(k)即为滤除噪声的核磁共振全波信 号。

野外采集的核磁共振全波数据含有FID的波形为：

$f\left(k\right)=1\times {10}^{-7}\exp (-\frac{k}{0.3\times 100000})\cos (\frac{2\pi \times 2326k}{100000}+\frac{\pi }{6})$

上式中，FID信号的初始幅度为100nV，平均弛豫时间T₂^*为300ms，拉莫 尔频率f0为2326Hz，初始相位为30°，接收机采样率为100kHz，信号采集时 间为500ms，采样点数为50000。不含噪声的FID信号波形见图2所示。

采集到的核磁共振全波数据中除了微弱的FID信号外，还包含工频及其谐波 噪声、随机噪声和尖峰噪声，噪声的表达式为：

n(k)＝10^-6square(t,50)+wag(1,N,-140)+10^-6pulstran(t,0.05:0.05:0.5,'gauspuls')

上式中，Matlab中square用于产生50Hz方波，来模拟工频工频及其谐波噪 声，方波的幅度为1μV；wgn(1,N,-140)用于产生标准差为100nV的零均值随机 噪声；pulstran函数用于产生尖峰噪声，尖峰噪声的幅度为1μV；这些噪声为加 性噪声叠加在FID信号上，核磁共振全波数据s(k)的表达式为：

s(k)＝f(k)+n(k)

该核磁共振全波数据的信噪比为-30dB。包含噪声的核磁共振全波数据波形 如图3所示。

在Matlab中使用函数fft函数对全波数据x(k)进行快速傅里叶变换，得到的 频谱图如图4所示，频谱数据S(ω)包含50000个频率所对应的复数，频率间隔为 2Hz。

从频谱数组S(ω)中找到第1164个元素为，对此复数取相角得到θ₀为31°。

根据拉莫尔频率2326Hz，得到最接近拉莫尔频率的两个工频谐波频率分别 为2250Hz和2350Hz，从频谱数组S(ω)中找到第1126个元素为1.3216e-09- 1.3472e-08i，对此复数取模得到a₁为27.73nV，θ₁为-1.5°；从频谱数组S(ω)中 找到第1176个元素为6.5767e-10-1.3849e-08i，对此复数取模得到a₂为27.07nV， θ₂为-1.4°。重构得到的这两个工频谐波波形为：

$n_1\left(k\right)=27.73\times {10}^9\cos (\frac{2\mathrm{\pi k}\times 2250}{100000}-\frac{1.5\pi }{180}),n_2\left(k\right)=27.07\times {10}^{-9}\cos (\frac{2\mathrm{\pi k}\times 2350}{100000}-\frac{1.4\pi }{180})$

从核磁全波数据s(k)中减去这两个工频谐波噪声得到x(k)，计算式如下：

x(k)＝s(k)-n₁(k)-n₂(k)

由于减去的这两个谐波噪声相对总体噪声而言很小，核磁共振全波数据的波 形几乎不会发生明显变化。将x(k)与参考信号cos(2πkf₀/f_s)按照如下公式进行互 相关运算，在Matlab中使用xcorr命令即可得到两个数组的互相关函数，得到的 R(m)波形如图5所示，可见互相关运算后的波形与原始的FID信号存在差异， 波形可分为两段：前一段呈现幅度按指数规律增加的正弦信号，后一段呈现幅度 按指数规律衰减的正弦信号。

利用希尔伯特变换获取R(m)的包络信号En(m)，将En(m)分为两段，每段 长度均为50000，即En₁(k)＝En(1:50000)，En₂(k)＝En(50000:99999)，分别使用 双指数函数拟合En₁(k)和En₂(k)，采用方法为最小二乘拟合，得到的拟合结果为：

a₁＝0.001364，b₁＝0.1022，c₁＝-0.001367，d₁＝-3.625；a₂＝-0.0001759， b₂＝0.4319，c₂＝0.001402，d₂＝-3.628。

使用双指数函数的拟合结果重构包络信号，Env(1:N)＝a₁·exp(b₁kf₀/f_s)+ c₁·exp(d₁kf₀/f_s)，Env(N:2N-1)＝a₂·exp(b₂kf₀/f_s)+c₂·exp(d₂kf₀/f_s)，拟合后的包络信号 将不含有噪声，图6是互相关波形的包络拟合曲线，包络拟合曲线分为两段。

构造不含噪声的互相关波形rf(m)＝Env(m)·cos(2πmf₀/f_s+θ₀)，Env(m)为拟合得 到的包络，图7是重构后得到的不含噪声的互相关波形图，重构后的互相关波形 比原互相关波形光滑。

使用参考信号cos(2πkf₀/f_s)对rf(k)进行解卷积运算，在Matlab是使用解卷积 命令deconv即可得到ff(k)：ff＝deconv(rf,cos(2πf₀t))，得到的解卷积结果ff(k)即 为滤除噪声的核磁共振全波信号，图8是解卷积得到的滤除噪声的核磁共振全波 数据波形图，对比图2和图8发现二者波形一致，说明基于互相关的核磁共振全 波信号噪声滤除方法能够恢复淹没在噪声中的FID信号。