核磁共振驰豫信号反演问题

关键字：核磁共振 驰豫时间 多指数反演 奇异值分解反演 差分进 化
1、 问题重述
1.1、问题的背景 核磁共振(NMR)在测井技术和岩心分析中已经得到广泛应用，正在
为油气资源的勘探和开发发挥重要作用。油气井中NMR技术提供的油气 藏流体特性和储集层参数，如地层孔隙度、孔径分布、束缚水与可动流 体孔隙体积、渗透率、以及流体的扩散系数和黏度等信息，都需要经过
4)、寻找最优解中最大的负项;
5)、将矩阵A 中对应的第k 列归零,转回步骤 1 重新计算,直到满足非
负约 束。
因此,当最优解中出现负项之后,归零A 中相应的列,使得最后的T2
谱上不出现与这些列对应的值,导致谱线不连续。另外,这个算法的计算 量主要集中在矩阵的奇异值分解上,其效率为。在T2 谱线布点密集时, 要不断重复地进行奇异值分解,其计算量很大。
（二）模型的求解 据此可由MATLAB编程算得矩阵： （全部数据见附录2） 再由常数阵 因此根据： 即可算出 ，结果列表1如下：
0.1 0.36 1.29 4.64 16.68 59.95 215.44 774.26 2782.6 10000
如下：
-1.0255e+009 3.4813e+006 1.7888e+003 -3.5789e+006 5.8481e+007 -2.0196e+008 -2.2378e+009 -1.7246e+009 4.9485e+010 -4.5376e+010
4.4445
5.5556
6.6667
7.7778
由于第一个数据太小，先忽略，即此时 同理去和去噪信号的9组数据，得到系数矩阵（全部数据参见附录2)如 下：
计算出如下表2： 1.0e+006 * 【 0.0641 0.9336 3.8016 1.7125】 做出曲线图：
4.2610
6.4140
2.8537 0.0014 1.7329
（1） 其中为第类孔隙在总孔隙中所占的份额，为第类孔隙的驰豫时间， 通常范围为，为回波间隔时间，为采集时间。 从物理意义上讲,方程的解 应该是非负的。这样, 的全部采样值构成的向 量可以写成如下矩阵形式：
（2） 其中：
（3） （4） （5）来自中是采样点数，是布点数, 是核磁共振总磁化强度信号在时刻的采样
： 含行列的矩阵
5、 模型的建立与求解
5.1、问题一模型的建立与求解 5.1.1对数均匀布点模型 （1） 模型的建立
先对原始信号进行分析，由题目中为第类孔隙的驰豫时间，通常 范围，利用MATLAB软件在内对数均匀选取10个点， 执行命令为： x=logspace(log10(0.1),log10(10000),10)； 得到驰豫时间/s分布点为：
3、 模型假设
1）假设附录一中做给的测量数据每组都是相互独立的。因此可随机选 取其中的数据参与计算。 2）假设核磁共振信号的弛豫时间甚至是连续分布的。 3）假设在信噪比极低时仍能较好的保持驰豫谱分布的真实性。
4、 符号说明
： 第类孔隙的驰豫时间 ： 磁化强度信号
： 第信号采集所用时间 ： 第类孔隙在总孔隙中所占的份额
经计算得：
(2) 改进的SVD 反演算法 可以从向量空间的角度分析非负约束 的求解问题。即将f 看成n
维向量空间的一个坐标点,对方程的求解,实际上是要在n 维空间中的各 坐标分量大于或等于0 的子区域（简称为I区）。寻找满足方程的点。 由于上述方程是严重病态的,通常该点不存在,所以改为寻求n 维空间内 使 最小的一个点。这个点通常有无限个,原因就在于系数矩阵A 的秩及 其运算时的数值秩都比n 要小得多,解的自由度较大。SVD 方法给出的 解是最小意义下的一个最优解, 但不一定满足非负约束条件, 即不一定 位于I 区。传统的SVD 算法通过缩减A 迭代求解, 等于降低了解的维 数,丢掉了振荡最厉害的坐标分量。正是由于缺失了这部分坐标分量,才 导致谱线的不连续。为了获得连续的T2 谱,现改为求解的最小二乘解。
都可以分解为正交矩阵U ,V 和非负对角阵Σ的乘积形式,即：
（6）
式中 r 是A 的秩, T 代表转置。不管A 是否奇异,式(6) 总可以进行
分解, 并且 Σ的形式是惟一的。通过奇异值分解, 式( 2 ) 满足最小
意义下的最优解, 即其最小二乘解可以表示为：
（7）
式中 和分别是U 和V 的第i 个列向量,，是 Σ 中主对角线上的第i
假设已知一个初始解, 令 , 这样原来的方程可以变为： （8）
若求得式（8）在最小意义下的最优解 ,则显然是 最小意义下的最 优解。由于是求解 ,因此在实现非负约束时,就不再缩小矩阵A, 而只须 将小于零的分量改为零,再重新迭代计算 即可,直到所有坐标分量满足 非负约束。由于上述迭代过程没有丢失坐标分量,因此反演出来的T2 谱 是连续的。
（2） 分析问题二，关键是从带有非负约束条件的式中反演出不同的单 弛豫组分所对应的,也就是T2 谱，这样的反演问题测量数据微小 的波动都可能导致解的剧烈变化。本文通过建立多指数驰豫模 型，用奇异值分解反演算法及改进的反演算法进行求解。并通过 计算机编程对两种算法进行检验，改进的反演算法计算稳定，可 满足实际应用需求。
对于问题二岩石核磁共振分析的关键是从带有非负约束条件的式中 反演出不同的单弛豫组分所对应的,也就是T2 谱。这样的反演问题是一 个典型的病态问题,或称非适定问题。测量数据微小的波动都可能导致 解的剧烈变化,尤其是在施加非负约束时,情况更加严重。而奇异值分解 算法是求解这类问题的有效数值解法之一。
对于问题三考虑到测量数据的误差，经过算例会给计算结果带来误 差，例如前面有噪和无噪时造成的结果误差，经查相关文献造成测量数 据误差有等待时间、回波间隔、轻烃或稠油、顺磁物质以及采集方法等 因素。而这些都是在测井的时候就已经存在了的，在这里无法避免，只 能通过改进现有的计算模型以及相关的算法来使得结果更加精确。这里 引进了考虑基线偏离的多指数反演优化模型（即将磁化强度信号描述 为：）和差分进化算法。
核磁共振驰豫信号反演问题
队号：35
队员：张阳（物理院）， 戚向波（物理院）， 杨金新（物理院）
第四轮培训论文
时间：2010年8月23日
核磁共振驰豫信号反演问题
摘要
根据核磁共振理论，单个空隙中的磁化强度信号的衰减满足单指数 衰减规律，但由于岩石内部是一系列大小不等的孔隙群体组成，所以在 岩石核磁共振中测得的中磁化强度信号为一系列单个孔隙磁化强度信号 的叠加，因此可描述为：。其中为第类孔隙在总孔隙中所占的份额，为 第类孔隙的驰豫时间，通常范围为，为回波间隔时间，为采集时间。
一个基本的反演处理，即NMR驰豫信号的多指数拟合，得到驰豫时间的 分布。
在核磁共振测井中, 一般采用CMPG方法测量自旋回波串，纵向驰豫 时间和横向驰豫时间都是用来描述自旋回波信号的驰豫特征。由于谱能 够提供许多岩石物性和流体特性的信息，越来越受到人们的关注。根据 核磁共振理论，单个空隙中的磁化强度信号的衰减满足单指数衰减规 律，但由于岩石内部是一系列大小不等的孔隙群体组成，所以在岩石核 磁共振中测得的中磁化强度信号为一系列单个孔隙磁化强度信号的叠 加，因此可描述为 其中为第类孔隙在总孔隙中所占的份额，为第类孔隙的驰豫时间，通常 范围为，为回波间隔时间，为采集时间。
个元素。
对于带非负约束条件的式（2）的求解, 传统的SVD 算法步骤如
下：
1)、对矩阵A 进行奇异值分解,得到U ,V 和；
2)、由于T2 分布的跨度达到,导致矩阵A 的条件数量太大,所以用信噪
比
（）限制矩阵条件数。忽略Σ中所有小
于 的奇异值,确定矩阵A 的 数值秩;
3)、根据新的,U 和V ,计算最优解和 ;
岩石核磁共振分析的关键是如何从上式中解出各类孔隙的驰豫时间 以及各类孔隙在总孔隙中所占的份额，称之为谱。 1.2、问题的提出
1) 预先给定驰豫时间分布，典型取法为在内对数均匀选取个 点，称之为驰豫时间布点，也可考虑用2的幂指数布点、线性 均匀布点等方式。基于上述几种布点方式建立求解谱的数学 模型
2) 假设不预先给定驰豫时间分布，建立求解谱的数学模型 3) 考虑到测量数据的误差，通过算例分析误差给前面两个问题
（1） 分析题目，建立以线性方程组（其中，）为基础的对数均与布点 模型。用MATLAB软件在内对数均匀选取10个点（），结合附录1 所给数据（），算出系数矩阵，进而得出谱中的分布（表1和图 A）。同理可分别建立2的幂指数布点模型和线性均匀模型（表2 和图B）。根据原始信号和去噪信号得到的不同结果进行对比检 验。
可见此时的均为正值，符合实际。 5.1.2 2的幂指数布点模型 （一）模型的建立与求解
利用MATLAB软件在内选取10个点，执行命令为：
5.1.3线性均匀布点模型 （一）模型建立与求解 利用MATLAB软件在内线性均匀选取10个点，执行命令为：
x=[linespace(0.0001,10,10)] 得到的10个值如下： 0.0001 1.1112 2.2223 3.3334 8.8889 10.0000 趋势如下：
得到 的图像大致
图A （3） 模型分析与改进
由上面得出的结果，可以看出各类孔隙所占总孔隙的份额含有大量 的负值，与实际情况有很大出入。这是因为该信号是原始信号，有许多 的干扰因素。 下面采用前10组去噪信号进行类似的模型求解：
中的没改变，此时
由MATLAB计算得到的值如下： 4.1966e+004 1.5717e+003 7.1043e+003 2.7801e+004 1.9839e+004 3.9737e+004 3.6850e+005 2.7741e+005 3.4944e+006 2.9405e+006 其曲线图如下：