本发明公开了一种基于张量分解稀疏约束的三维心脏磁共振成像方法，其采用三维径向采样轨迹实现全心脏三维K空间数据的欠采样；通过高阶张量分解实现磁共振图像的稀疏表示，实现全心脏三维磁共振数据的最优稀疏表示，提高磁共振成像的精度；结合张量分解的稀疏数据的L1范数与全变差变换的复合正则化项作为约束项，利用快速复合分裂算法实现欠采样三维心脏磁共振图像重构。本发明方法有效缩短磁共振成像扫描时间，提高心脏磁共振成像的速度，有利于消除心脏检测过程中产生的运动伪影，提高三维心脏磁共振成像的精度。
1.一种基于张量分解稀疏约束的三维心脏磁共振成像方法，包括如下步骤：(1)利用三维径向采样轨迹模式对人体心脏的磁共振K空间数据进行欠采样，得到心脏的K空间欠采样数据；所述的K空间欠采样数据分为N_z个采样层面，每个采样层面包含有N_p条投影线，每条投影线上包含有N_s个样本点；N_z、N_p和N_s均为大于1的自然数；所述的K空间欠采样数据中各样本点的三维坐标表示如下： G z ( p , s , k ) = 2 ( k - 1 ) N z - 1 - 1 ]]> G x ( p , s , k ) = p h c o s ( s θ + 2 k π N z ) G y ( p , s , k ) = p h s i n ( sθ + 2 k π N z ) ]]> θ = 2 π N s h = R N p ]]> 其中：G_x(p,s,k)、G_y(p,s,k)和G_z(p,s,k)分别为K空间欠采样数据中第k个采样层面的第p条投影线上第s个样本点在x轴、y轴和z轴上的坐标，R为磁共振K空间数据的最大径向半径，k、p和s均为自然数且1≤k≤N_z，1≤p≤N_p，1≤s≤N_s；(2)对所述的K空间欠采样数据进行傅里叶逆变换，得到初始三维心脏磁共振图像u₀；(3)利用张量分解的稀疏正则项结合全变差的稀疏正则项作为约束，建立重构三维心脏磁共振图像的目标函数如下： m i n u { | | A u - F | | 2 2 + α T V ( u ) + β H O S V D ( u ) } ]]> 其中：u为三维心脏磁共振图像，A为系统矩阵，F为K空间欠采样数据，||||₂为2-范数，HOSVD(u)表示对三维心脏磁共振图像u张量分解的稀疏正则项，TV(u)表示对三维心脏磁共振图像u全变差变换的稀疏正则项，α和β均为权重参数；所述的稀疏正则项HOSVD(u)的表达式如下： H O S V D ( u ) = u × 1 U 1 T × 2 U 2 T × 3 U 3 T ]]> 其中：×_n表示n-模乘法，U_n为左奇异向量且通过对三维心脏磁共振图像u进行n模矩阵化奇异值分解得到，n＝1、2或3，^T表示转置；(4)根据初始三维心脏磁共振图像u₀通过以下迭代方程组对上述目标函数进行最小化求解，重建得到三维心脏磁共振图像u； x t = r t - ▿ { f ( r t ) } f ( r t ) = 1 2 | | Ar t - F | | 2 2 ]]> a t = m i n u { α T V ( u ) + 1 2 | | u - x t | | 2 } b t = m i n u { β H O S V D ( u ) + 1 2 | | u - x t | | 2 } ]]> u t = 1 2 ( a t + b t ) ]]> w t + 1 = 1 + 1 + 4 w t 2 2 r t + 1 = u t + ( w t - 1 w t + 1 ) ( u t - u t - 1 ) ]]> 其中：u_t和u_t-1分别为第t次迭代和第t-1次迭代的三维心脏磁共振图像，▽{}表示微分算子，a_t和b_t均为第t次迭代的重构子图像，r₁＝u₀，w₁＝1，t为迭代次数。
2.根据权利要求1所述的三维心脏磁共振成像方法，其特征在于：所述的稀疏正则项TV(u)的表达式如下： T V ( u ) = Σ i = 1 N x Σ j = 1 N y Σ k = 1 N z D x , i j k ( u ) 2 + D y , i j k ( u ) 2 + D z , i j k ( u ) 2 ]]> D x , i j k ( u ) = u i + 1 , j , k - u i , j , k i < N x 0 i = N x ]]> D y , i j k ( u ) = u i , j + 1 , k - u i , j , k j < N y 0 j = N y ]]> D z , i j k ( u ) = u i , j , k + 1 - u i , j , k k < N z 0 k = N z ]]> 其中：N_x、N_y和N_z分别为三维心脏磁共振图像u在x轴、y轴和z轴上的维度，u_i,j,k为三维心脏磁共振图像u的第k层平面图像中第i行第j列像素的像素值，u_i+1,j,k为三维心脏磁共振图像u的第k层平面图像中第i+1行第j列像素的像素值，u_i,j+1,k为三维心脏磁共振图像u的第k层平面图像中第i行第j+1列像素的像素值，u_i,j,k+1为三维心脏磁共振图像u的第k+1层平面图像中第i行第j列像素的像素值，i、j和k均为自然数且1≤i≤N_x，1≤j≤N_y，1≤k≤N_z。
3.根据权利要求1所述的三维心脏磁共振成像方法，其特征在于：所述的重构子图像a_t和b_t通过收缩阈值化的快速迭代算法求解得到。
4.根据权利要求1所述的三维心脏磁共振成像方法，其特征在于：通过所述的迭代方程组进行迭代计算，使达到最大迭代次数或迭代收敛后的三维心脏磁共振图像即作为重建得到的三维心脏磁共振图像u；迭代收敛条件如下： | | u t - u t - 1 | | 2 | | u t | | 2 < 10 - 5 . ]]> 技术领域
本发明属于磁共振成像技术领域，具体涉及一种基于张量分解稀疏约束的三维心脏磁共振成像方法。
背景技术
心血管疾病是导致心脏突然停止跳动的主要原因。当前心血管病发病和死亡率居高不下，使心血管病防治负担加重，成为重要公共卫生问题，加强心血管病防治刻不容缓。心脏磁共振成像(CardiacMagneticResonanceImaging，CMR)是利用核磁共振原理进行人体心脏断层成像的技术，能准确地反映心脏的解剖结构、形态功能、血流特性和心肌活性，已迅速发展成为心脏疾病诊断中的主要工具。心脏磁共振成像在成像过程中具有较好的软组织对比度，没有任何放射性污染，分辨率高，可任意层面成像；而且由于参与磁共振成像的因素较多，得到的图像信息量大，优于现有的其它各种影像学成像技术，在心脏疾病诊断中有很大的优越性和应用潜力。在早中期的磁共振成像中，设备扫描需时较长，如一次心脏扫描需要大约数小时左右，甚至更长时间，这限制了心脏磁共振成像的应用范围。磁共振成像大多被用在静态部位的成像中，通常不适用于动态成像，这是因为在动态成像中，数据量的增加会大幅度延长数据获取时间，需时会更长，病人常常会因为长时间的等待数据扫描感觉不适；或在数据获取过程中发生自主和非自主的运动，导致图像中出现各种伪影。提高磁共振成像速度，缩短磁共振成像扫描时间的意义不仅仅在于提高磁共振设备的工作效率，减轻病人痛苦，更重要的是它有利于消除心脏运动以及呼吸等造成运动伪影的影响。目前所研究的磁共振图像重构方法主要两种：一种是多线圈并行成像技术，主要是利用相控阵线圈中单个接收线圈的空间敏感度差异来编码空间信息，降低成像所必需的梯度编码步数，其采用多线圈阵列同时采集信号，允许对K空间进行欠采样以减少相位编码步数，在保持图像空间分辨率不变的同时，能大幅度缩短扫描时间，提高成像速度；但多线圈并行成像技术涉及多线圈敏感度分布的估计，需要增加计算量。另一种基于压缩感知理论的磁共振成像重构方法，由于磁共振影像具有稀疏特性，可以采用压缩感知理论从随机欠采样的k空间数据进行图像重构，减少采样数据，提高成像速度，目前常用的奇异值分解、离散小波变换、离散余弦变换等稀疏变换方法，只考虑单层心脏磁共振图像的稀疏表示，而没有考虑三维心脏磁共振图像层与层之间的稀疏性。
发明内容
针对现有技术所存在的上述技术问题，本发明提供了一种基于张量分解稀疏约束的三维心脏磁共振成像方法，不仅考虑心脏磁共振图像层内的稀疏特性，也考虑心脏磁共振成像层间的稀疏特性，可提高磁共振成像精度。一种基于张量分解稀疏约束的三维心脏磁共振成像方法，包括如下步骤：(1)利用三维径向采样轨迹模式对人体心脏的磁共振K空间数据进行欠采样，得到心脏的K空间欠采样数据；(2)对所述的K空间欠采样数据进行傅里叶逆变换，得到初始三维心脏磁共振图像u₀；(3)利用张量分解的稀疏正则项结合全变差的稀疏正则项作为约束，建立重构三维心脏磁共振图像的目标函数如下： min u { | | Au - F | | 2 2 + αTV ( u ) + βHOSVD ( u ) } ]]> 其中：u为三维心脏磁共振图像，A为系统矩阵(该矩阵由傅立叶变换矩阵和K空间数据欠采样掩模构成)，F为K空间欠采样数据，||||₂为2-范数，HOSVD(u)表示对三维心脏磁共振图像u张量分解的稀疏正则项，TV(u)表示对三维心脏磁共振图像u全变差变换的稀疏正则项，α和β均为权重参数；(4)根据初始三维心脏磁共振图像u₀对上述目标函数进行最小化求解，重建得到三维心脏磁共振图像u。所述的K空间欠采样数据分为N_z个采样层面，每个采样层面包含有N_p条投影线，每条投影线上包含有N_s个样本点；N_z、N_p和N_s均为大于1的自然数。所述的K空间欠采样数据中各样本点的三维坐标表示如下： G z ( p , s , k ) = 2 ( k - 1 ) N z - 1 - 1 ]]> G x ( p , s , k ) = ph cos ( sθ + 2 kπ N z ) , G y ( p , s , k ) = ph sin ( sθ + 2 kπ N z ) ]]> θ = 2 π N s , h = R N p ]]> 其中：G_x(p,s,k)、G_y(p,s,k)和G_z(p,s,k)分别为K空间欠采样数据中第k个采样层面的第p条投影线上第s个样本点在x轴、y轴和z轴上的坐标，R为磁共振K空间数据的最大径向半径，k、p和s均为自然数且1≤k≤N_z，1≤p≤N_p，1≤s≤N_s。所述的稀疏正则项HOSVD(u)的表达式如下： HOSVD ( u ) = u × 1 U 1 T × 2 U 2 T × 3 U 3 T ]]> 其中：×_n表示n-模乘法，U_n为左奇异向量且通过对三维心脏磁共振图像u进行n模矩阵化奇异值分解得到，n＝1、2或3，^T表示转置。所述的稀疏正则项TV(u)的表达式如下： TV ( u ) = Σ i = 1 N x Σ j = 1 N y Σ k = 1 N z D x , ijk ( u ) 2 + D y , ijk ( u ) 2 + D z , ijk ( u ) 2 ]]> D x , ijk ( u ) = u i + 1 , j , k - u i , j , k i < N x 0 i = N x ]]> D y , ijk ( u ) = u i , j + 1 , k - u i , j , k j < N y 0 j = N y ]]> D z , ijk ( u ) = u i , j , k + 1 - u i , j , k k < N z 0 k = N z ]]> 其中：N_x、N_y和N_z分别为三维心脏磁共振图像u在x轴、y轴和z轴上的维度，u_i,j,k为三维心脏磁共振图像u的第k层平面图像中第i行第j列像素的像素值，u_i+1,j,k为三维心脏磁共振图像u的第k层平面图像中第i+1行第j列像素的像素值，u_i,j+1,k为三维心脏磁共振图像u的第k层平面图像中第i行第j+1列像素的像素值，u_i,j,k+1为三维心脏磁共振图像u的第k+1层平面图像中第i行第j列像素的像素值，i、j和k均为自然数且1≤i≤N_x，1≤j≤N_y，1≤k≤N_z。所述的步骤(4)中通过以下迭代方程组对目标函数进行最小化求解： x t = r t - ▿ { f ( r t ) } , f ( r t ) = 1 2 | | Ar t - F | | 2 2 ]]> a t = min u { α · TV ( u ) + 1 2 | | u - x t | | 2 } , b t = min u { β · HOSVD ( u ) + 1 2 | | u - x t | | 2 } ]]> u t = 1 2 ( a t + b t ) ]]> w t + 1 = 1 + 1 + 4 w t 2 2 , r t + 1 = u t + ( w t - 1 w t + 1 ) ( u t - u t - 1 ) ]]> 其中：u_t和u_t-1分别为第t次迭代和第t-1次迭代的三维心脏磁共振图像，ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="109"/>表示微分算子，a_t和b_t均为第t次迭代的重构子图像，r₁＝u₀，w₁＝1，t为迭代次数。所述的重构子图像a_t和b_t通过收缩阈值化的快速迭代算法求解得到。通过所述的迭代方程组进行迭代计算，使达到最大迭代次数或迭代收敛后的三维心脏磁共振图像即作为重建得到的三维心脏磁共振图像u；迭代收敛条件如下： | | u t - u t - 1 | | 2 | | u t | | 2 < 10 - 5 ]]> 本发明基于张量分解稀疏约束的三维心脏磁共振成像方法为三维心脏磁共振成像研究探索出新的方法，其基于三维径向采样模式实现心脏三维K空间数据的欠采样，可实现全心脏数据径向欠采样，加快磁共振数据采集速度，降低磁共振设备的扫描时间；其基于高阶奇异值分解的全心脏磁共振数据的稀疏表示，实现心脏磁共振采样数据的最优稀疏，提高磁共振成像的精度；其基于快速复合分裂算法的三维心脏磁共振图像重建算法研究，提高磁共振成像重构的速度。
附图说明
图1为三维径向欠采样轨迹模式的示意图。图2为张量分解的n-模矩阵化示意图。图3是本发明利用高阶奇异值分解实现三维心脏数据稀疏过程的示意图。图4为本发明三维心脏磁共振成像方法的流程示意图。
具体实施方式
为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。如图4所示，本实施方式的基于张量分解稀疏约束的三维心脏磁共振成像方法，具体实施步骤如下：(1)利用三维径向采样轨迹实现K空间数据的欠采样。一个三维的径向采样轨迹，其中含有N_z个采样层面，每个层面包含N_p条投影线，每条投影线上包含有N_s个样本点。每个采样层面都是由第一个层面围绕k_z轴作一定旋转得到，如图1所示，三维径向采样轨迹的采样点三维坐标(G_x,G_y,G_z)定义如下： G z ( p , s , k ) = 2 ( k - 1 ) N z - 1 - 1 ]]> G x ( p , s , k ) = ph cos ( sθ + 2 kπ N z ) , G y ( p , s , k ) = ph sin ( sθ + 2 kπ N z ) ]]> θ = 2 π N s , h = R N p ]]> 其中：R为磁共振K空间数据的最大径向半径；p为第p条投影线，p＝1,2,…,N_p；k为第k个层面，k＝1,2,…,N_z；s为第s个样本点，s＝1,2,…,N_s。三维径向采样轨迹的欠采样率由层面数目N_z、投影线数目N_p和样本点数目N_s确定。本实施方式设定N_z＝25，N_p＝64，N_s＝255。(2)利用张量分解实现三维心脏磁共振数据的稀疏表示。高阶奇异值分解是矩阵奇异值分解在高阶张量中的推广，可以将高阶张量分解成核心张量和各子空间矩阵的积。将心脏数据作为三阶张量ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="271"/>，如图2所示，对三维张量u按模型n矩阵化分别得到u(n)，再对u(n)进行奇异值分解，即：[U_n,S_n,V_n]＝SVD(u(n))n＝1,2,3式中：U_n为其左奇异向量且为I_n×I_n的正交子空间矩阵；利用高阶奇异值分解(HighOrderSingularValueDecomposition，HOSVD)，可以通过下列公式计算得到核心张量Z： Z = HOSVD ( u ) = u × 1 U 1 T × 2 U 2 T × 3 U 3 T ]]> 式中：×_n是n-模乘法，可以被看作是先将三维张量u按模型n矩阵化分别得到u(n)，再进行普通的矩阵相乘，然后将其重新排列成张量。通过高阶奇异值分解，可以将全心脏的三维数据张量u转换成一个三维的稀疏张量Z，如图3所示。(3)采用线性组合稀疏变换L1正则化项约束和全变差(TotalVariation，TV)稀疏项约束，也就是同时考虑磁共振图像的某一稀疏变换和全变差变换作为图像的稀疏表示，以提高磁共振成像的精度，通过三维全变差变换与高阶奇异值分解相结合实现三维心脏磁共振成像的稀疏表示，三维全变差变换定义如下： TV ( u ) = Σ i = 1 N x Σ j = 1 N y Σ k = 1 N z D x , ijk ( u ) 2 + D y , ijk ( u ) 2 + D z , ijk ( u ) 2 ]]> D x , ijk ( u ) = u i + 1 , j , k - u i , j , k i < N x 0 i = N x ]]> D y , ijk ( u ) = u i , j + 1 , k - u i , j , k j < N y 0 j = N y ]]> D z , ijk ( u ) = u i , j , k + 1 - u i , j , k k < N z 0 k = N z ]]> 式中，u为三维心脏磁共振图像数据，N_x、N_y、N_z分别为x，y，z轴的最大范围值。(4)三维心脏磁共振成像重构算法用公式表示如下： min u { | | Au - F | | 2 2 + αTV ( u ) + βHOSVD ( u ) } ]]> 其中，F表示K空间欠采样数据，A表示系统矩阵，式中设定α＝0.01和β＝0.005分别为全变差稀疏正则项和高阶奇异值分解稀疏正则项的正则化参数。将以上重构图像的目标函数分解成为两个简单的正则化子问题，分别为TV和HOSVD，定义如下： min u { f ( u ) + α · TV ( u ) } ]]> min u { f ( u ) + β · HOSVD ( u ) } ]]> 式中： f ( u ) = 1 2 | | Au - F | | 2 2 . ]]>迭代计算的详细过程如下：首先，设定初始状态；将欠采样的K空间数据直接进行傅里叶逆变换得到的图象u₀作为初始图像；r₁＝u₀，w₁＝1。然后，根据快速迭代的收缩-阈值化算法求解迭代子问题的解a_t和b_t，并线性组合得到重构图像u_t： x t = r t - ▿ { f ( r t ) } , f ( r t ) = 1 2 | | Ar t - F | | 2 2 ]]> a t = min u { α · TV ( u ) + 1 2 | | u - x t | | 2 } , b t = min u { β · HOSVD ( u ) + 1 2 | | u - x t | | 2 } ]]> u t = 1 2 ( a t + b t ) ]]> w t + 1 = 1 + 1 + 4 w t 2 2 , r t + 1 = u t + ( w t - 1 w t + 1 ) ( u t - u t - 1 ) ]]> 最后，判断迭代收敛，如果前后两次重构误差范围满足以下收敛条件，则停止迭代计算；否则重复迭代，直至达到最大迭代次数200为止。 | | u t - u t - 1 | | 2 | | u t | | 2 < 10 - 5 ]]> 以下我们采用3T飞利浦磁共振扫描仪全采样心脏电影成像的K空间数据，获取数据所用成像序列是稳态自由进动(SteadyStateFreePrecession，SSFP)序列，采用心脏门控技术。磁共振扫描参数TR/TE为3.95/1.97ms，翻转角度为60°，视野(FieldOfView，FOV)为256×256×10mm³，空间分辨率为1.3×1.3×1.3mm³。在该项目研究中，采用全采样K空间采样点重建心脏电影的幅度图像，用来作为衡量重建图像质量的金标准u₀。在三维径向欠采样模式中，设置K空间数间欠采样率20％，根据高阶奇异值分解和三维全变差方法实现欠采样磁共振数据的稀疏表示，利用迭代计算实现全心脏磁共振图像的重构。该项目所研究的三维心脏磁共振成像方法将与小波变换的稀疏表示方法重构的磁共振图像进行比较，将从多个方面评价磁共振重构图像的质量。首先利用信噪比、相对误差、重构时间三个定量指标评价重构算法的性能；其次对重建图像、以及误差图(重建图像和原始图像之间差的绝对值)从视觉角度评价重构质量。重构图像的信噪比(Signal-to-NoiseRatio，SNR)和相对误差(RelativeError，RE)定义如下： SNR = 10 log 10 Var ( u r ) MSE ( u - u r ) , RE = | | u - u r | | 2 | | u r | | 2 ]]> 式中：u为重构的三维心脏磁共振图像，u_r为参考标准图像，MSE(u-u_r)为重构图像与参考标准图像的均方误差，Var(u_r)为参考标准图像u_r的方差值。基于张量分解和小波分解稀疏变换重构的心脏磁共振图像结果及其性能分别见表1；从表中结果可以看出本实施方式相对现有基于小波分解稀疏变换的重构方法，从磁共振成像的精度和速度上都有显著提高。表1 ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/> ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/>