基于Backus-Gilbert理论的孔隙介质核磁共振弛豫反演

肖立志，张恒荣，廖广志，傅少庆，李 奎

油气资源与探测国家重点实验室中国石油大学，北京 102249

1 引 言

岩石的NMR弛豫测量能够提供岩石物理特性和流体特性的许多信息，如岩石孔径分布、孔隙度、束缚水、可动流体以及渗透率等.NMR测井从井下采集到的是一系列反映地层弛豫特性的回波串，需要经过回波串多指数拟合才能得到T2谱.回波串拟合是一个数学反演过程[1-3］.国内外学者提出了多种反演方法，如Butler等[4］基于非负约束的模平滑方法（BRD）；Prammer[5］基于奇异值分解的反演方法（SVD）；Dunn等[6］对比了模平滑方法和曲率平滑方法的优缺点，发展了奇异值分解反演方法；Liaw等[7］和 Miller等[8］分别提出基于B样条函数和Gamma函数的基函数方法；Borgia和Brown[9］设计了一个“曲率平滑”的复杂方法；Rafael和Boqin[10］等提出了蒙特卡罗最优反演方法.国内王才志[11］，姜瑞忠等[12］尝试了奇异值分解方法；王忠东等[13］尝试了基于整体迭代修正的反演方法；翁爱华[14］尝试了高分辨率反演方法，谭茂金等[15］尝试了基于遗传算法的反演方法，马建海等[16］尝试了迭代Tikhonov正则化反演方法.这些算法尝试从不同角度解决核磁共振弛豫测量的多指数反演问题，但通常会受低信噪比、正则化参数、迭代收敛等的影响.

本文基于Backus-Gilbert理论，尝试一种新的孔隙介质NMR弛豫反演方法，简称BGW，其核心是从解的非唯一性出发，构造出一系列解估计，然后通过引入解估计分辨率和解估计方差对解进行评价，最终在其中找出一个最佳折衷的解估计.实践表明，该方法不仅稳定，并且适应于低信噪比数据.

2 孔隙介质NMR弛豫反演问题描述

NMR测井采集包含数百个乃至数千个回波的CPMG自旋回波串信号.经过一些预处理，如相位旋转和累加等，然后进行多指数反演.回波测量值M（ti）可以写成：

式中M（ti）表示第i个回波幅度（i=1，…，n），ti是第n个相等回波间隔的时间，εi代表第n个回波的噪声，并有m个已知的预先选择的弛豫时间T2在对数刻度上等间隔分布，a（T2j）是由方程（1）解出的弛豫时间T2j对应的幅度，可以刻度成区间孔隙度.

从数学角度，求解a（T2j）分布是一个反拉普拉斯变换[17-18］.噪声的存在会导致反演结果的非唯一性.一个可行的方法是，利用最小二乘拟合使下面的求和最小化：

3 基于BG线性评价理论的反演方法

BG反演理论是地球物理反演的重要成果，其核心思想是从反问题解的非唯一性出发，对有限而精确数据的反演，反问题的解虽不唯一，但解的平均可以是唯一的；对有限带误差数据的反演，反问题的解可能不存在，即使存在也不唯一，这时要用分辨率和方差来评价解估计.由于解估计的方差和分辨率不可能同时达到极小，只能通过修改折衷参数的办法求在分辨率和方差之间取最佳折衷的解估计[19-21］.

BG反演理论讨论的是连续模型，与BG理论对应的离散模型反演由Wiggins和Jackson先后提出[22］.

NMR数据反演问题是离散非线性反问题，采集到的回波串是有限带误差的数据.方程（1）的数值解将该方程变换为线性方程组的形式如下：

其中Aij，fj和gi是它们各自相应的连续函数的离散值，gi为带误差的回波幅度值，Aij为预先设置基函数组成的系数矩阵，fj为待定值，即待求的T2分布.

根据奇异值分解定理，n×m实矩阵A可分解为n×n的正交阵U和m×m的正交阵V和一个对角阵S：

由于系数矩阵A的病态性，其中sq+1，…，sm很小，根据离散BG理论，模型的协方差矩阵为

此时数据分辨率矩阵为

对应分辨率矩阵的解估计展布为

从理论上讲，希望方差最小，分辨率最高.根据BG理论，由于解估计的方差和分辨率不可能同时达到极小和极大，只能通过修改折衷参数的办法求在分辨率和方差之间取最佳折衷的解估计，而保留奇异值个数q则为折衷参数.随着折衷参数q的变化，解估计以及解估计的方差和展布也随之变化，q从小到大逐渐变化，解估计的方差越来越大，展布越来越小，然后用（6）和（8）式计算方差和展布，以展布为横坐标，方差为纵坐标绘出最佳折衷曲线图，最佳折衷的q值应位于最接近原点的左下方，因为原点处展布和方差都为零.然后根据求出解估计中第k个分量，解估计中每一个分量都对应一个最佳折衷q值，而q值可能相同，也可能不相同.

传统正则化方法是对整个fj向量直接求取，在选择正则化参数时无论是人工还是自动都是对整个向量进行处理.例如，SVD方法是取一个截断奇异值应用到每个分量；BRD方法是取一个平滑值加到每个分量；SIRT是基于整个向量的迭代修改.而BGW方法是分步求取每个分量，对每个分量都构造出一系列解估计，然后引入解估计分辨率和解估计方差来评价所有解估计，找出最佳折衷解估计，因此本方法在精度上更胜一筹，只是由于每个分量都要依次解出，反演时间上会有所增加.

4 基于BG理论的NMR弛豫反演方法

首先给定一个具有双峰特征的T2分布模型，如图1中的实线，根据这个模型生成回波串，设回波间隔TE=1.2ms，回波个数NE=1024，对数均匀布点32个，T2=0.3～3000ms，加入信噪比为100的高斯分布随机噪声（信噪比定义为第一个回波幅度与噪声标准差的比值），用此模型考察BGW反演方法.

根据BG理论的反演方法，求出每一个弛豫分量的方差和展布，本文以第5个弛豫分量为例，将其方差和展布绘成折衷曲线，如图2所示.

图2中横坐标（线性）为解估计的展布，纵坐标（对数）为解估计的方差；最佳折衷参数为6，计算此点与原点的距离最小.32个弛豫分量可以依次计算出与最佳折衷参数对应的最佳折衷解，即可得到T2分布.T2反演结果如图1的圈实线.可见反演结果与模型吻合很好，从而证实了该方法的有效性.

下面通过构造不同模型，如大孔单峰、小孔单峰、大孔占优双峰和小孔占优双峰，并加不同噪声，来考察BGW方法的适应性.

从图3中看出，从单一孔隙结构到复杂孔隙结构，从信噪比高达100到信噪比低达10，BGW均能给出较好的结果.

5 方法对比

将BGW方法分别与SVD、BRD以及SIRT方法反演结果进行对比，仍然采用大孔占优T2分布模型，信噪比SNR分别为100，30，20，10（见图4）.

各个反演方法计算结果的相对误差如表1表示，BGW方法在各信噪比下的相对误差都是最小的.

图4 大孔占优模型不同信噪比不同方法反演结果（a）BGW；（b）SVD；（c）BRD；（d）SIRT.Fig.4 Inversion results of dominant large pore model with different SNR（a）BGW；（b）SVD；（c）BRD；（d）SIRT.

图5是模拟某地层NMR T2分布.第一道是模型，在此模型基础上产生回波串，并加入噪声，使其信噪比为10，然后分别用四种方法进行反演.第二道是BGW反演结果；第三道是SIRT反演结果；第四道是SVD反演结果；第五道是BRD反演结果.可见，在信噪比为10的情况下，BGW反演结果与模型最为接近，说明在低信噪比情况下，BGW方法能得到更为有效的反演结果.

表1 4种方法反演结果的相对误差比较（%）Table 1 Comparison of the relative error of inversion results（%）

图5 四种方法反演结果与地层模型的对比Fig.5 Comparison of formation model computed by four methods

6 BGW方法影响因素研究

由式（1）知，不同噪声、回波间隔、回波个数及T2布点都可能影响反演结果.

6.1 噪声

信噪比低是NMR测井的重要特征，所以，反演方法必须考虑它对噪声的依赖性.

噪声影响可以用 Hansen等[23］给出的离散Picard条件进行分析.离散Picard准则认为，若方程组Ax=y的傅里叶系数，〈ui，y〉趋于零的速度在平均意义下快于矩阵A的奇异值σi趋于零的速度，则称该方程组满足离散Picard准则（条件）.

Picard准则是检验模型的不适定性以及模型受噪声污染情况的重要准则.对于一个给定的模型，若其满足或部分满足Picard准则，则可以通过适当的处理求解之；当模型的不适定性较强，或者数据受噪声污染较严重，以至于完全违背这一准则，那么再好的方法对该模型也是徒劳的.

图6显示了信噪比分别为106、30、10、5、2.5时傅里叶系数的变化.可以看出，傅里叶系数随信噪比的降低而逐步偏离奇异值曲线，即傅里叶系数趋于零的速度逐步慢于奇异值趋于零的速度.从本文大孔占优双峰模型中，信噪比低于5后，已经不满足离散Picard准则，也就是说在这个模型下，由于噪声的影响，从数学角度已经不能反演出真实可靠信息，必须要对数据滤波或者重新采集.

图7是BGW方法反演结果，证明信噪比小于5以后反演结果不再具有实际意义.

6.2 回波间隔影响

NMR回波串数据在时间域采集，按回波间隔TE记录N 个离散回波幅度，总时间为T=N×TE.分辨T2谱的能力依赖于信噪比（SNR）和采集模式.定义一个有效采集时间ΔT（T2）为T2组分从信号最强衰减到噪声水平的时间间隔.则相应的有效回波个数NE（T2）=ΔT（T2）／TE，T2谱的分辨能力可以用下面公式表示[24］：

从式中可以看出：（1）信噪比越低，T2谱分辨能力越差；（2）相对于长弛豫组分，短弛豫组分更难分辨；（3）在信噪比和T2组分确定时，短的TE能提高T2谱分辨能力.

图8是不同信噪比下BGW方法反演结果，显示回波间隔对T2谱反演的影响，同时证明了短回波间隔能提高T2谱分辨能力，并且信噪比越低时效果越明显.

6.3 回波个数与布点影响

回波个数会影响数据采集的质量，反演的T2谱可能出现不收敛.引起这种不收敛的主要原因是由于所设置的回波间隔较小而采集的回波个数又不够，导致长弛豫组分的信号采集不充分.因此，信号采集时要设置合理的等待时间和足够的回波个数，使长弛豫组分的极化信号充分采集.

T2谱的布点个数，即反演时预先设置的基函数个数，对多指数反演也会有影响.当布点个数较少时，弛豫时间谱的分辨率较低，不能精细的反映储层孔隙结构；当布点个数较多时，可以提高弛豫时间谱对孔隙结构的分辨率，但是布点数越多计算速度就越慢.实验表明，岩心分析时，布点数在32～256个之间较合适[25］.

7 结 论

（1）根据BG反演理论，引入解估计分辨率和解估计方差，每一个解分量会对应一系列折衷参数，其中必定有一个最佳折衷参数使得分辨率和方差达到最佳平衡，此时该点在折衷曲线图中离原点最近.

（2）与以往正则化方法每个解分量都对应相同的正则化参数不同，BGW方法在求解每个解分量时，最佳折衷参数不一定都相同，使求解结果更精确.

（3）通过构造符合不同地层T2谱的模型，并与其它算法相比较，验证了BGW方法在孔隙介质NMR弛豫反演中的有效性.

（4）BGW方法在求解T2谱时，虽然不存在迭代收敛问题，但每个弛豫分量都是依次计算，因此整体计算速度有待优化.

（5）利用离散Picard准则分析了噪声对反演的影响，发现BGW方法在本文大孔占优双峰模型中信噪比低于5后，难以从回波串中提取有效信息.短回波间隔有利于提高T2谱分辨能力，采集回波个数尽量使长弛豫组分采集完全，T2布点需要包含长、短弛豫组分，尽量横跨三个数量级，使得反演结果更能反映粘土束缚水、毛管束缚水和可动流体信息.

（References）

[1]肖立志.核磁共振成像测井与岩石核磁共振及其应用.北京：科学出版社，1998.Xiao L Z.NMR Image Logging and NMR in Rock Experiments（in Chinese）.Beijing：Science Press，1998.

[2]Coates G R，Xiao L Z，Prammer M G.NMR Logging Principles and Applications.Texas：Gulf Publishing Company，1999.

[3]Dunn K J，Bergman D J，Latorraca G A.Nuclear Magnetic Resonance：Petrophysical and Logging Applications.Pergamon：Elsevier Science，2002.

[4]Butler J P，Reeds J A，Dawson S V.Estimating solutions of first kind integral equations with nonnegative constraints and optimal smoothing.SIAMJ Numer.Anal.，1981，18（3）：381-397.

[5]Prammer M G.NMR Pore Size Distributions and Permeability at the Well Site.SPE 28368SPE Annual Technical Conference and Exhibition，25-28September 1994，New Orleans，Louisiana.

[6]Dunn K J，Latorraca G A，Warner J L.On the calculation and interpretation of NMR relaxation time distribution.SPE28367，69th Annual SPE Technical Conference and Exhibition，New Orleans，1994：45-54.

[7]Liaw H K，Kulkarni R，Chen S.Characterization of fluid distributions in porous media by NMR techniques.AICHE J.，1994，42（2）：538-548.

[8]Miller A，Chen S H，Georgi D T，Vozoff K.A new method for estimating T2distributions from NMR measurements.Mag.Res.Imaging，1998，16（5-6）：617-619.

[9]Borgia G C，Brown R J，Fantazzini P.Uniform-penalty inversion of multiexponential decay data.Journal of Magnetic Resonance，1998，132（1）：65-77.

[10]Salazar-Tio R，Sun B Q.Monte Carlo Optimization-inversion Methods for NMR，SPWLA 50th Annual Logging Symposium，June 21-24，2009.

[11]王才志，尚卫忠.应用奇异值分解算法的核磁共振测井解谱方法.石油地球物理勘探，2003，38（1）：91-94.Wang C Z，Shang W Z.Application of singular value decomposition （SVD）in solution of T2relaxation spectra from nuclear magnetic resonance （NMR）log data.Oil Geophysical Prospecting （in Chinese），2003，38（1）：91-94.

[12]姜瑞忠，姚彦平，苗盛等.核磁共振T2谱奇异值分解反演改进算法.石油学报，2005，26（6）：57-59.Jiang R Z，Yao Y P，Miao S，et al.Improved algorithm for singular value decomposition inversion of T2spectrum in nuclear magnetic resonance.Acta Petrolei Sinica（in Chinese），2005，26（6）：57-59.

[13]王忠东，肖立志，刘堂宴.核磁共振弛豫信号多指数反演新方法及其应用.中国科学（G），2003，33（4）：323-332.Wang Z D，Xiao L Z，Liu T Y.New inversion method with multi-exponent and its application.Science in China （G）（in Chinese），2003，33（4）：323-332.

[14]翁爱华，李舟波，莫修文等.核磁共振测井数据高分辨反演方法研究.测井技术，2002，26（6）：455-459.Weng A H，Li Z B，Mo X W，et al.On high resolution inversion of NMR logging data.WLT （in Chinese），2002，26（6）：455-459.

[15]谭茂金，石耀霖，谢关宝.基于遗传算法的核磁共振T2分布反演.测井技术，2007，31（5）：413-416.Tan M J，Shi Y L，Xie G B.NMR T2distribution inversion based on genetic algorithm.Well Logging Technology（in Chinese），2007，31（5）：413-416.

[16]马建海，孙建孟，孙萌等.基于迭代Tikhonov正则化的核磁测井解谱方法研究.西南石油大学报（自然科学版），2009，31（1）：37-40.Ma J H，Sun J M，Sun M，et al.Tikhonov iteration method of echo-stripping T2relaxation spectra from nuclear magnetic resonance log data.Journal of Southwest Petroleum Institute （Science &Technology Edition） （in Chinese），2009，31（1）：37-40.

[17]Song Y Q，Venkataramanan L，Burcaw L.Determining the resolution of Laplace inversion spectrum.The Journal of Chemical Physics，2005，122：104104.

[18]Song Y Q.Resolution and uncertainty of Laplace inversion spectrum.Magnetic Resonance Imaging，2007，25（4）：445-448.

[19]Backus G E，Gilbert J F.The resolving power of gross earth data.Geophys.J.R.Astron.Soc.，1968，16（2）：169-205.

[20]Backus G E，Gilbert J F.Uniqueness in the inversion of inaccurate gross earth data.Philosophical Transactions，1970，266（1173）：123-192.

[21]杨文采.地球物理反演的理论和方法.北京：地质出版社，1996.Yang W C.Theory and Methods of Geophysical Inversion（in Chinese）.Beijing：Geological Publishing House，1996.

[22]Jackson D D.Interpretation of inaccurate，insufficient and inconsistent data.Geophys.，1972，28（2）：97-109.

[23]Hansen P C.Rank-Deficient and Discrete Ill-Posed Problems：Numerical Aspects of Linear Inversion.Philadelphia：SIAM，1998.

[24]Chen S.Improving Estimation of NMR Log T2cutoff Value.SCA，2000.

[25]廖广志，肖立志，谢然红等.孔隙介质核磁共振弛豫测量多指数反演影响因素研究.地球物理学报，2007，50（3）：932-938.Liao G Z，Xiao L Z，Xie R H，et al.Influence factors of multi-exponential inversion of NMR relaxation measurement in porous media.Chinese J.Geophys.（in Chinese），2007，50（3）：932-938.