CN103728581A - 基于离散余弦变换的speed快速磁共振成像方法 - Google Patents

Abstract

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

Description

基于离散余弦变换的SPEED快速磁共振成像方法

技术领域

本发明属于磁共振的图像成像领域，涉及一种基于离散余弦变换的SPEED快速磁共振成像方法。

背景技术

由于磁共振成像(Magnetic Resonance Imaging,MRI)对机体没有不良影响和较强的软组织分辨能力等优点，目前已在临床疾病检测中得到广泛应用。但是，MRI在临床应用中还常受到数据采集时间过长的限制。研究人员已通过提高MRI硬件性能、采用高效的k-空间（频率空间）数据采集轨迹、研制快速序列和并行数据采集等方式来极大地提高了MRI数据采集的速度，但是很多应用中，还是不能完全满足临床对快速数据采集的需求，例如脑功能成像，三维高分辨率磁敏感加权成像和动态数据采集等场合。MRI的数据采集是在k-空间进行的，在k-空间的相位编码（Phase Encoding,PE）方向上减少数据采集点数，可以有效地缩短数据采集时间。近年来，如何在不降低图像质量条件下减少数据采集总量成为进一步快速磁共振成像的研究热点，如压缩感知（Compressed Sensing,CS）技术，SPEED（SkippedPhase Encoding and Edge Deghosting）技术等。

CS技术利用信号自身或其在变换域中的稀疏性，通过随机采集少量的k-空间数据点和非线性迭代优化等算法来重建出高质量的图像（M Lustig,et.al.,Sparse MRI:The Application of Compressed Sensing for Rapid MRImaging,Magnetic Resonance in Medicine,58:1182-1195,2007）。然而，由于受到硬件扫描系统有规律往返轨迹的限制，在实际应用中难以满足CS技术要求完全随机数据采集，不能满足压缩感知重建时需要的不相关性（Incoherent），从而影响了图像的质量；此外，基于非线性迭代的重建方式使得CS技术的图像重建时间很长，制约了其临床普及应用。

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

目前已申请的关于快速成像方面的MRI专利有：快速磁共振成像方法及系统（申请号：201210013858.2），一种基于部分回波压缩感知的快速磁共振血管成像方法（申请号：201010272089.8）；一种基于CS压缩感知技术的快速磁共振成像方法（申请号：201010272095.4）等。基于小波域稀疏表示的SPEED快速磁共振成像方法（申请号：2013102071971），提出基于小波域的数据稀疏特性来提高常规的基于差分变换的SPEED成像质量。但是目前还未能查询到任何关于基于离散余弦变换的SPEED快速成像方面的授权发明专利或申请。

国内外已发表的SPEED快速成像方面的文章有：2013年，金朝阳和向清三提出了通用G-SPEED（General-SPEED）采样方法（Jin Z,Xiang QS.Accelerated MRI by SPEED with Generalized Sampling Schemes.MagneticResonance in Medicine.70:1674-1681，2013），突破了传统SPEED方法的采样间隔周期N必须是质数（例如：N=5、7、11）的限制，通过秩判据的方式，使得N不但可为质数，也可为合数（例如：N=2、4、6、8、9）。基于MRA数据本身就非常稀疏的特性，常征和向清三将SPEED的双层模型简化到单层模型（Chang Z and Xiang QS.Simplified skipped phaseencoding and edge deghosting(SPEED)for imaging sparse objects withapplications to MRA.Med Phys.34:3173-3182,2007），提出了S-SPEED（Simplified-SPEED）算法，适用于数据本身就非常稀疏的场合，例如暗背景亮信号的MRA应用（Chang Z,Xiang QS,Shen H and Yin FF.Acceleratingnon-contrast-enhanced MR angiography with inflow inversion recoveryimaging by skipped phase encoding and edge deghosting(SPEED).Journal ofMagnetic Resonance Imaging.31:757-765,2010）。2006年，常征和向清三将SPEED算法与并行成像技术进一步结合，提出了SPEED-ACE成像法（Chang Z and Xiang QS.Highly accelerated MRI by skipped phase encodingand edge deghosting with array coil enhancement(SPEED-ACE).Med Phys.33:3758-3766,2006），利用多个采集线圈共同采集k-空间欠采样数据，进一步提高成像速度。2009年，常征等人提出EMA-SPEED（Efficient MultipleAcquisition by SPEED）算法（Chang Z,Xiang QS,Ji J,and Yin FF.Efficientmultiple acquisitions by skipped phase encoding and edge deghosting(SPEED)using shared spatial information.Magnetic Resonance in Medicine.61:229-233,2009），通过共享多个采集间的相似空间信息进一步缩短了SPEED的数据采集时间，从而可获得比单次采集更高的加速比。

以上发表的关于SPEED快速磁共振成像方面的文章或已申请的发明专利，重建时是基于差分或小波变换的方式进行数据的稀疏表示，还未公开过任何基于离散余弦变换的SPEED快速成像方法。

发明内容

本发明针对现有SPEED技术的不足，将离散余弦变换用于SPEED快速磁共振成像中，提供了一种新的SPEED数据稀疏表示方法，提高了SPEED快速成像技术的图像质量。本发明主要包括八个步骤：k空间数据采集、填零傅立叶重建、离散余弦变换、确立离散余弦变换域双层稀疏鬼影模型、基于最小平方误差法（Least Square Error，LSE）的重叠鬼影分离、离散余弦逆变换、多个鬼影子图的配准求和、图像重建。

1、k空间数据采集

在k-空间的相位编码PE方向每隔N行采集一行数据，共采集三组，分别用S₁、S₂和S₃表示。用d₁,d₂,d₃表示每组欠采样数据在PE方向上的偏移量，采样方式用N(d₁,d₂,d₃)表示。根据图像大小，在k-空间中心区域PE方向分别采集16至64行数据。

2、填零傅立叶重建

对于三组欠采样的数据，k-空间中没有采集的数据用0表示，进行常规的填零傅立叶重建，重建后图像分别为I₁、I₂和I₃。

k-空间中每隔N行采集一行数据使得对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多可能有N个重叠的鬼影。例如，当N=5时，I₁、I₂和I₃上分别有5层重叠的鬼影。

3、离散余弦变换

对步骤2得到的图像I₁、I₂和I₃分成若干个大小为N_x×N_y的图像子块，对各个图像子块的实部和虚部分别施加离散余弦变换（DCT），例如，对任一图像子块I的DCT可表示为：

$E(k_x,k_y)={a_k}_x{a}_{k_y}\underset{y=0}{\overset{N_y-1}{\mathrm{\Σ}}}\underset{x=0}{\overset{N_x-1}{\mathrm{\Σ}}}I(x,y)\cos \left(\frac{\mathrm{\π}(2y+1)k_y}{{2N}_y}\right)\cos \left(\frac{\mathrm{\π}(2x+1)k_x}{{2N}_x}\right)---\left[1\right]$

其中， ${a_k}_x=\left\{\begin{array}{cc}\frac{1}{N_x},& k_x=0\\ \sqrt{\frac{2}{N_x}},& 1\≤k_x\≤N_x-1\end{array}\right.,$ ${a_k}_y=\left\{\begin{array}{cc}\frac{1}{N_y},& k_y=0\\ \sqrt{\frac{2}{N_y}},& 1\≤k_y\≤N_y-1\end{array}\right.$

并且，k_x＝0,1,2,3,…N_x-1,k_y＝0,1,2,3,…N_y-1，x和y为图像子块的像素点坐标，k_x和k_y为离散余弦域对应的像素点坐标。

每个图像子块经过离散余弦变换后得到大小为N_x×N_y的变换系数矩阵，将变换系数矩阵中对应于（0，0）位置的系数值设为零，从而去掉低频信号，得到离散余弦域的稀疏鬼影图E₁、E₂和E₃。

4、确立离散余弦变换域双层稀疏鬼影模型

步骤3得到的离散余弦变换域稀疏鬼影图E₁、E₂和E₃中，每个像素点上只有两层鬼影的重叠，可采用双层稀疏鬼影模型来描述E₁、E₂和E₃的每个像素：

$\begin{array}{c}{E}_1=P_{d1}^{n1}{G}_{n1}+P_{d1}^{n2}{G}_{n2}\\ E_2=P_{d2}^{n1}{G}_{n1}+P_{d2}^{n2}{G}_{n2}\\ E_3=P_{d3}^{n1}{G}_{n1}+P_{d3}^{n2}{G}_{n2}\end{array}---\left[2\right]$

公式[2]中

为相位因子，G_n1和G_n2分别为每个像素点上需要确定的不同阶的鬼影（不同的阶表示不同的鬼影位置），n₁和n₂表示不同的鬼影阶数。

定义为：

$\begin{array}{cc}{P}_d^n=e^{i(2\mathrm{\πdn}/N)},& n=\mathrm{0,1,2},\·\·\·,N-1---\left[3\right]\end{array}$

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

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

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

$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[5\right]$

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

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

G_LSE＝[P+P]^-1P+E   [6]公式[6]中上标符号“+”表示复数共轭，上标符号“-1”表示矩阵逆运算。

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

6、离散余弦逆变换

对步骤5得到的N个分离鬼影子图，按子块的实部和虚部进行离散余弦逆变换（IDCT），再将每个鬼影子图对应的逆变换子块系数矩阵分别合成一个复数图像，得到N个图像域的鬼影子图R₁，R₂，…，R_N。

7、多个鬼影子图的配准求和

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

8、图像重建

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

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

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

（2）本发明采用离散余弦变换代替常规SPEED的差分变换，离散余弦变换具有很强的"能量集中"特性：大多数的自然信号（包括声音和图像）的能量都集中在离散余弦变换后的低频部分，通过分离低频信号，离散余弦变换可以提供比差分变换更稀疏的数据表示，重建出质量更好的图像。

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

附图说明

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

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

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

具体实施方式

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

本发明主要包括八个步骤：k空间数据采集、填零傅立叶重建、离散余弦变换、确立离散余弦变换域双层稀疏鬼影模型、基于最小平方误差法LSE的重叠鬼影分离、离散余弦逆变换、多个鬼影子图的配准求和、图像重建。

1、k空间数据采集

如图1所示，在k-空间的相位编码PE方向每隔N行采集一行数据，共采集三组，分别用S₁、S₂和S₃表示。用d₁,d₂,d₃表示每组欠采样数据在PE方向上的偏移量，采样方式用N(d₁,d₂,d₃)表示。根据图像大小，在k-空间中心区域PE方向分别采集16至64行数据。

2、填零傅立叶重建

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

3、离散余弦变换

对步骤2得到的图像I1、I2和I3分成多个大小为N_x×N_y的图像子块，对各个图像子块的实部和虚部分别施加离散余弦变换（DCT），例如，对任一图像子块I的DCT可表示为：

$E(k_x,k_y)={a_k}_x{a}_{k_y}\underset{y=0}{\overset{N_y-1}{\mathrm{\Σ}}}\underset{x=0}{\overset{N_x-1}{\mathrm{\Σ}}}I(x,y)\cos \left(\frac{\mathrm{\π}(2y+1)k_y}{{2N}_y}\right)\cos \left(\frac{\mathrm{\π}(2x+1)k_x}{{2N}_x}\right)---\left[1\right]$

其中， ${a_k}_x=\left\{\begin{array}{cc}\frac{1}{N_x},& k_x=0\\ \sqrt{\frac{2}{N_x}},& 1\≤k_x\≤N_x-1\end{array}\right.,$ ${a_k}_y=\left\{\begin{array}{cc}\frac{1}{N_y},& k_y=0\\ \sqrt{\frac{2}{N_y}},& 1\≤k_y\≤N_y-1\end{array}\right.$

并且，k_x＝0,1,2,3,…N_x-1,k_y＝0,1,2,3,…N_y-1，x和y为图像子块的像素点坐标，k_x和k_y为离散余弦域对应的像素点坐标。

每个图像子块经过离散余弦变换后得到大小为N_x×N_y的变换系数矩阵，将变换系数矩阵中对应于（0，0）位置的系数值设为零，从而去掉低频信号，得到离散余弦域的稀疏鬼影图E₁、E₂和E₃，如图2所示。

4、确立离散余弦变换域双层稀疏鬼影模型

步骤3得到的离散余弦变换域稀疏鬼影图E₁、E₂和E₃中，每个像素点上只有两层鬼影的重叠，可采用双层稀疏鬼影模型来描述E₁、E₂和E₃的每个像素：

$\begin{array}{c}{E}_1=P_{d1}^{n1}{G}_{n1}+P_{d1}^{n2}{G}_{n2}\\ E_2=P_{d2}^{n1}{G}_{n1}+P_{d2}^{n2}{G}_{n2}\\ E_3=P_{d3}^{n1}{G}_{n1}+P_{d3}^{n2}{G}_{n2}\end{array}---\left[2\right]$

公式[2]中

为相位因子，G_n1和G_n2分别为每个像素点上需要确定的不同阶的鬼影（不同的阶表示不同的鬼影位置），n₁和n₂表示不同的鬼影阶数。定义为：

$\begin{array}{cc}{P}_d^n=e^{i(2\mathrm{\πdn}/N)},& n=\mathrm{0,1,2},\·\·\·,N-1---\left[3\right]\end{array}$

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

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

E＝PG，   [4]

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

$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[5\right]$

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

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

G_LSE＝[P+P]^-1P+E   [6]公式[6]中上标符号“+”表示复数共轭，上标符号“-1”表示矩阵逆运算。

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

6、离散余弦逆变换

对步骤5得到的N个分离鬼影子图，按子块的实部和虚部进行离散余弦逆变换（IDCT），再将每个鬼影子图对应的逆变换子块系数矩阵分别合成一个复数图像，得到N个图像域的鬼影子图R₁，R₂，…，R_N，如图2所示。

7、多个鬼影子图的配准求和

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

8、图像重建

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

以下结合人体脑部MRI图像的欠采集和重建进行实例说明。假设要采集的MRI图像的矩阵大小为k_x×k_y=256×256。首先以采样方式N(d₁,d₂,d₃)=5(0,1,2)进行数据采集，在k-空间的相位编码PE方向每隔N=5行采集一行k-空间数据，共采集三组，分别得到欠采样的k-空间数据S₁、S₂和S₃。在信息量集中的k-空间中心区域采集全部数据，共采集32行相位编码数据。接下来，对三组欠采样数据分别进行常规的填零傅立叶重建，重建后图像分别为I₁、I₂和I₃。然后对I₁、I₂和I₃的各个图像子块的实部和虚部分别施加离散余弦变换DCT，并将各个字块的离散余弦变换低频系数置为零，得到离散余弦变换域的稀疏重叠边缘鬼影图E₁、E₂和E₃。然后，确立离散余弦变换域双层稀疏鬼影模型，并基于最小平方误差法进行重叠鬼影的分离，得到5个分布位置不同的子鬼影图。接着，对得到的5个鬼影子图的实部和虚部分别进行基于子块的离散余弦逆变换，子块的大小设为8×8，合成对应的复数图像，得到5个图像域的鬼影图。通过像素点的移位配准、对齐这5个鬼影，5个鬼影求和后，可得到无重叠鬼影的图像E₀。最后，对E₀进行IDFT离散傅立叶逆变换，用实际采集数据替换，再经过离散傅立叶正变换DFT，重建出最终的SPEED图像I₀，如图3所示。图3（a）为全采样数据重建后的参考图像，图3（b）和3（c）分别为基于离散余弦变换和常规差分变换的SPEED重建结果。误差分析表明，基于离散余弦变换的SPEED重建图误差（3.94e-4）小于基于差分的SPEED重建误差（4.62e-4）。本发明采用离散小波变换进行SPEED欠采样数据的稀疏表示，相比SPEED所用的常规差分方法，提高了数据的稀疏性，从而有效地改善了重建图像质量。

Claims (1)

1.基于离散余弦变换的SPEED快速磁共振成像方法,其特征在于该方法包括以下步骤：

步骤1、k空间数据采集

在k-空间的相位编码PE方向每隔N行采集一行数据，共采集三组，分别用S₁、S₂和S₃表示；用d₁,d₂,d₃表示每组欠采样数据在PE方向上的偏移量，采样方式用N(d₁,d₂,d₃)表示；根据图像大小，在k-空间中心区域PE方向分别采集16至64行数据；

步骤2、填零傅立叶重建

对于三组欠采样的数据，k-空间中没有采集的数据用0表示，进行常规的填零傅立叶重建，重建后图像分别为I₁、I₂和I₃；

k-空间中每隔N行采集一行数据使得对应的填零傅立叶重建图像中有N层重叠的鬼影，每个像素点上最多可能有N个重叠的鬼影；

步骤3、离散余弦变换

对步骤2得到的图像I₁、I₂和I₃分成若干个大小为N_x×N_y的图像子块，对各个图像子块的实部和虚部分别施加离散余弦变换，对任一图像子块I的离散余弦变换可表示为：

$E(k_x,k_y)={a_k}_x{a}_{k_y}\underset{y=0}{\overset{N_y-1}{\mathrm{\Σ}}}\underset{x=0}{\overset{N_x-1}{\mathrm{\Σ}}}I(x,y)\cos \left(\frac{\mathrm{\π}(2y+1)k_y}{{2N}_y}\right)\cos \left(\frac{\mathrm{\π}(2x+1)k_x}{{2N}_x}\right)---\left[1\right]$ 其中， ${a_k}_x=\left\{\begin{array}{cc}\frac{1}{N_x},& k_x=0\\ \sqrt{\frac{2}{N_x}},& 1\≤k_x\≤N_x-1\end{array}\right.,$ ${a_k}_y=\left\{\begin{array}{cc}\frac{1}{N_y},& k_y=0\\ \sqrt{\frac{2}{N_y}},& 1\≤k_y\≤N_y-1\end{array}\right.$

并且，k_x＝0,1,2,3,…N_x-¹,k_y＝0,1,2,3,…N_y-¹，x和y为图像子块的像素点坐标，k_x和k_y为离散余弦域对应的像素点坐标；

每个图像子块经过离散余弦变换后得到大小为N_x×N_y的变换系数矩阵，将变换系数矩阵中对应于（0，0）位置的系数值设为零，从而去掉低频信号，得到离散余弦域的稀疏鬼影图E₁、E₂和E₃；

步骤4、确立离散余弦变换域双层稀疏鬼影模型

步骤3得到的离散余弦变换域稀疏鬼影图E₁、E₂和E₃中，每个像素点上只有两层鬼影的重叠，可采用双层稀疏鬼影模型来描述E₁、E₂和E₃的每个像素：

$\begin{array}{c}{E}_1=P_{d1}^{n1}{G}_{n1}+P_{d1}^{n2}{G}_{n2}\\ E_2=P_{d2}^{n1}{G}_{n1}+P_{d2}^{n2}{G}_{n2}\\ E_3=P_{d3}^{n1}{G}_{n1}+P_{d3}^{n2}{G}_{n2}\end{array}---\left[2\right]$

公式[2]中

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

定义为：

$\begin{array}{cc}{P}_d^n=e^{i(2\mathrm{\πdn}/N)},& n=\mathrm{0,1,2},\·\·\·,N-1---\left[3\right]\end{array}$

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

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

E＝PG，   [4]

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

$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[5\right]$

步骤5、基于最小平方误差法的重叠鬼影分离

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

G_LSE＝[P+P]^-1P+E   [6]公式[6]中上标符号“+”表示复数共轭，上标符号“-1”表示矩阵逆运算；

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

步骤6、离散余弦逆变换

对步骤5得到的N个分离鬼影子图，按子块的实部和虚部进行离散余弦逆变换，再将每个鬼影子图对应的逆变换子块系数矩阵分别合成一个复数图像，得到N个图像域的鬼影子图R₁，R₂，…，R_N；

步骤7、多个鬼影子图的配准求和

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

步骤8、图像重建

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

Images