CN110728624B - 一种高分辨率扩散加权图像重建方法 - Google Patents

Abstract

一种高分辨率扩散加权图像重建方法，涉及多激发扩散加权图像的重建方法。包括以下步骤：1)获取多激发傅里叶空间信号，根据对数据的要求设置欠采样方式和实验参数，此时被采样信号点的位置可以确定，进而得到采样模板；2)构建结构化汉克尔矩阵；3)构建基于结构化汉克尔矩阵的扩散加权图像重建模型，建立基于结构化汉克尔矩阵的扩散加权图像的重建模型；4)通过交替方向乘子法求解重建模型获得多激发图像；5)将步骤4)获得的多激发图像合成得到无伪影的高分辨率重建图像。无需额外采集导航回波信息，不仅减少了采样时间，还避免了利用导航回波重建时，图像与导航回波间不匹配的问题。重建出的高分辨率图像是无伪影的。

Description

一种高分辨率扩散加权图像重建方法

技术领域

本发明涉及多激发扩散加权图像的重建方法，尤其是涉及利用无导航回波的多激发扩散加权技术来重建出无伪影高分辨率图像的一种高分辨率扩散加权图像重建方法。

背景技术

扩散加权成像由Stejskal和Tanner于1965年提出(E.O.Stejskal andJ.E.Tanner,"Spin Diffusion Measurements:Spin Echoes in the Presence of aTime-Dependent Field Gradient,"Journal of Chemical Physics,vol.42,no.1,pp.288-292,1965.)，这是一种无入侵的检测组织内水分子的扩散运动的方式。传统扩散加权成像是基于一种单激发的平面回波成像序列(EPI序列)(D.Le Bihan,E.Breton,D.Lallemand,P.Grenier,E.Cabanis,and M.Laval-Jeantet,"MR imaging of intravoxelincoherent motions:application to diffusion and perfusion in neurologicdisorders,"Radiology,vol.161,no.2,pp.401-407,1986.)，这种序列采样速度快，受运动影响小，但是这种成像方式对不均匀场十分敏感，易使得图像在相位编码维变得扭曲，而且受限于采样带宽，难以采集高分辨率的图像(F.Farzaneh,S.J.Riederer,and N.J.Pelc,"Analysis of T2 limitations and off-resonance effects on spatial resolutionand artifacts in echo-planar imaging,"Magnetic Resonance in Medicine,vol.14,no.1,pp.123-139,1990.)。多激发平面回波成像方式减小了回波链长度，抵抗了不均匀场带来的影响(R.Bammer,R.Stollberger,M.Augustin,J.Simbrunner,H.Offenbacher,H.Kooijman,S.Ropele,P.Kapeller,P.Wach,and F.Ebner,"Diffusion-weighted imagingwith navigated interleaved echo-planar imaging and a conventional gradientsystem,"Radiology,vol.211,no.3,pp.799-806,1999.)，但是在扩散梯度的作用下，多次激发所采集到的图像间会有不同的相位，若直接将多激发所采集到的傅里叶空间数据直接合成，将会给图像带来伪影。典型的MUSE重建方法(N.K.Chen,A.Guidon,H.C.Chang,andA.W.Song,"A robust multi-shot scan strategy for high-resolution diffusionweighted MRI enabled by multiplexed sensitivity-encoding(MUSE),"Neuroimage,vol.72,no.2,pp.41-47,2013.)重建出的图像仍然有伪影(如图3所示)。为了解决伪影问题，受启发于约束单张图像的相位平滑，用于单张磁共振图像随机欠采样的重建的LORAKS方法(J.P.Haldar,"Low-rank modeling of local k-space neighborhoods(LORAKS)forconstrained MRI,"IEEE transactions on medical imaging,vol.33,no.3,pp.668-681,2013.)。

低秩重建方法在磁共振波谱欠采样重建中有很好的应用(X.Qu,M.Mayzel,J.F.Cai,Z.Chen,and V.Orekhov,"Accelerated NMR spectroscopy with low-rankreconstruction,"Angewandte Chemie International Edition,vol.54,no.3,pp.852-854,2015.；J.Ying,H.Lu,Q.Wei,J.-F.Cai,D.Guo,J.Wu,Z.Chen,and X.Qu,"Hankelmatrix nuclear norm regularized tensor completion for N-dimensionalexponential signals,"IEEE Transactions on Signal Processing,vol.65,no.14,pp.3702-3717,2017.；J.Ying,J.-F.Cai,D.Guo,G.Tang,Z.Chen,and X.Qu,"Vandermondefactorization of Hankel matrix for complex exponential signal recovery—Application in fast NMR spectroscopy,"IEEE Transactions on Signal Processing,vol.66,no.21,pp.5520-5533,2018.；H.Lu,X.Zhang,T.Qiu,J.Yang,J.Ying,D.Guo,Z.Chen,and X.Qu,"Low rank enhanced matrix recovery of hybrid time andfrequency data in fast magnetic resonance spectroscopy,"IEEE Transactions onBiomedical Engineering,vol.65,no.4,pp.809-820,2017.)。

但是，在多激发扩散加权图像重建中，尚无利用图像的相位平滑特性构建低秩汉克尔矩阵的方法。

发明内容

本发明的目的在于提供一种无导航回波多激发一种高分辨率扩散加权图像重建方法。

本发明包括以下步骤：

1)获取多激发傅里叶空间信号，根据对数据的要求设置欠采样方式和实验参数，此时被采样信号点的位置可以确定，进而得到采样模板；

2)构建结构化汉克尔矩阵；

3)构建基于结构化汉克尔矩阵的扩散加权图像重建模型，建立基于结构化汉克尔矩阵的扩散加权图像的重建模型；

4)通过交替方向乘子法求解重建模型获得多激发图像；

5)将步骤4)获得的多激发图像合成得到无伪影的高分辨率重建图像。

在步骤2)中，所述构建结构化汉克尔矩阵的具体方法可为：令[c,d]为图像坐标，则复数图像ρ[c,d]的相位为ω[c,d]＝∠ρ[c,d]；j表示虚数单位；定义：

h[c,d]＝exp(-jω[c,d]) (1)

其中，

和

分别表示ρ[c,d]和h[c,d]的傅里叶变换；傅里叶空间中半径为R的圆块用

表示，

表示卷积核的范围；

表示图像ρ第n次激发的傅里叶空间数据；Re(·)和Im(·)分别表示取实部和虚部；下标g表示傅里叶空间中取出的第g个圆块，g＝1,2,…,G，G表示傅里叶空间中圆块Λ_R的总数；下标m表示圆块Λ_R中第m个元素，m＝1,2,...,M，M表示圆块Λ_R中元素的个数，由R的大小决定，当R＝1,2,3时，分别对应的M＝5,13,29；定义第n次激发的结构化汉克尔矩阵中每个元素为：

然后，构建具有低秩特性的第n次激发的结构化汉克尔矩阵如下：

其中，

表示将矩阵转换为结构化汉克尔矩阵的算子，N₁×N₂表示输入矩阵的维度，G＝(N₁-2R-1)×(N₂-2R-1)。

在步骤3)中，所述构建基于结构化汉克尔矩阵的扩散加权图像重建模型，建立基于结构化汉克尔矩阵的扩散加权图像的重建模型的具体方法如下：

其中，

表示将矩阵转换为结构化汉克尔矩阵的算子，X＝[X₁,...,X_n,...,X_N]，X_n表示第n次激发的傅里叶空间数据，n＝1,2,...,N，N表示激发次数；C_q表示第q个通道的灵敏度系数矩阵，q＝1,2,...,Q，Q表示通道数；

表示傅里叶变换算子；

表示欠采样并在未采样点填零的算子；Y_q表示采集到的未采样位置已填零的第q个通道的傅里叶空间数据；λ表示正则化参数；||·||_F表示矩阵的弗罗贝尼乌斯范数；||·||*表示矩阵的核范数；右上标*表示伴随算子，如

表示

的伴随算子。

在步骤4)中，所述通过交替方向乘子法求解重建模型获得多激发图像的具体方法可为：在步骤3)的基础上提出基于结构化汉克尔矩阵的扩散加权图像重建模型的求解算法，采用交替乘子法求解式(4)中的最优化问题，引入中间变量

和拉格朗日乘子D，将式(4)中的优化问题转化为以下优化问题：

其中，＜·,·＞表示矩阵在希尔伯特空间的内积，即

表示取复数的实部，trace(·)表示矩阵的迹；参数ρ是正数；

引入交替乘子法对式(5)进行求解，式(5)的优化问题求解可通过交替求解以下问题得到：

求解式(6)中的问题，根据以下公式迭代更新变量：

当达到迭代停止准则时，迭代停止；迭代停止准则设定为达到最大迭代次数或相邻两次迭代中的误差||X^(k+1)-X^(k)||_F小于设置的大于0的阈值η；其中，X^(k+1)，Z^(k+1)和D^(k+1)分别表示变量X，Z和D在第k+1次迭代时的值；

表示奇异值收缩算子；正则化参数λ是正数；右上标H表示共轭转置，如

表示对C_q进行共轭转置操作。

在步骤5)中，所述合成的公式如下：

其中，θ表示重建的图像；

表示求解得到的第n次激发的图像。

本发明首先利用多激发图像的相位平滑特性，构建一种结构化汉克尔矩阵的低秩重建模型，然后通过交替方向乘子法求解重建模型获得多激发图像，最后将得到的多激发图像合成，重建出无伪影的高分辨扩散加权图像。本发明将LORAKS方法用于扩散加权成像中，约束了多激发图像间的相位平滑。这种重建方法无需额外采集导航回波信息，不仅减少了采样时间，还避免了利用导航回波重建时，图像与导航回波间不匹配的问题。本发明重建出的高分辨率图像是无伪影的。

附图说明

图1是4次激发采用的不同的采样模板图。

图2是4次采集的傅里叶空间信号直接合成后做傅里叶反变换得到的没有校正的图像。

图3是使用典型的MUSE方法重建出的有伪影的高分辨率图像。

图4是使用本发明方法重建出的无伪影的高分辨率图像。

具体实施方式

本发明实例是一个多激发傅里叶空间信号高分辨率扩散加权重建的具体过程，是对本发明所提出方法的详细描述。

具体实施过程如下：

第一步：获取多激发傅里叶空间信号

本实施案例中，包含了32通道，4次激发所采集到的傅里叶空间数据，其采样模板如图1所示，图中白色的点为采样点，表示该位置对应的数据被采样，黑色的点表示没有采样到的点，该位置对应的数据丢失；每次采集傅里叶空间1/4的数据量，并且使用了部分傅里叶采集技术，系数为0.56，因此，实际每次激发采集数据量约为14％。对4次激发的数据直接合并后做傅里叶反变换得到没有校正的图像，如图2所示。

第二步：构建结构化汉克尔矩阵

令[c,d]为图像坐标，则复数图像ρ[c,d]的相位为ω[c,d]＝∠ρ[c,d]；j表示虚数单位。

定义：

h[c,d]＝exp(-jω[c,d]), (1)

和

分别表示ρ[c,d]和h[c,d]的傅里叶变换；傅里叶空间中半径为R的圆块用

表示，

表示卷积核的范围，本方法中R＝3；

表示图像ρ第n次激发的傅里叶空间数据；Re(·)和Im(·)分别表示取实部和虚部；下标g表示傅里叶空间中取出的第g个圆块，g＝1,2,...,G，G表示傅里叶空间中圆块Λ_R的总数；下标m表示圆块Λ_R中第m个元素，m＝1,2,...,M，M表示圆块Λ_R中元素的个数，本方法中M＝29。定义第n次激发的结构化汉克尔矩阵中每个元素为：

那么可以构建具有低秩特性的第n次激发的结构化汉克尔矩阵如下：

其中，

表示将矩阵转换为结构化汉克尔矩阵的算子，N₁×N₂表示输入矩阵的维度，本方法中N₁＝248，N₂＝244，G＝(N₁-2R-1)×(N₂-2R-1)＝57117。

第三步：构建基于结构化汉克尔矩阵的扩散加权图像重建模型

构建基于结构化汉克尔矩阵的扩散加权图像重建模型如下：

其中，

表示将矩阵转换为结构化汉克尔矩阵的算子，X＝[X₁,...,X_n,...,X_N]，X_n表示第n次激发的傅里叶空间数据，n＝1,2,...,N，N表示激发次数，在本方法中N＝4；C_q表示第q个通道的灵敏度系数矩阵，q＝1,2,...,Q，Q表示通道数，在本方法中Q＝32；

表示傅里叶变换算子；

表示欠采样并在未采样点填零的算子；Y_q表示采集到的未采样位置已填零的第q个通道的傅里叶空间数据；λ表示正则化参数，在本方法中λ＝10³；||·||_F表示矩阵的弗罗贝尼乌斯范数；||·||*表示矩阵的核范数；右上标*表示伴随算子，如

表示

的伴随算子。

第四步：提出构建基于结构化汉克尔矩阵的扩散加权图像重建模型的求解算法

采用交替乘子法求解公式(4)中的最优化问题。引入中间变量

和拉格朗日乘子D，公式(4)中的优化问题转化为以下优化问题：

其中＜·,·＞表示矩阵在希尔伯特空间的内积，即

表示取复数的实部，trace(·)表示矩阵的迹；在本方法中参数ρ＝1。

引入交替乘子法对公式(5)进行求解，公式(5)的优化问题求解可以通过交替求解以下问题得到：

求解(6)中的问题，根据以下公式迭代更新变量：

当达到迭代停止准则时，迭代停止。迭代停止准则设定为达到最大迭代次数或相邻两次迭代中的误差||X^(k+1)-X^(k)||_F小于设置的大于0的阈值η，本方法中η＝10^-6。其中，X^{(k +1)}，Z^(k+1)和D^(k+1)分别表示变量X，Z和D在第k+1次迭代时的值；本方法初始化时，设置k＝0，Z⁽⁰⁾为随机矩阵，D⁽⁰⁾为各位置全为1的矩阵；

表示奇异值收缩算子；右上标H表示共轭转置，如

表示对C_q进行共轭转置操作。

第五步：将求解得到的多激发图像合成无伪影的高分辨率重建图像

将求解得到的多次激发图像合成无伪影的高分辨率重建图像。将第四步求解得到的多次激发图像，对其进行合成即可得到无伪影的高分辨率重建图像如图4所示，合成公式如下：

其中，θ表示重建的图像；

表示求解得到的第n次激发的图像。

从图3可看出，典型的MUSE方法的重建图像有明显伪影。而从图4可看出，本发明重建出的高分辨率图像是无伪影的。

Claims (3)

1.一种高分辨率扩散加权图像重建方法，其特征在于包括以下步骤：

1)获取多激发傅里叶空间信号，根据对数据的要求设置欠采样方式和实验参数，此时被采样信号点的位置可以确定，进而得到采样模板；

2)构建结构化汉克尔矩阵；

3)构建基于结构化汉克尔矩阵的扩散加权图像重建模型，建立基于结构化汉克尔矩阵的扩散加权图像的重建模型，具体方法如下：

其中，

表示将矩阵转换为结构化汉克尔矩阵的算子，X＝[X₁,...,X_n,...,X_N]，X_n表示第n次激发的傅里叶空间数据，n＝1,2,...,N，N表示激发次数；C_q表示第q个通道的灵敏度系数矩阵，q＝1,2,...,Q，Q表示通道数；

表示傅里叶变换算子；u表示欠采样并在未采样点填零的算子；Y_q表示采集到的未采样位置已填零的第q个通道的傅里叶空间数据；λ表示正则化参数；||·||_F表示矩阵的弗罗贝尼乌斯范数；||·||_*表示矩阵的核范数；右上标*表示伴随算子，如

表示

的伴随算子；

4)通过交替方向乘子法求解重建模型获得多激发图像，具体方法为：在步骤3)的基础上提出基于结构化汉克尔矩阵的扩散加权图像重建模型的求解算法，采用交替乘子法求解式(4)中的最优化问题，引入中间变量

和拉格朗日乘子D，将式(4)中的优化问题转化为以下优化问题：

其中，＜·,·＞表示矩阵在希尔伯特空间的内积，即

表示取复数的实部，trace(·)表示矩阵的迹；参数ρ是正数；

引入交替乘子法对式(5)进行求解，式(5)的优化问题求解可通过交替求解以下问题得到：

求解式(6)中的问题，根据以下公式迭代更新变量：

当达到迭代停止准则时，迭代停止；迭代停止准则设定为达到最大迭代次数或相邻两次迭代中的误差||X^(k+1)-X^(k)||_F小于设置的大于0的阈值η；其中，X^(k+1)，Z^(k+1)和D^(k+1)分别表示变量X，Z和D在第k+1次迭代时的值；

表示奇异值收缩算子；正则化参数λ是正数；右上标H表示共轭转置，如

表示对C_q进行共轭转置操作；

5)将步骤4)获得的多激发图像合成得到无伪影的高分辨率重建图像。

2.如权利要求1所述一种高分辨率扩散加权图像重建方法，其特征在于在步骤2)中，所述构建结构化汉克尔矩阵的具体方法为：令[c,d]为图像坐标，则复数图像x[c,d]的相位为ω[c,d]；j表示虚数单位；定义：

h[c,d]＝exp(-jω[c,d]) (1)

对x[c,d]和h[c,d]做傅里叶变换，分别记为X[e,f]和H[p,q]；将X[e,f]中的第n次激发的傅里叶空间数据取出，记为X_n[e,f]，n＝1,2,.. ,N，其中N表示总共激发次数；傅里叶空间中半径为R的圆块用

表示，

表示卷积核的范围；Re(·)和Im(·)分别表示取实部和虚部；下标g表示傅里叶空间中取出的第g个圆块，g＝1,2,...,G，G表示傅里叶空间中圆块Λ_R的总数；下标m表示圆块Λ_R中第m个元素，m＝1,2,...,M，M表示圆块Λ_R中元素的个数，由R的大小决定，当R＝1,2,3时，分别对应的M＝5,13,29；定义第n次激发的结构化汉克尔矩阵中每个元素为：

然后，构建具有低秩特性的第n次激发的结构化汉克尔矩阵如下：

再将每次激发得到的数据串联为矩阵S，即S＝[S₁,...,S_n,...,S_N]；

为便于后面的权利要求描述，将上述从磁共振图像的傅里叶空间数据X中构建结构化汉克尔矩阵S的操作记为算子

即

3.如权利要求1所述一种高分辨率扩散加权图像重建方法，其特征在于在步骤5)中，所述合成的公式如下：

其中，θ表示重建的图像；

表示求解得到的第n次激发的图像。

Images