一种动态磁共振成像重建方法与流程

本发明涉及动态磁共振成像，尤其是涉及同时利用三维动态磁共振图像时间维和空间维冗余信息并且可区分时间维稀疏和空间维稀疏各自重要性的三维稀疏约束重建来加快高空间时间分辨率的一种动态磁共振成像重建方法。

背景技术：

动态磁共振成像在医学临床诊断中有着广泛的应用，如功能磁共振成像、心脏磁共振成像、动态对比度增强磁共振成像等。动态磁共振成像需要快速的数据采集速度来保证数据具有合适的空间分辨率和时间分辨率。为了加速动态磁共振成像，出现了一系列加速动态磁共振成像的技术，如多线圈并行成像技术sense(pruessmannkp,weigerm,scheideggermb,boesigerp,"sense:sensitivityencodingforfastmri,"magneticresonanceinmedicine,vol.42,pp.952-962,1999)，多线圈并行成像技术利用单个磁共振信号接收线圈的空间敏感度差异来编码空间信息，降低了磁共振成像所需的梯度编码步数以加快扫描速度，对于不同的线圈而言由于其与组织的距离不同，采集到的信号强度也就不一致。靠近线圈的组织信号强，远离线圈的组织信号弱，视野以外的组织将卷折到图像对侧。并行成像将各个线圈采集的数据与敏感度信息相结合，去除单个线圈的卷折伪影，重建完整的图像。然而并行成像的加速技术受限于数据的信噪比和采样线圈的设计，在很高的加速时会导致数据难以重建。黄金角放射线径向采样技术(linw,huangf,liy,reykowskia,"grappaoperatorforwiderradialbands(growl)withoptimallyregularizedself-calibration,"magneticresonanceinmedicine,vol.64,pp.757-766,2010)对k空间中心过采样，其它地方的采集很少的数据来加速动态磁共振成像，并且黄金角放射线径向采样的数据采集方式与其它采样方式相比能更好的抵抗运动伪影。基于动态磁共振三维图像在笛卡尔空间和傅里叶域内的稀疏特性的k-t加速方法，如radialk-tspirit(santellic,schaefftert,kozerkes."radialk-tspirit:autocalibratedparallelimagingforgeneralizedphase-contrastmri,"magneticresonanceinmedicine,vol.72,pp.1233-1245,2014)。此外，还有利用人工稀疏的方法，如：arts-growl(chenz,xial,liuf,etal,"animprovednon-cartesianpartiallyparallelimagingbyexploitingartificialsparsity,"magneticresonanceinmedicine,vol.78,pp.271-279,2017)，k-tarts-growl(z.chenetal,"technicalnote:sequentialcombinationofparallelimaginganddynamicartificialsparsityframeworkforrapidfree-breathinggolden-angleradialdynamicmri:k-tarts-growl,"medicalphysics,vol.45,no.1,pp.202-213,2018)。

近年来结合并行成像、非笛卡尔黄金角放射线采样方式和压缩感知的加速动态磁共振成像的技术得到发展。在得到非笛卡尔黄金角放射线采样的动态磁共振数据后，将一定数量(通常选斐波那契数列中的数)的放射线排成一个时间帧，然后联合多通道的稀疏性结合压缩感知和并行成像可以重建出高空间时间分辨率的动态磁共振三维图像。如igrasp(fengl,grimmr,blockkt,etal,"golden-angleradialsparseparallelmri:combinationofcompressedsensing,parallelimaging,andgoldenangleradialsamplingforfastandflexibledynamicvolumetricmri,"magneticresonanceinmedicine,vol.72,no.3,pp.707-717,2014)，但是此方法只利用到了动态磁共振三维图像时间维的冗余信息，而图像空间维的冗余信息未被利用。

非均匀采样和压缩感知满足快速扫描的迫切需求，只需要采样少量数据点，通过约束图像的稀疏性，利用非线性优化算法恢复丢失数据实现磁共振图像的重建。设计基于压缩感知的磁共振图像重建方法关键在于磁共振图像的稀疏表示，目前磁共振图像的稀疏表示主要分为两种类型，正交稀疏表示(j.m.p.m.lustig,d.donoho,"sparsemri:theapplicationofcompressedsensingforrapidmrimaging,"magneticresonanceinmedicine,vol.58,no.6,pp.1182-1195,2007)和冗余稀疏表示,qu等(x.qu,d.guo,b.ning,y.hou,y.lin,s.cai,andz.chen,"undersampledmrireconstructionwithpatch-baseddirectionalwavelets,"magneticresonanceimaging,vol.30,no.7,pp.964-977,2012)、(x.qu,y.hou,f.lam,d.guo,j.zhong,andz.chen,"magneticresonanceimagereconstructionfromundersampledmeasurementsusingapatch-basednonlocaloperator,"medicalimageanalysis,vol.18,no.6,pp.843-856,2014)、zhan等(z.zhan,j.cai,d.guo,y.liu,z.chen,andx.qu,"fastmulticlassdictionarieslearningwithgeometricaldirectionsinmrireconstruction,"ieeetransactionsonbiomedicalengineering,vol.63,no.9,pp.1850-1861,2016)。并且，y.liu等人针对单通道的笛卡尔采样磁共振图像在紧标架下的重建模型提出了一种快速算法pfista(y.liu,z.zhan,j.-f.cai,d.guo,z.chen,andx.qu,"projectediterativesoft-thresholdingalgorithmfortightframesincompressedsensingmagneticresonanceimaging,"ieeetransactionsonmedicalimaging,vol.35,no.9,pp.2130-2140,2016)。

在动态磁共振三维图像重建问题中，需要保证重建结果的高时间空间分辨率，提高重建图像的质量以方便医生诊断。对于动态对比度增强磁共振图像重建还需要保证重建结果的时间-信号强度曲线(也叫动态增强曲线，指的是感兴趣区域信号强度均值随时间的变化)符合血液动力学模型，即重建结果的动态增强曲线的上升期、信号最大增强点、下降期和延长期要能够正确反映造影剂进入人体后组织血液代谢情况。此外，并行成像和非笛卡尔采样的引入也导致数值算法求解时间长。

技术实现要素：

本发明的目的在于提供一种同时利用三维动态磁共振图像时间维和空间维冗余信息并且可区分时间维稀疏和空间维稀疏各自重要性的三维稀疏约束重建方法以实现三维动态磁共振图像的快速重建。

本发明包括以下步骤：

1)提出一个时间空间联合加权稀疏约束的重建模型；

2)通过迭代运算重建三维动态磁共振图像；

3)利用重建得到的三维动态磁共振图像对结果进行修正，加快三维动态磁共振图像重建的收敛速度。

在步骤1)中，所述提出一个时间空间联合加权稀疏约束的重建模型的具体方法可为：

将测量得到的非笛卡尔欠采样动态多通道多帧的磁共振成像数据(也称为k空间)重排为列向量m，可以表示为m＝ed，其中，d表示待重建的三维动态磁共振图像经重排后的列向量，e表示对待重建图像d进行非均匀傅里叶变换和多通道灵敏度编码，可以表示为e＝fs其中f表示非均匀傅里叶变换，s表示灵敏度编码矩阵；用表示对重排后的三维动态磁共振图像的冗余稀疏变换，其中rt表示对重排后的三维动态磁共振图像的时间维施加的冗余稀疏变换，rs表示对重排后的三维动态磁共振图像的空间维施加的冗余稀疏变换；表示变换a的复共轭，其中，表示变换rt的复共轭，表示变换rs的复共轭；表示变换a的反变换；此外引入一个权值矩阵w＝diag(wt,…,wt,ws,…,ws)，其中diag(a1,a2,...al)表示一个除主对角线之外的元素全为0的矩阵，a1,a2,...al为主对角线上的元素，该权值矩阵为一个对角矩阵，对角线上的元素为权值，可以表示为wi，其中，wi＝wt表示时间维稀疏变换域系数的权重，wi＝ws表示空间维稀疏变换域系数的权重，通过调节wt，ws的取值来调节三维动态磁共振图像经时间维稀疏变换和空间维稀疏变换后的稀疏系数各自的重要性；提出的时间空间联合加权稀疏约束重建模型为：

其中，表示向量的二范数的平方，||·||1表示向量的一范数，正则化参数λ用于权衡||wad||1和两项的重要性。

在步骤2)中，所述通过迭代运算重建三维动态磁共振图像的具体方法可为：

将待重建三维动态磁共振图像d用稀疏反变换φ与稀疏变换域系数α的乘积来表示，即d＝φα；三维动态磁共振图像d可以按照下式的迭代运算来获得：

其中，dn和dn+1是在第n和第n+1次迭代得到的三维动态磁共振图像，e^*表示利用帧内采样轨迹和线圈灵敏度分布将非笛卡尔采样的k空间数据转换为笛卡尔的三维动态磁共振图像的变换，可以表示为e^*＝s^*f^*，其中f^*表示非均匀傅里叶变换f的复共轭，s^*表示s(灵敏度编码矩阵)的复共轭，γ是步长，其取值范围为0＜γ≤1；步长γ与正则化参数λ、权值wi的乘积为γλwi；表示对给定向量的所有元素进行软阈值操作，软阈值操作的定义为：给定向量z，若向量z的下标为i的元素zi的绝对值|zi|≤γλwi，则zi＝0，若|zi|＞γλwi，则zi＝sgn(zi)(|zi|-γλwi)，其中sgn(zi)是符号函数。

在步骤3)中，所述利用重建得到的三维动态磁共振图像对结果进行修正，加快三维动态磁共振图像重建的收敛速度的具体方法可为：

引入一组因子{tn}(其中n是大于等于0的整数，初始化的因子t0＝1)；已知这组因子中的第n项因子tn，通过表达式可以计算第n+1项因子tn+1；在第n+1次迭代重建后，利用第n和n+1次重建得到的三维动态磁共振图像dn和dn+1对结果进行修正，作为第n+1次重建的结果，从而加快三维动态磁共振图像重建的收敛速度；

修正的方法是：

上式中的符号表示将上式中右边的结果赋值给上式左边的dn+1；

迭代过程中，当相邻两次的重建三维动态磁共振图像dn和dn+1的l2范数误差||dn+1-dn||2/||dn||2缩小到指定的范围(小于10^-4)时，停止迭代运算；将最后一步迭代的重建结果dn+1输出为最终的重建结果。

本发明涉及高空间时间分辨率三维动态磁共振图像重建方法，尤其是涉及同时利用三维动态磁共振图像时间维和空间维冗余信息并且可区分时间维稀疏和空间维稀疏各自重要性的三维稀疏约束重建来加快高空间时间分辨率动态磁共振成像时的一种三维图像重建方法。本发明实现了动态磁共振图像的快速重建，并且保证重建结果具有较高的空间分辨率和时间分辨率。

附图说明

图1是本发明实施例连续的欠采样示意图。

图2是本发明实施例的重建流程图。

图3是本发明实施例迭代过程中目标函数的变化曲线。

图4是使用本发明方法重建出的28帧动态对比度增强三维磁共振图像。

图5是使用非均匀傅里叶变换重建出的动态对比度增强三维磁共振图像即初始输入，在打药前、动脉期、静脉期和延长期的4帧图像。

图6是使用本发明方法重建出的动态对比度增强三维磁共振图像在打药前、动脉期、静脉期和延长期的4帧图像。

具体实施方式

以下实施例将结合附图对本发明作进一步的说明。

本发明实施例使用磁场强度为3特斯拉的磁共振成像仪对志愿者的腹部进行成像。本实施例使用的采样序列为可自由呼吸的三维放射线快速低角度(flash)脉冲序列，在扫描志愿者时为志愿者静脉注射30ml生理盐水(20ml)和造影剂(10ml)的混合物。采样得到的数据的参数如下：序列的重复时间tr为3.83ms，回波时间te为1.71ms，选层厚度δkz为3mm，角度δθ为111.25度，视野fov大小为370mm×370mm，层数：40层，每层数据放射线的条数：600，每条放射线上的采样点数：768，采样时间：90s。

具体实施过程如下：

第一步：提出时间空间联合加权稀疏约束重建模型

对于使用连续的三维欠采样(如图1所示)，采样得到的实施例中的数据使用所提方法重建了其代表层的数据，代表层数据的维度为768×600×12，其中，768为每条放射线上的采样点数，600为代表层数据包含的放射线条数，12为通道数，根据代表层的数据结合采样轨迹利用非均匀傅里叶变换可估计出重建所需要的线圈灵敏度分布。将代表层的数据以每21条放射线为一个时间帧划分为4维的多通道多帧的非笛卡尔的k空间数据(768×21×12×28)并将其重排为一个列向量m，再用e^*(利用帧内采样轨迹和线圈灵敏度分布将非笛卡尔的k空间数据转换为笛卡尔的动态对比度增强三维磁共振图像)乘以上述非笛卡尔的k空间数据m后，就可以得到空间时间待重建的动态对比度增强三维磁共振图像d，即d＝e^*m。上述待重建的动态对比度增强三维磁共振图像也将作为算法的初始输入d0，然后对待重建的动态对比度增强三维磁共振图像d的时间维和空间维同时施加冗余的稀疏变换，其中时间维采用循环移位小波变换rt，采用多贝西小波基，分解级数为5，空间维为sidwt(移不变小波变换)rs，采用多贝西小波基，分解级数为4。时间维的稀疏变换rt和空间维的稀疏变换rs共同组成了对动态对比度增强三维磁共振图像的冗余稀疏变换此外，引入一个权值矩阵w＝diag(wt,...,wt,ws,…,ws)，其中，diag(a1,a2,...al)表示一个除主对角线之外的元素全为0的矩阵，a1,a2,...al为主对角线上的元素，该权值矩阵为一个对角矩阵，对角线上的元素为权值，可以表示为wi，其中，wi＝wt表示时间维稀疏变换域系数的权重，wi＝ws表示空间维稀疏变换域系数的权重，通过调节wt，ws的取值来调节动态对比度增强三维磁共振图像经时间维稀疏变换和空间维稀疏变换后的稀疏系数各自的重要性。在本实施例中设置wt＝1，ws＝0.09。提出的时间空间联合加权稀疏约束重建模型为：

其中，表示向量的二范数的平方，||·||1表示向量的一范数，正则化参数λ用于权衡||wad||1和两项的重要性，e表示对待重建的动态对比度增强三维磁共振图像d进行非均匀傅里叶变换与多通道灵敏度编码。实施例动态对比度增强三维磁共振图像的重建中，设定正则化参数λ＝0.06a0(a0为初始输入d0所有元素中模值最大的元素的模值，等于2e-4)，本发明实施例的重建流程图如图2所示。

第二步：重建动态对比度增强三维磁共振图像

待重建的动态对比度增强三维磁共振图像d通过下式的迭代运算来获得

其中，dn和dn+1是在第n和n+1次迭代时得到的动态对比度增强三维磁共振图像，表示变换a的反变换，表示变换a的复共轭，其中，表示变换rt的复共轭，表示变换rs的复共轭，e^*表示反多通道非均匀傅里叶变换(将非笛卡尔的k空间数据转换为笛卡尔的动态对比度增强三维磁共振图像)，γ是步长，实施例中设置迭代步长γ＝1，设定正则化参数λ＝0.06a0(a0为初始输入d0所有元素中模值最大的那个元素的模值)，设定时间维稀疏系数的权重wt＝1，空间维稀疏系数的权重ws＝0.09。步长γ与正则化参数λ、权值wi的乘积为γλwi。表示对给定向量的所有元素进行软阈值操作，软阈值操作的定义为：给定向量z，若向量z的下标为i的元素zi的绝对值|zi|≤γλwi，则zi＝0，若|zi|＞γλwi，则zi＝sgn(zi)(|zi|-γλwi)，其中sgn(zi)是符号函数。

第三步：加快动态对比度增强三维磁共振图像重建的速度

引入一组因子{tn}(其中n是大于等于0的整数，初始化的因子t0＝1)。已知这组因子中的第n项因子tn，通过表达式可以计算第n+1项因子tn+1。在第n+1次迭代重建动态对比度增强三维磁共振图像后，利用第n和n+1次迭代时得到的动态对比度增强三维磁共振图像dn和dn+1对迭代的动态对比度增强三维磁共振图像进行修正，从而加快动态对比度增强三维磁共振图像重建的收敛速度；

修正的方法是：

上式中的符号表示将上式中右边的结果赋值给上式左边的dn+1中。

迭代过程中，当第n和n+1次的重建动态对比度增强三维磁共振图像dn和dn+1的l2范数误差||dn+1-dn||2/||dn||2小于10^-4时重建图像稳定，停止迭代运算；将最后一步迭代的重建结果dn+1输出为最终的重建结果，本实施例的目标函数收敛曲线如图3所示，从图3可以看出，723s后目标函数值不再下降。

本实施例重建出的动态对比度增强三维磁共振图像如图4所示，本发明实施例重建结果的时间分辨率为：3.21秒/帧，空间分辨率为：1.4×1.4×3mm³，本发明实施例的初始输入在打药前、动脉期、静脉期和延长期的4帧图像如图5所示，本发明实施例的重建结果在打药前、动脉期、静脉期和延长期的4帧图像如图6所示。