简介面向单像素成像与压缩感知重建学习者的Matlab算法对照包聚焦在同一实验框架下比较不同SPI重建思路的表现适合需要理解算法差异、复现实验或做课程与科研验证的本科生、研究生及工程人员。包内共10个文件以9个.m函数脚本和1个.txt说明文件为主压缩包约12KB代码按重建方法拆分另含统一演示入口与误差评估函数便于逐项调用和横向比较。内容覆盖差分鬼成像、梯度下降、共轭梯度下降、泊松最大似然、交替投影、稀疏表示压缩感知和全变分压缩感知等典型方法。已有2762人学习下载读者可据此快速搭建比较实验观察各算法在重建质量与收敛特性上的差异并在此基础上修改采样矩阵、参数或评价指标作为单像素成像入门复现与算法选型的参考。1. 单像素成像选算法这件事先别急着抄代码用 Matlab 做单像素成像算法比较最容易翻车的地方往往不是算法选错而是比较时变量没控住。我第一次把共轭梯度下降和一组全变分重建并排跑结论是「CG 不行」半年后回去查发现是测量矩阵行序和投影图案时序错开一位谁上都得糊。从那以后我给自己立了条规矩任何一次压缩感知重建的对比先把采样率、测量矩阵类型、噪声模型写死在脚本头部再谈谁好谁坏。这篇讲的是在单像素成像这条链路上共轭梯度下降、正交匹配追踪、全变分和 FISTA 各自适合什么采样率与图案Matlab 实现里哪些细节会让排名反转以及怎么用一周时间跑出属于你自己硬件条件的结论。2. 从桶探测到图像单像素成像的重建链路拆给你看单像素成像的物理过程很朴素用空间光调制器把一组预先设计好的图案打到目标上每个图案对应一次桶探测器的总光强读数最后拿「图案矩阵 读数向量」反解出图像。真正决定成败的环节不在探测器而在反解这一步用什么算法、带什么先验。很多人在这一步直接抄一段压缩感知的重建代码跑通了就算完事结果换一组图案、换一个采样率指标就掉得没法看——不是代码写错了是链路里每个环节的耦合关系没搞清楚。2.1 一条完整的重建链路里算法到底插在哪把链路拆开看它其实只有五段。第一段是图案生成用随机二值、高斯随机、Hadamard 排序或者傅里叶基底生成 m 个 n×n 的图案第二段是调制与投影图案被 DMD 或 LCD 加载目标经过调制后总光强落在单点探测器上第三段是采集探测器输出模拟电压ADC 量化成一个标量读数第四段是把 m 个读数拼成向量 y把 m 个图案拉直拼成测量矩阵 A得到方程 y A x n第五段才是重建算法登场从欠定的 y 和 A 里恢复 x。算法只负责第五段但它能做的事情被前四段死死限制住。A 的质量决定了解的上限图案如果不是理想的 ±1 正交基底或者 DMD 的二值化把灰度图案截断失真矩阵里就混进了系统性误差再好的算法也只是在拟合误差。探测器噪声模型同样如此光子计数低的时候是泊松噪声读出电路主导的时候近似高斯噪声两种噪声下同一套正则参数的取值能差一个数量级。所以「算法比较」这个说法本身要拆细。我一般会把它拆成三件套来比先验是什么无先验、稀疏先验、全变分先验、算子是什么显式矩阵、函数句柄、快速变换、求解器是什么共轭梯度、贪婪追踪、近端梯度、增广拉格朗日。同一个 TV 先验用 ADMM 解和用原始-对偶解收敛速度可以差两三倍同一个共轭梯度作用在显式矩阵上和作用在快速变换算子上内存占用差一个量级。只看算法名字做对比等于只看发动机型号不看变速箱和车重。2.2 比算法之前先把三类前置变量钉死我见过太多对比实验两组结果的差异其实来自测量矩阵不同却被归因到算法头上。要让比较有意义下面三类变量必须先固定且固定在脚本最顶部用一个结构体统一管理换值只改一处。前置变量常见取值对算法排名的影响我的固定建议测量矩阵高斯随机、伯努利、Hadamard 排序、小波、傅里叶基底相干性低的矩阵让 OMP 类贪婪算法占便宜结构化矩阵让快速算子占便宜基线统一用归一化高斯随机行向量除以 sqrt(m)采样率5% 到 50%按 N 的百分比取整低于 15% 时无先验的共轭梯度基本失效TV 和 L1 才有优势基线取 25%再扫 5/10/15/30/50 五个点噪声模型无噪仿真、加性高斯、泊松光子噪声泊松噪声下 L1 类算法对低光子区域的重建偏差明显更大先跑无噪验证链路再加 30dB 高斯最后上泊松这三类固定之后才轮到图像内容。我习惯准备四张测试图一张分段常数的块状图、一张自然纹理图、一张二值稀疏点阵、一张带渐变背景的目标。分段常数图用来暴露 TV 的优势和阶梯效应自然纹理图用来暴露 TV 的细节丢失二值点阵用来验证 OMP 在强稀疏下的表现渐变图用来检验各算法在低频分量上的还原能力。四张图跑完算法的大致画像就出来了比只跑一张 Lena 可靠得多。还有一件事容易漏图案的排序。Hadamard 基底如果按自然序取前 m 行低频分量集中在前面的行里重建结果会严重偏向低频看起来「平滑但糊」按随机序或者按 sequency 排序取频率分布均匀得多。这个差别会让同一个算法在同样的采样率下 PSNR 差 3 到 5dB足以让排名颠倒。2.3 判断一个算法好不好我只看四个指标指标不能只看重建质量单像素成像的落地场景里速度和可调性往往比多 1dB 更重要。我固定看四项。第一项是重建质量但必须看成曲线而不是单点。固定其他变量扫采样率画 PSNR 和 SSIM 随采样率变化的曲线。两条曲线的交叉点比单点数值有信息量得多很多算法在 10% 采样率下差距巨大到 40% 以上就收敛到差不多的水平。第二项是收敛成本不能只报总耗时要拆成迭代次数乘以单次迭代的矩阵向量乘次数再乘以单次乘法的耗时。共轭梯度每次迭代两次算子作用FISTA 两次TV 的 ADMM 每次迭代四到六次这个系数在比较时经常被忽略。第三项是参数敏感度也就是正则参数在多大范围内波动时PSNR 下降不超过 1dB。敏感度高的算法在真实系统里没法用因为真实系统的噪声水平每次都在变你不可能每次都手工调参。第四项是内存与可扩展性能不能从 64×64 撑到 128×128 甚至 256×256这一项直接决定算法能不能上真实样机。四项指标的权重取决于你的场景。做离线验证的质量权重最高做实时成像的收敛成本优先做产品化封装的参数敏感度是生死线。3. 四类重建算法逐项对比CG、OMP、TV、FISTA 谁更适合你3.1 共轭梯度下降快、稳但别指望它自带稀疏先验共轭梯度下降在单像素成像里最常见的用法是解最小二乘问题 min ||A x - y||²。它不直接对正规方程做显式求逆而是在 Krylov 子空间里构造一组共轭方向理论上 n 步内收敛实际因为舍入误差和条件数限制通常几十步就够。Matlab 里对应的入口是 pcg 和 lsqr前者解对称正定的正规方程后者直接处理长方矩阵数值稳定性更好。它的最大优点是干净没有正则参数只有最大迭代次数和容差两个旋钮调起来没什么玄学。每次迭代的算子作用次数是两次一次正向 A p一次反向 A r内存占用是 O(N)不需要存中间的大矩阵。在采样率比较高、测量矩阵条件数可控的场景里它是极好的基线和快速预览工具。问题也很明确最小二乘本身不含任何稀疏先验。当 m 远小于 n 时方程欠定解不唯一共轭梯度会收敛到最小范数解附近噪声被完整地搬进重建图像。表现就是图整体偏亮偏糊高频全是噪点采样率越低越明显。我做过一组对照25% 采样率下纯共轭梯度的 PSNR 比带小波软阈值的版本低 6dB 左右到 10% 采样率直接不可用。真要用共轭梯度做压缩感知重建得把它嵌进迭代软阈值的框架里每一轮先做一次共轭梯度下降逼近数据保真项再在小波域做一次软阈值收缩交替进行。这样既保留了共轭梯度的收敛速度又把稀疏先验加了进去代价是多了阈值参数和内外两层循环调参复杂度上来了。我的经验是如果只是想在写正式算法之前先验证测量链路通不通纯共轭梯度足够了如果要拿它当最终方案必须加先验。3.2 正交匹配追踪稀疏度高时漂亮字典一相干就崩正交匹配追踪属于贪婪类算法逻辑直白每一轮从字典里挑出与当前残差最相关的那一列把它的系数用最小二乘更新然后从残差里减掉这一列的贡献重复 K 轮K 就是预设的稀疏度。它的计算成本主要在相关计算和最小二乘更新上K 通常远小于 n所以总耗时在四类算法里最低。在单像素成像里OMP 有一个非常合适的场景目标是强稀疏的比如二值点阵、稀疏荧光标记、远距离小目标。这时候如果测量矩阵本身就是稀疏基或者目标在某个变换域里只有几十个非零系数OMP 在 5% 到 10% 的采样率下就能给出干净的结果速度快到可以边采集边重建。但自然图像是它的死穴。自然图像在小波域里是近似稀疏而不是严格稀疏系数呈幂律衰减尾部一大堆小系数。更麻烦的是单像素成像的等效字典是测量矩阵和小波基的乘积随机矩阵和小波基之间的相干性虽然不算最差但也足够让 OMP 在迭代早期选错原子。一旦选错后面每一轮都在错误的支撑集上做最小二乘误差越滚越大最后的结果是一堆分散的亮点加一块糊掉的背景。稀疏度 K 的设定也是个大坑。设小了重建欠拟合细节全丢设大了等于把噪声也当成信号选进来。我一般先用小波系数的能量累积曲线估一个 K 的量级再在附近扫三五个值取 PSNR 最高的那个。这个做法在仿真里可行在真实系统里不可行因为你事先不知道目标的稀疏度。这也是 OMP 在工程落地里用得比论文里少的原因。3.3 全变分与 FISTA两类正则化求解器的分工全变分正则化解决的是另一类图像分段常数、块状目标、二值图案、边缘锐利的工业检测对象。它的目标函数是数据保真项加上全变分项全变分项惩罚相邻像素差的 L1 范数等价于鼓励图像分块平滑。求解上常用增广拉格朗日或者 ADMM把问题拆成几个可以交替求解的子问题TVAL3 是这类思路里比较常见的实现叫法。TV 在块状图上表现非常好采样率 10% 就能把方块的边缘还原得干干净净PSNR 能比共轭梯度高出 10dB 以上。代价是阶梯效应渐变区域会被还原成一层层的台阶纹理细节被抹平人脸皮肤这类缓变区域会出现明显的块状伪影。这是 L1 类先验的通病不是实现问题调参也救不回来。正则参数 λ 的敏感度也比较高λ 大了过度平滑λ 小了噪点全回来可用的窗口通常在半个数量级以内。FISTA 走的是另一条路数据保真项加小波域 L1 范数用近端梯度下降求解软阈值算子做收缩再用 Nesterov 动量做加速收敛速率从 O(1/k) 提升到 O(1/k²)。它更适合自然图像这类在小波域近似稀疏的目标纹理保留比 TV 好但边缘锐利度不如 TV。步长取 1/LL 是 A 的 Lipschitz 常数也就是 AA 的最大特征值实际用幂迭代估十几轮就够准。步长取大了会发散取小了收敛慢这个参数不能拍脑袋。这两类算法在实现上有一个共同点容易被忽略都不需要显式存 A只要能用函数句柄算出 A x 和 A y 就行。这意味着它们可以直接对接快速变换把 O(mN) 的矩阵乘法降到 O(N log N)在 256×256 以上的分辨率上是唯一可行的路子。3.4 逐项对比表与选型结论把四类算法放在同一张表里差异就清楚了。下面的表是在 25% 采样率、归一化高斯测量矩阵、无噪条件下测的换条件排名会变但相对关系大体稳定。算法先验关键参数每次迭代算子调用适用采样率优势翻车场景Matlab 入口共轭梯度无可外挂软阈值maxit、tol240% 以上无参数、速度快、内存小低采样率下噪声放大pcg、lsqrOMP稀疏度 KK1 次相关 1 次最小二乘5% 到 15%强稀疏目标极快、结果可解释自然图像、字典相干性高自写主循环TVADMM全变分λ、μ、maxit4 到 610% 到 30%块状图、二值目标表现最好纹理丢失、阶梯效应自写 ADMMFISTA小波 L1λ、步长、maxit210% 到 30%自然图像均衡、收敛有加速λ 敏感、步长估错会发散自写主循环结论按人群分。如果你在搭第一套验证链路先用共轭梯度跑通全流程确认图案、时序、读数都对得上再换算法这一天的投入能省掉后面一周的排查。如果你的目标明确是稀疏点阵或者二值图案采样率还想压到 10% 以下OMP 和 TV 都值得试前者快后者稳块状目标优先 TV。如果你做的是自然场景成像采样率在 20% 左右FISTA 是综合最平衡的选择纹理和边缘都能兼顾。如果你的系统要求实时别在算法上死磕先把测量矩阵换成结构化快速变换把单次算子作用从 O(mN) 降到 O(N log N)这比换算法带来的提速大得多。4. Matlab 实现层的性能拆解同一算法为什么跑出两种结果4.1 测量矩阵存 sparse 还是 full先把内存算清楚单像素成像的测量矩阵尺寸是 m×Nm 是测量次数N 是像素总数。这个矩阵是稠密的稀疏存储没有意义所以问题就变成 full 存储能不能扛住。按 double 类型每个元素 8 字节算不同分辨率的账很好算。图像分辨率像素数 N25% 采样 mfull double 内存full single 内存建议64×644096102432 MB16 MBfull 存储随便跑128×128163844096512 MB256 MB用 single或改算子256×25665536163848 GB4 GB必须用函数句柄512×51226214465536128 GB64 GB结构化快速变换 分块这张表是我踩过坑之后贴在显示器边上的。64×64 的时候一切都很美好代码随手写成 A randn(m, N)跑得飞快然后直接把分辨率改成 128×128Matlab 卡死三分钟之后报内存不足。更隐蔽的情况是内存够但速度崩了单个矩阵乘法涉及 4096×16384 次浮点运算每轮迭代做两次几百轮迭代下来时间全花在等内存带宽上。正确的做法是从一开始就把 A 写成算子而不是矩阵。这不只是为了省内存更是为了让同一份重建代码能在不同分辨率下无缝切换。转 single 也要谨慎单精度在矩阵乘法里会累积误差共轭梯度这类依赖内积精度的算法迭代几十轮之后残差就可能不再下降表现为「跑到一半卡住」。4.2 算子写成函数句柄别写成显式矩阵Matlab 的 lsqr、pcg、bicgstab 都支持传函数句柄签名是 afun(x, trans_flag)其中 trans_flag 取 notransp 或 transp。把正向和反向作用包在一个子函数里外部用匿名函数捕获图案参数就能做到图案常驻内存、矩阵不落地。% 用一个子函数承载正反向算子作用外部只传句柄给 lsqr/pcg function y opA(P, x, transFlag) if strcmp(transFlag, notransp) y P * x; % 正向图像向量 - 测量向量 else y P * x; % 反向测量向量 - 图像向量 end end % 主脚本片段 n 64; N n^2; rate 0.25; m round(rate*N); rng(2024); % 固定种子保证每次比较用同一组图案 P randn(m, N, single) / sqrt(m); % 归一化避免内积量级失控 y P * double(phantom(n)(:)); % 无噪仿真测量 Afun (x, t) opA(P, x, t); % 句柄捕获 P不复制矩阵 [x, flag, relres, iter] lsqr(Afun, y, 1e-6, 200, [], [], zeros(N,1)); img reshape(x, n, n);这段代码里几个参数值得说明。rng 的种子必须固定否则每次跑出来的图案都不一样两个算法的比较就没有意义。图案矩阵除以 sqrt(m) 是归一化让 AA 的对角元素接近 1这样共轭梯度的收敛行为和正则参数的取值在不同采样率之间才可比。lsqr 的第三个参数是容差第四个是最大迭代次数第五到第七个是阻尼和初始值这里阻尼取空表示不做 Tikhonov 正则初始值全零。返回的 flag 为 1 表示达到容差为 2 表示达到最大迭代次数为 3 表示残差不再下降——看到 flag 为 3 就要回头查矩阵是不是秩亏或者量级失衡。如果图案是 Hadamard 基底根本不用生成 m×N 的矩阵直接调用快速沃尔什-哈达玛变换把向量变换后取前 m 个系数即可复杂度从 O(mN) 降到 O(N log N)。这个改动在 256×256 以上是数量级的差别值得单独花半天实现。4.3 停止准则与正则参数最容易被忽略的调参旋钮共轭梯度的停止准则用相对残差别用绝对残差。abs(norm(A*r)) tol 这种写法在图像动态范围变化时会完全失效换成 norm(r)/norm(y) tol 才稳定。tol 取 1e-6 是常见起点但压缩感知的问题里没必要追求这么高1e-4 通常就够了因为再往下迭代下降的往往是噪声而不是信号。正则参数 λ 是 TV 和 FISTA 的核心旋钮。我的做法是先算一个参考量 scale norm(A*y, inf)然后令 λ α * scaleα 从 1e-3 扫到 1e-1取五个点画 PSNR 随 α 变化的曲线。曲线通常有一个宽顶宽顶越宽说明算法对参数越不敏感。如果曲线是尖峰说明这个算法在你的数据上不稳落地要谨慎。还有一个现象值得单独说残差降到很低但 PSNR 反而下降。这是过拟合噪声的典型信号早停本身就是一种正则化。做法是每迭代十轮算一次验证指标取指标最高那一轮的结果而不是最后一轮。这个技巧在真实系统里几乎是必用的因为真实噪声水平未知靠 λ 调不准靠早停更稳。4.4 并行与 GPU 的边界什么时候值得上parfor 只在「多重独立重建」的场景下有用比如扫采样率扫参数每个组合独立跑一次重建这时候把外层循环改成 parfor 能线性加速。单次重建内部没法用 parfor因为共轭梯度和 FISTA 都是严格串行的迭代强行拆分只会带来数据同步开销。GPU 的收益取决于矩阵规模。64×64 的时候把图案传到 GPU 的传输开销比计算本身还大实测反而更慢。128×128 开始有微弱收益256×256 以上收益明显但前提是图案一直留在显存里只在初始化时上传一次。判断标准很简单单次矩阵乘法的耗时要超过 50 毫秒GPU 才划算低于这个量级传输开销吃掉全部收益。另外 GPU 上的 single 精度和 CPU 上的行为不完全一致收敛曲线会有细微差别做严格对比实验时要么全用 CPU要么全用 GPU别混着来。5. 避坑与排查单像素成像跑不出结果时的五条血泪经验5.1 图糊成条纹或棋盘格先查图案同步与 ±1 差分现象重建结果不是噪点而是规则的横条纹、竖条纹或者棋盘格PSNR 在 10dB 上下换算法也没用。原因最常见的是测量矩阵的行序和实际投影图案的时序错位。图案按某个顺序加载读数按另一个顺序拼接A 的第 i 行和 y 的第 i 个元素对应的不是同一个图案。其次是 ±1 图案的差分没做对Hadamard 基底的理论取值是 1 和 -1而 DMD 只能显示 0 和 1正确做法是投影 1 和 -1 两帧做差有些人图省事只投一帧然后把 0 映射成 -1实际光强不是线性的矩阵里就混进了系统偏差。解决写一个自检脚本用第一个图案单独投影一次比对读数和 A 的第一行和已知图像的内积误差在千分之一以内才算同步。±1 差分必须实打实采两帧别用数学映射偷懒。5.2 指标好看图难看评价流程里的三个陷阱现象PSNR 报出来 32dB肉眼看重建图明显糊边缘全是锯齿。原因第一评价前对重建结果和真值分别做了 min-max 归一化等于把动态范围的误差抹掉了实际显示时又用原始范围当然对不上。第二重建结果里有负值和超过上限的值显示时被自动裁剪但算 PSNR 时没裁剪两种处理给出的数值不一致。第三只报了一个采样率下的数值恰好在那个采样率上算法表现好换个采样率就露馅。解决评价函数统一入口先对重建结果做 clip 到真值范围再用同一个函数算 PSNR 和 SSIM绝不各自归一化。指标必须报曲线至少五个采样率点。5.3 跑得慢、跑崩、结果飘性能和数值上的三条排查线现象一同一个算法换台机器跑耗时分差三倍以上。原因图案矩阵每次运行都重新生成随机矩阵的生成成本不低更麻烦的是没固定种子导致每次测的其实是不同的数据。解决图案预生成后存成 mat 文件测量时直接 load种子写进配置结构体。现象二迭代到某一步残差突然往上跳然后一路震荡不收敛。原因步长或正则参数和数据的量级不匹配最典型的是图案矩阵没做归一化AA 的最大特征值远大于 1FISTA 的步长 1/L 就取得过大。解决先跑幂迭代估出 L或者干脆把图案按列归一化让每个原子的二范数为 1。现象三64×64 跑得好好的改成 128×128 内存直接爆。原因显式存了测量矩阵或者中间变量用了 double 且没预分配。解决改成函数句柄算子中间变量预分配并统一用 single迭代里的临时变量用完就 clear。现象四重建结果整体偏暗或偏亮形状对但灰度不对。原因桶探测器的暗电流没扣或者归一化时用了不同批次的参考光强。解决每次采集前采一帧全暗和一帧全亮做参考测量向量先减暗再除以亮暗差这一步在算法之前做完别指望重建算法去补。6. 落到你手上的这份比较怎么用一周跑出自己的结论拿到一份算法比较最没用的做法是照抄结论最有用的做法是照抄实验框架然后把自己的硬件参数填进去。我给一个一周的节奏第一天把测量链路自检跑通用共轭梯度确认图案、时序、读数三者对齐PSNR 在无噪仿真下能到 40dB 以上第二天把四张测试图和五个采样率组合跑一遍纯共轭梯度建立基线第三天换 TV 和 FISTA各扫五个正则参数画出 PSNR 曲线第四天换 OMP扫稀疏度 K第五天加噪声重跑前三天的关键配置第六天看残差曲线和耗时分解把每个算法的算子调用次数乘出来第七天整理成自己的表按你的场景加权选型。下面这段是可以直接拿去改的主循环把正则项留成插槽换算法只改一行。n 64; N n^2; rate 0.25; m round(rate*N); rng(2024); P randn(m, N, single) / sqrt(m); xTrue double(phantom(n)); xTrue xTrue(:); y P * xTrue; % 无噪声基线先验证链路 % 共轭梯度解正规方程 AA x Ay用 CGNR 形式保证残差单调下降 x zeros(N, 1, double); r y - P * x; % 初始残差 z P * r; p z; rho z * z; maxit 200; tol 1e-5; for k 1:maxit w P * p; % 正向作用一次 alpha rho / (w * w); % 步长由 p 方向的曲率决定 x x alpha * p; r r - alpha * w; z P * r; % 反向作用一次 rhoNew z * z; if sqrt(rhoNew) tol * norm(P * y), break; end p z (rhoNew / rho) * p; % 共轭方向更新 rho rhoNew; end imgCG reshape(x, n, n); fprintf(iter %d, relres %.2e\n, k, sqrt(rhoNew) / norm(P * y));这段代码里alpha 的表达式是整个算法的核心它保证每一步都在当前共轭方向上取最优步长rho 的递推则保证方向之间共轭。停止准则用 sqrt(rhoNew) 和初始量的比值等价于相对残差比绝对残差稳。maxit 设 200 是保守值实际通常几十轮就满足容差把 k 打出来看一眼就知道收敛快不快。要换成 FISTA把这段循环体换成「梯度下降 小波软阈值 动量更新」三步正则参数从插槽里读其余结构不变。回过头看我在这个方向上最大的一次教训是把一次「算法 A 比算法 B 好 4dB」的结论写进报告三个月后被人复现不出来查到最后是两次实验的图案排序不同。从那以后我每次做算法对比都强制走一遍三件事固定种子并把种子写进结果文件名、把前置变量集中在配置结构体里、把残差曲线和指标曲线一起存档。这三件事花不了半小时但能省掉事后几个月的扯皮。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/26 20:21:02 · 阅读 1.2k 单像素成像Matlab算法选型：CG、OMP、TV与FISTA对比避坑
简介面向单像素成像与压缩感知重建学习者的Matlab算法对照包聚焦在同一实验框架下比较不同SPI重建思路的表现适合需要理解算法差异、复现实验或做课程与科研验证的本科生、研究生及工程人员。包内共10个文件以9个.m函数脚本和1个.txt说明文件为主压缩包约12KB代码按重建方法拆分另含统一演示入口与误差评估函数便于逐项调用和横向比较。内容覆盖差分鬼成像、梯度下降、共轭梯度下降、泊松最大似然、交替投影、稀疏表示压缩感知和全变分压缩感知等典型方法。已有2762人学习下载读者可据此快速搭建比较实验观察各算法在重建质量与收敛特性上的差异并在此基础上修改采样矩阵、参数或评价指标作为单像素成像入门复现与算法选型的参考。1. 单像素成像选算法这件事先别急着抄代码用 Matlab 做单像素成像算法比较最容易翻车的地方往往不是算法选错而是比较时变量没控住。我第一次把共轭梯度下降和一组全变分重建并排跑结论是「CG 不行」半年后回去查发现是测量矩阵行序和投影图案时序错开一位谁上都得糊。从那以后我给自己立了条规矩任何一次压缩感知重建的对比先把采样率、测量矩阵类型、噪声模型写死在脚本头部再谈谁好谁坏。这篇讲的是在单像素成像这条链路上共轭梯度下降、正交匹配追踪、全变分和 FISTA 各自适合什么采样率与图案Matlab 实现里哪些细节会让排名反转以及怎么用一周时间跑出属于你自己硬件条件的结论。2. 从桶探测到图像单像素成像的重建链路拆给你看单像素成像的物理过程很朴素用空间光调制器把一组预先设计好的图案打到目标上每个图案对应一次桶探测器的总光强读数最后拿「图案矩阵 读数向量」反解出图像。真正决定成败的环节不在探测器而在反解这一步用什么算法、带什么先验。很多人在这一步直接抄一段压缩感知的重建代码跑通了就算完事结果换一组图案、换一个采样率指标就掉得没法看——不是代码写错了是链路里每个环节的耦合关系没搞清楚。2.1 一条完整的重建链路里算法到底插在哪把链路拆开看它其实只有五段。第一段是图案生成用随机二值、高斯随机、Hadamard 排序或者傅里叶基底生成 m 个 n×n 的图案第二段是调制与投影图案被 DMD 或 LCD 加载目标经过调制后总光强落在单点探测器上第三段是采集探测器输出模拟电压ADC 量化成一个标量读数第四段是把 m 个读数拼成向量 y把 m 个图案拉直拼成测量矩阵 A得到方程 y A x n第五段才是重建算法登场从欠定的 y 和 A 里恢复 x。算法只负责第五段但它能做的事情被前四段死死限制住。A 的质量决定了解的上限图案如果不是理想的 ±1 正交基底或者 DMD 的二值化把灰度图案截断失真矩阵里就混进了系统性误差再好的算法也只是在拟合误差。探测器噪声模型同样如此光子计数低的时候是泊松噪声读出电路主导的时候近似高斯噪声两种噪声下同一套正则参数的取值能差一个数量级。所以「算法比较」这个说法本身要拆细。我一般会把它拆成三件套来比先验是什么无先验、稀疏先验、全变分先验、算子是什么显式矩阵、函数句柄、快速变换、求解器是什么共轭梯度、贪婪追踪、近端梯度、增广拉格朗日。同一个 TV 先验用 ADMM 解和用原始-对偶解收敛速度可以差两三倍同一个共轭梯度作用在显式矩阵上和作用在快速变换算子上内存占用差一个量级。只看算法名字做对比等于只看发动机型号不看变速箱和车重。2.2 比算法之前先把三类前置变量钉死我见过太多对比实验两组结果的差异其实来自测量矩阵不同却被归因到算法头上。要让比较有意义下面三类变量必须先固定且固定在脚本最顶部用一个结构体统一管理换值只改一处。前置变量常见取值对算法排名的影响我的固定建议测量矩阵高斯随机、伯努利、Hadamard 排序、小波、傅里叶基底相干性低的矩阵让 OMP 类贪婪算法占便宜结构化矩阵让快速算子占便宜基线统一用归一化高斯随机行向量除以 sqrt(m)采样率5% 到 50%按 N 的百分比取整低于 15% 时无先验的共轭梯度基本失效TV 和 L1 才有优势基线取 25%再扫 5/10/15/30/50 五个点噪声模型无噪仿真、加性高斯、泊松光子噪声泊松噪声下 L1 类算法对低光子区域的重建偏差明显更大先跑无噪验证链路再加 30dB 高斯最后上泊松这三类固定之后才轮到图像内容。我习惯准备四张测试图一张分段常数的块状图、一张自然纹理图、一张二值稀疏点阵、一张带渐变背景的目标。分段常数图用来暴露 TV 的优势和阶梯效应自然纹理图用来暴露 TV 的细节丢失二值点阵用来验证 OMP 在强稀疏下的表现渐变图用来检验各算法在低频分量上的还原能力。四张图跑完算法的大致画像就出来了比只跑一张 Lena 可靠得多。还有一件事容易漏图案的排序。Hadamard 基底如果按自然序取前 m 行低频分量集中在前面的行里重建结果会严重偏向低频看起来「平滑但糊」按随机序或者按 sequency 排序取频率分布均匀得多。这个差别会让同一个算法在同样的采样率下 PSNR 差 3 到 5dB足以让排名颠倒。2.3 判断一个算法好不好我只看四个指标指标不能只看重建质量单像素成像的落地场景里速度和可调性往往比多 1dB 更重要。我固定看四项。第一项是重建质量但必须看成曲线而不是单点。固定其他变量扫采样率画 PSNR 和 SSIM 随采样率变化的曲线。两条曲线的交叉点比单点数值有信息量得多很多算法在 10% 采样率下差距巨大到 40% 以上就收敛到差不多的水平。第二项是收敛成本不能只报总耗时要拆成迭代次数乘以单次迭代的矩阵向量乘次数再乘以单次乘法的耗时。共轭梯度每次迭代两次算子作用FISTA 两次TV 的 ADMM 每次迭代四到六次这个系数在比较时经常被忽略。第三项是参数敏感度也就是正则参数在多大范围内波动时PSNR 下降不超过 1dB。敏感度高的算法在真实系统里没法用因为真实系统的噪声水平每次都在变你不可能每次都手工调参。第四项是内存与可扩展性能不能从 64×64 撑到 128×128 甚至 256×256这一项直接决定算法能不能上真实样机。四项指标的权重取决于你的场景。做离线验证的质量权重最高做实时成像的收敛成本优先做产品化封装的参数敏感度是生死线。3. 四类重建算法逐项对比CG、OMP、TV、FISTA 谁更适合你3.1 共轭梯度下降快、稳但别指望它自带稀疏先验共轭梯度下降在单像素成像里最常见的用法是解最小二乘问题 min ||A x - y||²。它不直接对正规方程做显式求逆而是在 Krylov 子空间里构造一组共轭方向理论上 n 步内收敛实际因为舍入误差和条件数限制通常几十步就够。Matlab 里对应的入口是 pcg 和 lsqr前者解对称正定的正规方程后者直接处理长方矩阵数值稳定性更好。它的最大优点是干净没有正则参数只有最大迭代次数和容差两个旋钮调起来没什么玄学。每次迭代的算子作用次数是两次一次正向 A p一次反向 A r内存占用是 O(N)不需要存中间的大矩阵。在采样率比较高、测量矩阵条件数可控的场景里它是极好的基线和快速预览工具。问题也很明确最小二乘本身不含任何稀疏先验。当 m 远小于 n 时方程欠定解不唯一共轭梯度会收敛到最小范数解附近噪声被完整地搬进重建图像。表现就是图整体偏亮偏糊高频全是噪点采样率越低越明显。我做过一组对照25% 采样率下纯共轭梯度的 PSNR 比带小波软阈值的版本低 6dB 左右到 10% 采样率直接不可用。真要用共轭梯度做压缩感知重建得把它嵌进迭代软阈值的框架里每一轮先做一次共轭梯度下降逼近数据保真项再在小波域做一次软阈值收缩交替进行。这样既保留了共轭梯度的收敛速度又把稀疏先验加了进去代价是多了阈值参数和内外两层循环调参复杂度上来了。我的经验是如果只是想在写正式算法之前先验证测量链路通不通纯共轭梯度足够了如果要拿它当最终方案必须加先验。3.2 正交匹配追踪稀疏度高时漂亮字典一相干就崩正交匹配追踪属于贪婪类算法逻辑直白每一轮从字典里挑出与当前残差最相关的那一列把它的系数用最小二乘更新然后从残差里减掉这一列的贡献重复 K 轮K 就是预设的稀疏度。它的计算成本主要在相关计算和最小二乘更新上K 通常远小于 n所以总耗时在四类算法里最低。在单像素成像里OMP 有一个非常合适的场景目标是强稀疏的比如二值点阵、稀疏荧光标记、远距离小目标。这时候如果测量矩阵本身就是稀疏基或者目标在某个变换域里只有几十个非零系数OMP 在 5% 到 10% 的采样率下就能给出干净的结果速度快到可以边采集边重建。但自然图像是它的死穴。自然图像在小波域里是近似稀疏而不是严格稀疏系数呈幂律衰减尾部一大堆小系数。更麻烦的是单像素成像的等效字典是测量矩阵和小波基的乘积随机矩阵和小波基之间的相干性虽然不算最差但也足够让 OMP 在迭代早期选错原子。一旦选错后面每一轮都在错误的支撑集上做最小二乘误差越滚越大最后的结果是一堆分散的亮点加一块糊掉的背景。稀疏度 K 的设定也是个大坑。设小了重建欠拟合细节全丢设大了等于把噪声也当成信号选进来。我一般先用小波系数的能量累积曲线估一个 K 的量级再在附近扫三五个值取 PSNR 最高的那个。这个做法在仿真里可行在真实系统里不可行因为你事先不知道目标的稀疏度。这也是 OMP 在工程落地里用得比论文里少的原因。3.3 全变分与 FISTA两类正则化求解器的分工全变分正则化解决的是另一类图像分段常数、块状目标、二值图案、边缘锐利的工业检测对象。它的目标函数是数据保真项加上全变分项全变分项惩罚相邻像素差的 L1 范数等价于鼓励图像分块平滑。求解上常用增广拉格朗日或者 ADMM把问题拆成几个可以交替求解的子问题TVAL3 是这类思路里比较常见的实现叫法。TV 在块状图上表现非常好采样率 10% 就能把方块的边缘还原得干干净净PSNR 能比共轭梯度高出 10dB 以上。代价是阶梯效应渐变区域会被还原成一层层的台阶纹理细节被抹平人脸皮肤这类缓变区域会出现明显的块状伪影。这是 L1 类先验的通病不是实现问题调参也救不回来。正则参数 λ 的敏感度也比较高λ 大了过度平滑λ 小了噪点全回来可用的窗口通常在半个数量级以内。FISTA 走的是另一条路数据保真项加小波域 L1 范数用近端梯度下降求解软阈值算子做收缩再用 Nesterov 动量做加速收敛速率从 O(1/k) 提升到 O(1/k²)。它更适合自然图像这类在小波域近似稀疏的目标纹理保留比 TV 好但边缘锐利度不如 TV。步长取 1/LL 是 A 的 Lipschitz 常数也就是 AA 的最大特征值实际用幂迭代估十几轮就够准。步长取大了会发散取小了收敛慢这个参数不能拍脑袋。这两类算法在实现上有一个共同点容易被忽略都不需要显式存 A只要能用函数句柄算出 A x 和 A y 就行。这意味着它们可以直接对接快速变换把 O(mN) 的矩阵乘法降到 O(N log N)在 256×256 以上的分辨率上是唯一可行的路子。3.4 逐项对比表与选型结论把四类算法放在同一张表里差异就清楚了。下面的表是在 25% 采样率、归一化高斯测量矩阵、无噪条件下测的换条件排名会变但相对关系大体稳定。算法先验关键参数每次迭代算子调用适用采样率优势翻车场景Matlab 入口共轭梯度无可外挂软阈值maxit、tol240% 以上无参数、速度快、内存小低采样率下噪声放大pcg、lsqrOMP稀疏度 KK1 次相关 1 次最小二乘5% 到 15%强稀疏目标极快、结果可解释自然图像、字典相干性高自写主循环TVADMM全变分λ、μ、maxit4 到 610% 到 30%块状图、二值目标表现最好纹理丢失、阶梯效应自写 ADMMFISTA小波 L1λ、步长、maxit210% 到 30%自然图像均衡、收敛有加速λ 敏感、步长估错会发散自写主循环结论按人群分。如果你在搭第一套验证链路先用共轭梯度跑通全流程确认图案、时序、读数都对得上再换算法这一天的投入能省掉后面一周的排查。如果你的目标明确是稀疏点阵或者二值图案采样率还想压到 10% 以下OMP 和 TV 都值得试前者快后者稳块状目标优先 TV。如果你做的是自然场景成像采样率在 20% 左右FISTA 是综合最平衡的选择纹理和边缘都能兼顾。如果你的系统要求实时别在算法上死磕先把测量矩阵换成结构化快速变换把单次算子作用从 O(mN) 降到 O(N log N)这比换算法带来的提速大得多。4. Matlab 实现层的性能拆解同一算法为什么跑出两种结果4.1 测量矩阵存 sparse 还是 full先把内存算清楚单像素成像的测量矩阵尺寸是 m×Nm 是测量次数N 是像素总数。这个矩阵是稠密的稀疏存储没有意义所以问题就变成 full 存储能不能扛住。按 double 类型每个元素 8 字节算不同分辨率的账很好算。图像分辨率像素数 N25% 采样 mfull double 内存full single 内存建议64×644096102432 MB16 MBfull 存储随便跑128×128163844096512 MB256 MB用 single或改算子256×25665536163848 GB4 GB必须用函数句柄512×51226214465536128 GB64 GB结构化快速变换 分块这张表是我踩过坑之后贴在显示器边上的。64×64 的时候一切都很美好代码随手写成 A randn(m, N)跑得飞快然后直接把分辨率改成 128×128Matlab 卡死三分钟之后报内存不足。更隐蔽的情况是内存够但速度崩了单个矩阵乘法涉及 4096×16384 次浮点运算每轮迭代做两次几百轮迭代下来时间全花在等内存带宽上。正确的做法是从一开始就把 A 写成算子而不是矩阵。这不只是为了省内存更是为了让同一份重建代码能在不同分辨率下无缝切换。转 single 也要谨慎单精度在矩阵乘法里会累积误差共轭梯度这类依赖内积精度的算法迭代几十轮之后残差就可能不再下降表现为「跑到一半卡住」。4.2 算子写成函数句柄别写成显式矩阵Matlab 的 lsqr、pcg、bicgstab 都支持传函数句柄签名是 afun(x, trans_flag)其中 trans_flag 取 notransp 或 transp。把正向和反向作用包在一个子函数里外部用匿名函数捕获图案参数就能做到图案常驻内存、矩阵不落地。% 用一个子函数承载正反向算子作用外部只传句柄给 lsqr/pcg function y opA(P, x, transFlag) if strcmp(transFlag, notransp) y P * x; % 正向图像向量 - 测量向量 else y P * x; % 反向测量向量 - 图像向量 end end % 主脚本片段 n 64; N n^2; rate 0.25; m round(rate*N); rng(2024); % 固定种子保证每次比较用同一组图案 P randn(m, N, single) / sqrt(m); % 归一化避免内积量级失控 y P * double(phantom(n)(:)); % 无噪仿真测量 Afun (x, t) opA(P, x, t); % 句柄捕获 P不复制矩阵 [x, flag, relres, iter] lsqr(Afun, y, 1e-6, 200, [], [], zeros(N,1)); img reshape(x, n, n);这段代码里几个参数值得说明。rng 的种子必须固定否则每次跑出来的图案都不一样两个算法的比较就没有意义。图案矩阵除以 sqrt(m) 是归一化让 AA 的对角元素接近 1这样共轭梯度的收敛行为和正则参数的取值在不同采样率之间才可比。lsqr 的第三个参数是容差第四个是最大迭代次数第五到第七个是阻尼和初始值这里阻尼取空表示不做 Tikhonov 正则初始值全零。返回的 flag 为 1 表示达到容差为 2 表示达到最大迭代次数为 3 表示残差不再下降——看到 flag 为 3 就要回头查矩阵是不是秩亏或者量级失衡。如果图案是 Hadamard 基底根本不用生成 m×N 的矩阵直接调用快速沃尔什-哈达玛变换把向量变换后取前 m 个系数即可复杂度从 O(mN) 降到 O(N log N)。这个改动在 256×256 以上是数量级的差别值得单独花半天实现。4.3 停止准则与正则参数最容易被忽略的调参旋钮共轭梯度的停止准则用相对残差别用绝对残差。abs(norm(A*r)) tol 这种写法在图像动态范围变化时会完全失效换成 norm(r)/norm(y) tol 才稳定。tol 取 1e-6 是常见起点但压缩感知的问题里没必要追求这么高1e-4 通常就够了因为再往下迭代下降的往往是噪声而不是信号。正则参数 λ 是 TV 和 FISTA 的核心旋钮。我的做法是先算一个参考量 scale norm(A*y, inf)然后令 λ α * scaleα 从 1e-3 扫到 1e-1取五个点画 PSNR 随 α 变化的曲线。曲线通常有一个宽顶宽顶越宽说明算法对参数越不敏感。如果曲线是尖峰说明这个算法在你的数据上不稳落地要谨慎。还有一个现象值得单独说残差降到很低但 PSNR 反而下降。这是过拟合噪声的典型信号早停本身就是一种正则化。做法是每迭代十轮算一次验证指标取指标最高那一轮的结果而不是最后一轮。这个技巧在真实系统里几乎是必用的因为真实噪声水平未知靠 λ 调不准靠早停更稳。4.4 并行与 GPU 的边界什么时候值得上parfor 只在「多重独立重建」的场景下有用比如扫采样率扫参数每个组合独立跑一次重建这时候把外层循环改成 parfor 能线性加速。单次重建内部没法用 parfor因为共轭梯度和 FISTA 都是严格串行的迭代强行拆分只会带来数据同步开销。GPU 的收益取决于矩阵规模。64×64 的时候把图案传到 GPU 的传输开销比计算本身还大实测反而更慢。128×128 开始有微弱收益256×256 以上收益明显但前提是图案一直留在显存里只在初始化时上传一次。判断标准很简单单次矩阵乘法的耗时要超过 50 毫秒GPU 才划算低于这个量级传输开销吃掉全部收益。另外 GPU 上的 single 精度和 CPU 上的行为不完全一致收敛曲线会有细微差别做严格对比实验时要么全用 CPU要么全用 GPU别混着来。5. 避坑与排查单像素成像跑不出结果时的五条血泪经验5.1 图糊成条纹或棋盘格先查图案同步与 ±1 差分现象重建结果不是噪点而是规则的横条纹、竖条纹或者棋盘格PSNR 在 10dB 上下换算法也没用。原因最常见的是测量矩阵的行序和实际投影图案的时序错位。图案按某个顺序加载读数按另一个顺序拼接A 的第 i 行和 y 的第 i 个元素对应的不是同一个图案。其次是 ±1 图案的差分没做对Hadamard 基底的理论取值是 1 和 -1而 DMD 只能显示 0 和 1正确做法是投影 1 和 -1 两帧做差有些人图省事只投一帧然后把 0 映射成 -1实际光强不是线性的矩阵里就混进了系统偏差。解决写一个自检脚本用第一个图案单独投影一次比对读数和 A 的第一行和已知图像的内积误差在千分之一以内才算同步。±1 差分必须实打实采两帧别用数学映射偷懒。5.2 指标好看图难看评价流程里的三个陷阱现象PSNR 报出来 32dB肉眼看重建图明显糊边缘全是锯齿。原因第一评价前对重建结果和真值分别做了 min-max 归一化等于把动态范围的误差抹掉了实际显示时又用原始范围当然对不上。第二重建结果里有负值和超过上限的值显示时被自动裁剪但算 PSNR 时没裁剪两种处理给出的数值不一致。第三只报了一个采样率下的数值恰好在那个采样率上算法表现好换个采样率就露馅。解决评价函数统一入口先对重建结果做 clip 到真值范围再用同一个函数算 PSNR 和 SSIM绝不各自归一化。指标必须报曲线至少五个采样率点。5.3 跑得慢、跑崩、结果飘性能和数值上的三条排查线现象一同一个算法换台机器跑耗时分差三倍以上。原因图案矩阵每次运行都重新生成随机矩阵的生成成本不低更麻烦的是没固定种子导致每次测的其实是不同的数据。解决图案预生成后存成 mat 文件测量时直接 load种子写进配置结构体。现象二迭代到某一步残差突然往上跳然后一路震荡不收敛。原因步长或正则参数和数据的量级不匹配最典型的是图案矩阵没做归一化AA 的最大特征值远大于 1FISTA 的步长 1/L 就取得过大。解决先跑幂迭代估出 L或者干脆把图案按列归一化让每个原子的二范数为 1。现象三64×64 跑得好好的改成 128×128 内存直接爆。原因显式存了测量矩阵或者中间变量用了 double 且没预分配。解决改成函数句柄算子中间变量预分配并统一用 single迭代里的临时变量用完就 clear。现象四重建结果整体偏暗或偏亮形状对但灰度不对。原因桶探测器的暗电流没扣或者归一化时用了不同批次的参考光强。解决每次采集前采一帧全暗和一帧全亮做参考测量向量先减暗再除以亮暗差这一步在算法之前做完别指望重建算法去补。6. 落到你手上的这份比较怎么用一周跑出自己的结论拿到一份算法比较最没用的做法是照抄结论最有用的做法是照抄实验框架然后把自己的硬件参数填进去。我给一个一周的节奏第一天把测量链路自检跑通用共轭梯度确认图案、时序、读数三者对齐PSNR 在无噪仿真下能到 40dB 以上第二天把四张测试图和五个采样率组合跑一遍纯共轭梯度建立基线第三天换 TV 和 FISTA各扫五个正则参数画出 PSNR 曲线第四天换 OMP扫稀疏度 K第五天加噪声重跑前三天的关键配置第六天看残差曲线和耗时分解把每个算法的算子调用次数乘出来第七天整理成自己的表按你的场景加权选型。下面这段是可以直接拿去改的主循环把正则项留成插槽换算法只改一行。n 64; N n^2; rate 0.25; m round(rate*N); rng(2024); P randn(m, N, single) / sqrt(m); xTrue double(phantom(n)); xTrue xTrue(:); y P * xTrue; % 无噪声基线先验证链路 % 共轭梯度解正规方程 AA x Ay用 CGNR 形式保证残差单调下降 x zeros(N, 1, double); r y - P * x; % 初始残差 z P * r; p z; rho z * z; maxit 200; tol 1e-5; for k 1:maxit w P * p; % 正向作用一次 alpha rho / (w * w); % 步长由 p 方向的曲率决定 x x alpha * p; r r - alpha * w; z P * r; % 反向作用一次 rhoNew z * z; if sqrt(rhoNew) tol * norm(P * y), break; end p z (rhoNew / rho) * p; % 共轭方向更新 rho rhoNew; end imgCG reshape(x, n, n); fprintf(iter %d, relres %.2e\n, k, sqrt(rhoNew) / norm(P * y));这段代码里alpha 的表达式是整个算法的核心它保证每一步都在当前共轭方向上取最优步长rho 的递推则保证方向之间共轭。停止准则用 sqrt(rhoNew) 和初始量的比值等价于相对残差比绝对残差稳。maxit 设 200 是保守值实际通常几十轮就满足容差把 k 打出来看一眼就知道收敛快不快。要换成 FISTA把这段循环体换成「梯度下降 小波软阈值 动量更新」三步正则参数从插槽里读其余结构不变。回过头看我在这个方向上最大的一次教训是把一次「算法 A 比算法 B 好 4dB」的结论写进报告三个月后被人复现不出来查到最后是两次实验的图案排序不同。从那以后我每次做算法对比都强制走一遍三件事固定种子并把种子写进结果文件名、把前置变量集中在配置结构体里、把残差曲线和指标曲线一起存档。这三件事花不了半小时但能省掉事后几个月的扯皮。希望帮到你。本文还有配套的精品资源点击获取