Matlab单像素成像重建算法选型：五类方法对比与避坑

简介：面向单像素成像学习与算法复现的Matlab代码包，围绕同一成像框架对比差分鬼成像、梯度下降、共轭梯度下降、泊松最大似然、交替投影、稀疏表示压缩感知与全变分压缩感知七类重建方法，适合具备Matlab基础、希望理解压缩感知与迭代重建差异的学生和研究人员用于课程实验、论文复现或方案选型。zip压缩包共10个文件，以9个.m函数脚本和1个txt说明文档为主，体积约12KB，函数按算法拆分便于逐项调用和替换比较。已有2762人学习下载，说明该主题在单像素成像入门与算法对照中具有一定参考需求。读者可借助统一演示入口快速跑通各重建流程，观察不同先验与优化策略在重建效果、误差表现上的区别，并据此扩展自己的采样矩阵、目标图像或评价指标，形成可复用的实验脚本与排错思路。

1. 单像素成像算法选型：25% 采样率下差 10 dB，多半不是算法本身的锅

同一组单像素成像观测数据，我在 Matlab 里换过五套重建脚本，PSNR 从 15.2 dB 一路跳到 26.4 dB。乍看像算法碾压，可把归一化、迭代停止条件和测量矩阵的随机种子统一之后，差距缩到 4 dB 左右，其余部分全是评测口径造成的偏差。这件事之后，我每次做算法对比都先固定评测流程，再谈谁好谁坏。下面这份对比围绕 Matlab 单像素成像里最常被拿来比较的几类重建方法——线性反投影、贪婪类 OMP、压缩感知的 FISTA/TwIST、共轭梯度下降（CGLS）以及 TV 正则——给出评判维度、逐项对比和能直接照抄的脚本骨架。做重建选型的人能拿到一张决策表，正在啃压缩感知代码的人能对上参数含义，被「跑不出别人论文数字」卡住的人，多半能在避坑那一章找到自己的现场。

2. 评判维度先定死：单像素成像算法在 Matlab 里怎么比才不算白比

我见过太多“A 算法比 B 算法好”的结论，翻开源码一看，A 跑了 300 次迭代，B 只跑了 50 次；A 的测量矩阵按行做了归一化，B 没做；A 的重建结果裁到了 [0,1]，B 还带着负数。这种比法得出的排序没有参考价值。单像素成像的特殊性在于，它的“采样”是物理层面的——每个测量模式对应一次探测器读数，采多少、采什么模式，直接决定了后面所有算法的上限。所以评判维度必须围绕“给定测量预算下能做到多好”来设计，而不是围绕某个算法在某个数据集上的单点成绩。我一般从三条线切进去：质量指标怎么算、采样率与耗时怎么扫、实现复杂度能不能在 Matlab 里改得动。这三条线定完，剩下的对比才有意义。

2.1 重建质量：PSNR 要和 SSIM 配对看

PSNR 对整体亮度偏移极其敏感。单像素成像里如果重建结果做了 min-max 归一化、参考图却按原始灰度算，PSNR 能凭空掉好几 dB，而人眼看结构几乎没差别。反过来，某些平滑性很强的正则会把图像“抹”得很干净，MSE 降下去了，PSNR 上去了，但边缘全糊——这时候 PSNR 会骗你，SSIM 不会。我的做法是两个都算，并且规定 PSNR 必须在数据范围归一化到 1 的口径下计算，也就是10*log10(1/MSE)，而不是用 255 或 65535 当峰值。别小看这一个常数，它能让两份对比表差出 48 dB 的假象。评价函数的骨架我固定成下面这样，所有算法共用一份，谁也别想在自己的脚本里偷偷改口径。

% 统一评价口径：先裁范围，再算指标，所有算法共用这一个函数 function [psnr_db, ssim_val] = eval_recon(x_hat, x_ref, N) x_hat = min(max(real(x_hat(:)), 0), 1); % 卡到 [0,1]，剔除虚部与越界 x_ref = min(max(x_ref(:), 0), 1); mse = mean((x_hat - x_ref).^2); psnr_db = 10 * log10(1 / mse); % 峰值取 1，不是 255 ssim_val = ssim(reshape(x_hat, N, N), reshape(x_ref, N, N)); end

函数里real()那一步不是多余的：共轭梯度下降在条件数差、迭代次数又给得多的时候，数值误差会累积出微小虚部，直接进ssim会报类型错误。SSIM 依赖图像处理工具箱，窗口用默认的高斯窗即可，但要注意两张图尺寸和数据类型必须一致，double对uint8混着传是另一个经典翻车点。

2.2 采样率-耗时曲线比单点耗时更有价值

只看“某个算法跑一张 64×64 图要几秒”没有决策价值，因为耗时随采样率 m/n 变化很大。共轭梯度下降的每次迭代成本是两次矩阵-向量乘，与 m 成正比；FISTA 还要多一次 DCT 和一次逆变换，开销随 N 增长；TV 类因为要解内层子问题，耗时对正则权重相当敏感。我的做法是把采样率从 5% 扫到 50%，步长取 5%，每个点跑三次取中位数，画两条曲线：PSNR 随采样率的变化、耗时随采样率的变化。真正的选型依据是“要达到 25 dB，需要多少采样、多少秒”，而不是“谁在 25% 采样率下高 0.3 dB”。这条曲线还能暴露一个隐藏特性：反投影的 PSNR 曲线在采样率低于 40% 时会急剧下降，20% 附近几乎不可用；而压缩感知类在 15%～25% 区间仍然有明显的上升趋势；共轭梯度下降介于两者之间，斜率最平缓，也就意味着它对采样预算的变化最不敏感。

2.3 实现复杂度：能不能在 Matlab 里改得动

这一条最容易被忽略，却决定你三个月后还愿不愿意用这套代码。核心判断是：算法能不能写成矩阵自由（matrix-free）的形式，也就是只需要提供A*p和A'*r两个算子，而不必显式存下整个测量矩阵。算一下内存就明白了——64×64 图像、25% 采样率，测量矩阵是 1024×4096，约 33 MB，随便存；128×128 时变成 4096×16384，约 537 MB，开着 Matlab 再开个浏览器就开始卡；256×256 直接冲到 4 GB 以上，普通笔记本当场翻车。反投影和 OMP 基本要求显式矩阵，共轭梯度下降和 FISTA 可以完全用函数句柄实现，TV 类介于中间。所以在选型表里我专门留了一列“是否需要显式测量矩阵”，这一列的权重在图像尺寸超过 128×128 之后会迅速上升，甚至压过 PSNR 那 0.5 dB 的差距。

3. 五类重建算法逐项对比：反投影、OMP、FISTA、共轭梯度下降与 TV

把维度定死之后，逐项看每个算法在单像素成像这个具体场景下的脾气。需要说明的是，下面所有结论都基于 64×64 灰度图、高斯随机测量矩阵、稀疏基取 DCT 这一套设定；换成部分 Hadamard 或者小波基，排序会有变化，但每个算法的“性格”不会变。我把判断依据写在每一项里，你可以据此判断自己的场景落在哪一档。

3.1 线性反投影：快得离谱，但采样率是硬门槛

反投影（关联成像）的一行公式就是x = Phi' * y，本质是把每次测量结果按对应的模式加权回投到像素上。它的优点没有对手：一次矩阵-向量乘，64×64 图上是毫秒级，不需要任何正则参数，不需要迭代，也不需要调参。缺点同样明确——它隐含假设测量模式之间近似正交，只有在采样接近满采（m 接近 n）时才成立。在 20% 采样率下，反投影出来的图基本上是噪声上浮着一点轮廓，PSNR 掉到 15 dB 上下，人眼只能看出“这里有个亮块”。所以在我的选型表里，反投影唯一的适用场景是实时预览：先投几十个模式，用反投影出一张糊图，确认视场和目标位置对不对，再决定要不要花时间做迭代重建。把它当成最终重建结果去和别人比 PSNR，是很多对比表失真的起点。

3.2 贪婪类 OMP：稀疏度已知时才值得拿出来

OMP 的逻辑是每次挑一个与当前残差最相关的原子，解一次最小二乘，更新残差，重复 K 次。它天然适合压缩感知，因为压缩感知的理论前提就是信号在某基下 K-稀疏。问题在于 K 从哪来——自然图像在 DCT 下的系数并不是严格稀疏的，nnz算出来几千个都不止，你只能人为设一个 K 当超参数。K 设小了丢细节，设大了计算量上去还引入噪声原子。计算成本也不友好，每轮都要算一次Phi' * r，总体是 O(K·m·n)，64×64 上跑 K=200 大概几秒，128×128 就要以分钟计。我一般只在两类情况下用它：一是被测目标本身就是稀疏的（比如几个点目标、二值图案），二是要拿它当“稀疏先验到底有多大用”的对照组。日常重建选型，它很少是最终答案。

3.3 压缩感知 FISTA/TwIST：低采样率下的主力，正则参数带点玄学

FISTA 是近端梯度法，每次迭代做一次梯度下降再加一次软阈值，收敛速度比 ISTA 快一个量级。它解决的是min 0.5*||Phi*x - y||^2 + tau*||Psi*x||_1，其中Psi是稀疏基（DCT 或小波），tau控制稀疏强度。这是单像素成像在低采样率下最常被引用的方案，15%～25% 采样率仍能出结构完整的图。它的麻烦全在参数上：tau的量级和测量矩阵的谱范数绑定，换一个归一化方式，最优tau就变一个数量级；迭代次数给少了没收敛，给多了会在最优解附近来回蹭。我一般先估L = norm(Phi)^2定步长，然后按tau ∈ {1e-3, 3e-3, 1e-2, 3e-2} * L扫四个点，取 SSIM 最高的那个，这套流程看着不优雅，但比拍脑袋定参数可靠得多。至于它和 TwIST 的差别，在单像素成像这种测量数远小于像素数的场景里往往被正则参数的影响盖过去，先把这个参数调明白比换算法更划算。

3.4 共轭梯度下降 CGLS：不假设稀疏，靠早停和阻尼当正则

共轭梯度下降在单像素成像里常被误解。严格说它不是一种“重建算法”，而是一个线性方程组求解器——用来解最小二乘的正规方程(Phi'*Phi + lambda*I) x = Phi'*y。如果只做无正则的最小二乘，欠定条件下解不唯一，直接跑满迭代会过拟合噪声，这就是半收敛现象：PSNR 先升后降。所以真正好用的形态是 CGLS 加 Tikhonov 阻尼，再加一个残差范数触发的早停。它的优势在于完全不假设稀疏基，对纹理丰富、DCT 下并不稀疏的图像反而比 FISTA 更合适；每次迭代只要两次矩阵-向量乘，64×64 上跑 50 次迭代通常在 0.1 秒量级，比 FISTA 快一个数量级。还有一个容易被忽略的事实：TV 正则的内层子问题，主流实现里也是用共轭梯度下降或者它的变体去解的，所以你在很多论文的“TV 重建”代码里翻到 CG 循环，一点都不奇怪。

3.5 对比表与逐项结论

下面这张表是我在 64×64、高斯随机测量矩阵、25% 采样率下的量级参考，用来判断排序和比例关系，绝对值会随图像内容和随机种子变化，不要当成标准答案抄。耗时列是单帧重建的中位数量级，机器不同浮动很大。

从这张表里能读出四条结论。第一，反投影和其余四类根本不在一个赛道，它只解决“有没有”的问题，不解决“好不好”的问题。第二，低采样率下 TV 和 FISTA 处在同一梯队，TV 通常高 1～2 dB，代价是十几倍的耗时。第三，CGLS 是性价比的最好选择：PSNR 只差 2 dB，速度快十倍，而且它对参数的敏感程度远低于 FISTA，换个场景不用重调。第四，OMP 在这套设定下几乎不占优势，除非被测目标天生稀疏。

4. Matlab 复现：从测量矩阵到重建脚本的完整骨架

对比要能复现才叫对比。这一章给出四个可直接跑的最小骨架，覆盖观测生成、共轭梯度下降重建、压缩感知对照和批量扫采样率。所有代码都假设图像边长 N 为 2 的幂，图像先归一化到 [0,1]，测量矩阵按行归一化。

4.1 统一观测生成：固定随机种子和归一化方式

观测生成本身不复杂，坑在于“每个算法各自生成一遍”。只要有一处rng没固定、一处归一化方式不同，后面的对比全部作废。所以我把观测生成单独抽成一个函数，任何算法都只读它输出的Phi和y。

% 生成统一的单像素观测：固定种子 + 按行归一化，所有算法共用 function [Phi, y, x_true] = make_obs(N, ratio, seed) n = N * N; m = max(1, round(ratio * n)); % 采样次数 = 采样率 * 像素数 rng(seed); % 固定种子，保证可比 Phi = randn(m, n) / sqrt(m); % 高斯随机矩阵，按行能量归一化 img = imresize(imread('cameraman.tif'), [N N]); x_true = double(img(:)) / 255; % 统一到 [0,1] y = Phi * x_true; % 单像素探测器的一维读数 end

Phi除以sqrt(m)是为了让Phi'*Phi的量级接近单位阵，这样 FISTA 的步长和 TV 的正则权重在不同采样率下不用重新标定。如果你打算换成部分 Hadamard 测量矩阵，注意 Hadamard 的行是 ±1，探测器读数会出现负值，实际光学系统里要靠差分测量实现，仿真阶段可以直接算，但报告结果时要说明这是理想差分。x_true必须和重建结果走同一个裁剪函数，否则 PSNR 没有可比性。

4.2 共轭梯度下降重建：CGLS 的循环、停止条件和参数怎么改

下面这段是我用了很久的 CGLS 实现，阻尼项直接写进正规方程，早停靠残差范数的相对变化判断。

function [x, res_hist] = cgls_solve(Phi, y, max_iter, lambda) % CGLS：解 min ||Phi*x - y||^2 + lambda*||x||^2 % max_iter: 迭代上限，通常 30~80 就够 % lambda: 阻尼系数，取 1e-3~1e-1 * trace(Phi'*Phi)/n 量级 n = size(Phi, 2); x = zeros(n, 1); r = y; % 初始残差 = y（因为 x0 = 0） s = Phi' * r - lambda * x; % 正规方程的梯度 p = s; gamma = s' * s; res_hist = zeros(max_iter, 1); for it = 1:max_iter q = Phi * p; alpha = gamma / (q' * q + lambda * (p' * p)); x = x + alpha * p; r = r - alpha * q; s = Phi' * r - lambda * x; gamma_new = s' * s; res_hist(it) = norm(r); if sqrt(gamma_new) < 1e-6 * sqrt(gamma) % 相对收敛判据 res_hist = res_hist(1:it); break; end p = s + (gamma_new / gamma) * p; gamma = gamma_new; end end

三个参数的实际影响：max_iter不是越大越好，前面说的半收敛在这里就是体现在res_hist上——残差单调下降，但 PSNR 会在某个迭代点之后开始掉，我一般画一次 PSNR-迭代曲线，把峰值位置记下来当作该采样率下的默认迭代数。lambda控制阻尼强度，取 0 就是纯最小二乘，欠定时前几步就发散；取太大等于把解往零拉，图像整体变暗。经验值是从1e-3开始按 10 倍往上试。相对收敛判据里的1e-6不用调得太细，因为这个循环通常是某个外层算法（比如 TV 的 ADMM）的内层，内层解得过精反而拖慢整体。

4.3 压缩感知对照实现：FISTA 加 DCT 软阈值

对照实现要保证和 CGLS 读同一份观测、用同一个评价函数，否则前面所有的对比工作白做。

function x = fista_dct(Phi, y, N, tau, max_iter) % FISTA：解 min 0.5*||Phi*x-y||^2 + tau*||DCT(x)||_1 % tau 与测量矩阵谱范数绑定，建议按 L 的 1e-3~3e-2 倍扫 L = norm(Phi)^2; % 谱范数平方，用于定步长 x = zeros(N*N, 1); z = x; t = 1; for it = 1:max_iter g = Phi' * (Phi * z - y); % 数据保真项的梯度 v = z - g / L; c = reshape(dct2(reshape(v, N, N)), [], 1); c = sign(c) .* max(abs(c) - tau / L, 0); % 软阈值 = L1 近端算子 x_new = reshape(idct2(reshape(c, N, N)), [], 1); t_new = (1 + sqrt(1 + 4 * t^2)) / 2; z = x_new + ((t - 1) / t_new) * (x_new - x); x = x_new; t = t_new; end end

dct2是正交变换，所以逆变换直接idct2不需要额外缩放，这一点如果用dctmtx手动构造要注意别漏了转置。阈值是tau/L而不是tau，因为软阈值算子的输入是v而不是梯度，步长已经被1/L吸收进去了，漏掉这个除法会让tau的实际含义差好几个数量级，这也是很多人“照论文设 tau 结果一片黑”的真实原因。max_iter取 100～300，再多收益很小。

4.4 批量扫采样率：把对比表自动跑出来

手动一个个改采样率再抄数字，既慢又容易抄错。写个双层循环，结果直接存表。

N = 64; ratios = 0.05:0.05:0.50; seeds = [11 22 33]; rows = {}; for ri = 1:numel(ratios) for sd = seeds [Phi, y, x_true] = make_obs(N, ratios(ri), sd); x_bp = Phi' * y; % 反投影基线 x_cg = cgls_solve(Phi, y, 60, 1e-2); x_fi = fista_dct(Phi, y, N, 1e-2 * norm(Phi)^2, 200); p1 = eval_recon(x_bp, x_true, N); p2 = eval_recon(x_cg, x_true, N); p3 = eval_recon(x_fi, x_true, N); rows(end+1, :) = {ratios(ri), sd, p1, p2, p3}; %#ok<SAGROW> end end

seeds取三个以上再对结果取中位数，能显著压掉随机测量矩阵带来的抖动。我试过只用一个种子，同一算法在相邻两个采样率点上的 PSNR 能出现非单调的跳变，看着像 bug，其实是种子导致的。另外这段循环里的三个算法共用同一个Phi，这是对比有效性的前提，别图省事在每个算法内部各生成一次。

5. 避坑与排查：单像素成像重建里最容易翻车的五个现场

这一章记录的五个现场，基本覆盖了我这几年在 Matlab 里做单像素成像重建时被绊倒的全部类型。每条按现象、原因、解决三段来写，照着对号入座通常能省掉半天时间。

5.1 现象：换个算法 PSNR 差 8 dB，图像看起来却差不多

原因几乎总是归一化口径不一致。反投影出来的结果天然带着测量矩阵的尺度因子，Phi'*y的数值范围可能是 [-0.3, 1.7]，直接和 [0,1] 的参考图算 MSE，误差主要来自整体偏移而不是结构差异；而 FISTA 的结果因为软阈值的存在，天然被压到接近 [0,1]。两边不裁剪就比，等于拿两把不同的尺子量同一块布。解决办法是把裁剪写进评价函数本身，而不是指望每个重建函数自己收尾——eval_recon里min(max(...))那两行就是干这个的。我后来还额外加了一条：先对重建结果做一次线性回归把增益和偏置对齐，再算 PSNR，这样连尺度因子都不影响结论。

5.2 现象：CGLS 迭代到 200 次还在震荡，PSNR 反而比 50 次时低

这是半收敛，不是代码写错了。欠定最小二乘的迭代解会先快速拟合信号的主要成分，然后开始拟合观测噪声，残差范数一路下降但真实误差开始上升。判断方法很简单：把每次迭代的 PSNR 存下来画一条曲线，看到明显的单峰就是它。解决办法有三条，可以叠加用：一是早停，把峰值处的迭代数当默认值；二是加 Tikhonov 阻尼，lambda取1e-2量级通常能把曲线压平；三是换用带稀疏先验的方法，从根上避免这个问题。我一般先加阻尼，因为改一个参数就能解决，比重新选算法便宜。

5.3 现象：10% 采样率下 TV 或 FISTA 出现明显块状伪影

块状伪影有两个来源，得分开排查。如果伪影是 8×8 的方格，那是 DCT 分块导致的，说明你用的稀疏基和图像内容不匹配，换成小波基或者整图 DCT 会好转。如果伪影是大尺度的平滑块、边缘被抹成台阶，那是正则过强，tau或 TV 权重给大了，把细节当成噪声压掉了。排查顺序建议先扫正则参数：固定其他条件，把tau按 10 倍量级从1e-3*L扫到1e-1*L，看 SSIM 的峰值在哪。如果扫完所有tau都救不回来，那才是稀疏基的问题。另外 10% 采样率本身已经接近这类方法的极限，25 dB 以下的 PSNR 别指望能出干净的边缘。

5.4 现象：耗时统计完全不可比，同一段代码两次跑差三倍

Matlab 的首次函数调用要加载、要 JIT，计时器精度也有限。第一次跑cgls_solve可能 0.8 秒，第二次同样的输入 0.2 秒，这不是算法变快了。解决办法是预热：正式计时前先空跑一次，把函数和内存池都热起来；计时用timeit而不是tic/toc，它会自动多次采样；如果装了并行工具箱，注意 BLAS 的多线程会让小矩阵运算的耗时抖动变大，做对比时把线程数固定住。还有一点，别在计时循环里做imshow或者save，I/O 的开销会把你测的东西全盖掉。

5.5 现象：复现别人的曲线，整体低 3～5 dB，形状却对得上

形状对得上说明算法逻辑没问题，整体下移说明输入的统计特性不一样。三个最常见的差异点：测量矩阵类型（高斯、伯努利、部分 Hadamard 三者的重建难度差异在低采样率下可以到 2 dB）、随机种子（单个种子的偏差能到 1 dB 以上）、测试图像是否做了对比度拉伸。我的做法是把这三件事写成一个配置结构体，跑任何对比实验都从同一个结构体读参数，换图像、换种子、换矩阵类型只改这一处。这样即使复现不出别人的绝对值，至少自己的曲线之间是可比的，而选型决策需要的恰恰是相对排序。

6. 结论与上手路径：按采样率分三档选，再加一个 30 秒自检

把前面所有内容收敛成一张决策依据，按你手上的采样预算分三档。

选完之后还有一个 30 秒自检，我每次交付重建代码前都会走一遍：拿同一份观测，把重建结果和参考图并排imshow一次，看整体亮度有没有系统性偏移；再把残差y - Phi*x_hat的范数打出来，和norm(y)比一下，如果残差范数接近观测范数，说明这个“重建结果”基本没拟合上数据，PSNR 再好看也是归一化凑出来的；最后把max_iter翻倍再跑一次，看 PSNR 变化是否小于 0.2 dB，如果还在明显上升，说明之前根本没收敛，那张对比表得重做。

最后这个自检习惯是被坑出来的。早些年我做一份对比表，反投影那一栏的 PSNR 高得反常，查了两天才发现是重建结果忘了裁剪，负值在 MSE 里被平方放大后反而“凑”出了一个好看的数。从那以后我每次跑完重建，都强制把残差范数和裁剪后的图像范围打印出来看一眼，再谈指标。血泪经验谈不上，但这一步确实帮我省掉了好几次返工。希望帮到你。

本文还有配套的精品资源，点击获取