CN103443643A - 用于执行并行磁共振成像的方法 - Google Patents

Abstract

一种身体的并行磁共振成像的方法，包括：从具有已知的或估算的灵敏度分布和噪声协方差矩阵的各个接收天线中获取所述身体的一组基础磁共振图像，所述基础图像在k-空间中欠采样；以及执行所述身体的磁共振图像的正则化重建，其中，所述执行磁共振图像的正则化重建的步骤是无人监管的或在离散的框架空间中进行。一种执行身体的动态的并行磁共振成像的方法，包括：从具有已知的或估算的灵敏度分布和噪声协方差矩阵的各个接收天线中获取所述身体的一组时间序列的基础磁共振图像，所述基础图像在k-空间中欠采样；以及执行所述身体的时间序列的磁共振图像的正则化重建。

Description

用于执行并行磁共振成像的方法

技术领域

本发明涉及一种用于执行人体的并行磁共振成像（pMRI）的方法，其包括并行的、动态的（时间分辨的）磁共振成像，例如功能性磁共振成像（fMRI）。

背景技术

在医学磁共振成像（MRI）中，尤其当涉及如fMRI的动态成像时，主要关心的是减少整体获取时间。实际上，短获取时间允许改进所获取的fMRI数据的空间分辨率/时间分辨率，这导致更有效的统计分析。此外，通过减少整体成像时间，可避免由患者移动所导致的一些额外的伪像。出于该原因，已经开发了并行成像系统：采用围绕处于下方的物体或身体而放置的多个具有互补敏感轮廓的接收器表面线圈，以同时采集频域内（即所谓的k-空间）的以沿着至少一个空间方向（即相位编码方向）的且比奈奎斯特（Nyquist）采样速率低R倍的速率所采样的数据，R通常称作“缩减因数”。因此，总获取时间比传统的非并行成像短R倍。然后，通过展开由各个接收器所获取的欠采样的“基础”图像而执行重建步骤以建立全视野（FOV）图像。由于与欠采样率有关的混叠伪像所导致的并行MRI（pMRI）中的低信噪比（SNR），这些混叠伪像由获取过程期间的噪音以及在线圈敏感场（sensitivity map）估计中误差的存在所引起，因此该重建是一项具有挑战性的任务。

空间谐波同步采集（SMASH）[Sodickson等，1997]是第一种重建方法，其在k-空间域操作。它使用预估算的线圈敏感场的线性组合来产生丢失的相位编码步骤。

还提出了一些其他的基于k-空间的重建技术，如GRAPPA（全局自动校准部分并行采集）[Griswold等，2002]和SENSE（灵敏度编码）[Pruessmann等，1999]。SENSE是首先依赖于缩小的FOV图像的重建、其次依赖于空间展开技术的两步式过程，其相当于加权最小平方估算法。该技术要求使用参考扫描（通常为2D梯度回波（GRE））的线圈敏感场的精确估算。目前，它是最频繁利用的pMRI技术，尤其被应用于脑成像和心脏成像。

关于pMRI中的重建方法的综述参见[Hoge等，2005]。

在具有无噪声数据和完美的线圈敏感场知识的情况下，SENSE通常被期望实现精确的重建，这对于全部上文所提到的方法也是可行的。然而，实际上，数据中噪声的存在以及线圈敏感场估算的不精确性是无法避免的且使重建问题情况较差。

由于图像重建是不适定的逆问题，因此正则化技术通常用于更好地估算全FOV图像。这些方法中的大多数工作在图像域中；对于Tikhonov正则化[Ying等，2004]尤其是这种情况，其使用二次惩罚项以促进平滑约束或说明重建图像和先验的参考图像之间的平方差。然而，尽管使用正则化，由于重建的图像受严重的混叠伪像影响，故当使用低的磁场强度（直到1.5Tesla）时，高的缩减因数（超过R=2的值）通常被认为是不可行的。

在等，2008]和

等，2009]中，发明人已经描述了一种用并行MRI来执行正则化图像重建的方法，其使用基于小波变换的正则化方案，从而允许增加缩减因数R。

发明内容

本发明目的在于提供所述方法的多种改进，包括使所述方法延伸至动态成像（例如fMRI）且使所述方法全部或部分自动校准（或“无人监管”）。

因此，本发明的一个目的在于一种身体的并行磁共振成像的方法，包括：

-从具有已知的或估算的灵敏度分布和噪声协方差矩阵的各个接收天线中获取所述身体的一组基础磁共振图像，所述基础图像在k-空间中欠采样；以及

-执行所述身体的磁共振图像的正则化重建；

其中，通过最小化成本函数而在离散的框架空间中进行所述执行磁共振图像的正则化重建的步骤，该成本函数包括：

-误差项，在给定所述所获取的基础图像的情况下，该误差项表示重建图像的似然；和

-框架惩罚项，该框架惩罚项表示所述重建图像的框架系数的实际统计分布与所述系数的先验分布之间的偏差；

基于所述身体的辅助的磁共振图像估算重建图像的框架系数的先验分布。

根据所述方法的不同实施方式：

-在给定所述所获取的基础图像的情况下，所述误差项可表示所述重建图像的负对数似然（neg-log-likelihood）。

-在给定所述所获取的基础磁共振图像和所述框架系数的先验分布的情况下，可通过在所述框架空间中最大化限定身体的图像的一组框架系数的完全后验分布而进行所述执行所述身体的磁共振图像的正则化重建的步骤。

-可根据所述所获取的基础磁共振图像重建所述身体的所述辅助的磁共振图像。更具体地，可使用灵敏度编码-SENS算法，根据所述所获取的基础磁共振图像重建所述身体的辅助的磁共振图像。甚至更具体地，可使用选自未正则化的SENSE算法、在图像空间中正则化的SENSE算法和在k-空间中正则化的SENSE算法中的一种算法，根据所述所获取的基础磁共振图像重建所述身体的辅助的磁共振图像。

-可假设所述框架系数的广义高斯-拉普拉斯先验统计分布，使用最大似然估算量或后验均值估算量，基于所述身体的所述辅助的磁共振图像估算所述分布的参数。

-所述误差项可为二次平均误差项。

-所述所获取的基础图像可为三维图像，且可在离散的三维框架空间中进行所述执行磁共振图像的正则化重建的步骤。更具体地，可以通过叠加待成像的对象的切片的二维基础图像来获得所述所获取的三维基础图像。

-所述执行磁共振图像的正则化重建的步骤可基于冗余小波框架表示。可替选地，所述执行磁共振图像的正则化重建的步骤可基于非冗余小波表示。

-所述成本函数可以还包括选自重建图像的全变差范数和凸约束中的至少一个空间域惩罚项。

本发明的另一目的在于一种执行身体的动态的并行磁共振成像的方法，该方法包括：

-从具有已知的或估算的灵敏度分布和噪声协方差矩阵的各个接收天线中获取所述身体的一组时间序列的基础磁共振图像，所述基础图像在k-空间中欠采样；以及

-执行所述身体的时间序列的磁共振图像的正则化重建；

其中，所述执行时间序列的基础磁共振图像的正则化重建的步骤通过最小化成本函数而进行，该成本函数包括：

-误差项，在给定相应的所获取的基础图像的情况下，所述误差项表示每个重建图像的似然；和

-时间惩罚项，该时间惩罚项表示所述时间序列的连续图像之间的像素间的差异或体素间的差异。

根据所述方法的不同实施方式：

-所述时间惩罚项可基于边缘保持函数，更具体地基于凸边缘保持函数，且甚至更具体地基于L_p范数，其中p≥1且优选地1≤p<1.5。

-可通过第一部分时间惩罚项和第二部分时间惩罚项的总和而给出所述时间惩罚项，其中：所述第一部分时间惩罚项表示序列的每个偶数图像与前一个奇数图像之间的像素间或体素间的差异；以及所述第二部分时间惩罚项表示序列的每个奇数图像与前一个偶数图像之间的像素间或体素间的差异；通过使用所述第一部分时间惩罚项和第二部分时间惩罚项的邻近算子最小化所述成本函数。

-可在离散的框架空间中进行所述执行磁共振图像的正则化重建的步骤，且所述成本函数可以还包括框架惩罚项，其表示每个重建图像的框架系数的统计分布与所述系数的先验分布之间的偏差；基于所述身体的辅助的磁共振图像估算重建图像的框架系数的所述先验分布。在这种情况下，所述基础图像可为三维图像，且所述离散的框架空间可为离散的三维框架空间。

-本方法可以还包括使用最大似然估算量来自动确定所述时间惩罚项的加权参数的步骤。更确切地说，本方法可包括估算待重建的图像的针对每个像素或体素、或者邻近像素或体素的组的时间惩罚项的所述加权参数。此外，所述时间惩罚项可基于L_p范数，且可使用最大似然估算量来共同地确定时间惩罚项的所述加权参数和参数p。特别地，在约束p≥1的情况下，可使用最大似然估算量来共同地确定时间惩罚项的所述加权参数和参数p。

-所述误差项可取决于几何参数，该几何参数限定了所述基础磁共振图像中每个图像相对于作为参照的基础磁共振图像的刚性变换，且其中，还通过相对于所述几何参数最小化所述函数来进行所述执行时间序列的基础磁共振图像的正则化重建的步骤。

-可通过回波平面成像获取所述基础图像。

-所述基础图像可以欠采样，且缩减因数大于或等于4。

附图说明

从结合附图而进行的以下描述中本发明的附加特征和优点将变得显而易见，附图示出：

-图1为人脑的磁共振图像的小波系数的实部和虚部的经验性的柱状图；

-图2和图3为使用现有技术中已知的Tikhonov-正则化的SENSE pMRI重建方法而获得的人脑的九个解剖的轴向切片，其中R分别等于2和4；

-图4为使用现有技术中已知的TV-正则化的SENSE pMRI重建方法而获得的同一人脑的解剖的轴向切片，其中R等于2（左面板）和R等于4（右面板）；

-图5和图6为使用根据本发明的实施方式的基于自动校准的2D小波变换的正则化方案（UWR-SENSE）而获得的同一人脑的九个解剖的轴向切片，其中R分别等于2和4；

-图7和图8为使用根据本发明的另一实施方式的基于受约束的、自动校准的2D小波变换的正则化方案（CWR-SENSE）而获得的同一人脑的九个解剖的轴向切片，其中R分别等于2和4；

-图9为使用根据本发明的另一实施方式的自动校准的组合的小波-全变差正则化方案且对非冗余标准正交的小波基使用图像分解而获得的同一人脑的九个解剖的轴向切片，其中R等于4；

-图10为使用根据本发明的另一实施方式的自动校准的组合的小波-全变差正则化方案且对由两个标准正交基的联合而构成的冗余小波框架使用的图像分解而获得的同一人脑的九个解剖的轴向切片，其中R等于4；

-图11为使用传统的非并行MRI（顶行）、R等于4的mSENSE并行MRI重建（中间行）和根据本发明另一实施方式的基于自动校准的3D小波变换的正则化方案（3D-UWR-SENSE）（底行）而获得的同一人脑的三个解剖的轴向切片；

-图12为使用传统的非并行MRI（顶行）、R等于4的mSENSE并行MRI重建（中间行）和根据本发明另一实施方式的基于自动校准的3D小波变换的正则化方案（3D-UWR-SENSE）（底行）而获得的同一人脑的三个解剖的矢状切片；

-图13为使用TV正则化（左）、根据本发明的实施方式的对非冗余小波基使用图像分解的组合小波-全变差正则化方案（中央）以及根据本发明的实施方式对由两个标准正交基的联合构成的冗余小波框架使用图像分解的组合小波-全变差正则化方案（右）而获得的同一人脑的三个解剖的轴向切片（上行）及其放大的细节（下行）；

-图14为使用回波平面（EPI）fMRI获取的，使用现有技术中已知的mSENSE方法和根据本发明的实施方式的基于自动校准的2D小波变换的正则化方案（4D-UWR-SENSE）重建的人脑的轴向切片、冠状切片和矢状切片；

-图15为使用EPI fMRI所检测的aC-aS对照的叠加到解剖图像的对象级的student-t图片；使用mSENSE、UWR-SENSE和4D-UWR-SENSE方法重建数据，其中R分别等于2（图的顶部）和4（图的底部）；示出矢状图、冠状图和轴向图；

-图16为分别使用mSENSE、UWR-SENSE和4D-UWR-SENSE方法且以R等于2重建的两个对象的aC-aS对照的对象级的student-t图片；示出矢状图、冠状图和轴向图；

-图17为使用EPI fMRI所检测到的Lc-Rc对照的叠加到解剖图像的对象级的student-t图片；使用mSENSE、UWR-SENSE和4D-UWR-SENSE方法重建数据，其中R分别等于2（图的顶部）和4（图的底部）；示出矢状图、冠状图和轴向图；

-图18为分别使用mSENSE、UWR-SENSE和4D-UWR-SENSE方法且以R等于4重建的两个对象的Lc-Rc对照的对象级的student-t图片；示出矢状图、冠状图和轴向图；

-图19为aC-aS对照的组级的student-t图片，其中使用mSENSE、UWRSENSE和4D-UWR-SENSE重建数据，且R等于2和4；示出矢状图、冠状图和轴向图；以及

-图20为Lc-Rc对照的组级的student-t图片，其中使用mSENSE、UWRSENSE和4D-UWR-SENSE重建数据，且R等于2和4；示出矢状图、冠状图和轴向图。

具体实施方式

在详细描述本发明之前，需要回顾关于pMRI（尤其SENSE）的一些基本情况。

MRI是一种成像技术，其可根据相关的RF脉冲设计以2D（二维）或直接以3D（三维）进行。在2D情况下，使用相邻的切片覆盖体积。由于该方法使k-空间扫描更快速，因此PMRI可适用于这两种情况，不论其清晰度如何（即以2D或3D）。为了简洁，在下文中仅示出2D情况。

利用一排C线圈来测量进入被研究的对象（例如在患者头部中的脑）中的自旋密度ρ。图像获取基于具体的成像序列；在本发明的示例性实施方式中，解剖的MRI基于3D MPRAGE序列，而所涉及的功能性MRI则使用2D回波平面成像（EPI）来执行。下面的描述将聚焦于2D情况；那么，每个线圈l所接收的信号d_c（1≤c≤C）为在k-空间中的一些位置k（k=(k_y,k_x)^T）（上标^T是指对换）上所估算的、由线圈灵敏度轮廓s_c加权的所期望的2D场ρ的傅里叶变换。因此，所接收的信号

由采样方案限定：

${\tilde{d}}_c\left(k_r\right)=\&Integral;\stackrel{\&OverBar;}{\mathrm{\&rho;}}\left(r\right)s_c\left(r\right)e^{-v2\mathrm{\&pi;}{k}_r^{T}r}\mathrm{dr}+{\tilde{n}}_c\left(k_r\right)---\left(1\right)$

其中，是加性高斯白噪声（AWGN）的体现，且r=(y,x)^T是图像域中的空间位置。为了简洁，通常采用笛卡尔坐标系。在其最简单的形式中，SENSE成像意味着解决了由于傅里叶变换的可分离性而导致的一维反演问题。使Δy=Y/R为采样时段，其中Y为沿着相位编码方向的视野（FOV）的尺寸，使y为在图像域中沿着相位编码方向的位置，x为在图像域内沿着频率编码方向的位置，且R≤L为缩减因数。2D傅里叶逆变换允许在空间域中恢复所测量的信号。通过用R来解释k-空间的欠采样，（1）可用以下矩阵形式重新表示：

$d\left(r\right)=S\left(r\right)\stackrel{\&OverBar;}{\mathrm{\&rho;}}\left(r\right)+n\left(r\right)---\left(2\right)$

其中：

和

在等式（2）中，(n(r))_r是循环的零均值高斯的复值随机向量的序列。这些噪声向量为i.i.d.（独立的、恒等分布的）且在空间上不受尺寸为C×C的协方差矩阵Ψ约束。实际上，通过从所有没有射频脉冲的线圈获取数据来估算Ψ，且其对应于两个线圈c₁和c₂之间的协方差的通用项Ψ(c₁,c₂)由下式给出：

$\mathrm{\&Psi;}(c_1,c_2)=\frac{1}{Y\&times;X}\underset{(y,x)}{\mathrm{\&Sigma;}}{\underset{\&OverBar;}{d}}_{c_1}(y,x){\underset{\&OverBar;}{d}}_{c_2}^{*}(y,x),\&ForAll;(c_1,c_2)\&Element;{\{1,\&CenterDot;\&CenterDot;\&CenterDot;,C\}}^2---\left(4\right)$

其中，(·)^*代表复共轭。注意，统计模型可被假定为矩阵Ψ（例如线圈间的统计独立性）或最邻近的统计相关模型。在第一种情况下，矩阵Ψ变成对角的且取决于线圈的方差可使用[Donoho,1995]所描述的平均绝对偏差（MAD）来估算，其依赖于加性高斯白噪声（AWGN）的假设。在最简单的情况下，Ψ可简单地等于单位矩阵。在存在非零的非对角项的情况下，在等式（4）中所给出的经验性协方差用于估算这些项。

因此，重建步骤包括使（2）反向和在空间位置r=(y,x)^T处从d（r）恢复

注意，数据(d_c)_1≤l≤C和未知图像

为复值的，尽管仅被认为用于可视化目的。

一种简单的重建方法，也称作SENSE方法[Pruessmann等，1999]，基于加权最小二乘法（WLS）准则的最小限度。

目的在于在每个空间位置r处找到向量

使得：

其中，

(.)^H代表转置的复共轭，(.)^#代表伪逆，且 ${\left|\right|\&CenterDot;\left|\right|}_{{\mathrm{\&Psi;\&psi;}}^{-1}}=\sqrt{{(\&CenterDot;)}^H\mathrm{\&Psi;\&psi;}(\&CenterDot;)}$ 限定

上的范数。

如上所述，该逆问题通常是不适定的且需要正则化，例如Tikhonov正则化。正则化过程通常包括计算

作为以下惩罚加权最小二乘法（PWLS）准则的最小值：

其中，I_R是R维的单位矩阵。正则化参数κ>0确保了数据的贴近度和惩罚项之间的平衡，其控制来自给定的参考向量ρ_r(r)的偏差。解允许以下封闭式表达：

${\widehat{\mathrm{\&rho;}}}_{\mathrm{PWLS}}\left(r\right)={\mathrm{\&rho;}}_r\left(r\right)+{(S^H\left(r\right){\mathrm{\&Psi;}}^{-1}S\left(r\right)+\mathrm{\&kappa;}{I}_R)}^{-1}{S}^H\left(r\right){\mathrm{\&Psi;}}^{-1}(d\left(r\right)-S\left(r\right){\text{\&rho;}}_r\left(r\right))---\left(7\right)$

注意，解的精确度取决于参考向量和正则化参数κ的选择。

已知Tikhonov正则化会引入图像中的模糊。

为了克服该限制，已经提出了“边缘保持”惩罚项，其被应用在图像域中且通过限制模糊效应和保持图像边缘而使正则化更加有效。然而，这些项的引入可导致在数学上不一定易于解决的不可微分的优化问题。

如上所述，在等，2008]和

等，2009]中，发明人已经提出了一种新的基于小波变换（WT）的正则化方案，将在这里概述该方案。

事实上，在基于SENSE的重建图像中，在空间上很好定位的伪像呈现为具有极高强度或极低强度的变形的曲线，且WT已经被认为是促使有用信息的良好的空间和频率定位的强大工具。

[Pesquet-Popescu,Pesquet]提出了一种小波和小波变换的大体介绍。

在下文中，T代表WT算子。它对应于在分辨率水平j_max下所进行的针对可分离的2D M-带小波基的离散分解。尺寸为Y×X的目标图像

可被视为欧几里得空间

的元素，其中K=Y×X且赋予标准内积<·|·>和范数||·||。

使(e_k)_1≤k≤K为所考虑的空间

的离散小波基。小波分解算子T被限定为线性算子：

用作重建目的的伴随算子T*被限定为双射线性算子：

所形成的目标图像函数ρ的小波系数场通过ζ=((ζ_a,k)_{1≤k≤Kjmax},(ζ_o,j,k)_{1≤j≤jmax,1≤k≤Kj})限定，其中，K_j=KM^-2j是在分辨率j下给定子带中的小波系数的数目（假设Y和X是M^jmax的倍数），并且通过某种方式重新索引这些系数，在这种方式中，ζ_a,k表示在分辨率水平j_max下的近似系数，ζ_o,j,k表示在分辨率水平j下且定向o∈O={0,...,M-1}²\{(0，0)}的细节系数。

在二元的情况下（M=2），具有对应于水平方向、垂直方向或对角线方向的三个取向。当考虑标准正交的小波基时，伴随算子T*减小至逆WT算子T^-1且T的算子范数||T||等于1。

通过重建小波算子T*生成目标图像

的估算量。使

为未知的小波系数，使得

目的在于从观察值d建立系数

的向量的估算量

这基于贝叶斯方法，其依赖于针对小波系数的合适的先验。

假定等式（2）中的观察模型以及关于噪声的假设（i.i.d.圆高斯具有零均值和线圈之间的相关矩阵Ψ），则似然函数对Y×X的FOV中的像素进行因式分解：

其中，ρ=T^*ζ且

使f为小波域中的图像的在先的概率密度函数（pdf）。在此假设小波系数的实部和虚部是独立的。还假设小波系数的实部（或虚部）在每个子带中是独立的且恒等分布的（i.i.d.）。然而，它们的统计特征可以在两个截然不同的子带之间变化。此外，通过查看所考虑的小波系数的实部和虚部的经验分布，本发明人已经注意到它们的经验柱状图与“广义Gauss-Laplace”（GGL）分布很好地匹配，这说明单个模式及其形状在高斯密度和拉普拉斯算子密度之间变化。对应的pdf为：

其中，

和

为待估算的超参数。图1示出了使用利用长度为8的多贝西（Daubechies）滤波器的二元（M=2）小波分解在第一分辨率水平下的水平的具体子带的实部和虚部的经验柱状图。该图还示出所采用的GGL分布比广义高斯（GG）pdf更好地匹配经验柱状图。

在最粗略的分辨率水平j_max下，由于近似系数的实部和虚部都属于低频子带，因此近似系数的实部和虚部的分布可通过高斯分布来进行建模。

由于MAP（最大后验概率）估算量的相似性和简明性，MAP估算量将被用于估算目的，尽管可能有不同的选择。基于上文给出的先验和似然，通过最大化全部后验分布或最小化其负对数似然函数来计算MAP估算量：

或者等效地，通过最小化以下准则来计算：

其中：

${\mathrm{\&Phi;}}_a\left({\mathrm{\&zeta;}}_{a,k}\right)=\frac{{(\mathrm{Re}\left({\mathrm{\&zeta;}}_{a,k}\right)-{\mathrm{\&mu;}}^{\mathrm{Re}})}^2}{{\left(\sqrt{2}{\mathrm{\&sigma;}}^{\mathrm{Re}}\right)}^2}+\frac{{(\mathrm{Im}\left({\mathrm{\&zeta;}}_{a,k}\right)-{\mathrm{\&mu;}}^{\mathrm{Im}})}^2}{{\left(\sqrt{2}{\mathrm{\&sigma;}}^{\mathrm{Im}}\right)}^2},$

${\mathrm{\&Phi;}}_{o,j}\left({\mathrm{\&zeta;}}_{o,j,k}\right)={\mathrm{\&alpha;}}_{o,j}^{\mathrm{Re}}\left|\mathrm{Re}\left({\mathrm{\&zeta;}}_{o,j,k}\right)\right|+\frac{{\mathrm{\&beta;}}_{o,j}^{\mathrm{Re}}}{2}{\left|\mathrm{Re}\left({\mathrm{\&zeta;}}_{o,j,k}\right)\right|}^2+{\mathrm{\&alpha;}}_{o,j}^{\mathrm{Im}}\left|\mathrm{Im}\left({\mathrm{\&zeta;}}_{o,j,k}\right)\right|$ $+\frac{{\mathrm{\&beta;}}_{o,j}^{\mathrm{Im}}}{2}{\left|\mathrm{Im}\left({\mathrm{\&zeta;}}_{o,j,k}\right)\right|}^2.---\left(13\right)$

在上文中，Re(·)和Im(·)(或·^Re和·^Im)分别代表实部和虚部。先验参数（超参数）

和都是未知的且需要被估算。超参数的估算是本发明的重要部分且将在下文中对其展开讨论。

当J_L利用李普希茨（Lipschitz）连续梯度是可微的时，J_P是不可微的。因此，尽管罚函数J_WT是凸面的，但是优化过程不能够依赖于传统的凸优化技术，例如同伪共轭梯度。

可使用在[Daubechies等，2004；Chaux等，2007]中所开发的迭代优化过程的广义形式来进行优化，其基于向前-向后（FB）算法和邻近算子的概念，被概括成复变函数的情形。

对于函数

其中，φ^Reandφ^Im是

中的函数，Re(x)(或者Im(x))是的分量的实部（或者虚部）的向量，邻近算子被限定为：

通过将[Chaux等，2007]中的算法扩展到复数情况，可根据以下算法（算法1）迭代地计算J_WT的最小值，其中最先计算J_L的梯度，然后更新框架系数。可注意到，λ_n和γ_n分别对应于松弛参数和步长参数。

算法1：

使(γ_n)_n>0和(λ_n)_n>0为正实部的序列。

1：设定n=0且ε≥0。初始化ζ⁽ⁿ⁾且设定

2：重复。

3：通过设定ρ⁽ⁿ⁾=T^*ζ⁽ⁿ⁾而重建图像。

4：计算图像μ⁽ⁿ⁾使得：

${\&ForAll;}_r\&Element;\{1,...,Y/R\}\&times;\{1,...,X\},$

μ⁽ⁿ⁾(r)=2S^H(r)Ψ^-1(S(r)ρ⁽ⁿ⁾(r)-d(r))，

其中，向量μ⁽ⁿ⁾(r)以与相同的方式由μ⁽ⁿ⁾限定，由

限定（参见等式（2.17））。

5：确定μ⁽ⁿ⁾的小波系数

6：更新重建图像的ρ⁽ⁿ⁺¹⁾的邻近系数：

$\&ForAll;k\&Element;\{1,...,K_{j_{\max }}\},{\mathrm{\&zeta;}}_{a,k}^{(n+1)}={\text{\&zeta;}}_{a,k}^{\left(n\right)}+{\mathrm{\&lambda;}}_n({\mathrm{prox}}_{{\mathrm{\&gamma;}}_n{\mathrm{\&Phi;}}_n}({\mathrm{\&zeta;}}_{a,k}^{\left(n\right)}-r_n{v}_{a,k}^{\left(n\right)})-{\mathrm{\&zeta;}}_{a,k}^{\left(n\right)}).$

7：更新重建图像的ρ⁽ⁿ⁺¹⁾的细节系数：

8：计算

9：n←n+1。

10：直到

11:返回ρ⁽ⁿ⁾=T^*ζ⁽ⁿ⁾。

对于每个r∈{1,…,Y}×{1,…,X}，使0_r≥0为厄米半正定矩阵S^H(r)Ψ^-1S(r)的最大本征值且使θ=maxr∈{1,…,Y}×{1,···,X}θ_r＞0。为了保证算法1的收敛，步长参数和松弛参数服从于以下条件:

(i)inf_n＞0γ_n＞0和 ${\sup }_{n>0}{\mathrm{\&gamma;}}_n\text{<}\frac{1}{\mathrm{\&theta;}{\left|\right|T\left|\right|}^2},$

(ii)inf_n＞0λ_n＞0和sup_n＞0λ_n≤1.

如

等，2009]中所述，优化问题（11）、（12）可被修改成包括额外的（优选地，凸的）约束，例如图像强度值的上界和下界。

本发明人已经发现使用如上文所述小波正则化允许保持图像细节同时使重建伪像平滑，但可能在图像的均匀区域中引起少许不规则。另一方面，已知其它的（不基于小波的）正则化方案适用于正则化平滑区域，但以图像细节的过度平滑为代价。尤其在[Raj等，2007]和[Block等，2007]中所描述的“全变差”TV正则化正是这种情况。

等，2008]和

等，2009]的方法的第一种改进，构成本发明的一个方面，其在于在联合正则化框架中使WT和TV相结合以利用它们的不同属性，从而允许它们缓解彼此的缺陷。联合小波变换-全变差(WTTV)的正则化可归纳为以下惩罚准则的最优化：

其中，J_L和J_P如上文所限定，κ₁>0和κ₂>0为正则化参数，以及T^*是小波伴随算子。图像

的离散全变差范数由下式给出：

${\left|\right|\mathrm{\&rho;}\left|\right|}_{\mathrm{TV}}=\underset{y=1}{\overset{Y}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{X}{\mathrm{\&Sigma;}}}\sqrt{{\left|\left({\&dtri;}_1\mathrm{\&rho;}\right)(y,x)\right|}^2+{\left|{\left({{\&dtri;}_1\mathrm{\&rho;}}^T\right)}^T(y,x)\right|}^2}---\left(16\right)$

其中，对于每个

是水平平滑的梯度算子，其由下式限定：

${\&dtri;}_1\left(\mathrm{\&rho;}\right)=\frac{1}{2}(\mathrm{\&rho;}(y+1,x+1)-\mathrm{\&rho;}(y,x+1)+\mathrm{\&rho;}(y+1,x)-\mathrm{\&rho;}(y,x){)}_{1\&le;y\&le;Y,1\&le;x\&le;X}---\left(17\right)$

为了说明，假设图像ρ是周期性的或等效地，图像边界是环形的。

最小化最优性准则（15）比最小化式（11）至式（13）所给出的仅空间的准则更加困难，因为涉及超过两项，然而上文所讨论的向前-向后（FB）算法仅适用于包含两项的准则的最小化。包含超过两项的不可微分的凸函数的迭代最小化可使用所谓的并行邻近算法（PPXA）[Combettes and Pesquet, 2008]来执行，其也需要计算三个所涉及的项中每项的邻近算子。

困难源于||·||_TVοT^*的邻近算子的计算。为了克服该困难，等式(15)中的TV惩罚项可被分成如[Combettes and Pesquet, 2008]中的四项。因此，TV 惩罚项为：

${\left|\right|\mathrm{\&rho;}\left|\right|}_{\mathrm{TV}}=\underset{i=0}{\overset{3}{\mathrm{\&Sigma;}}}{\mathrm{tv}}_i\left(\mathrm{\&rho;}\right)$

其中，对于每个(q,r)∈{0,1}²，

${\mathrm{tv}}_{2q+r}\left(\mathrm{\&rho;}\right)=\underset{y=1}{\overset{Y/2}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{X/2}{\mathrm{\&Sigma;}}}\sqrt{{\left|\left({\&dtri;}_1\mathrm{\&rho;}\right)(2y+q,2x+r)\right|}^2+{\left|{\left({\&dtri;}_1{\mathrm{\&rho;}}^T\right)}^T(2y+q,2x+r)\right|}^2}$

对于{0,1}中的每个q和r，使↓q,r为抽取算子，其由下式限定：

且U_q+2r为以下算子：

其中，对于每个算子

和

由下式限定：

${\&dtri;}_0\left(\mathrm{\&rho;}\right)=\frac{1}{2}{(\mathrm{\&rho;}(y+1,x+1)+\mathrm{\&rho;}(y,x+1)+\mathrm{\&rho;}(y+1,x)+\mathrm{\&rho;}(y,x))}_{1\&le;y\&le;Y,1\&le;x\&le;X}$ 和

${\&dtri;}_2\left(\mathrm{\&rho;}\right)=\frac{1}{2}{(\mathrm{\&rho;}(y+1,x+1)-\mathrm{\&rho;}(y,x+1)-\mathrm{\&rho;}(y+1,x)+\mathrm{\&rho;}(y,x))}_{1\&le;y\&le;Y,1\&le;x\&le;X}.$

也使h为通过以下公式，根据所限定的函数：

$h\left(\mathrm{\&rho;}\right)=\underset{y=1}{\overset{Y/2}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{X/2}{\mathrm{\&Sigma;}}}\sqrt{{\left|\left({\&dtri;}_1\mathrm{\&rho;}\right)(y,x+X/2)\right|}^2+{\left|{\left({\&dtri;}_1{\mathrm{\&rho;}}^T\right)}^T(y+Y/2,x)\right|}^2}$

从以上等式可得出，

tv_i=hοU_i。因此，在等式（15）中的优化问题可被重写成：

正则化的WT-TV重建方法用算法2概括，其中，PPXA算法[Combettes和Pesquet，2008]用于最小化等式（18）中的最优性准则，其由J=6个凸函数构成。

算法2

设定γ∈[0,+∞]，n=0，(ω_i)_1≤i≤6∈[0,1]⁶，使得

其中，对于每个i∈{1,...,6}，

也设定ε≥0，初始化 ${\mathrm{\&zeta;}}^{\left(n\right)}={\mathrm{\&Sigma;}}_{i=1}^6{\mathrm{\&omega;}}_i{\mathrm{\&zeta;}}_i^{\left(n\right)}$ 和

1：重复。

2：计算图像μ⁽ⁿ⁾使得：

${\&ForAll;}_r\&Element;\{1,...,Y/R\}\&times;\{1,...,X\},$

${\mathrm{\&mu;}}^{\left(n\right)}\left(r\right)={(I_R+\frac{2\mathrm{\&gamma;}}{{\mathrm{\&omega;}}_1}{S}^H\left(r\right){\mathrm{\&Psi;}}^{-1}S\left(r\right))}^{-1}({\mathrm{\&rho;}}_1^{\left(n\right)}\left(r\right)+\frac{2\mathrm{\&gamma;}}{{\mathrm{\&omega;}}_1}{S}^H\left(r\right){\mathrm{\&Psi;}}^{-1}d\left(r\right)),$ 其中 ${\mathrm{\&rho;}}_1^{\left(n\right)}=T^{*}{\mathrm{\&zeta;}}_1^{\left(n\right)}.$

3：计算μ⁽ⁿ⁾的小波系数p₁ ⁽ⁿ⁾=Tu⁽ⁿ⁾。

4：计算

5：对于每个i∈{0,1,2,3}，计算 $p_{i+3}^{\left(n\right)}=\left({\mathrm{prox}}_{\mathrm{\&gamma;}{K}_2{\mathrm{\&omega;}}_i/{\mathrm{\&omega;}}_i}\right({\mathrm{\&zeta;}}_{i+3}^{\left(n\right)}).$

6：设定 $P^{\left(n\right)}={\mathrm{\&Sigma;}}_{i=1}^6{\mathrm{\&omega;}}_i{p_i}^{\left(n\right)}.$

7：设定λ_n∈|0,2|。

8：对于i=1至6，进行

$9:{\mathrm{\&zeta;}}_i^{\left(n\right)}={\mathrm{\&zeta;}}_i^{\left(n\right)}+{\mathrm{\&lambda;}}_n(2P^{\left(n\right)}-{\mathrm{\&zeta;}}^{\left(n\right)}-p_i^{\left(n\right)}).$

10：比对

11：ζ⁽ⁿ⁺¹⁾=ζ⁽ⁿ⁾+λ_n(P⁽ⁿ⁾-ζ⁽ⁿ⁾)。

12：计算

13：n←n+1。

14：直到

15：返回ρ⁽ⁿ⁾=T^*ζ⁽ⁿ⁾。

在

等，2008]和等，2009]的方法中,通过叠加正则化的二维（2D）图像而进行三维（3D）成像。该方法的另一改进，构成本发明的另一方面，其在于执行直接的正则化的3D图像重建。在3D情况下，图像重建问题仍可被写成：

$d\left(r\right)=S\left(r\right)\stackrel{\&OverBar;}{\mathrm{\&rho;}}\left(r\right)+n\left(r\right)---\left(19\right)$

如在等式（12）中；然而，现在r=(y,x,z)是三维的空间位置，z∈{1,...,Z}是沿着第三方向（切片选择）的位置。此外，惩罚项J_P将取决于3D二元小波变换的小波系数的分布。

该方法的另一改进，构成本发明另一方面，其在于用动态的MRI执行四维（4D）的时空正则化。本发明的该实施方式可被应用至例如脑的fMRI，其中，全部脑体积不得不被获取数次，以产生4D数据集。在传统的fMRI中，3D图像被假定为独立的，尽管它们属于相同的fMRI时段。然而，实际上，3D时间图像是以某种方式而独立的，由于它们属于涉及相同实验范例的相同的fMRI时段。在脑的fMRI时段中，BOLD（血氧水平依赖）时间序列和获取噪声实际上是时间相关的。出于该原因，考虑到3D图像之间的时间依赖性有助于增加穿过所获取的图像的SNR，因此增强fMRI中的统计分析的可靠性。然而，由于在动态的MRI中，成像的对象几何形状通常在获取期间发生变化，因此将重建过程与时间正则化结合是非常困难的。

为了处理N_r（(N_r通常是偶数）个三维图像的时间序列的4D重建，等式（12）和等式（19）的观察模型将被重写成如下形式：

d_t(r)=S(r)ρ_t(r)+n_t(r)  (20)

其中，t∈{1,…,N_r}，是获取时间。使用二元的3D小波算子T，系数将被重新索引，使得

其中

说明额外的时间正则化项l_p，4D“体积”的（即3D图像或体积的时间序列的）重建通过以下最优性准则的最小化而进行：

其中，κ>0是正则化参数且J_P被限定为在式（13）中。

在等式（21）中，鉴于对应的所获取的基础图像， $\underset{t=1}{\overset{N_r}{\mathrm{\&Sigma;}}}\underset{r\&Element;\{1,...,Y/R\}\&times;\{1,...,X\}\&times;\{1,...,Z\}}{\mathrm{\&Sigma;}}{\left|\right|d_t\left(r\right)-S\left(r\right)\left(T^{*}{\mathrm{\&zeta;}}^t\right)\left(r\right)\left|\right|}_{{\mathrm{\&Psi;}}^{-1}}^2$ 是表示每个重建的图像的似然的“误差项”；

是时间惩罚项，其表示序列的连续图像之间的像素间或体素间的差异，以及

是用于空间正则化的小波惩罚项。应当理解，也可使用不同的（例如不基于小波的）空间正则化项。

于是将小波伴随算子T^*应用至ζ的每个分量ζ^t，以通过考虑与其他所获取的图像的时间依赖性而获得在获取时间t处所重建的3D图像ρ^t。

注意到，时间惩罚项不基于小波，但表示序列的连续图像之间的像素间或体素间的差异。更确切地说，在等式（21）中，时间惩罚项基于L_p范数；本发明人已经发现p必须大于或等于1，且p应优选地满足：1≤p<1.5。可使用其它形式的时间惩罚项，尤其基于序列的连续图像之间的像素间或体素间的差异的边缘保持函数（优选地为凸函数）。

最小化最优性准则（21）比最小化式（11）至式（13）所给出的仅空间的准则更加困难，因为涉及多于两项，而上文所述的向前-向后（FB）算法仅适用于包含两项的准则的最小化。包含超过两项的不可微分的凸函数的迭代最小化可使用已经引证的并行邻近算法（PPXA）[Combettes and Pesquet,2008]而执行，该算法也需要计算三个所涉及的项中每项的邻近算子。由于等式（21）的两个第一项相对于时间变量t和空间位置r是独立的，因此该任务对于他们是非常简单的。然而，对于时间惩罚项（等式（21）中的第三项）则并非如此，其使得对应的邻近算子的计算变得稍微复杂。因此，建议通过将时间惩罚项分解成相对于时间变量t（在第一项或第二项中涉及给定的获取时间t）独立的两项而以J_ST重写最优性准则，且对于该两项，邻近算子容易计算。更具体地，将时间惩罚项J_T表达为第一部分时间惩罚项J_T ¹和第二部分时间惩罚项J_T ²的总和，其中：第一部分时间惩罚项表示序列的每个偶数号的图像和前一个奇数号的图像之间的像素间或体素间的差异；第二部分时间惩罚项表示序列的每个奇数号的图像和前一个偶数号的图像之间的像素间或体素间的差异：

由于J_T ¹和J_T ²相对于时间变量t是可分离的，因此可以基于所涉及的项中的每项的邻近算子而很容易地计算对应的邻近算子。考虑以下函数：

其中，且H为由下式限定的线性算子：

因此，它的相关的伴随算子H^*由下式给出：

那么，Φ的邻近算子由下式给出：

在下文中给出用于时空最优性准则J_ST的最小化的形成算法（算法3）。

算法3：

设定γ∈[0,+∞]，n=0，(ω_i)_1≤i≤4∈[0,1]⁴，使得

其中，对于每个i∈{1,...,4}和t∈{1,...,N_r}，

和

也设定ε≥0，初始化 ${\mathrm{\&zeta;}}^{\left(n\right)}={\mathrm{\&Sigma;}}_{i=1}^4{\mathrm{\&omega;}}_i{\mathrm{\&zeta;}}_i^{\left(n\right)}$ 和

1：重复。

2：设定 $p_4^{1,\left(n\right)}={\mathrm{\&zeta;}}_4^{1,\left(n\right)}.$

3：对于t=1至Nr，进行

4：计算图像μ_t ⁽ⁿ⁾使得：

${\&ForAll;}_r\&Element;\{1,...,Y/R\}\&times;\{1,...,X\}\&times;\{1,...,Z\},$

${{\mathrm{\&mu;}}_t}^{\left(n\right)}\left(r\right)={(I_R+\frac{2\mathrm{\&gamma;}}{{\mathrm{\&omega;}}_1}{S}^H\left(r\right){\mathrm{\&Psi;}}^{-1}S\left(r\right))}^{-1}({\mathrm{\&rho;}}_{t,1}^{\left(n\right)}\left(r\right)+\frac{2\mathrm{\&gamma;}}{{\mathrm{\&omega;}}_1}{S}^H\left(r\right){\mathrm{\&Psi;}}^{-1}d\left(r\right)),$ 其中 ${\mathrm{\&rho;}}_{t,1}^{\left(n\right)}=T^{*}{\mathrm{\&zeta;}}_1^{t,\left(n\right)}.$

5：计算小波系数

6：计算

7：如果t是偶数，则

8：计算 $(p_3^{t,\left(n\right)},p_3^{t-1,\left(n\right)})={\mathrm{prox}}_{\mathrm{\&gamma;k\&Phi;}/{\mathrm{\&omega;}}_3}\left(({\mathrm{\&zeta;}}^{t,\left(n\right)},{\mathrm{\&zeta;}}^{t-1,\left(n\right)})\right),$

9：如果t是奇数且t>1，则

10：计算 $(p_4^{t,\left(n\right)},p_4^{t-1,\left(n\right)})={\mathrm{prox}}_{\mathrm{\&gamma;k\&Phi;}/{\mathrm{\&omega;}}_4}\left(({\mathrm{\&zeta;}}^{t,\left(n\right)},{\mathrm{\&zeta;}}^{t-1,\left(n\right)})\right).$

11：结束条件

12：如果t>1，则

13：设定 $P^{t-1,\left(n\right)}={\mathrm{\&Sigma;}}_{i=1}^4{\mathrm{\&omega;}}_i{p}_i^{t-1,\left(n\right)}.$

14：结束条件

15：比对

16：设定 $p_4^{N_r,\left(n\right)}={\mathrm{\&zeta;}}_4^{N_r,\left(n\right)}.$

17：设定 $P^{N_r,\left(n\right)}={\mathrm{\&Sigma;}}_{i=1}^4{\mathrm{\&omega;}}_i{p}_i^{N_r,\left(n\right)}.$

18：设定 $p_1^{\left(n\right)}=(p_1^{1,\left(n\right)},p_1^{2,\left(n\right)},...,p_1^{N_r,\left(n\right)}),$ $p_2^{\left(n\right)}=(p_2^{1,\left(n\right)},p_2^{2,\left(n\right)},...,p_2^{N_r,\left(n\right)}),$ $p_3^{\left(n\right)}=(p_3^{1,\left(n\right)},p_3^{2,\left(n\right)},...,p_3^{N_r,\left(n\right)}),$ $p_4^{\left(n\right)}=(p_4^{1,\left(n\right)},p_4^{2,\left(n\right)},...,p_4^{N_r,\left(n\right)})$ 和 $p^{\left(n\right)}=(p^{1,\left(n\right)},p^{2,\left(n\right)},...,p^{N_r,\left(n\right)}).$

19：设定λ_n∈[0,2]。

20：对于i=1至4，进行

$21:{\mathrm{\&zeta;}}_i^{\left(n\right)}={\mathrm{\&zeta;}}_i^{\left(n\right)}+{\mathrm{\&lambda;}}_n(2P^{\left(n\right)}-{\mathrm{\&zeta;}}^{\left(n\right)}-p_i^{\left(n\right)}).$

22：比对

23：ζ⁽ⁿ⁺¹⁾=ζ⁽ⁿ⁾+λ_n(P⁽ⁿ⁾-ζ⁽ⁿ⁾)。

24：计算

25：n←n+1。

26：直到

27：设定 $\widehat{\mathrm{\&zeta;}}={\mathrm{\&zeta;}}^{\left(n\right)}.$

28：对于每个t∈{1,...,N_r}，返回

应当理解，时空正则化也可适用于二维图形的时间序列的情况，即使这样不太常见。此外，空间惩罚项可以还包括全变差分量和/或约束量，如在纯的空间2D或3D正则化的情况下。

值得注意的是，在4D fMRI数据集中，由于成像主体的体积被获取多次，因此所述主体（即脑）可在连续扫描之间轻微移动。对于fMRI分析，固有的运动伪像可以是引人注目的。在这种背景下，由于标准的fMRI研究涉及通过应用刚性变换的运动校正步骤，因此提出了估算方法的另一延伸。该方法通过在估算步骤中引入刚性变换参数δr来解释最终的运动伪像。变换首先被应用至将用于估算步骤的标准的SENSE数据。还提出了上文所述方法的延伸，其不是对标准的SENSE重建数据，而是对二次正则化（QR）版本的SENSE（即QR-SENSE）进行估算步骤，QR-SENSE使用具有平滑矩阵的Tikhonov正则化。QR-SENSE对于改进估算步骤的稳健性比在SENSE解决方案上引入额外的滤波步骤表现得更有效率，且由于QR-SENSE允许闭型表达式，因此其保持了计算上的高效率。为了进一步考虑发生在fMRI中的运动伪像，所提出的4D重建方法（算法3）也可被扩展为解释重建过程期间的刚性变换参数。因此，在等式（22）中的最优性准则可被重写为以下形式：

其中，

来源于在应用运动校正变化之后的ζ^t，且假设：

和

由于，

或者

相对于时间变量t是可分离的，因此其邻近算子可基于在上文的总和中所涉及的项的每项的邻近算子而被很容易地计算。

也相对于这些参数，刚性变换参数δr可通过最小化等式（22bis）所表达的

而被发现。可替代地，可以基于传统的SENSE重建，单独地计算所述刚性变换参数。

本文中已经描述的全部正则化方法（基于小波的2D空间正则化、基于小波的3D空间正则化、基于混合小波/全变差的2D和3D空间正则化、3D和4D时空正则化）涉及必须使用合适的估算算法手动设定或自动设定的数值参数。这些参数包括小波系数的统计分布的先验参数或“超参数”：

和它们也可包括时间正则化参数κ，所谓的“形状参数”，即在时间正则化中所使用的限定L_p范数||·||_p的“p”的值（参见等式21），以及确定了在混合正则化中小波和全变差惩罚项的相对权重的参数κ₁、κ₂。

已引证的文献

等，2008]和

等，2009]没有教导如何确定有关的参数。这些文献所描述的方法的进一步改进，构成本发明的另一方面，其在于提供自动校准算法，从而允许全部或部分的无监督的基于小波的正则化。两个不同的框架被认为根据分解的类型而实现该目的。在标准正交的小波分解的情况下，可以使用最大似然法（参见下文）来精确地估算超参数。相比之下，当借助于冗余框架分解时，随机采样方法似乎更为有效。在使用随机采样法的情况下，使用最小均方差（MMSE）或者等效的后验均值估算量（参见下文）来估算超参数。后验均值估算量也可与标准正交的小波分解一起使用，如同这是框架的特例。

首先，将考虑确定限定用于2D或3D空间正则化的小波系数的先验统计分布的超参数。

用于确定超参数

的第一种可能值是根据以下等式最大化完整的似然：

$\widehat{\mathrm{\&Theta;}}=\underset{\mathrm{\&Theta;}}{\arg \max }p(d;\mathrm{\&Theta;})=\underset{\mathrm{\&Theta;}}{\arg \max }\&Integral;p\left(d\right|T^{*}\mathrm{\&zeta;})f(\mathrm{\&zeta;};\mathrm{\&Theta;})\mathrm{d\&zeta;}---\left(24\right)$

由于ζ是未知的，因此最大化等式（24）是遗漏数据问题。它需要整合出所寻找的图像分解ζ且使用[Dempster等，1977]所描述的密集型EM算法在图像重建和超参数估算之间进行迭代。为了减轻计算负担，通过假设参考的或“辅助的”全FOV图像

是可用的来不同地进行是有利的，因此其小波分解

也是如此。实际上，辅助图像

可以按照图像空间中或者k-空间中的传统的方式（例如使用Tikhonov正则化），使用R值相同的未正则化的或二次正则化的1D-SENSE重建而获得。

因此，最大似然（ML）估算方法包括假设该辅助图像是全部先验分布的实现，从而基于其直接拟合Θ：

在统计学上，该方案被称为完整数据最大似然，其与上文提到的遗漏数据最大似然估算量相反。该方法可在两个独立的步骤中进行分解，第一步骤包括设定隶属于近似系数

的高斯先验参数（μ，σ），第二步骤涉及估算来自对应的细节系数的GGL先验参数

另一方面，通过经验均值和标准偏差而明确地给出ML估算量

${\widehat{\mathrm{\&mu;}}}^{\mathrm{Re}}=\frac{1}{k_{j\max }}\underset{k=1}{\overset{K_{j\max }}{\mathrm{\&Sigma;}}}\mathrm{Re}({\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_a,k),{\widehat{\mathrm{\&sigma;}}}^{\mathrm{Re}}=\sqrt{\frac{1}{K_{j\max }}\underset{k=1}{\overset{K_{j\max }}{\mathrm{\&Sigma;}}}{(\mathrm{Re}\left({\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_{a,k}\right)-{\widehat{\mathrm{\&mu;}}}^{\mathrm{Re}})}^2}$

${\widehat{\mathrm{\&mu;}}}^{\mathrm{Im}}=\frac{1}{k_{j\max }}\underset{k=1}{\overset{K_{j\max }}{\mathrm{\&Sigma;}}}\mathrm{Im}\left({\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_{a,k}\right)\mathrm{and}{\widehat{\mathrm{\&sigma;}}}^{\mathrm{Im}}=\sqrt{\frac{1}{K_{j\max }}\underset{k=1}{\overset{K_{j\max }}{\mathrm{\&Sigma;}}}{(\mathrm{Re}\left({\stackrel{\&OverBar;}{\mathrm{\&zeta;}}}_{a,k}\right)-{\widehat{\mathrm{\&mu;}}}^{\mathrm{Im}})}^2}.---\left(25\right)$

对于每个分辨率水平j和定向o，按照下式从估算

和

超参数

通过在等式（26）中用Im(·)替代Re(·)而以相同方式被估算。

这种两维最小化问题并不允许闭型解决方案。因此，使用数值最优化方法，例如直接搜索法（例如，罗森布鲁克（Rosenbrock）法，参见[Bertsekas,1995,第一章,p159-165]），来计算ML参数。基于蒙特卡罗（Monte Carlo）法或斯坦（Stein）原理的可替代方案也可被考虑，但以增加计算负担为代价。

当考虑时间正则化时，必须确定两个额外的参数：时间正则化参数κ和形状参数p(参见等式21和等式22)。

实际上，在对象级的分析中，仅处理为一个对象所记录的数据。出于该原因，考虑到合理的数据量，手动地固定正则化参数κ和形状参数p可能是适宜的，即使这使所提出的方法并不完全地自动化且在某种程度上保持用户依赖性。当考虑组级的分析时，所处理的数据变得更大，且由于许多对象被考虑，因此因对象间差异会导致同时为所有对象手动设定的时间正则化参数肯定是不理想的。该差异从解剖学观点或功能性观点来看可能是显著的。为每个分开的对象手动设定该参数是可能的但却是麻烦的，且考虑到在组级的分析中所涉及的对象的数目（通常涉及15个对象），这会使所提出的方法变得更具有用户依赖性。

甚至在组级的分析中，针对每个分开的对象进行估算。另外，优选地，针对3D体积的每个体素估算参数值，这意味着每个体素独立于其他体素而处理。尽管如此，出于计算效率的原因，同一形状参数p被认为用于整个体积。除非另有说明，为每个体素确定各自的参数κ_r，而不是确定唯一的时间正则化参数κ。

再次，执行所述参数的最大似然估算。使ρ¹(r)…ρ^Nr(r)为形成4D数据集的3D图像，且将V_r限定为：

$V_r=({\mathrm{\&rho;}}^2\left(r\right)-{\mathrm{\&rho;}}^1\left(r\right),{\mathrm{\&rho;}}^3\left(r\right)-{\mathrm{\&rho;}}^2\left(r\right),...,{\mathrm{\&rho;}}^{2N_r}\left(r\right)-{\mathrm{\&rho;}}^{N_{r-1}}\left(r\right))$

因此，在ML意义上的参数p和参数κ_r的联合估算可通过最小化以下准则而执行：

$L(p,k)\&Proportional;-\log \left(k_r\right)+k_r{\mathrm{\&Sigma;}}_{t=1}^{N_r-1}{\left|V_r\left(t\right)\right|}^p---\left(27\right)$

优选地，约束条件p≥1。

显然，给定的体素应当更可能地与其相邻的体素行为相似，与解剖学观点和功能性观点一样。可替代的方法包括估算脑的每个功能区（即具有相似功能作用的一组相邻的体素）的参数值。该估算基于脑的先验分类步骤以识别不同的功能区。由于在给定的估算步骤中，较多的样品是可用的，因此该延伸使估算过程更加稳健。实际上，与所有属于同一功能区的体素有关的时间信号被集中在同一向量中。如果用S_ι指示功能区ι的尺寸，则对包括S_ι倍的更多样品的信号进行最大似然估算。因此，所产生的估算对于非典型的观察值（异常值）更加稳健，其可发生在fMRI时段中。

等，2008]和

等，2009]的图像重建方法基于小波基上所观察到的体积的分解。构成本发明的另一方面的另一改进，在于将所述方法推广至对框架进行分解的情形。

框架概括了小波、曲波、带波（bandlet）等的概念。框架是函数族，信号（或图像，图像仅仅是多维信号）可被分解到该函数族上。通常，所形成的信号表示是冗余的；当框架不是冗余的时，框架被称为基（basis）。框架通常被用于协助捕获处于研究的图像的更多的几何细节。出于该原因，在pMRI重建中使用框架有助于追踪重建伪像。

更确切地说，认为长度为L的实值的数字信号是欧几里得空间

的元素，该空间被赋予通常表示为<·|·>的数积和表示为||·||的范数。

使K为大于或等于L的整数。当]0,+∞[中存在常数μ时，在有限维空间

中的向量族(e_k)_1≤k≤K是框架，使得：

如果不等式（28）变成等式，则被称为(e_k)_1≤k≤K紧框架。有界的线性框架分析算子F和伴随合成框架算子F*被限定为：

注意，F是单射的，而F*是满射的。当F^-1=F*时，(e_k)_1≤k≤K是标准正交基。冗余框架的简单示例是标准正交基的联合。在该情况下，框架是紧密的且μ=M，因而F*F=MI，其中I是恒等算子。

在下文中，将解决复杂性增大的两个框架估算问题，第一个框架估算问题对应于降噪模式，第二个框架估算问题对应于涉及退化算子的扩展例如在并行MRI重建中的灵敏度矩阵S。首先解决降噪问题，其针对pMRI的应用与这样的情况相匹配，在这种情况下，下文中（例如SENSE重建）的表示为y的噪声参考图像可用于估算框架表示的超参数。在示出算法4后将讨论第二个问题且第二个问题对应于估算直接源自pMRI背景下的缩小的FOV图像的框架分解的超参数。

根据涉及系数x∈R^K的框架表示(FR)，可将所观察的信号y=F^*x+n写成如下形式：

y=F^*x+n   (30)

其中，n是所观察的信号与其FR F^*x之间的误差。该误差通过强加属于以下闭凸集的部分而被建模：

其中，δ∈[0,∞[是某个误差界限，N可以是上的任何范数。

在信号/图像恢复问题中，n是破坏所测量的数据的加性噪声。以下进展将聚焦于由均匀噪声构成的有界观察误差的情况。通过采用概率处理法，y和x被假定为随机向量Y和X的实现。在这种背景下，目的在于通过考虑某种参数概率模型和估算相关的超参数而描述X|Y的概率分布的特征。估算问题比在非冗余基（例如小波基）上的分解的情况下更加困难，因为F^*不是双射的。

在贝叶斯框架中，为框架系数限定先验分布是必要的。例如，该先验可被选择以便促进表示的稀疏性。在下文中，f(x,θ)表示框架系数的概率密度函数，其取决于未知的超参数向量θ，且f(θ)为针对超参数向量θ的先验pdf。按照约束C_δ，n被假定为均匀地分布在球

上，且f(y|x)为位于闭型凸面球上的均匀的pdf:

超参数向量θ也被认为是随机向量Θ的实现，且给定Y的情况下（X,Θ）的条件pdf可被写为：

f(x,θ|y)∝f(y|x)f(x|θ)f(θ)  (31)

为了简洁，假设框架系数为独立于边际GG（广义高斯）分布的先验，即使如上文所述，GGL分布更适合应用于pMRI。这导致以下框架系数先验：

$f\left(x_k\right|{\mathrm{\&alpha;}}_k,{\mathrm{\&beta;}}_k)=\frac{{\mathrm{\&beta;}}_k}{2{\mathrm{\&alpha;}}_k\mathrm{\&Gamma;}\left(\frac{1}{{\mathrm{\&beta;}}_k}\right)}\exp (-\frac{{\left|x_k\right|}^{{\mathrm{\&beta;}}_k}}{{\mathrm{\&alpha;}}_k^{{\mathrm{\&beta;}}_k}})---\left(32\right)$

其中，α_k>0,β_k>0（k∈{1,...,K}）为与x_k相关的尺度参数和形状参数，x_k是框架系数向量的分量，且Γ(·)是伽马函数。

通过引入

在框架先验可被重写为：

$f\left(x_k\right|{\mathrm{\&gamma;}}_k,{\mathrm{\&beta;}}_k)=\frac{{\mathrm{\&beta;}}_k}{2{\mathrm{\&gamma;}}_k^{1/{\mathrm{\&beta;}}_k}\mathrm{\&Gamma;}\left(\frac{1}{{\mathrm{\&beta;}}_k}\right)}\exp (-\frac{{\left|x_k\right|}^{{\mathrm{\&beta;}}_k}}{|{\mathrm{\&gamma;}}_k})---\left(33\right)$

框架系数的分布通常从一系数到另一系数是不同的。然而，一些框架系数可具有非常相似的分布，其可由相同的超参数来限定。结果，提出了将框架系数分成G个不同的组。第g组将通过表示为

的唯一的超参数向量而参数化。在该情况下，框架先验可被表示为：

$f\left(x\right|\mathrm{\&theta;})$ $=\underset{g=1}{\overset{G}{\mathrm{\&Pi;}}}\left[{\left(\frac{{\mathrm{\&beta;}}_g}{2{\mathrm{\&gamma;}}_g^{1/{\mathrm{\&beta;}}_g}\mathrm{\&Gamma;}\left(\frac{1}{{\mathrm{\&beta;}}_g}\right)}\right)}^{n_g}\exp (-\frac{1}{{\mathrm{\&gamma;}}_g}\underset{k\&Element;S_g}{\mathrm{\&Sigma;}}{\left|x_k\right|}^{{\mathrm{\&beta;}}_g})\right]---\left(34\right)$

其中，求和项覆盖了包含n_g个元素和θ_g=(θ₁...θ_G)的第g组的元素的指数集S_g。例如，每个组可对应于给定的小波子带。当使用多尺度的框架表示时，考虑到在给定分辨率水平下的所有框架系数都属于同一个组，较粗的分类可能是不明智的。

用于框架分解的分级贝叶斯模型可通过以下不适当的超先验来完成：

其中，如果ξ∈A，则1_A(ξ)是通过1_A(ξ)=1在

上所限定的函数，否则1_A(ξ)=0。

使用这种先验的动机如下：

●区间[0,3]覆盖了实际应用中所有可能遇到的β_g的值，且不存在关于参数的额外信息。

●用于参数γ_g的先验是杰弗里（Jeffrey）分布，其反映了关于该参数的知识的缺乏。

因此，所形成的后验分布由下式给出：

与后验分布（35）相关的贝叶斯估算量[例如，最大后验概率（MAP）或也被称为后验平均估算量的最小均方差(MMSE)估算量]不具有简单的闭型表达式。

在下文中，随机方法被提出用以计算MMSE估算量，其依赖于混合马尔可夫链蒙特卡尔(MCMC)算法。该思想在于使用合适的算法来生成根据后验分布（35）而分布的样本。在收敛后，所生成的样本被分别用于计算未知的模型参数的MMSE的估算量和超参数向量x和θ。在pMRI环境下，x是参考图像的搜索框架表示。

当所考虑的框架为标准正交M基的联合且N(·)是欧几里得范数时，可使用众所周知的吉布斯采样法（GS）[Geman,1984],其迭代地生成根据与目标分布相关的条件分布而分布的样本。更确切地说，基本的GS迭代地生成根据f(x|θ,y)和f(θ|x,y)而分布的样本以模拟全部后验的实现f(x,θ|y)。

直接的计算产生以下条件分布：

$f\left(x\right|\mathrm{\&theta;},y)\&Proportional;1_{C_{\mathrm{\&delta;}}}\left(x\right)\underset{g=1}{\overset{G}{\mathrm{\&Pi;}}}\exp (-\frac{1}{{\mathrm{\&gamma;}}_g}\underset{k\&Element;S_g}{\mathrm{\&Sigma;}}{\left|x_k\right|}^{{\mathrm{\&beta;}}_g})---\left(37\right)$

该条件分布是在C_δ上截断GG分布的产物。实际上，根据该截断分布的采样并不总是易于执行，因为伴随框架算子F_通常尺寸大。然而，在下文中详细介绍了两种可替换的采样策略。

通过根据以下独立GG分布的采样而进行简单采样：

$\underset{g=1}{\overset{G}{\mathrm{\&Pi;}}}\exp (-\frac{1}{{\mathrm{\&gamma;}}_g}\underset{k\&Element;S_g}{\mathrm{\&Sigma;}}{\left|x_k\right|}^{{\mathrm{\&beta;}}_g})$

然后，仅在

时，接受所提出的候选x。该方法可用于任何框架分解和任何范数。然而，由于非常低的接受率，尤其当δ采用小值时，该方法可能是极其低效的。

吉布斯采样法被设计成：当所考虑的框架是M个标准正交基的联合且N(·)是欧几里得范数时，从等式（36）中的条件分布中更高效地进行采样。在该情况下，分析框架算子和相应的伴随框架算子可被分别写成F=[F₁...F_M]^T和F^*=[F^* ₁...F^* _M]^T，其中，

F_m是第m次标准正交基上的分解算子，例如F^* _mF_m=F_mF^* _m=Id。在下文中，每个

（K=ML）被分解成

其中，

用于生成框架系数的GS根据约束条件N(y-F^*x)≤δ下的条件分布f(x_n|x_-n,y,θ)绘制向量，其中，x_-n是通过去除第n个向量x_n而由x构建的次元的减小的尺寸向量。如果N(·)是欧几里得范数，则

其中， $c_n=F_{n}y-\underset{m\&NotEqual;n}{\mathrm{\&Sigma;}}{F}_n{F}_m^{*}{x}_m.$

为了采集每个x_n，提出使用MH步骤，其所建议分布被支撑在球B_cn,δ上,B_cn,δ由下式限定：

现在将讨论来自基于B_cn,δ限定的pdf q_δ的随机生成。

首先，将考虑如何在

的单元l_p球（p∈]0,+∞]）中采样向量。在特殊情况下（p=+∞），这可通过沿着根据区间[-1,1]上的分布的每个空间坐标独立地进行采样而很容易地被执行。不应将这里所考虑的参数“p”与用在时间正则化中的形状参数混淆，参见等式21和等式22。

因此，该部分聚焦于与p的有限数值相关的更困难的问题。在下文中，||·||_p表示l_p范数。

使A=[A₁,...,A_L’]^T为i.i.d.分量的随机向量，i.i.d.分量具有以下GG(p^1/p,p)pdf：

使U=[U₁,...,U_L']^T=A/||A||_p。然后，可以示出（[Gupta等，1997]），随机向量U均匀地分布在

的l_p单位球的表面上，且U₁,…,U_L’-1的联合pdf为：

$f(u_1,...\&CenterDot;,u_{L^{\&prime;}-1})=\frac{p^{L^{\&prime;}-1}\mathrm{\&Gamma;}(L^{\&prime;}/p)}{2^{L^{\&prime;}-1}{\left(\mathrm{\&Gamma;}(1/p)\right)}^{L^{\&prime;}}}{(1-\underset{k=1}{\overset{L^{\&prime;}-1}{\mathrm{\&Sigma;}}}{\left|u_k\right|}^p)}^{(1-p)/p}{1}_{D_{p,L^{\&prime;}}}(u_1,...,u_{L^{\&prime;}-1})$

其中，

在

的单位l_p球上的均匀分布将由u(L’,p)表示。因此，U=[U₁,...,U_L']^T～u(L',p)。

对于每个∈{1,...,L’-1}，可以示出([Gupta等，1997])，V=[U₁,...,U_L]^T的pdf可由下式给出：

$q_1(u_1,...,u_L)\frac{p^L\mathrm{\&Gamma;}(L^{\&prime;}/p){(1-{\mathrm{\&Sigma;}}_{k=1}^L{\left|u_k\right|}^p)}^{(L^{\&prime;}-L)/p-1}}{2^L{\left(\mathrm{\&Gamma;}(1/p)\right)}^L\mathrm{\&Gamma;}((L^{\&prime;}-L)/p)}{l}_{D_{p,L+1}(u_1,...,u_L)}---\left(40\right)$

特别地，如果且L’=L+p，则等式（40）提供在R^L’的单位l_p球上的均匀分布。来自在半径满足δ≥η≥0的l_p球上的任何分布q_η的采样可通过缩放比例V而直接地被推导出。

具有该pdf的闭型表达式是重要的，以能够计算MH移动的接受率。考虑到在前次迭代（i-1）中所获得的x_n ^(i-1)的值，然而优选地，可以选择支撑在包含x_n ^(i-1)的半径η∈]0,δ[的受限球上的建议分布。该策略，与随机游动MH相似，导致与条件分布f(x|θ,y)的较大值相关的区域的更好的探索。更确切地说，提出了选择在

上所限定的建议分布，其中

${\hat{x}}_n^{(i-1)}=P(x_n^{(i-1)}-c_n)+c_n---\left(41\right)$

且P是到球B_0,δ-η上的投影，被限定成：

球的中心的这种选择保证了

此外，在

中连续绘制后可达到的任一点。半径η必须被调整以确保

的良好探索。实际上，也希望固定足够小的η值（与δ相比）以提高接受率。

取代根据f(θ|x,y)采样超参数向量θ，提出了根据f(γ_g|β_g,x,y)和f(β_g|γ_g,x,y)迭代地采样超参数向量θ。直接的计算产生下列结果：

$f\left({\mathrm{\&beta;}}_g\right|{\mathrm{\&gamma;}}_g,x,y)\&Proportional;\frac{{\mathrm{\&beta;}}_g^{n_g}{1}_{[0,3]}\left({\mathrm{\&beta;}}_g\right)}{{\mathrm{\&gamma;}}_g^{n_g/{\mathrm{\&beta;}}_g}{\left[\mathrm{\&Gamma;}(1/{\mathrm{\&beta;}}_g)\right]}^{n_g}}\exp \left(\underset{k\&Element;S_g}{\mathrm{\&Sigma;}}\frac{-{\left|x_k\right|}^{{\mathrm{\&beta;}}_g}}{{\mathrm{\&gamma;}}_g}\right)---\left(44\right)$

因此，f(γ_g|β_g,x,y)是易于采样的逆伽马分布的pdf。

相反地，根据截断的pdff(β_g|γ_g,x,y)采样是更加困难的。这通过使用MH移动而实现，MH移动的建议

是在区间[0,3]上被截断的高斯分布，其标准偏差为

[Dobigeon等，2007]、[Robert,1995]。该分布的模式是在前次迭代（i-1）中的参数

的值。所形成的方法是算法4中所概括的混合GS。

如果想要在同一贝叶斯框架中将重建步骤和超参数估算步骤整合在一起，则可将等式（30）中的观察模型扩展至y=SF^*x+n，其中，S为在等式（3）中所限定的灵敏度矩阵。因此，必须相应地重新进行推论。

算法4

1：利用一些

和x⁽⁰⁾∈C_δ进行初始化，且设定i=1。

2：重复。

3：采样x：

4：对于n=1到M，进行

5：计算 $c_n^{\left(i\right)}=F_n(y-{\mathrm{\&Sigma;}}_{m<n}{F}_m^{*}{x}_m^{\left(i\right)}-{\mathrm{\&Sigma;}}_{m>n}{F}_m^{*}{x}_m^{(i-1)})$ 和 ${\hat{x}}_n^{(i-1)}=P(x_n^{(i-1)}-c_n^{\left(i\right)})+c_n^{\left(i\right)}.$

6：如下模拟

·生成

其中在B_0,η上限定q_η，

·计算比例

$r({\tilde{x}}_n^{\left(i\right)},x_n^{(i-1)})=\frac{f\left({\tilde{x}}_n^{\left(i\right)}\right|{\mathrm{\&theta;}}^{(i-1)},{\left(x_m^{\left(i\right)}\right)}_{m<n},{\left(x_m^{(i-1)}\right)}_{m>n},y)q_{\mathrm{\&eta;}}(x_m^{(i-1)}-P({\tilde{x}}_n^{\left(i\right)}-c_n^{\left(i\right)})-c_n^{\left(i\right)})}{f\left(x_n^{(i-1)}\right|{\mathrm{\&theta;}}^{(i-1)},{\left(x_m^{\left(i\right)}\right)}_{m<n},{\left(x_m^{(i-1)}\right)}_{m>n},y)q_{\mathrm{\&eta;}}({\tilde{x}}_n^{\left(i\right)}-{\tilde{x}}_n^{(i-1)})}$

且根据概率

接受所提出的候选。

7：比对

8：采样θ：

9：对于g=1到G，进行

10：生成

11：如下模拟

·生成 ${\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}~q\left({\mathrm{\&beta;}}_g\right|{\mathrm{\&beta;}}_g^{(i-1)}),$

·计算比例

$r({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)},{\mathrm{\&beta;}}_g^{(i-1)})=\frac{f\left({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}\right|{\mathrm{\&gamma;}}_g^{\left(i\right)},x^{\left(i\right)},y)q\left({\mathrm{\&beta;}}_g^{(i-1)}\right|{\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)})}{f\left({\mathrm{\&beta;}}_g^{(i-1)}\right|{\mathrm{\&gamma;}}_g^{\left(i\right)},x^{\left(i\right)},y)q\left({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}\right|{\mathrm{\&beta;}}_g^{(i-1)})}$

根据概率

接受所提出的候选。

12：比对

13：设定i←i+1。

14：直到收敛。

尽管该算法是直观的且易于实现，但必须指出，该算法是在所考虑的框架是M个标准正交基的联合的限制性假设下被推导出的。当这些假设不成立时，如下文所讨论的另一算法允许通过利用框架的代数性质来采样框架系数和有关的超参数。实际上，由于根据f(x|θ,y)的样本的直接生成通常是不可能的，因此提出了可替代的方法，该方法用MH移动替换Gibbs移动。该MH移动目的在于根据建议分布整体采样候选x。利用标准的MH接受率，该候选被接受或被拒绝。MH移动的效率在很大程度上取决于对于x的建议分布的选择。

将x⁽ⁱ⁾表示为算法的第i次接受的样本，且建议q(x|x^(i-1))表示为用于产生在迭代i时的候选。对于选择q(x|x^(i-1))的主要困难源于以下事实：必须保证x∈C_δ，同时产生已处理的表达式q(x^(i-1)|x)/q(x|x^(i-1))。

出于该原因，提出了利用框架表示的代数性质。更确切地说，任何框架系数向量可被分解成x=x_H+x_H⊥，其中x_H和x_H⊥为分别在H=Ran(F)和H^⊥=[Ran(F)]^⊥=Null(F*)中取值的随机向量的实现。检索到F的范围为且F*的零空间为Null(F^*)={x∈R^K|F^*x=0}。

这里所使用的建议分布允许生成样本x_H∈H和x_H⊥,∈H^⊥。更确切地说，建议pdf的以下的可分离的形式被认为：

$q\left(x\right|x^{\left(i\right)})=q\left(x_H\right|x_H^{(i-1)})q\left(x_{H^{\&perp;}}\right|x_{H^{\&perp;}}^{(i-1)})---\left(45\right)$

其中， $x_H^{(i-1)}\&Element;H,$ $x_{H^{\&perp;}}^{(i-1)}\&Element;H^{\&perp;}$ 且 $x^{(i-1)}=x_H^{(i-1)}+x_{H^{\&perp;}}^{(i-1)}.$ 换句话说，将执行x_H和x_H⊥的独立采样。

如果考虑分解x=x_H+x_H⊥，则在C_δ中采样x相当于采样

其中，

实际上，x_H=Fλ，其中，λ∈R^L，并且，由于x_H⊥∈Null(F^*)，因此F^*x=F^*Fλ。可以很容易地实现在

中采样λ，例如，通过从球B_y,δ上的分布产生u且通过使λ=(F*F)^-1u。

为了使在迭代i时的x_H的采样更为有效，考虑到在前次迭代时的采样值 $x_H^{(i=1)}={\mathrm{F\&lambda;}}^{(i-1)}=F{\left(F^{*}F\right)}^{-1}{u}^{(i-1)}$ 可能是引人注意的。

与随机游动生成技术相似，以

生成u，其中η∈]0,δ[且

这允许绘制向量u，使得x_H=F(F^*F)^-1u∈C_δ且N(u-u^(i-1))≤2η。

假定N(·)是l_p范数且p∈[1,+∞]，则可根据等式（40）执行u的生成。

一旦已经模拟x_H=Fλ∈H∩C_δ(这确保了x在C_δ内)，则x_H⊥必须被采样作为H^⊥的元素。由于y=F^*x+n=F^*x_H+n，因此没有关于x_H⊥的y的信息。因此，提出了通过根据高斯分布绘制z且通过将z投射到H^⊥上而采集x_H，即 $x_{H^{\&perp;}}={\mathrm{\&Pi;}}_{H^{\&perp;Z}},$ ${\mathrm{\&Pi;}}_{H^{\&perp;}}=I-F{\left(F^{*}F\right)}^{-1}{F}^{*}$ 是到H^⊥上的正交投影算子。

因此，可以示出建议pdf的表达式：

$\frac{q\left(x^{(i-1)}\right|x)}{q\left(x\right|x^{(i-1)})}=\frac{q_{\mathrm{\&eta;}}(u^{(i-1)}-P(u-y)-y)}{q_{\mathrm{\&eta;}}(u-{\hat{u}}^{(i-1)})}---\left(46\right)$

在K=L（产生x_H⊥=0）的退化情况下，该表达式保持有效。最终，重要的是需要注意，如果q_η是球B_0,η上的均匀分布，则上述比率减小至1，这简化了MH接受率的计算。最终的算法被概括在算法5中。值得注意的是，如同对于算法4的混合GS一样，执行超参数向量的采样。

算法5

1：利用一些

和u⁽⁰⁾∈B_y,δ进行初始化，且设定x⁽⁰⁾=F(F^*F)^-1u⁽⁰⁾且i=1。

2：重复。

3：采样x：

·计算 ${\hat{u}}^{(i-1)}=P(u^{(i-1)}-y)+y.$

·生成

其中q_η被限定B_0,η上。

·计算 ${\tilde{x}}_H^{\left(i\right)}=F{\left(F^{*}F\right)}^{-1}{\tilde{u}}^{\left(i\right)}.$

·生成 $z^{\left(i\right)}~N(x^{(i-1)},{\mathrm{\&sigma;}}_x^{2}I).$

·计算 ${\tilde{x}}_{H^{\&perp;}}^{\left(i\right)}={\mathrm{II}}_{H^{\&perp;}}{z}^{\left(i\right)}$ 和 ${\tilde{x}}^{\left(i\right)}={\tilde{x}}_H^{\left(i\right)}+{\tilde{x}}_{H^{\&perp;}}^{\left(i\right)}.$

·计算比率 $f({\tilde{x}}^{\left(i\right)},x^{(i-1)})=\frac{f\left({\tilde{x}}^{\left(i\right)}\right|{\mathrm{\&theta;}}^{(i-1)},y)q_{\mathrm{\&eta;}}(u^{(i-1)}-P({\tilde{u}}^{\left(i\right)}-y)-y)}{f\left(x^{(i-1)}\right|{\mathrm{\&theta;}}^{(i-1)},y)q_{\mathrm{\&eta;}}({\tilde{u}}^{\left(i\right)}-{\hat{u}}^{(i-1)})}$ 且根据概率

接受建议的候选

4：采样θ：

5：对于g=1到G，进行

6：生成 ${\mathrm{\&gamma;}}_g^{\left(i\right)}~\mathrm{Ig}(\frac{n_g}{{\mathrm{\&beta;}}_g^{(i-1)}},{\mathrm{\&Sigma;}}_{k\&Element;S_g}{\left|x_k^{\left(i\right)}\right|}^{{\mathrm{\&beta;}}_g^{(i-1)}}).$

7：如下模拟

·生成 ${\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}~q\left({\mathrm{\&beta;}}_g\right|{\mathrm{\&beta;}}_g^{(i-1)}),$

·计算比例 $r({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)},{\mathrm{\&beta;}}_g^{(i-1)})=\frac{f\left({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}\right|{\mathrm{\&gamma;}}_g^{\left(i\right)},x^{\left(i\right)},y)q\left({\mathrm{\&beta;}}_g^{(i-1)}\right|{\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)})}{f\left({\mathrm{\&beta;}}_g^{(i-1)}\right|{\mathrm{\&gamma;}}_g^{\left(i\right)},x^{\left(i\right)},y)q\left({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}\right|{\mathrm{\&beta;}}_g^{(i-1)})}$ 且根据概率

接受建议的候选。

8：比对

9：设定i←i+1。

10：直到收敛。

等式（31-35）的分级贝叶斯模型可被扩展成包括取决于待重建图像的全变差（TV）的先验中的附加项。如同冗余框架表示，使用TV先验还导致超参数估算问题。

假定指数形状，则新的先验可被表达为：

$f\left(x\right|\mathrm{\&theta;})=\frac{1}{Z\left(\mathrm{\&theta;}\right)}\exp (-\mathrm{\&kappa;}{\left|\right|F^{*}x\left|\right|}_{\mathrm{TV}})\underset{g=1}{\overset{G}{\mathrm{\&Pi;}}}\left[{\left(\frac{1}{{\mathrm{\&gamma;}}_g^{1/{\mathrm{\&beta;}}_g}}\right)}^{n_g}\exp (-\frac{1}{{\mathrm{\&gamma;}}_g}\underset{k\&Element;S_g}{\mathrm{\&Sigma;}}{\left|x_k\right|}^{{\mathrm{\&beta;}}_g})\right]---\left(47\right)$

其中，θ=(θ₁,...,θ_G,k)是新的超参数向量，其中κ>0，||·||_TV是TV半范数以及Z(θ)是标准化常数。用于框架分解的新的分级贝叶斯模型可通过以下不适当的超先验来完成：

其中，κ_max是待确定的正实数（在优选实施方式中，κ_max=10）。

可以示出，Z(θ)相对于γ_g是均匀界定的，因而当γ_g→+∞时，等式（47）中的超先验f(θ)具有稳定的渐进特性。

因此，所形成的新的后验分布由下式给出：

与等式（48）中的后验分布相关的贝叶斯估算量不具有简单的闭型表达式。出于该原因，将应用上文所讨论的相同的采样策略，且将如在算法5中一样采样框架系数。然而，对于超参数向量，直接的计算示出针对超参数γ_g、β_g和κ的后验分布将被分别表达为：

$f\left({\mathrm{\&beta;}}_g\right|{\mathrm{\&gamma;}}_g,\mathrm{\&kappa;},x,y)\&Proportional;\exp (\underset{k\&Element;S_g}{\mathrm{\&Sigma;}}-\frac{1}{{\mathrm{\&gamma;}}_g}{\left|x_k\right|}^{{\mathrm{\&beta;}}_g})1_{\left[\mathrm{0,3}\right]}\left({\mathrm{\&beta;}}_g\right)---\left(51\right)$

以及

$f\left(\mathrm{\&kappa;}\right|{\mathrm{\&gamma;}}_1,...,{\mathrm{\&gamma;}}_G,{\mathrm{\&beta;}}_1,...,{\mathrm{\&beta;}}_G,x,y)\&Proportional;\exp (-\mathrm{\&kappa;}{\left|\right|F^{*}x\left|\right|}_{\mathrm{TV}})1_{[0,{\mathrm{\&kappa;}}_{\max }]}\left(\mathrm{\&kappa;}\right)---\left(52\right)$

因此，f(γ_g|β_g,k,x,y)是逆伽马

的pdf。因此，将如在算法4中一样精确地采样γ_g。相反地，根据f(β_g|γ_g,k,x,y)和f(k|γ₁,...,γ_G,β₁,...,β_G,x,y)的采样更加困难。通过使用两次MH移动而实现该任务，两次MH移动的建议分布和q(k|k^(i-1))是在区间[0,3]和[0,κ_max]上被截断的高斯分布，其标准偏差分别为σ_βg=0.05和σ_κ=0.01。已基于实际观察结果确定了这些标准偏差值。所形成的根据等式（49）中的后验分布的采样方法被概括在算法6中。

算法6

1：利用一些 ${\mathrm{\&theta;}}^{\left(0\right)}=({\left({\mathrm{\&theta;}}_g^{\left(0\right)}\right)}_{1\&le;g\&le;G},k^{\left(0\right)})=({({\mathrm{\&gamma;}}_g^{\left(0\right)},{\mathrm{\&beta;}}_g^{\left(0\right)})}_{1\&le;g\&le;G},k^{\left(0\right)})$ 和u⁽⁰⁾∈B_y,δ进行初始化，且

设定x⁽⁰⁾=F(F^*F)^-1u⁽⁰⁾且i=1。

2：重复。

3：采样x：

·计算 ${\hat{u}}^{(i-1)}=P(u^{(i-1)}-y)+y.$

·生成 ${\tilde{u}}^{\left(i\right)}~q_{\mathrm{\&eta;}}(u-{\hat{u}}^{(i-1)}),$ 其中q_η被限定B_0,η上。

·计算 ${\tilde{x}}_H^{\left(i\right)}=F{\left(F^{*}F\right)}^{-1}{\tilde{u}}^{\left(i\right)}.$

·生成 $z^{\left(i\right)}~N(x^{(i-1)},{\mathrm{\&sigma;}}_x^{2}I).$

·计算 ${\tilde{x}}_{H^{\&perp;}}^{\left(i\right)}={\mathrm{\&Pi;}}_{H^{\&perp;}}{z}^{\left(i\right)}$ 和 ${\tilde{x}}^{\left(i\right)}={\tilde{x}}_H^{\left(i\right)}+{\tilde{x}}_{H^{\&perp;}}^{\left(i\right)}.$

·计算比率 $r({\tilde{x}}^{\left(i\right)},x^{(i-1)})=\frac{f\left({\tilde{x}}^{\left(i\right)}\right|{\mathrm{\&theta;}}^{(i-1)},y)q_{\mathrm{\&eta;}}(u^{(i-1)}-P({\tilde{u}}^{\left(i\right)}-y)-y)}{f\left(x^{(i-1)}\right|{\mathrm{\&theta;}}^{(i-1)},y)q_{\mathrm{\&eta;}}({\tilde{u}}^{\left(i\right)}-{\hat{u}}^{(i-1)})}$ 且根据概率

接受建议的候选

4：采样θ：

5：对于g=1到G，进行

6：生成 ${\mathrm{\&gamma;}}_g^{\left(i\right)}~\mathrm{IG}(\frac{n_g}{{\mathrm{\&beta;}}_g^{(i-1)}},{\mathrm{\&Sigma;}}_{k\&Element;\mathrm{Sg}}{\left|x_k^i\right|}^{{\mathrm{\&beta;}}_g^{(i-1)}}).$

7：如下模拟

·生成 ${\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}~q\left({\mathrm{\&beta;}}_g\right|{\mathrm{\&beta;}}_g^{(i-1)}),$

·计算比例

$r({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)},{\mathrm{\&beta;}}_g^{(i-1)})=\frac{f\left({\widetilde{\mathrm{\&beta;}}}_f^{\left(i\right)}\right|{\mathrm{\&gamma;}}_g^{\left(i\right)},{\mathrm{\&kappa;}}^{\left(i\right)},x^{\left(i\right)},y)q\left({\mathrm{\&beta;}}_g^{(i-1)}\text{|}{\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}\right)}{f\left({\mathrm{\&beta;}}_g^{(i-1)}\right|{\mathrm{\&gamma;}}_g^{\left(i\right)},{\mathrm{\&kappa;}}^{\left(i\right)},x^{\left(i\right)},y)q\left({\widetilde{\mathrm{\&beta;}}}_g^{\left(i\right)}\right|{\mathrm{\&beta;}}_g^{(i-1)})}$ 且根据概率接受建议的候选。

8：比对

9：如下模拟k⁽ⁱ⁾:

·生成 ${\tilde{k}}^{\left(i\right)}~q\left(k\right|k^{(i-1)}),$

·计算比例 $r({\widetilde{\mathrm{\&kappa;}}}^{\left(i\right)},{\mathrm{\&kappa;}}^{(i-1)})=\frac{f\left({\widetilde{\mathrm{\&kappa;}}}^{\left(i\right)}\right|{\mathrm{\&gamma;}}_1^{\left(i\right)},...,{\mathrm{\&gamma;}}_G^{\left(i\right)},{\mathrm{\&beta;}}_1^{\left(i\right)},...,{\mathrm{\&beta;}}_G^{\left(i\right)},x^{\left(i\right)},y)q\left({\mathrm{\&kappa;}}^{(i-1)}\right|{\widetilde{\mathrm{\&kappa;}}}^{\left(i\right)})}{f\left({\mathrm{\&kappa;}}^{(i-1)}\right|{\mathrm{\&gamma;}}_1^{\left(i\right)},...,{\mathrm{\&gamma;}}_G^{\left(i\right)},{\mathrm{\&beta;}}_1^{\left(i\right)},...,{\mathrm{\&beta;}}_G^{\left(i\right)},x^{\left(i\right)},y)q\left({\widetilde{\mathrm{\&kappa;}}}^{\left(i\right)}\right|{\mathrm{\&kappa;}}^{(i-1)})}$ 且根据概率

接受建议的候选。

10：设定i←i+1。

11：直到收敛。

下面将参照图2至图20讨论本发明方法的不同实施方式的技术结果。

在下列条件下获得图2至图10以及图13中所报道的解剖结果。对包括利用0.93x0.93x8mm³的空间分辨率解剖的256x256x14的梯度回波（GE）的实数集进行了实验。GE解剖图像在TE/TR=10/500ms且BW=31.25kHz下进行获取。还注意到，在具有八通道头部线圈的±1.5特斯拉（Tesla）的GE保健扫描器上使用加速度因子R=2和R=4获取这些图像。有趣的是，在非并行成像中，解剖数据的扫描时间持续5分钟，然而，在R=2和R=4的并行成像中，获取持续时间分别减少至3分钟10秒和2分钟20秒。

图2和图3示出使用具有Tikhonov正则化的SENSE所获得的人脑的9张解剖切片。Tikhonov参数已经被手动固定。对于图2，缩减因数R=2，对于图3，缩减因数R=4。尽管正则化，但仍可以看到一些曲线形式的混叠伪像，特别是在图3上。还可看到图像的在一定程度上的过度平滑。

在相同的实验条件下，使用具有手动调谐的正则化参数κ的TV正则化（左面板：R=2；右面板：R=4，示出单切片）获得了图4。可看到极强的混叠伪像和一些“楼梯”形缺陷；这些缺陷可通过增加κ值而减轻，但以图像的信息内容为代价。

在相同的实验条件下，使用如上所述的基于自动校准的2D小波变换的正则化方案获得了图5和图6；该方法被被称作2D-UWR-SENSE，其中“UWR”代表“无约束的小波正则化”。此外，对于图5，R=2，对于图6，R=4。

更确切地说，与长度为8的滤波器相关的二元（M=2）对称标准正交小波基被用在j_max=3的分辨率等级下。至于小波系数，已经利用等式（10）的先验。

使用上文描述的贝叶斯方法估算有关的超参数（一对超参数适合于每个子带的实部/虚部，即在每个分辨率等级和定向下的每个近似系数和细节系数）。因而，使用基于小波的正则化进行全部FOV图像重建。

在利用Tikhonov正则化的图2和图3中所观察到的平滑效应不再存在于图5和图6的WT正则化图像中，其中，在脑遮蔽范围内进行非常精确的重建，而不会引入利用TV正则化所观察到的楼梯效应（参见图4）。

可以通过将额外约束并入上述方法中而获得图像质量的进一步改进，以便更好地正则化伪像区域。因而产生的算法被称作CWR-SENSE，其中，“CWR”代表“受约束的小波正则化”。

该方法暗示了在伪像区域中将局部的下界和上界施加到图像强度值，而不考虑伪像区域的形状和/或位置。这些边界限定非空的闭凸集：

其中，利用下式来建模在位置r∈2{1,…,Y/R}×{1,…,X}处在范围值上所引入的约束：

其中，

且那么，最优性准则变成：

其中，

且

是闭凸集C^*的指标函数，由下式限定：

形态学梯度[Serra,1982]可被用于检测伪像区域（在梯度图像中的极低/极高的转变），在伪像区域上应用额外的凸约束。由于限定凸集C_r的上界和下界依赖于r，因此它们在空间上是不同的。因此，可使用应用到基本SENSE重建图像（辅助图像或参考图像）的形态学开闭运算来计算这些上界和下界，以便删除极低和极高的强度。

图7和图8示出了分别对于R=2和R=4的CWR-SENSE方法的结果。可以看出，图5和图6中存在的伪像现已由于使用额外的凸约束的各向异性平滑而被去除。

从定量角度看，与基本SENSE重建和Tikhonov重建相比较，通过UWR/CWR-SENSE算法实现了显著的改进。下面的表1针对图3、图6和图8（R=4）中所示的解剖的脑体积的示出切片，报道了对应于基本SENSE、Tikhonov正则化和所提出的UWR/CWR-SENSE技术的以dB表示的信号-噪声比（SNR）。通常，从所提出的受约束的正则化策略所得到的增益达到了1.08dB和较好的视觉质量。

表1

	SENSE	Tikhonov	UWR-SENSE	CWR-SENSE
					切片#1	14.36	14.63	14.77	14.95
切片#2	11.55	11.62	12.01	12.53
					切片#3	12.95	13.44	14.02	14.22
切片#4	9.24	9.48	10.01	10.30
					切片#5	11.50	11.79	12.06	12.25
切片#6	9.68	9.87	10.12	10.32
					切片#7	11.00	11.56	11.85	12.00
切片#8	12.16	12.49	13.00	13.36
					切片#9	13.78	14.77	15.83	16.04
平均体积	11.80	12.18	12.63	12.88

还研究了选择小波基的影响。更具体地，考虑了四个不同基：二元Symmlet8、二元多贝西（Daubechies）8、二元Haar和具有M=4带的Meyer[Chaux等，2006b]。前三个基产生了非常相似的结果（除了一些仅利用Haar基才发生的方块效应外），二元Symmlet8导致略微较高的SNR。相反地，Meyer的4带小波基导致显著较低的SNR。

图9示出了通过在j_max=3分辨率级别下，使用与长度为8的滤波器相关的二元（M=2）Symmlet标准正交小波基的组合小波-全变差正则化（算法2）而获得的结果，其中，缩减因数R=4。图10示出了通过使用具有Symmlet8滤波器和Symmlet4滤波器的2个标准正交基的联合（U2OB）所构成的冗余小波框架的组合小波-全变差正则化而获得的结果。使用上文所描述的方法估算小波先验的超参数和TV正则化参数。这些图示出，所重建的图像比使用UWR-SENSE算法所重建的图像（图6）具有更好的正则性。然而，从视觉角度看，使用冗余WT未必会导致更好的重建质量。从定量角度看，在下文表2中提供了以dB表示的SNR值。与表1中的SNR值的比较证实了在联合正则化框架中将WT和TV组合的有效性。该表示出，当使用Symmlet小波时，与UWR-SENSE算法相比，获得了0.02dB的微小改进。然而，当使用U2OB冗余小波框架时，获得了更加显著的改进（0.24dB）。于是证明了，即使使用标准正交小波基和冗余小波框架所获得的重建性能似乎是相同的，但是SNR值指示出使用冗余框架效果更好。

图13示出了使用TV正则化（左）、具有由两个标准正交小波基的联合所构成的小波框架的3D-UWR-SENSE（中）以及混合框架-TV正则化方法（右）所重建的图像。顶行示出了整个图像，底行示出了放大的细节。两个标准正交基的联合构成了用于重建的小波框架（中列和右列），这两个标准正交基是长度为4的Daubechies基和长度为8的移位的Daubechies基。使用了三种分辨率级别，这意味着考虑了G=20组的小波系数。为了表示框架，使用上文所述的混合MCMC算法选择超参数。

表2

	OB	U2OB
			切片#1	14.85	15.16
切片#2	12.18	12.35
			切片#3	13.85	14.32
切片#4	9.78	9.85
			切片#5	12.17	12.62
切片#6	10.20	10.17
			切片#7	11.75	12.04
切片#8	13.15	13.44
			切片#9	15.94	15.88
平均体积	12.65	12.87

在下列条件下获得了在图11和图12中所提到的解剖结果。使用3D T1-加权的MP-RAGE脉冲序列和由32个接收通道构成的矩阵阵列头部线圈，在3TTim Trio西门子扫描仪上进行解剖的MRI扫描。选择扫描参数如下：切片定向=矢状的（右>左)、切片厚度=1.1mm；TE=2.98ms、TR=2300ms；TI=900ms、翻转角：9°、BW=61kHz，单次激发。视野为256x240x176mm³，256x240x160的矩阵尺寸对应于1x1x1.1mm³的各向异性的分辨率。使用传统的MRI的扫描时间为TA_conv=9分14秒，即不需要加速（既不是pMRI也不是部分傅里叶）。使用6/8部分傅里叶方案能够使扫描时间减少至7分46秒。收集并行MRI数据，而不需使用加速度因数R=2或R=4的部分傅里叶。这些测试的各自的扫描时间是5分钟3秒和2分钟59秒，而不是所期望的TA_conv/R。其原因在于西门子k-空间采样策略。制造商对24条中央线采用全采样方案，从而导致实际的加速度因数小于规定的加速度因数（1.83而不是2以及3.09而不是4）。

图11至图12示出了来自使用矩阵阵列32通道的头部线圈在3Tesla（西门子Tim Trio）处所获取的单个对象T1-加权的MRI数据的不同的MRI重建算法。图11和图12分别示出了对于同一对象的轴向视图和冠状视图。在这些图中的顶行示出了地面实况，即根据全k-空间获取所执行的重建。因此，在噪声背景下（即R=4），比较mSENSE（中间行）算法和3D-UWR-SENSE（底行）算法。对于后一算法，已针对J=2分辨率级别和Daubechies小波（非冗余）获得了结果。此外，如在上文所述，已使用完整数据最大似然法估算小波表示的超参数。在灰质边界上呈现为白质斑点的MSENSE重建伪像由两个视图轮廓中的白色环所指示。很明显，因为这些伪像没有出现在相同位置处，因此3D-UWR-SENSE算法胜过mSENSE技术。此外，通过本发明的基于小波的算法，在很大程度上过滤了在轴向视图的中央内的椭圆形伪像。冠状视图证实了使用所提出的发明在脑干和皮质下的区域中的实质性的降噪。

这里未示出使用2D小波基所获得的结果，因为该结果在质量上接近于图11和图12中所报道的结果。然而，就信噪比（SNR）来说，基于3D小波的重建在数量上胜过基于2D小波的重建。注意，相对于2D技术，3D重建产生了1.3dB的SNR改进。

到这里为止，仅仅讨论了关于静态的解剖学成像的情况的结果。然而，如上文所述，本发明还适于例如通过回波平面成像（EPI）所获得的4D fMRI图像序列的时空正则化。图14至图20示出了相应的结果。

出于验证目的，在认知定位器[Pinel等，2007]拟定期间，使用梯度回波EPI（GE-EPI）序列（TE=30ms、TR=2400ms、切片厚度=3mm、横向定向、FOV=192mm²），在3T西门子Trio磁体上获取fMRI数据。该实验被设计成映射听觉的、视觉的和运动的脑功能以及较高的认知任务，例如数字处理和语言理解（听和读）。它由包括N_r=128次扫描的单个时段构成。该范例为包含六十个听觉的、视觉的和运动的刺激的与快速事件相关的设计，这些刺激被限定在十种实验条件下（听觉的和视觉的句子、听觉的和视觉的计算、左/右听觉的和视觉的点击、水平的和垂直的棋盘）。使用L=32通道的线圈以促使并行成像。

对于15个对象中的每个对象，使用不同的缩减因数（R=2或R=4），在2×2mm²空间平面分辨率下收集fMRI数据。基于扫描仪所传输的原始数据文件，使用以下两种具体的处理方式重建缩小的FOV EPI图像：

i）在读出渐变坡度期间，重绘k-空间网格，以解释非均匀的k-空间采样，这发生在快速MRI序列（如GE-EPI）中；

ii）在EPI图像的k-空间获取期间，校正奈奎斯特重影，以去除奇-偶回波的不一致性。

然后，在最后步骤中利用4D-UWR-SENSE算法和已经讨论的3D-UWR-SENSE（简称为“UWR-SENSE”)算法来重建全FOV EPI图像，并将这两种算法与mSENSE方案（mSENSE是利用西门子扫描仪所实现的未正则化的SENSE算法）进行比较。

图14比较了两种pMRI重建算法，以在轴向切片、冠状切片、矢状切片上示出如何使用4D-UWR-SENSE方法去除掉mSENSE重建伪像。mSENSE重建图像实际上呈现出位于脑的感官区域和认知区域（颞叶、额叶皮质和运动皮质）中的中央和边界处的大型伪像；在图上，伪像由椭圆画出轮廓。这导致SNR损失，从而可能对这些脑区域中的活化检测产生显著影响。

不考虑重建管线，于是使用SPM5软件（http://www.fil.ion.ucl.ac.uk/spm/software/spm5/）预处理全FOV fMRI图像：预处理所涉及的重新排列、运动的校正、切片获取时间上的差异、空间标准化、以及利用在半极大值时的4mm全宽的各向同性高斯核的平滑。通过具有利用32通道头部线圈所获取的解剖的T1扫描的功能性图像的配准而执行解剖标准化到MNI空间。用于标准化到MINI空间的参数通过将该扫描标准化成SPM5所提供的T1MINI模板来估计，并且随后将这些参数应用于全部的功能性图像。

为了进行对象级的分析，一般线性模型（GLM）被构造成捕获与刺激有关的BOLD响应。设计矩阵依赖于十种实验条件且因此由对应于粘贴函数的21个回归量构成，粘贴函数与规范的血流动力学响应函数（HRF）及其第一时间导数相卷积，最后的回归量构成基线。因此，该GLM适于相同的获取图像且使用西门子重建器、UWR-SENSE和4D-UWR-SENSE进行重建。

针对运动反应和较高认知功能（计算、语言）的对照估计图像进一步经历对象级和组级的分析。由于所期望的激活位于不同的脑区域中，且该两种对照可能被重建伪像不同地破坏，因此这些两种对照是互补的。

更确切地说，研究了：

·听觉计算与听觉语句（aC-aS）的对照，其被认为引起额叶和顶叶中的诱发活性，因为解决心算任务涉及工作存储器，更具体地涉及顶内沟；

·左点击与右点击（Lc-Rc）的对照，对于该对照，期望在右运动皮质（中央前回、额中回）中诱发的活性。事实上，Lc-Rc对照限定了涉及以视觉或听觉形态所呈现的两种运动刺激的混合比较。因此，该比较的目的在于检测运动皮质中的单侧化效应。

当查看感官区（视觉区、听觉区和运动区）或者涉及更高的认知功能（阅读、计算）的区域时，因为这两个对照概括了针对该范例而发生的完全不同的情形（大激活群与小激活群、分布的激活图案与聚集的激活图案、双向活性与单向活性），因此选择该两个对照。

图15示出了针对aC-aS对照的叠加在解剖的MRI上的对象级的student-t图片。分别使用mSENSE、UWR-SENSE和4D-UWR-SENSE重建数据，其中R=2(图的顶部)和R=4（图的底部）。采用神经学惯例（左是左）。交叉表示最大激活峰。

对于最显著的切片且R=2，全部pMRI重建算法在左顶叶皮质和额叶皮质中，更确切地说是在顶下小叶和额中回中，成功地发现了诱发的活性。然而，对于R=4，仅有UWR-SENSE和4D-UWR-SENSE，优选后者，能够检索到通过心算所引发的可靠的额活性，其通过mSENSE算法找不到。从定量角度看，所提出的4D-UWR-SENSE算法发现较大的群，其局部最大值比使用mSENSE和UWR-SENSE所获得的局部最大值更为显著，如表3中所示出的。关于对于R=2的最显著的群，不论什么重建算法，峰位置都保持稳定。然而，检验其显著性水平，可以首先在将UWR-SENSE结果与mSENSE结果进行比较时测量基于小波的正则化的优点，然后在查看4D-UWR-SENSE结果时测量时间正则化和3D小波分解的额外正效应。对于R=4，也可展示出这些效益。表3示出了针对aC-aS对照的对象级的显著统计结果（在α=0.05的显著性水平下对多重比较进行校正，这意味着如果p值小于或等于α，则拒绝虚假设)。

不应当将统计的显著性试验的“p值”与时间正则化（参见等式21和等式22）中所使用的“形状参数”相混淆，或与用于阐明MCMC算法的参数“p”相混淆。

表3

图16示出了在R=2时针对aC-aS对照的所检测的激活的对象间的差异性。事实上，当比较使用不同管线（R=2）所重建的对象级的student-t图片时，可以观察到，mSENSE算法不能检测到在第二对象的期望区域中的任何激活群。相反，4D-UWR-SENSE方法检索到更多相干的激活，但对于第一对象来讲并不恰好在相同的位置处。

图17示出了针对Lc-Rc对照的叠加在解剖的MRI上的对象级的student-t图片。分别使用mSENSE、UWR-SENSE和4D-UWR-SENSE重建数据。

可以看出，所有的重建方法都能够在右侧中央前回中检索到所期望的激活。然而，当更为仔细地查看统计结果（参见表4）时，UWR-SENSE算法且更优选的4D-UWR-SENSE算法在右侧额中回中检索到额外的群。关于R=4时所获取的数据，相同的Lc-Rc对照引发了类似的激活，即在相同的区域中。在该图的底部可以看出，当根据本发明的方法进行pMRI重建时，该活性增强。

在表4中的定量结果在数字上证实了可在图中观察到：对于R=2和R=4，都使用4D-UWR-SENSE算法来检测具有较高的局部t-值的较大集群。更确切地说，表4示出了针对Lc-Rc对照的处于对象级的显著统计结果（在α=0.05的显著水平下对多重比较进行校正，这意味着如果p值小于或等于α，则拒绝虚假设）。

表4

图18指出了对于该运动对照，提出的pMRI管线对于对象间差异的稳健性。由于感官功能被期望生成较大的BOLD效应（较高的SNR）且表现得更稳定，因此在R=4时进行本发明的比较。比较使用不同的pMRI算法所重建的两个对象级的student-t图片。对于第二对象，可以观察到mSENSE算法不能检测在右侧运动皮质中的任何激活集群。相反，对于该第二对象，4D-UWR-SENSE方法在所期望的区域中检索到更多相干的激活。

总之，关于这两个对照，在统计显著性（集群数目、集范围、峰值…）以及稳健性方面，4D-UWR-SENSE算法始终胜过其它重建方法。

由于对象间的解剖变异和功能变异，组级分析是必需的，以便在总体水平上导出稳健且可重复的结论。为了该验证，对先前研究的对象级的对照图片进行了涉及15个健康对象的随机效应分析（RFX）。更确切地说，使用SPM5对对象级的对照图像（例如Lc-Rc图像、aC-aS图像）进行样本的Student-t试验。

图19示出了针对aC-aS对照的组级的student-t图片，其中，使用mSENSE、UWRSENSE和4D-UWR-SENSE重建数据（R=2和R=4）。使用了神经学惯例。箭头指示全局最大激活峰值。

这些图片说明，不考虑重建方法，给定更好的SNR，则在利用R=2所获取的数据集上发现更大且更为显著的激活。其次，对于R=2，目测证实了仅有4D-UWR-SENSE算法允许在顶叶皮质（参见轴向MIP切片）和较大群中检索到显著的双向激活，以及针对横跨不同的重建器的稳定的集群检索到显著的增益。对于R=4，当查看图的底部时，可得到类似的结论。对于R=2和R=4，在表5中可得到互补的结果且这些结果在数值上证实了该视觉比较：

·不论使用何种重建方法，使用R=2时，统计性能是更为显著的，尤其在集群级下，因为集群范围降低了一个数量级。

·使用4D-UWR-SENSE方法，而不是mSENSE或UWR-SENSE，增强了体素级结果和集群级的结果。

表5-针对aC-aS对照的组级的显著统计结果（在p=0.05时对多重比较进行校正）。

图20报道了关于Lc-Rc对照的R=2和R=4的类似的组级的MIP结果。该图示出，无论使用何种加速度因数R，管线都能够在运动皮质中检测到更多的在空间上延伸的激活区域。当将使用4D-UWR-SENSE方法所检测到的集群与那些通过mSENSE而发现的集群相比较时（再次不考虑R），在表6中从数量上证实了该目测。最终，4D-UWR-SENSE算法胜过UWR-SENSE算法，这证实了所提出的时空正则化方案的优点。

表6-针对Lc-Rc对照的组级的显著统计结果（在p=0.05时对多重比较进行校正）。

参考文献

[Bertsekas,1995]Bertsekas,D.P.(1995).Nonlinear programming,Second Edition.Athena Scientific,Belmont，美国，特别是第159-165页。

[Block等，2007]Block,K.T.,Uecker,M.和Frahm,J.(2007).Undersampled radialMRI with multiple coils.Iterative image reconstruction using a total variationconstraint.Magnetic Resonance in Medicine,56(7):1086–1098。

等，2008]:L.

J.-C.Pesquet,A.Benazza-Benyahia,P.Ciuciu,Autocalibrated Parallel MRI Reconstruction in the Wavelet Domain,in:IEEEInternational Symposium on Biomedical Imaging,法国，巴黎，2008，第756-759页。

等，2009]:L.

J.-C.Pesquet,Ph.Ciuciu和A.Benazza-Benyahia.AnIterative Method for Parallel MRI SENSE-based Reconstruction in the WaveletDomain,.arXiv:0909.0368v1[math.OC]。

[Chaux等，2007]:Chaux,C.,Combettes,P.,Pesquet,J.-C.和Wajs,V.(2007).Avariational formulation for frame-based inverse problems.Inverse Problems,23(4):1495–1518。

[Combettes和Pesquet,2008]:Combettes,P.L.和Pesquet,J.-C.(2008).A proximaldecomposition method for solving convex variational inverse problems.InverseProblems,24(6):27。

[Daubechies等，2004]:I.Daubechies,M.Defrise,C.DeMol,An iterativethresholding algorithm for linear inverse problems with a sparsity constraint,Communications on Pure and Applied Mathematics57(11)(2004)1413–1457。

[Dempster,1997]:Dempster,A.,Laird,A.和Rubin,D.(1977).Maximum likelihoodfrom incomplete data via the EM algorithm(with discussion).Journal of the RoyalStatistical Society,Series B,39:1–38。

[Dobigeon等，2007]:Dobigeon,N.和Tourneret,J.-Y.(2007).Truncated multivariateGaussian distribution on a simplex.Technical report,University of Toulouse。

[Donoho,1995]:D.Donoho.De-noising by soft-thresholding.IEEE Transactions onInformation Theory,41(3):613–627,1995。

[Geman,1984]:S.Geman和D.Geman,“Stochastic relaxation,Gibbs distributionand the Bayesian restoration of image,”IEEE Trans.Pattern Anal.Mach.Intell.,第6卷，第721–741页，1984。

[Griswold等，2002]:M.A.Griswold,P.M.Jakob,R.M.Heidemann,M.Nittka,V.Jellus,J.Wang,B.Kiefer,A.Haase,Generalized autocalibrating partially parallelacquisitions GRAPPA,Magnetic Resonance in Medicine47(6)(2002)1202–1210。

[Gupta等，1997]:Gupta,A.K.和Song,D.(1997).Journal of Statistical Planningand Inference.第60卷，第241-260页。

[Hoge等，2005]:W.S.Hoge,D.H.Brooks,B.Madore,W.E.Kyriakos,A tour ofaccelerated parallel MR imaging from a linear systems perspective,Concepts inMagnetic Resonance27A(1)(2005)17–37。

[Pesquet-Popescu,Pesquet]Ondelettes et applications,Techniques de l’Ingénieur,traitéTélécoms,TE5215。

[Pinel等，2007]:P.Pinel,B.Thirion,S.Mériaux,A.Jobert,J.Serres,D.Le Bihan,J.-B.Polin和S.Dehaene,“Fast reproducible identification and large-scaledatabasing of individual functional cognitive networks”BMC Neurosci.,第8卷，

第1期，第91页，2007年10月。

[Pruessmann等，1999]:K.P.Pruessmann,M.Weiger,M.B.Scheidegger,P.Boesiger,SENSE:sensitivity encoding for fast MRI,Magnetic Resonance inMedicine42(5)(1999)952–962。

[Raj等，2007]:Raj,A.,Singh,G.,Zabih,R.,Kressler,B.,Wang,Y.,Schuff,N.和Weiner,M.(2007).Bayesian parallel imaging with edge-preserving priors.MagneticResonance in Medicine,57(1):8–21。

[Robert,1995]:“Simulation of truncated normal variables,”Statistics andComputing,第5卷，第121–125页，1995。

[Serra,1982]Serra,J.(1982).Image Analysis and Mathematical Morphology.Academic Press,伦敦。

[Sodickson等，1997]:Sodickson,D.K.和Manning,W.J.(1997).Simultaneousacquisition of spatial harmonics(SMASH):fast imaging with radiofrequency coilarrays.Magnetic Resonance in Medicine,38(4):591–603。

[Ying等，2004]:L.Ying,D.Xu,Z.-P.Liang,On Tikhonov Regularization for imagereconstruction in parallel MRI,IEEE Engineering in Medicine and Biology Society,San Francisco,美国，2004，第1056–1059页。

Claims (27)

1.一种身体的并行磁共振成像的方法，包括：

-从具有已知的或估算的灵敏度分布和噪声协方差矩阵的各个接收天线中获取所述身体的一组基础磁共振图像，所述基础图像在k-空间中欠采样；以及

-执行所述身体的磁共振图像的正则化重建；

其中，执行磁共振图像的正则化重建的所述步骤通过最小化成本函数而在离散的框架空间中进行，所述成本函数包括：

-误差项，在给定所获取的基础图像的情况下，所述误差项表示重建图像的似然；和

-框架惩罚项，所述框架惩罚项表示所述重建图像的框架系数的实际统计分布与所述系数的先验分布之间的偏差；

基于所述身体的辅助的磁共振图像估算所述重建图像的框架系数的先验分布。

2.根据权利要求1所述的方法，其中，在给定所述所获取的基础图像的情况下，所述误差项表示所述重建图像的负对数似然。

3.根据权利要求1或2所述的方法，其中，在给定所述所获取的基础磁共振图像和所述框架系数的先验分布的情况下，通过在所述框架空间中最大化限定所述身体的图像的一组框架系数的完全后验分布，而进行所述执行所述身体的磁共振图像的正则化重建的步骤。

4.根据前述任一项权利要求所述的方法，其中，根据所述所获取的基础磁共振图像重建所述身体的所述辅助的磁共振图像。

5.根据权利要求4所述的方法，其中，使用灵敏度编码-SENSE算法，根据所述所获取的基础磁共振图像重建所述身体的所述辅助的磁共振图像。

6.根据权利要求5所述的方法，其中，使用选自以下算法中的一种算法，根据所述所获取的基础磁共振图像重建所述身体的所述辅助的磁共振图像：

-未正则化的SENSE算法；

-在图像空间中正则化的SENSE算法；和

-在k-空间中正则化的SENSE算法。

7.根据前述任一项权利要求所述的方法，其中，假设所述框架系数的广义Gauss-Laplace先验统计分布，使用最大似然估算量或后验均值估算量，基于所述身体的所述辅助的磁共振图像估算所述分布的参数。

8.根据前述任一项权利要求所述的方法，其中，所述误差项为二次平均误差项。

9.根据前述任一项权利要求所述的方法，其中，所述所获取的基础图像为三维图像，且其中，在离散的三维框架空间中进行所述执行磁共振图像的正则化重建的步骤。

10.根据权利要求9所述的方法，其中，通过叠加待成像的对象的切片的二维基础图像来获得所述所获取的三维基础图像。

11.根据前述任一项权利要求所述的方法，其中，所述执行磁共振图像的正则化重建的步骤基于冗余小波框架表示。

12.根据权利要求1至10中任一项所述的方法，其中，所述执行磁共振图像的正则化重建的步骤基于非冗余小波表示。

13.根据前述任一项权利要求所述的方法，其中，所述成本函数还包括至少一个选自以下项的空间域惩罚项：

-所述重建图像的全变差范数；和

-凸约束。

14.一种执行身体的动态的并行磁共振成像的方法，包括：

-从具有已知的或估算的灵敏度分布和噪声协方差矩阵的各个接收天线中获取所述身体的一组时间序列的基础磁共振图像，所述基础图像在k-空间中欠采样；以及

-执行所述身体的时间序列的磁共振图像的正则化重建；

其中，通过最小化成本函数进行所述执行时间序列的基础磁共振图像的正则化重建的步骤，所述成本函数包括：

-误差项，在给定相应的所获取的基础图像的情况下，所述误差项表示各重建图像的似然；和

-时间惩罚项，所述时间惩罚项表示在所述序列的连续图像之间的像素间或体素间的差异。

15.根据权利要求14所述的方法，其中，所述时间惩罚项基于边缘保持函数。

16.根据权利要求15所述的方法，其中，所述时间惩罚项基于凸边缘保持函数。

17.根据权利要求16所述的方法，其中，所述时间惩罚项基于L_p范数，其中，p≥1且优选地1≤p<1.5。

18.根据权利要求14至17中任一项所述的方法，其中，通过第一部分时间惩罚项和第二部分时间惩罚项的总和而给出所述时间惩罚项，其中：

-所述第一部分时间惩罚项表示所述序列的每个偶数图像与前一个奇数图像之间的像素间或体素间的差异；以及

-所述第二部分时间惩罚项表示所述序列的每个奇数图像与前一个偶数图像之间的像素间或体素间的差异；

通过使用针对所述第一部分时间惩罚项和第二部分时间惩罚项的邻近算子最小化所述成本函数。

19.根据权利要求14至18中任一项所述的执行身体的动态的并行磁共振成像的方法，其中，在离散的框架空间内进行所述执行磁共振图像的正则化重建的步骤，且所述成本函数还包括框架惩罚项，所述框架惩罚项表示每个重建图像的框架系数的统计分布与所述系数的先验分布之间的偏差；基于所述身体的辅助的磁共振图像估算所述重建图像的框架系数的所述先验分布。

20.根据权利要求19所述的方法，其中，所述基础图像为三维图像，且其中，所述离散的框架空间为离散的三维框架空间。

21.根据权利要求14至20中任一项所述的方法，还包括使用最大似然估算量来自动确定所述时间惩罚项的加权参数的步骤。

22.根据权利要求21所述的方法，其中，所述方法包括：估算待重建的图像的针对每个像素或体素、或者成组的邻近像素或体素的时间惩罚项的所述加权参数。

23.根据权利要求21或22所述的方法，其中，所述时间惩罚项基于L_p范数，且其中，使用所述最大似然估算量来共同地确定所述时间惩罚项的所述加权参数和参数p。

24.根据权利要求23所述的方法，其中，在约束p≥1的情况下，使用所述最大似然估算量来共同地确定所述时间惩罚项的所述加权参数和所述参数p。

25.根据权利要求14至23中任一项所述的方法，其中，所述误差项取决于几何参数，所述几何参数限定所述基础磁共振图像的每个图像相对于作为参照的基础磁共振图像的刚性变换，且其中，还通过相对于所述几何参数最小化所述函数来进行所述执行时间序列的基础磁共振图像的正则化重建的步骤。

26.根据前述任一项权利要求所述的方法，其中，通过回波平面成像获取所述基础图像。

27.根据前述任一项权利要求所述的方法，其中，所述基础图像欠采样，且缩减因数大于或等于4。

Images