01 压缩感知原理和建模

传统的数据采样和重构需要遵循Nyquist采样定律，即采样频率必须大于信号频率带宽的2倍，才能完整的重建信号。如果采样频率低于2倍的频率带宽，信号在频域频谱搬移后就会发生混叠，产生伪影。

压缩感知（Compressed Sensing）理论提出：如果一个信号是稀疏的，或者在其某个变换域是稀疏的，那么信号可以从远低于Nyquist采样定律的采样频率中重建出来，即稀疏或可压缩信号从其投影到某个子空间的少量值上通过非线性重建。

假设一个信号 ，存在正交变换 ，在此变换域中的稀疏表示 是稀疏的，那么可以通过一个观测矩阵 ，产生一个观测向量 ，满足，即采样值（已知）。压缩感知问题就是在已知测量值和测量矩阵的基础上，求解欠定方程组得到原信号，。

可以将该问题建模成一个优化问题，在数据保真项的约束下，通过控制的稀疏程度（实际应用中更通常使用范数），求解以下反问题，可以从少量的信号采样中，得到最终的信号重建结果：

想要通过CS理论重构出精确的原始信号，需要满足以下三个要求[2]：

①信号的稀疏性：待重建信号需要在某一个变换域（包括恒等变换）是稀疏的，即可压缩的；

②非相关性/等距约束性：非相关性要求欠采样造成的伪影和信号在稀疏表示下是不相关（noise-like）的，这样可以通过非线性算法从带有混叠伪影的信号中重构出原始信号。这一条件约束的是测量矩阵，要求观测矩阵满足有限等距性质（RIP），RIP的等价条件是观测矩阵和稀疏表示基不相关（incoherent）。独立同分布的高斯随机矩阵、随机部分傅里叶矩阵都可以作为普适的测量矩阵；

③非线性重构算法：从维测量向量重构长度为的原始信号是一个NP-hard问题，无法直接求得最优解，只能通过高度非线性的最优化算法求得其次最优解。

02 压缩感知在MRI上的应用

MR图像的变换稀疏性和MR采集的编码性质是使CS在MRI中实现的两个关键特性。

磁共振成像由于其K空间采样和磁共振图像的稀疏性，以及临床上快速成像的要求，非常适合压缩感知的应用。

大部分的磁共振图像是稀疏的（如血管造影）或者在特定的变换域如小波变换、有限差分变换下具有很好的稀疏表示。自然的满足了压缩感知的第一条要求， 稀疏性是利用压缩感知重建磁共振图像的基础 。

MRI在K空间即频率空间采集数据，再通过反傅里叶变换得到空间域坐标值，傅里叶变换矩阵作为观测矩阵可以增加非相干性，前一小节中提到随机部分傅里叶矩阵满足等距约束性，可以作为观测矩阵产生不相干伪影。然而囿于磁共振采集硬件和人体生理上的限制，采样轨迹必须是相对平滑的直线或曲线，完全的随机采样是无法实现的。如何设计合理的欠采样方案，使得重建出的图像在稀疏域上的混叠伪影尽可能少，是快速磁共振成像的一个关键问题。

除此之外，MR图像在k空间中的能量分布并不是均匀的，MR I中的大部分能量集中在靠近 k-space 的中心，并向四周衰减，因此欠采样方案的设计应该具有可变密度，在k空间中心附近采样更加密集，同时采样轨迹要尽可能的不规则（随机）以满足观测矩阵的不相干性要求，满足以上条件的欠采样越快越好。

此时压缩感知模型中的观测矩阵对应于部分傅里叶变换矩阵，其中是采样矩阵（对应于我们常说的mask），由0和1组成，是傅里叶变换矩阵，图中给出了一些欠采样方案例子。

满足稀疏性和非相干性要求的MRI重建过程可以建模成以下优化问题：

其中为感知（观测）矩阵，单通道数据情况下，多通道数据情况下，是线圈敏感度矩阵，用来合并和恢复多通道数据。

有时还会加入全变分（TV）算子作为一个惩罚项，TV被定义为图像中绝对变化的总和，目标函数变为：

其中是参数，用于权衡稀疏性约束和TV约束。

03 非线性重建算法

上述优化问题可以转化成拉格朗日形式：

或者加入TV算子：

其中和都是参数

解决该优化问题的主要难度来自于范数和变分算子等正则项的不可微性，导致无法直接通过线性最小二乘算法求解，转而寻求更多的非线性迭代算法。下面简单介绍几种常见的非线性算法：

①非线性共轭梯度法

Lustig等人[1]将小波变换作为稀疏变换应用于CS-MRI，建立如下图像重建模型，并采用非线性共轭梯度下降算法和回溯线搜索进行求解：

非线性共轭梯度法将线性共轭梯度（CG）法扩展到非线性求解问题，用CG的方式计算优化方向，然后用回溯线搜索算法计算一个可接受的步长，用梯度替代残差

对于非线性共轭梯度法的关键在于如何去计算梯度，由于范数是绝对值的和，并不是一个光滑可微的函数（这也是求解该优化问题的难点），Lustig等人用一个平滑函数$xapproxsqrt{x*x+mu } f(m)$的近似梯度：

其中是对角矩阵，对角元素

算法流程：

②迭代软阈值算法（Iterative Shrinkage Thresholding Algorithm，ISTA）

在凸优化问题中，对于可微分的目标函数，可以通过梯度下降法（gradient descent）迭代求得最优解；对于不可微分的目标函数，引入次梯度（subgradient）也可以求解，但速度比较缓慢。为此，针对一些整体不可微分，但却可以分解的目标函数，引入近端梯度法(Proximal Gradient Descent)，用来解决包含不可导项的优化问题。

对于一个目标函数可以分解成如下形式的：

其中，是凸函数且可以微分，也是凸函数但可能不可微分，仍然可以用平滑部分的二次近似来定义向最优解下降的一步。

定义近端函数为

迭代软阈值算子为近端算子的一种形式，当为范数的形式时：

为软阈值算子，

Ricardo Otazo等人[4]用L+S分解模型求解动态磁共振成像问题时用到了迭代软阈值方法，L+S模型为：

和是参数，用来权衡数据一致性与核范数和范数之间的关系。

在迭代求解L和S时用到了软阈值算子（也称为收缩算子），每一步迭代更新L和S，最后利用得到的L和S计算矩阵M

算法流程：

③增广拉格朗日和ADMM

由于范数不可微，可以采用变量分离法（Variable Splitting Method），如增广拉格朗日算法等，通过引入辅助变量，将上面带TV算子约束项的优化问题替换为一个包含辅助变量的等式约束最小化问题，Junfeng Yang等人[5]引入辅助变量和，问题建模为：

其中表示待重建图像，为到的采集信号，被定义为全变分（TV）的离散化。

增广拉格朗日方法通过增加罚项，将不光滑不可微的项转变为光滑的表示，

引入函数和分别对应项和项，

这样不可微的项就变成了一个二次凹函数，上述约束问题对应的增广拉格朗日函数为：

增广拉格朗日函数的求解可以通过下面的迭代流程：

虽然解决了不可微的问题，但迭代求解第一项仍然是很麻烦的，为了保证ALM的收敛性，在乘子迭代更新之前，需要将每个最小化子问题求解到一定的高精度，并且第一项要达到的联合最小化。此时可以利用交替乘子迭代法（ADMM），利用中变量的可分离结构，每次迭代时仅执行一次交替最小化，然后立即更新乘子，大大降低迭代成本。

算法流程：

④原始-对偶法

通过构建原始变量的对偶变量，在原始变量和对偶变量交替迭代更新[6]