CN103323805B - 基于小波域稀疏表示的speed快速磁共振成像方法 - Google Patents

Abstract

本发明公开了一种基于小波域稀疏表示的SPEED快速磁共振成像方法。本发明主要包括两大步骤：有规律的欠采样数据采集、基于小波域数据稀疏表示的SPEED图像重建。有规律的欠采样数据采集包括k-空间有规律的欠采样和k-空间中心区域全采样两个步骤。基于小波域数据稀疏表示的SPEED图像重建包括填零傅立叶重建、离散小波变换、确立小波域双层稀疏鬼影模型、基于最小平方误差法的重叠鬼影分离、离散小波逆变换、多个鬼影子图的配准求和、图像重建七个步骤。采用本发明方法可以将离散小波变换用于SPEED的数据稀疏表示，通过离散小波变换来提高图像的稀疏性，进而提高SPEED重建图像的质量。

Description

基于小波域稀疏表示的SPEED快速磁共振成像方法

技术领域

本发明属于磁共振的图像成像领域，涉及一种基于小波域数据稀疏表示的SPEED快速磁共振成像方法。

背景技术

磁共振成像(Magnetic Resonance Imaging, MRI) 是一种利用体外测量的核磁共振信号产生体内器官物理和化学特性的断层图像成像技术。由于MRI对机体没有不良影响，具有较高的软组织分辨能力等优点，目前已在临床疾病检测中得到广泛应用。但是，MRI在临床应用中还常受到数据采集时间过长的限制。长期以来，研究人员通过提高MRI梯度性能、研制快速序列和高效的k-空间（频率空间）数据采集轨迹、重建算法和并行数据采集等方式来提高MRI数据采集速度，取得了很多研究成果，但是还存在很多问题。例如，过高的梯度值和快速的梯度切换会产生病人的周围神经刺激，病人的生理限制了梯度系统性能的发挥。

MRI的数据采集是在k-空间进行的。k-空间编码由相位编码（Phase Encoding, PE）和频率编码（Frequency Encoding, FE）组成。在相位编码方向上减少数据采集点数，可以有效地缩短数据采集时间。因此，如何在不降低图像质量条件下减少数据采集总量成为进一步快速磁共振成像的关键。

近年来，压缩感知（Compressed Sensing, CS）技术突破了香浓-奈奎斯特（Shannon-Nyquist）采样定理关于采样速率必须高于2倍信号带宽的极限，从而在MRI等多个领域成为新的研究热点（M Lustig, et. al., Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging, Magnetic Resonance in Medicine, 58:1182-1195, 2007）。CS技术利用信号自身或其在变换域中的稀疏性，只需随机采集少量的数据点，就可通过非线性重建等算法恢复出原始的信号来。然而，MRI受硬件扫描系统有规律往返轨迹的限制，不能做到完全的随机采集，难以确保压缩感知重建时需要的不相关性（Incoherent），从而影响了图像的质量；而且随机的采集方式也不易于临床实现；此外，基于非线性迭代的重建方式使得CS技术的图像重建时间很长，尤其在MRI三维或动态数据采集时，更难以实现实时重建。

SPEED（Skipped Phase Encoding and Edge Deghosting）也是一种通过减少数据采集总量来缩短数据采集时间的MRI快速成像方法（QS Xiang, Accelerating MRI by skipped phase encoding and edge deghosting (SPEED)，Magnetic Resonance in Medicine, 53:1112-1117, 2005）。和CS类似，SPEED也利用了信号在变换域的稀疏特性，但是，与CS技术不同，SPEED通过简单有规则的欠采样来采集数据，然后基于解析法来重建图像，相比CS技术中的非线性迭代重建，SPEED解析求解过程非常快速。由于SPEED采集数据时无需CS那种复杂的随机采样，只需简单有规律的欠采样，不但易于实现，也易于和现有的采集方式结合，是一种很有应用潜力的成像方式。

目前已申请的关于快速成像方面的MRI专利有：一种基于部分回波压缩感知的快速磁共振血管成像方法（申请号：201010272089.8）；一种基于CS压缩感知技术的快速磁共振成像方法（申请号：201010272095.4）；快速磁共振成像方法及系统（申请号：201210013858.2）。目前还未能查询到任何关于SPEED快速成像方面的授权发明专利或申请。国内外已发表的SPEED快速成像方面的文章有：

2006年，常征和向清三将SPEED算法与并行成像技术进一步结合，提出了SPEED-ACE成像法（Chang Z and Xiang QS. Highly accelerated MRI by skipped phase encoding and edge deghosting with array coil enhancement (SPEED-ACE). Med Phys. 33:3758-3766, 2006），利用多个采集线圈共同采集k-空间欠采样数据，进一步提高成像速度。2007，基于MRA数据本身就非常稀疏的特性，常征和向清三将SPEED的双层模型简化到单层模型（Chang Z and Xiang QS. Simplified skipped phase encoding and edge deghosting (SPEED) for imaging sparse objects with applications to MRA. Med Phys. 34:3173-3182, 2007），提出了S-SPEED（Simplified-SPEED）算法，在保持图像质量不变的情况下进一步减少了数据采集的总量，缩短了数据采集时间。S-SPEED适用于数据本身就非常稀疏的场合，例如暗背景亮信号的MRA应用（Chang Z, Xiang QS, Shen H and Yin FF. Accelerating non-contrast-enhanced MR angiography with inflow inversion recovery imaging by skipped phase encoding and edge deghosting (SPEED). Journal of Magnetic Resonance Imaging. 31:757-765, 2010）。2009年，常征等人提出EMA-SPEED（Efficient Multiple Acquisition by SPEED）算法（Chang Z, Xiang QS, Ji J, and Yin FF. Efficient multiple acquisitions by skipped phase encoding and edge deghosting (SPEED) using shared spatial information. Magn Reson Med. 61:229-233, 2009.），通过共享多个采集间的相似空间信息进一步缩短了SPEED的数据采集时间，从而可获得比单次采集更高的加速比。

以上发表的关于SPEED快速磁共振成像的文章，重建时都是基于差分的方式进行数据的稀疏表示的，还未公开过基于离散小波变换的数据稀疏表示方法进行SPEED快速成像。

发明内容

本发明针对现有SPEED技术的不足，将基于小波域数据稀疏表示用于SPEED快速磁共振成像中，提供了一种新的SPEED数据稀疏表示方法，提高了SPEED快速成像技术的图像质量。

本发明主要包括两大步骤：有规律的欠采样数据采集、基于小波域数据稀疏表示的SPEED图像重建。

有规律的欠采样数据采集包括k-空间有规律的欠采样和k-空间中心区域全采样两个步骤。

1-1、k-空间有规则的欠采样

在k-空间的相位编码PE方向每隔N行采集一行k-空间数据，共采集三组，得到三组欠采样的k空间数据S ₁、S ₂ 和S ₃。用 d ₁, d ₂, d ₃表示每组欠采样数据在PE方向上的偏移量，采样方式用N (d ₁, d ₂, d ₃)表示。

1-2、k-空间中心区域全采样

在信息量较集中的k-空间中心区域采集全部数据，跟据需采集数据矩阵的大小，相应地采集8至32行k-空间中心区域的相位编码数据。

基于小波域数据稀疏表示的SPEED图像重建包括填零傅立叶重建、离散小波变换、确立小波域双层稀疏鬼影模型、基于最小平方误差法（Least Square Error，LSE）的重叠鬼影分离、离散小波逆变换、多个鬼影子图的配准求和、图像重建七个步骤。

2-1、填零傅立叶重建

对于三组欠采样的数据，k-空间中没有采集的数据用0表示，进行常规的填零傅立叶重建，重建后图像分别为I ₁、I ₂和I ₃。k-空间中每隔N行采集一行数据使得对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多可能有N个重叠的鬼影。例如，当N = 6时， I ₁、I ₂和I ₃上分别有6层重叠的鬼影。

2-2、离散小波变换

对步骤2-1得到的图像I ₁、I ₂和I ₃分别施加离散小波变换Ψ，得到小波域的稀疏重叠鬼影图E ₁、E ₂和E ₃。

2-3、确立小波域双层稀疏鬼影模型

步骤2-2得到的小波域图像E ₁、E ₂和E ₃中，信号非常稀疏，每个像素点上只有两层鬼影的重叠，采用双层稀疏鬼影模型来描述E ₁、E ₂和E ₃的每个像素：

                              [1]

公式[1]中为相位因子，G _n1和G _n2分别为每个像素点上需要确定的不同阶的鬼影（不同的阶表示不同的鬼影位置），n₁和n₂表示不同的鬼影阶数。定义为：

                   [2]

公式[2]中d表示每组欠采样数据在PE方向上的偏移量d ₁，d ₂和d ₃，为固定值， n为鬼影阶数。

方程组[1]也可用矩阵表示为：

，                                [3]

公式[3]中，E = (E ₁, E ₂, E ₃) ^T 和G = (G _n1, G _n2) ^T 都是矢量，P为公式[1]中的系数矩阵：

                               [4]

2-4、基于最小平方误差法的重叠鬼影分离

公式[1]或[3]有三个等式，两个未知数G _n1和G _n2，属于过定（Over Determined）问题。对于每一对可能出现的鬼影阶数（n₁, n₂），公式[3]存在最小平方解向量G _LSE ：

                       [5]

公式[5]中上标符号“+”表示复数共轭，上标符号“-1”表示矩阵逆运算。

对于任一采样间隔N，可能出现的（n₁, n₂）组合个数是确定的，通过求取公式[5]中的最小平方误差，选取平方误差最小的（n₁, n₂）组合为像素点正确的（n₁, n₂）组合。每个像素点对应的（n₁, n₂）确定后，按各自取值进行归类，可得到N个不同阶的分离鬼影子图。

2-5、离散小波逆变换

对步骤2-4得到的N个分离鬼影子图，分别进行小波逆变换，得到N个图像域的鬼影子图R₁，R₂，…，R_N。

2-6、多个鬼影子图的配准求和

步骤2-5得到图像域鬼影子图R₁，R₂，…， R_N中，鬼影分布位置不同，通过像素点的移位和对齐来配准多个鬼影子图；对应像素点求和后可得到无重叠鬼影的图像R ₀。

2-7、图像重建

对R ₀进行离散傅立叶逆变换（IDFT），然后用实际采集数据替换IDFT后的数据，且k空间的中心区域用实际全采样数据替代，再经过离散傅立叶正变换DFT，重建出最终的SPEED图像I ₀。

采用本发明方法可以将离散小波变换用于SPEED的数据稀疏表示，通过离散小波变换来提高图像的稀疏性，进而提高SPEED重建图像的质量，同时本发明具有以下特点：

（1）本发明采取简单有规律方式进行数据的欠采样，无需更改MRI硬件，就能提高数据采集的速度，且易于与常规的临床数据采集方式集成。

（2）本发明采用离散小波变换代替常规SPEED的差分变换，对于对比度差、结构较复杂的图像，离散小波变换可以提供比差分变换更稀疏的数据表示，重建出质量更好的图像。

（3）采用离散小波变换替换差分变换后，省去了常规SPEED重建的逆滤波过程，有效地避免了逆滤波时可能出现的极值状态。

（4）本发明没有采用CS那种长时间的非线性迭代来求解，而是基于解析法来求解，重建速度较快。

附图说明

图1是SPEED数据采集方式的示意图；

图2是SPEED重建过程部分中间数据的示意图；

图3是采用本发明进行SPEED数据采集和重建实例的结果图。

具体实施方式

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

本发明主要包括二大步骤：有规律的欠采样数据采集、基于小波域数据稀疏表示的SPEED图像重建。

有规律的欠采样数据采集包括k-空间有规律的欠采样和k-空间中心区域全采样两个步骤。

1-1、k-空间有规则的欠采样

在k-空间的相位编码PE方向每隔N行采集一行k-空间数据，共采集三组（如图1所示），得到三组欠采样的k空间数据S ₁、S ₂ 和S ₃（如图2所示）。用 d ₁, d ₂, d ₃表示每组欠采样数据在PE方向上的偏移量，采样方式用N (d ₁, d ₂, d ₃)表示。

1-2、k-空间中心区域全采样

在信息量较集中的k-空间中心区域采集全部数据，跟据需采集数据矩阵的大小，相应地采集8至32行k-空间中心区域的相位编码数据。

基于小波域数据稀疏表示的SPEED图像重建包括填零傅立叶重建、离散小波变换、确立小波域双层稀疏鬼影模型、基于最小平方误差法的重叠鬼影分离、离散小波逆变换、多个鬼影子图的配准求和、图像重建七个步骤。

2-1、填零傅立叶重建

对于三组欠采样的数据，k-空间中没有采集的数据用0表示，进行常规的填零傅立叶重建，重建后图像分别为I ₁、I ₂和I ₃（如图2所示）。k-空间中每隔N行采集一行数据使得对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多可能有N个重叠的鬼影。例如，当N = 6时， I ₁、I ₂和I ₃上分别有6层重叠的鬼影（如图2所示）。

2-2、离散小波变换

对步骤2-1得到的图像I ₁、I ₂和I ₃分别施加离散小波变换Ψ，得到小波域的稀疏重叠鬼影图E ₁、E ₂和E ₃（如图2所示）。

2-3、确立小波域双层稀疏鬼影模型

步骤2-2得到的小波域图像E ₁、E ₂和E ₃中，信号非常稀疏，每个像素点上只有两层鬼影的重叠，采用双层稀疏鬼影模型来描述E ₁、E ₂和E ₃的每个像素：

                              [1]

公式[1]中为相位因子，G _n1和G _n2分别为每个像素点上需要确定的鬼影（不同的阶表示不同的鬼影位置），n₁和n₂表示不同的鬼影阶数。定义为：

                   [2]

公式[2]中d表示每组欠采样数据在PE方向上的偏移量d ₁，d ₂和d ₃，为固定值， n为鬼影的阶数。

方程组[1]也可用矩阵表示为：

，                                [3]

公式[3]中，E = (E ₁, E ₂, E ₃) ^T 和G = (G _n1, G _n2) ^T 都是矢量，P为双层稀疏鬼影公式[1]中的系数矩阵：

                               [4]

2-4、基于最小平方误差法的重叠鬼影分离

公式[1]或[3]有三个等式，两个未知数G _n1和G _n2，属于过定（Over Determined）问题。对于每一对可能出现的鬼影阶数（n₁, n₂），公式[3]存在最小平方解向量G _LSE ：

                       [5]

公式[5]中上标符号“+”表示复数共轭，上标符号“-1”表示矩阵逆运算。

对于任一采样间隔N，可能出现的（n₁, n₂）组合个数是确定的，通过求取公式[5]中的最小平方误差，选取平方误差最小的（n₁, n₂）组合为像素点正确的（n₁, n₂）组合。每个像素点对应的（n₁, n₂）确定后，按各自取值进行归类，可得到N个不同阶的分离鬼影子图。

2-5、离散小波逆变换

对步骤2-4得到的N个分离鬼影子图，分别进行小波逆变换，得到N个图像域的鬼影子图R₁，R₂，…， R_N（如图2所示，图中N=5）。

2-6、多个鬼影子图的配准求和

步骤2-5得到图像域鬼影子图R₁，R₂，…， R_N中，鬼影分布位置不同，通过像素点的移位和对齐来配准多个鬼影子图；对应像素点求和后可得到无重叠鬼影的图像R ₀（如图2所示）。

2-7、图像重建

对R ₀进行离散傅立叶逆变换（IDFT），然后用实际采集数据替换IDFT后的数据，且k-空间的中心区域用实际全采样数据替代，再经过离散傅立叶正变换DFT，重建出最终的SPEED图像I ₀（如图2所示）。

以下结合人体关节部位MRI图像的欠采集和重建进行实例说明。假设要采集的MRI图像的矩阵大小为kx×ky=256×256。首先以采样方式N (d ₁, d ₂, d ₃)=5 (0, 1, 2) 进行数据采集，在k-空间的相位编码PE方向每隔N=5行采集一行k-空间数据，共采集三组，分别得到欠采样的k-空间数据S ₁、S ₂和S ₃。在信息量集中的k-空间中心区域采集全部数据，共采集32行相位编码数据。接下来，基于小波域数据稀疏表示进行SPEED图像重建。首先对三组欠采样数据分别进行常规的填零傅立叶重建，重建后图像分别为I ₁、I ₂和I ₃。然后对I ₁、I ₂和I ₃分别施加离散小波变换Ψ，得到小波域的稀疏重叠边缘鬼影图E ₁、E ₂和E ₃。然后，确立小波域双层稀疏鬼影模型，并基于最小平方误差法进行重叠鬼影的分离，得到5个分布位置不同的子鬼影图。接着，对得到的5个鬼影子图分别进行小波逆变换，得到5个图像域的鬼影图。通过像素点的移位配准、对齐这5个鬼影，5个鬼影求和后，可得到无重叠鬼影的图像E ₀。最后，对E ₀进行IDFT离散傅立叶逆变换，用实际采集数据替换，再经过离散傅立叶正变换DFT，重建出最终的SPEED图像I ₀，如图3所示。图3（a）为全采样数据重建后的参考图像，图3（b）和3（c）分别为基于小波稀疏和常规差分稀疏的SPEED重建结果。误差分析表明，基于小波稀疏表示的SPEED重建图误差（0.720）小于基于差分的SPEED重建误差（0.819）。本实例在个人PC机（CPU为2.6GHz，4G内存）的Matlab环境下运行时间为1.13秒。本发明采用离散小波变换进行SPEED欠采样数据的稀疏表示，相比SPEED所用的常规差分方法，提高了数据的稀疏性，从而有效地改善了重建图像质量；此外，基于解析求解的过程使得SPEED的重建速度较快。

Claims (1)

1.基于小波域稀疏表示的SPEED快速磁共振成像方法,其特征在于该方法包括两大步骤，分别是有规律的欠采样数据采集、基于小波域数据稀疏表示的SPEED图像重建；

所述的有规律的欠采样数据采集包括k-空间有规律的欠采样和k-空间中心区域全采样两个步骤，具体是：

1-1、k-空间有规则的欠采样

在k-空间的相位编码PE方向每隔N行采集一行k-空间数据，共采集三组，得到三组欠采样的k空间数据S₁、S₂和S₃；用d₁,d₂,d₃表示每组欠采样数据在PE方向上的偏移量，采样方式用N(d₁,d₂,d₃)表示；

1-2、k-空间中心区域全采样

在信息量较集中的k-空间中心区域采集全部数据，跟据需采集数据矩阵的大小，相应地采集8至32行k-空间中心区域的相位编码数据；

所述的基于小波域数据稀疏表示的SPEED图像重建包括填零傅立叶重建、离散小波变换、确立小波域双层稀疏鬼影模型、基于最小平方误差法的重叠鬼影分离、离散小波逆变换、多个鬼影子图的配准求和、图像重建七个步骤；

2-1、填零傅立叶重建

对于三组欠采样的数据，k-空间中没有采集的数据用0表示，进行常规的填零傅立叶重建，重建后图像分别为I₁、I₂和I₃；k-空间中每隔N行采集一行数据使得对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多可能有N个重叠的鬼影；

2-2、离散小波变换

对步骤2-1得到的图像I₁、I₂和I₃分别施加离散小波变换Ψ，得到小波域的稀疏重叠鬼影图E₁、E₂和E₃；

2-3、确立小波域双层稀疏鬼影模型

步骤2-2得到的小波域图像E₁、E₂和E₃中，信号非常稀疏，每个像素点上只有两层鬼影的重叠，采用双层稀疏鬼影模型来描述E₁、E₂和E₃的每个像素：

$E_1=P_{d1}^{n1}{G}_{n1}+P_{d1}^{n2}{G}_{n2}$

$E_2=P_{d2}^{n1}{G}_{n1}+P_{d2}^{n2}{G}_{n2}---\left[1\right]$

$E_3=P_{d3}^{n1}{G}_{n1}+P_{d3}^{n2}{G}_{n2}$

公式[1]中为相位因子，G_n1和G_n2分别为每个像素点上需要确定的不同阶的鬼影，n₁和n₂表示不同的鬼影阶数；定义为：

$P_d^n=e^{i(2\mathrm{\πdn}/N)},n=\mathrm{0,1,2},...,N-1---\left[2\right]$

公式[2]中d表示每组欠采样数据在PE方向上的偏移量d₁，d₂和d₃，为固定值，n为鬼影阶数；

公式[1]也可用矩阵表示为：

E＝PG，  [3]

公式[3]中，E＝(E₁,E₂,E₃)^T和G＝(G_n1,G_n2)^T都是矢量，P为公式[1]中的系数矩阵：

$P=\left[\begin{array}{cc}{P}_{d1}^{n1}& P_{d1}^{n2}\\ P_{d2}^{n1}& P_{d2}^{n2}\\ P_{d3}^{n1}& P_{d3}^{n2}\end{array}\right]---\left[4\right]$

2-4、基于最小平方误差法的重叠鬼影分离

公式[1]或[3]有三个等式，两个未知数G_n1和G_n2，属于过定问题；对于每一对可能出现的鬼影阶数(n₁,n₂)，公式[3]存在最小平方解向量G_LSE：

G_LSE＝[P⁺P]^-1P⁺E  [5]

公式[5]中上标符号“+”表示复数共轭，上标符号“-1”表示矩阵逆运算；

对于任一采样间隔N，可能出现的(n₁,n₂)组合个数是确定的，通过求取公式[5]中的最小平方误差，选取平方误差最小的(n₁,n₂)组合为像素点正确的(n₁,n₂)组合；每个像素点对应的(n₁,n₂)确定后，按各自取值进行归类，可得到N个不同阶的分离鬼影子图；

2-5、离散小波逆变换

对步骤2-4得到的N个分离鬼影子图，分别进行小波逆变换，得到N个图像域的鬼影子图R₁，R₂，…，R_N；

2-6、多个鬼影子图的配准求和

步骤2-5得到图像域鬼影子图R₁，R₂，…，R_N中，鬼影分布位置不同，通过像素点的移位和对齐来配准多个鬼影子图；对应像素点求和后可得到无重叠鬼影的图像R₀；

2-7、图像重建

对R₀进行离散傅立叶逆变换，然后用实际采集数据替换IDFT后的数据，且k空间的中心区域用实际全采样数据替代，再经过离散傅立叶正变换DFT，重建出最终的SPEED图像I₀。