一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法与流程

本发明涉及一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，属于医学磁共振成像技术领域。

背景技术：

并行磁共振(magneticresonance，mr)成像是一种众所周知的加速成像方法。其优点是通过多线圈阵列接收器接收空间灵敏度信息来减少mr成像的采样时间。在过去的二十年里，许多并行成像重构方法被提出，这些方法因灵敏度信息的使用不同方式而有所不同；如：灵敏度编码(sensitivityencoding，sense)方法是利用显式灵敏度信息来进行重构的。这一方法在实际的应用中最主要的限制是很难准确的测量灵敏度信息。另一类是一般自校准部分并行采样(generalizedautocalibaratingpartiallyparallelacquisition，grappa)和迭代自一致性并行成像重构(iterativeself-consistentparallelimagingreconstruction，spirit)方法，它们利用隐式灵敏度信息进行重构，避免了难以准确测量灵敏度信息的限制。除此之外，还有一些不需要校准灵敏度信息的mr成像重构方法，如：无校准多线圈mr图像重构。但是相比于基于sense和spirit的重构方法，无校准的mr图像重构方法在使用相同正则项的情况下，重构效果通常不理想。

最近，一些研究人员利用稀疏正则项模型来改善grappa校准内核的重构质量。基于grappa自校准方法，spirit是一种在k空间网格中为每一个点(采样区和未采样区)执行校准一致性的重构方法。为了提高重构质量，一些研究者提出了联合l1(l1,2)正则项来促进多线圈图像在一些变换域中的联合稀疏性，并且在频率域或图像域中使用基于插值算子的凸集投影(projectionoverconvexsets，pocs)方法求解spirit重构问题。weller等引入联合总变分(jointtotalvariation，jtv)和l1,2正则项，提出了非笛卡尔spirit重构方法，并通过基于先验条件的变量分离(variablesplitting，vs)和交替方向乘子法(alternatingdirectionmethodofmultipliers，admm)进行求解。在前期工作中，提出了含jtv和l1,2正则项的笛卡尔spirit重构问题的快速求解方法，即：基于spirit的含联合稀疏的有效算子分裂方法(efficientoperatorsplittingalgorithmforjointsparsity-regularized，eosjs)。

为了进一步改善重构质量，一些研究工作者提出了基于稀疏变换(或字典)学习的重构方法。通过对参考图像进行全采样，提出了基于块的k-svd(singularvaluedecomposition，svd)字典的变换学习方法。然而，在一些新图像的新特征中，从原图像获取的字典学习块通常不是有效稀疏的，从而不能有效的对图像进行重构。为此，otaza等研究了一种在mr成像中利用k-svd的一维(onedimensional，1d)字典的重构方法，能有效的将学习字典块进行稀疏化，但是在图像的一些二维局部结构中，1d字典法也不能使学习的字典块有效的稀疏。ravishankar等提出了一种有效的字典学习mr成像(dictionarylearningmrimaging，dlmri)重构方法，它是通过从亚采样k空间方案和同时采样的mr图像中自适应学习的稀疏变换字典，该字典是利用k-svd方法得到的。为了更有效地重构mr成像，ravishankar等进一步提出了一种稀疏变换学习mr成像(transformlearningmrimaging，tlmri)重构方法。这种方法也是一种可以从训练数据中通过学习条件良好的平方稀疏变换的新重构方法。该方法可以更新学习的稀疏变换，从而得到问题的近似解。

相比于传统的固定稀疏变换的压缩感知mr成像，tlmri方法可以有效的提高图像的质量。另外，相比于字典是利用迭代k-svd方法求解的dlmri方法的合成字典法，tlmri方法是通过引入一个变换学习字典来对求出的近似解进行分析，从而能在单线圈mr成像重构中有效的节约采样时间。虽然，tlmri方法相比于其他方法能有效的提高图像的质量和处理速度，但该方法在多线圈mr成像中的图像质量和处理速度上还有一定的改善空间。

技术实现要素：

本发明要解决的技术问题是提供了一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，在于克服现有技术的不足，提高mr成像的质量。

本发明采用的技术方案是：一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，所述方法的具体步骤如下：

s0：初始化，令γ^-1＝((g-i)^h(g-i)+μ1i)^-1，x⁰＝a^hy，k＝1；

其中γ是一个具有对角结构的矩阵，能够直接进行求逆，g表示从自校准线中获得的图像域中的k空间自一致性卷积，i表示一个(n·nc)×(n·nc)的单位矩阵，μ1＞0是正则项参数，(·)^h表示共轭转置；是待重构的所有线圈图像，x＝(x1,...,xc,...,xnc)，表示x的第c列，也就是第c个线圈图像的向量化表示，n＝nx×ny表示单幅图像的像素点个数，nx和ny分别表示单幅图像的行数和列数，nc表示线圈图像的总个数，x⁰表示待重构的所有线圈图像的初始值；a是系统编码矩阵，即部分傅里叶变换矩阵，且是由n×n的单位矩阵in的m行构成的亚采样矩阵，m表示单位矩阵in中的m行，即采样数据的点数，是一个n×n的矩阵，分别是nx点和ny点的离散傅里叶变换矩阵，表示kronecker积；是所有线圈图像的部分傅里叶测量数据，y＝[y1,...,yc,...,ync]，表示y的第c列，也就是第c个线圈图像的部分傅里叶测量数据；按照重叠步幅为1且边界环绕方式，以单线圈二维图像的每个像素为左上顶点生成n个的重叠二维图像块，n表示二维图像块的总像素点数，1≤j≤n表示与重叠二维图像块相关的索引变量，表示第j个的二维图像块抽取矩阵，pj由n×n的单位矩阵in的n行构成，pj与的向量化单线圈图像相乘，则可抽取得到的二维图像块的n点列向量表示；pj与的向量化多线圈图像相乘，则可抽取得到nc个的二维图像块的列向量表示组成的n×nc矩阵，(·)^t表示矩阵转置；w表示图像块的方形稀疏变换，w^k表示第k次迭代的图像块的方形稀疏变换，k为迭代次数，w⁰表示图像块的方形稀疏变换的初始值；由于按照重叠步幅为1且边界环绕方式将二维图像分割为重叠二维图像块，故dc是对角矩阵，c表示线圈图像的索引变量；表示点离散余弦变换矩阵；为辅助变量z对应的对偶变量，表示uz的第c列，为待重构的所有线圈图像x的辅助变量，z＝(z1,...,zc,...,znc)，表示z的第c列；为辅助变量，b＝[b1,:,...,bj,:,...,bn,:]，表示b的第j个分块，表示bj,:的第c列；为辅助变量b对应的对偶变量，ub＝[(ub)1,:,...,(ub)j,:,...,(ub)n,:]，表示ub的第j个分块，表示(ub)j,:的第c列；为辅助变量z对应的对偶变量的初始值；为辅助变量b对应的对偶变量的初始值；

s1：计算辅助变量z^k，计算公式如下：

其中z^k表示第k次迭代的待重构的所有线圈图像x^k的辅助变量，x^k-1表示第k-1次迭代的待重构的所有线圈图像，是第k-1次迭代的辅助变量z^k-1的对偶变量；vec(·)将矩阵操作数的所有列堆积成一个列向量，unvec(·)是将vec(·)转化成的列向量转化成原来形式的矩阵；

s2：初始化，j＝1；

s3：计算辅助变量的第j个分块公式如下：

式中，w^k-1表示第k-1次迭代的图像块的方形稀疏变换，表示从第k-1次迭代的待重构的所有线圈图像x^k-1中抽取nc个图像块并组成n×nc的矩阵；表示第k次迭代的辅助变量b^k的第j个分块，与相对应，表示第k-1次迭代的辅助变量b^k-1的对偶变量，表示第k-1次迭代的对偶变量的第j个分块，μ2＞0是正则项参数；hj(b,θ)是硬阈值处理函数，其第一个输入参数b＝[b1,...,bc,...,bnc]，为b的第c列，第二个参数θ＞0为一个常数，hj(b,θ)定义为：

其中，表示矩阵的逐点乘积，为一个中间变量，计算公式为：repmat(s,[1,nc])是matlab中的函数，功能为将第一个参数s按照第二个参数[1,nc]的维度进行复制拓展，拓展之后得到的矩阵；

s4：判断是否完成b的所有图像块分块的更新，若是，则进入步骤s5；否则，对j进行加一操作后返回s3；

s5：更新第k次迭代的图像块的方形稀疏变换w^k，计算公式如下：

其中：in表示n×n的单位矩阵，是一个分解因子，l＝(x^k-1(x^k-1)^h+0.5λin)^1/2，是一个中间变量矩阵，其值为x^k-1＝[p1x^k-1,...,pjx^k-1,...,pnx^k-1]，x^k-1表示第k-1次迭代的待重构所有线圈图像，且的完整奇异值分解为uσv^h，是奇异值分解因子，b^k表示第k次迭代的辅助变量，λ＞0表示正则项参数；

s6：初始化，c＝1；

s7：计算待重构的第c个线圈图像公式如下：

其中：表示第k次迭代的待重构的第c个线圈图像，表示第k次迭代的辅助变量z^k的第c列，表示第k-1次迭代的对偶变量的第c列，表示第k次迭代的辅助变量b^k的第j个分块的第c列，表示第k-1次迭代的对偶变量的第j个分块的第c列，α＞0表示正则项参数；

s8：判断是否完成所有线圈，若是，进入步骤s9；否则，对c进行加一操作后返回s7；

s9：计算对偶变量公式如下：

其中，是第k次迭代的辅助变量z^k对应的对偶变量，x^k表示第k次迭代的待重构的所有线圈图像；

s10：计算对偶变量公式如下：

其中，表示第k次迭代的辅助变量b^k对应的对偶变量，是一个中间变量矩阵，其值为x^k＝[p1x^k,...,pjx^k,...,pnx^k]，表示第k次迭代的待重构的所有线圈图像x^k中抽取nc个图像块并组成n×nc的矩阵；

s11：判断是否达到最大迭代次数k，k是最大迭代次数，若是，输出重构的所有线圈图像x；否则，对k进行加一操作后返回步骤s1；

s12：得到重构的所有线圈图像x后，使用平方和的平方根方法来形成幅度图像，计算公式如下：

本发明的有益效果是：在本发明中，提出了一种将自适应稀疏变换学习和联合稀疏正则项与笛卡尔spirit并行mr成像相结合的方法，以提高重建质量。通过使用vs和admm技术，可以将重建问题转换为具有近似形式解的一些子问题。对两个活体数据集的实验仿真表明，所提出的方法优于其他比较的方法，能有效的提高重构图像的视觉质量。

附图说明

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

图2在data1下，使用admmjtvli方法重构的图像。

图3在data1下，使用dlspirit方法重构的图像。

图4在data1下，使用jtlspirit方法重构的图像。

图5在data1下，使用admmjtvli方法重构的误差图像。

图6在data1下，使用dlspirit方法重构的误差图像。

图7在data1下，使用jtlspirit方法重构的误差图像。

图8在data2下，使用admmjtvli方法重构的图像。

图9在data2下，使用dlspirit方法重构的图像。

图10在data2下，使用jtlspirit方法重构的图像。

图11在data2下，使用admmjtvli方法重构的误差图像。

图12在data2下，使用dlspirit方法重构的误差图像。

图13在data2下，使用jtlspirit方法重构的误差图像。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案：一种基于变换学习和联合稀疏性的迭代自一致性并行成像重构方法，所述方法的具体步骤如下：

假设表示待重构的所有线圈图像，x＝(x1,...,xc,...,xnc)，表示x的第c列，也就是第c个线圈图像的向量化表示，n＝nx×ny表示单幅图像的像素点个数，nx和ny分别表示单幅图像的行数和列数，nc表示线圈图像的总个数。所有线圈图像的部分傅里叶测量数据或k空间亚采样数据由下式给定：

y＝ax(1)

其中，表示y的第c列，也就是第c个线圈图像的部分傅里叶测量数据；a是一个部分傅里叶变换矩阵，且是由n×n的单位矩阵in的m行构成的亚采样矩阵，m表示单位矩阵in中的m行，即采样数据的点数。是一个n×n的矩阵，分别是nx点和ny点的离散傅里叶变换矩阵，表示kronecker积。

spirit内核一致性项是由通过在图像域中的矩阵乘积和相应k空间中的卷积块的基形成的。在图像域中的spirit校准一致性方程如下所示：

vec(x)＝g·vec(x)(2)

其中g表示从自校准线中获得的图像域中的k空间自一致性卷积，vec(·)将矩阵操作数的所有列堆积成一个列向量，unvec(·)是将vec(·)转化成的列向量转化成原来形式的矩阵。

考虑到线圈图像之间小波变换域的相似性，murphy等引入了联合稀疏来提高重建质量。因此，并行成像重建问题可以表述为以下最小化问题：

其中ψ是逐线圈小波变换，且||ω||1,2定义为：

其中ω表示一个中间变量，ωrc表示第c个线圈的第r个k空间位置的中间变量。c是线圈图像索引变量，r是k空间位置索引变量。

murphy等提出了一种有效的pocs方法求解问题(3)。lustig等和vasanawala等利用了k空间自一致性卷积，也通过使用pocs方法来解决类似问题。为了提高并行mr图像的质量，weller等提出了利用jtv和l1,2正则项来代替问题(3)中的l1,2正则项来进行图像重构，相应的无约束问题如下：

其中表示frobenius范数，||·||jtv表示联合总变分正则项。i表示一个(n·nc)×(n·nc)的单位矩阵，α1，α2和γ均是正则化参数。weller等提出了利用vs和admm方法求解问题(5)。

为了进一步改善图像重构的质量，本发明提出了一种引入多线圈图像的稀疏变换学习(transformlearning，tl)和联合稀疏正则项模型，将重构问题转换为：

其中为中间变量，表示w的第c列。表示线圈图像的方形稀疏变换，按照重叠步幅为1且边界环绕方式，以单线圈二维图像的每个像素为左上顶点生成n个的重叠二维图像块，n表示二维图像块的总像素点数，1≤j≤n表示与重叠二维图像块相关的索引变量，表示第j个的二维图像块抽取矩阵，pj由n×n的单位矩阵in的n行构成，pj与的向量化单线圈图像相乘，则可抽取得到的二维图像块的n点列向量表示；pj与的向量化多线圈图像相乘，则可抽取得到nc个的二维图像块的列向量表示组成的n×nc矩阵，表示从待重构的所有线圈图像x中抽取nc个图像块并组成n×nc的矩阵。函数是w变换的一个正则项。与之前的变换学习的方法类似，定义在目标函数中作为正则项函数，det(·)表示矩阵的行列式。α，λ均表示正则项参数。

引入辅助变量z＝x和b＝[b1,:,...,bj,:,...,bn,:]＝[wp1x,...,wpjx,...,wpnx]，以及相应的对偶变量uz和ub。表示从多线圈图像x抽取nc个图像块组成的n×nc的矩阵的方形稀疏变换，表示b的第j个分块，表示bj,:的第c列。利用admm求解问题(6)，子问题的迭代解如下：

在(7)-(12)中，为待重构的所有线圈图像x的辅助变量，z＝(z1,...,zc,...,znc)，表示z的第c列，z^k表示第k次迭代的所有线圈图像x的辅助变量，k为迭代次数；为辅助变量z对应的对偶变量，表示uz的第c列，是第k-1次迭代的辅助变量z^k-1对应的对偶变量；x^k-1表示第k-1次迭代的待重构的所有线圈图像。表示第k次迭代的辅助变量b^k的第j个分块，与相对应；w^k-1表示第k-1次迭代的线圈图像的方形稀疏变换；表示第k-1次迭代的待重构的所有线圈图像x^k-1中抽取nc个图像块并组成n×nc的矩阵；为辅助变量b对应的对偶变量，表示ub的第j个分块，表示的第c列，表示第k-1次迭代的对偶变量的第j个分块。w^k表示第k次迭代的线圈图像的方形稀疏变换。是一个中间变量矩阵，其值为b^k表示第k次迭代的辅助变量。表示第k次迭代的待重构的第c个线圈图像，表示第k次迭代的辅助变量z^k的第c列，表示第k-1次迭代的对偶变量的第c列，表示从向量化的第c个线圈图像xc中抽取一个向量化为n点的图像块，表示第k次迭代的辅助变量b^k的第j个分块的第c列，表示第k-1次迭代的对偶变量的第j个分块的第c列。是第k次迭代的辅助变量z^k对应的对偶变量，x^k表示第k次迭代的待重构所有线圈图像。表示第k次迭代的辅助变量b^k对应的对偶变量，是一个中间变量矩阵，其值为表示第k次迭代的待重构的所有线圈图像x^k中抽取nc个图像块并组成n×nc的矩阵。μ1，μ2均表示大于零的正则项参数。

设γ＝(g-i)^h(g-i)+μ1i，子问题(7)可以通过下式进行求解：

式中，γ是一个具有对角结构的矩阵，(·)^h表示共轭转置，由于spirit一致性算子(g-i)的对角块结构仅在线圈之间耦合，因此可以直接对γ中的每一个nc×nc块进行求逆。

问题(8)的解如下：

其中hj(b,θ)是硬阈值处理函数，其第一个输入参数为b的第c列，第二个参数θ＞0为一个常数，hj(b,θ)定义为：

其中，表示矩阵的逐点乘积，为一个中间变量，计算公式为：repmat(s,[1,nc])是matlab中的函数，功能为将第一个参数s按照第二个参数[1,nc]的维度进行复制拓展，拓展之后得到的矩阵。

子问题(9)可以由下式求解：

其中in表示n×n的单位矩阵，是一个分解因子，l＝(x^k-1(x^k-1)^h+0.5λin)^1/2，且具有uσv^h的全部奇异值分解(singularvaluedecomposition，svd)，是奇异值分解因子。

由于按照重叠步幅为1且边界环绕方式将二维图像分割为重叠二维图像块，故是一个对角矩阵，其中(·)^t表示矩阵转置。因此，子问题(10)能有效的求解：

现在，得到了子问题(7)-(10)的求解方法(如公式(13)-(17)所示)，一个基于联合稀疏和变换学习的迭代自一致性并行mr成像重建方法(jointsparsityandtransformlearningbasedspiritparallelmrimagingreconstructionalgorithm，jtlspirit)，具体流程如图1所示，其中步骤如下：

s0：初始化，令γ^-1＝((g-i)^h(g-i)+μ1i)^-1，x⁰＝a^hy，

x⁰表示待重构的所有线圈图像的初始值；w⁰表示图像块的方形稀疏变换的初始值；表示离散余弦变换矩阵；为辅助变量z对应的对偶变量的初始值；为辅助变量b对应的对偶变量的初始值。

s1：计算辅助变量z^k，计算公式如(13)；

s2：初始化，j＝1；

s3：计算辅助变量的第j个分块公式如(14)；

s4：判断是否完成b的所有图像块分块更新，若是，则进入步骤s5；否则，对j进行加一操作后返回s3；

s5：更新第k次迭代的图像块的方形稀疏变换w^k，计算公式如(16)；

s6：初始化，c＝1；

s7：计算待重构的第c个线圈图像公式如(17)；

s8：判断是否完成所有线圈，若是，进入步骤s9；否则，对c进行加一操作后返回s7；

s9：计算对偶变量公式如(11)；

s10：计算对偶变量公式如(12)；

s11：判断是否达到最大迭代次数k，k是最大迭代次数，若是，输出重构的所有线圈图像x；否则，对k进行加一操作后返回步骤s1；

s12：得到重构的所有线圈图像x后，使用平方和的平方根方法来形成幅度图像，计算公式如下：

在以下实验中，比较了所提出的方法jtlsplrit与使用jtv和l1正则项的spirit并行成像重构方法(admmjtvl1)和基于k-svd字典学习的spirit并行成像重构方法(dlspirit)。所有比较的方法都在matlab中实现。

本发明使用了两个由八通道线圈获得的完整视场(fieldofview，fov)的体内大脑数据集，即data1和data2。为了生成测试数据集，采用加速因子为r×(不包括自校准采样(self-calibratedsampling，acs)线)和24×24acs线的笛卡尔泊松盘亚采样模式对完全采样的数据集进行人工亚采样。此外，在所有比较方法的实验中使用5×5spirit校准内核，并且从完全采样数据集中使用公式(18)计算参考图像。

对于dlspirit，使用大小为6×6(n＝36)的图像块，并且利用从14400个随机选择块中用k-svd迭代方法(10次迭代)来学习1.5倍的超完备合成词典。对于jtlspirit，使用6×6(n＝36)的图像块。手动调整方法的正则化参数以获得更好的重建质量。当x^k和x^k-1之间的相对误差低于公差tol＝2e-4时，所有比较的方法都被终止。本发明所有的实验都是在配备intelcorei7-7500u@2.7ghz处理器，8gb内存和windows10操作系统(64位)的笔记本电脑上进行的。

图2-13展示了采用加速因子为6×和24×24acs线的泊松盘亚采样模式，在数据集data1和data2下所有比较方法的并行mr成像的重构图像及其相应的误差图。在图2中，可以明显看出admmjtvl1方法的重构的图像出现了较多的伪影。在图3中，dlspirit方法可以稍微缓解图像重构的伪影，但在过平滑的区域中，其重构的图像存在一些细节上的缺失。相比之下，在图4中，提出的jtlspirit方法有效的抑制了图像重构中的伪像，并且保留了其他方法中缺失的细节。图5-7分别是admmjtvli，dlspirit和jtlspirit的重建误差图像，从图7中可以看出，提出的jtlspirit重构图像的误差算子较小(黑色的部分越多，误差算子越小)，能够重构出图像的更多细节。因此jtlspirit在所有比较方法中能得到最佳的视觉质量。对于图8-13，也可以得出同样的结论。

在实验仿真中，还对所有比较的方法和本发明提出的方法的运行时间进行了比较。对于data1，admmjtvl1,dlspirit和jtlspirit的平均运行时间分别为20秒，1600秒和260秒。对于data2，admmjtvl1,dlspirit和jtlspirit的平均运行时间分别为68秒，1950秒和410秒。综上可以看出，与admmjtvl1相比，提出的jtlspirit算法尽管需要更多的重构时间，但重构质量有较大的提高。而与dlspirit相比，jtlspirit重构图像的snr有所增加，视觉质量有较大改善，而且重构速度是dlspirit的5倍以上。

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