CN108120944B - 一种加权迭代的低场核磁共振t2谱反演算法 - Google Patents

Abstract

本发明涉及一种加权迭代的低场核磁共振T2谱反演算法，包括如下步骤：(1)读取低场核磁共振设备采集得到的原始数据文件；(2)对原始数据进行预处理操作得到反演核心矩阵K以及原始数据中各回波波峰的回波时刻的信号幅值组成的向量m；(3)计算采集数据的信噪比、曲率、斜率，进而获取反演权重矩阵；(4)利用反演权重矩阵进行加权迭代求解式：m＝Ks，得到s的最优解，s表示横向弛豫时间所对应的物质的含量组成的向量；(5)根据最优解绘制T2谱。与现有技术相比，本发明对短弛豫组分进行加权，更贴近原始信号的稀疏性，计算精度高，鲁棒性好，在不同信噪比数据情况下，能够得到稳定的反演结果。

Description

一种加权迭代的低场核磁共振T2谱反演算法

技术领域

本发明涉及一种核磁共振信号处理技术，尤其是涉及一种加权迭代的低场核磁共振T2谱反演算法。

背景技术

在世界范围内，核磁共振技术发展迅猛，在很多领域得到了广泛的发展，比如高场强核磁共振技术应用于对人体的临床诊断等。随着核磁共振技术的不断发展，人们发现核磁共振技术不仅可以应用于临床，还可以在其他领域(如食品科学、农业、石油能源、材料科学、纺织化工等)发挥其他科学仪器所不能发挥的作用。很多科学研究已经证明，在这些领域中应用核磁共振分析技术，可以解决现有的其他科学仪器所不能解决的问题，对这些领域的科学进步起到了举足轻重的推动作用。与用于临床的高场核磁共振技术不同的是，在食品科学、农业、石油能源、材料科学、纺织化工等领域采用低场核磁共振技术就可以很好的解决相关问题。此外，低场核磁共振分析仪器由于采用低场强磁体，可以使所研发的分析仪器体积较小，安装、调试、维护、操作都很方便。因此，低场核磁共振分析仪器得到了科学界的关注。

CPMG(Carr-Purcell-Meiboom-Gill)序列速度快，是低场核磁共振中最常用的 T2序列之一。研究者常常利用CPMG等序列的原始数据和样品横向弛豫时间、纵向弛豫时间的分布特点，进行时域谱反演的相关研究。但是核磁共振采集到的原始信号并不能得到直观的样品结构信息，需要通过反演技术才能得到易于理解的时域谱。一般，一维时域谱就可以为研究者分析样品的组成、性质等提供重要的依据。但是，随着低场核磁共振应用的不断深入，研究者们发现实验得到的T2谱在反演时长弛豫组分会把短弛豫组分吸收掉，进而不能反映短弛豫组分的真实分布。

国内外通常使用三种快速T2反演的思路：第一种是截断方法，比如奇异值截断法；第二种正则化方法，如标准正则化法、L曲线法；第三种迭代法，比如SIRT 法等。现有正则化因子的选择算法需要一定的人为干预，还待完善。

CPMG序列中短弛豫组分衰减速度比较快，现有的算法都是平等对待各个数据，这样使得反演计算时长弛豫组分容易将短弛豫组分吸收掉合成一个大的谱峰。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种加权迭代的低场核磁共振T2谱反演算法。

本发明的目的可以通过以下技术方案来实现：

一种加权迭代的低场核磁共振T2谱反演算法，包括如下步骤：

(1)读取低场核磁共振设备采集得到的原始数据文件；

(2)对原始数据进行预处理操作得到反演核心矩阵K以及原始数据中各回波波峰的回波时刻的信号幅值组成的向量m；

(3)计算采集数据的信噪比、曲率、斜率，进而获取反演权重矩阵；

(4)利用反演权重矩阵进行加权迭代求解式：m＝Ks，得到s的最优解，s表示横向弛豫时间所对应的物质的含量组成的向量；

(5)根据最优解绘制T2谱。

步骤(2)具体为：

(21)提取低场核磁共振设备测定的各回波波峰的回波时刻，组成向量τ₂，并获取各回波波峰的回波时刻的信号幅值组成向量m；

(22)对横向弛豫时间进行布点，将横向弛豫时间组成向量T₂；

(23)根据向量τ₂和m计算反演核心矩阵K，反演核心矩阵中第i行第j列元素为K^i,j：

其中，τ_2,i表示向量τ₂中的第i个元素，T_2,j表示向量T₂中的第j个元素， i＝1,2……I，j＝1,2……J，I表示回波波峰的个数，J表示横向弛豫时间布点个数。

步骤(3)具体为：

(31)计算原始数据的信噪比，绘制τ₂和m构成的曲线L，计算曲线L的曲率和斜率；

(32)构造反演权重矩阵W，反演权重矩阵对角线上元素为权重，其余元素均为0，具体地：

确定权重个数：权重个数等于衰减时间在10ms之内的数据个数；

确定反演权重矩阵上对角线元素的具体分布：最小权重为1，根据信噪比、曲率、斜率确定最大权重为w_max，然后在以2为底的指数空间中产生由最大权重到最小权重等间距的权重系数，最后按照权重系数由大到小分布在反演权重矩阵对角线上。

最大权重为w_max具体为：

w_max＝α_snr×(α_cC_max+α_pP_max)，

α_snr为信噪比的倒数，C_max为曲线L的曲率最大值，P_max为曲线L的斜率最大值，α_c和α_p均为常数。

步骤(4)具体为：

(41)对m和K进行加权得到K^w和m^w：K^w＝K^TWK，m^w＝K^TWm，W为反演权重矩阵；

(42)令迭代次数q＝0，给定s的初始值s^q；

(43)令s＝s^q带入公式：m＝K^ws，计算得到m^q；

(44)求取Δm^q＝m^w-m^q；

(45)根据Δm^q求取误差分配矩阵Δs^q；

(46)计算s^(q+1)＝s^q+Δs^q；

(47)令q＝q+1，返回步骤(43)进行迭代计算直到满足迭代终止条件得到s 的最优解。

步骤(45)Δs^q中第j个元素为：

nnz(K^w(j))表示矩阵K^w中第j列中的非零元素个数，k_ij表示矩阵K^w中第i行第j列的元素，表示向量Δm^q中第i个元素，i＝1,2……I，I表示回波波峰的个数，j＝1,2……J，J表示横向弛豫时间布点个数。

与现有技术相比，本发明具有如下优点：

本发明引入反演权重矩阵，对短弛豫组分进行加权，更贴近原始信号的稀疏性，得到的解更能体现真实的谱分布；计算精度高；鲁棒性好，在不同信噪比数据情况下，能够得到稳定的反演结果。

附图说明

图1为本发明加权迭代的低场核磁共振T2谱反演算法的流程框图。

图2的(a)-图2的(d)是信噪比为1000的仿真实验结果示意图：其中，图 2的(a)是采用TSVD方法的反演结果；图2的(b)是采用SIRT方法的反演结果；图2的(c)是采用BRD方法的反演结果；图2的(d)是采用本发明方法的反演结果。

图3的(a)-图3的(d)是信噪比为100的仿真实验结果示意图：其中，图3的(a)是采用TSVD方法的反演结果；图3的(b)是采用SIRT方法的反演结果；图3的(c)是采用BRD方法的反演结果；图3的(d)是采用本发明方法的反演结果。

图4的(a)-图4的(d)是信噪比为10的仿真实验结果示意图：其中，图4 的(a)是采用TSVD方法的反演结果；图4的(b)是采用SIRT方法的反演结果；图4的(c)是采用BRD方法的反演结果；图4的(d)是采用本发明方法的反演结果。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

实施例

如图1所示，一种加权迭代的低场核磁共振T2谱反演算法，包括如下步骤：

(1)读取低场核磁共振设备采集得到的原始数据文件；

(2)对原始数据进行预处理操作得到反演核心矩阵K以及原始数据中各回波波峰的回波时刻的信号幅值组成的向量m；

(3)计算采集数据的信噪比、曲率、斜率，进而获取反演权重矩阵；

(4)利用反演权重矩阵进行加权迭代求解式：m＝Ks，得到s的最优解，s表示横向弛豫时间所对应的物质的含量组成的向量；

(5)根据最优解绘制T2谱。

A、步骤(2)的原理基于如下：

一维反演问题就是求解如式(1)所示的具有一个核的Fredholm积分方程，τ₂表示使用低场核磁共振序列测定的各回波波峰的回波时刻组成的向量，s表示横向弛豫时间所对应的物质的含量组成的向量，s(T₂)表示某个特定的横向弛豫时间所对应的物质的含量，m表示各回波波峰的回波时刻的信号幅值组成的向量，m(τ₂)表示向量τ₂中某个特定时刻的采样数据的信号幅值。T_2max、T_2min分别表示T₂的最大布点值与最小布点值，一般采用以10为底的对数空间内均匀布点。

将公式(1)中的弛豫时间、离散化后可以用矩阵的形式表示如公式 (2)所示：

m＝Ks， (2)

其中，K表示反演核心矩阵，而K可以通过向量τ₂和m计算得到，因此，这样转换之后，反演问题就变成一个已知m和K求s的问题，其中矩阵K称为反演的核心矩阵。为了得到符合真实情况的谱分布，要对T₂进行布点，布点数目根据具体实验情况而定，实验精度要求越高，布点数目越多，一般为了节约计算时间将布点数目设为64。

基于上述原理，步骤(2)具体为：

(21)提取低场核磁共振设备测定的各回波波峰的回波时刻，组成向量τ₂，并获取各回波波峰的回波时刻的信号幅值组成向量m；

(22)对横向弛豫时间进行布点，将横向弛豫时间组成向量T₂；

(23)根据向量τ₂和m计算反演核心矩阵K，反演核心矩阵中第i行第j列元素为K^i,j：

其中，τ_2,i表示向量τ₂中的第i个元素，T_2,j表示向量T₂中的第j个元素， i＝1,2……I，j＝1,2……J，I表示回波波峰的个数，J表示横向弛豫时间布点个数。

B、步骤(3)的原理基于如下：

在T2反演中，常用的权重施加方式如下：

K^TWm＝K^TWKs， (3)

W就是反演权重矩阵，权重由它的对角元素决定。然后令K^w＝K^TWK， m^w＝K^TWm，这样公式(2)可以转化为：

m^w＝K^ws。 (4)

基于上述原理，步骤(3)具体为：

(31)计算原始数据的信噪比，绘制τ₂和m构成的曲线L，计算曲线L的曲率和斜率；

(32)构造反演权重矩阵W，反演权重矩阵对角线上元素为权重，其余元素均为0，具体地：

确定权重个数：权重个数等于衰减时间在10ms之内的数据个数；

确定反演权重矩阵上对角线元素的具体分布：最小权重为1，根据信噪比、曲率、斜率确定最大权重为w_max，然后在以2为底的指数空间中产生由最大权重到最小权重等间距的权重系数，最后按照权重系数由大到小分布在反演权重矩阵对角线上。

最大权重为w_max具体为：

w_max＝α_snr×(α_cC_max+α_pP_max)，

α_snr为信噪比的倒数，C_max为曲线L的曲率最大值，P_max为曲线L的斜率最大值，α_c和α_p均为常数。

C、步骤(4)的原理如下：

首先给定出s的初始值s'，根据式(4)计算m'，然后计算模型信号与实际信号差值：

Δm＝m^w-m′，

那么待求解的s的误差Δs同样满足式(4)：

Δm＝K^w·Δs，

它们之间的关系可以表述如下：

Δs_j为向量Δs中第j个元素，k_ij表示矩阵K^w中第i行第j列的元素，Δm_i表示向量Δm中第i个元素，i＝1,2……I，I表示回波波峰的个数，j＝1,2……J，J表示横向弛豫时间布点个数。设置Δs_j＝λk_ijΔm_i，联立这些上式可以得到：

然后就得到了误差分配公式：

nnz(K^w(j))表示矩阵K^w中第j列中的非零元素个数，经过上面的计算可以得到 j从1取到J的误差向量Δs，s^q＝s'+Δs，q表示迭代次数，接下来将得到的s^q作为新的s'继续执行上述过程完成迭代。

因此，根据上述原理，步骤(4)具体为：

(41)对m和K进行加权得到K^w和m^w：K^w＝K^TWK，m^w＝K^TWm，W为反演权重矩阵；

(42)令迭代次数q＝0，给定s的初始值s^q；

(43)令s＝s^q带入公式：m＝K^ws，计算得到m^q；

(44)求取Δm^q＝m^w-m^q；

(45)根据Δm^q求取误差分配矩阵Δs^q；

(46)计算s^(q+1)＝s^q+Δs^q；

(47)令q＝q+1，返回步骤(43)进行迭代计算直到满足迭代终止条件得到s 的最优解，迭代终止条件为迭代次数q到达最大迭代次数或者拟合误差小于根据样品设定的限定值。

步骤(45)Δs^q中第j个元素为：

nnz(K^w(j))表示矩阵K^w中第j列中的非零元素个数，k_ij表示矩阵K^w中第i行第j列的元素，表示向量Δm^q中第i个元素，i＝1,2……I，I表示回波波峰的个数，j＝1,2……J，J表示横向弛豫时间布点个数。

D、根据最优解绘制T2谱

设s的最优解为s″，s″包括J个元素，分别表示J个横向弛豫时间布点处所对应的物质的含量，T₂表示横向弛豫时间组成向量，由此，根据T₂和s″绘制出T2谱。

本发明加权迭代的低场核磁共振T2谱反演算法的效果能够通过以下实验进一步说明。

1.仿真实验：

实验首先构造一个中心分别位于T₂＝0.2ms和T₂＝100ms处的高斯峰作为理想的T2谱，然后向正演结果中加入一定程度的高斯白噪声，得到不同信噪比的仿真数据。

2.仿真实验结果

图2、图3、图4分别是信噪比为1000、信噪比为100和信噪比为10时采用三种不同的方法进行反演得到的仿真结果，这三幅图中(a)图显示的采用TSVD 方法反演结果；(b)图是采用SIRT方法的反演结果；(c)图是采用BRD方法的反演结果；(d)图是采用本发明算法的反演结果。

为了进一步验证本发明算法的抗噪性能，表1列举了信噪比分别为1000、500、300、200、100、50、10时四种方法仿真数据的处理结果。

表1不同信噪比下的4种算法的S/L和相对拟合误差比较 (S/L为短弛豫组分与长弛豫组分积分值之比)

3.仿真实验分析

通过以上图形及表1的结果，我们可以得到以下结论：在高信噪比下，4种方法都能得到反演结果，并且TSVD的结果更好，但是这种算法鲁棒性差，信噪比降低变化太大；SIRT算法在同等迭代次数下，显示短弛豫组分的能力最差，在不同信噪比下几乎都无法显示；BRD与本发明方法对于短弛豫组分显示都具有一定的优越性，在拟合误差方面，BRD算法优于本发明的算法，但是在S/L方面，本发明的算法明显好于BRD算法；在低信噪比下，TSVD算法和BRD算法都容易产生虚假小峰。综上所述，本发明提出的算法具有很高的计算精度，对于短弛豫组分显示具有极大优势，并且鲁棒性好，抗干扰能力强，能够得到更加符合真实谱分布的T2谱。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员，在不脱离本发明方法的前提下，还可以做出若干改进和补充，这些改进和补充也应视为本发明的保护范围。

Claims (1)

1.一种加权迭代的低场核磁共振T2谱反演算法，其特征在于，包括如下步骤：

(1)读取低场核磁共振设备采集得到的原始数据文件；

(2)对原始数据进行预处理操作得到反演核心矩阵K以及原始数据中各回波波峰的回波时刻的信号幅值组成的向量m；

(3)计算采集数据的信噪比、曲率、斜率，进而获取反演权重矩阵；

(4)利用反演权重矩阵进行加权迭代求解式：m＝Ks，得到s的最优解，s表示横向弛豫时间所对应的物质的含量组成的向量；

(5)根据最优解绘制T2谱；

步骤(2)具体为：

(21)提取低场核磁共振设备测定的各回波波峰的回波时刻，组成向量τ₂，并获取各回波波峰的回波时刻的信号幅值组成向量m；

(22)对横向弛豫时间进行布点，将横向弛豫时间组成向量T₂；

(23)根据向量τ₂和m计算反演核心矩阵K，反演核心矩阵中第i行第j列元素为K^i,j：

其中，τ_2,i表示向量τ₂中的第i个元素，T_2,j表示向量T₂中的第j个元素，i＝1,2……I，j＝1,2……J，I表示回波波峰的个数，J表示横向弛豫时间布点个数；

步骤(3)具体为：

(31)计算原始数据的信噪比，绘制τ₂和m构成的曲线L，计算曲线L的曲率和斜率；

(32)构造反演权重矩阵W，反演权重矩阵对角线上元素为权重，其余元素均为0，具体地：

确定权重个数：权重个数等于衰减时间在10ms之内的数据个数；

确定反演权重矩阵上对角线元素的具体分布：最小权重为1，根据信噪比、曲率、斜率确定最大权重为w_max，然后在以2为底的指数空间中产生由最大权重到最小权重等间距的权重系数，最后按照权重系数由大到小分布在反演权重矩阵对角线上；

最大权重为w_max具体为：

w_max＝α_snr×(α_cC_max+α_pP_max)，

α_snr为信噪比的倒数，C_max为曲线L的曲率最大值，P_max为曲线L的斜率最大值，α_c和α_p均为常数；

步骤(4)具体为：

(41)对m和K进行加权得到K^w和m^w：K^w＝K^TWK，m^w＝K^TWm，W为反演权重矩阵；

(42)令迭代次数q＝0，给定s的初始值s^q；

(43)令s＝s^q带入公式：m＝K^ws，计算得到m^q；

(44)求取Δm^q＝m^w-m^q；

(45)根据Δm^q求取误差分配矩阵Δs^q；

(46)计算s^(q+1)＝s^q+Δs^q；

(47)令q＝q+1，返回步骤(43)进行迭代计算直到满足迭代终止条件得到s的最优解；

步骤(45)Δs^q中第j个元素为：

nnz(K^w(j))表示矩阵K^w中第j列中的非零元素个数，k_ij表示矩阵K^w中第i行第j列的元素，表示向量Δm^q中第i个元素，i＝1,2……I，I表示回波波峰的个数，j＝1,2……J，J表示横向弛豫时间布点个数。