简介：凸集投影法（POCS）超分辨重建的MATLAB实现，面向具备一定图像处理与信号恢复基础、希望理解迭代投影思路的学习者，可用于复现低分辨率图像提升至高分辨率图像的关键流程。压缩包内共1个文件，为1个m脚本，体积约1KB，轻量紧凑，核心代码便于逐行阅读和调试。已有3872人学习，说明该主题在相关课程与实践中关注度较高。代码以初始化、投影操作、迭代更新和停止条件为主线，体现非负性、能量或频谱等凸集约束下的逐步逼近过程，也便于进一步结合频域扩展、正则化策略理解超分辨重建的实现逻辑。读者可据此梳理POCS算法的代码结构、迭代流程与约束投影方式，并在此基础上开展MATLAB实验、参数调整和算法对照学习。

1. 从一张糊掉的监控截图说起：POCS 超分辨重建到底在补什么

监控画面里那辆车只有三十几像素宽，放大两倍后车牌是一整片模糊。双三次插值放大，边缘发虚；想上深度学习模型，手上没有配对训练数据，显卡也不够折腾。凸集投影法（POCS）超分辨重建算法 MATLAB 实现这条路正好卡在中间：它不需要训练，只靠多帧之间的亚像素位移提供的互补信息，一轮一轮往凸约束集里投影，几十行 MATLAB 就能跑起来，中间每一步都看得见。

它适合谁？手里有连拍或多帧序列、能估出亚像素位移、又不想引入模型依赖的人。不适合单帧——一帧没有互补信息，POCS 会退化成带约束的锐化，效果和普通反卷积拉不开差距。

2. 观测模型与凸集：把「超分辨」翻译成三个可投影的约束

POCS 的骨架其实就两句话：先写清楚低分辨图是怎么从高分辨图退化出来的，再把「高分辨图应该满足的常识」写成凸集，然后反复往这些凸集上投影。退化模型写成公式是y_k = D·H·F_k·x + n_k，其中x是待求的高分辨图像，F_k是第 k 帧相对参考帧的亚像素平移，H是光学点扩散函数（PSF）卷积，D是降采样，n_k是噪声。四个算子各有物理含义：F_k来自相机或物体的微小移动，H来自镜头与传感器的模糊，D就是你想要的放大倍数。

2.1 观测模型里每一项对应的物理量

把公式拆开对照实物，调参时才知道该动哪个旋钮。F_k的单位是「高分辨像素」，典型值在 ±r 像素以内（r 是放大倍数），超出这个范围说明帧间位移已经大到可以用整像素配准解决，反而更好处理。H一般用高斯核近似，标准差 0.8～1.5 之间比较常见，它同时承担了镜头模糊和传感器像元积分两部分，所以不要照搬光学设计值，按实测边缘的过渡宽度反推更靠谱。D是理想的隔点采样，真实传感器还有积分窗，严格讲应该是先做 r×r 的盒式滤波再抽样，这一点在 r ≥ 3 时影响明显，r = 2 时可以先忽略。

噪声那一项决定了后面数据保真约束的容差 δ。如果用仿真数据，噪声标准差是你自己加的，直接取 2～3 倍即可；用实拍数据，就在平坦区域取一块算标准差，别凭感觉填 0。

2.2 三个必须实现的凸集及其几何含义

这三个集合都是凸的，交集也是凸的，所以交替投影才有收敛保证。幅度集最便宜，几乎不耗时；数据保真集是主力，信息全靠它进；梯度集是保险丝，专门压振铃和噪点。工程上再加一个可选的参考帧集——把某一帧当作最可信的观测，直接约束它降质后与原始帧的差异——在帧数少于 6 张时能明显稳住结果。

2.3 松弛迭代：为什么不直接用严格投影

严格交替投影在理论上收敛，实际跑起来抖得厉害，尤其数据保真集的投影方向在相邻两轮之间会来回摆。常见做法是加松弛因子 λ：x^{n+1} = x^n + λ·(P(x^n) − x^n)。λ = 1 是严格投影，λ 在 1.0～1.6 之间通常收敛最快，超过 2 一定发散。这个结论不是玄学，它对应的是算子范数约束：只要 λ 落在 (0, 2/K) 区间内迭代就是压缩的，K 是复合算子D·H·F的范数平方。

% pocs_sr.m —— 松弛形式的交替投影主循环 function [X, hist] = pocs_sr(LR, cfg, X0) X = X0; K = estimate_opnorm(cfg); % 复合算子范数平方的估计值 hist = zeros(cfg.maxIter, 1); for it = 1:cfg.maxIter Xp = X; Xp = proj_data_fidelity(Xp, LR, cfg, K); % C_D：多帧残差反投影 Xp = proj_gradient(Xp, cfg.tvT, cfg.beta); % C_T：压制过冲梯度 Xp = proj_amplitude(Xp, cfg.lo, cfg.hi); % C_A：最后截断，保证合法值域 X = X + cfg.lambda * (Xp - X); % 松弛更新 hist(it) = norm(X(:) - Xp(:)) / numel(X); % 投影残差，用于判断是否早停 end end

这段循环里，三个投影的先后顺序在理论上不影响收敛点，实践上把幅度投影放最后更稳——前两个投影都可能把值推出 [0,1]，最后截一刀成本最低。hist记录的是本轮迭代「被投影改动了多少」，它单调下降到某个平台期就说明约束基本都满足了，继续迭代只剩噪声在动。这一步千万别省，后面判断早停全靠它。

estimate_opnorm用幂迭代在小尺寸上估一个标量就行，尺寸对算子范数影响不大，但必须估：手工填 λ 又不定 K，是新手最容易翻的车。

function K = estimate_opnorm(cfg) % 幂迭代估 ||D·H·F||^2，忽略边界效应 v = randn(32); v = v / norm(v(:)); for t = 1:20 w = blur2(v, cfg); w = downsample2(w, cfg.r); w = upsample2(w, cfg.r); w = blur2(w, cfg); nrm = norm(w(:)); v = w / nrm; end K = nrm; end

幂迭代收敛很快，20 次足够。nrm就是主特征值的近似，也就是 K。注意这里用的upsample2是零插值（严格伴随），不是双线性插值，两者数值上差一个 r² 的因子，混用会让 K 偏大或偏小，λ 的可用区间跟着漂移。

3. MATLAB 工程骨架：五个函数搭出可跑通的迭代闭环

算法本身不到 100 行，但能不能一次跑通，取决于外围那几件事：仿真数据生成、上采样初始化、算子的小心实现。我习惯先把观测模型单独写成可复现的仿真模块，用已知真值跑通，再换真实数据。这样出问题时能立刻分清是算法错了还是数据/配准错了——这一步省不得。

3.1 目录结构与函数职责

建议按这个结构铺，后面换数据类型、换 PSF 都不用动主循环：

pocs_sr/ ├── main_pocs_sr.m % 入口：配置 + 调用 + 出图 ├── sim_observation.m % 仿真观测：模糊-平移-降采样-加噪 ├── pocs_sr.m % 迭代主循环 ├── proj_data_fidelity.m % 数据保真投影 ├── proj_gradient.m % 梯度约束投影 ├── proj_amplitude.m % 幅度约束投影 ├── shift_subpixel.m % 亚像素平移（傅里叶/插值两种） ├── blur2.m % 高斯 PSF 卷积 ├── downsample2.m % 隔点降采样 ├── upsample2.m % 零插值上采样（D 的伴随） ├── estimate_opnorm.m % 幂迭代估算子范数 └── eval_metrics.m % PSNR / SSIM

把算子拆成独立文件的好处是每一个都能单测：先测shift_subpixel平移 1 个整像素和circshift对不对得上，再测平移 0.5 像素后反平移能不能还原，这两个过了，后面九成的怪现象都能排除掉。我自己就吃过亏——傅里叶平移的频率轴排错了一半，结果整幅图在迭代中缓慢漂移，查了两天才定位到那三行。

3.2 仿真观测生成：把每一步都写成可关掉的开关

function [LR, GT, cfg] = sim_observation(cfg) % 退化链：模糊 -> 亚像素平移 -> 降采样 -> 加噪 GT = double(imread('cameraman.tif')) / 255; GT = imresize(GT, [256 256]); [H, W] = size(GT); cfg.sy = zeros(1, cfg.L); cfg.sx = zeros(1, cfg.L); LR = zeros(H/cfg.r, W/cfg.r, cfg.L); for k = 1:cfg.L cfg.sy(k) = (rand - 0.5) * 2 * cfg.r; % 高分辨像素为单位的亚像素位移 cfg.sx(k) = (rand - 0.5) * 2 * cfg.r; Z = blur2(GT, cfg); % 光学模糊 Z = shift_subpixel(Z, cfg.sy(k), cfg.sx(k), ... cfg.shiftMode); % 亚像素平移 Z = downsample2(Z, cfg.r); % 降采样 Z = Z + cfg.noiseSig * randn(size(Z)); % 高斯噪声 LR(:,:,k) = Z; end % 第一帧作为参考帧，位移置零，避免后面重复补偿 cfg.sy(1) = 0; cfg.sx(1) = 0; end

位移量把rand-0.5乘2*r，是为了保证位移落在 ±r 个高分辨像素以内、且不全挤在整数附近——这是验证算法能否真正利用亚像素信息的关键。如果把位移量化到整像素，POCS 的优势立刻消失，PSNR 提升会掉到 0.3 dB 以内，很多人以为算法不行，其实是仿真设计错了。

cfg.shiftMode留了'fourier'和'interp'两个选项：前者精确但有循环边界，后者用interp2的'spline'，边界干净但会引入轻微模糊。两边都保留，可以在排查时交叉验证。

function Z = shift_subpixel(X, dy, dx, mode) % dy/dx 以高分辨像素为单位，正数表示向下、向右 [M, N] = size(X); if strcmp(mode, 'fourier') u = (-floor(N/2):ceil(N/2)-1) / N; % 列方向归一化频率 v = (-floor(M/2):ceil(M/2)-1) / M; % 行方向 [U, V] = meshgrid(u, v); F = fftshift(fft2(X)) .* exp(-1i*2*pi*(V*dy + U*dx)); Z = real(ifft2(ifftshift(F))); else [xx, yy] = meshgrid(1:N, 1:M); Z = interp2(xx, yy, X, xx + dx, yy + dy, 'spline', 0); end end

傅里叶分支里fftshift之后频率向量按-floor(N/2) : ceil(N/2)-1排，偶数长度时中间那个 Nyquist 点落在负频率一侧，这个约定必须和ifftshift配对，否则相位斜着走、图像整体平移量只有一半。'interp'分支的两个 meshgrid 注意是(xx, yy)对应(列, 行)，写反了结果是转置的，症状是重建图看起来「像旋转了 90 度」。

3.3 初始估计与算子归一化

初始值只影响收敛速度，不影响最终收敛点（凸集交非空的前提下）。用双三次上采样第一帧最省事，比全零初始化少迭代十几次。

function X0 = init_upsample(y, r) X0 = imresize(y, r, 'bicubic'); X0 = min(max(X0, 0), 1); % 先落进幅度集，避免首轮投影跳变 end

另外两个伴随算子写清楚，后面才能谈归一化：

function Y = downsample2(X, r) Y = X(1:r:end, 1:r:end); % 理想隔点采样 end function Y = upsample2(X, r) % D 的严格伴随：零插值，不乘 r^2 [M, N] = size(X); Y = zeros(M*r, N*r); Y(1:r:end, 1:r:end) = X; end

严格伴随就是零插值、不缩放。有些实现里乘了 r²，那是为了补偿零插值带来的能量损失，属于人为归一化，会让 K 的估计值差 r⁴，λ 就得跟着改。两种写法都能用，但在同一个工程里必须统一，混着写就会出现「明明 λ 设了 1.2 还是发散」这种情况。

3.4 入口脚本：一张图看清所有参数

% main_pocs_sr.m clear; clc; close all; rng(20240501); % 固定种子，保证每次实验可比 cfg.r = 2; % 放大倍数 cfg.L = 12; % 低分辨帧数 cfg.psfSig = 1.2; % 高斯 PSF 标准差 cfg.maxIter = 60; % 迭代上限 cfg.lambda = 1.4; % 松弛因子 cfg.tvT = 0.06; % 梯度约束阈值 cfg.beta = 0.25; % 梯度投影步长 cfg.lo = 0; % 幅度下界 cfg.hi = 1; % 幅度上界 cfg.noiseSig = 0.005; % 仿真噪声标准差 cfg.shiftMode = 'fourier'; [LR, GT, cfg] = sim_observation(cfg); X0 = init_upsample(LR(:,:,1), cfg.r); [X, hist] = pocs_sr(LR, cfg, X0); eval_metrics(X, GT);

rng固定种子这件事别嫌麻烦。调参时每次都换一组随机位移，PSNR 波动能有 1 dB，你根本分不清是参数变好了还是运气好。固定种子之后，改一个参数看一次结果，因果链才清晰。

4. 投影算子逐个写：幅度、残差、梯度三类的实现与参数

三个投影里，幅度投影几乎没有实现成本，梯度投影最容易写歪，数据保真投影决定成败。下面按「信息量从大到小」的顺序排。

4.1 数据保真投影：残差反投影与逐帧累加

这是唯一一个把观测信息注入进来的算子。思路直白：当前估计降质后与真实低分辨帧的差，沿原路反投回去。

function X = proj_data_fidelity(X, LR, cfg, K) % 逐帧做一次残差反投影，K 为 ||D·H·F||^2 的估计 for k = 1:size(LR, 3) y = LR(:,:,k); % 前向：把当前估计退化到第 k 帧的观测空间 yh = blur2(shift_subpixel(X, cfg.sy(k), cfg.sx(k), cfg.shiftMode), cfg); yh = downsample2(yh, cfg.r); r = y - yh; % 低分辨残差 % 反向：上采样 -> 反卷积（同核）-> 反平移 g = upsample2(r, cfg.r); g = blur2(g, cfg); g = shift_subpixel(g, -cfg.sy(k), -cfg.sx(k), cfg.shiftMode); X = X + cfg.lambda * g / K; % 松弛 + 归一化 end end

逐句说：yh是「用当前猜的高分辨图，按照第 k 帧的位移和模糊重新拍一张低分辨图」，r是这张重拍图和真实图的差。g是把差沿原路送回高分辨网格——上采样对应 D 的伴随，blur2对应 H 的伴随（高斯核是对称的，正反一样），反平移对应 F 的伴随。除以 K 是保证更新步长落在稳定区内。

帧数多的时候不必每轮都过所有帧。我一般每轮随机抽 3～5 帧，洗牌顺序，跑 60 轮下来每帧平均被用到 15 次以上，效果和全量遍历差不多，时间省一半。

提示：cfg.sy/cfg.sx里第一帧必须是 0。如果参考帧本身带位移，所有帧的相对位移都要减去它，否则整幅重建结果会稳定地偏一个方向。

4.2 梯度约束投影：软化阈值而不是硬砍

严格的 Stark-Oskoui 梯度投影要做逐点积分，实现麻烦、边界还容易出环。工程上我常用一个近似版本：只对超过阈值的那部分梯度做反向扩散。

function X = proj_gradient(X, T, beta) [gx, gy] = gradient(X); mag = hypot(gx, gy); excess = max(mag - T, 0); % 超出门限的量 ratio = excess ./ (mag + eps); % 归一化到 [0,1) lap = 4 * del2(X); % 离散拉普拉斯 X = X + beta * ratio .* lap / 4; % 只在过冲像素上做反向扩散 X = min(max(X, 0), 1); end

ratio是软权重：梯度刚过阈值时权重接近 0，几乎不动；远超阈值时权重趋近 1，按拉普拉斯方向拉回来。这个设计的好处是不会在平滑区域制造假的纹理——真正对应锋利边缘的像素才被处理。

阈值 T 的取法有讲究。图像归一化到 [0,1] 之后，自然图像的梯度幅度中位数一般在 0.02～0.08，T 取在 90 分位数附近比较合适，先跑一遍统计再定，别拍脑袋填。

g = hypot(gradient(X0), gradient(X0)); T0 = prctile(g(:), 90); % 用初始估计的梯度分布定阈值

beta控制投影强度，0.1～0.3 之间。它和 λ 会互相影响：λ 调大之后单轮残差注入变强，振铃也跟着变强，这时候要么把beta提一点，要么把 T 降一点。

4.3 幅度投影与参数联动表

幅度投影就是截断，但上下界不能瞎填：

function X = proj_amplitude(X, lo, hi) X = min(max(X, lo), hi); end

关键是lo/hi要和数据本身量纲一致。低分辨序列如果是uint8读进来的，先除以 255 归一化，否则 PSNR 计算和阈值全部错位，超分辨结果看起来「整体发黑」，你会以为是约束太强，实际是量纲没统一。

调参顺序我的习惯是：先固定 λ = 1.0 把psfSig和tvT调到 PSNR 不再上升，再回头把 λ 拉到 1.4 左右看是否还能涨。反过来先调 λ 会掩盖 PSF 建模的误差——你以为是松弛因子的功劳，其实是模糊核错了。

5. 排查清单：POCS 跑不出效果的五个典型翻车现场

下面五条基本覆盖了我遇到过的绝大多数问题，症状、成因、验证方法都写清楚，照着对照比盲调参数快得多。

5.1 现象：迭代到十几轮后 PSNR 反而开始下降

原因：数据保真投影用的是严格投影（δ = 0），噪声被当成必须满足的约束，算法开始拟合噪声。低分辨帧信噪比不高时，这一条几乎必然发生。

解决：把严格投影改成带容差的版本，δ取噪声标准差的 2～3 倍。实现上很简单，算完残差r后先做一次软阈值再反投影：

r = y - yh; r = sign(r) .* max(abs(r) - cfg.delta, 0); % 容差内残差不反投影

同时把hist曲线和 PSNR 曲线画在一起，两条线一个平台期就早停，别死磕maxIter。

5.2 现象：重建结果有重影，边缘分裂成两条

原因：配准误差。仿真时位移是自己给的不会有这问题，换真实序列就冒出来。用的是整像素配准、或者相位相关没做上采样细化，亚像素位移估偏了 0.3 像素以上，POCS 会把错误位移当成硬约束，反而把高频抹平。

解决：配准单独做，用相位相关拿到峰值后用 5 点抛物线插值细化到 0.01 像素；拿不准就先用已知位移的仿真数据对照一遍，确认算法本身没问题，再怀疑配准。

5.3 现象：图像四周出现一圈暗边或环状波纹

原因：傅里叶平移的循环边界，加上blur2卷积时的边界补零。两者都会在边缘制造不连续，POCS 每一轮都在放大这个不连续。

解决：imfilter明确指定'symmetric'，傅里叶平移前先做 8～16 像素的镜像扩展，平移完再裁回去。图省事也可以对整幅图乘一个 Tukey 窗，代价是边缘一小圈不参与评价。

5.4 现象：跑着跑着 MATLAB 提示内存不足

原因：LR存了全部全分辨率帧，每帧还要做 FFT，L = 30、r = 4的时候双精度内存直接翻十几倍。

解决：LR转single；FFT 分块处理；每轮只抽 3～5 帧。这三条一起上，内存占用能压到原来的四分之一以内，PSNR 损失不到 0.05 dB。

5.5 现象：结果整体偏灰，对比度明显塌陷

原因：三条并发的路都会导致这个症状——数据没归一化（uint8直接进算法）、幅度集上下界写成[0, 255]但数据已经归一化、或者upsample2里顺手乘了 r² 导致残差被系统性放大。

解决：入口处统一X = double(X) / 255，之后所有阈值、下界都按 [0,1] 走；upsample2坚持零插值不缩放；最后出图前imshow(X, [])只是显示拉伸，别拿它当真实动态范围去判断。

6. 验证与进阶：PSNR 曲线、早停判据和一个收敛加速技巧

评测本身不复杂，麻烦的是评价口径要统一：

function m = eval_metrics(X, GT) X = min(max(X, 0), 1); e = X - GT; m.psnr = 10 * log10(1 / mean(e(:).^2)); % 数据归一化到 [0,1] m.ssim = ssim(X, GT); fprintf('PSNR = %.2f dB, SSIM = %.4f\n', m.psnr, m.ssim); end

双三次上采样作为基线一起算，比较才有意义。r = 2、12 帧、噪声 0.005 的仿真条件下，双三次一般在 27 dB 上下，POCS 跑到收敛通常能到 31～33 dB，SSIM 从 0.78 提到 0.90 左右。如果只涨了 0.5 dB，先回去检查配准和 PSF，别急着改算法。

早停判据我用留一法：从 L 帧里留一帧不参与重建，每 5 轮用这一帧算一次残差，残差开始回升就停。这比盯着hist靠谱，因为hist只反映投影幅度，反映不了泛化。

加速方面有个便宜好用的招：给松弛更新加动量项，形式上是X = X + λ(Xp − X) + μ(X − Xprev)，μ 取 0.3～0.5。它等价于在投影方向上做一次轻度外推，实测能把达到同样 PSNR 的迭代次数砍掉三分之一。

% 主循环里替换单步更新 Xnew = X + cfg.lambda * (Xp - X) + cfg.mu * (X - Xprev); Xprev = X; X = min(max(Xnew, 0), 1); % 动量可能把值推出去，补一刀截断

μ 超过 0.7 之后震荡明显，尤其是梯度约束开得比较弱的时候。我现在的习惯是 μ 从 0.3 起步，每加 0.1 看一次hist曲线，出现锯齿状抖动就退回去。从那以后每次开新数据集，我都会先用已知位移的仿真跑一遍全流程、把 PSNR 曲线和hist曲线并排画出来贴在脚本注释里，再动真数据——这一步多花二十分钟，能省掉后面两天的自我怀疑。希望帮到你。

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