读 Lustig、Donoho 和 Pauly 2007 年的 Sparse MRI 论文

论文要解决什么问题

扫描时间与采集的 k 空间样本数成正比，而 Nyquist 采样要求的样本数决定了常规扫描有多慢。论文的出发点是：MR 图像在合适的变换下往往是稀疏或可压缩的，大部分系数接近零。只要满足三个前提，就可以少采样本，再把图像算回来。

三个前提是：（a）图像在已知变换域中有稀疏表示；（b）k 空间欠采样造成的混叠在该变换域中不相干，像噪声一样；（c）用非线性重建同时保证表示的稀疏性和与采集数据的一致性。论文的图 1 画出了图像、k 空间和变换域三者之间的关系。

来源：M. Lustig · MIT 6.3000 Signal Processing

不相干：点扩散函数和式 [1]

图 2 用一个简单的例子说明思路：一个稀疏信号被 8 倍欠采样。等间隔欠采样并补零得到一组平移副本，原信号与副本同样可能，无法区分；随机欠采样得到的伪影像随机噪声，但它们其实是每个非零系数的能量漏到别处。强的分量高出干扰，可以用阈值找出来，算出它们造成的干扰并减掉，再找较弱的分量。

为了量化，论文定义了点扩散函数 PSF 和变换点扩散函数 TPSF（式 [2]）：一个变换系数经过欠采样以后，有多少漏到第 $j$ 个系数上。从 $D$ 个网格点里随机取 $N$ 个时，点外的值的实部和虚部各自很像零均值的白高斯噪声，标准差由式 [1] 给出，$p=D/N$ 是缩减倍数。

来源：M. Lustig · MIT 6.3000 Signal Processing

笛卡尔采样的三种情形和采样设计

论文比较了三种笛卡尔采样。单层二维扫描只能在相位编码方向欠采样，干扰只沿一个方向扩散，按式 [1] 比纯二维随机采样大 $D^{1/4}$ 倍，所以只能期望不大的加速。多层二维扫描可以让每层取不同的相位编码行。作者偏好三维扫描：在 $k_y$–$k_z$ 平面上随机欠采样，能达到二维的高不相干性。

采样设计用 Monte Carlo 方法：按期望的分辨率和视野定网格，构造一个随离原点距离按幂次下降的概率密度，从中随机抽取样本；因为抽样是随机的，可能碰上不好的样本，所以重复多次，选 TPSF 峰值副瓣最小的那一组。变密度采样的理由是图像能量集中在 k 空间原点附近。

来源：M. Lustig · MIT 6.3000 Signal Processing

重建：约束优化和共轭梯度

重建是式 [3]：在与测量数据相差小于 $\varepsilon$ 的图像中，找 $\|\Psi m\|_1$ 最小的那个，$\varepsilon$ 通常设得低于预期噪声。附录 A 把它写成拉格朗日形式 [A1]，选 $\lambda$ 让残差约等于 $\varepsilon$，用带回溯线搜索的非线性共轭梯度求解；文中例子需要 80 到 200 次迭代。他们的 Matlab 实现重建一幅 $480\times92$ 的血管造影（5 倍加速）大约迭代 150 次，用时 30 秒。

稀疏变换因应用而异：小波、有限差分（即全变差），或两者结合。脑部例子用 Daubechies 4 小波加全变差惩罚。

来源：M. Lustig · MIT 6.3000 Signal Processing

结果和作者承认的局限

模拟：$100\times100$ 的体模上随机放 18 个特征，6 种大小、3 种强度，575 个非零像素。8 倍和 12 倍加速（12 倍只限变密度采样）时，压缩感知重建完全恢复；加速继续增大时，丢掉的是低对比度特征，而不是通常的分辨率。

实测：多层二维快速自旋回波的脑部扫描，每层从 192 条相位编码中按二次变密度随机取 80 条，加速 2.4 倍；三维腹部血管造影事后欠采样 5 倍，相当于把扫描时间从 22 秒缩到 4.4 秒。

讨论部分写得很直接：随着加速增大，压缩感知最明显的伪影不是通常的分辨率下降或混叠干扰增加，而是低对比度特征消失。所以它特别适合有高分辨率、高对比度特征又需要快速成像的场合。

来源：M. Lustig · MIT 6.3000 Signal Processing

完整例题：复算式 [1] 和 16 倍

附录 B 的论证只有几行。把点扩散函数的总能量记作 $N/D$，主瓣能量是 $(N/D)^2$，点外的能量就是两者之差。除以点外像素数和主瓣能量，再令 $p=D/N$，就得到式 [1]。若点外像素数按精确的 $D-1$ 计算，得到下面右边的式子。

取 $D=256$、$N=32$：式 [1] 给出 0.165359，精确式给出 0.165683。本课程的基准程序抽了 400 组随机样本，测得 0.165683，与精确式一致；实部的超额峰度是 −0.057，接近高斯分布的 0，与论文“很像白高斯噪声”的说法相符。

再复算 16 倍。$256\times256$ 的图像有 $D=65\,536$ 个像素。纯二维随机采样时干扰分摊到全部像素上；单层二维扫描只能沿一个方向的 256 个像素分摊。两者之比 $\sqrt{65\,536/256}=16=D^{1/4}$，与论文一致。

来源：M. Lustig · MIT 6.3000 Signal Processing

本课程的实验与论文有什么不同

k 空间重建实验沿用了论文的思路：随机选相位编码行、中心 8 行全采、其余按变密度抽取，再解拉格朗日形式的目标函数。但它在几处做了简化，结果不能直接与论文比较：矩阵只有 $64\times64$，单层、单线圈，加速 2 倍；稀疏变换是三层 Haar 小波，没有全变差；求解用加速软阈值迭代 50 步，不用共轭梯度；权重由留出行选择，不按噪声水平选。

在这些简化下，补零重建的误差 23.9%、稀疏重建 13.3%，小圆点对比度从 71% 回升到 87.3%，方向与论文一致：稀疏重建去掉了大部分不相干伪影，但对比度没有完全恢复。规则跳行时稀疏重建同样无效，误差停在 49.8%，这对应论文图 2 中等间隔欠采样无法恢复的情形。

来源：M. Lustig · MIT 6.3000 Signal Processing