李 辉， 肖忠祥， 刘晓娟

(西安石油大学，陕西 西安 710065)

1946年哈佛大学的Purcell 和斯坦福大学的Bloch 两位教授分别用吸收法和感应法发现核磁共振(NMR)现象以来，经过近60年的发展，核磁共振T2 谱已经能够提供许多的岩石储集物性参数及流体物性参数(肖立志，1998)，如岩石孔隙度、有效孔隙度、孔隙尺寸分布、渗透率、可动流体与束缚流体饱和度、粘度、流体扩散系数等。而这些参数都是由测量的岩石样品及流体样品的核磁共振弛豫信号经过多指数反演后的T1、T2谱计算得到的。这些岩石物理信息为储层评价、产量预测等提供了重要的理论依据。

目前，国内外先后提出的核磁共振T2谱反演方法主要有三种:罚函数法(BRD)(Butler et al.，1981)、奇异值分解法(SVD)(王忠东等，2003)以及联合迭代重建算法(SIRT)。

早期解谱方法多为单指数或多指数拟合，Whitall 等(1989)提出NNLS 算法，美国Numar 公司提出的商业化的基于正则化方法的UPEN 方案，Gy.Steinbrecher 提出的基于Laplas 变化的反卷积法(姚绪刚等，2003)，这些方法都存在着各种缺点，反演结果不理想，需要做进一步改进。

本文通过对传统奇异值算法进行改良，实现低信噪比的可操作性。

1 核磁共振T2 谱反演原理

1.1 基本原理

核磁共振测井的原始数据是幅度随时间衰减的回波信号，这些回波信号隐含重要的岩石物理和油气信息，并且能够反映横向弛豫过程，必须经过一个基本的多指数变化，得到反映孔径特征的T2谱分布，这一过程称为核磁共振测井的解谱。如果核磁共振检测到的信号很微弱，便会出现信噪比低的情况，在解谱过程中必须考虑这一问题，提高信噪比，突出有用信号。

理论和实验研究都表明，单个孔隙内磁化强度信号的横向弛豫过程是按照单指数规律递减的(Coates et al.，2007)。测井过程中，由于探测区域存在一系列大小不同的孔隙群体，孔隙中存在多种弛豫组分，测量得到的信号就是这些不同弛豫组分信号的迭加，即总磁化强度(Total Magnetization)是一系列指数衰减信号的迭加(肖立志等，2010)(图1，2)。

1.2 计算过程

图1 指数衰减信号曲线Fig.1 Exponential decay signal curve

图2 NMR 原始信号—回波串Fig.2 The original NMR signals-Echo string

T2谱横向弛豫信号分布满足Fredholm 第一类积分方程，但实际过程中，由于噪声的存在，弛豫过程满足下面的方程:

其离散化模型为:

其中，Z(tj)为tj观测到的回波幅度;T2i为第i个弛豫分量的横向弛豫时间;T2min和T2max分别为测量的自由感应衰减信号所能辨别的最短和最长弛豫时间;xi为第i 种弛豫分量零时刻的信号大小;n 为弛豫分量个数;tj为时间(通常为回波间隔的整数倍);εi为噪声信号。

解谱过程就是如何从(2)中分解出T2i和xi。

弛豫分量的横向弛豫时间T2i一般取2i+1(i =1，2，…)ms。假设弛豫分量的个数为n，观测的回波个数为m，则由式(2)写成矩阵形式如下:

其中

在式(3)左右两边同时左乘矩阵的广义逆矩阵A+，可得

其中，A+为n × m 矩阵，且由最小二乘法知A+=(ATA)-1AT。

根据矩阵理论的奇异值定理(肖立志，1998)可知，对于A ∈，r(A)= r，存在m 阶酉矩阵U和n 阶酉矩阵V，使

其中，∑r= diag(δ1，δ2，…，δr)，δ1，δ2，…，δr是A 的全部非零奇异值，且呈递减排列。为m × n矩阵，于是有，

不论A 是否奇异，总可以分解成(5)形式，并且∑r是唯一的。但通常情况下，由式(7)解出的X 值不稳定，且常有负数出现。根据核磁共振的解谱结果的物理意义可知，xi为第i 种弛豫分量在零时刻的初始信号，不能为负值，而一般纯碎的数学过程不能保证这一点。

1.3 数值解稳定性控制

本文采用对系数矩阵A 加阻尼项的方法使数值解稳定，并结合迭代算法进行非负性约束。

1.3.1 阻尼方法

根据式(3)构造函数

其中，λ2称为阻尼因子。为求解F 最小值，令F 梯度为零，有

所以

又因

所以

对式(12)求逆矩阵得

综上可得

另外，给定x 一个初始解X = X0，得到Z0=AX0，与Z = AX 相减，可得Z - Z0= A(X - X0)，记为

与上述推理类似，能够得出

同时，需对ΔX 进行截断处理，得

信噪比SNR 是由原始资料通过下述公式确定的

而由X = X0+ ΔX 即可得到X。

最终，当SNR 降低时，对角线上保留较少项数，反之，SNR 提高时，对角线上保留较多项数。从而减少了噪声对结果的影响，增加了解的稳定性，在不同信噪比情况下都可以得到相对稳定的T2谱。

1.3.2 非负性约束实现

(1)给定初始值X(k)，令k = 0;

(4)计算X(k+1)= X(k)+ ΔX(k);

(5)若X(k+1)中分量全不为负，则终止，输出X(k+1)，否则转(6);

(6)令X(k+1)中负值为0，k = k +1，转(2)再次迭代，直到满足条件(5)。

2 反演效果检验

为验证该方法效果，用微机编制相关程序，并对其进行数值模拟。图2 中实验T2谱的布点数为32，其余实验中T2谱的布点数为64，布点区间为0.1 ～10 000 ms。

反演结果和构造T2谱对比结果(图3，4)。其中图3 为32个布点的反演无噪回波数据T2谱与构造谱对比图，图3 为64个布点反演无噪回波数据T2谱与构造谱对比图。可以看出，反演的曲线和数据点与构造谱的曲线和数据点基本重合。这表明，本文的奇异值分解法可以较准确反映真实的T2谱。

分别加入信噪比SNR =30，SNR =60，SNR =100 的高斯白噪声后的数据反演结果对比(图5)。从图中可看出，三种情况下反演结果与构造谱分布趋势一致，谱线光滑并连续，具有较高吻合度。同时还可以看出，信噪比越高，反演T2谱与构造T2谱的符合率越高。该算法具有较高的稳定性和收敛性。

用该方法处理某岩心实验室NMR 数据，解谱结果与实验室解谱结果对比(图6)。可以看出，岩心奇异值分解法反演结果与实验室反演结果吻合性较好，能够得到连续的T2谱。

3 结论

此种SVD 法能够获得连续性较好的反演T2谱，算法较稳定，收敛性好，非负约束条件容易实现。反演T2谱与构造T2谱吻合性较高。由NMR实验室数据反演结果可知，该方法可以应用于实验室岩心的T2谱反演。

图3 无噪声反演结果对比(32个布点)Fig.3 Inversion results of NMR data without noise and origin T2 spectrum(32 stationing)

图4 无噪声反演结果对比(64个布点)Fig.4 Inversion results of NMR data without noise and origin T2 spectrum(64 stationing)

图5 有噪声反演结果对比Fig.5 Inversion results of three groups of NMR data with noise

图6 岩心反演结果对比Fig.6 Inversion results of NMR data from core using Improved SVD and Laboratory solution spactrum

王忠东，肖立志，刘堂宴.2003.核磁共振驰豫信号多指数反演新方法及其应用［J］.中国科学(G 辑)，33(4):323-324.

肖立志，柴细元，孙宝喜，等.2010.核磁共振测井资料解释与应用导论［M］.北京:石油工业出版社:15-38.

肖立志.1998.核磁共振成像测井与岩石核磁共振及其应用［M］.北京:科学出版社:11-18.

姚绪刚，王忠东.2003.一种新的核磁共振弛豫谱反演算法［J］.测井技术，27(5):373-376.

Butler J P，Reeds J A，Dawson S V.1981.Estimating Solutions of First Kind Integral Equations with Nonnegative Constraints and Optimal Smoothing［J］.SIAMJ Numer Anal.，18(3):381-397.

Coates G，肖立志，Prammer M. 2007. 核磁共振测井原理与应用［M］.北京:石油工业出版社:39-41.

Whitall K P，Mackay A I.1989.Quantitative Interpretation of NMR Relaxation Data［J］. Journal of Megnetic Resonance，4(5):134-152.