CN114004764B

CN114004764B - 一种基于稀疏变换学习的改进灵敏度编码重建方法

Info

Publication number: CN114004764B
Application number: CN202111296938.9A
Authority: CN
Inventors: 段继忠; 李玺兰
Original assignee: Kunming University of Science and Technology
Current assignee: Kunming University of Science and Technology
Priority date: 2021-11-03
Filing date: 2021-11-03
Publication date: 2024-03-15
Anticipated expiration: 2041-11-03
Also published as: CN114004764A

Abstract

本发明涉及一种基于稀疏变换学习的改进灵敏度编码重建方法，属于磁共振成像技术领域。灵敏度编码(SENSE)是一种显式利用多个接收线圈灵敏度信息来减少扫描时间的技术。为了提高磁共振重建质量、减少重建的伪影，本发明将数据驱动的自适应稀疏变换学习TL引入SENSE模型中，提出了一种基于灵敏度编码结合变换学习正则项的并行磁共振成像重建算法(Transform Learning‑Sensitivity Encoding,TL‑SENSE)。该模型利用交替方向乘子法(Alternating Direction Method Of Multipliers,ADMM)进行求解，通过变换更新、硬阈值去噪和图像更新三步实现磁共振成像重建。仿真实验结果表明，与基于全变分(Total Variation，TV)和Lp伪范数全变分(LpTV)正则项的灵敏度编码重建模型相比，提出的模型对图像去噪和修复具有较好的效果。

Description

一种基于稀疏变换学习的改进灵敏度编码重建方法

技术领域

本发明涉及一种基于稀疏变换学习的改进灵敏度编码重建方法，属于磁共振成像技术领域。

背景技术

磁共振成像(Magnetic Resonance Imaging,MRI)技术是目前临床医学重要的影像学工具之一。如何加快磁共振成像扫描速度及提高重建质量广受关注。

并行成像(Parallel Imaging,PI)和压缩感知(Compressed Sensing,CS)是加速磁共振成像的常用技术。并行成像技术通过使用多个灵敏度不同的接收线圈同时采集信号。灵敏度编码(SENSitivity Encoding,SENSE)模型显式利用灵敏度信息进行并行磁共振成像的重建，可显著提高磁共振成像的性能。2001年，Pruessmann等人结合了网格划分原理和共轭梯度迭代法提出针对任意k空间轨迹的SENSE技术，减少了非笛卡尔SENSE重构所需的时间。并行成像有效的加速了磁共振成像，且提高了成像的鲁棒性，但是边缘的保留完整度较差，灵敏度估计也存在误差。

压缩感知理论表明：若信号具有有效稀疏性，则可以从欠采样数据中精确恢复出原信号。为了弥补传统线圈灵敏度估计和并行磁共振图像重建技术的不足，Ying等人于2007年提出了一种联合估计线圈灵敏度和重建所需图像的方法，该方法从线圈灵敏度的初始估计开始，在每次迭代中交替更新重建图像和线圈灵敏度。2008年Liu等人将压缩感知与SENSE模型结合，引入全变分(Total Variation,TV)正则项，即SparseSENSE，该方法能有效的提高重构图像质量。2011年，Ramani等人提出利用增广拉格朗日方法求解含有全变分(Total Variation,TV)正则项、小波变换的L1范数正则项的SENSE模型重构问题，即TV-SENSE，该算法有较快的收敛速度。2020年，鲍中文等人提出基于Lp伪范数全变分(LpTV)正则项的图像重建算法，利用变量分裂技术将复杂的LpTV-SENSE问题解耦成多个容易求解的子问题进行求解。由此可见，正则项的引入使SENSE模型在抑制噪声和去除伪影方面表现出良好的性能。

基于灵敏度编码技术进行线圈灵敏度估计重建图像需要正则化来抑制噪声和减少混叠效应。目前常见的基于正则化器的灵敏度编码技术，主要有基于TV正则项的TV-SENSE算法和基于Lp伪范数全变分(LpTV)正则项的LpTV-SENSE算法。但这些算法处理的图像仍存在问题：重建图像存在阶梯伪影和模糊伪影，细节和边缘保留完整度较低，目标图像与原始图像一致性较低。

发明内容

本发明的目的在于克服现有技术的不足，提供基于稀疏变换学习的改进灵敏度编码重建方法，能进一步提高磁共振成像重建图像的质量。

本发明的目的是通过以下技术方案来实现的：基于稀疏变换学习的改进灵敏度编码重建方法，它包括以下几个步骤：

S0：初始化，令k＝0，x⁰＝(RFS)^Hy，X⁰＝[P₁(x⁰),...,P_J(x⁰)]，B⁰＝H(W⁰X⁰,α/μ₂)；

其中，k是循环变量，表示数据的第k次迭代，上标“⁰”表示初始值；表示列向量化的待重建图像，N＝m×n，m和n分别是图像的行数和列数，x⁰表示x的初始值；表示列向量化的欠采样多线圈k空间数据，M为单线圈k空间数据实际采样的点数，且M＜＜N，L表示并行成像所使用的接收线圈个数，上标“H”表示复共轭转置。y的第l个线圈数据用y_l表示，l＝1,...,L表示线圈索引，例如，y₁和y_L分别表示列向量化的欠采样多线圈k空间数据y的第1个线圈数据和第L个线圈数据。/> 是对角矩阵，S_l对应第l个接收线圈的灵敏度，例如：S₁表示第1个线圈的灵敏度图，S_L表示第L个线圈的灵敏度图；/> 和/>分别为n×n和m×m的傅里叶变换矩阵，I_L为L×L的单位矩阵，/>表示克罗内克积；/>是一个欠采样矩阵，其中/>表示从单线圈的k空间网格中选择采样点位置的矩阵。

P_j(·):表示第j个图像块提取线性算子，从单个线圈图像提取大小为的图像块再列向量化成/>的列向量，例如：P₁(·)表示第1个图像块提取线性算子，P_J(·)表示第J个图像块提取线性算子。/>表示/>图像块的稀疏变换矩阵，W⁰表示W的初始值，/>表示/>点离散余弦变换矩阵。/>表示WP_j(x)对应的辅助变量，将所有B_j按顺序水平拼接成矩阵/>B_j表示B的第j个分量，J表示图像块的个数，例如：B₁表示B的第1分量，B_J表示B的第J个分量，B⁰表示B的初始值。中间变量/>表示从x中抽取J个/>的图像块并水平拼接成的矩阵，P₁(x)和P_J(x)分别表示从单线圈图像x中提取的第1个图像块和第J个图像块，X⁰＝[P₁(x⁰),...,P_J(x⁰)]表示X＝[P₁(x),...,P_J(x)]的初始值。/>是硬阈值函数，/> 表示输入矩阵，θ表示阈值，α,μ₂为参数，且α＞0,μ₂＞0。/>为FSx对应的辅助变量，/>为z对应的拉格朗日乘子，/>表示u_z的初始值。

S1：对X^k(B^k)^H进行奇异值分解，得到X^k(B^k)^H＝UΣV^H；

其中，B^k是辅助变量B的第k次迭代，X^k是中间变量X的第k次迭代，U,Σ,V是奇异值分解因子。

S2：变换更新，计算第k+1次迭代的稀疏变换W^k+1，计算公式如下：

W^k+1＝VU^H。

S3：硬阈值去噪，计算第k+1次迭代的辅助变量B^k+1，计算公式如下：

B^k+1＝H(W^k+1X^k,α/μ₂)。

S4：计算第k+1次迭代的辅助变量z^k+1，计算公式如下：

其中，x^k表示列向量化的待重建图像x的第k次迭代，R^T表示从k空间网格中选择采样点位置的矩阵的转置，表示拉格朗日乘子u_z的第k次迭代，μ₁＞0，(·)^-1表示求逆算子。

S5：计算第k+1次迭代的待重建图像x^k+1，计算公式如下：

其中，上标“^*”是伴随算子，是P_j(·)的伴随算子，/>是对角矩阵，其对角元素等于对应的像素点出现在所有图像块中的次数，μ₁＞0,μ₂＞0，/>是B_j的第k+1次迭代。

S6：计算第k+1次迭代的中间变量X^k+1，计算公式如下：

X^k+1＝[P₁(x^k+1),...,P_J(x^k+1)]。

S7：更新第k+1次迭代的拉格朗日乘子

S8：计算x^k+1和x^k之间的相对误差(RE)，计算式为：

S9：判断是否满足算法停止条件，若满足RE＜tol，则进入S10，否则令k＝k+1，返回S1。

S10：输出x＝x^k+1，得到重建图像x。

本发明的有益效果是：SENSE是一种并行磁共振成像的重建技术，通过求解关于线圈灵敏度的线性系统来重建图像。本发明将数据驱动自适应稀疏变换学习(TransformLearning,TL)引入灵敏度编码(SENSE)技术中，进一步加强了灵敏度编码技术的性能，确保MRI重建图像的质量。该模型利用ADMM技术将TL-SENSE模型分解为多个容易求解的子问题，然后通过变换更新，硬阈值去噪和图像更新三步实现MRI重建。理论分析和实验结果表明，与其他基于灵敏度编码(SENSE)技术的算法(如TV-SENSE和LpTV-SENSE)相比，本发明提出的新算法去模糊效果显著，边缘保留较完整，重构图像清晰。

附图说明

图1为本发明方法流程图；

图2为使用12通道的头部线圈进行完全采样来获取受试者的体内人脑切片数据，取第95个片层，尺寸为218×170(即dataset 1)；

图3为具有3倍加速因子和24×24中心全采样自校准区域的二维泊松圆盘欠采样掩模图；

图4～6分别为使用TV-SENSE算法、LpTV-SENSE算法和TL-SENSE算法从含有3倍加速因子和24×24中心全采样自校准区域的二维泊松圆盘欠采样掩模进行欠采样的dataset1数据重建出的图像；

图7～9分别为图4～6重建图像对应的误差图，也就是由L1-SENSE算法、LpTV-SENSE算法和TL-SENSE算法重建得到的误差图。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案。

实施例1：本发明是基于SENSE框架提出的一种高效的重建方法。

1)SENSE框架：

假设表示列向量化的待重建图像，N＝m×n，m和n分别是图像的行数和列数；/>表示列向量化的欠采样多线圈k空间数据，M表示单个线圈k空间欠采样数据点个数，M为单线圈k空间数据实际采样的点数，M＜＜N，L表示并行成像所使用的接收线圈个数，上标"H"表示共轭转置。y的第l个线圈数据用y_l表示，l＝1,...,L表示线圈索引，例如，y₁和y_L分别表示列向量化的欠采样多线圈k空间数据y的第1个线圈数据和第L个线圈数据。

然后，得到SENSE重建模型为：

y＝RFSx (1)

其中，是对角矩阵，S_l对应第l个接收线圈的灵敏度，例如：S₁表示第1个线圈的灵敏度图，S_L表示第L个线圈的灵敏度图。/> 和/>分别为n×n和m×m的傅里叶变换矩阵，I_L为L×L的单位矩阵，/>表示克罗内克积；是一个欠采样矩阵，其中/>表示从单线圈的k空间网格中选择采样点位置的矩阵。

在SENSE中，关于灵敏度算子的计算式为：

其中s_r是的特征值为“1”的特征向量，r表示k空间位置索引变量；/>表示每个位置r的半定矩阵值的卷积，定义为/>定义为/>是针对向量化的多线圈图像进行二维傅里叶变换的矩阵，/>分别为m点和n点傅里叶变换矩阵，I_L是C×C的单位矩阵，/>表示克罗内克积。/>是一个具有矩阵值内核的卷积，定义为：

其中V_||由k空间校准矩阵得到，若校准矩阵定义为A，则矩阵A的奇异值分解为：

A＝UΣV^H (4)

其中Λ是A的奇异值对角矩阵，U和V都是酉矩阵，其列分别是奇异值的左、右奇异向量。将V分解可以得到A的零空间跨度V_⊥和A的行空间跨度V_||。由于V的列是矩阵A的行的基，所有的自校准信息都位于V_||中。

2)算法推导：

尽管结合了LpTV正则项方法的SENSE算法可以一定程度上改善重建质量，但是仍有较大的改进余地。为了进一步提高重建质量，本发明将数据驱动自适应稀疏变换学习(TL)引入SENSE模型中，提出TL-SENSE算法，得到的优化问题表示如下：

式中，α为正则化参数，J表示图像块的个数。表示第j个图像块提取线性算子，从单个线圈图像提取大小为/>的图像块再列向量化成/>的列向量。/>表示/>图像块的稀疏变换矩阵。/>表示WP_j(x)对应的辅助变量，将所有B_j组合成即B＝[B₁,...,B_J]。式(5)利用交替方向乘子法(ADMM)进行求解。

引入辅助变量并引入拉格朗日乘子/>则优化问题(5)可以转化为：

式(7)是第k+1次迭代的辅助变量z^k+1对应的拉格朗日乘子，x^k+1表示第k+1次迭代的待重建图像。

使用ADMM技术，则问题(6)可以分解为以下子问题：

X^k+1＝[P₁(x^k+1),...,P_J(x^k+1)] (13)

在(8)～(13)中，变量的上标“k+1”和“k”表示算法的第k+1次迭代与第k次迭代。为FSx对应的辅助变量，/>表示第k+1次迭代的辅助变量，/>是第k+1次迭代z^k+1对应的拉格朗日乘子，X^k+1表示第k+1次迭代的中间变量，W^k+1表示第k+1次迭代的稀疏变换矩阵，B^k+1表示第k+1次迭代的辅助变量，μ₁＞0,μ₂＞0。x^k表示第k次迭代的待重构图像，x^k ⁺¹表示第k+1次迭代的待重构图像。然后通过更新变换、硬阈值去噪和图像更新分别求出满足上述子问题的解。

更新变换：该步骤求解W，令表示第k次迭代的中间变量x^k中抽取J个图像块组成n×J的矩阵，例如：P₁(x^k)和P_J(x^k)分别表示第k次迭代的中间变量从x^k中提取的第1个图像块和第J个图像块。则问题(8)重写为：

对X^k(B^k)^H进行奇异值分解得到X^k(B^k)^H＝UΣV^H，则(14)的解是唯一的：

W^k+1＝VU^H (15)

其中，U，∑，V是奇异值分解因子。

硬阈值去噪：子问题(9)重写为：

子问题(16)可用硬阈值法求解，其解表示为：

B^k+1＝H(W^k+1X^k,α/μ₂) (17)

其中，是硬阈值函数，/> 表示输入矩阵，θ表示阈值，且α＞0,μ₂＞0。

图像更新：图像更新求解辅助变量z和待重建图像x。更新z时，控制其他变量不变，令(10)的目标函数对z的导数为0，得：

其中R^TR为对角矩阵，表示第k次迭代的待重建图像，/>是第k次迭代的辅助变量z^k对应的拉格朗日乘子，μ₁＞0，则可以求出z^k+1，表示如下：

更新x，令其他变量不变，令(11)的目标函数对x的导数为0，得：

其中W^HW＝I、均是对角矩阵，则：

其中，(·)^-1表示求逆算子，μ₁＞0,μ₂＞0；上标"H"表示共轭转置。为辅助变量，/>属于辅助变量。上标“*”是伴随算子，/>是P_j(·)的伴随算子，是对角矩阵，其对角元素等于对应的像素点出现在所有图像块中的次数，是B_j的第k+1次迭代。

综上所述，所有的子问题均可有效求解。获得了一种新的具有TL正则项的SENSE并行MRI重建算法，称为TL-SENSE算法。当相对误差(RE)低于容差tol，算法停止。x^k+1和x^k的RE定义为：

其中x^k+1表示第k+1次重建图像。

具体流程图如图1所示，其中步骤如下：

其中上标“⁰”表示初始值，x⁰表示待重建图像x的初始值；W⁰是W的初始值，表示/>点离散余弦变换矩阵；B⁰是B的初始值；/>表示中间变量从x中抽取J个图像块并组成n×J的矩阵，例如：P₁(x^k)和P_J(x^k)分别表示第k次迭代的中间变量从x^k中提取的第1个图像和第J个图像，X⁰＝[P₁(x⁰),...,P_J(x⁰)]表示X＝[P₁(x),...,P_J(x)]的初始值；/>是对应于z的拉格朗日乘子，/>表示它的初始值。

S1：对X^k(B^k)^H进行奇异值分解，得到X^k(B^k)^H＝UΣV^H；

S2：更新变换，计算W^k+1，计算公式如(15)；

S3：硬阈值去噪，计算第k+1次迭代的辅助变量B^k+1，计算公式如(17)；

S4：计算第k+1次迭代的辅助变量z^k+1，计算公式如(19)：

S5：计算第k+1次迭代的待重建图像x^k+1，计算公式如(21)：

S6：计算第k+1次的中间变量X^k+1，计算公式如(13)；

S7：更新第k+1次迭代的拉格朗日乘子计算公式如(12)；

S8：计算RE，公式如(22)；

S9：判断是否满足算法停止标准，若满足RE＜tol，则进入S10，否则令k＝k+1，返回S1。

S10：输出x＝x^k+1，得到重建图像x。

实验结果：

在下面的实验中，为了验证TL-SENSE算法的重建性能，我们将提出的TL-SENSE算法分别与TV-SENSE和LpTV-SENSE方法进行比较。所有算法都能得到重建图像，所有实验都在Matlab(MathWorks，Natick，MA)上实现。所有实验均在配置为英特尔酷睿i5-4210U@2.4GHz处理器、8GB内存和Windows10操作系统(64位)的笔记本电脑上进行。

为了对各算法的性能进行比较，本发明使用1张人类脑部切片图进行仿真实验，命名为dataset 1，如图2所示，dataset 1为12通道头部线圈进行完全采样获取受试者的体内人脑切片图，取第95个片层，尺寸为218×170。使用加速因子为AF＝3和24×24中心全采样自校准区域的二维泊松圆盘欠采样掩模，如图3所示，对全采样的数据集进行人工采样生成测试数据集。

所有算法的参数均针对信噪比(Signal Noise Ratio,SNR)最优进行调整，所有实验均用SNR来定量评判重建质量，SNR越高重建效果越好，在感兴趣区域计算SNR值。

SNR定义如下：

其中，MSE表示重建图像和参考图像x之间的均方误差，Var表示x的方差。

首先，对数据集dataset 1下的三个重建算法进行视觉比较，本发明选择加速因子为3时的重建图像进行比较，如图4～6所示。图4～6分别依次展示了算法TV-SENSE、算法LpTV-SENSE和算法TL-SENSE的重建图像。图4是使用算法TV-SENSE进行重建得到的图像，可以看出，该图像存在模糊伪影、部分细节没有重建出来。图5是使用算法LpTV-SENSE进行重建得到的图像，该算法的重建质量优于TV-SENSE算法，但是也存在一些伪影。图6是使用算法TL-SENSE进行重建得到的图像，该算法能有效的去除模糊伪影，有着较清楚的纹理细节和边缘轮廓信息，且重建图像与原始图像的一致性较高。

为了进一步说明本发明提出的新算法TL-SENSE在重建性能方面优于算法TV-SENSE和LpTV-SENSE，图7～9展示了dataset 1在三种算法下的误差图(误差图越白表明误差越大)。由误差图可以看出，TV-SENSE的重建图像存在较大误差，其次是LpTV-SENSE，TL-SENSE的误差图白点最少，具有较小的重建误差，所提出的算法的重建性能较好。

另外，实验通过SNR值对三种重建算法的性能进行了定量评估，本发明提出的算法TL-SENSE对重建图像的SNR有明显提高。

综上所述，对选用的数据集dataset 1在不同的重建算法中进行实验，TL-SENSE算法在评估指标和视觉上都表现出了自身的优势，视觉效果明显优于TV-SENSE和LpTV-SENSE算法。

以上结合附图对本发明的具体实施方式作了详细分析说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims

1.一种基于稀疏变换学习的改进灵敏度编码重建方法，包括以下步骤：

其中，k是循环变量，表示数据的第k次迭代，上标“⁰”表示初始值；表示列向量化的待重建图像，N＝m×n，m和n分别是图像的行数和列数，x⁰表示x的初始值；/>表示列向量化的欠采样多线圈k空间数据，M为单线圈k空间数据实际采样的点数，且M＜＜N，L表示并行成像所使用的接收线圈个数，上标“^H”表示复共轭转置，y的第l个线圈数据用y_l表示，l＝1,...,L表示线圈索引，y₁和y_L分别表示列向量化的欠采样多线圈k空间数据y的第1个线圈数据和第L个线圈数据，/> 是对角矩阵，S_l对应第l个接收线圈的灵敏度，S₁表示第1个线圈的灵敏度图，S_L表示第L个线圈的灵敏度图；和/>分别为n×n和m×m的傅里叶变换矩阵，I_L为L×L的单位矩阵，/>表示克罗内克积；/>是一个欠采样矩阵，其中/>表示从单线圈的k空间网格中选择采样点位置的矩阵；

表示第j个图像块提取线性算子，从单个线圈图像提取大小为/>的图像块再列向量化成/>的列向量，P₁(·)表示第1个图像块提取线性算子，P_J(·)表示第J个图像块提取线性算子，/>表示/>图像块的稀疏变换矩阵，W⁰表示W的初始值，表示/>点离散余弦变换矩阵，/>表示WP_j(x)对应的辅助变量，将所有B_j按顺序水平拼接成矩阵/>B_j表示B的第j个分量，J表示图像块的个数，B₁表示B的第1分量，B_J表示B的第J个分量，B⁰表示B的初始值，中间变量/>表示从x中抽取J个/>的图像块并水平拼接成的矩阵，P₁(x)和P_J(x)分别表示从单线圈图像x中提取的第1个图像块和第J个图像块，X⁰＝[P₁(x⁰),...,P_J(x⁰)]表示X＝[P₁(x),...,P_J(x)]的初始值，/>是硬阈值函数，/> 表示输入矩阵，θ表示阈值，α,μ₂为参数，且/>为FSx对应的辅助变量，/>为z对应的拉格朗日乘子，/>表示u_z的初始值；

S1：对X^k(B^k)^H进行奇异值分解，得到X^k(B^k)^H＝UΣV^H；

其中，B^k是辅助变量B的第k次迭代，X^k是中间变量X的第k次迭代，U,Σ,V是奇异值分解因子；

W^k+1＝VU^H；

B^k+1＝H(W^k+1X^k,α/μ₂)；

S4：计算第k+1次迭代的辅助变量z^k+1，计算公式如下：

其中，x^k表示列向量化的待重建图像x的第k次迭代，R^T表示从k空间网格中选择采样点位置的矩阵的转置，表示拉格朗日乘子u_z的第k次迭代，μ₁＞0，(·)^-1表示求逆算子；

S5：计算第k+1次迭代的待重建图像x^k+1，计算公式如下：

其中，上标“^*”是伴随算子，是P_j(·)的伴随算子，/>是对角矩阵，其对角元素等于对应的像素点出现在所有图像块中的次数，μ₁＞0,μ₂＞0，/>是B_j的第k+1次迭代；

S6：计算第k+1次迭代的中间变量X^k+1，计算公式如下：

X^k+1＝[P₁(x^k+1),...,P_J(x^k+1)]；

S7：更新第k+1次迭代的拉格朗日乘子

S8：计算x^k+1和x^k之间的相对误差RE，计算式为：

S9：判断是否满足算法停止条件，若满足RE＜tol，则进入S10，否则令k＝k+1，返回S1；

S10：输出x＝x^k+1，得到重建图像x。