CN103064046A - 一种基于稀疏采样的核磁共振成像的图像处理方法 - Google Patents

Abstract

本发明提供了一种基于稀疏采样的核磁共振成像的图像处理方法，其利用多个相控线圈同时以高密度采集K空间中心区域的数据、以高斯分布对K空间中心区域四周随机降采样K空间数据；将每一个线圈所采集的Ｋ空间数据进行傅里叶变换转换到图像空间，根据每个线圈的敏感度信息对所述图像空间的信号进行线性拟合，形成重建的自旋密度图像ρ，并利用PMRI降采样后的频域信号获得参考图像，从而基于局部噪声方差信息对正则约束项的参数进行定量计算来实现准确的图像重建。

Description

一种基于稀疏采样的核磁共振成像的图像处理方法

【技术领域】

本发明涉及基于稀疏采样的核磁共振成像的图像处理方法。

【背景技术】

核磁共振成像(MRI)利用人体组织中氢原子核（质子）在磁场中受到射频脉冲的激励而发生核磁共振现象，产生磁共振信号，经过电子计算机处理，重建出人体某一层面的图像的成像技术。但是，该成像技术基于理论和硬件等条件，使得它的成像速度相对慢于其它成像方法，并且对于动态的对象重建的图像较为模糊。于是，为了加速磁共振成像采样速度，许多研究者提出基于稀疏采样理论来实现快速磁共振成像。然而，基于稀疏采样理论的快速磁共振成像方法都是采用非现象约束重建方法来实现，无法满足不同对象、不同图像获取方式下对图像的后续处理，使得快速获得高分辨率的图像。

【发明内容】

鉴于上述状况，本发明提供了一种基于稀疏采样的核磁共振成像的图像处理方法，其具体技术方案如下所示：

一种基于稀疏采样的核磁共振成像的图像处理方法，其包括，

第一步，利用多个相控线圈同时以高密度采集K空间中心区域的数据、以高斯分布对K空间中心区域四周随机降采样K空间数据；将每一个线圈所采集的Ｋ空间数据进行傅里叶变换转换到图像空间，根据每个线圈的敏感度信息对所述图像空间的信号进行线性拟合，形成重建的自旋密度图像ρ，具体采用以下公式1：

$d_l\left(k\right)=\∫\mathrm{\ρ}\left(\stackrel{\→}{r}\right)s_l\left(\stackrel{\→}{r}\right)e^{-i\stackrel{\→}{k}\·\stackrel{\→}{r}}d\left(\stackrel{\→}{r}\right)$ （公式1）

l＝1,2,…N

其中，ρ是自旋密度图像；s_l是第l个线圈灵敏度矩阵；d_l是第l个线圈所获得的欠采数据矩阵；

是图像空间中的位置坐标向量（x,y,z）；k是K空间中的位置坐标向量（kx,ky,kz）；N为利用并行磁共振成像（PMRI）采集预扫描数据的线圈总个数，l是当前采集数据的线圈序数；

其中

是通过N个通道加权求均值后的K空间数据;

第二步，结合压缩感知成像（CS）稀疏约束来重建自旋密度图像ρ，获得重构图像ρ′，其重建方法包括：

按照下述公式2利用非线性迭代方法求解重构图像ρ′，

$\left\{\begin{array}{c}{\mathrm{\&rho;}}^{\&prime;}=\arg \min \mathrm{\&lambda;}{\left|\right|\mathrm{W\&rho;}\left|\right|}_1+{\mathrm{\&lambda;}}_2\mathrm{TV}\left(\mathrm{\&rho;}\right)\\ {\left|\right|\left(\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_n\end{array}\right)-\left(\begin{array}{c}\mathrm{PF\&rho;}\&CenterDot;s_1\\ \mathrm{PF\&rho;}\&CenterDot;s_2\\ .\\ .\\ .\\ {\mathrm{PF\&rho;}\&CenterDot;s}_n\end{array}\right)\left|\right|}_2<\mathrm{\&epsiv;}\end{array}\right.$ （公式2）

上述公式2中，F是傅里叶变换；P是欠采样因子；W为稀疏变换操作符；ε为噪声因子；TV(ρ)为约束条件；λ为局部正则化参数、λ₂为全局正则项参数；ρ是自旋密度图像；d_l是第l个线圈所获得的欠采数据矩阵，l＝1,2,…N；s_l是第l个线圈灵敏度矩阵，l＝1,2,…N；

上述约束条件TV(X)函数采用下述公式3：

$\mathrm{TV}\left(X\right)=\underset{i,j}{\mathrm{\&Sigma;}}\sqrt{{(X_{i,j}-X_{i-1,j})}^2+{(X_{i,j}-X_{i,j-1})}^2}---\left(3\right)$

其中，X即表示上述求解得到的自旋密度图像ρ向量，X_i,j表示ρ向量中i行、j列的元素值；

在上述公式2的迭代过程中，以第一步获得的重建图像作为引入了噪声的降质图像，并根据降质图像的噪声分布情况选取当前迭代次数下对应的局部正则化参数λ；

第三步，迭代到达设定值，输出重构图像。

其中，所述全局正则项参数λ₂的获取采用以下公式8：

${\mathrm{\&lambda;}}_2=\frac{{\left|\right|\left(\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_n\end{array}\right)-\left(\begin{array}{c}\mathrm{PF\&rho;}\&CenterDot;s_1\\ \mathrm{PF\&rho;}\&CenterDot;s_2\\ .\\ .\\ .\\ {\mathrm{PF\&rho;}\&CenterDot;s}_n\end{array}\right)\left|\right|}_2^2}{{\left|\right|\mathrm{W\&rho;}\left|\right|}_1+\mathrm{\&epsiv;}}$ （公式8）

上式中，F是傅里叶变换；P是欠采样因子；W为稀疏变换操作符；ε为噪声因子；ρ是自旋密度图像，即上述第一步获得的降质图像；d_l是第l个线圈所获得的欠采数据矩阵，l＝1,2,…N；s_l是第l个线圈灵敏度矩阵，l＝1,2,…N。其中，根据降质图像的噪声分布情况选取当前迭代次数下对应的局部正则化参数λ的方法如下所示：

第一步，根据降质图像g经均匀算子d卷积后方差的变化量，估计获得降质图像g的噪声方差矩阵；

第二步，根据下述公式获取局部正则化参数：

λ(i,j)＝α(i,j)·λ₂                            (公式7)

其中，α(i,j)称为正则化权系数；λ(i,j)为图像中任一像素(i,j)处局部正则化参数，λ₂表示全局正则项参数;

在估计出降质图像的噪声方差矩阵后，对于具有最大估计噪声

的局部区域，取一个较大值V_max作为该区域正则化权系数α(i,j)，而对于具有最小估计噪声

的局部区域，取一个较小的值V_min，作为该局部区域的正则化权系数α(i,j),基于上述获得的区域正则化权系数α(i,j)，建立一个由 $\left[\begin{array}{cc}{\mathrm{\&sigma;}}_{n\max }^2& {\mathrm{\&sigma;}}_{n\min }^2\end{array}\right]$ 至[V_maxV_min]的线性映射函数，对于其他区域处的正则化权系数α(i,j)通过该线性映射函数在[V_max V_min]之间选取。

其中，所述降质图像g的噪声方差矩阵的获取方法如下：

对降质图像g经均匀算子d卷积后获得第一图像z，计算公式如下所示：

z＝d·g                (公式4)

其中，d＝[0.5 0.5]或[0.25 0.25;0.25 0.25]；

对第一图像z进行平滑处理，利用下述公式6估计降质图像g的噪声方差

${\mathrm{\&sigma;}}_n^2=r{\mathrm{\&sigma;}}_g^2-{\mathrm{\&sigma;}}_z^2$ （公式6）

上式中，修正因子r选取范围为0<r<1；

基于上述噪声方差的估计结果利用窗口卷积得出所述降质图像g中每一个像素在图像窗口大小的噪声方差，从而获得所述降质图像g的噪声方差矩阵。

其中，所述线性拟合采用最小二乘法。

其中，所述迭代采用共轭梯度法。

其中，修正因子r选取0.5。

本发明提供了一种利用并行磁共振成像（Parallel MRI，简称PMRI）中的先验信息，基于局部噪声方差信息对正则约束项的参数进行定量计算来实现准确的图像重建的方法。因为对频域数据随机降采样后再经傅里叶反变换得到的信号与等间隔降采样的情况截然不同，随机降采样重建的信号中出现极不相似（incoherent）的伪影，这些伪影与加性噪声非常相似。所以，利用PMRI降采样后的频域信号获得参考图像，从而基于局部噪声方差信息对正则约束项的参数进行定量计算来实现准确的图像重建。

【附图说明】

图1为本发明方法的原理图；

图2（a）-图2（c）为测试结果对比图；

图3(A)-图3（B）为等间隔采样A与随机降采样B重建结果对比。

【具体实施方式】

为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的较佳的实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。

如图1所示，本发明实施例提供了一种基于稀疏采样的核磁共振成像的图像处理方法，其主要包括以下步骤：

第一步，利用多个相控线圈同时以高密度采集K空间中心区域的数据、以高斯分布对K空间中心区域四周随机降采样K空间数据，将每一个线圈所采集的Ｋ空间数据进行傅里叶变换转换到图像空间，根据每个线圈的敏感度信息对所述图像空间的信号进行线性拟合，形成重建的自旋密度图像ρ，具体采用以下公式1：

$d_l\left(k\right)=\&Integral;\mathrm{\&rho;}\left(\stackrel{\&RightArrow;}{r}\right)s_l\left(\stackrel{\&RightArrow;}{r}\right)e^{-i\stackrel{\&RightArrow;}{k}\&CenterDot;\stackrel{\&RightArrow;}{r}}d\left(\stackrel{\&RightArrow;}{r}\right)$ （公式1）

l＝1,2,…N

其中，ρ是自旋密度图像；s_l是第l个线圈灵敏度矩阵；d_l是第l个线圈所获得的欠采数据矩阵；

是图像空间中的位置坐标向量（x,y,z）；k是K空间中的位置坐标向量（kx,ky,kz）；N为利用并行磁共振成像（PMRI）采集预扫描数据的线圈总个数，l是当前采集数据的线圈序数；

其中

是通过N个通道加权求均值后的K空间数据。

这一步中，主要是基于并行磁共振成像（Parallel MRI，简称PMRI）原理，利用多个相控线圈同时欠采样获取K空间数据，对于未采集的数据进行估计，组合成全K空间数据后进行傅里叶变换重建出图像。具体又可以分为基于K空间（K空间是寻常空间Rⁿ在傅里叶变换下的对偶空间，主要应用在磁振造影的成像分析，其他如磁振造影中的射频波形设计，以及量子计算中的初始态准备亦用到K空间的概念。）处理的自校正并行获取(Generalized autocalibratingpartially parallel acquisitions，简称GRAPPA)算法以及基于图像域处理的敏感度编码(Sensitivity encoding，简称Sense)成像，本发明实施例在此优选采用Sense成像方法。基于MRI的多相控阵线圈敏感度编码(Sense)成像过程中，由每一个线圈采集的K空间数据经傅里叶变换可以重建出一幅含混迭伪影的图像，根据每个线圈的敏感度信息通过线性拟合方式（比如最小二乘法）消除混迭伪影重建图像，具体成像公式1所示。

对于多个相控线圈笛卡尔采样轨迹而言，由于能量主要集中在磁共振数据K空间的低频部分，所以本发明实施例的采样轨道中，需对K空间中心区域的低频部分进行高密集采样，而对中心区域四周的高频部分按高斯分布对并行磁共振成像（Parallel MRI，简称PMRI）K空间进行相同轨迹随机降采样。本步骤中，对于并行磁共振成像（PMRI）K空间随机欠采样时，由于磁共振成像中，K空间的中心部分包含了图像的低频信号，所以采用伪随机（Paseudo-random）的变密度采样形式，以高概率采集K空间中心数据，因此选择N个通道的加权求均值作为K空间数据。

依据上述公式1，在预先获知线圈灵敏度矩阵的情况下，将依据上述公式1获得的不同线圈的采集数据进行组合，利用最小二乘法线性拟合，求解公式1中的自旋密度图像ρ，即为利用并行磁共振成像（PMRI）采集预扫描对象欠采集的K空间数据。

第二步，结合压缩感知成像（CS）稀疏约束来重建自旋密度图像ρ，获得重构图像ρ′。这里的压缩感知成像（CS），又称压缩采样、压缩传感，它通过开发信号的稀疏特性，在远小于奈奎斯特（Nyquist）采样率的条件下，用随机采样获取信号的离散样本，然后通过非线性重建技术获得重建信号。这一步骤中，根据压缩感知成像（CS）被引入到并行磁共振成像（PMRI）之后，很多先验信息都可以作为先验约束在图像重建中。具体如下述重建方法：

按照下述公式2利用非线性迭代方法求解重构图像ρ′，

$\left\{\begin{array}{c}{\mathrm{\&rho;}}^{\&prime;}=\arg \min \mathrm{\&lambda;}{\left|\right|\mathrm{W\&rho;}\left|\right|}_1+{\mathrm{\&lambda;}}_2\mathrm{TV}\left(\mathrm{\&rho;}\right)\\ {\left|\right|\left(\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_n\end{array}\right)-\left(\begin{array}{c}\mathrm{PF\&rho;}\&CenterDot;s_1\\ \mathrm{PF\&rho;}\&CenterDot;s_2\\ .\\ .\\ .\\ {\mathrm{PF\&rho;}\&CenterDot;s}_n\end{array}\right)\left|\right|}_2<\mathrm{\&epsiv;}\end{array}\right.$ （公式2）

上述公式2中，F是傅里叶变换；P是欠采样因子；W为稀疏变换操作符，代表对于某一幅图像做小波变换；ε为噪声因子，一般设为10的负10次方左右；TV(ρ)为约束条件；λ为局部正则化参数、λ₂为全局正则项参数；ρ是自旋密度图像；d_l是第l个线圈所获得的欠采数据矩阵，l＝1,2,…N；s_l是第l个线圈灵敏度矩阵，l＝1,2,…N；

上述约束条件TV(X)函数采用下述公式3：

$\mathrm{TV}\left(X\right)=\underset{i,j}{\mathrm{\&Sigma;}}\sqrt{{(X_{i,j}-X_{i-1,j})}^2+{(X_{i,j}-X_{i,j-1})}^2}$ （公式3）

其中，X即表示上述求解得到的自旋密度图像ρ向量，X_i,j表示ρ向量中i行、j列的元素值；

在上述公式2的迭代过程中，以第一步获得的重建图像作为引入了噪声的降质图像，并根据降质图像的噪声分布情况选取当前迭代次数下对应的局部正则化参数λ。

第三步，迭代到达设定值，输出重构图像。

比如上述公式（2）中采用TV(X)作为约束条件增强MR图像在特定的变换域和有限差分变换域的稀疏性，并在运算过程加入噪声序列设置限制条件共同重构图像。求解上述公式2可以有很多种非线性迭代方法，比如最速下降法，共轭梯度法，高斯牛顿法等。在一个迭代搜索方法中，最关键的是在每一步迭代中找到搜索方向，共轭梯度法是一种线性搜索方法，是目前求解线性方程中最重要的一个方法，它可以很快的线性收敛到局部最优解。因此在求解非线性最优化公式2时，本发明实施例优选采用共轭梯度法来求解。

而对于上述提到的根据降质图像的噪声分布情况自适应选取局部正则化参数λ是为了加强图像在变换域的稀疏表示，进一步提供图像重构的精度。

经过上述处理后，因为对频域数据随机降采样后再经傅里叶反变换得到的信号与等间隔降采样的情况截然不同，随机降采样重建的信号中出现极不相似（incoherent）的伪影，这些伪影与加性噪声非常相似（见附图1）。所以，上述根据第一步中随机降采样重建的自旋密度图像ρ可以看作加噪后的降质图像，具体自适应选取局部正则化参数的方法如下所示：

第一步，根据降质图像g经均匀算子d卷积后方差的变化量，估计获得降质图像g的噪声方差矩阵。具体步骤如下所示：

首先，对降质图像g经均匀算子d卷积后获得第一图像z，计算公式如下所示：

z＝d·g                        (公式4)

为了保证降质图像g中模糊部分在图像平均前后方差近似相等，在选取均匀算子d时，其尺寸应尽量小，本发明实施例中选取d＝[0.5 0.5]或[0.25 0.25;0.250.25]。这样小尺寸的均匀算子d的取法，可使降质图像g中模糊部分的每一个像素值在均匀算子d的尺寸范围内的像素值都比较接近。

然后，对第一图像z进行平滑处理，平滑处理的具体处理方法如下所示：

利用降质图像g经均匀算子d卷积后方差的变化量

将其认为是平滑了噪声的结果，也就是，降质图像g的方差

减去图像z的方差

即认为是降质图像g的噪声方差

如下式（5）所示： ${\mathrm{\&sigma;}}_n^2={\mathrm{\&sigma;}}_g^2-{\mathrm{\&sigma;}}_z^2$ （公式5）

考虑到实际运用中，降质图像经均匀算子d卷积后方差也会有所下降，噪声也不可能完全得到平滑，因此，在利用上述公式5估计降质图像g的噪声方差

作为平滑处理结果时，时引入修正因子r，以便提高平滑处理的效果，具体可参见下面的公式6估计降质图像g的噪声方差， ${\mathrm{\&sigma;}}_n^2=r{\mathrm{\&sigma;}}_g^2-{\mathrm{\&sigma;}}_z^2$ （公式6）

上式中，修正因子r为降质图像g的自相关系数，降质图像g模糊程度越高时，相关性越高，r越接近于1，选取范围为0<r<1，本发明实施例中r优选0.5。

利用上述有关降质图像g的噪声方差

的估计公式，选取合适大小的图像窗口（这里的图像窗口大小经验取得，一般选取类似3*3的窗口），基于上述噪声方差

的估计结果利用窗口卷积得出降质图像g中每一个像素在图像窗口大小的噪声方差，得到一个估计的降质图像g的噪声方差矩阵，即为降质图像g噪声分布的一个估计结果。

第二步，获取自适应局部正则化参数。在这一步中，本发明根据图像中的不同区域的信息，可自适应的选取正则化参数，相比原有基于稀疏采样的PMRI算法中正则化参数对于图像中每一个像素值都是同一个值，而本发明得到的正则化参数可以对应图像中不同像素选取会有不同的正则化参数值，从而有效的地保持了图像不同区域平滑噪声与保留细节之间的平衡。具体过程如下所示：

对于图像中任一像素(i,j)处局部正则化参数λ(i,j)的选取，包含下述公式中的两个部分，具体如下式（7）所示，

λ(i,j)＝α(i,j)·λ₂                        (公式7)

其中，α(i,j)称为正则化权系数，其值由图像局部噪声方差确定，λ₂表示全局正则项参数。

对于局部噪声方差大的地方，α(i,j)取较大的值以平滑噪声，而局部噪声方差小的地方，α(i,j)可取较小的值以利于边缘恢复。

在估计出降质图像的噪声方差矩阵后，对于具有最大估计噪声

的局部区域，取一个较大值V_max作为该区域正则化权系数α(i,j)，而对于具有最小估计噪声

的局部区域，取一个较小值V_min，作为该局部区域的正则化权系数α(i,j)。这里的较大较小值V_max、V_min均可按照经验值取得，而且大小是相对于V_max、V_min两者相比而言的，比如这里的V_max必须是大于V_min的值。

基于上述获得的区域正则化权系数α(i,j)，建立一个由 $\left[\begin{array}{cc}{\mathrm{\&sigma;}}_{n\max }^2& {\mathrm{\&sigma;}}_{n\min }^2\end{array}\right]$ 至[V_maxV_min]的线性映射函数，对于其他区域处的正则化权系数α(i,j)通过该线性映射函数在[V_max V_min]之间选取。

上述步骤中，基于线性映射函数自适应求取正则化权系数α(i,j)，可以确保求得的局部正则项参数大小只取决于图像局部噪声方差的相对大小，而不是图像局部噪声方差的绝对值。

正则化参数λ(i,j)的第二部分全局正则项参数λ₂，本发明实施例中采取Kang的全局正则化参数计算方法，即下述公式8所示的： ${\mathrm{\&lambda;}}_2=\frac{{\left|\right|\left(\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_n\end{array}\right)-\left(\begin{array}{c}\mathrm{PF\&rho;}\&CenterDot;s_1\\ \mathrm{PF\&rho;}\&CenterDot;s_2\\ .\\ .\\ .\\ {\mathrm{PF\&rho;}\&CenterDot;s}_n\end{array}\right)\left|\right|}_2^2}{{\left|\right|\mathrm{W\&rho;}\left|\right|}_1+\mathrm{\&epsiv;}}$ （公式8）

上式中，F是傅里叶变换；P是欠采样因子；W为稀疏变换操作符，表示对于图像或者矩阵做小波变换；ε为噪声因子，靠经验值，一般设为10的负10次方左右；ρ是自旋密度图像，即上述第一步获得的降质图像；d_l是第l个线圈所获得的欠采数据矩阵，l＝1,2,…N；s_l是第l个线圈灵敏度矩阵，l＝1,2,…N；||Wρ||₁表示Wρ向量的1范数，同理，|| ||₂表示向量的2范数。

全局正则项参数λ₂在每一次迭代时的这个图像空间上是不变的，在迭代恢复算法中它只随迭代每一步的ρ的改变而改变。

从上述方法可以看出，本发明实施例主要是采用PMRI中的先验信息，基于局部噪声方差信息对正则约束项的参数进行定量计算来实现准确的图像重建，这样在保证快速采集的同时不损失图像信噪比。本发明提供的利用PMRI中的先验信息，基于局部噪声方差信息对正则约束项的参数进行定量计算来实现准确的图像重建的方法，主要采用的是对K空间中间低频部分进行密集采样，而对两边高频部分按高斯(Gaussian)分布对PMRI K空间进行相同轨迹随机降采样，具体如图3A和B所示，本发明设计的采样轨道中，需因为对频域数据随机降采样后再经傅里叶反变换得到的信号与等间隔降采样的情况截然不同，随机降采样重建的信号中出现极不相似（incoherent）的伪影，这些伪影与加性噪声非常相似。所以利用PMRI降采样后的频域信号获得参考图像，从而基于局部噪声方差信息对正则约束项的参数进行定量计算来实现准确的图像重建，保证快速采集的同时不损失图像信噪比，从而有效的地保持了图像不同区域平滑噪声与保留细节之间的平衡，使得图像分辨率和清晰度都相对提高，具体如图2所示，图2中，（a）为模板参考图像；(b）给出了本发明实施例的重建图像；（c）全局同一个正则化参数。通过对比同样可以看出，本发明实施例无论在图像质量及边缘的细节分辨都具有明显优势。

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

Claims (7)

1.一种基于稀疏采样的核磁共振成像的图像处理方法，其特征在于，包括：

第一步，利用多个相控线圈同时以高密度采集K空间中心区域的数据、以高斯分布对K空间中心区域四周随机降采样K空间数据；将每一个线圈所采集的Ｋ空间数据进行傅里叶变换转换到图像空间，根据每个线圈的敏感度信息对所述图像空间的信号进行线性拟合，形成重建的自旋密度图像ρ，具体采用以下公式1：

$d_l\left(k\right)=\&Integral;\mathrm{\&rho;}\left(\stackrel{\&RightArrow;}{r}\right)s_l\left(\stackrel{\&RightArrow;}{r}\right)e^{-i\stackrel{\&RightArrow;}{k}\&CenterDot;\stackrel{\&RightArrow;}{r}}d\left(\stackrel{\&RightArrow;}{r}\right)$ （公式1）

l=1,2,…N

其中，ρ是自旋密度图像；s_l是第l个线圈灵敏度矩阵；d_l是第l个线圈所获得的欠采数据矩阵；

是图像空间中的位置坐标向量（x,y,z）；k是K空间中的位置坐标向量（kx,ky,kz）；N为利用并行磁共振成像采集预扫描数据的线圈总个数，l是当前采集数据的线圈序数；

其中

是通过N个通道加权求均值后的K空间数据;

第二步，结合压缩感知成像稀疏约束来重建自旋密度图像ρ，获得重构图像ρ′，其重建方法包括：

按照下述公式2利用非线性迭代方法求解重构图像ρ′，

$\left\{\begin{array}{c}{\mathrm{\&rho;}}^{\&prime;}=\arg \min \mathrm{\&lambda;}{\left|\right|\mathrm{W\&rho;}\left|\right|}_1+{\mathrm{\&lambda;}}_2\mathrm{TV}\left(\mathrm{\&rho;}\right)\\ {\left|\right|\left(\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_n\end{array}\right)-\left(\begin{array}{c}\mathrm{PF\&rho;}\&CenterDot;s_1\\ \mathrm{PF\&rho;}\&CenterDot;s_2\\ .\\ .\\ .\\ {\mathrm{PF\&rho;}\&CenterDot;s}_n\end{array}\right)\left|\right|}_2<\mathrm{\&epsiv;}\end{array}\right.$ （公式2）

上述公式2中，F是傅里叶变换；P是欠采样因子；W为稀疏变换操作符；ε为噪声因子；TV(ρ)为约束条件；λ为局部正则化参数、λ₂为全局正则项参数；ρ是自旋密度图像；d_l是第l个线圈所获得的欠采数据矩阵，l=1,2,…N；s_l是第l个线圈灵敏度矩阵，l=1,2,…N；

上述约束条件TV(X)函数采用下述公式3： $\mathrm{TV}\left(X\right)=\underset{i,j}{\mathrm{\&Sigma;}}\sqrt{{(X_{i,j}-X_{i-1,j})}^2+{(X_{i,j}-X_{i,j-1})}^2}$ （公式3）

其中，X即表示上述求解得到的自旋密度图像ρ向量，X_i,j表示ρ向量中i行、j列的元素值；

在上述公式2的迭代过程中，以第一步获得的重建图像作为引入了噪声的降质图像，并根据降质图像的噪声分布情况选取当前迭代次数下对应的局部正则化参数λ；

第三步，迭代到达设定值，输出重构图像。

2.如权利要求1所述的基于稀疏采样的核磁共振成像的图像处理方法，其特征在于，所述全局正则项参数λ₂的获取采用以下公式8：

${\mathrm{\&lambda;}}_2=\frac{{\left|\right|\left(\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_n\end{array}\right)-\left(\begin{array}{c}\mathrm{PF\&rho;}\&CenterDot;s_1\\ \mathrm{PF\&rho;}\&CenterDot;s_2\\ .\\ .\\ .\\ {\mathrm{PF\&rho;}\&CenterDot;s}_n\end{array}\right)\left|\right|}_2^2}{{\left|\right|\mathrm{W\&rho;}\left|\right|}_1+\mathrm{\&epsiv;}}$ （公式8）

上式中，F是傅里叶变换；P是欠采样因子；W为稀疏变换操作符；ε为噪声因子；ρ是自旋密度图像，即上述第一步获得的降质图像；d_l是第l个线圈所获得的欠采数据矩阵，l=1,2,…N；s_l是第l个线圈灵敏度矩阵，l=1,2,…N。

3.如权利要求1或2所述的基于稀疏采样的核磁共振成像的图像处理方法，其特征在于，根据降质图像的噪声分布情况选取当前迭代次数下对应的局部正则化参数λ的方法如下所示：

第一步，根据降质图像g经均匀算子d卷积后方差的变化量，估计获得降质图像g的噪声方差矩阵；

第二步，根据下述公式获取局部正则化参数：

λ(i,j)=α(i,j)·λ₂                (公式7)

其中，α(i,j)称为正则化权系数；λ(i,j)为图像中任一像素(i,j)处局部正则化参数，λ₂表示全局正则项参数；

在估计出降质图像的噪声方差矩阵后，对于具有最大估计噪声

的局部区域，取一个较大值V_max作为该区域正则化权系数α(i,j)，而对于具有最小估计噪声

的局部区域，取一个较小值V_min，作为该局部区域的正则化权系数α(i,j),基于上述获得的区域正则化权系数α(i,j)，建立一个由 $\left[\begin{array}{cc}{\mathrm{\&sigma;}}_{n\max }^2& {\mathrm{\&sigma;}}_{n\min }^2\end{array}\right]$ 至[V_max V_min]的线性映射函数，对于其他区域处的正则化权系数α(i,j)通过该线性映射函数在[V_max V_min]之间选取。

4.如权利要求3所述的基于稀疏采样的核磁共振成像的图像处理方法，其特征在于，所述降质图像g的噪声方差矩阵的获取方法如下：

对降质图像g经均匀算子d卷积后获得第一图像z，计算公式如下所示：

z=d·g                        (公式4)

其中，d＝[0.5 0.5]或[0.25 0.25;0.25 0.25]；

对第一图像z进行平滑处理，利用下述公式6估计降质图像g的噪声方差

${\mathrm{\&sigma;}}_n^2=r{\mathrm{\&sigma;}}_g^2-{\mathrm{\&sigma;}}_z^2$ （公式6）

上式中，修正因子r选取范围为0<r<1；

基于上述噪声方差的估计结果利用窗口卷积得出所述降质图像g中每一个像素在图像窗口大小的噪声方差，从而获得所述降质图像g的噪声方差矩阵。

5.如权利要求1所述的基于稀疏采样的核磁共振成像的图像处理方法，其特征在于，所述线性拟合采用最小二乘法。

6.如权利要求1所述的基于稀疏采样的核磁共振成像的图像处理方法，其特征在于，所述迭代采用共轭梯度法。

7.如权利要求4所述的基于稀疏采样的核磁共振成像的图像处理方法，其特征在于，修正因子r选取0.5。

Images