二维核磁共振识别与评价储层流体是当今测井技术的发展方向(谭茂金和邹友龙，2012; 赵文杰和谭茂金，2008).目前，普遍应用的核磁共振测井只基于横向弛豫时间T₂分布做油藏估计和储层评价，称为“一维核磁共振测井”.然而，实际的地层环境中，储层孔隙结构复杂、原油成分多样、黏度变化范围大，造成不同流体信号在T₂分布上重叠(谢然红和肖立志，2009；谢然红等，2005)，特别是短弛豫时间的黏土束缚水和毛细管束缚水与黏度较大的原油在T₂信号上重叠，导致目前一维核磁共振测井不能彻底区分来自水和稠油的信号.二维核磁共振可从扩散系数、弛豫时间两个维度上展现流体的性质，将在一维空间上混叠的信息分开，因此，在快速确定流体类型、原油黏度及含油饱和度方面，具有理论优势(Sun and Dunn， 2004).二维核磁共振测井给岩石物理勘探开辟了广阔的领域，它利用梯度磁场中增强 T₂ 弛豫以及油、水扩散系数的差异生成能明显区分油、水的二维NMR图，称为(T₂，D)分布图(肖亮，2006).

目前，世界上主流石油勘探公司的核磁仪器，如Baker Atlas的MR-Explorer和Schlumberger的 MR-Scanner 都可进行二维核磁核磁测井.最早的二维 核磁共振测井脉冲序列由Freeman等(Freedman et al., 2001; Freedman et al., 2002)提出采用不同的回波间隔T_E，测量多组CPMG自旋回波，经反演得到(T₂，D)二维分布.Schlumberger公司Doll研究中心的Hurlimann(2002)和Venkatarama(2002，2003)以及美国Chevron石油公司的BoqinSun和Dunn(2005)分别提出利用2个窗口改进自旋回波脉冲序列CPMG，实现了(T₂，D)分析的2D NMR测井，极大地提高了流体识别能力和饱和度计算精度.国内，谢然红、肖立志等(2009)对含有顺磁物质的人造砂岩和天然泥质砂岩饱和水，进行2D NMR实验测量，给出了岩石横向弛豫时间-内部磁场梯度(T₂，G)的二维分布图，研究结果对分析陆相沉积地层复杂岩性具有重要指导意义.

在反演方法方面，目前主要采用基于截断奇异值分解(Truncated Singular Value Decomposition，TSVD)的算法(顾兆斌和刘卫，2007)实现了不同流体模型的二维谱反演.然而，这种方法的评价效果很大程度上取决于保留奇异值个数，且后者的选取需基于先验信息(所测地层流体类型与分布特性的经验资料)，此条件在实际测井应用中往往无法满足(Venkataramanan et al., 2002).此外，TSVD反演算法需要对大型矩阵做奇异值分解，运算量大，严重影响二维核磁共振测井的评价效果.

基于此，本文在阐述二维核磁共振测井原理与数理模型的基础上，提出一种基于SVD(Singular Value Decomposition)和BRD(Butler et al., 1981)的正则化反演方法.该算法通过压缩、平滑、迭代优化等步骤获取流体二维分布信息，较大幅度减少计算量，无需先验信息，自动寻求最佳优化因子，为二维核磁共振测井数据实时处理提供必要的基础.数值模拟与在自行研发的二维核磁共振测井仪器上的初步实验表明，本文所设计的反演算法可有效取流体的(T₂，D)分布. 2 理论分析 2.1 二维核磁共振测井时序

理论上，获取地层流体分布的测量时序较多.例如，采用反转恢复(Inversion recovery)时序可对纵向弛豫时间T₁编码；利用CPMG序列可测量横向弛豫时间T₂；扩散信息则通过改变回波间隔或外部磁场梯度获取.在实际的应用中，受测井速度与固定外部磁场梯度的限制，目前通常采用图 1所示的两种测量时序对流体(T₂，D)二维分布编码.

图 1

Fig. 1

图 1(T2，D)二维核磁共振测量时序 (a)固定梯度改变回波间隔测量多组CPMG；(b)固定梯度两窗口改进型CPMG. Fig. 1(T2，D)2D NMR Logging Sequence (a)Using different echo spacing and a constant field gradient in a regular CPMG for encoding， and (b)constant field and two windows CPMG for encoding.

图 1a中使用了固定梯度、标准CPMG、不同回波间隔(t_E)编码的方法.为彻底将油和水分开，这种方法中需采用很长的回波间隔，但当回波间隔过大时，受仪器精度限制，仅能测得前面几个回波，无法获取长扩散时间的信号.图 1b所示为两个窗口改进型CPMG脉冲序列，该时序中，第一个窗口采用不同长度的回波间隔给扩散信息编码，第二个窗口记录最小回波间隔下的弛豫信息.第一个窗口中的回波个数较少(一般为两个)，回波间隔相对较长，称为长回波间隔t_EL.在此序列中，能提供扩散率信息的独立变量是前两个回波的间隔.由于扩散的影响，增加时序中头两个回波的间隔将减小余下回波的幅度.余下的序列的获取采用最小回波间隔(t_ES)的脉冲序列，目的是为了将扩散的影响降到最低.在一系列不同的长回波间隔t_EL序列中，头两个回波的幅度衰减反映了流体扩散率的信息；而余下回波幅度的衰减形态提供驰豫时间分布特性.图 2描述了两窗口改进型CPMG序列中采用3组不同长回波间隔时信号幅度衰减情况.

图 2

Fig. 2

图 2 两窗口改进型CPMG序列信号幅度衰减 Fig. 2 The echo decay string in two windows CPMG

两窗口改进型CPMG序列，克服了“回波间隔过长”缺陷，极大地提高了流体识别效果和饱和度计算的精度，其回波幅度的衰减方程如下式所示：

令

表 1(Table 1)

表 1 测量时序中各参数个数 Table 1 The parameter number in logging sequence 长回波个数 标准回波个数 扩散系数布点数 弛豫时间布点数 NtEL=30 NE=1000 NT2=128 ND=128

表 1 测量时序中各参数个数 Table 1 The parameter number in logging sequence

为实现二维谱的实时反演，摆脱对被测流体分布信息的依赖，本文提出了一种效率较高的非负性约束反演算法.根据以上分析，可采用Tikhonov正则化方法，将反演问题转化为寻求满足式(11)的密度函数的优化问题.

(1)数据压缩

K₁ 的奇异值分解为

(2)非负性正则化反演

同前一样，式(14)可变形为

(3)确定最佳平滑因子

式(24)中描述的是在给定正则化因子 α 下取得密度函数的最优解.为进一步提高反演精度，需确定最佳平滑因子.常用的平滑因子选择方法有，Morozov残差准则，L曲线(Hansen，1992)和GCV方法等，本文采用BRD算法的扩展方法.令 f_T 为未知的扩散-弛豫分布的真实值， f_α 为利用步骤(1)、(2)估算出的分布函数值.根据BRD算法， α 最优值的确定通过最小化式(25)中的 S(α) 获得.

给定平滑因子 α 的初始值，利用(23)式计算出 c， 再通过(31)式更新平滑因子的取值，如此反复迭代，使残差达到期望水平，直至获得优化后的平滑因子 α_opt. 3 数值模拟与实验 3.1 数值模拟

为验证所设计时序和反演算法性能，构造具有双峰结构的二维图谱模拟含有束缚水和稠油的地层，其中稠油的弛豫时间 T₂ 为4 ms，扩散系数 D 为2×10^－11cm²·s^-1；束缚水的弛豫时间 T₂ 为10 ms，扩散系数 D 为2×10^－9cm²·s^-1.选择固定磁场梯度g=13.2 G/cm；短回波间隔t_ES=0.2 ms，回波个数为1000；长回波间隔t_EL在区间1~30 ms上线性均匀布点，布点数为30；横向弛豫时间 T₂ 在区间0.001~10 s上对数均匀布点，布点数为128；扩散系数 D 在区间10^－8~10^－12cm²·s^-1上对数均匀布点，布点数为128. (T₂，D) 分布采用高斯函数模型，利用(1)式产生模拟数据，在Matlab中，绘制模拟流体的 (T₂，D) 的二维谱分布的等高线、三维分布图，扩散系数和弛豫时间的一维分布，如图 3a，同时在图 3b中绘制了其中5条对应的回波衰减曲线.由图 3可以看出，模拟的稠油和束缚水在 T₂ 分布上几乎是重叠在一起.

图 3

Fig. 3

图 3 模拟流体的(T2，D)二维分布与回波串信号 (A)模拟流体(T2，D)二维分布模型;(B)模拟流体回波串信号. Fig. 3 The(T2，D)distribution and echo signal of fluid model (A)The(T2，D)distribution of fluid model;(B)The Echo signal of fluid model.

在Matlab中，分别编程实现TSVD与本文设计的非负性正则化反演算法处理模拟数据，并在程 序中加入计时模块，记录程序执行时间.将模型原始数据 X _mod 与反演结果 X _inv 代入(32)式，计算相对误差R_X.

TSVD与扩展型正则化算法反演结果分别如图 4和图 5所示.TSVD算法的处理时间为1824.35 s，反演曲线偏离模型曲线幅度较大，经计算，横向弛豫 时间与扩散系数的相对误差分别为6.66%与19.02%. 由于原始数据有效压缩，且仅通过14次迭代便求解出最佳平滑因子，本文设计的改进型反演算法整个处理过程在143.27 s内完成，而且反演曲线几乎与模型曲线重合，扩散系数与横向弛豫时间的相对误差仅为0.08%与0.00%.可见，基于SVD与BRD 的正则化算法相对传统的TSVD方法，在处理速度与计算精度上均有了大幅提高.

图 4

Fig. 4

图 4 TSVD反演结果 (A)二维反演谱;(B)反演结果一维对比. Fig. 4 TSVD inversion result (A)2D inversion spectrum;(B)Inversion result comparison in 1D.

图 5

Fig. 5

图 5 基于SVD与BRD的Tikhonov正则化反演结果 (A)二维反演谱;(B)反演结果一维对比. Fig. 5 The inversion result of Tikhonov regulation algorithm based on SVD and BRD (A)2D inversion spectrum;(B)Inversion result comparison in 1D.

向模拟的模型数据中混入不同幅度的加性高斯白噪声，得到不同信噪比的原始数据，并利用本文设计的正则化算法分别对其反演.当信噪比SNR分别为150、100、50时得到的 (T₂，D) 二维分布如图 6所示.当原始数据SNR=150时，反演的二维谱与原始值的基本相同，稠油和束缚水明显分离，信号较小的地方出现小幅度误差，导致分布二维分布中出现一些毛边，弛豫时间与扩散系数分布基本与模拟数据吻合，相对差分别为2.63%和3.67%；SNR=100时，稠油和束缚水清晰可辨，反演的二维图谱出现一些形变且毛边较为严重，弛豫时间与扩散系数相对误差分别为3.47%和7.52%；SNR=50时，反演的二维谱变形较为严重，稠油与束缚水分离界限模糊，弛豫时间与扩散系数相对误差较大，分别为6.00%和10.99%.

图 6

Fig. 6

图 6 不同信噪比下基于SVD与BRD的Tikhonov正则化算法的反演结果 (A)SNR=150，二维反演谱;(B)SNR=150，反演结果一维对比;(C)SNR=100，二维反演谱;(D)SNR=100，反演结果一维对比; (E)SNR=50，二维反演谱;(F)SNR=50，反演结果一维对比. Fig. 6 The inversion result of Tikhonov regulation algorithm based on SVD and BRD with different (A)SNR=150，2D inversion spectrum;(B)SNR=150，Inversion result comparison in 1D;(C)SNR=100，2D inversion spectrum;(D)SNR=100， Inversion result comparison in 1D;(E)SNR=50，2D inversion spectrum;(F)SNR=50，Inversion result comparison in 1D.

为进一步验证上述时序与反演算法的性能，配制一定浓度的CuSO₄溶液灌入刻度筒中，并在自行 研发的核磁共振测井仪EMRT(Elis Magnet Resonance Tool)上实现CPMG-DE脉冲序列，采集CuSO₄溶液二维磁共振信号.EMRT为中海油田服务股份有限公司(China Oilfield Services Limited，COSL)自主研发的国内首支现场成功应用的核磁共振测井仪器(目前仅限于一维核磁共振测井).它可挂接于中海油服的ELIS测井地面系统上，主要由探头、电路和储能三个短节构成，并配有泥浆排除器与偏心器，以保证其作业质量.仪器结构与实验平台分别如图 7a和图 7b所示所示.测试前，通过岩芯分析仪测量得所配制CuSO₄溶液的扩散系数与弛豫时间分别为1.205×10^－6cm²·s^-1和61.67 ms.利用CPMG- DE序列，采用表 2中的5组时序参数，测量仅含CuSO₄溶液的样品.实验测得回波串衰减曲线如图 7c所示.随着各组时序参数中t_EL逐渐增大，扩散和弛豫对回波的影响也越来越严重，对应回波幅度衰减 加剧，相同位置的回波幅度随t_EL增大而依次降低.

图 7

Fig. 7

图 7 基于EMRT的(T2,D)测量实验 (a) EMRT仪器框图；(b) EMRT 实验仪器； (c) 回波衰减曲线. Fig. 7 (T2,D) distribution measurement experiment on EMRT (a) EMRT structure graph, (b) EMRT experimental apparatus, (c) Echo decay curves.

采用基于SVD与BRD的正则化反演算法，处理CPMG-DE获取的回波串数据，由于被测样品处在屏蔽较好的刻度箱中，且样品成分单一，信噪比较实际测井时大，所以测量时序仅设计了5组不同的t_EL(模拟实验采用t_EL个数为30).仅经过5次迭代，耗时10.74 s便得到CuSO₄水溶液的扩散系数与弛豫时间的分布，如图 8a所示.可以看出，所得反演分布在两个维度上(T₂，D)均成单峰性，且随着与特性峰值点距离增大，信号幅度迅速衰减，符合样品仅含CuSO₄溶液的特性.通过计算，扩散系数与横向弛豫时间的估计值分别为1.176×10^-6 cm^-2·s^-1和59.23 ms，相对误差分别仅为2%和4%.同时，本文也采用TSVD算法处理实验数据，作为对比.由于未采用压缩算法，直接对大型矩阵做奇异值分解，TSVD算法的处理时间为60.73 s，其处理结果如图 8b所示.由于舍弃奇异值损失信号，反演结果出现一定误差，在T₂一维分布图中，除特性峰值外，在6.542 ms 处也出现了一个峰；在D一维分布中，特性峰值点以外的信号未得到充分抑制.TSVD计算出扩散系数与横向弛豫时间分别为1.311×10^-6 cm^-2·s^-1和53.15 ms，相对误差分别为14.8%和13.8%.由此可见，本文提出的算法在处理效率和性能上，均具有较大优势.

图 8

Fig. 8

图 8 CuSO4水溶液的(T2,D)分布 (A) 基于SVD和BRD的正则化反演算法处理结果; (B) TSVD反演结果. Fig. 8 The CuSO4 solution′s (T2,D) distribution (A) The result of Tikhonov regulation algorithm based on SVD and BRD; (B) The result of TSVD.

二维(T₂, D)核磁共振方法能有效识别储层流体，但国内的相关研究才刚刚开始，大量的应用基础问题需要探索.本文深入研究二维核磁共振测井理论，系统分析二维脉冲序列CPMG-DE测量扩散系数与弛豫时间的方法，结合核磁共振二维谱数理模型，提出一种基于SVD与BRD的非负性正则化反演算法.通过数值分析与实验验证得到以下结论：

(1)当油水同层，束缚水与稠油的T₂谱上重叠，但利用(T₂，D)二维分布，二者可在扩散系数上区分；

(2)本文所设计的反演算法通过压缩、优化、迭 代等步骤降低对先验信息的依赖、提高运算效率、减小相对误差，且在原始数据信噪比低至50时，仍可有效获 取流体(T₂，D)二维分布，较传统的反演算法(如TSVD)在二维核磁共振实时数据解释中具有较大优势；

(3)在自主研发的核磁共振测井仪上，测量刻度箱中CuSO₄溶液(T₂，D)分布.结果显示，该算法对弛豫时间和扩散系数的反演误差分别为2%和4%，但对实际油井的测量性能，尚待进一步验证.