CN104777442A - 一种核磁共振测深fid信号噪声抑制方法 - Google Patents

Abstract

本发明涉及一种核磁共振测深FID信号的噪声抑制方法。包括以下步骤：对磁共振测深系统检测到的信号进行频谱分析，利用归一化正交检测技术将检测到的信号分解为同向分量X、正交分量Y，并进行硬件滤波处理得到低频的FID信号；对采集数据采用非线性能量算子算法分别剔除FID信号中X、Y分量的尖峰噪声；基于主成分分析方法对X、Y分量分别进行初步的信噪分离；基于经验模态分解方法对PCA处理后的X、Y分量进一步分解提取信号趋势项；对EMD处理后的X、Y分量分别叠加求平均后获得e指数曲线。本发明全摒弃了传统滤波手段容易对信号成分造成损失等问题，实现了对核磁共振测深FID信号中包含各种复杂噪声的有效抑制。

Description

一种核磁共振测深FID信号噪声抑制方法

技术领域

本发明涉及一种核磁共振领域中信号处理，特别涉及一种基于主成分分析与经验模态分解的核磁共振测深FID信号噪声抑制方法。

背景技术

核磁共振测深(Magnetic resonance sounding，MRS)作为一种直接探测地下水的地球物理方法是利用水中氢质子的弛豫特性实现对地下水的探测的，并且已在地下水勘查、水文地质调查以及滑坡稳定性评价、堤坝渗漏、煤矿老空水探测和隧道突水超前探测预警等方面均获得了应用。伴随着科学技术的进步，经过十多年的发展，目前已经商品化的仪器系统包括法国IRIS公司的NUMIS系统、NUMISplus系统、NUMISpoly系统，美国Vista Clara公司的GeoMRI系统以及中国吉林大学JLMRS-I型找水仪及其联用仪、JLMRS-III型灾害水源探测仪。

然而，上述仪器在野外工程实验中均暴露出极易受到环境电磁噪声影响的问题，尤其在中国，电网覆盖率达到96％以上，存在大量的电力线和生活用电，这将导致引入的强工频谐波等电磁干扰会将核磁共振测深信号淹没，使得核磁共振探水仪的探测准确度和灵敏度大大降低。

为了解决强电磁干扰下核磁共振测深信号有效检测难题，国内外专家学者开展了一系列研究。在国外，2002年Plata and Rubio利用采集数据的多次叠加来提高信噪比，但这种方法最致命的弱点在于过多的叠加次数将导致系统工作效率的降低；Legchenko and Valla采用同步检波技术将信号分解成两个正交的分量，使用低通滤波器滤除高频噪声的方法在NUMIS系统被采用；2003年Legchenko and Valla采用区块对消法，正弦对消法和陷波滤波器法进行电力线谐波干扰去除；2006年Strehl等提出采用小波变换手段对尖峰噪声进行处理。国内的研究起步较晚，2006年，中国地质大学在潘玉玲教授带领下研究人员采用引进的NUMIS系统进行找水实验，李振宇等曾发表文章探讨应用小波分析的方法提高核磁共振找水信号信噪比，但在信噪比较高时去噪效果不明显；2009年，王中兴根据地面核磁共振信号的特征和奇异噪声的来源与特点，提出差阈值替代叠加是对奇异噪声进行削弱与抑制的有效方法以及采用4N倍采样率进行核磁共振FID信号提取与噪声抑制，可以部分抵消随机噪声与工频谐波干扰；2009-2011年，蒋川东采用统计叠加对磁共振FID信号中尖峰噪声以及部分随机噪声进行滤除，采用自适应陷波滤除工频谐波干扰。

上述研究均是针对磁共振测深FID信号噪声去除展开，并从时域或频域信号和噪声的特性出发，进而实现信噪分离的，但是由于环境噪声的复杂性、不稳定性和无规律性，使得算法的应用效果受到一定限制，从而影响到结果的精确性。

发明内容

本发明所要解决的技术问题在于提供一种核磁共振测深FID信号噪声抑制方法，旨在提供一种全新的现代滤波手段抑制核磁共振测深FID信号的噪声，以解决强电磁干扰下核磁共振测深信号的有效提取难题。

本发明是这样实现的，一种核磁共振测深FID信号噪声抑制方法，包括以下步骤：

a、对磁共振测深系统检测到的信号进行频谱分析，利用归一化正交检测技术将检测到的信号分解为同向分量X、正交分量Y，并进行硬件滤波处理得到低频的FID信号；所得到的同向分量X和正交分量Y表示为：

E(t)＝x(t)+i·y(t)

其中，x(t)为同向分量X的表达式，y(t)为正交分量Y的表达式，E₀为初始振幅，为平均衰减时间，δω是接收角频率与拉莫尔角频率ω₀的差值，而ε_x(t)和ε_y(t)分别为同向分量X和正交分量Y上的噪声成分；

b、对采集数据采用非线性能量算子算法分别剔除FID信号中同向分量X、正交分量Y的尖峰噪声；

c、基于主成分分析方法对步骤b处理后的同向分量X、正交分量Y分别进行初步的信噪分离，将与有用信号具有强相关性的噪声分离出来并消除；

d、基于经验模态分解方法对步骤c处理后的同向分量X、正交分量Y进一步分解提取信号趋势项；

e、对步骤d处理后的同向分量X、正交分量Y分别叠加求平均后获得e指数曲线。

根据本发明的优选实施例，步骤b中将检测到的核磁信号的X分量x(t)、Y分量y(t)分别离散化成x(n)和y(n)，经过非线性能量算子算法运算后，得到：

选择绝对中位差MAD方法设定阈值门限，设一组数据为z₁,z₂,…z_n表示x(n)或y(n)，则其MAD定义为：

MAD＝median_i{|z_i-median_j{z_j}|}

其中z_j表示在同一时刻下的全部数据，z_j表示每一组数据，定义阈值的上下门限为：

$\left\{\begin{array}{c}\mathrm{low}=\mathrm{median}\×0.01;\\ \mathrm{high}=\mathrm{median}+b\×\mathrm{MAD}.\end{array}\right.$

其中b为可根据信噪比以及多次重复试验确定的常系数，经过阈值限定后，取其范围

在上述公式中范围内的信号被保留，而尖峰噪声则被剔除。

根据本发明的优选实施例，步骤c中，首先，对于一维的时间序列x(t)或y(t)，用z(t)表示，定义一个N×M的实对称矩阵

求取矩阵Q的协方差矩阵C_Q，并获得其特征值λ₀,λ₁,…,λ_N-1和特征向量A₀,A₁,…,A_N-1,由此构造正交矩阵

A＝[A₀,A₁,…,A_N-1]^T

使得

I＝AQ

转化后的矩阵I的各个列向量依次为矩阵Q的第一主分量，第二主分量，以此类推，各个列向量呈现方差递减的趋势，若方差累计贡献率a满足：

$a=\underset{m=1}{\overset{i}{\mathrm{\Σ}}}{\mathrm{\λ}}_m/\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_m\×100\%>95\%$

其中，n为所有特征值的个数，i为所选取的主成分的个数；

由选取的主成分个数i对矩阵I和A降维，保留矩阵I和A的前i行，分别得到矩阵I_i和A_i，重构信号矩阵为

$Q_i=A_i^T{I}_i$

由重构信号矩阵得到消掉具有强相关性的噪声成分的时间序列x(t)或y(t)。

根据本发明的优选实施例，步骤c中采用多次循环主成分分析方法，若满足方差累计贡献率a>95％且主成分个数为1，或者是主成分分析方法适用性检验KMO<0.7，则判定循环过程停止。

根据本发明的优选实施例，步骤d，经验模态分解方法包括：通过拟合出x(t)或y(t)的极大值包络线e_max(t)和极小值包络线e_min(t)，选用G.Rilling sift criterion判据来获取IMF分量以及分解余项，该判据的评价函数为：

$\mathrm{\&sigma;}\left(t\right)=\frac{|e_{\max }\left(t\right)+e_{\min }\left(t\right)|}{|e_{\max }\left(t\right)-e_{\min }\left(t\right)|}$

设定三个阈值常量θ₁、θ₂和α，若满足：

①σ(t)中小于θ₁的比例不低于1-α，即

#{t∈D|σ(t)<θ₁}/#{t∈D}≥1-α

其中，D为信号持续的总时间，#{A}表示在集合A中元素的个数；

②对于整个时间范围内的每个时刻t，满足：

σ(t)＜θ₂

则认为经验模态分解得到的IMF分量满足要求；

当分解余项r_n(t)是单调函数或常量时，经验模态分解过程停止，最后经过经验模态分解后得到的信号x(t)或y(t)为：

$x\left(t\right)=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{c}_j\left(t\right)+r_k\left(t\right);$ 或 $y\left(t\right)=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{c}_j\left(t\right)+r_k\left(t\right)$

其中，c_j(t)表示每一个IMF分量，r_k(t)表示分解余项，即信号的趋势项；分别得到n组相应的同向分量X和正交分量Y，表示为：

x_i(t)＝r_xki(t)

y_i(t)＝r_yki(t)；(i＝1，2…n)

根据本发明的优选实施例，将同向分量X和正交分量Y各自叠加求取平均值，得到：

$x\left(t\right)=\frac{x_1\left(t\right)+x_2\left(t\right)+...+x_n\left(t\right)}{n}$

$y\left(t\right)=\frac{y_1\left(t\right)+y_2\left(t\right)+...+y_n\left(t\right)}{n}$

最终求出去噪后的e(t)指数曲线：

e(t)＝x(t)+i·y(t)

本发明与现有技术相比，有益效果在于：本发明鉴于核磁共振找水仪在同一环境下工作时，所采集FID信号中包含的噪声具有较大相关性的特点，提出基于PCA方法，将一组相关变量通过线性变换转换成另一组互不相关的变量，实现核磁共振测深FID信号中噪声的有效去除；同时，针对PCA处理后仍然含有少量不相关的工频干扰和随机噪声的问题，提出进一步采用EMD方法进行处理；这两种方法的结合，完全摒弃了传统滤波手段容易对信号成分造成损失等问题，实现了对核磁共振测深FID信号中包含各种复杂噪声的有效抑制，为后续数据反演获取精确的水文地质信息提供了技术保障。

附图说明

图1为采用本发明核磁共振测深FID信号噪声抑制方法的流程图；

图2为理想的核磁信号X分量、Y分量分解图，(A)为X分量分解图，(B)为Y分量分解图；

图3为含噪声的核磁信号X分量、Y分量分解图，(A)为X分量分解图，(B)为Y分量分解图；

图4为NEO去尖峰噪声之后的X分量、Y分量，(A)为X分量，(B)为Y分量；

图5为PCA信噪分离后的X分量、Y分量，(A)为X分量，(B)为Y分量；

图6为X分量、Y分量的EMD分解信号趋势项，其中，(1)、(2)、(3)、(4)表示的是X分量的EMD分解信号趋势项，(5)、(6)、(7)和(8)表示的是Y分量的EMD分解信号趋势项；

图7为叠加求平均后的e(t)曲线图。

具体实施方式

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

基于主成分分析与经验模态分解的核磁共振测深FID信号的噪声抑制方法，参照图1的流程图，包括以下步骤：

a、首先对磁共振测深系统检测到的信号进行频谱分析，比较信号接收频率与电流发射频率的关系，确认所检测到的信号是地下水中的磁共振信号E(t)，然后利用归一化数字正交检测技术将该信号分解为同向分量X、正交分量Y，通过硬件滤波滤掉其高频分量，从而得到低频的FID信号，如图2、图3所示。所得到的同向分量X、正交分量Y可以表示为：

E(t)＝x(t)+i·y(t)

其中，E₀为初始振幅，为平均衰减时间，δω是接收角频率与拉莫尔角频率ω₀的差值，而ε_x(t)和ε_y(t)分别为同向分量X、正交分量Y上的噪声成分。

b、利用非线性能量算子(Non-linear energy operator，NEO)算法，以下简称NEO算法，检测FID信号中的尖峰噪声并剔除。

将检测到的核磁信号的同向分量X、正交分量Y的表达式x(t)和y(t)分别离散化成x(n)和y(n)，经过NEO算法计算后，得到：

由公式(2)可知，相比于核磁信号，尖峰噪声的幅值会突增，因此只要进行合理的阈值设定，便可剔除尖峰噪声。鉴于较弱的尖峰噪声可能会同核磁信号的幅值相近，选择绝对中位差(median absolute deviation，MAD)方法设定阈值门限。以同向分量X为例，检测到多组数据x₁,x₂,…x_n，则其MAD定义为：

MAD＝median_i{|x_i-median_j{x_j}|}                (3)

其中，median为中位数，x_j表示在同一时刻下的全部数据，x_i表示每一组数据。定义阈值的上下门限为：

$\left\{\begin{array}{c}\mathrm{low}=\mathrm{median}\&times;0.01;\\ \mathrm{high}=\mathrm{median}+b\&times;\mathrm{MAD}.\end{array}\right.---\left(4\right)$

其中，b为可根据信噪比以及多次重复试验确定的常系数。经过阈值限定后，取其范围

在公式(5)中范围内的信号被保留，而尖峰噪声则被剔除，信号被剔除的部分由同一时刻其它组数据中的非坏值替代。去除尖峰噪声后的信号如图4所示。正交分量Y同上述处理过程相同。

c、基于主成分分析(Principal Component Analysis，PCA)方法，以下简称PCA消噪算法，分别对同向分量X、正交分量Y进行初步的信噪分离。

以同向分量X为例，PCA消噪算法的基本步骤如下：

首先，对于一维的时间序列x(t)来说，需要将其构造成一定维数的数据矩阵Q，才能执行PCA消噪算法。定义一个N×M的实对称矩阵

求取矩阵Q的协方差矩阵C_Q，并获得其特征值λ₀,λ₁,…,λ_N-1和特征向量A₀,A₁,…,A_N-1,由此构造正交矩阵

A＝[A₀,A₁,…,A_N-1]^T              (7)

使得

I＝AQ                     (8)

转化后的矩阵I的各个列向量依次为矩阵Q的第一主分量，第二主分量……。而且它们互不相关，并呈现方差递减的趋势。若方差累计贡献率a满足：

$a=\underset{m=1}{\overset{i}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_m/\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_m\&times;100\%>95\%---\left(9\right)$

其中，n为所有特征值的个数，i为所选取的主成分的个数。

即可由选取的主成分个数i对矩阵I和A降维，即保留矩阵I和A的前i行，分别得到矩阵I_i和A_i，重构信号矩阵为

$Q_i=A_i^T{I}_i---\left(10\right)$

单次的PCA消噪算法并不能够完全消掉具有强相关性的噪声成分，因此采用多次PCA循环消噪算法，若满足方差累计贡献率a>95％且主成分个数为1，或者是PCA适用性检验KMO<0.7，则判定PCA循环过程停止，否则循环将多次执行，直到满足条件为止。KMO与PCA适用性关系如表1所示，表1 KMO值与PCA适用性关系

KMO取值区间	(0,0.5]	(0.5,0.7]	(0.7,0.8]	(0.8,0.9]	(0.9,1]
						PCA适合度	极不适合	不大适合	一般	适合	非常适合

经过PCA处理后的信号如图5所示。对比图4和图5，可以看出，原先复杂的具有强相关性的噪声已经被大量消除，且经过PCA处理后的信号具有明显的理想核磁信号的曲线趋势。

d、基于经验模态分解(Empirical mode decomposition，EMD)，简称EMD方法分别对X分量、Y分量进行分解并提取趋势项，提取结果如图6所示，可以看出，无论是X分量，还是Y分量，均分解出了信号的趋势项。

经过PCA消噪算法处理过的FID信号噪声成分明显减少，但仍残留互不相关的工频干扰成分。因此对每一个X分量、Y分量分别进行EMD分解，提取趋势项。同样以X分量的一组数据x(t)为例，通过三次样条函数分别拟合出x(t)的极大值包络线e_max(t)和极小值包络线e_min(t)，根据法国学者Rilling提出的判据：

$\mathrm{\&sigma;}\left(t\right)=\frac{|e_{\max }\left(t\right)+e_{\min }\left(t\right)|}{|e_{\max }\left(t\right)-e_{\min }\left(t\right)|}---\left(11\right)$

设定三个阈值常量θ₁＝0.05，θ₂＝0.5，α＝0.05；相应的若满足以下两个条件，就认为得到的IMF分量满足要求。

①σ(t)中小于θ₁的比例不低于1-α，即

#{t∈D|σ(t)<θ₁}/#{t∈D}≥1-α                 (12)

其中，D为信号持续的总时间，#{A}表示在集合A中元素的个数。

②对于整个时间范围内的每个时刻t，满足：

σ(t)＜θ₂                            (13)

则认为EMD分解得到的IMF分量满足要求。将x(t)减去σ(t)就得到一个去掉低频的新信号得到：

$h_1^1\left(t\right)=x\left(t\right)-\mathrm{\&sigma;}\left(t\right)---\left(14\right)$

一般不是一个平稳信号，不能够满足上述两个判定条件，因此重复步骤d，假定经过k次之后(k小于10)，满足IMF分量的判定条件，则原信号x(t)的一阶IMF分量为：

$c_1\left(t\right)=\mathrm{im}{f}_1\left(t\right)=h_1^k\left(t\right)---\left(15\right)$

将x(t)减去c₁(t)，得到一个去掉高频成分的新信号r₁(t)，则：

r₁(t)＝x(t)-c₁(t)                     (16)

对去掉高频成分的r₁(t)重复处理信号x(t)的过程，得到第二个IMF分量c₂(t)，如此反复进行，直到第k阶IMF分量c_k(t)或其余项r_k(t)小于预设值；或当残余分量r_k(t)是单调函数或常量时，EMD分解过程停止。最后x(t)经EMD分解后得到：

$x\left(t\right)=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{c}_j\left(t\right)+r_k\left(t\right)---\left(17\right)$

其中，c_j(t)表示每一个IMF分量，r_k(t)表示分解余项，即信号的趋势项。

e、对EMD处理后的X分量、Y分量分别叠加求平均后获得e(t)曲线，如图7所示。从图7可以看出，经过此方法对检测到的核磁FID信号进行消噪处理后得到的e指数曲线与理想FID信号的e指数曲线十分吻合，验证了本发明方案的优越性。

核磁信号经过上述消噪算法处理后，分别得到n组相应的X分量和Y分量，即

x_i(t)＝r_xki(t)

y_i(t)＝r_yki(t)；(i＝1，2…n)               (18)

将其各自叠加求取平均值，得到：

$x\left(t\right)=\frac{x_1\left(t\right)+x_2\left(t\right)+...+x_n\left(t\right)}{n}$

$y\left(t\right)=\frac{y_1\left(t\right)+y_2\left(t\right)+...+y_n\left(t\right)}{n}---\left(19\right)$

最终可以求出e(t)曲线：

e(t)＝x(t)+i·y(t)                  (20)

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

Claims (6)

1.一种核磁共振测深FID信号噪声抑制方法，其特征在于，包括以下步骤：

a、对磁共振测深系统检测到的信号进行频谱分析，利用归一化正交检测技术将检测到的信号分解为同向分量X、正交分量Y，并进行硬件滤波处理得到低频的FID信号；所得到的同向分量X和正交分量Y表示为：

E(t)＝x(t)+i·y(t)

其中，x(t)为同向分量X的表达式，y(t)为正交分量Y的表达式，E₀为初始振幅，为平均衰减时间，δω是接收角频率与拉莫尔角频率ω₀的差值，而ε_x(t)和ε_y(t)分别为同向分量X和正交分量Y上的噪声成分；

b、对采集数据采用非线性能量算子算法分别剔除FID信号中同向分量X、正交分量Y的尖峰噪声；

c、基于主成分分析方法对步骤b处理后的同向分量X、正交分量Y分别进行初步的信噪分离，将与有用信号具有强相关性的噪声分离出来并消除；

d、基于经验模态分解方法对步骤c处理后的同向分量X、正交分量Y进一步分解提取信号趋势项；

e、对步骤d处理后的同向分量X、正交分量Y分别叠加求平均后获得e指数曲线。

2.如权利要求1所述的核磁共振测深FID信号噪声抑制方法，其特征在于，步骤b中将检测到的核磁信号的X分量x(t)、Y分量y(t)分别离散化成x(n)和y(n)，经过非线性能量算子算法运算后，得到：

选择绝对中位差MAD方法设定阈值门限，设一组数据为z₁,z₂,…z_n表示x(n) 或y(n)，则其MAD定义为：

MAD＝median_i{|z_i-median_j{z_j}|} 

其中z_j表示在同一时刻下的全部数据，z_j表示每一组数据，定义阈值的上下门限为：

其中b为可根据信噪比以及多次重复试验确定的常系数，经过阈值限定后，取其范围

在上述公式中范围内的信号被保留，而尖峰噪声则被剔除。

3.如权利要求1所述的核磁共振测深FID信号噪声抑制方法，其特征在于，步骤c中，首先，对于一维的时间序列x(t)或y(t)，用z(t)表示，定义一个N×M的实对称矩阵

求取矩阵Q的协方差矩阵C_Q，并获得其特征值λ₀,λ₁,…,λ_N-1和特征向量A₀,A₁,…,A_N-1,由此构造正交矩阵

A＝[A₀,A₁,…,A_N-1]^T

使得

I＝AQ

转化后的矩阵I的各个列向量依次为矩阵Q的第一主分量，第二主分量，以此类推，各个列向量呈现方差递减的趋势，若方差累计贡献率a满足：

其中，n为所有特征值的个数，i为所选取的主成分的个数；

由选取的主成分个数i对矩阵I和A降维，保留矩阵I和A的前i行，分别得到矩阵I_i和A_i，重构信号矩阵为

由重构信号矩阵得到消掉具有强相关性的噪声成分的时间序列x(t)或y(t)。

4.如权利要求3所述的核磁共振测深FID信号噪声抑制方法，其特征在于，步骤c中采用多次循环主成分分析方法，若满足方差累计贡献率a>95％且主成分个数为1，或者是主成分分析方法适用性检验KMO<0.7，则判定循环过程停止。

5.如权利要求1所述的核磁共振测深FID信号噪声抑制方法，其特征在于，步骤d，经验模态分解方法包括：通过拟合出x(t)或y(t)的极大值包络线e_max(t)和极小值包络线e_min(t)，选用G.Rilling sift criterion判据来获取IMF分量以及分解余项，该判据的评价函数为：

设定三个阈值常量θ₁、θ₂和α，若满足：

①σ(t)中小于θ₁的比例不低于1-α，即

#{t∈D|σ(t)<θ₁}/#{t∈D}≥1-α

其中，D为信号持续的总时间，#{A}表示在集合A中元素的个数；

②对于整个时间范围内的每个时刻t，满足：

σ(t)<θ₂

则认为经验模态分解得到的IMF分量满足要求；

当分解余项r_n(t)是单调函数或常量时，经验模态分解过程停止，最后经过经验模态分解后得到的信号x(t)或y(t)为：

或

其中，c_j(t)表示每一个IMF分量，r_k(t)表示分解余项，即信号的趋势项；分别得到n组相应的同向分量X和正交分量Y，表示为：

x_i(t)＝r_xki(t)

y_i(t)＝r_yki(t)；(i＝1，2…n)。

6.如权利要求5所述的核磁共振测深FID信号噪声抑制方法，其特征在于，将同向分量X和正交分量Y各自叠加求取平均值，得到：

最终求出去噪后的e(t)指数曲线：

e(t)＝x(t)+i·y(t) 。