c题 核磁共振驰豫信号反演问题
承 诺 书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 中南大学 参赛队员 (打印并签名) ：1. 杨金风 指导教师或指导教师组负责人 (打印并签名) ：
日期： 2010年 08月 23 日
赛区评阅编号（由赛区组委会评阅前进行编号）：
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
核磁共振驰豫信号反演问题
摘要
核磁共振在测井技术和岩心分析中已经得到广泛应用。本文主要讨论核磁共振驰豫信号反演问题。针对横向驰豫时间T 2的不同布点方式，建立数学模型，通过NNLS 算法和差分进化算法进行T 2谱的多指数反演，并对结果进行相关讨论。
考虑到原始数据含有噪声，在前两问中均采用去噪信号进行求解，从而使所得T 2谱更符合要求；第三问中，给定一个双峰T 2谱，并利用matlab7.0产生随机噪声生成模拟信号，通过对比反演结果，总结测量误差对前两问算法结果的影响。
在预先给定弛豫信号分布的情况下，在不同布点方式中，根据T 2的时间分布区间，分别合理确定m 值进行取点，此时T 2j 均确定。求解f j 的过程转化为曲线拟合问题，且需保证f j >0，建立以曲线拟合残差为目标函数的非负优化模型，采用非负最小二乘（NNLS ）算法，利用matlab7.0求解出T 2谱和其反演曲线与原曲线对比图（见图4.1至4.11）。由各T 2谱的数据表（见表4.1，表4.2，表4.3）及对比图，得出如下结论：利用NNLS 算法进行多指数反演时，最佳的布点方式是2的幂指数均匀布点，次之是线性均匀布点方式，而且随着布点数增多，二者反演曲线与原曲线的拟合度提高；该算法不太适用于对数均匀布点方式，拟合曲线的残差较大。
上述T 2谱反演方法中对f j 非负约束的处理，采用非负迭代的方法。由于T 2分布事先给定，如果所测量的T 2组分离散且分布较宽时，反演得到的T 2谱分辨率会较低。于是考虑未预先给定T 2j 分布的情况下，以T 2j 和f j 为待求量，利用差分优化算法，求解此带非负约束的优化问题，同时解出T 2j 和f j ，即T 2谱。此处令m =5,10得出两组T 2谱。
考虑到原题中测量误差对结果的影响，本文采用计算机产生模拟信号来分析误差对于前两个问题的影响。首先给定一个具有双峰特性的T 2谱分布y ，由该频谱可计算不同t i 对应得信号强度y i 。由matlab7.0随机生成四组噪声，与y 叠加得到四组模拟线号y ，再由公式（2）计算出各模拟信号信噪比SNR 。对于四组模拟信号，在预先给定驰豫时间分布T 2的情况下，采用NNLS 算法和差分进化算法反演出T 2频谱。将反演结果与原信号y 对比，得出如下结论：NNLS 算法和差分进化算法在SNR 较高的情况下计算结果均比较准确。在SNR 较低时，NNLS 算法反演结果与实际信号差距较大，而差分进化算法精确度仍然比较高。
根据上述分析可知，首先对原始信号进行去噪处理，若信号SNR 较高采用NNLS 算法或者差分进化算法，若信号SNR 较低则采用差分进化算法进行驰豫信号反演，从而减小误差带来的不良影响。
关键字：核磁共振 多指数反演 NNLS 算法 差分进化算法 SNR
~
一、问题重述
1.1 问题背景
核磁共振(NMR)在测井技术和岩心分析中已经得到广泛应用，正在为油气资源的勘探和开发发挥重要作用。油气井中NMR 技术提供的油气藏流体特性和储集层参数，如地层孔隙度、孔径分布、束缚水与可动流体孔隙体积、渗透率、以及流体的扩散系数和黏度等信息，都需要经过一个基本的反演处理，即NMR 驰豫信号的多指数拟合，得到驰豫时间的分布。
在核磁共振测井中, 一般采用CMPG 方法测量自旋回波串，纵向驰豫时间T 1
和横向驰豫时间T 2都是用来描述自旋回波信号的驰豫特征。由于T 2谱能够提供许多岩石物性和流体特性的信息，越来越受到人们的关注。 1.1 问题提出
(1)预先给定T 2驰豫时间分布T 2j ，典型取法为在(T 2m in , T 2m ax ) 内对数均匀选取
m 个点，称之为驰豫时间布点，也可考虑用2的幂指数布点、线性均匀布点等方
式。基于上述几种布点方式快慢建立求解T 2谱的数学模型
(2)假设不预先给定T 2驰豫时间分布T 2j ，建立求解T 2谱的数学模型
(3)考虑到测量数据的误差，通过算例分析误差给前面两个问题带来的影响；并考虑如何克服误差带来的不良影响
二、模型假设
（1） 假设本模型中的自旋回波串是无噪声(SNR=∞) 理想回波串 （2） 假设研究对象是匀质样品，NMR 信号通常表现为单指数衰减 （3） 假设每个点的测量误差是随机独立的，且标准差都一样
三、符号说明
f j T 2j
： 第j 类孔隙在总孔隙中所占的份额 ： 第j 类孔隙的T 2驰豫时间
： 回波间隔时间
t
： 采集时间
四、模型的建立与求解
4.1 问题一：已知弛豫时间T 2分布求解T 2谱 4.1.1 问题分析
在岩心核磁共振实验分析与核磁共振测井解释分析中，关键是由方程
m
-t T 2
j
y (t ) =
∑
j =1
f j ⋅e
, t =n ⋅τ
求解样品各个流体单元(岩石孔隙) 的T 2j 弛豫时间和T 2j
及其在总流体(岩石总孔隙) 中的相应份额f j ，即求解出样品的T 2谱。
横向弛豫时间T 2的测量一般采用自旋回波，或者多回波(CPMG)脉冲序列。后者由于测量速度快而被更广泛地采用．对于单个弛豫时间的数据分析较简单，采用线性拟合算法就可以得到较精确的结果。当样品中包含多个离散的弛豫时间时，一般采用非线性拟合算法来分析处理。而在样品的弛豫时间呈连续分布时，较有效的处理方法是多指数反演算法。多指数反演算法在处理多个离散(≥4) 甚至连续分布的弛豫数据时较有效，但是其计算量相对较大，因此基(横坐标的点) 的选取不能过多。
实际上，T 2j 各组分的幅度不会出现负值，这就要求保证反演结果的非负性。根据上述分析，采用非负最小二乘算法。 4.1.2 模型建立
原始时域信号y (t ) 可描述为
m
-t T 2
j
y (t ) =
∑
j =1
f j ⋅e
, t =n ⋅τ
（1）
其中f j 为第j 类孔隙在总孔隙中所占的份额，T 2j 为第j 类孔隙的T 2驰豫时间，通常范围为0.1m s ≤T 2j ≤10000m s ，τ为回波间隔时间，t 为采集时间。、 定义
-t i T 2
j
A ij =e
假设每个点的测量误差是随机独立的，且标准差都一样。问题一中已假设每个点的测量误差是随机独立的，且标准差都一样。问题一中预先给定T 2驰豫时间分布T 2j ，且0.1m s ≤T 2j ≤10000m s ，因此当m 一定时，采用对数均匀布点、2的幂指数布点和线性均匀布点任意一种布点方式，T 2j 均为固定值。求解f j 过程
实质上就变成了一个最优曲线拟合过程。曲线拟合时残差越小，表示拟合度越好，因此以曲线拟合的残差为目标函数。由于T 2频谱的幅度必正数因此得到约束条件
f
j
≥0。最终建立具有非负约束的优化模型：
n
min R =
∑
i =1
⎡
⎢y i -⎣
m
∑
j =1
⎤f j A ij ⎥
⎦
2
s . t .
f j ≥0
j =1, 2, m
其中R 表示曲线拟合残差，y i 表示实测信号采样强度，f i 为第j 类孔隙在总孔隙中所占的份额。 4.1.3 模型求解及结果
下面，通过不同的布点方式，将上述非负最小二乘算法进行具体步骤分析。具体步骤流程如下：
Step 1：对数布点方式
（1）对数均匀取点
在对数取点方式中，由于T 2=10000ms 过大，在取点计算时无法得出合适的数据，因此在该取点方式取中T 2的单位为秒。为了研究在对数布点方式中不同取点数对结果的影响，此处分别取m =5,m =15,m =50。
（2）计算不同取点时f j
利用matlab7.0编程，在m 取不同值的情况下，通过非负最小二乘拟合曲线得出对应的f j 。当m =5,15,50时，由于曲线拟合结果始终相同，在此均用表4.1表示结果：
表4.1 m=a(a=5,15,50)时对数取点拟合结果f j 及残差表
结论：对数均匀取点时，时间t 单位为s ，当m =5, m =15, m =50时，T 2谱上始终只有最后一个分量，即T 2谱始终相同，y (t ) =107440e -t /10。
（3）反演曲线与去噪信号对比图和T 2谱
图4.1 m=5时对数均匀布点反演曲线图和T 2谱
图4.2 m=15时对数均匀布点反演曲线图和T 2谱
图4.3 m=50时对数均匀布点反演曲线图和T 2谱
结论：在所给时间样本的约束下，在对数布点方式中，不同的m 值对最终的结果并没有影响。比较实测曲线和反演曲线可知，曲线拟合存在一定的误差。由MATLAB7.0计算出残差为2.77E+09。 Step 2：2的幂指数布点方式 （1）2的幂指数均匀取点
在该取点方式取中T 2的单位为毫秒。为了研究在2的幂指数布点方式中不同取点数对结果的影响，此处分别取m =5,m =15,m =20,m =50。
（2）计算不同取点时f j
利用matlab7.0编程，在m 取不同值的情况下，通过非负最小二乘拟合曲线得出对应的f j 。当m =5，15，20，50时，拟合结果f j 如表4.2所示：
表4.2 2的幂指数取点拟合结果f j 及残差表
结论：2的幂指数均匀取点时，时间t 单位为m s ，由上表可知，随着布点点数的增多，反演曲线与去噪信号曲线的相似度增大。同时，残差r 从，依次减小。说明随着布点数目增多，拟合程度越来越高。
（3）反演曲线与去噪信号对比图和T 2谱
图4.4 m=5时2的幂指数均匀布点反演曲线图和T 2谱
图4.5 m=15时2的幂指数均匀布点反演曲线图和T 2谱
图4.6 m=20时2的幂指数均匀布点反演曲线图和T 2谱
图4.7 m=50时2的幂指数均匀布点反演曲线图和T 2谱
Step 3：线性布点方式
（1）线性均匀取点
在该取点方式取中T 2的单位为毫秒。为了研究在2的幂指数布点方式中不同取点数对结果的影响，此处分别取m =100,m =300,m =600,m =1000。
（2）计算不同取点时f j
利用matlab7.0编程，在m 取不同值的情况下，通过非负最小二乘拟合曲线得出对应的f j 。当m =100，300，600，1000时，拟合结果f j 如表4.3所示：
表4.3 线性均匀取点拟合结果f j 及残差表
结论：线性均匀取点时，时间t 单位为m s ，随着布点数目的增多，反演曲线
与原曲线近似程度明显增加。因此，在T 2跨度为0.1ms 到10000ms 之间时，布点数达到600至1000点时，反演曲线与原曲线相似度就很高了。
（3）反演后的曲线对比图和T 2谱
图4.8 m=100时线性均匀布点反演曲线图和T 2谱
图4.9 m=300时线性均匀布点反演曲线图和T 2谱
图4.10 m=600时线性均匀布点反演曲线图和T 2谱
图4.11 m=1000时线性均匀布点反演曲线图和T 2谱
总结：由不同布点方式得出的结果对比，对数均匀布点方式得出的反演曲线与原曲线差距最大，而且只有一个结果；而线性均匀布点方式，在取点数较少时，得出的反演曲线与原曲线差异较大。综合来说，2的指数幂布点方式得出的反演曲线与原曲线最接近。
4.2 问题二：求T 2谱的一般模型 4.2.1 问题分析
横向弛豫时间T 2分布跨度大，通常范围从0.1ms 到10000ms ，导致方程的
测量数据很小的偏差也导致最后计算结果偏差y (t )求解是个典型的不适定问题，
很大。现有各种T 2谱反演方法分布都是事先采取某种对数均匀布点得到，然后由方程y (t )得到一个线性方程组，再利用NNLS 算法求解。而对f j 非负约束的处理，都是采用非负迭代的方法处理。由于T 2分布事先给定，导致这些方法都存在一个普遍的问题：就是在测量T 2组分离散且分布较宽时，反演得到的T 2谱分辨率较低。而且这种基于迭代处理非负约束的优化问题，有时候收敛的比较慢，大大增
加了计算量。
为避免这些问题考虑在未预先给定T 2j 分布的情况下，通过求解带非负约束的优化问题，反演出T 2j 和f j ，即T 2谱。考虑到原问题的病态性质和测量的误差，一般通过曲率平滑的方法，将原问题转化为带非负约束的优化问题。 4.2.2 模型建立
根据上述分析，建立模型如下： 首先利用范数平滑平滑公式：
n
min
∑
i =1
⎛ y - i ⎝
m -
t i T 2
j
∑
j =1
f i ∙e
m ⎫
2⎪+αf j ∑⎪j =1
⎭
2
s . t .
f j ≥0, 1≤j ≤m
其中f j 为第j 类孔隙在总孔隙中所占的份额，T 2j 为第j 类孔隙的T 2驰豫时间，通常范围为0.1m s ≤T 2j ≤10000m s ，τ为回波间隔时间，t 为采集时间。 接着利用曲率平滑公式：
n
min
∑
i =1
⎛ y - i ⎝
m -
t i T 2
j
∑
j =1
f i ∙e
m -1⎫
2⎪+α(f j +1-2f j +f j -1) ∑⎪j =1
⎭
2
s . t .
f j ≥0, 1≤j ≤m
其中，α为正则化参数，它的选取与原始数据的误差水平相关。 4.2.3 模型求解
对本文模型中的非线性优化问题，采取差分进化算法求解。差分进化算法的步骤如下：为求非线性函数f (x 1, x 2,..., x D )的最小值，进化过程中的第G 代利用
NP
个D 维参数向量构成种群{x i , G , i =1, 2,..., NP }，其中种群大小NP 在整个进
化过程中保持不变。 Step1：种群初始化
在D 维空间，随机产生满足约束条件的NP 个向量：
x ji , 0=rand j ∙(b j , U -b j , L )+b j , L 构成初始种群，其中rand j 为区间(0 1)上的随机数。b j , U , b j , L 量x j 的上限和下限。 Step2：变异
对每一个目标向量，x i , G 在{1, 2,..., NP }范围内随机选择三个互易整数r 1, r 2, r 3, 且使得r j ≠i (1≤i ≤3)，得到变异变量
v i , G +1=x r , G +F ∙(x r
1
2
分j 个变
, G
-x r 3, G ) ，
其中，变异因子F ∈(0 , 2)为一常数，控制偏差变量的放大程度。 Step3：交叉
交叉是为了增加种群的多样性，通过如下方式得到试验向量：
u ji , G +1
⎧⎪v ji , G +1, rand j ≤CR or =⎨⎪⎩x ji , G +1, rand j ≥CR or
j =rnbr i j ≠rnbr i
其中j =1, 2,..., D ，式中rand j 为(0 , 1)之间的随机数，rnbr ∈{1, 2,..., D }为随机选择j 的序列，用它来确保u i , G +1从v i , G +1至少获取一个分量，交叉因子CR ∈[0, 1]。 Step4：选择
u i , G +1的目标函数值小于的 x i , G +1差分进化的选择方案基于局部竞争机理，如果
x i , G +1等于 x i , G +1等于 x i , G 。 u i , G +1；反之，令 目标函数值，则令
Step5：边界条件处理
在有边界约束的问题中，必须确保产生新个体的参数值落在问题的可行域中，一个简单方法是将不符合边界约束的新个体用在可行域中随机产生的参数向量代替。
计算结果表明，基于曲率光滑的正则化优化算法能更好地抑制f (T 2)的振荡，对正则化参数的选取不如范数平滑时敏感，对此类多指数反演问题更加有效。 4.3 问题三：减小测量数据的误差的影响 4.3.1 问题分析
由于测量方式和测量工具的差异，测量数据必然会存在误差。前两问中考虑到测量数据的误差，均通过分析去噪信号，求解测量信号的T 2频谱。在实际检测计算过程中，需要考虑测量误差对结果的影响，因此需要在数据处理和建模求解过程中，尽可能地降低测量误差带来的影响。
测量误差可视为测量过程各种干扰导致的结果。为了衡量测量误差程度，此处，引入信噪比概念（SNR ）[4]，信噪比SNR 按如下方式定义：
SNR =20log
||y ||
10
~2
（2）
2
||y -y ||
~
其中y 表示有效信号，即去噪信号，y 表示原信号，即未经去噪处理的信号。SNR 数值越大，表明信号受干扰越少，即测量误差越小。
本文采用计算机模拟来分析测量误差对于前面两个问题带来的影响。首先给定一个具有双峰特性的T 2谱分布y (t ) =200e -t /16+800e -t /1024（单位：ms ），由该频谱可计算出不同时刻t i 可计算出信号强度y i ，再对y 加上一定大小的噪音得到
y ，由数据可计算出受干扰信号的信噪比，模拟信号的噪声由
~
~
matlab7.0随机生
成。对于模拟所得的信号y ，在预先给定驰豫时间分布T 2j 和未预先给定驰豫时间分布T 2j 两种情况下，分别采用非负最小二乘算法和差分进化算法反演T 2频谱。最终通过对比得出测量误差对于这两种情况的T 2频谱计算结果的影响。 4.3.2 预先给定驰豫分布的反演结果
由问题一的结果分析可得，采用2的幂指数均匀布点方式曲线拟合相似度较高，因此本问采用2的幂指数均匀布点方式，取T 2j =2j (j =1, 2, 15) （单位：ms ）。分别在信号y (t ) =200e -t /16+800e -t /1024上加以不同程度的噪声，计算模拟信号的信噪比SNR ，在利用问题一中的NNLS 算法对模拟信号进行反演，得出反演曲线的残差和对应的f j 。此处用反演结果中的f j 与原信号f j 的方差作为评判反演结果好坏的指标。
（1）加为强度为0-1的随机噪音
利用matlab7.0计算结果如表4.4和图4.10所示
图4.12 强度为0-1随机噪音信号反演结果和T2谱
（2）加为强度为0-5的随机噪音
利用matlab7.0计算结果如表4.5和图4.11所示
图4.13 强度为0-5随机噪音信号反演结果和T2谱
（3）加为强度为0-20的随机噪音
利用matlab7.0计算结果如表4.6和图4.12所示
图4.14
强度为0-20随机噪音信号反演结果和T2谱
（4）不加随机噪音
利用matlab7.0计算结果如表4.7和图4.13所示
表4.7 不加随机噪音信号反演结果
图4.15 不加随机噪音信号反演结果和T2谱
结论：通过对比不同程度噪音干扰下的反演结果可知，NNLS 算法在信噪比较高的情况下计算比较准确，但在测量误差比较大的情况下，反演得到的结果与实际信号差距很大，已经不能正确计算出较好的结果。
总结：根据上述分析可知，要减小测量误差带来的不良影响，可从两方面下手。一方面在进行驰豫信号反演前，对测量数据进行去噪处理。另一方面，在去噪处理后，，估算去噪信号的SNR ，如果信号SNR 较高，进行驰豫信号反演时采用NNLS 算法或者差分进化算法均可，如果信号SNR 较低，则应采用差分进化算法进行驰豫信号反演。通过上述两种方法，可以减小测量误差为驰豫信号反演结果带来的不良影响。
五、模型评价与推广
5.1 模型的评价
本文采用的非负最小二乘算法（NNLS 法）的特点是数字稳定性好、运算速度快、容易实现，从计算结果来看，NNLS 算法的可信度较高。NNLS 法虽然得不到精确的T 2值，但是从它的结果中可以基本知道弛豫数据的T 2组分数．而且NNLS 法算出的各组分的T 2值和幅度与理想结果偏差不大，作为非线性拟合的初值比较理想从问题三的结果可以看出，NNLS 算法在信噪比比较低的情况下误差较小；但在信噪比比较大的情况下，该算法就体现出了其局限性。因此，接着采用差分进化算法，能克服NNLS 的上述问题。
基于差分进化的优化算法不依赖于初值选取，计算稳定；更容易实现非负约束，计算速度快。该算法非常适合T 2组分离散且分布较宽，信噪比较低的T 2弛豫信号反演。
5.2 模型的推广
本文中建立的模型均可用来解决多指数反演问题，例如热中子寿命谱的计算，PNN 测井时间谱的计算和地表温度和发射率的计算。
其中主要应用了NNLS 算法和差分进化算法。NNLS 算法几乎适用于任何具有非负约束条件的线性和非线性拟合问题，也可以变相地解决具有非负约束条件的优化问题。而差分进化算法在约束优化计算、聚类优化计算、非线性优化控制、神经网络优化、滤波器设计、阵列天线方向图综合及其它方面得到了广泛的。
参考文献
[1] 刘卫国编著，MATLAB 程序设计教程，北京：中国水利水电出版社，2005.3； [2] 姚绪刚，王忠东，一种新的核磁共振弛豫谱反演算法，测井技术，2003，27(5)：373~376；
[3] 王 鹤，李鲠颖，反演与拟合相结合处理核磁共振弛豫数据的方法，物理学报，第54卷第3期2005年3月；
[4] 潘克家，陈华，谭永基，基于差分进化算法的核磁共振T2谱多指数反演反，物理学报，第57卷第9期2008年9月；
附录
附录一 问题一的求解
附录1.1 对数均匀分布NNLS 算法matlab 源程序 x=linspace(10^0.0001,10^10,n);
t2=log10(x);
A = exp( -kron(te, 1./ t2) );
A = [A ones(size(A,1), 1)];
[x,resnorm]= lsqnonneg(A, y_e);
baseline = x(end);
amplitudes = x(1:end-1);
y_recon = A * x;
figure(1),clf;
subplot(1,2,1);
hold on;
plot(te, y_e', 'rx',te, y_recon, 'b-');
grid on;
xlabel('采样时间/s');
ylabel('磁化信号强度');
legend('实测曲线',' 反演曲线')
title('对数均匀布点m=n');
subplot(1,2,2);
semilogx(t2, amplitudes);
xlabel('T2 (s)');
ylabel('f2');
title('对应T2谱');
grid on;
附录1.2 2的幂指数均匀分布NNLS 算法matlab 源程序 x=linspace(log2(0.1),log2(10000),n);
t2=2.^x;
A = exp( -kron(te, 1./ t2) );
A = [A ones(size(A,1), 1)];
[x,resnorm]= lsqnonneg(A, y_e);
baseline = x(end);
amplitudes = x(1:end-1);
y_recon = A * x;
figure(1),clf;
subplot(1,2,1);
hold on;
plot(te, y_e', 'rx',te, y_recon, 'b-');
grid on;
xlabel('采样时间/ms');
ylabel('磁化信号强度');
legend('实测曲线',' 反演曲线')
title('2的幂指数均匀布点m=n');
subplot(1,2,2);
semilogx(t2, amplitudes);
xlabel('T2 (ms)');
ylabel('f2');
title('对应T2谱');
grid on;
附录1.3 线性均匀分布NNLS 算法matlab 源程序 t2 =linspace(0.1,10000,n);
A = exp( -kron(te, 1./ t2) );
A = [A ones(size(A,1), 1)];
[x,resnorm]= lsqnonneg(A, y_e);
baseline = x(end);
amplitudes = x(1:end-1);
y_recon = A * x;
figure(1),clf;
subplot(1,2,1);
hold on;
plot(te, y_e', 'rx',te, y_recon, 'b-');
grid on;
xlabel('采样时间/ms');
ylabel('磁化信号强度');
legend('实测曲线',' 反演曲线')
title('线性均匀布点m=n');
subplot(1,2,2);
semilogx(t2, amplitudes);
xlabel('T2 (ms)');
ylabel('f2');
title('对应T2谱');
grid on;
附录二 问题二的求解
附录三 问题三的求解
附录3.1 模拟T2频谱及反演过程的matlab 源程序 te = [10:10:320]';
y = 200 * exp(-te ./ 16 ) + 800 * exp(-te ./1024); noise_factor = 20;
e = noise_factor * randn(size(y));
y_e = y +e;
a=1:15;
t2 =2.^a;
A = exp( -kron(te, 1./ t2) );
A = [A ones(size(A,1), 1)];
[x,resnorm]= lsqnonneg(A, y_e);
baseline = x(end);
amplitudes = x(1:end-1);
subplot(1,2,2)
semilogx(t2, amplitudes,'r:');
c=[16,1024];d=[200,800];
hold on;
stem(c,d,'b');
legend('反演信号',' 模拟信号')
xlabel('T2 (ms)');
ylabel('fj');
title('模拟信号T2谱及其反演T2谱');
grid on;
y_recon = A * x;
subplot(1,2,1)
plot(te, y_e', 'rx',te, y_recon, 'b-');
legend('模拟信号',' 反演信号')
grid on;
xlabel('采样时间(ms)');
ylabel('信号磁化强度');
title('模拟信号及其反演信号');
附录3.2 计算SNR 的matlab 源程序
function y=snr(x1,x2);%x1是原始信号，x2是去噪信号 N=length(x1);
y1=sum(x1.^2);
y2=sum((x1-x2).^2);
y=10*log((y1/y2));