CN108717171A

CN108717171A - 一种压缩感知低场磁共振成像算法

Info

Publication number: CN108717171A
Application number: CN201810511209.2A
Authority: CN
Inventors: 柴青焕; 苏冠群; 侯学文; 常晓; 聂生东
Original assignee: University of Shanghai for Science and Technology
Current assignee: University of Shanghai for Science and Technology
Priority date: 2018-05-24
Filing date: 2018-05-24
Publication date: 2018-10-30

Abstract

本发明涉及一种压缩感知低场磁共振成像算法，对采集的低场核磁共振原始数据进行欠采样处理，得到随机变密度K空间数据；欠采样的K空间数据与压缩感知理论相结合建模，建模问题成一个包含数据保真项、稀疏先验项和全变分项的线性组合最小化问题，将双树小波变换与小波树稀疏联合作为压缩感知理论中的稀疏变换；利用双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法对其进行求解，得到欠采样的变密度K空间数据重建出的低场磁共振图像。引入双树小波变换作为CS‑MRI模型中的稀疏变换，克服了传统小波变换具有的敏感性、移变性、缺乏方向性和重构图像带来混叠伪影的缺陷。计算精度高，鲁棒性好，提高成像速度的同时，提高重建图像的信噪比。

Description

一种压缩感知低场磁共振成像算法

技术领域

本发明涉及一种核磁共振成像技术，特别涉及一种双树小波变换与小波树稀疏联合的压缩感知低场磁共振成像算法。

背景技术

在世界范围内，磁共振技术发展迅猛，在很多领域都得到了广泛的发展。根据背景强度的不同，磁共振技术可分为高场、中场、低场三类，高场磁共振强度≥3T、中场磁共振强度为0.5T＜中场＜3T、低场磁共振强度≤0.5T(T，特拉斯)。高场磁共振技术的磁源多采用超导磁体产生高强磁场，采集到的信号幅值强，信噪比高。但超导磁体对设备环境要求高(如必须在液态氮冷却条件下使用)且购买、安装和维护非常昂贵。中场磁共振技术的磁源采用超导磁体、永磁体或电磁，采集到的信号信噪比较高，中场技术一般用于医用磁共振成像系统等。低场磁共振技术的磁源一般采用永磁体或电磁(室温下可使用)，采集到的信号幅值很弱、信噪比极低。但是低场磁共振设备与中高场相比具有便携、廉价、制造维护应用成本低，且测量方便、快速、无损、无副作用等优点。低场磁共振技术已被广泛应用在食品化工、在线无损检测、聚合物材料、灾害防治、生物医药、能源勘探、石油测井、岩心分析等各种领域；也逐渐成为一种必要的分析检查手段，在科研和生产中也发挥了无可替代的作用。近年来随着电磁技术、计算机技术的发展，低场磁共振成像(Magnetic Resonance Imaging，MRI)技术也迅速发展起来，高性能的低场开放型永磁MRI系统开始受到人们的青睐。

压缩感知(Compressive Sensing，CS)是一个新的信号采样和压缩理论。它突破了Nyquist采样定理，利用信号在特定域上的稀疏性及可压缩性，在远小于Nyquist采样率的条件下，从少量的信号中通过非线性重建算法恢复出原始信号。对MRI而言，大部分MR图像都具有稀疏性且其采样的原始数据K空间为抽象的频率空间，存储的是图像的傅里叶数据。传统的MRI技术需要很长的扫描时间，使重建图像的分辨率不理想并且需要花费昂贵的成本。而利用CS理论进行MRI(CS-MRI)能够极大地减少傅里叶变换域的采样数据，降低扫描时间。CS-MRI在不显著降低图像质量的情况下克服了MRI成像速度慢的缺点；验证了从少量稀疏测量数据重建出图像的可能性。低场CS-MRI重建技术是通过借鉴现已发展成熟的中高场CS-MRI技术，吸收其优点并规避其缺点来实现的。

对于CS-MRI，有两个关键点需要进一步研究。第一：稀疏变换，基于CS理论的MRI技术能从欠采样的稀疏图像数据中重建出原始图像，而MR图像本身不具有稀疏性，因此需要将图像在稀疏域进行稀疏表示。传统CS-MRI的稀疏变换基有小波变换、全变分(TotalVariation，TV)等，这些稀疏基可能不提供足够的稀疏表示。如离散小波变换(DiscreteWavelet Transform，DWT)有三个缺点：具有转移敏感性、缺乏方向性、缺乏相位信息且DWT不能稀疏表示二维信号的一维奇异点，很难捕捉到规则的图像轮廓。因此，可以通过改善，减少DWT的缺点。如文献中使用双树小波变换(Dual-Tree Complex Wavelet Transform，DT-CWT)代替小波变换作为稀疏基，弥补了传统小波变换的一些缺点。第二：重建算法，由于重建过程是一个病态问题，我们应该使用正则约束项去求得病态问题的稳定解。近年来，已经提出了大量的算法用于求解CS-MRI模型。如：共轭梯度(Conjugate Gradient，CG)算法求解全变分(Total Variation，TV)和小波稀疏联合的正则化目标函数问题；快速复合分裂算法(Fast Composite Splitting Algorithm，FCSA)将正则化问题分裂为两个简单的子问题，每个子问题再通过快速收缩算法求解，使正则化问题更有效的得到解决；一种先进的小波树稀疏MRI(Wavelet Tree Sparsity MRI，WaTMRI)重建算法，它在TV和小波稀疏联合正则项基础上增加了小波树稀疏正则项；与先前的算法相比使测量的复杂性由O(K+Klogn)减少到O(K+logn)。

发明内容

本发明是针对低场磁共振技术存在的问题，提出了一种压缩感知低场磁共振成像算法，改善低场MRI的图像质量及提高低场MRI的成像速率。

本发明的技术方案为：一种压缩感知低场磁共振成像算法，具体包括如下步骤：

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

2)原始数据作为全采样K空间数据，进行欠采样处理操作，得到模拟的随机变密度K空间数据：

利用变密度随机采样矩阵与全采样K空间数据矩阵相乘得到随机变密度欠采样K空间数据，变密度随机采样矩阵为由0和1组成的大小为m行n列的笛卡尔采样矩阵，定义欠采样K空间数据的欠采样率为m/n；

3)将双树小波变换与小波树稀疏联合作为压缩感知理论中的稀疏变换；

4)将步骤2)得到的欠采样K空间数据与步骤3)压缩感知理论相结合构建模型，将磁共振图像重建问题建模成一个包含数据保真项、稀疏先验项和全变分项的线性组合最小化问题；

5)根据步骤4)建立的模型，提出一种双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法，并用算法对模型进行求解，得到欠采样的变密度K空间数据重建出的低场磁共振图像。

所述双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法具体是：在步骤4)所建CS-MRI模型中稀疏变换采用双树小波变换，树结构被建模为一个重叠的组正则项，TV作为惩罚项被用在重建公式中来抑制噪声，TV惩罚项和L1范数约束项作为树形结构模型的补充。

本发明的有益效果在于：本发明压缩感知低场磁共振成像算法，引入双树小波变换作为CS-MRI模型中的稀疏变换，克服了传统小波变换具有的敏感性、移变性、缺乏方向性的缺陷，并克服小波变换作为稀疏变换对复杂信号不适用及给重构图像带来混叠伪影的问题。本发明提出的基于双树小波变换和小波树稀疏联合的压缩感知低场磁共振图像重建模型；计算精度高，鲁棒性好，在提高成像速度的同时，提高了重建图像的信噪比。

附图说明

图1为本发明压缩感知低场磁共振成像算法流程图；

图2a为采用FCSA算法对matlab自带phantom仿真数据欠采样进行重建的结果图；

图2b为采用WaTMRI算法对matlab自带phantom仿真数据欠采样进行重建的结果图；

图2c为采用本发明DT-WaTMRI算法对matlab自带phantom仿真数据欠采样进行重建的结果图；

图3a为本发明采用FCSA算法对低场条件下采集的哈密瓜数据欠采样进行重建的结果图；

图3b为采用WaTMRI算法对低场条件下采集的哈密瓜数据欠采样进行重建的结果图；

图3c为采用本发明DT-WaTMRI算法对低场条件下采集的哈密瓜数据欠采样进行重建的结果图；

图4a为采用FCSA算法对低场条件下采集的小鼠数据欠采样进行重建的结果图；

图4b为采用WaTMRI算法对低场条件下采集的小鼠数据欠采样进行重建的结果图；

图4c为采用本发明DT-WaTMRI算法对低场条件下采集的小鼠数据欠采样进行重建的结果图；

图5a为本发明对phantom仿真数据重建结果用RMSE评价指标的客观评价图；

图5b为本发明对phantom仿真数据重建结果用PSNR评价指标的客观评价图；

图5c为本发明对phantom仿真数据重建结果用SSIM评价指标的客观评价图；

图6a为本发明对采集的低场哈密瓜数据重建结果用RMSE评价指标的客观评价图；

图6b为本发明对采集的低场哈密瓜数据重建结果用PSNR评价指标的客观评价图；

图6c为本发明对采集的低场哈密瓜数据重建结果用SSIM评价指标的客观评价图；

图7a为本发明对采集的低场小鼠数据重建结果用RMSE评价指标的客观评价图；

图7b为本发明对采集的低场小鼠数据重建结果用PSNR评价指标的客观评价图；

图7c为本发明对采集的低场小鼠数据重建结果用SSIM评价指标的客观评价图。

具体实施方式

如图1所示，一种双树小波变换与小波树稀疏联合的压缩感知(低场)磁共振成像算法，包括如下步骤：

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

2、原始数据作为全采样K空间数据，进行模拟欠采样处理操作，得到随机变密度欠采样K空间数据；

3、磁共振图像本身或其在某个变换域内越稀疏，其重建质量就越好，将双树小波变换与小波树稀疏联合作为压缩感知理论中的稀疏变换；

4、将步骤2得到的欠采样K空间数据与步骤3压缩感知理论相结合构建模型，将磁共振图像重建问题建模成一个包含数据保真项、稀疏先验项和全变分项的线性组合最小化问题。

5、根据步骤4建立的模型，提出一种双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法，并用算法对模型进行求解，得到欠采样的变密度K空间数据重建出的低场磁共振图像。

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

MR图像的数据空间是K空间，在K空间内存储在不同位置的MR信号对图像的贡献有所差异，K空间中心部分，所对应的MR信号空间频率低、幅度大，主要形成图像的对比度；K空间的外围部分，所对应的MR信号空间频率高、幅度低，主要形成图像的分辨力。由于K空间的冗余特性采用任意变密度的采样模式，都能缩短采样时间、提高成像速度、减少运动伪。但是，不同采样脉冲序列产生不同类型的K空间模型，而不同的K空间模型需选用不用的随机变密度采样模式。本发明中实验所用低场磁共振采样脉冲序列为多回波SE序列，所得到的K空间为笛卡尔单向平行轨迹，因此本发明实验中采用一维随机变密度采样方式，来进行磁共振图像重建。

随机测量矩阵能在很大概率上满足有限等距RIP性质。RIP性质是CS理论完整恢复出原始信号的保证。但是随机测量矩阵硬件实现困难，所以当前实际应用中测量矩阵一般选用伪随机确定性测量矩阵。因此，基于CS理论的MRI对K空间数据欠采样也选用变密度的随机测量矩阵。CS-MRI模型对K空间以随机变密度采样方式在K空间中心部分的采样权重多于K空间外围部分，仅需采集Nyquist采样定理要求数据量的30％甚至更少的K空间数据，来重建原始图像，即加快了采样速率又保证了成像质量，而且也降低了存储要求。

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

2.1利用变密度随机采样矩阵与全采样K空间数据矩阵相乘得到随机变密度欠采样K空间数据；

2.2变密度随机采样矩阵是由0和1组成的大小为m行n列的笛卡尔采样矩阵；

2.3我们定义欠采样K空间数据的欠采样率为m/n。

B、步骤3的原理如下：

稀疏变换在CS理论中扮演着至关重要的角色，CS-MRI重建的结果很大程度上取决于稀疏变换。一般来说，信号越稀疏重建所需的测量值越少。除此之外，当存在测量不足和/或噪声时，重构图像中伪影的消除也依赖于稀疏转换。例如，有限差分(即全变分)变换对分段常量信号重建实施稀疏非常有效。但是，当测量值不足时，有限差分会导致CS-MRI模型重建图像的块模糊和细节丢失。同样地，CS-MRI中经常使用的另一稀疏变换——DWT也有一些缺点，第一：DWT系数常常在奇异点(边缘)周围振荡使这种奇异点的提取变得复杂；从而，无法恢复奇异点周围所有的DWT系数导致CS-MRI重建图像产生边缘伪影。第二：DWT具有移变性，这意味着信号微小的偏移会严重扰乱其小波系数。第三：在DWT中非理想的低通滤波器和高通滤波器产生了相当大的混叠，当这种混叠在理想条件下逆DWT被取消时，DWT系数的任何扰动(如噪声、阈值等)会导致重构图像中的混叠伪影。第四：DWT缺乏方向性不能有效地表示曲线。而DT-CWT是对DWT的一种增强，它通过使用复值小波和尺度函数弥补了DWT的缺点。

除了小波系数的稀疏性，树结构也为压缩感知成像提供了好的技术优势。小波系数之间的统计依赖性通常是用如隐藏的马尔可夫树(Hidden Markov Tree，HMTs)这样的图形模型表示，即小波系数的树结构。根据CS理论如果一个信号是稀疏的或者能够被稀疏表示出来，那么重构该信号所需的采样率比Nyquist-Shannon抽样定理要小得多。此外，如果我们事先知道原始信号的结构，如组或图结构，则可以进一步减少采样率。当前提出的一些算法已经开始利用小波系数的树结构来改进CS重构，这些算法很好的利用了树结构的设想，使用图解模型去表示小波树结构，小波树结构的父值决定每个系数的分配，如：树结构中的父-子关系，父系数有一个大的或小的值，那么他的子系数值往往与其保持一致。以上的先验知识如果能够被完全利用，仅仅O(K+logn)的观测数据就能恢复出满足树稀疏结构的原始数据；比CS理论中标准K稀疏数据需要O(K+Klog，n)的采样数据才足够恢复出长度为n的原始数据减少了接近一倍的采样量。当n无限大时利用树结构的优势就更加显著了。目前有文献验证了小波树结构也能帮助提高磁共振成像速度。

基于以上理论，本发明将双树小波变换与小波树稀疏联合作为压缩感知理论中的稀疏变换。以便在提高CS-MRI技术的成像速度的同时提高成像质量。

C、步骤3的原理基于如下：

磁共振应用压缩感知理论重建，就是将欠采样得到的少量K空间数据，进行某种稀疏变换，然后再结合L1范数凸优化问题，最终重构出图像。CS-MRI的数学模型可表示为式(1)

其中x是待重建的磁共振图像，Ψ表示线性稀疏变换，F_S表示欠采样傅里叶变换。y是从MRI扫描仪欠采样测量得到的K空间数据。ε控制重建数据和测量数据的保真度，代表测量数据的噪声水平。利用拉格朗日乘子法将(1)式转换为无约束优化问题形式，如式(2)所示

其中λ为正则化参数，对(2)式进行求解即可重建出磁共振图像。近年来提出的各种CS-MRI模型都是基于改变(2)式的稀疏变换基或加入不同的正则项演化而来。目前，CS-MRI模型的重建多采用全变分和小波稀疏约束项线性结合的方式来重建MR图像，如式(3)所示。

其中α和β是两个正则参数，Ψ_d是任意稀疏变换，和分别为图像水平和垂直方向的梯度算子。

D、提出一种双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法，并对其进行求解得到欠采样的变密度K空间数据重建出的低场磁共振图像。

DT--WaTMRI与已有CS-MRI重建算法的不同主要是将稀疏变换基由小波变换改变为双树小波变换。算法的求解步骤阐述如下：

提出一种小波树结构稀疏、双树复小波稀疏变换和梯度(TV)稀疏联合的新的CS-MRI模型来进行低场磁共振的成像以便在提高其成像速度的同时提高其成像质量。该CS-MRI模型中稀疏变换采用双树小波变换，树结构被建模为一个重叠的组正则项，TV作为惩罚项也被用在重建公式中来抑制噪声。TV惩罚项和L1范数约束项实际上是作为树形结构模型的补充，使我们新提出CS-MRI模型对实际低场MR图像更具鲁棒性。实现以上条件的目标函数如下(4)式所示。

式中Φ表示双树复小波变换系数，G表示所有父-子系数组的集合，x_g是属于母-子系数组的所有数据。

(4)式所表示的树结构模型算法引入了TV项和L1项与之前的树结构模型算法大有不同，但是不能有效的求解，因此引入变量z约束x重叠结构，使(4)式变为非重叠凸优化问题，令GΦx＝z，则(4)变为下式：

公式(5)最后一项是正则化项也即罚范数二范数，正则化项本质上是一种先验信息，该项对模型中向量进行’惩罚’，从而避免最小二乘问题的过拟合问题。

对x子问题我们可以联合(5)式中第一个和最后一个二次项惩罚项来求解，而z子问题通过分组软阈值来解决，剩下部分与FCSA算法类似，可以通过迭代方案有效地解决。z子问题能被写作如下(6)式所示：

其中gi表示第i组，s表示全部组数。z_gi与(GΦx)_gi都表示的是对待重建的磁共振图像x进行双小波稀疏变换后的系数，这些系数之间的统计依赖关系可以用马尔可夫树这样的图解模型表示，即使用图解模型去表示小波树结构，而小波树结构的父值决定每个子系数的分配。这是增加到目标函数中的一个二次惩罚函数，控制解的范数或模，λ＞0权衡测量保真度，λ越大，越小，可能为零。

Z子问题有一个封闭式的软阈值解决方案如下(7)式所示：

式中ri＝(GΦx)_gi。令这是一个有着莱布尼兹常数(Lipschitz，L_f)的凸平滑函数，再令g₁(x)＝α||x||_TV，g₂(x)＝β||Φx||₁，g₁(x)和g₂(x)是凸的非平滑函数。x子问题能通过FCSA算法有效的解决。为方便起见，我们将(7)式表示为：

我们提出的算法总结如下表1所示DT-WaTMRI算法主要步骤：

表1

算法中L_f表示为Lipschitz常数(L_f＞0)、N表示最大迭代次数、x表示待重建的图像、r迭代时存放相应值的中间变量、t表示适当的迭代步长、表示函数f在点rk处的梯度，函数其中F_S ^T表示部分逆傅里叶转换。

任意标量ρ＞0时，近端映射(Proximal map，prox)归一化连续凸函数g(x)的定义如下(8)式所示：

公式(9)是对x的正交投影操作，将图像的像素值规范化到[l，u]范围内，且u＞l之0，若不执行该操作，则重建图像由于存在负数像素数值而出现伪影。

本发明一种双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法的效果能够通过以下实验进一步说明。

1、仿真实验：

实验测试数据分为三组，分别为MATLAB自带phantom仿真数据和使用苏州纽迈分析仪器有限公司的低场磁共振成像实验仪，对实验样本哈密瓜及小鼠进行实验采集所得的fid格式的大小为256*128的全K空间数据。对以上三种样本的全K空间数据欠采样得到变密度K空间数据，结合CS-MRI模型，利用双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法求解凸优化问题，重建出欠采样数据下的MR图像。

2、仿真实验结果

图2a-2c、图3a-3c、图4a-4c分别是对MATLAB自带phantom仿真数据、低场磁共振采集的哈密瓜及小鼠全K空间数据，采用三种不同的方法进行重建得到实验结果结果，这三幅图中2a、3a、4a图显示的采用FCSA算法重建结果；2b、3b、4b图是采用WaTMRI算法的重建结果；2c、3c、4c图是采用本发明算法DT-WaTMRI方法的重建。

图5a-5c是对phantom仿真数据重建结果用RMSE、PSNR、SSIM评价指标的客观评价。图6a-6c是对采集的低场哈密瓜数据重建结果用RMSE、PSNR、SSIM评价指标的客观评价。图7a-7c是对采集的低场小鼠数据重建结果用RMSE、PSNR、SSIM评价指标的客观评价。

为了进一步验证本发明算法提高了重建效率，表2列举了不同算法对不同样本的重建时间结果，重建时间单位为秒/s。

表2

3、仿真实验分析

通过以上图形及表2的结果，表2中对所有实验样本来说，稀疏变换采用传统小波变换的FCSA算法重建出MR图像所用时间比稀疏变换采用传统小波变换与结构稀疏联合的WaTMRI算法少，但是WaTMRI算法与FCSA算法相比牺牲了重建时间来换取重建质量。而稀疏变换采用双树小波变换与结构稀疏联合的DT-WaTMRI算法重建出MR图像所用时间比WaTMRI算法少，在保证重建质量的同时提高了成像效率。图3是不同全采样数据经一维变密度随机欠采样后简单填零的重建结果，由结果可知对K空间数据欠采样后直接傅里叶变换进行简单填零重建所得到的重建图像伪影比较严重。图4是利用不同的CS-MRI重建算法对欠采样数据进行重建的结果，由重建结果可知本文提出的DT-WATMRI算法与FCSA算法和WaTMRI算法相比能显著改善图像的重建质量。图5是用评价图像质量的客观指标RMSE、PSNR、SSIM对不同算法的重建结果进行评价的结果。由评价结果可知对任何实验样本来说，随着迭代次数的增加每种算法都趋于收敛。而本文所提出的DT-WATMRI算法的重建效果与FCSA算法和WaTMRI算法相比，其RMSE值远小于、PSNR值、SSIM值远大于FCSA和WaTMRI两种算法。因此结合理论和实践我们可以得到以下结论：本文提出的DT-WATMRI算法改进了现有的重建算法，且能在降低采样时间的情况下，提高磁共振图像的重建质量。

Claims

1.一种压缩感知低场磁共振成像算法，其特征在于，具体包括如下步骤：

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

2）原始数据作为全采样K空间数据，进行欠采样处理操作，得到模拟的随机变密度K空间数据：

3）将双树小波变换与小波树稀疏联合作为压缩感知理论中的稀疏变换；

4）将步骤2)得到的欠采样K空间数据与步骤3)压缩感知理论相结合构建模型，将磁共振图像重建问题建模成一个包含数据保真项、稀疏先验项和全变分项的线性组合最小化问题；

5）根据步骤4）建立的模型，提出一种双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法，并用算法对模型进行求解，得到欠采样的变密度K空间数据重建出的低场磁共振图像。

2.根据权利要求1所述压缩感知低场磁共振成像算法，其特征在于，所述双树小波变换与小波树稀疏联合的压缩感知磁共振成像算法具体是：在步骤4）所建CS-MRI模型中稀疏变换采用双树小波变换，树结构被建模为一个重叠的组正则项，TV作为惩罚项被用在重建公式中来抑制噪声，TV惩罚项和L1范数约束项作为树形结构模型的补充。