核磁共振弛豫信号的多指数反演
王武蕾王梦佳刘⽂英
2010年8⽉24⽇
摘要
本⽂针对核磁共振信号处理问题，系统地研究了多指数反演的⽅法。
⾸先，本⽂构造了⼀个具有双峰特征的模拟弛豫时间谱，并在不同信噪⽐下，进⾏多指数正演计算，得到反演所需的原始回波数据，正演模拟结果如图2所⽰。 为求解2T时间谱，本⽂从第⼀类积分⽅程的反演求解⼊⼿，在预先知道2T弛豫时T分布的情况下，推导出适合于核磁共振弛豫时间多指数反演的两种算法——奇异间2j 值分解⽅法和阻尼最⼩⼆乘算法。对算法的具体实现过程进⾏了详细的论述，从理论和实例处理两⽅⾯分析讨论了两种算法的优缺点。针对这两种算法，分析了噪声对其解谱的影响，并确定了反演的条件。在进⾏核磁共振信号的多指数反演处理时，应优先
≥时，可选⽤奇异值分解选⽤阻尼最⼩⼆乘算法，只有当原始数据的信噪⽐SNR40
反演算法。
若预先给定弛豫时间分布，由于布点间断、不连续导致所得结果分辨率较低，结果不是很理想。本⽂使⽤基于差分进化算法对核磁共振弛豫信号进⾏多指数反演，在将反演问题转化为带⾮负约束的⾮线性拟合优化问题的基础上，根据⾮负最⼩⼆乘⽅法(NNLS)所确定的弛豫时间的⼤致组分数，直接利⽤差分进化算法进⾏反演进化，从⽽避免了对弛豫时间2T分布的确定。仍利⽤所构造的模拟时间谱，在不同的信噪⽐下使⽤该⽅法进⾏指数反演，从计算精度、抗噪能⼒和计算速度上说，所得结果都⽐上述算法优越，尤其是在信噪⽐较低的情况下，仍保证了弛豫时间谱的真实性。
在对反演算法进⾏研究的基础上，针对提供的实际核磁共振数据，分别使⽤上述三种算法对原始信号和去噪信号进⾏多指数反演，反演结果见图11。根据结果分析，阻尼最⼩⼆乘算法和差分进化算法具有较好的抗噪能⼒和较⾼的计算精度，奇异值分解算法只有在信噪⽐较⾼时才能得到精确的2T时间谱，这与问题⼀、⼆中的分析结果保持⼀致。同时，本⽂使⽤对数均匀分布和2的幂指数分布对弛豫时间布点，⼆者所得结果基本相同，对反演计算没有产⽣重要影响，具体数值见表4。
由于奇异值分解算法只适⽤于⾼信噪⽐数据的反演，在此给出⼀个算法改进的展望，如若可以得到奇异值最佳保留个数和信噪⽐的关系，同样可以得到适⽤于低信噪⽐问题的求解算法。并且对⾮负约束
可构造更快速的迭代算法。基于约束阻尼最⼩⼆乘法算法反演迅速，极为适合多孔介质的核磁共振2T谱反演问题。差分进化算法⽆论信噪⽐的⼤⼩如何，都可有效地求解出核磁共振2T谱反演问题，同时该算法亦可推⼴应⽤于其他⾮指数多项式⾮线性拟合问题的求解中去。
本⽂算法简单有效，思路清晰，对核磁共振岩⼼分析和核磁共振测井⼯作具有⼀定参考价值。⽂中给出的差分进化算法对核磁共振2T谱反演问题的求解虽然有效，但⽬标函数的正则化参数到⽬前为⽌并没有较好的选取⽅法，还有待改进。
关键词：核磁共振弛豫信号多指数反演数值模拟差分进化算法
1. 问题重述
孔隙介质核磁共振(NMR)弛豫信号的多指数反演，在NMR 岩⼼分析和测井技术中起着关键作⽤。⼀般采⽤CMPG ⽅法测量得到反映地下岩⽯孔隙中所含流体的核磁共振信号的原始数据，即⾃旋回波串。⽽岩⽯物性参数以及流体属性等信息才是实际需要的，这些信息需从原始回波信号中提取出来，所以提取这些信息的关键就在于对岩⽯⾃旋回波信号进⾏多弛豫反演。
纵向弛豫时间1T 和横向弛豫时间2T 是描述岩⽯核磁共振信号弛豫特征的两个参数，由于1T 的测量时间很长，在岩⽯核磁共振中⼀般进⾏2T 的测量，故对岩⽯核磁共振信号的多指数反演即指对横向弛豫时间2T 分布谱进⾏反演处理。
岩⽯孔隙是由⼀系列⼤⼩不等的孔隙体组成，单个孔隙的磁化强度信号的衰减满⾜单指数衰减规律，NMR 测得的总弛豫信号()y t 是⼀系列单个孔隙弛豫信号的叠加
()1,j
t
m
T j j y t f e t n τ?==?=
∑ 其中，j f 为第j 类孔隙在总孔隙中所占份额，2j T 为第j 类孔隙的2T 弛豫时间，通常范围为20.110000j ms T ms ≤≤，τ为回波间隔时间，t 为采集时间。岩⽯核磁共振分析的关键就是求解2T 谱，谱的求解本质上为多指数反演问题，即从上式中解出各类孔隙的2T 驰豫时间2j T 以及各类孔隙在总孔隙中所占的份额j f 。
需要求解的问题涉及三个⽅⾯：
1) ⾸先在2T 弛豫时间的分布预先给定的情况下，主要考虑选⽤对数均匀布点、2的幂
尼罗河宝石
指数布点、线性均匀布点等⽅式，建⽴基于上述布点⽅式的、2T 谱求解模型； 2) 然后在不预先给定2T 驰豫时间分布2j T 的情况下，建⽴求解2T 谱的数学模型； 3) 最后考虑测量数据的误差即随机噪声，利⽤题⽬所给数据，进⾏实例计算，分析误
差对模型求解的影响；并给出如何消除或减⼩误差带来的不良影响的具体⽅法。
2. 问题分析
岩⽯核磁共振分析的关键是如何得到各类孔隙的2T 弛豫时间2j T 以及各类孔隙在总孔隙中所占的份额j f ，称之为2T 谱。核磁共振弛豫信号的解谱本⾝是⼀个⾼度病态的反演问题，各种解谱⽅法在应⽤过程中都有各⾃的使⽤条件[1]。
求解2T 谱时要考虑两种情况，⼀种即为2T 弛豫时间2j T 的分布预先给定，另⼀种就是不预先给出2T 弛豫时间2j T 的分布。
2T 弛豫时间2j T 分布的构造⽅法可以是在()2min 2max ,T T 内进⾏对数均匀布点或是⽤
2的幂指数布点和线性均匀布点等⽅式进⾏布点。进⾏反演计算，需要⼀定量的原始回波信号数据，由于附表所给数据有限，因此本⽂构造⼀个具有双峰特征的2T 谱分布，构造的2T 谱分为10类。考虑在⽆噪声条件和有噪声条件下利⽤奇异值分解算法、约束最⼩⼆乘法和约束阻尼最⼩⼆乘法进⾏反演
计算。对⽐各种算法所得结果，⽐较计算速度，选取最佳的⽅法进⾏核磁弛豫信号反演。
因2T 跨度⽐较⼤，导致布点数⾮常多，布点不连续，反演所得2T 谱的分辨率就⽐较低。当2T 弛豫时间2j T 分布未知时，采⽤⾮负最⼩⼆乘法和⾮线性拟合相结合的思想，将2T 布点当做未知数，采⽤⾮线性拟合求解。由于该⽅法对初值⽐较敏感，可能陷⼊问题的局部极值。⽽差分进化算法则⽆需设置初值，且原理简单，易于理解和实现，故考虑运⽤差分进化算法求解。
数据测量时会产⽣⼀定的误差，利⽤上述给出的反演算法，对所得结果进⾏分析，考察随机噪声在不同算法下对反演结果的影响。发现信噪⽐的⼤⼩⾄关重要，当信噪⽐所在范围不同时，选⽤的算法也要随之改变。为了克服误差所带来的影响，针对不同的信噪⽐，选择不同的算法，可以避免误差影响。
3. 模型假设与符号说明
3.1. 模型假设
1. 以构造谱正演所得到的回波信号，作为进⾏多指数反演的回波信号数据；
2. 正演所得的回波信号为⽆噪声的回波信号，即SNR =+∞；
3. 以在正演所得的回波信号中逐点加⼊⼀定⼤⼩的信噪⽐的⽅式，
构造有噪声的回波信号； 4. 测量数据误差—随机噪声为⾼斯随机⽩噪声，服从正态分布； 5. 不考虑回波间隔时间的⼤⼩对2T 谱反演精度的影响。
3.2. 符号说明
2j T 第j 类空隙的2T 弛豫时间
j f 第j 类孔隙在总孔隙中所占份额 SNR 从测量数据中估计出的信噪⽐
()X t m 维时间衰减向量
α正则化参数
D 问题空间维数λ加权因⼦
()t ε测量误差或随机噪声
其他变量符号在⽂中使⽤时给出说明。
4. 问题⼀：确定分布下弛豫时间2T 谱的求解
训练日剧情
核磁共振(NMR)在测井技术和岩⼼分析中已得到⼴泛应⽤，在岩⼼分析中，NMR 技术提供的岩⽯物性参数以及流体属性等信息，都需要经过⼀个基本的反演处理，即NMR 弛豫信号的多指数反演，得到弛豫时间的分布[2]。
孔隙介质核磁共振(NMR)弛豫信号的多指数反演在NMR 岩⼼分析与测井解释中起着关键作⽤。⼀般采⽤CMPG ⽅法，经过测量得到反映地下岩⽯孔隙中所含流体的核磁共振信号的原始数据，即⾃旋回波信号。⾃旋回波信号的弛豫特征由纵向弛豫时间1T 和横向弛豫时间2T 进⾏描述。⽽为了得到岩⽯的物性参数以及流体属性等信息，需从原始回波信号中提取，所以提取这些信息的关键就在于对岩⽯⾃旋回波信号进⾏多弛豫反演。多弛豫反演的关键是解出各类孔隙的2T 驰豫时间2j T 以及各类孔隙在总孔隙中所占的份额j f ，称之为2T 谱[3]。
4.1. 岩⽯核磁共振信号的多指数特征
纵向弛豫时间1T 和横向弛豫时间2T 是描述岩⽯核磁共振信号弛豫特征的两个参数，由于1T 的测量时间很长，在岩⽯核磁共振中⼀般进⾏2T 的测量，故对岩⽯核磁共振信号的多指数反演即指对横向弛豫时间2T 分布谱进⾏反演处理。
岩⽯孔隙是由⼀系列⼤⼩不等的空隙体组成，单个孔隙的磁化强度信号的衰减满⾜单指数衰减规律。
()
2
0t T s t S e ?=
(1)
⽽NMR 测得的总弛豫信号()y t 是⼀系列单个孔隙弛豫信号的叠加
()1
,j
t
m
T j j y t f e t n τ?==?=
∑ (2)
其中，j f 为第j 类孔隙在总孔隙中所占份额，2j T 为第j 类空隙的2T 弛豫时间，通常范围为20.110000j ms T ms ≤≤，τ为回波间隔时间，t 为采集时间。
求解上式实际上是解第⼀类fredholm 积分⽅程，下⾯给出求解上述⽅程的具体模型[4] 。
4.2. 模型数据的准备
为进⾏反演计算，需要得到⼀定量的原始回波信号数据，由于附表所给数据量有限。为此，先构造⼀个具有双峰特性的2T 谱分布(图1)，构造的2T 谱分为10类，即10m =，
具体分布见附录表 3。
通过对构造谱进⾏正演计算，基于式(2)可构造出⼀个⽆噪声的回波衰减信号()y t 。利⽤该衰减信号，从中选取1024个信号强度数据，作为模型中反演计算的原始回波数据。在上述⽆噪声的回波衰减信号中，加⼊信噪⽐SNR=80,40,20的噪声后，相应的信号强度分布如图2所⽰。
4.3. 基于奇异值分解(SVD)的反演模型
4.3.1. 奇异值分解算法
奇异值分解(SVD)算法可以⽤来求解⼤多数的线性最⼩⼆乘法问题，SVD 算法基于如下分解定理：对任意的矩阵m n A ×，当其⾏数m ⼤于等于列数n 时，可以分解为正交矩阵n n W ×以及正交矩阵n n V ×的转置的乘积，即：
()T m n m n j n n n n
A U diag W V ×××× =?? (3)
其中0j W ≥，U V ，为正交矩阵，即满⾜：
1
m
ij ik
jk i u u
δ==∑，1
m
锵锵三人行2008
ij ik jk i ννδ==∑
当m n <，SVD 也可以执⾏，在这种情况下，奇异值0j W =，并且U 中相应列的值半角字符
都为零，这时上式仅在,j k m ≤时成⽴。故不管矩阵A 是否是奇异，总可以对其进⾏分解，⽽且这个分解是唯⼀的。SVD 分解明确地构造了矩阵零空间和值域的正交标准化基，特别地，对U 的列，若与其标号相同的元素j W 为零，则其列为值域的⼀个正交标准化的基础⽮量；对于V 列，若与其标号相同
的元素j W 为零，则其列为零空间的⼀个正交标准化基。
4.3.2. 基于SVD 的反演模型的建⽴
对于多指数衰减2T 模型
图1构造的弛豫时间分布
图 2 模拟的回波信号串
亚星奔驰y M f =?
其中，()
12,,,T
n y y y y = 为正演计算得到的⾃旋回波衰减信号，
2/i j t T ij n m
n m
M
课堂内外初中版m e ?×× =
= ，()12,,,T m f f f f = 为弛豫时间2j T 对应的各点幅度值。
4.3.2.1.
对数均匀分布条件下的M 矩阵的建⽴给定2T 弛豫时间分布2j T 为对数均匀分布，即在()2min 2max ,T T 内对数均匀选取m 个点，满⾜
10211021021021
log log 1log log j j j j T T T T +??=?。由此给出M 矩阵为：
11121222m
222221
22
2122
2m T T T T T T m n n n t t t t t t t t t e e
e e e e M e
e
e =
4.3.2.2. 2的幂指数分布条件下的M 矩阵的建⽴
给定2T 弛豫时间分布2j T 为2的幂指数分布，即在()2min 2max ,T T 内对数均匀选取m 个点，满⾜22,1,2,j j
T j m ==
，由此可得到相应的M 矩阵。
4.3.3. 基于SVD 的反演模型的求解
4.3.3.1.
⽆噪声条件下模型的求解采⽤SVD 分解法，系数矩阵()T
n m n m j m m m m
M U diag W V ×××× =?? ，这⾥U V ，为正交矩阵，()j diag W 为对⾓矩阵，其对⾓元呈递减排列，则很容易地求得最⼩⼆乘意义下的解
()12111
,,,,0,,0T n f V diag U y w w w ∧
=
(4)
运⽤MATLAB 编程进⾏反演计算，求解得到2T 谱分布。 4.3.3.2.
具有⼀定信噪⽐条件下模型的求解考虑到正演计算所得数据是⽆噪声的，⽽在实际测量过程中，由于受到电⼦部件
和环境的影响，不可避免地在测量过程中产⽣随机误差。为此在进⾏SVD 分解时，加⼊随机噪声这⼀因素，得到了在噪声条件下的最⼩⼆乘解
()1211,,,,0,,0T n SNR
f V dia
g U y w w w ∧
=
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议