岩石磁共振多指数反演算法研究

岩石磁共振多指数反演算法研究

磁共振石油岩心分析和钻探技术可以提供丰富的岩石和流量特征参数，例如岩石中的空隙大小分布；有效孔隙度、渗透率、可动流和不可动流之间的关系等。通过差谱法和移谱法，可以定性或定量地提供地层含油气的指示。通过测量扩散系数，我们可以获得各种重要信息，如岩层的润湿度和流粘度。这些独特的岩石和材料的物理信息为地质勘探和油田开发工程师以及制定储层评价、产量预测和开发方案提供了重要依据。磁共振技术在油田勘探和开发中的应用引起了国内外的高度重视。核磁共振石油应用技术是一个正在发展的新技术,本身还不完善,近年来,我国花费了巨额外汇引进多台美国NUMAR公司核磁共振成像测井仪,但应用表明:国外核磁测井解释软件由于反演指数少、弛豫谱精度差、无法反映裂缝及溶洞等复杂的地质情况,应用效果不好.要解决这个问题,必须研究更精确的多弛豫反演算法,并形成适合于我国复杂陆相储层地质条件的岩石核磁测量数据的解释标准、解释模型.在自然科学中,很多参数都是采用间接方法获得的,也就是说,与观测数据相关的预测函数或预测矢量是由某种算符确定的,这种算符一般可采用线性积分算符来表示.而在观测数据中通常带有未知的噪声,问题关键是怎样从含有噪声的数据中估计或预测出我们所期望获得的函数或矢量.现代核磁共振石油工业应用测试技术也是一种间接测量技术,它所采集到的原始数据是地下岩石孔隙中所含流体的核磁共振信号,即所谓的自旋回波信号.我们所需的岩石物性参数以及流体属性等信息是从原始回波信号中提取出来的,提取这些信息的关键在于对岩石自旋回波信号进行多弛豫反演,这个反演实际上是求解第一类Fredholm线性积分方程,这种积分方程通常的解不是惟一的.在对算符求逆时,积分核作用到函数上,起着光滑的作用,输入原始数据中的小变化(如核磁共振测量过程中的噪声等)对解的结果有较大的影响.光滑的结果常常是丢失了信息,而逆运算又无法找回它,因此,我们必须利用解的一些先验知识(如解的连续性、数据的信噪比、解的非负限制等)来最大限度地恢复我们所需的信息.本文从第一类积分方程的反演求解入手,推导出适合于岩石核磁弛豫多指数反演的两种算法:奇异值分解(singularvaluedecomposition,SVD)方法和变换反演算法,分析了两种算法优缺点,给出了最优反演模型的选取原则.1孔隙磁化强度信号有两个参数可用来描述核(在岩石核磁共振中指氢核)自旋信号的弛豫特征:纵向弛豫时间T1和横向弛豫时间T2.由于T1的测量时间很长,所以在岩石核磁共振中一般进行T2测量.根据核磁共振理论分析,单个孔隙中的磁化强度信号的衰减满足单指数衰减规律,但由于岩石内部是由一系列大小不等的孔隙群体组成的,所以在岩石核磁共振中测得的总磁化强度信号y(t)是一系列单个孔隙磁化强度信号的叠加,同时在实际测量过程中,不可避免地产生随机噪声,因此y(t)可描述为y(t)=∑ifi⋅exp(-tΤ2i)+ε(t),t=n⋅τ,(1)y(t)=∑ifi⋅exp(−tT2i)+ε(t),t=n⋅τ,(1)其中fi为第i类孔隙在总孔隙中所占的份额,T2i为第i类孔隙的T2弛豫时间,τ为回波间隔时间,ε(t)为随机噪声序列.岩石核磁共振分析的关键是如何从(1)式中求解出各类孔隙的T2弛豫时间T2i以及各类孔隙在总孔隙中所占的份额fi,我们有时也称之为T2谱.2不同f2弛豫时间分布函数的估计结果求解方程(1),实际上是解第一类Fredholm积分方程,这是一个非适定问题,也就是说,在允许的测量误差ε条件下,存在不同的f(T2)弛豫时间分布函数都能相当好地拟合原始回波衰减曲线,这种相容性增大了反演过程的难度.2.1svd分解的明确构造SVD算法可用来求解大多数的线性最小二乘法问题.SVD算法基于如下分解定理:对任意的矩阵Am×n,当其行数m大于等于列数n时,可以分解为正交矩阵Um×n,非负对角矩阵Wn×n以及正交矩阵Vn×n的转置的乘积,即Am×n=Um×n⋅[diag(wj)]n×n⋅VΤn×n,(2)Am×n=Um×n⋅[diag(wj)]n×n⋅VTn×n,(2)其中wj≥0(1≤j≤n);U,V为正交矩阵,即满足m∑i=1uijuik=δjk‚n∑i=1vijvik=δjk(1≤j,k≤n),(3)∑i=1muijuik=δjk‚∑i=1nvijvik=δjk(1≤j,k≤n),(3)当m<n时,SVD也可以执行,在这种情况下,奇异值wj=0(m+1≤j≤n),并且U中相应的列都是零,这时(3)式仅对j,k≤m时成立.故不管矩阵A是否是奇异,(2)式的分解总可以进行,而且这个分解几乎是惟一的.也就是说,其分解形式惟一到:对矩阵U的列、W的元素和V的列能做相同的置换,或者矩阵U和V的任意列的线性组合,在W中对应的元素仍恰好完全相同.SVD分解明确地构造了矩阵零空间和值域的正交标准化基.特别地,对U的列,若与其标号相同的元素wj为非零元,则其列为值域的一个正交标准化的基础矢量;对V的列,若与其标号相同的wj为零,则其列为零空间的一个正交标准化基.对于如下的多指数衰减T2模型,有y=Μ⋅f,(4)y=M⋅f,(4)其中y=(y1,y2,…,yn)T为测量的自旋回波衰减信号,M=[mij]n×m=[exp(-ti/T2j)]n×m;f=(f1,f2,…,fm)T为弛豫时间T2j对应的各点的幅度值,T2j(j=1,2,…,m)为预先指定的T2时间分布系列,典型的取法为在(T2min,T2max)区间内对数均匀地选取m个点,我们称为弛豫时间布点,也可采用2的幂指数布点、线性均匀布点等方式.矩阵条件数的定义为矩阵的最大特征值与最小特征值的比值.若矩阵的条件数为无穷大,则该矩阵奇异;若矩阵的条件数太大,即其倒数超出了机器的浮点精度,则称该矩阵为病态的矩阵.若直接采用Gauss分解求上式,几乎是不可能的,原因是矩阵M的条件数相当大,例如:若回波间隔τ=1.2×10-3s,T2在0.1×10-3~10000×10-3s内对数均匀地取50个点,则矩阵M的条件数可达1016数量级,很明显,矩阵M是高度病态的.采用SVD分解法来求解上式,系数矩阵Mn×m=Un×m·[diag(wj)]m×m·VTm×m,这里U,V为正交矩阵,diag(wj)为对角矩阵,其对角元递减排列,则我们就可以很容易地求得最小二乘意义下的解为ˆf=V⋅[diag(1w1,1w2,⋯‚SΝRw1,0,⋯‚0)]⋅(UΤ⋅y).(5)fˆ=V⋅[diag(1w1,1w2,⋯‚SNRw1,0,⋯‚0)]⋅(UT⋅y).(5)这里给出了矩阵条件数小于等于SNR的限制,避免了解的不稳定性.其中SNR为从测量数据中估计出的信噪比.SNR定义为第1个回波的幅度值除以误差矢量r(r=y-Μ⋅ˆf)的标准差σ.2.2im1,,首先给定如下的目标函数:χ2=n∑i=1[yi-m∑j=1(fj⋅mij)]2+λ⋅m∑j=1f2j=∥y-Μf∥2+λ∥f∥2,(6)这里M=[mij]=[exp(-ti/T2j)],λ为平滑因子.对幅度f=(f1,f2,…,fm)T的第k分量求极值并令其等于0,则有∂χ2∂fk=-2n∑i=1(yi-m∑j=1fj⋅mij)⋅mik+2λ⋅fk=0.(7)交换求和顺序,并移项整理,可得m∑j=1[fj⋅n∑i=1(mik⋅mij)]+λfk=n∑i=1(mik⋅yi).(8)很容易验证,k=1,2,…,m的m个等式组成的方程组满足(采用矩阵形式给出)(ΜΤΜ)⋅f+λΙm×m⋅f=ΜΤ⋅y,(9)上式中Im×m为m×m单位矩阵.我们对方程(9)作如下线性变换,令f=ΜΤ⋅c,(10)注意(10)式中未知变量c=(c1,c2,…,cn)T为n×1维的,而不是m×1维的.采用上述变换我们将m维T2域空间的解变换到n维时域空间来求解.将(10)式代入(9)式,则有ΜΤ⋅(Μ⋅ΜΤ+λΙn×n)⋅c=ΜΤ⋅y,(11)则原问题的解f可以通过求解出方程(Μ⋅ΜΤ+λΙn×n)⋅c=y(12)的解c,再通过线性变换(10)式回代而获得.选择合适的λ,以保证矩阵(M·MT+λIn×n)的可逆,则我们就可以很容易地求得方程(1)的最小二乘解ˆf=ΜΤ⋅(Μ⋅ΜΤ+λΙn×n)-1⋅y.(13)平滑因子λ的选取:增大λ会使解的残差r=∥y-Μˆf∥远离最小的情形,趋向于使‖ˆf‖最小;反之,减小λ又会遇到矩阵(Μ⋅ΜΤ+λΙn×n)病态的问题,使我们获得的解变得不可靠.λ的选取应从主观和客观两个方面加以判断,客观上要求λ尽可能地小,主观上要求解尽可能地光滑连续.根据先验信息,若我们已知测量误差σ,则比较理想的平滑因子为λ=√n⋅σ∥c∥,(14)还可以选择不同平滑项,以实现不同的平滑方式.2.3svd算法的改进若对上述两种算法不加任何限制条件,则求解出的T2谱的幅度ˆf中很可能存在小于零的分量,这在物理意义上来说是不合逻辑的.故在求解模型中,还必须加上非负性限制条件,即fj≥0(j=1,2,…,m).两种算法实现非负性限制的方法是不一样的.SVD算法采用迭代法逐渐消去负的组分,即在我们求出的ˆf中,找出具有最大负数的组分,假定为第k项,说明第k项最不合理,应该在原始模型中去除该组分,令与fk有关的所有系数为零,消去第k项,然后重新进行反演计算,直到满足所有的fj≥0(j=1,2,…,m)为止.这种方法的缺点是会破坏T2分布的连续性,造成T2谱的畸形.变换反演算法也是采用迭代法逐渐消去负数的组分,但该方法通过T2域到时域的变换后,将第k项相邻点的T2分布信息叠加到负数的T2k组分中,重新进行反演计算,直到满足所有的fj≥0(j=1,2,…,m)为止.这种方法的优点是可以保持T2分布的连续性,使反演结果更符合实际.3t2反演实验变换反演算法采用变换的方式,将T2域空间的解变换到时域空间进行求解,优点是算法在实现非负性限制条件时,可利用相邻点的T2分布信息来保持T2分布的连续性,并且其设计矩阵(MMT+λI)的大小固定不变,总是n×n的,算法稳定且容易实现.SVD算法是对设计矩阵直接进行SVD分解,采用原始数据的信噪比SNR作为矩阵条件数的限制,且SNR可以在拟合过程中动态的获取,避免了变换反演算法中提到的平滑参数λ的选取问题.但SVD算法在非负限制条件的实现过程中,有可能造成T2谱的畸形,特别是当原始数据信噪比较低时,反演结果的分辨率比较低,解会出现不规则的跳动,其原因在于低信噪比的数据会使解的自由度增大.我们可用如下的实例来说明这个问题.假设核磁测量的回波间隔τ=1.2×10-3s,T2弛豫时间在0.1×10-3~10000×10-3s内对数均匀地选取50个点,采用SVD算法进行反演时,在不同信噪比下所获得的解的自由度见表1.解的自由度越大,解中所包含的线性无关的信息含量就越少,也就是说,解的分辨率越低,我们所获得的信息量就越小.同样,解的自由度太小,也会造成严重的误差,原因是解对数据的适应性差,表明我们所选的模型不合适,如T2布点数太少等.目前的核磁测井数据的信噪比一般为30左右,采用SVD算法反演出的T2谱中所含有的线性无关的信息量不超过4,T2谱的分辨率比较低.这就是美国Numar核磁解释处理软件所反演出的T2谱大多呈现一个峰的缘故.同时,SVD算法计算处理相当费时,尤其是T2布点数较多时就更为明显.我们的计算结果表明当原始数据的信噪比大于80时,SVD算法可以给出好的结果.为了便于对比分析,我们构造了一个具有双峰结构的典型石油岩心T2谱(图1),假定回波间隔时间τ=1.2×10-3s,回波数NE=1024,可构造出无噪声的回波衰减信号.在4种不同T2对数均匀布点数(100,40,15和8)条件下,用SVD法进行反演计算,图1给出了反演结果的对比.从图中可以看出,T2布点数设为100和40时所得的结果差异较小,而T2布点数为15和8的谱分辨率明显变差.因此为保证T2分布的真实性,在反演过程中应保持一定的T2点数和覆盖范围.我们的研究结果表明,T2布点数为30~50较为合适,布点区间应跨越3个数量级.为了对两种解法进行直观对比分析,在上述构造的无噪声回波衰减信号中,逐点加上随机Gauss白噪声,构造出不同信噪比的回波串.分别对SNR=∞,100,50,30,20,10的6个回波串(其中SNR=∞为无噪声的回波串)进行反演,反演中均采用对数均匀布点,布点区间为0.1×10-3~10000×10-3s,布点数为50.通过反演共获得12个不同的T2谱.图2给出了SVD反演的6个T2谱的对比.信噪比≥50的3条T2谱线具有双峰结构,与真实T2谱线非常接近;SNR为30,20的2条T2谱线只具有单峰结构,其单峰仍能较好地反映出真实谱的大峰结构;而SNR为10的T2谱线与真实谱线相差较远,且存在不连续的现象.从表1中可以发现上述几种情况的解所包含的线性无关信息量是不同的,这就是它们的谱与真实谱差异的原因.图3给出了变换反演算法反演的结果对比.SNR≥20的T2谱线均呈现双峰态结构,与真实谱吻合较好,而SNR为10的T2谱线偏离真实谱线较远.4精细解释及其应用作为一个应用实例,这里介绍本文变换反演算法在青海柴达木盆地狮子沟构造狮32井中的应用情况,狮子沟储层岩石孔隙结构比较复杂,既存在基质孔,也存在裂缝和微裂缝,还存在少量溶蚀孔,其T2谱大部分呈三峰态.青海油田采用美国NUMAR公司生产的C型核磁共振成像测井仪对该井进行了测井,并采用该公司的解释软件对核磁回波串进行了处理,NUMAR公司所采用的是SVD算法,T2布点方式为2的幂次方,布点数只有8个,因而很难反映复杂孔隙结构的细节,解释结果无法用于石油储层评价和储量计算.针对上述问题,我们以本文变换反演算法为基础,编写了一套核磁测井数据精细解释软件,对该井的核磁测