CN101975935A - 一种基于部分回波压缩感知的快速磁共振成像方法 - Google Patents

Abstract

本发明公开了一种基于部分回波压缩感知的快速磁共振成像方法。现有的成像方法速度比较慢，硬件成本高。本发明方法首先采集随机变密度部分回波数据；具体方法为在k-空间中心区域高密度采集数据，四周随机稀疏采集，生成二维随机掩模；将频率编码轴上需采集的数据点处都添加二维随机掩模，形成三维随机掩模；按照生成的三维随机掩模采集k-空间数据；然后基于小波域软域值去噪的凸集投影方法重建；最后基于有限差分变换的最小化L1范数的非线性重建，具体方法对图像空间信号x进行稀疏变换，确立优化目标，求解优化目标。本发明方法将部分回波和压缩感知技术结合起来应用在MRI的数据采集中，在减小回波时间的同时，缩减了数据采集时间。

Description

一种基于部分回波压缩感知的快速磁共振成像方法

技术领域

本发明属于磁共振的图像成像领域，涉及一种基于部分回波压缩感知的快速磁共振成像方法。

背景技术

磁共振成像MRI(Magnetic Resonance Imaging)是一种利用体外测量的核磁共振信号产生体内器官物理和化学特性的断层图像成像技术。由于MRI对机体没有不良影响，具有较高的软组织分辨能力等优点，目前已在临床疾病检测中得到广泛应用。但是，MRI在应用中还经常受到数据采集时间过长的限制。

MRI的数据采集是在k-空间(频率空间)进行的。k空间编码由kz方向的层编码、ky方向的相位编码和kx方向的频率编码组成。在相位和层编码平面(ky_kz)上减少数据采集点数，可以有效地缩短数据采集时间。在磁共振血管成像MRA(Magnetic Resonance Angiography)中常用部分回波(Partial Echo)采集来缩短回波时间，减小高阶流动伪影。部分回波采集，就是只采集频率编码kx方向的一部分数据。

近年来，压缩感知(Compressed Sensing)技术突破了香浓(Shannon)采样定理关于采样速率必须高于2倍信号带宽的极限，从而在信号处理等多个领域成为新的研究热点(Candès et al.，Robust Uncertainty Principles：ExactSignal Reconstruction from Highly Incomplete Frequency Information，2004)。压缩感知技术利用信号自身或其在变换域中的稀疏性，只需随机采集少量的数据点，就可通过非线性重建等算法恢复出原始的信号来。然而，MRI受硬件扫描系统有规律往返轨迹的限制，不能做到完全的随机采集。

2007年，M.Lustig等人提出了适合MRI硬件扫描系统的伪随机数据采集方式，首次将压缩感知技术成功应用于MRI领域(M.Lustig，et.al.，SparseMRI：The Application of Compressed Sensing for Rapid MR Imaging，MagneticResonance in Medicine，58：1182-1195，2007)。他们利用变密度随机采集的方式采集数据，采用小波变换等方式作为图像的稀疏变换，用基于L1范数最小化的非线性优化方法重建图像。此后，有人将压缩感知技术应用到MRI的多个应用领域。

目前已申请的有关压缩感知技术方面的发明专利有：宽带移动通信中利用压缩感知减少导频数的信道估计方法(200910079441.4)，提出一种利用压缩感知技术的原理，降低系统估计信道时所需的导频符号个数和保证信道估计性能的方法。一种稀疏信号的重构方法(200910023785.3)提出一种稀疏信号的重构方法，主要解决从观测向量重构原始稀疏信号速率低的问题。基于压缩感知技术的分布式信源编码的方法(200910242622.4)，提出利用压缩感知技术的优点和视频图像的稀疏特性，形成一种新的分布式信源编码方法。基于正交匹配追踪的具有压缩感知过程的信号检测方法(201010032485.4)，提出一种基于正交匹配追踪的具有压缩感知过程的信号检测方法，解决了利用现有的压缩感知实现信号检测时需要重构信号所带来的资源浪费的问题。传感器网络物理信号协同压缩感知系统及感知方法(200910198056.1)，提出一种传感器网络物理信号协同压缩感知系统及感知方法。以上专利主要应用于宽带移动、无线电等信号检测和处理领域。

以上申请的发明专利表明，还没有一种专门针对MRI领域的压缩感知技术方面的专利。

发明内容

部分回波技术是有规律的采集k-空间，而压缩感知技术则是随机的采集，本发明针对现有技术的不足，将部分回波和压缩感知技术结合起来应用在MRA的数据采集中，在缩短回波时间，抑制高阶流动伪影的同时，缩减数据采集的时间，提供一种适用于硬件的随机变密度部分回波数据采集方法，并提供一种基于小波域软域值去噪的凸集投影方法POCS(Projection onto Convex Sets)重建和基于有限差分变换的最小化L1范数的非线性重建相结合的方法来重建图像。

本发明主要包括三大步骤：随机变密度部分回波数据采集、基于小波软域值去噪的POCS重建、基于有限差分变换的最小化L1范数的非线性重建。

随机变密度部分回波数据采集包括产生三维随机掩模、根据三维随机掩模采集数据两个步骤：

1-1、产生三维随机掩模

在相位编码和层编码平面产生二维的随机掩模，具体是：在信息量集中的k-空间中心区域高密度采集数据，在中心区域四周进行随机稀疏采集，生成二维随机掩模；

频率编码轴上需采集的数据点处都添加二维随机掩模，由此形成基于多层二维随机掩模的三维随机掩模。需采集的数据点在三维随机掩模中以“1”表示，不采集的数据点以“0”表示。

1-2、根据三维随机掩模采集数据

按照步骤1-1生成的三维随机掩模采集k-空间数据。

基于小波域软域值去噪的POCS重建包括初始化、利用低频信息产生相位约束、傅立叶反变换至图像空间、离散小波变换、小波域的软域值法去噪、离散小波反变换、图像空间数据合成、经过傅立叶变换转换到k-空间和k-空间数据合成：

2-1、初始化

初始化序号i和k-空间数据y_i，使得i＝0，y_i＝y；其中y代表采集到的欠采样数据；

2-2、利用低频信息产生相位约束

取步骤1-1的k-空间数据y_i的中心低频部分，做傅立叶反变换得到图像空间信号，取出该图像空间信号中的相位信号φ作为三维POCS迭代的约束条件。

2-3、傅立叶反变换至图像空间

对步骤2-1的y_i做傅立叶反变换到图像空间，形成图像空间信号，以x_i表示。

2-4、离散小波变换

用正交小波对步骤2-3的x_i进行离散小波变换。

2-5、小波域的软域值法去噪

用软域值(Soft Thresholding)的方法对步骤2-4生成的小波系数进行去噪处理。

$\mathrm{ST}({\mathrm{\ψx}}_i,\mathrm{\λ})=\frac{{\mathrm{\ψx}}_i}{\left|{\mathrm{\ψx}}_i\right|}\left(\right|{\mathrm{\ψx}}_i|-\mathrm{\λ})\_+$ (公式1)

$\mathrm{\λ}=\mathrm{\α}\×\frac{\mathrm{median}\left(\right|{\mathrm{\ψx}}_i\left|\right)}{0.6745}\×\sqrt{\ln (\mathrm{xres}\×\mathrm{yres})}$ (公式2)

其中，ψ表示步骤2-4中的小波变换，λ为域值，median表示中值滤波，α是常数，xres×yres为图像的像素点数。

2-6、离散小波反变换

将步骤2-5中去噪后的小波系数经过离散小波反变换至图像空间。

2-7、图像空间数据合成

取出步骤2-6产生的图像空间信号x_i的幅值信号|x_i|，与步骤2-2产生的相位信号φ组合成新的图像空间数据x_i′。

x_i ′＝|x_i|e^iφ                    (公式3)

2-8、经过傅立叶变换转换到k-空间

将步骤2-7产生的x_i′数据，经傅立叶变换到k-空间形成数据w_i。

2-9、k-空间数据合成

用实际采集到的y信号替换步骤(8)产生的w_i中相应位置的信号，得到信号y_i′。

2-10、重复迭代

将步骤2-9中生成的y_i′信号作为y_i送至步骤2-3，重复步骤2-3至步骤2-10n次，其中n为自然数且n≥4。

2-11、将y_i作傅立叶反变换至图像空间形成图像空间信号x。

基于有限差分变换的最小化L1范数的非线性重建包括有限差分变换(Finite Difference)、确立最小化L1范数优化目标和优化目标的求解三个步骤。

3-1、有限差分变换

采用有限差分变换对图像空间信号x进行稀疏变换。

3-2、确立最小化L1范数优化目标

假设F_u代表欠采样数据的傅立叶变换，非线性凸优化重建需解决以下问题：

$\hat{x}=\underset{x}{\min }{\left|\right|\mathrm{\&psi;x}\left|\right|}_1,s.t.{\left|\right|{y_i}^{\&prime;}-F_{u}x\left|\right|}_2<\mathrm{\&epsiv;}$ (公式4)

其中，需重建的图像由优化变量x表示，ψ为有限差分变换，优化目标为x在有限差分变换域的L1范数最小化，域值参数ε是期望的噪声水平，用以控制重建的保真度。

3-3、优化目标的求解

采用复共轭梯度法进行优化目标的求解，得到磁共振图像。

采用本发明方法可以将部分回波和压缩感知技术结合起来应用在MRI的数据采集中，在减小回波时间的同时，缩减数据采集时间，同时本发明具有以下特点：

(1)将压缩感知技术中的随机数据采集与有规律的部分回波采集相结合，解决完全随机采样不适用于实际硬件的矛盾。通过频率编码方向的欠采样缩短回波时间，减小高阶的流动伪影，通过相位编码与层编码平面的随机欠采样，缩短数据采集时间。

(2)用已采数据和已采数据的低频相位约束，通过POCS算法在图像域和k-空间域反复迭代，可以恢复大部分未采集的k-空间数据。

(3)利用小波域软域值去噪的方法抑制随机采样带来的变换域中的噪声状伪影。

(4)针对MRA图像在有限差分域中具有较高的稀疏表现，采用基于有限差分变换的最小化L1范数的非线性重建进一步优化MRA图像重建质量。

附图说明

图1是相位编码和层编码平面的二维随机掩模示意图；

图2是三维随机掩模在相位编码和频率编码截面上的示意图；

图3是POCS算法的二维图解；

图4(a)是常规全采样参考数据经过傅立叶重建的图像；

图4(b)是欠采样MRA数据的傅立叶重建结果图；

图4(c)是欠采样MRA数据用本发明方法的重建结果图。

具体实施方式

以下结合附图对本发明作进一步说明。

一种基于部分回波压缩感知的快速磁共振成像方法包括三大步骤：随机变密度部分回波数据采集、基于小波软域值去噪的POCS重建、基于有限差分变换的最小化L1范数的非线性重建。

随机变密度部分回波数据采集包括产生三维随机掩模、根据三维随机掩模采集数据两个步骤：

1-1、产生三维随机掩模

在相位编码和层编码平面产生二维的随机掩模，具体是：在信息量集中的k-空间中心区域高密度采集数据，在中心区域四周进行随机稀疏采集，生成二维随机掩模；

频率编码轴上需采集的数据点处都添加二维随机掩模，由此形成基于多层二维随机掩模的三维随机掩模。需采集的数据点在三维随机掩模中以“1”表示，不采集的数据点以“0”表示。

1-2、根据三维随机掩模采集数据

按照步骤1-1生成的三维随机掩模采集k-空间数据。

基于小波域软域值去噪的POCS重建包括初始化、利用低频信息产生相位约束、傅立叶反变换至图像空间、离散小波变换、小波域的软域值法去噪、离散小波反变换、图像空间数据合成、经过傅立叶变换转换到k-空间和k-空间数据合成：

2-1、初始化

初始化序号i和k-空间数据y_i，使得i＝0，y_i＝y；其中y代表采集到的欠采样数据；

2-2、利用低频信息产生相位约束

取步骤1-1的k-空间数据y_i的中心低频部分，做傅立叶反变换得到图像空间信号，取出该图像空间信号中的相位信号φ作为三维POCS迭代的约束条件。

2-3、傅立叶反变换至图像空间

对步骤2-1的y_i做傅立叶反变换到图像空间，形成图像空间信号，以x_i表示。

2-4、离散小波变换

用正交小波对步骤2-3的x_i进行离散小波变换。

2-5、小波域的软域值法去噪

用软域值(Soft Thresholding)的方法对步骤2-4生成的小波系数进行去噪处理。

$\mathrm{ST}({\mathrm{\&psi;x}}_i,\mathrm{\&lambda;})=\frac{{\mathrm{\&psi;x}}_i}{\left|{\mathrm{\&psi;x}}_i\right|}\left(\right|{\mathrm{\&psi;x}}_i|-\mathrm{\&lambda;})\_+$ (公式1)

$\mathrm{\&lambda;}=\mathrm{\&alpha;}\&times;\frac{\mathrm{median}\left(\right|{\mathrm{\&psi;x}}_i\left|\right)}{0.6745}\&times;\sqrt{\ln (\mathrm{xres}\&times;\mathrm{yres})}$ (公式2)

其中，ψ表示步骤2-4中的小波变换，λ为域值，median表示中值滤波，α是常数，xres×yres为图像的像素点数。

2-6、离散小波反变换

将步骤2-5中去噪后的小波系数经过离散小波反变换至图像空间。

2-7、图像空间数据合成

取出步骤2-6产生的图像空间信号x_i的幅值信号|x_i|，与步骤2-2产生的相位信号φ组合成新的图像空间数据x_i′。

x_i′＝|x_i|e^iφ                     (公式3)

2-8、经过傅立叶变换转换到k-空间

将步骤2-7产生的x_i′数据，经傅立叶变换到k-空间形成数据w_i。

2-9、k-空间数据合成

用实际采集到的y信号替换步骤(8)产生的w_i中相应位置的信号，得到信号y_i′。

2-10、重复迭代

将步骤2-9中生成的y_i′信号作为y_i送至步骤2-3，重复步骤2-3至步骤2-10n次，其中n为自然数且n≥4。

2-11、将步骤2-10产生的y_i′(此处的y_i′＝y_i)信号作傅立叶反变换至图像空间形成图像空间信号x。

基于有限差分变换的最小化L1范数的非线性重建包括有限差分变换(Finite Difference)、确立最小化L1范数优化目标和优化目标的求解三个步骤。

3-1、有限差分变换

采用有限差分变换对图像空间信号x进行稀疏变换。

3-2、确立最小化L1范数优化目标

假设F_u代表欠采样数据的傅立叶变换，非线性凸优化重建需解决以下问题：

$\hat{x}=\underset{x}{\min }{\left|\right|\mathrm{\&psi;x}\left|\right|}_1,s.t.{\left|\right|{y_i}^{\&prime;}-F_{u}x\left|\right|}_2<\mathrm{\&epsiv;}$ (公式4)

其中，需重建的图像由优化变量x表示，ψ为有限差分变换，优化目标为x在有限差分变换域的L1范数最小化，域值参数ε是期望的噪声水平，用以控制重建的保真度。

3-3、优化目标的求解

采用复共轭梯度法进行优化目标的求解。

以下结合MRA图像的部分回波压缩感知欠采集和重建进行实例说明。假设要采集的MRA图像的三维矩阵大小为kx×ky×kz＝512×384×64。首先，产生部分回波随机变密度三维掩模，如图1所示，三维掩模在相位编码与层编码平面二维的随机变密度欠采样率为70％，即该平面只采集384×64×30％＝7372个点。如图2所示，在频率编码方向部分回波的欠采样率为30％，只采集362个点。欠采样总量高达79％。

用三维梯度回波序列按三维掩模采集MRA数据。采集的数据先进行三维POCS重建，图3为三维POCS重建的二维示意图。如图3所示：初始化时，i＝0，令y_i＝y；取y_i的低频信号，在其周围填零后，做傅立叶反变换，取出变换后的相位信号φ作为约束条件；同时对y_i做傅立叶反变换，变换到图像空间，以x_i表示；再对x_i进行基于正交小波基的小波变换；在小波域，用公式(1)所示的软域值算法对小波系数去噪；去噪后的小波系数经过小波反变换到图像空间；取出图像空间信号的幅值信号|x_i|，与相位约束信号φ组合成新的x_i′信号；x_i′信号经过傅立叶变换到k-空间，形成w_i信号；然后用实际采集到的y信号替换w_i信号中相应位置的信号，形成y_i′，将y_i′作为下一次迭代的y_i信号；重复迭代5次，对y_i′作傅立叶反变换形成图像空间信号x。

POCS重建后恢复出的k-空间数据y_i′作为参考数据，按照公式(4)所示，进行基于有限差分变换的最小化L1范数的非线性重建。需重建的图像为非线性重建中需优化的变量x，以x在有限差分域中的L1范数最小为优化目标，以x在k-空间中相对参考数据y_i′的连续性为约束条件，用复共轭梯度的方法求解得到最优的重建图像

图4(a)为常规全采样参考数据采集经过傅立叶重建的图像，图4(b)为欠采样的MRA数据经过直接傅立叶重建的图像，图4(c)为本发明所用方法的重建图像。从放大的图中可以看到，欠采样的MRA数据经过直接傅立叶重建的图像质量较差，如箭头所示处，很多中小血管信息丢失了，而它们却在图4(c)中却清晰可见。

相比直接傅立叶重建方法，用本发明方法重建的图像质量更好。在高达79％的欠采样下，相比参考图像，采用本发明方法重建图像中的小动脉可见性仅略有下降。高欠采样带来数据采集时间的缩短，因此，部分回波压缩感知方法为临床MRA数据采集缩减时间提供了可能。

Claims (1)

1.一种基于部分回波压缩感知的快速磁共振成像方法，其特征在于该方法包括如下步骤：

步骤(1)采集随机变密度部分回波数据；具体方法为：

1-1、在信息量集中的k-空间中心区域高密度采集数据，在k-空间中心区域四周随机稀疏采集，生成二维随机掩模；将频率编码轴上需采集的数据点处都添加二维随机掩模，由此形成基于多层二维随机掩模的三维随机掩模；

1-2、按照步骤1-1生成的三维随机掩模采集k-空间数据；

步骤(2)基于小波域软域值去噪的凸集投影方法重建；具体方法为：

2-1、初始化序号i和k-空间数据y_i，使得i＝0，y_i＝y；其中y为采集到的欠采样数据；

2-2、取步骤2-1的k-空间数据y_i的中心低频部分，做傅立叶反变换得到图像空间信号，取出该图像空间信号中的相位信号φ作为三维凸集投影方法迭代的约束条件；

2-3、对y_i做傅立叶反变换到图像空间，形成图像空间信号x_i；

2-4、用正交小波对x_i进行离散小波变换；

2-5、用软域值法对步骤2-4中所生成的小波系数进行去噪处理；

2-6、将去噪后的小波系数经过离散小波反变换至图像空间；

2-7、将步骤2-6产生的图像空间信号x_i的幅值信号|x_i|与相位信号φ组合成新的图像空间数据x_i′，x_i′＝|x_i|e^iφ；

2-8、将x_i′经傅立叶变换到k-空间形成数据w_i；

2-9、用采集到的y信号替换步骤2-8产生的w_i中相应位置的信号，得到信号y_i′；

2-10、令y_i＝y_i′，重复步骤2-3至步骤2-10n次，其中n为自然数且n≥4；

2-11、将y_i作傅立叶反变换至图像空间形成图像空间信号x；

步骤(3)基于有限差分变换的最小化L1范数的非线性重建；具体方法为：

3-1、采用有限差分变换对图像空间信号x进行稀疏变换；

3-2、确立最小化L1范数优化目标，优化目标为x在有限差分变换域的L1范数最小化，即

其中ψ为有限差分变换，ε为域值参数，F_u为欠采样数据的傅立叶变换，为磁共振图像；

3-3、采用复共轭梯度法求解优化目标，得到磁共振图像。

Images