引言

核磁共振(NMR)T₂谱可用于获取岩石孔隙结构、可动流体等信息^[1]，但实测T₂弛豫信号较弱，容易受到噪声影响^[2]，因而对NMR去噪处理较为普遍。采用小波变换对NMR信号去噪处理方法的关键在于小波基和阈值的选取，通常是根据噪声差异选择小波阈值或采用多尺度处理来压制噪声，并尽量保留图像边缘信息，但很可能破坏信号本身的数据特征^[3-5]。基于经验模态分解(EMD)的去噪方法是一种处理非线性非平稳信号的自适应方法^[6-7]，将EMD与小波变换相结合进行NMR去噪得到了更高的信噪比^[8-9]。2010年，蔡剑华等^[10]去除第一个本征模态函数，对第二和第三个本征模态函数设阈值，较好地保护了有用信号，使得去噪后弛豫信号反演T₂分布更精确，孔隙度计算结果一致性更高^[11-12]。然而，目前EMD核磁去噪方案很少考虑NMR噪声自身特性，信噪分离准则的科学性不足，设置阈值的方式也破坏了EMD算法的自适应性优势。

本文通过研究NMR噪声特性，提出一种基于EMD的NMR去噪方法，并使用岩芯分析和测井NMR衰减信号进行实验，评价算法的优越性。

1 核磁共振噪声特性 1.1 NMR噪声类型

使用NIUMAG-Ⅱ型低场核磁共振信号测量仪(主频2~MHz)，采用CPMG脉冲序列测量NMR横向弛豫(T₂)衰减数据。测量条件设定为：温度25 ℃，回波时间0.2 ms，等待时间6 s，回波个数2 000，叠加次数32。

在岩芯室空载时测量背景T₂衰减数据，即得到NMR噪声信号

$ S'({t_i}) = S({t_i}) - \frac{1}{P}\sum\limits_{i = 1}^P {S({t_i})} $

(1)

式中：S'(t_i)-在t_i时刻提取的NMR噪声；

S(t_i)-没有样品时在t_i时刻测量的衰减数据；

i-第i个回波；

P-回波数。

典型的核磁噪声信号的幅度概率分布密度图及高斯拟合如图 1所示。

图 1中，高斯分布拟合的复相关系数高达0.959，显示NMR噪声服从高斯分布。根据功率谱定义^[13]可得到NMR噪声信号的功率谱，如图 2所示。NMR噪声功率谱在噪声频率范围内保持稳定，符合白噪声的定义^[14]。

1.2 采集参数对NMR噪声的影响

在核磁数据采集中，主要设定参数包括回波间隔、迭代次数和等待时间。不同回波间隔核磁实验显示(图 3a)，噪声信号量与回波间隔存在幂函数关系

$ S = 26.476T_{\rm{E}}^{{\rm{ - 1}}.{\rm{560}}} $

(2)

式中：

S-噪声信号量；

T_E-回波间隔，ms。

回波间隔的减小扩大了核磁噪声信号量，但并不影响NMR噪声高斯分布的特征，如图 3b所示。不同叠加次数核磁实验结果显示，噪声信号量与采集参数呈正比(图 3c)

$ S=4.987N_{\rm S}+52.749 $

(3)

式中：N_S-叠加次数。

不同等待时间数核磁实验结果显示，NMR噪声信号量与等待时间关系不大(图 3d)。综合考虑各种采集参数与NMR噪声信号量的关系，可以对NMR噪声进行初步估计，见式(4)。

$ S = 4.988{N_{\rm{S}}} + 32.276T_{\rm{E}}^{ - 1.159} - 162.371 $

(4)

2 基于EMD核磁共振去噪方法 2.1 EMD分解原理与步骤

NMR噪声为典型高斯白噪声，可利用EMD方法进行去噪^[15-16]。EMD是一种适应性较强的时域分析方法，通过自适应筛选极值方式，将信号分解为一系列不同尺度的本征模态函数(IMF)，分解过程从高频到低频，最低频分量(余项)一般认为是信号的趋势，高频分量认为是噪声反映。EMD分解的一般步骤如下。

(1) 找出输入信号X(t)的所有局部极大值与极小值。

(2) 将输入信号的所有极大值点，使用三次样条插值拟合方法^[17]，求出其上包络线的均值，再使用同样的方法确定其下包络线的均值。

(3) 计算上下包络线的均值

$ {M_1}(t) = \frac{{{E_{\rm{u}}}(t) + {E_{\rm{l}}}(t)}}{2} $

(5)

式中：

M₁(t)-t时刻上下包络线的均值；

E_u(t)-t时刻上包络线的均值；

E₁(t)-t时刻下包络线的均值。

(4) 求取信号与均值的差

$ {H_{\rm{1}}}(t) = X(t) - {M_{\rm{1}}}(t) $

(6)

式中：

H₁(t)-t时刻信号与均值之差；

X(t)-t时刻的输入信号。

若在数据长度内，极值点数目与零点数目相等或至多相差一个，并且上下包络线的平均值等于0，则将H₁(t)作为一个本征模态函数。

(5) 从原始信号中减去本征模态函数，即

$ {R_1}\left( t \right) = X\left( t \right) - {H_1}\left( t \right) $

(7)

式中：R₁(t)-信号的剩余分量。

若R₁(t)是非单调或极值点数目为1，则分解结束；否则，继续执行步骤(1)~(5)，逐渐得到2~n个本征模态函数及信号的余项R_ES。

原始输入信号可以表示为

$ X(t) = \sum\limits_1^n {{H_i}(t) + {R_{{\rm{ES}}}}} $

(8)

2.2 核磁噪声及含噪信号的EMD分解

岩样实测NMR弛豫信号的离散程度较大，尖峰和无用毛刺较多，见图 4。对实测NMR弛豫信号进行EMD分解，见图 5，图 6。

由图 6可见，从IMF1到IMF5噪声(图中黑色曲线)和含噪信号(图中红色曲线)的本征模态函数一致性较强，从IMF6开始两种信号差异明显增大，有用信号在IMF6之后的本征模态函数中仍有分布，而噪声则几乎没有。

2.3 信噪分离准则 2.3.1 基于曲线趋势法的信噪分离准则

将余项(R_ES)与各个本征模态函数叠加，并与余项对比^[17]，见图 7。余项代表衰减数据的趋势，反映有用的弛豫信息。将余项与各个本征模态函数叠加，若某本征模态函数代表噪声，由于噪声的零均值分布特性，那么叠加后的曲线趋势应类似于余项^[18]；若某本征模态函数代表有用信号，那么叠加后的曲线趋势在信号起始处应与余项不同。

余项与IMF6~IMF9叠加的曲线趋势不同于余项；余项与IMF1~IMF5叠加的曲线除局部发生离散和跳变外，整体的趋势与余项基本一致，因此，认为IMF6~IMF9是有用信号而IMF1~IMF5是噪声成分。尽管曲线趋势法可以确定信噪分离准则，但分离过程主观性较强且无法做到程序化处理。

2.3.2 基于改进的过零点率曲线的信噪分离准则

传统的过零点率检测法^[19-21]以过零点数变化速度最快的IMF为信噪分离IMF，但实际上过零点数变化速度突变的IMF才是信噪分离IMF，因为他代表了噪声信号贡献的突然降低。因此，可以求出曲线在相邻两个IMF处的斜率，根据前后斜率比值反映噪声的突变情况，斜率比最大的IMF即为信噪分离IMF(图 8)。图 8中，最大斜率比对应IMF6，因此，判断IMF6~IMF9为有用信号。将IMF6~IMF9与余项叠加后可以得到NMR衰减数据的去噪信号，如图 9。

2.4 去噪效果分析

本文方法的去噪效果去噪前后各IMF1~IMF9、R_ES的相关系数分别为0.265，0.207，0.153，0.376，0.672，0.034，0.019，0.075，0.048和0.013；原始噪声信号量、EMD分解去除噪声信号量及利用采集参数估计的噪声信号量分别为286.278，205.687，275.786；小波阈值法、EMD-小波阈值法和本文方法的信噪比(SNR)分别为21.238，27.659，36.620。可见，噪声主要分布在前5个IMF，去除的噪声信号接近原始信号，本文方法去噪效果好于小波阈值法及EMD-小波阈值法^[22-25]。

表 1给出了去噪前后不同叠加次数下的信噪比，叠加次数低于16时，出现负的信噪比，说明数据质量较差，通过去噪处理后，使信噪比明显提高。叠加次数为32时，去噪后信号信噪比已达到36.6，此信噪比下的信号已达到生产要求。

进一步选取5块含泥质砂岩样品，对样品进行适当的预处理之后，测得NMR衰减信号并进行去噪处理，分别使用去噪前后的数据确定核磁孔隙度，结果见表 2。由表 2可见，去噪后效果明显，计算的孔隙度明显增大且与气测及水测孔隙度较为接近。

3 核磁测井信号处理

对XX井的NMR测井原始数据运用本去噪方法进行处理，处理井段为2 490.0~2 505.0 m，去噪前后的反演结果如图 10所示。在图 10a中，第1道是岩性曲线，第2道为深度曲线，第3道是长等待时间(12 s)短回波间隔(0.9 ms)的回波串反演结果，第4道是第3道回波串去噪后的反演结果。2 494.5 m处岩芯薄片、饱和水T₂谱(等待时间6 s，回波间隔0.2 ms，128次)、压汞孔隙半径分布图分别见图 10b，图 10c，图 10d。

试水资料显示，2 493.0~2 497.0 m处为一纯水层，NMR测井反演结果能较为准确地反映地层的孔隙结构。薄片分析显示岩芯孔隙类型为残余粒间孔，矿物颗粒分选性较好，缝洞率仅为2%，孔隙结构单一；岩芯饱和水T₂谱及由压汞测试得到的孔隙半径均呈明显的单峰形态。去噪后T₂谱呈单峰分布，与实验分析结果对应较好，这说明去噪后反演结果能够更准确地反映地层孔隙结构的真实情况。

4 结论

(1) NMR噪声是一种高斯白噪声，在已知叠加次数和回波间隔的前提下，可以用文中提供的方法对噪声信号量进行估计。

(2) 基于EMD的去噪方法被成功地运用在核磁共振衰减数据中，使用改进的过零点率曲线作为信噪分离准则，分离结果与曲线趋势法的结论及噪声的特性对应良好，可以在压制核磁噪声的同时提高数据的准确度。

(3) 针对EMD算法本身的改进研究是必要的，同时需要注意各种去噪方法在清除核磁噪声的同时是否消除了数据中的细节信息。在进一步的研究中，除改进去噪方法外，去除掉的细节信息应当得到关注。