---
title: "读 Lustig、Donoho 和 Pauly 2007 年的 Sparse MRI 论文"
description: "带着第五、六章的公式读压缩感知 MRI 的奠基论文：三个前提、点扩散函数与式 [1]、采样设计、重建和结果，复算干扰公式和 16 倍的结论。"
canonical_url: "https://sandboxphysics.com/zh/stories/mri-spatial-encoding/reading-lustig-2007/"
language: "zh-CN"
page_type: "story"
alternate_languages:
  en: "https://sandboxphysics.com/stories/mri-spatial-encoding/reading-lustig-2007/"
last_modified: "2026-10-29"
---

[磁共振空间编码](https://sandboxphysics.com/zh/stories/mri-spatial-encoding/)

L4 · 第 8 章

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

Sandbox Physics · 2026-10-29

第六章用到的随机跳行和稀疏重建，来自 Lustig、Donoho 和 Pauly 2007 年发表在 Magnetic Resonance in Medicine 上的论文 Sparse MRI。这一章按论文的顺序读下去，每一步都问：他们的论点靠的是哪个公式或哪组数字，哪些我们能复算，k 空间重建实验又在哪里与论文不同。

## 开始前需要什么

需要第五章的折叠和分辨率、第六章的点扩散函数、式 [1] 和稀疏重建目标函数。不需要读过论文，但手边有原文会更有收获。

- 说出论文提出的三个前提，以及它们分别对应论文的哪一部分。
- 由论文附录的能量论证复算式 [1]，并复算单层二维扫描干扰大 16 倍的结论。
- 指出论文的主要结果、作者自己承认的局限，以及本课程实验与论文的差别。

## 论文要解决什么问题

扫描时间与采集的 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$ 是缩减倍数。

$$
\sigma_{\mathrm{SPR}}=\sqrt{\frac{p-1}{D}}
$$

论文式 [1]。正文说推导在附录 II，实际在附录 B。

来源：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 小波加全变差惩罚。

$$
\min_m\ \|F_u m-y\|_2^2+\lambda\|\Psi m\|_1
$$

论文式 [A1]。

来源：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}$，与论文一致。

$$
\sigma^2=\frac{N/D-(N/D)^2}{(D-1)(N/D)^2}=\frac{D-N}{N(D-1)}\;\approx\;\frac{p-1}{D}
$$

当 $D\gg1$ 时回到式 [1]。

来源：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

## 到实验里试一试

### [在实验里重做论文图 2 的对比](https://sandboxphysics.com/zh/labs/mri-k-space-reconstruction/)

1. 选“隔行采样：头像折叠”，再切换到稀疏小波拟合。预测：折叠不会消失。读出误差约 49.8%。
2. 选“随机跳行，缺口补零”，记下误差 23.9%；切换到“随机跳行，稀疏拟合补全”，误差降到 13.3%。两次采的行数相同。
3. 在“拟合权重、中心行数与种子”里把中心始终保留的行数调到 0。预测：不保留中心行时误差变大。读数应从 13.3% 升到 18.8%，补零重建则从 23.9% 升到 30.9%。

## 练习与解答

1. 从 192 条相位编码中取 80 条，缩减倍数是多少？若把它看作一维随机采样，按式 [1] 干扰的标准差是多少？

   展开解答

   $p=192/80=2.4$，$\sigma=\sqrt{1.4/192}=0.085$。实际论文用的是变密度采样，干扰分布不均匀，这个数只是量级。

2. 论文的模拟中加速倍数增大时先丢掉的是什么？这对临床应用意味着什么？

   展开解答

   先丢掉的是低对比度特征，它们被干扰淹没后无法恢复。所以压缩感知适合血管造影这类高对比度、需要快速成像的场合；对靠细微对比度诊断的病变要谨慎。

3. 为什么作者偏好三维扫描而不是单层二维扫描？

   展开解答

   三维扫描可以在 $k_y$–$k_z$ 平面上随机欠采样，干扰分摊到二维，不相干性更高；单层二维扫描只能沿一个方向分摊，干扰大 $D^{1/4}$ 倍，$256\times256$ 的图像是 16 倍。

[上一章：扩散编码：一对梯度怎样测出水分子走了多远](https://sandboxphysics.com/zh/stories/mri-spatial-encoding/diffusion-encoding/)

## 来源与复算范围

1. [MIT 6.3000 Signal Processing, Magnetic Resonance Imaging lecture handout (6 May 2025), data and examples from E. Adalsteinsson](https://sigproc.mit.edu/_static/spring25/lectures/MRI-handout.pdf)

   通读全文。第 9 页：质子旋磁比取 42.58 MHz/T。第 11—13 页：梯度场、相位积分，以及把 k 定义为梯度对时间的积分。第 15—18 页：k 空间数据、复图像及其幅值。第 19 页：256 行、每行 2 s，约 8.5 分钟。第 20—21 页：只保留 k 空间的偶数列。第 24—36 页：两个线圈展开折叠，包括灵敏度无法区分折叠像素的那一对线圈。本课程沿用讲义中傅里叶核的符号约定。

2. [M. Lustig, D. Donoho and J. M. Pauly, Sparse MRI: the application of compressed sensing for rapid MR imaging, Magnetic Resonance in Medicine 58, 1182–1195 (2007), doi:10.1002/mrm.21391](https://doi.org/10.1002/mrm.21391)

   通读全文。图 1 与压缩感知的三个条件；图 2（256 点稀疏信号、8 倍欠采样）；式 1（干扰的标准差）及附录 B 中的推导（正文称为 Appendix II）；单层 2DFT 的干扰大 D^{1/4} 倍的说明；式 2（变换点扩散函数）；式 3—4（带约束的重建）；附录 A 式 A1—A3（拉格朗日形式与非线性共轭梯度，80 到 200 次迭代）；结果与讨论（脑部 2.4 倍、血管造影 5 倍从 22 s 缩到 4.4 s、低对比结构丢失）。在体图像只做描述，不复现。