【论文阅读】Calibrationless Parallel Imaging Reconstruction Based on Structured Low-Rank Matrix Completion

【论文阅读】Calibrationless Parallel Imaging Reconstruction Based on Structured Low-Rank Matrix Completion

提出了一种新的算法 ：SAKE ，同步校正和k空间估计

本篇博客只重点介绍关于论文中所使用的算法，具体实验和附录公式的推导请查询原文。
采样不足的多通道数据及被重建成一种单一的数据矩阵，然后该重建问题被表述为一个低秩矩阵的补全问题，重建结果也证明了回顾性和前瞻性采样不足，多通道的笛卡尔数据无校准信号，此外，本篇文章也提出了一种非笛卡尔数据的重构方法，最后将SAKE算法与基于小波的压缩感知所结合，证明了改进后图像的质量。

理论部分

文章中首先定义了一种被称为“数据矩阵”的增强数据结构，并描述了它的结构化、低秩性。

然后，从数据矩阵的角度讨论了类似于grappa的自动标定方法。最后对所提出的SAKE方法进行了说明。

结构化数据矩阵

将数据矩阵结构化的算法步骤如下图所示：

设原始数据的矩阵维度为 $N_x\times N_y$ ，线圈的数量为$N_c$，那么我们可以通过滑动一个大小为$w\times w\times N_c$的窗口来构造出一个$w^2{N}_c\times (N_x-w+1)(N_y-w+1)$的二维矩阵。

文中的附录中，假设线圈灵敏度具有紧凑的k空间支持，作者证明了块Hankel矩阵形式的数据矩阵通过选择合适的$w$的值被结构化成为一个秩亏矩阵，并且能够进一步证明这个矩阵的秩满足$rank(A)<(w+s-1)$，其中$A$表示的是数据矩阵，$s$表示的是k空间像素测量的线圈带宽。

一旦有了一个低秩的矩阵，我们就可以使用奇异值分解将MR数据集中的信息进行分解，将其分解成一个信号空间和噪声子空间。

基于k空间的并行成像重建

基于k空间的自动校准：GRAPPA或者是SPIRiT，这两种算法通过拟合GRAPPA kernel weight 来估计k空间样本中的线性依赖关系。在本节的重建过程中，假设k空间中所有位置的依赖关系都是相同的，通过对所有线圈上附近的k空间的点施加线性权重来合成未采样的数据点。

将多通道的ACS数据组织成数据矩阵(有时被称为校准矩阵)，这样就很容易估计出GRAPPA/SPIRiT核的线性权值。这样我们假设获得了足够多的ACS（自动校准信号），设$A^{ACS}$是标定矩阵，然后这个估计线性权值GRAPPA标定的过程可以表述为：
$$g_{ir}^H{A}^{ACS}=e_i^H{A}^{ACS}$$

在公式(1)中，$g_{ir}^H$是第i个通道的一个GRAPPA核，它在适当的位置包含线性权值和0值。$r$值表示在特定的采样模式下的索引。向量$e_i$是一个来自规范基的向量，它是在矩阵$A^{ACS}$中选择的一行数据。我们使用$g_{ir}^H$表示$g_{ir}$的复转置矩阵，在SPIRiT成像中，无论采样模式如何，都能找到周围样本的线性系数，所以公式(1)的$r$被忽略，公式(1)被重新写成：
$$(g_{ir}-e_i{)}^H{A}^{ACS}=0$$
从公式(2)中可以看出GRAPPA/SPIRiT核减去$e_i$后就是校准矩阵$A^{ACS}$的左空向量，所以我们可以将校准步骤看成在$A^{ACS}$的噪声子空间中寻找一组具有代表性的零向量的过程 ，噪声子空间就是矩阵$A^{ACS}$的左零空间。如果选择合适的核窗口大小的话，就能够使得$A^{ACS}$低秩，并且总是具有一个非平凡的左零空间。

GRAPPA算法中的假设是从ACS中估计出的线性依赖关系应该在整个k空间中保持不变，所以公式(2)能够进一步被写成：
$$(g_{ir}-e_i{)}^{H}A=0$$
这里的$A$矩阵是k空间组成的数据矩阵。公式(3)构成了GRAPPA/SPIRiT中最基本的机制，为获取数据的重构提供了基础。公式(3)说明了在k空间中的任意被向量化的数据块都能被向量$(g_{ir}-e_i)$通过内积变为空，这进一步说明了在k空间中的任意缺失的数据都能够按照满足要求的方式被合成出。

基于结构低秩矩阵补全的并行成像重建

前面讨论的方法都假设我们有ACS，并且能够从中获取额外的子空间信息。然而当采样不足没有ACS时，我们无法估计校正矩阵的子空间。我们现在所知道的先验信息是数据矩阵是一个低秩矩阵，所以SAKE算法的表述可以是：当k空间的欠采样只给出矩阵$A$的一个子集时，我们来恢复结构化的低秩矩阵$A$。

我们先来定义线性算子：

• 从多通道的k空间数据集到数据矩阵：
  $$H_W:{\mathbb{C}}^{Nx\cdot N_y\cdot Nc\times 1}\to {\mathbb{C}}^{w^2{N}_c\times \left(N_x-w+1\right)\left(N_y-w+1\right)}$$

• 从数据矩阵到响应的k空间数据集：
  $$H_w^{†}:{\mathbb{C}}^{w^2{N}_c\times \left(N_x-w+1\right)\left(N_y-w+1\right)}\to {\mathbb{C}}^{Nx\cdot N_y\cdot Nc\times 1}$$

$†$ 表示的是伪逆算子，而$H^{†}$的作用首先是将数据矩阵中来自相同k空间位置的多个反对角线上的元素求平均值来强制Hankel矩阵结构，一旦数据矩阵有着这种结构，$H^{†}$就将计算出的平均值储存在k空间中。

有文献详细讨论了如何在矩阵中实现线性算子$H$和$H^{†}$：

在实践中，由于数据量大，计算伪逆算子复杂，建议避免显式实现矩阵。相反，我们将伪逆实现为一个运算符，以更有效地完成相同的计算。作者在附带的源代码中提供了一个示例(http://www.eecs.berkeley.edu/?mlustig/)。

通过公式(5)，并行成像重建能够表述成一个结构化低秩矩阵补全问题：
$$\begin{gathered} minimizerank(A) \\ subjecttox=H^{†},{\Vert Dx-y\Vert }^2<\epsilon \end{gathered}$$
公式(6)中的$D$是一个线性算子，它能够将重建的k空间数据$x$和获取的数据$y$联系起来，$\epsilon$是噪声的界限。总结下来就是，我们找一个低秩的矩阵$A$，将其转化为k空间的数据$x$，与根据$D$定义的采样机制获取的$y$数据一致，其中$A$转化为$x$的过程依据的是结构一致性，而根据$D$得到$x$的过程依据的是数据一致性。公式(6)中的数据一致性在笛卡尔采样和非笛卡尔采样中是一致的，在本文的附录中有证明过程。

当公式(6)中$A$的秩已知时，公式(6)就能够被重新写为：
$$\begin{gathered} minimize{\Vert Dx-y\Vert }^2 \\ subjecttorank(A)=k,x=H^{†} \end{gathered}$$
其中$A$矩阵的$k$值需要使用Cadzow算法(参考文献17，26)进行估计，在实现整个算法时，文章定义了一下的投影操作符，对应于公式(7)中的每一个约束：

• Low-rankness projection: 由k空间数据的当前估计构造的数据矩阵的硬阈值奇异值。
• 结构一致性投影：将数据矩阵投射到分块Hankel矩阵的空间上。
• 数据一致性投影：将k空间数据的当前估计投影到问题$Dx=y$的最小二乘解集上。附录中详细介绍了笛卡尔采样和非笛卡尔采样情况下实现投影的详细信息

在每个循环迭代的过程中，这些投影被依次利用，直到算法收敛。算法迭代的过程如下图所示：

伪代码如下所示：

参考文献

Shin P J, Larson P E Z, Ohliger M A, et al. Calibrationless parallel imaging reconstruction based on structured low‐rank matrix completion[J]. Magnetic resonance in medicine, 2014, 72(4): 959-970.