CN104323776A - 脑功能磁共振成像方法和系统 - Google Patents

Abstract

本发明提供了一种新的脑功能磁共振成像方法和系统，其方法包括：获取基于多通道并行且稀疏采样获得的每一个相控阵线圈对应的动态数据图像和第一参考图像；根据所述第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像；根据所述基准图像，获得高阶广义级模型的基函数；利用L2范数对高阶广义级模型的级系数进行约束，基于所述动态数据图像和所述基函数估计获得所述级系数；根据所述级系数的估计值和所述基准图像，基于高阶广义级模型重建出每个通道各个时间点对应的中间图像；基于相控阵线圈的敏感度矩阵，对所述中间图像进行并行重建，获得重建后的动态脑功能磁共振图像。其解决了现有技术成像技术中存在图像几何畸变以及奈奎斯采样伪影问题。

Description

脑功能磁共振成像方法和系统

技术领域

本发明涉及脑功能磁共振成像技术，特别是涉及一种基于广义级模型的脑功能磁共振成像方法和系统。

背景技术

脑功能磁共振成像(function magnetic resonance imaging，fMRI)是在大脑神经活动情况下测量大脑的血液动力学反应。目前常用的EPI(Echo-planar imaging)成像能够在几百毫秒内采集一个层面，空间分辨率可以达到2-3毫米，然而随着功能成像研究的不断深入和发展，基于EPI的脑功能成像由于其本身硬件条件的限制，很难进一步的提高空间和时间分辨率。

BOLD效应fMRI是基于神经元功能活动对局部氧耗量和脑血流影响程度不匹配所导致的局部磁场性质变化的原理。由于这种信号强度较弱，且神经元活动引起局部力流增加是短暂的，普通的MRI成像速度较慢，很难捕捉到神经电活动引起的这种瞬时变化，目前各种各样的fMRI成像方法相继提出来进一步的提高分辨率和信号强度。例如，并行成像技术中的SMASH(SMASH技术主要利用相阵控线圈敏感度的线性组合代替由梯度产生的相位编码的空间调节，从而减少了图像采集时间。)等，可以联合EPI技术，在保证一定的时间分辨率的前提下提高空间分辨率，由于读出时间的减少，减少了磁敏感度伪影，并且EPI噪声也由于减少了梯度间的频繁切换而随之减少。然而在高场fMRI情况下，T2*衰减时间很短，在单次激发很难在保证图像质量的前提下采集大量的相位编码线，为了解决这些问题，研究者们提出采用分段的k空间获取，部分傅里叶采集(partial Fourier sampling)，外部容积压制(outer volume suppression)和内部容积选择(inner volume selection)等。除此之外，基于单击发多回波的EPI和非笛卡尔采样原理提出了多种快速磁共振BOLD功能成像序列。如2008年Rabrait等提出的单激发体回波成像(single shot echo-volumar imaging,EVI)技术、2009年Grotz等提出的降采样投影成像技术(projection imaging)、2010年Feinberg等提出的multiplexed EPI技术、2011年Zahneisen等提出的三维rosette轨迹采样成像以及2012年Stefan等提出的multi-slab echo-volumar imaging技术等。这些技术能将BOLD功能成像的时间分辨率从常规的2～3s提高至几百毫秒。然而由于这些方法均以EPI序列为基础，依然存在着EPI固有的成像缺点，需要矫正Nyquist伪影和几何形变，其时间分辨率的提高是以增加图像几何畸变或牺牲空间分辨率为代价换来的。

在现有的磁共振硬件条件下，目前脑功能磁共振成像方法主要是采用最快的EPI序列采集方式进行成像，由于该序列本身对磁场硬件梯度的要求过高，以及图像空间分辨率较低，容易出现图像几何畸变以及Nyquist(奈奎斯采样)伪影。

发明内容

基于此，有必要针对现有技术成像技术中存在图像几何畸变以及奈奎斯采样伪影问题，提供一种新的脑功能磁共振成像方法和系统。

本发明的一种脑功能磁共振成像方法，其包括：

获取基于多通道并行且稀疏采样获得的每一个相控阵线圈对应的动态数据图像和第一参考图像；

根据所述第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像；

根据所述基准图像，获得高阶广义级模型的基函数；

利用L₂范数对高阶广义级模型的级系数进行约束，基于所述动态数据图像和所述基函数估计获得所述级系数；

根据所述级系数的估计值和所述基准图像，基于高阶广义级模型重建出每个通道各个时间点对应的中间图像；

基于相控阵线圈的敏感度矩阵，对所述中间图像进行并行重建，获得重建后的动态脑功能磁共振图像。

在其中一个实施例中，所述第一参考图像包括：针对每一个相控阵线圈在采集所述动态数据图像前后分别采集的至少两个K空间欠采样图像。

在其中一个实施例中，所述动态数据图像为K空间欠采样、且边缘未采集的图像。

在其中一个实施例中，所述根据所述基准图像获得高阶广义级模型的基函数的过程包括：

给定控制参数μ；

调用下述公式(1)计算获得所述高阶广义级模型的基函数：

其中，表示高阶广义级模型的基函数；表示第l个相控阵线圈在第q个时间点对应的基准图像；l＝1,...,L，L表示相控阵线圈的个数，每一个相控阵线圈采集的数据对应一个输出通道；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

在其中一个实施例中，所述根据所述级系数、基于高阶广义级模型重建出每个通道各个时间点对应的中间图像的过程包括：

获取控制参数μ；

调用下述公式(3)计算每个通道各个时间点对应的中间图像：

${\widehat{\mathrm{\ρ}}}_{q,l}\left(x\right)=\underset{n=-M/2}{\overset{M/2-1}{\mathrm{\Σ}}}{c}_n({\mathrm{\ρ}}_{q,l}^{\mathrm{ref}}\left(x\right)+\mathrm{\μ})e^{i2\mathrm{\πn\Λkx}}$                    公式(3)

其中，c_n表示所述级系数；表示第l个相控阵线圈在第q个时间点对应的中间图像；M表示采集所述动态数据图像的相位编码数，n表示相位编码数变量，其取值范围为 表示第l个相控阵线圈在第q个时间点对应的基准图像；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

在其中一个实施例中，所述利用L₂范数对高阶广义级模型的级系数进行约束、基于所述动态数据图像和所述基函数估计获得所述级系数的过程包括：

根据获得的基函数生成矩阵向量；

依据欠采样因子构成的欠采样矩阵乘以所述矩阵向量获得中间矩阵向量；

调用下述公式(4)计算所述级系数：

${\hat{c}}_n={(A^{H}A+\mathrm{\λI})}^{-1}{A}^{H}D$                    公式(4)

其中，表示所述级系数的估计值；A表示所述中间矩阵向量，A^H表示矩阵A的共轭矩阵；λ表示正则化参数；D表示所述动态数据图像，即针对第l个相控阵线圈的T个采样时间帧对应的K空间欠采样图像D_1,l(k),...,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数；是单位算子。

在其中一个实施例中，所述利用灵敏度编码技术对所述中间图像进行并行重建的过程包括：

根据针对每一个相控阵线圈采集的一个或多个全采样的参考图像，获得敏感度矩阵；

结合所述敏感度矩阵和所述中间图像，采用最小二乘法估计每个时间点对应的初始重建图像；

基于高阶广义级模型，利用L₂范数对所述初始重建图像进行约束，估计获得每个时间点对应的重建后的脑功能磁共振图像；

汇总各个时间点对应的重建后的脑功能磁共振图像，形成所述重建后的动态脑功能磁共振图像。

在其中一个实施例中，所述利用L₂范数对所述初始重建图像进行约束、估计获得每个时间点对应的重建后的脑功能磁共振图像的步骤通过调用下述公式(5)执行计算获得：

其中，表示每个时间点对应的重建后的脑功能磁共振图像；表示所述初始重建图像；S表示所述敏感度矩阵，S^H表示S的共轭矩阵；γ表示正则化参数，γ可在[γ_min,γ_max]之间取值，其中，γ_min＝10^-4σ_M，σ_M为S的最大特征值，γ_max＝σ_M；是单位算子； 表示所述中间图像。

基于上述方法，本发明还提供了一种脑功能磁共振成像系统，其包括：

数据提取模块，用于获取基于多通道并行且稀疏采样获得的每一个相控阵线圈对应的动态数据图像和第一参考图像；

参考图像提取模块，用于根据所述第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像；

基函数计算模块，用于根据所述基准图像，获得高阶广义级模型的基函数；

级系数计算模块，用于利用L₂范数对高阶广义级模型的级系数进行约束，基于所述动态数据图像和所述基函数估计获得所述级系数；

第一图像重建模块，用于根据所述级系数的估计值和所述基准图像，基于高阶广义级模型重建出每个通道各个时间点对应的中间图像；及

第二图像重建模块，用于基于相控阵线圈的敏感度矩阵，对所述中间图像进行并行重建，获得重建后的动态脑功能磁共振图像。

在其中一个实施例中，所述基函数计算模块包括：

设定单元，用于给定控制参数μ；和

基函数计算单元，用于调用下述公式(1)计算获得所述高阶广义级模型的基函数：

其中，表示高阶广义级模型的基函数；表示第l个相控阵线圈在第q个时间点对应的基准图像；l＝1,...,L，L表示相控阵线圈的个数，每一个相控阵线圈采集的数据对应一个输出通道；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

在其中一个实施例中，所述第一图像重建模块包括：

参数提取单元，用于获取控制参数μ；和

中间图像获取单元，用于调用下述公式(3)计算每个通道各个时间点对应的中间图像：

${\widehat{\mathrm{\ρ}}}_{q,l}\left(x\right)=\underset{n=-M/2}{\overset{M/2-1}{\mathrm{\Σ}}}{c}_n({\mathrm{\ρ}}_{q,l}^{\mathrm{ref}}\left(x\right)+\mathrm{\μ})e^{i2\mathrm{\πn\Λkx}}$                    公式(3)

其中，c_n表示所述级系数；表示第l个相控阵线圈在第q个时间点对应的中间图像；M表示采集所述动态数据图像的相位编码数，n表示相位编码数变量，其取值范围为 表示第l个相控阵线圈在第q个时间点对应的基准图像；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

在其中一个实施例中，所述级系数计算模块包括：

特征向量提取单元，用于根据获得的基函数生成矩阵向量；

中间向量计算单元，用于依据欠采样因子构成的欠采样矩阵乘以所述矩阵向量获得中间矩阵向量；及

级系数提取单元，用于调用下述公式(4)计算所述级系数：

${\hat{c}}_n={(A^{H}A+\mathrm{\λI})}^{-1}{A}^{H}D$                    公式(4)

其中，表示所述级系数的估计值；A表示所述中间矩阵向量，A^H表示矩阵A的共轭矩阵；λ表示正则化参数；D表示所述动态数据图像，即针对第l个相控阵线圈的T个采样时间帧对应的K空间欠采样图像D_1,l(k),...,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数；是单位算子。

在其中一个实施例中，所述第二图像重建模块包括：

敏感度矩阵提取单元，用于根据针对每一个相控阵线圈采集的一个或多个全采样的参考图像，获得敏感度矩阵；

初始重建单元，用于结合所述敏感度矩阵和所述中间图像，采用最小二乘法估计每个时间点对应的初始重建图像；

再次重建单元，用于基于高阶广义级模型，利用L₂范数对所述初始重建图像进行约束，估计获得每个时间点对应的重建后的脑功能磁共振图像；及

输出单元，用于汇总各个时间点对应的重建后的脑功能磁共振图像，形成所述重建后的动态脑功能磁共振图像。

本发明基于传统的梯度回波序列，采用一种联合并行成像的高阶广义级模型(Higher-order Generalized Series，HGS)来实现高时空分辨的脑功能磁共振成像。相比于EPI成像，该方法能够在进一步提高时空分辨率的前提下，减少了图像的几何变形和Nyquist伪影。

附图说明

图1为本发明图像采集的示意图；

图2为本发明方法中并行重建的图像示意图；

图3为本发明的脑功能磁共振成像方法的流程示意图；

图4为本发明方法的一优选实施例的流程图；

图5为本发明方法的另一优选实施例的流程图；

图6为本发明的脑功能磁共振成像系统的结构示意图；

图7为本发明系统的一优选实施例的结构示意图。

具体实施方式

基于动态脑功能磁共振成像技术，本发明提供一种脑功能磁共振成像方法和系统，其基于传统的梯度回波序列，联合并行采集，采用高阶的广义级模型进行重建，该方法可以避免图像的几何畸变和伪影的前提下实现高分辨的fMRI成像。GS模型(Generalized Series:GS)早先已被应用到动态的MR成像，例如弥散加权成像，心脏成像，只需采集少量的相位编码线，远远的小于Nyquist采样就可以重建出高时间分辩的图像，但是这个模型不能捕捉到动态图像的细节特征，因此本发明结合并行成像的HGS模型，利用高频项来重建出局部的细节信息来提高空间分辨率，基于本发明的方法和系统只需稀疏采集k空间中间某些相位线，利用中间的数据进行HGS重建，得到一系列卷着的混叠图像，然后进一步结合线圈灵敏度矩阵进行并行重建来得到最终的高时空分辨的脑功能图像，本发明的方法和系统由于采用的是梯度回波序列采集，避免的EPI采集所产生的几何形变和Nyquist伪影。以下将结合附图详细说明本发明的各个具体实施例。

如图1、图2和图3所示，本实施例提供了一种脑功能磁共振成像方法，见图3，该方法包括以下步骤100至步骤600。

步骤100，获取基于多通道并行且稀疏采样获得的每一个相控阵线圈对应的动态数据图像和第一参考图像。这里的第一参考图像优选包括：针对每一个相控阵线圈在采集动态数据图像前后分别采集的至少两个K空间欠采样图像。而动态数据图像是指针对每一个相控阵线圈采集的T个采样时间帧对应的K空间欠采样图像D_1,l(k),...,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数。第一参考图像和动态数据图像均为K空间欠采样图像，但是动态数据图像是在k空间欠采集、且边缘未采集的图像。

本实施例的步骤100采集的图像数据优选采用如图1所示的方式，基于高阶广义级模型(HGS)，对于第l个线圈(l＝1,...,L)，获取一系列欠采集的动态数据图像以及第一参考图像，如图1所示，优选采集两个第一参考图像，即D_head,l(k)在获取动态数据图像集前获得，D_tail,l(k)在获取动态数据图像集后获得，HGS的第一参考图像是以欠采样因子R采集获得。如图1所示，多通道并行采集图像数据的过程依据M条相位编码线(图中实线)采集k空间中间某些行，形成K空间欠采样图像。欠采样因子R是已知给定的，比如一副全采集的图像，R＝1，如果R＝2，那么就只采集一副图像的1/2数据。那么采用过程中通过调节欠采样因子R和M条相位编码线可以获得不同的类型的K空间欠采样图像。

步骤200，根据第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像，其可以根据一个简单的指数模型得到。如图2所示，基于步骤100采集的第一参考图像，可以获得图2中a所表示的HGS参考图像，即上述基准图像，其对应每个相控阵线圈在每个动态数据图像获取时刻的图像数据。

步骤300，根据基准图像，获得高阶广义级模型的基函数。

步骤400，利用L₂范数对高阶广义级模型的级系数进行约束，基于动态数据图像和上述基函数估计获得级系数。

步骤500，根据级系数的估计值和基准图像，基于高阶广义级模型重建出每个通道各个时间点对应的中间图像。

步骤600，基于相控阵线圈的敏感度矩阵对中间图像进行并行重建，获得重建后的动态脑功能磁共振图像。上述过程主要是利用灵敏度编码技术(SENSE技术)对中间图像进行并行重建。

上述方法中，由于传统的GS模型采用k空间的低频信息进行重建，阶数较低，由于高阶的GS模型对fMRI成像来说，能够捕捉到一些细微的功能信号，所以高模阶的GS(即HGS模型)对脑功能成像显得很重要，但是由于高的模阶会产生较大的高频伪影，所以，本实施例中利用L₂范数约束高模阶的GS级系数c_n，可以去除掉高频伪影，同时改进空间分辨率、而且还可以信号的完整性。

基于上述实施例，本实施例中步骤400利用L₂范数对高阶广义级模型的级系数进行约束，估计获得级系数。这里在对高阶广义级模型的级系数进行约束时优选采用下述公式(2)的方式。

${\hat{c}}_n=\arg \underset{c_n}{\min }{\left|\right|D-\mathrm{\Ω\Ψ}{c}_n\left|\right|}_2^2+\mathrm{\λ}{\left|\right|c_n\left|\right|}_2^2$                    公式(2)

其中，表示级系数的估计值；是数据约束项；为二范数约束项；D表示动态数据图像(如图2中b所表示的示意图像)，即针对第l个相控阵线圈的T个采样时间帧对应的K空间欠采样图像D_1,l(k),...,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数；Ω表示欠采样因子构成的欠采样矩阵；Ψ表示上述基函数构成的矩阵向量；c_n表示高阶广义级模型的级系数；λ表示正则化参数。

上述步骤200至步骤500主要是基于高阶广义级模型(HGS)对图像进行初步重建，获得混叠伪影图像，如图2所示的d表示的示意图像。基于上述实施例，如图4所示，本实施例给出了上述步骤200至步骤500的最优实现方式，具体如下所示。

首先，执行步骤200，根据第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像该基准图像可以根据一个简单的指数模型得到。

然后，优选下面两个步骤执行上述步骤300。

执行步骤301，给定控制参数μ；

执行步骤302，调用下述公式(1)计算获得高阶广义级模型的基函数

其中，表示高阶广义级模型的基函数；表示第l个相控阵线圈在第q个时间点对应的基准图像；l＝1,...,L，L表示相控阵线圈的个数，每一个相控阵线圈采集的数据对应一个输出通道；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

其次，执行上述步骤400。当步骤400中优选上述公式(2)进行求解时，则最优方式是利用下述公式(4)进行计算，所以，本实施例中步骤400利用L₂范数对高阶广义级模型的级系数进行约束估计获得级系数的过程优选包括以下几个步骤：

步骤401，根据获得的基函数生成矩阵向量；

步骤402，依据欠采样因子构成的欠采样矩阵乘以矩阵向量获得中间矩阵向量；

步骤403，调用下述公式(4)计算级系数：

${\hat{c}}_n={(A^{H}A+\mathrm{\λI})}^{-1}{A}^{H}D$                    公式(4)

其中，表示级系数的估计值；A表示中间矩阵向量，A^H表示矩阵A的共轭矩阵，A＝ΩΨ，其中Ω表示欠采样因子构成的欠采样矩阵，Ψ表示根据上述基函数构成的矩阵向量；λ表示正则化参数；D表示动态数据图像，即针对第l个相控阵线圈的T个采样时间帧对应的K空间欠采样图像D_1,l(k),...,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数；是单位算子。

最后，采用以下两个步骤上述执行步骤500。

执行步骤501，获取控制参数μ；

执行步骤502，根据级系数的估计值和基准图像，调用下述公式(3)计算获得重建出的每个通道各个时间点对应的中间图像

${\widehat{\mathrm{\ρ}}}_{q,l}\left(x\right)=\underset{n=-M/2}{\overset{M/2-1}{\mathrm{\Σ}}}{c}_n({\mathrm{\ρ}}_{q,l}^{\mathrm{ref}}\left(x\right)+\mathrm{\μ})e^{i2\mathrm{\πn\Λkx}}$                    公式(3)

其中，c_n表示级系数；表示第l个相控阵线圈在第q个时间点对应的中间图像；M表示采集动态数据图像的相位编码数，n表示相位编码数变量，其取值范围为 表示第l个相控阵线圈在第q个时间点对应的基准图像；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

上述给出的步骤300、步骤400和步骤500的具体实现方式可以单独提出来用于另一个实施例，也可以自由组合，并且上述步骤400中仅给出了一种最优求解上述公式(2)的方法，但本发明不限于此，若采用其他方式执行步骤400获得级系数，也可以将其替换掉上述步骤401至步骤403的过程，而形成新的技术方案。

基于上述各个实施例，如图5所示，本实施例中步骤600利用灵敏度编码技术对中间图像进行并行重建的过程包括以下步骤：

步骤601，根据针对每一个相控阵线圈采集的一个或多个全采样的参考图像(参见图2中的c所表示的示意图像)，获得敏感度矩阵；

步骤602，结合敏感度矩阵和中间图像，采用最小二乘法估计每个时间点对应的初始重建图像。具体可参见下述公式(6)所示。

${\widehat{\mathrm{\ρ}}}_q\left(x\right)=\frac{{\mathrm{\Σ}}_{l=1}^L{\widehat{\mathrm{\ρ}}}_{q,l}\left(x\right)s_l^{*}\left(x\right)}{{\mathrm{\Σ}}_{l=1}^L{s}_l\left(x\right)s_l^{*}\left(x\right)}$                    公式(6)

其中，表示上述初始重建图像；l表示相控阵线圈序号，l＝1,...,L，L表示相控阵线圈的个数；表示上述重建出的每个通道各个时间点对应的中间图像；s_l(x)表示第l个相控阵线圈的敏感度矩阵，表示s_l(x)的伴随矩阵。

步骤603，基于高阶广义级模型，利用L₂范数对初始重建图像进行约束，估计获得每个时间点对应的重建后的脑功能磁共振图像(参见图2中的e所表示的示意图像)。这里利用L₂范数对初始重建图像进行约束时优选采用下述公式(7)所示模型。

其中，γ表示正则化参数，γ可在[γ_min,γ_max]之间取值，其中，γ_min＝10^-4σ_M，σ_M为S的最大特征值，γ_max＝σ_M；S表示敏感度矩阵；|| ||₂表示取L₂范数；为上述步骤500重建出的每个通道各个时间点对应的中间图像；表示上述初始重建图像 表示每个时间点对应的重建后的脑功能磁共振图像。

步骤604，汇总各个时间点对应的重建后的脑功能磁共振图像，形成重建后的动态脑功能磁共振图像。

基于上述实施例，上述步骤603可以考虑采用下述公式公式(5)来简单求解上述公式(7)。所以，本实施例中步骤603利用L₂范数对初始重建图像进行约束、估计获得每个时间点对应的重建后的脑功能磁共振图像的步骤优选通过调用下述公式(5)执行计算获得：

其中，表示每个时间点对应的重建后的脑功能磁共振图像；表示初始重建图像；S表示敏感度矩阵，S^H表示S的共轭矩阵；γ表示正则化参数，γ可在[γ_min,γ_max]之间取值，其中，γ_min＝10^-4σ_M，σ_M为S的最大特征值，γ_max＝σ_M；是单位算子； 表示上述中间图像。

以上各个实施例在具体说明中仅只针对相应步骤的实现方式进行了阐述，然后在逻辑不相矛盾的情况下，上述各个实施例是可以相互组合的而形成新的技术方案的，而该新的技术方案依然在本具体实施方式的公开范围内。同时，本发明也不限制上述各个步骤之间的先后顺序。上述方法中仅只公开了一种求解方式，本发明也不限于只保护上述最优实施例。

基于上述方法，本实施例还提供了一种脑功能磁共振成像系统700，如图6所示，其包括：

数据提取模块710，用于获取基于多通道并行且稀疏采样获得的每一个相控阵线圈对应的动态数据图像和第一参考图像；

参考图像提取模块720，用于根据第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像；

基函数计算模块730，用于根据基准图像，获得高阶广义级模型的基函数；

级系数计算模块740，用于利用L₂范数对高阶广义级模型的级系数进行约束，基于上述动态数据图像和基函数估计获得级系数；

第一图像重建模块750，用于根据级系数的估计值和基准图像，基于高阶广义级模型重建出每个通道各个时间点对应的中间图像；及

第二图像重建模块760，用于基于相控阵线圈的敏感度矩阵，对中间图像进行并行重建，获得重建后的动态脑功能磁共振图像。

上述各个功能模块的具体实现方式可以参见上述方法中关于步骤100至步骤600的相关说明。

基于上述实施例，如图7所示，本实施例中基函数计算模块730优选包括：

设定单元731，用于给定控制参数μ；和

基函数计算单元732，用于调用下述公式(1)计算获得高阶广义级模型的基函数：

其中，表示高阶广义级模型的基函数；表示第l个相控阵线圈在第q个时间点对应的基准图像；l＝1,...,L，L表示相控阵线圈的个数，每一个相控阵线圈采集的数据对应一个输出通道；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

基于上述实施例，如图7所示，本实施例中第一图像重建模块750优选包括：

参数提取单元751，用于获取控制参数μ；和

中间图像获取单元752，用于调用下述公式(3)计算每个通道各个时间点对应的中间图像：

${\widehat{\mathrm{\ρ}}}_{q,l}\left(x\right)=\underset{n=-M/2}{\overset{M/2-1}{\mathrm{\Σ}}}{c}_n({\mathrm{\ρ}}_{q,l}^{\mathrm{ref}}\left(x\right)+\mathrm{\μ})e^{i2\mathrm{\πn\Λkx}}$                    公式(3)

其中，c_n表示级系数；表示第l个相控阵线圈在第q个时间点对应的中间图像；M表示采集动态数据图像的相位编码数，n表示相位编码数变量，其取值范围为 表示第l个相控阵线圈在第q个时间点对应的基准图像；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

基于上述实施例，如图7所示，本实施例中级系数计算模块740优选包括：

特征向量提取单元741，用于根据获得的基函数生成矩阵向量；

中间向量计算单元742，用于依据欠采样因子构成的欠采样矩阵乘以上述矩阵向量获得中间矩阵向量；及

级系数提取单元743，用于调用下述公式(4)计算级系数：

${\hat{c}}_n={(A^{H}A+\mathrm{\λI})}^{-1}{A}^{H}D$                    公式(4)

其中，表示级系数的估计值；A表示中间矩阵向量，A^H表示矩阵A的共轭矩阵，A＝ΩΨ，其中Ω表示欠采样因子构成的欠采样矩阵，Ψ表示根据基函数构成的矩阵向量；λ表示正则化参数；D表示动态数据图像，即针对第l个相控阵线圈的T个采样时间帧对应的K空间欠采样图像D_1,l(k),...,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数；是单位算子。

基于上述实施例，如图7所示，本实施例中第二图像重建模块760优选包括：

敏感度矩阵提取单元761，用于根据针对每一个相控阵线圈采集的一个或多个全采样的参考图像，获得敏感度矩阵；

初始重建单元762，用于结合敏感度矩阵和中间图像，采用最小二乘法估计每个时间点对应的初始重建图像；

再次重建单元763，用于基于高阶广义级模型，利用L₂范数对初始重建图像进行约束，估计获得每个时间点对应的重建后的脑功能磁共振图像；及

输出单元764，用于汇总各个时间点对应的重建后的脑功能磁共振图像，形成重建后的动态脑功能磁共振图像。

上述系统中的各个功能模块或单元的具体实现方式可以参见上述方法步骤的相关说明，在此不累述。比如上述再次重建单元763可参见上述公式(7)和公式(5)的相关说明，再次重建单元763优选采用上述关于公式(5)的执行步骤。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个非易失性计算机可读存储介质(如ROM、磁碟、光盘)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的系统结构和方法。

综上所述，本发明提出了一种高时空分辨脑功能磁共振成像方法，该方法在前期已经成功的运用到弥散加权成像，心脏磁共振成像，完美的实现了高时空分辨的动态成像。本发明根据BOLD脑功能成像特点，在原有的HGS模型的基础上，进一步的做了改进来捕捉更细微的功能信号变化，其以降采样因子R并行采集获得至少两个HGS模型重建的参考图像。同时，采集一个或者是多个全采样的参考图像(采集欠采样的动态数据之前或者是之后)作为并行重建的参考图像，用于求解灵敏度矩阵。稀疏采集低频部分，高频部分未采得到随时间变换的一组动态数据。在传统的并行成像中，欠采样因子R经常要小于接受线圈L，但是在HGS模型中，R可以等于L，甚至可以大于L。然后利用HGS模型进行图像重建，得到一系列混叠伪影图像，其GS系数可以通过L2范数约束来获得，避免高阶的GS模型容易产生高频伪影。再次，对HGS重建获得的图像进一步结合并行重建的SENSE成像技术，获得高时空分辨的脑功能图像。

本发明的关键点和欲保护点主要有以下几方面：(1)本发明对现有的GS模型进行了改进，采用高阶的GS模型联合并行成像使其应用到脑功能成像中，可以获得高时空分辨的大脑激活图。(2)对HGS模型中的级系数c_n采用L₂范数约束重建，获得较为鲁棒的级系数。(3)充分结合参考图像及图像的灵敏度信息，更进一步地重建HGS混叠图像，获得高时空分辨的脑功能图像和刺激图。本发明基于高阶广义级模型，结合并行成像来实现脑功能成像，由于该模型是采用梯度回波序列，能够避免传统的EPI脑功能成像容易使得图像几何变形并且容易产生Nyquist伪影。因此，本发明的方法和系统能够实现高时空分辨的结构图像的前提下捕捉到大脑的细微激活信号，相比现有技术更加有效，获得的图像更加清晰、分辨率更高，并且减少了伪影，提高了图像质量。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (13)

1.一种脑功能磁共振成像方法，其包括：

获取基于多通道并行且稀疏采样获得的每一个相控阵线圈对应的动态数据图像和第一参考图像；

根据所述第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像；

根据所述基准图像，获得高阶广义级模型的基函数；

利用L₂范数对高阶广义级模型的级系数进行约束，基于所述动态数据图像和所述基函数估计获得所述级系数；

根据所述级系数的估计值和所述基准图像，基于高阶广义级模型重建出每个通道各个时间点对应的中间图像；

基于相控阵线圈的敏感度矩阵，对所述中间图像进行并行重建，获得重建后的动态脑功能磁共振图像。

2.根据权利要求1所述的脑功能磁共振成像方法，其特征在于，所述第一参考图像包括：针对每一个相控阵线圈在采集所述动态数据图像前后分别采集的至少两个K空间欠采样图像。

3.根据权利要求1所述的脑功能磁共振成像方法，其特征在于，所述动态数据图像为K空间欠采样、且边缘未采集的图像。

4.根据权利要求1所述的脑功能磁共振成像方法，其特征在于，所述根据所述基准图像获得高阶广义级模型的基函数的过程包括：

给定控制参数μ；

调用下述公式(1)计算获得所述高阶广义级模型的基函数：

   公式(1)

其中，表示高阶广义级模型的基函数；表示第l个相控阵线圈在第q个时间点对应的基准图像；l＝1,...,L，L表示相控阵线圈的个数，每一个相控阵线圈采集的数据对应一个输出通道；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

5.根据权利要求1所述的脑功能磁共振成像方法，其特征在于，所述根据所述级系数、基于高阶广义级模型重建出每个通道各个时间点对应的中间图像的过程包括：

获取控制参数μ；

调用下述公式(3)计算每个通道各个时间点对应的中间图像：

${\widehat{\mathrm{\ρ}}}_{q,l}\left(x\right)=\underset{n=-M/2}{\overset{M/2-1}{\mathrm{\Σ}}}{c}_n({\mathrm{\ρ}}_{q,l}^{\mathrm{ref}}\left(x\right)+\mathrm{\μ})e^{i2\mathrm{\πn\Λkx}}$    公式(3)

其中，c_n表示所述级系数；表示第l个相控阵线圈在第q个时间点对应的中间图像；M表示采集所述动态数据图像的相位编码数，n表示相位编码数变量，其取值范围为表示第l个相控阵线圈在第q个时间点对应的基准图像；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

6.根据权利要求1所述的脑功能磁共振成像方法，其特征在于，所述利用L₂范数对高阶广义级模型的级系数进行约束、基于所述动态数据图像和所述基函数估计获得所述级系数的过程包括：

根据获得的基函数生成矩阵向量；

依据欠采样因子构成的欠采样矩阵乘以所述矩阵向量获得中间矩阵向量；

调用下述公式(4)计算所述级系数：

${\hat{c}}_n={(A^{H}A+\mathrm{\λI})}^{-1}{A}^{H}D$    公式(4)

其中，表示所述级系数的估计值；A表示所述中间矩阵向量，A^H表示矩阵A的共轭矩阵；λ表示正则化参数；D表示所述动态数据图像，即针对第l个相控阵线圈的T个采样时间帧对应的K空间欠采样图像D_1,l(k),…,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数；是单位算子。

7.根据权利要求1所述的脑功能磁共振成像方法，其特征在于，所述利用灵敏度编码技术对所述中间图像进行并行重建的过程包括：

根据针对每一个相控阵线圈采集的一个或多个全采样的参考图像，获得敏感度矩阵；

结合所述敏感度矩阵和所述中间图像，采用最小二乘法估计每个时间点对应的初始重建图像；

基于高阶广义级模型，利用L₂范数对所述初始重建图像进行约束，估计获得每个时间点对应的重建后的脑功能磁共振图像；

汇总各个时间点对应的重建后的脑功能磁共振图像，形成所述重建后的动态脑功能磁共振图像。

8.根据权利要求1所述的脑功能磁共振成像方法，其特征在于，所述利用L₂范数对所述初始重建图像进行约束、估计获得每个时间点对应的重建后的脑功能磁共振图像的步骤通过调用下述公式(5)执行计算获得：

${\stackrel{\→}{\mathrm{\ρ}}}_r={\stackrel{\→}{\mathrm{\ρ}}}_q+{(S^{H}S+{\mathrm{\γ}}^{2}I)}^{-1}{S}^H\mathrm{\Δ}\stackrel{\→}{d}$    公式(5)

其中，表示每个时间点对应的重建后的脑功能磁共振图像；表示所述初始重建图像；S表示所述敏感度矩阵，S^H表示S的共轭矩阵；γ表示正则化参数，γ可在[γ_min,γ_max]之间取值，其中，γ_min＝10^-4σ_M，σ_M为S的最大特征值，γ_max＝σ_M；是单位算子； $\mathrm{\Δ}\stackrel{\→}{d}=\stackrel{\→}{d}-S{\stackrel{\→}{\mathrm{\ρ}}}_r,\stackrel{\→}{d}$ 表示所述中间图像。

9.一种脑功能磁共振成像系统，其特征在于，所述系统包括：

数据提取模块，用于获取基于多通道并行且稀疏采样获得的每一个相控阵线圈对应的动态数据图像和第一参考图像；

参考图像提取模块，用于根据所述第一参考图像，获得每个相控阵线圈在每个动态数据图像获取时刻处的基准图像；

基函数计算模块，用于根据所述基准图像，获得高阶广义级模型的基函数；

级系数计算模块，用于利用L₂范数对高阶广义级模型的级系数进行约束，基于所述动态数据图像和所述基函数估计获得所述级系数；

第一图像重建模块，用于根据所述级系数的估计值和所述基准图像，基于高阶广义级模型重建出每个通道各个时间点对应的中间图像；及

第二图像重建模块，用于基于相控阵线圈的敏感度矩阵，对所述中间图像进行并行重建，获得重建后的动态脑功能磁共振图像。

10.根据权利要求9所述的脑功能磁共振成像系统，其特征在于，所述基函数计算模块包括：

设定单元，用于给定控制参数μ；和

基函数计算单元，用于调用下述公式(1)计算获得所述高阶广义级模型的基函数：

   公式(1)

其中，表示高阶广义级模型的基函数；表示第l个相控阵线圈在第q个时间点对应的基准图像；l＝1,...,L，L表示相控阵线圈的个数，每一个相控阵线圈采集的数据对应一个输出通道；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

11.根据权利要求9所述的脑功能磁共振成像系统，其特征在于，所述第一图像重建模块包括：

参数提取单元，用于获取控制参数μ；和

中间图像获取单元，用于调用下述公式(3)计算每个通道各个时间点对应的中间图像：

${\widehat{\mathrm{\ρ}}}_{q,l}\left(x\right)=\underset{n=-M/2}{\overset{M/2-1}{\mathrm{\Σ}}}{c}_n({\mathrm{\ρ}}_{q,l}^{\mathrm{ref}}\left(x\right)+\mathrm{\μ})e^{i2\mathrm{\πn\Λkx}}$    公式(3)

其中，c_n表示所述级系数；表示第l个相控阵线圈在第q个时间点对应的中间图像；M表示采集所述动态数据图像的相位编码数，n表示相位编码数变量，其取值范围为表示第l个相控阵线圈在第q个时间点对应的基准图像；μ是个控制参数，Λk满足奈奎斯采样率的采样间隔。

12.根据权利要求9所述的脑功能磁共振成像系统，其特征在于，所述级系数计算模块包括：

特征向量提取单元，用于根据获得的基函数生成矩阵向量；

中间向量计算单元，用于依据欠采样因子构成的欠采样矩阵乘以所述矩阵向量获得中间矩阵向量；及

级系数提取单元，用于调用下述公式(4)计算所述级系数：

${\hat{c}}_n={(A^{H}A+\mathrm{\λI})}^{-1}{A}^{H}D$    公式(4)

其中，表示所述级系数的估计值；A表示所述中间矩阵向量，A^H表示矩阵A的共轭矩阵；λ表示正则化参数；D表示所述动态数据图像，即针对第l个相控阵线圈的T个采样时间帧对应的K空间欠采样图像D_1,l(k),…,D_T,l(k)，l＝1,...,L，L表示相控阵线圈的个数；是单位算子。

13.根据权利要求9所述的脑功能磁共振成像系统，其特征在于，所述第二图像重建模块包括：

敏感度矩阵提取单元，用于根据针对每一个相控阵线圈采集的一个或多个全采样的参考图像，获得敏感度矩阵；

初始重建单元，用于结合所述敏感度矩阵和所述中间图像，采用最小二乘法估计每个时间点对应的初始重建图像；

再次重建单元，用于基于高阶广义级模型，利用L₂范数对所述初始重建图像进行约束，估计获得每个时间点对应的重建后的脑功能磁共振图像；及

输出单元，用于汇总各个时间点对应的重建后的脑功能磁共振图像，形成所述重建后的动态脑功能磁共振图像。