简介一份基于MATLAB与CVX库的图像上采样凸优化实现面向具备基础图像处理知识、希望用优化方法替代传统插值提升画质的研究者或工程师。资源围绕凸优化建模与PSNR评估两条主线展开清晰演示如何将插值权重求解转化为凸问题借助CVX自动求解从而抑制最近邻、双线性等方法易产生的伪影与边缘模糊。压缩包大小约22KB平台未给出文件总数与类型明细但从标题与描述看内容以MATLAB脚本.m为主代码结构紧凑适合直接阅读、调试和二次改造。目前已有797人学习浏览。通过该资源可掌握CVX定义、求解和验证优化问题的基本流程理解如何设计目标函数与约束条件来获得更优插值权重并学会用PSNR量化比较上采样结果相关思路还能迁移到图像去噪、超分辨率重建等任务中为后续研究提供可复用的建模参考与实验起点。1. 图像上采样别只想到插值用 CVX 做凸优化结果更稳我第一次把放大倍数调到 4 的时候双三次插值出来的边缘全是锯齿换了好几种插值核也压不下去。后来我把图像上采样当成一个“由低分辨率观测反推高分辨率原图”的反问题用 Matlab 里的 CVX 工具箱直接求解带正则项的凸优化目标反而拿到了稳定、可控的结果。这套流程不需要深度学习框架不依赖训练数据换一台机器结果完全一致特别适合做算法基线、对比实验和工业可复现场景。下面我从退化模型怎么建、目标函数怎么写、CVX 代码怎么落地三个角度展开最后再把实际跑图时踩过的坑和调参经验一并交代清楚。2. 建立模型之前为什么图像上采样是反问题又为什么凸优化能解2.1 上采样的数学本质是把丢失信息当成约束估计大多数成像过程都可以写成退化模型y D B x n其中 x 是待求的高分辨率图像B 表示光学模糊核D 表示下采样矩阵n 是噪声y 是我们手里拿到的低分辨率图像。上采样的任务不是把像素填出来而是根据 y 反推 x。这个反问题天然不适定下采样已经丢失了高频信息满足“下采样后等于 y”的高分辨率图像有无数张。插值只是从中选了一张看起来平滑的没有任何物理约束。凸优化的思路是加先验。我们把问题改写成带正则项的约束优化minimize 0.5 * ||D B x - y||² λ * R(x)前半部分是数据保真项强制重建结果在退化后与观测一致后半部分 R(x) 是正则项把“边缘清晰、局部平滑”这类先验写进去。这个形式最大的好处是当 R(x) 选 L1 范数或总变分范数时整个问题保持凸能保证收敛到全局最优不依赖初始化。从这个角度看上采样和去噪、去模糊是同一个框架下的三个实例。区别只是 D 是不是下采样矩阵、B 是不是已知的模糊核。后面第 6 章我会专门说这个迁移先记住这个模型整篇文章都是围绕它在展开。2.2 双三次插值在放大倍数变大时的边界三个具体表现双三次插值本质是局部核方法用邻域像素的加权组合来估计未知像素。它在 1.5 到 2 倍放大时表现尚可一旦超过 3 倍问题就很明显第一是边缘出现锯齿和振铃。因为高频信息已经物理性丢失局部核只能猜而核的形状固定遇到陡峭边缘就会过冲。第二是纹理区域糊成一片。头发丝、草地、栅格这类周期性纹理双三次核会把相邻像素平均掉细节全部抹平。第三是误差是结构性的不是随机噪声。这个很关键插值误差集中在边缘附近且方向一致后续做配准、测量、检测时会产生系统性偏差不是加一点平滑就能补救的。我并不是说双三次插值没有用。在做初值、显示预览、放大倍数小于 2 的场景下它依然是最快的选择。但当你需要拿上采样结果去做定量分析或者放大倍数逼近 3 以上就必须换一个能约束全局一致性的方法。凸优化正好补上这个空档。2.3 凸优化和 CVX 的适用边界不是所有能量函数都能直接丢进去CVX 是一个基于 Matlab 的建模框架它要求你的目标函数符合 DCP 规则也就是只能由凸函数、凹函数和仿射函数按固定规则组合。很多东西随手写出来是不符合规则的比如我自己就犯过把二次型展开成 (Ax-y)(A*x-y) 直接写进 CVX 的错结果直接报 DCP 错误。图像反问题里常见的几类算子都能在 CVX 里安全表达数据保真项L2 用 sum_square(Ax-b)L1 用 norm(Ax-b,1)都是标准写法正则项L1 范数、各向同性或各向异性 TV都可以拆成若干范数的组合约束条件像素范围 0x255、能量守恒、线性等式约束全部支持因此这套方法很适合做原型验证和中小尺寸图像的正式实验。代价是 CVX 内部会做大幅度的问题转换变量一多内存和求解时间会急剧上升。这决定了它的应用边界整图做适合几百像素见方的小图大图必须拆 patch。第 4 章我会给出具体做法。2.4 什么场景值得选 CVX 上采样我的判断标准有三条第一你手里没有现成的深度学习环境或者你想先跑一个不吃数据的强基线第二你的图像尺寸可控或者可以接受分块重建第三你需要在论文或项目文档里写清楚“我的模型是凸的解是全局最优的”而不是“网络训练后效果变好了”。凡是满足其中两条的都值得把 CVX 方案作为首选或对照方法。反过来如果你的低分辨率图是 1000 万像素、要求 4 倍放大且必须在秒级出结果那凸优化不是好选择patch 拼接带来的缝合问题会让你耗尽耐心。这种情况我一般直接劝退去做深度学习或传统插值加锐化的组合别硬用 CVX。3. 搭建一个能复现的上采样目标退化矩阵、正则项和 CVX 最小实现3.1 准备一对有真值可对比的测试图像做算法实验不能上来就处理真实低分辨率图因为没有高分辨率真值你根本分不清结果是变好了还是变难看了。我习惯先用一张高分辨率灰度图人为加上模糊、下采样和噪声生成一张低分辨率观测图。这样后面可以用真实高分辨率图算 PSNR。clear; clc; s 2; % 放大倍数先不要超过 2 I_hr im2double(imread(cameraman.tif)); % 模拟光学模糊高斯核尺寸和方差按相机 PSF 估计 h fspecial(gaussian, 5, 1.2); I_blur imfilter(I_hr, h, conv, replicate); % 块平均降采样等价于后面 D 矩阵定义的采样过程 I_lr imresize(I_blur, 1/s, box); % 加一点高斯噪声模拟传感器噪声 I_lr I_lr 0.003 * randn(size(I_lr)); imwrite(I_lr, obs.png); [m, n] size(I_lr); M m * s; N n * s; % 高分辨率尺寸这里有一个关键选择低分辨率观测的生成方式必须和你后面构造的下采样矩阵 D 保持一致。我用了 box 降采样它在每个 s×s 块内取平均。如果你在这里用双三次降采样而后面 D 矩阵还是块平均那模型和观测就对不上结果一定出问题。这一点在真实图像上尤其重要后面避坑章节还会展开。3.2 构造下采样矩阵 D块平均模型与向量化写法D 的尺寸是 mn 行MN 列。每一行对应低分辨率上的一个像素行内非零位置是该像素对应的高分辨率 s×s 块数值取 1/s²也就是做块内平均。这个矩阵在数学上定义得很清楚但直接写双循环在中等尺寸下也够慢所以我把索引部分改成了向量化写法。D sparse(m*n, M*N); row 0; for i 1:m for j 1:n row row 1; base ((j-1)*s) * M (i-1)*s; % 按列主序生成块内线性偏移顺序为 (行偏移, 列偏移) offset kron((0:s-1), ones(s,1)) * M ... kron(ones(s,1), (0:s-1)); D(row, base offset 1) 1 / (s*s); end end B speye(M*N); % 模糊矩阵先设为单位阵跑通后再换成真实 PSF代码里最需要注意的是 MATLAB 的线性索引按列主序排列。高分辨率图像第 i 行、第 j 列的像素在线性向量中的位置是 (j-1)*M i。所以一个块的基准点不是 (i-1)*s (j-1)*s 这种行优先写法必须先乘以 M。我在这边上过一次当构造出的 D 矩阵转置后怎么都对不上最后把第 20 行的非零位置打印出来才看出索引错位。B 矩阵当前是单位阵不代表后面不用换。它存在的意义是让模型保留“先模糊后降采样”的表达能力。等你在模拟数据上把流程验证通了再把 B 替换成真实成像系统的 PSF 对应的循环矩阵。实际情况中块平均的下采样模型本身已经隐含了一点模糊效果所以很多人直接让 B 保持单位阵也能得到可接受结果。3.3 构造差分矩阵并选择正则项为什么先用各向异性 TVTV 正则的解释是图像在局部应该是分片平滑的但边缘处允许突变。实现上需要两个稀疏差分矩阵分别计算水平方向和垂直方向的相邻像素差。Tv sparse(M*N, M*N); % 垂直方向差分 Th sparse(M*N, M*N); % 水平方向差分 for ii 1:M for jj 1:N idx (jj-1)*M ii; if ii M Tv(idx, idx) 1; Tv(idx, idx1) -1; end if jj N Th(idx, idx) 1; Th(idx, idxM) -1; end end end这里用的是各向异性 TV也就是垂直和水平方向分开算绝对值再相加|Tvx| |Thx|。另一个常见选择是各向同性 TV形式是 sqrt((Tvx)²(Thx)²)边缘更平滑但表达式稍复杂。我推荐先用各向异性版本理由很实际它在 CVX 里的表达简单不增加求解器负担而且第一次写代码时不容易出错。等你把保真项、约束和 λ 都调顺了再换各向同性版本对比效果不迟。还有一个很容易混淆的点为什么要用 L1 而不是 L2 来做正则项L2 正则也叫 Tikhonov它惩罚的是像素差平方大梯度会被压得很狠结果就是边缘变软。L1 范数对小幅值惩罚相对温和对大幅值的惩罚与幅值成正比能保留陡峭边缘。图像上采样最需要保的就是边缘所以默认 TV 而不是 Tikhonov。3.4 核心CVX 最小上采样代码下面这段就是整个流程的主干。把观测图像拉成向量声明高分辨率变量写出保真项和正则项的凸组合。y I_lr(:); lambda 0.02; % 正则化权重经验区间 0.005~0.05 cvx_begin quiet variable x(M*N, 1) minimize( 0.5 * sum_square(D * B * x - y) ... lambda * (norm(Tv * x, 1) norm(Th * x, 1)) ) subject to x 0; x 1; cvx_end I_rec reshape(x, M, N);这个目标函数里的三项各有分工。0.5 * sum_square(DBx - y) 是 L2 保真项假设噪声是高斯分布它的梯度是平滑的不会引入椒盐噪声。norm(Tvx,1) 和 norm(Thx,1) 是 TV 正则控制重建图的分片平滑程度。λ 是两者的权重平衡“贴合观测”和“符合先验”。边界约束 x 0 和 x 1 是很多人会忘记的一步。图像是 0 到 1 的归一化量如果不加约束优化器可能为了极小化保真项把某些像素推到负值或超过 1结果图虽然 PSNR 看起来还行但显示出来全是过曝或死黑的点。加了约束后问题依然是凸的CVX 处理起来没有任何障碍。我特别强调一下如果你的输入图是 0 到 255 的 uint8请在进入 CVX 之前统一转成 double 并归一化到 0 到 1。否则 λ 的取法、约束边界、TV 项的尺度全都要跟着变调参就是一场灾难。这也是我最早把代码写得像玄学的原因之一。3.5 目标函数的三个细节能量守恒、初始化与求解器第一是能量守恒。构造 D 时行和要严格等于 1代表块内求平均后能量不变。验证方法是在 Matlab 里执行 full(sum(D,2))如果某个块因为索引越界少写了像素行和会小于 1重建结果就会出现整体变暗或局部发灰。第二是初始化。CVX 本质上不使用我们的初始化做热启动它的求解器从零向量开始。所以你在外面给 x 赋任何初值都不会影响求解结果唯一能做的是通过 precompute 减少建模时间。第三是求解器选择。CVX 默认的 SDPT3 在中小规模问题上比较稳如果碰到 L1 项导致模型转换后求解过慢可以换成 SeDuMi 试试同一个目标函数有时能快一倍但精度略低。4. 工程落地参数设置、patch 拼接和客观评估4.1 关键参数一览放大倍数、模糊核、λ 的经验区间参数作用经验取值s放大倍数2 最稳妥3 以上需配合更真实的 PSF模糊核 B光学退化建模先设单位阵跑通有 PSF 测量值再替换λ正则化权重0.005~0.05图像归一化到 0~1 时有效图像尺寸决定变量个数整图建议 M*N 不超过 30000大图拆 patch求解精度cvx_precision原型验证用默认批量实验用 low 加速λ 这个参数没有万能值。我的习惯是先固定一张小图把 λ 按 0.005、0.01、0.02、0.05、0.1 扫一遍画出 PSNR 曲线。正常情况下曲线是单峰的峰值对应的 λ 就是当前场景下的最优值。如果曲线一直单调上升或下降说明保真项和正则项的量纲不一致先检查 D 是否能量守恒再检查图像是否归一化而不是去继续调 λ。4.2 大图跑不动时的 patch 拼接方案CVX 的变量规模上去后问题转换阶段就会耗尽内存。我常用的上限是 M*N 不超过 30000 个变量也就是大约 173×173 的高分辨率图。实际工程中哪有这么小的图所以必须拆块。把低分辨率图切成若干小块每个小块独立做上采样再拼回全图。拼接时用重叠和加权平均消除接缝。function I_rec upsample_cvx_patch(I_lr_patch, s, lambda) % 输入一个低分辨率小块返回对应高分辨率块 [pm, pn] size(I_lr_patch); pM pm * s; pN pn * s; % 这里省略 D、Tv、Th 的构造直接复用第 3 章的代码 y I_lr_patch(:); cvx_begin quiet variable x(pM*pN, 1) minimize( 0.5*sum_square(D*x - y) ... lambda*(norm(Tv*x,1) norm(Th*x,1)) ) subject to x 0; x 1; cvx_end I_rec reshape(x, pM, pN); end调用时按重叠步长滑动窗口重建结果先累加最后除以权重图。patch_size 48; % 低分辨率 patch 尺寸 overlap 8; % 重叠像素 step patch_size - overlap; acc zeros(M, N); cnt zeros(M, N); lambda 0.02; for i 1:step:(m-patch_size1) for j 1:step:(n-patch_size1) sub I_lr(i:ipatch_size-1, j:jpatch_size-1); rec upsample_cvx_patch(sub, s, lambda); rows (i-1)*s (1:size(rec,1)); cols (j-1)*s (1:size(rec,2)); acc(rows, cols) acc(rows, cols) rec; cnt(rows, cols) cnt(rows, cols) 1; end end % 处理不够一个 patch 的右缘和下缘 if i patch_size - 1 m for j 1:step:(n-patch_size1) sub I_lr(m-patch_size1:m, j:jpatch_size-1); rec upsample_cvx_patch(sub, s, lambda); rows (m-patch_size)*s (1:size(rec,1)); cols (j-1)*s (1:size(rec,2)); acc(rows, cols) acc(rows, cols) rec; cnt(rows, cols) cnt(rows, cols) 1; end end I_rec_patch acc ./ max(cnt, 1);重叠区域被多个 patch 重建了多次取平均可以抑制接缝。如果接缝仍然可见把重叠加大到 16 像素或者对每个 patch 的贡献加一个以中心为高值的窗函数再加权融合。代价是计算量会上升所以我在大多数场景下用重叠平均就够了。4.3 用 PSNR 和边缘响应做客观评估PSNR 是图像重建领域最常用的指标但它对边缘的刻画很弱两个 PSNR 相同的重建图可能在视觉上差别巨大。我建议同时计算 PSNR 和边缘响应差值。function p psnr_db(x, ref) % 输入图像范围 0~1 x double(x); ref double(ref); mse_val mean((x(:) - ref(:)).^2); p 10 * log10(1 / (mse_val eps)); end I_bic imresize(I_lr, s, bicubic); fprintf(bicubic PSNR: %.2f dB\n, psnr_db(I_bic, I_hr)); fprintf(cvx PSNR: %.2f dB\n, psnr_db(I_rec, I_hr)); gx fspecial(sobel); gy fspecial(sobel); edge_cvx abs(imfilter(I_rec, gx, replicate)) abs(imfilter(I_rec, gy, replicate)); edge_bic abs(imfilter(I_bic, gx, replicate)) abs(imfilter(I_bic, gy, replicate)); % 边缘强度均值差值数值越小代表边缘形态越接近真值 delta_edge mean(abs(edge_cvx(:) - edge_bic(:)));拿到这几个数值后判断标准是PSNR 比双三次高至少 0.5 dB同时边缘差值显著下降。如果你跑出来 PSNR 高了但边缘差值反而变大多半是 λ 太小放大噪声。如果 PSNR 和双三次持平但视觉上更干净也说明方法本身有效只是指标没体现出来。这种情况我会再截一块纹理区放大对比主观判断在工程报告里同样要有据可依。5. 五个避坑记录从翻车里总结的 CVX 上采样经验5.1 报错“Disciplined convex programming error”不见得是你目标函数写错现象CVX 在解析阶段直接报错提示大意是“Cannot perform the operation {convex} ./ {convex}”。第一反应是目标函数选错了但问题是这段代码明明在别人那里能跑。原因CVX 只认它自己的原子函数组合。你把二次型写成了(D*x-y)*(D*x-y)或者把 L2 范数平方写成了norm(D*x-y,2)^2这两种写法在普通 Matlab 里合法在 CVX 里不被识别为凸表达。DCP 规则不允许直接对凸函数做任意幂次运算。解决二次型一律用sum_square(D*x-y)L1 项一律用norm(...,1)不要自己展开矩阵乘积。我还在目标里写过(D*x-y).^2这种逐元素平方同样报错换成sum_square立刻通过。记住一条原则CVX 的表达式只能用原子函数搭积木不能按数值计算的习惯展开。5.2 小图能跑大图直接内存爆炸现象64×64 的低分辨率图跑得挺顺利换到 256×256 的图CVX 在建模阶段就卡住后来直接提示内存不足。原因变量数等于高分辨率像素数。256×256 的低分辨率图放大 2 倍后有 512×512262144 个变量。CVX 处理 L1 项时会把问题转成一个大规模 SOCP内部变量还会成倍扩张内存消耗远不是稀疏矩阵表面看起来那么小。解决不要在整图上直接跑把输入切成 48×48 或 64×64 的低分辨率 patch逐块求解再拼接。如果业务不允许拼接就接受现实换成其他方法。CVX 适合做中小尺寸和原型验证不适合做大图整张重建这是方法论边界不是代码优化能解决的。5.3 重建图整体发灰亮度对不上现象求解完成后图像确实比双三次插值更清楚但整体发灰白色区域被压到了 0.8 左右黑色区域抬到 0.2 左右。原因最常见的是 D 矩阵能量不守恒。比如你用双三次下采样生成低分辨率观测但 D 矩阵却按块平均构造这会导致 D*x 的值域和 y 的值域有系统偏差。另一个原因是输入是 uint8 的 0 到 255直接放进了约束 0x255但 λ 还是按 0~1 范围设的正则项相对过强把像素往中间值拉。解决先执行full(sum(D,2))检查每一行确认非零和为 1。再统一数据范围输入图在 im2double 转成 0~1 后约束写 0x1λ 按 0.005~0.05 扫描。如果问题依旧把生成低分辨率观测的方式改成和 D 相同的 block 平均模型先做到自洽再谈真实退化。5.4 输出图像像马赛克拼接痕迹和网格状纹理明显现象重建结果边缘锐利但能隐约看到低分辨率像素块的网格像用积木拼出来的画。放大倍数越大网格越明显。原因这是 D 矩阵块平均模型带来的伪影。块平均相当于用一个矩形的低通滤波器它和现实中的光学 PSF 差别很大。当块边界和图像边缘恰好重合时重建结果会强化这种格子感。另一个诱因是 λ 过大平滑项把块内部的差异压平块之间又靠边缘检测拉回形成方块轮廓。解决把模糊矩阵 B 从单位阵换成高斯核对应的循环矩阵让退化模型变成“先模糊、后块平均”这样网格伪影会明显减弱。具体做法是用fspecial(gaussian,[5 5],1.0)生成核再用稀疏 Toeplitz 结构构造 B。如果暂时不想引入 B可以把 λ 调小到 0.01 量级观察网格是否缓解但这治标不治本。5.5 跑完比双三次还差PSNR 掉 0.5 dB 以上现象CVX 老老实实收敛了可拿到 PSNR 一看还不如imresize(...,bicubic)。求解时间花了十几分钟结果出来是负优化。原因大概率是正则项权重 λ 和噪声水平不匹配。λ 太小时重建会把下采样引入的噪声放大λ 太大时边缘和纹理被过度平滑丢失双三次还能保留的局部反差。另一个原因是观测图像本身已经被处理过比如你拿来的是一个已经增强或锐化过的低分辨率图它和 y DBxn 的假设不符凸优化没有任何高频信息可以恢复。解决先把 λ 在 0.005 到 0.05 之间做一次小扫描固定一张小图画出 PSNR 曲线。如果曲线最高点仍低于双三次那就检查退化模型。用合成数据进行闭环验证从一张清晰图生成 y再上采样确认流程能超过双三次然后再应用到真实图上。真实图若效果始终不理想说明该退化场景下凸优化的假设不成立放弃这条路不要硬调参。这个判断标准也回答了很多人的疑问凸优化上采样不是万能的但它值得做因为当你匹配好退化模型后效果是稳定且可解释的。6. 进阶用边缘自适应权重把 TV 正则局部化再把同一套框架迁移到去噪固定 λ 的 TV 正则有一个尴尬平滑区需要大 λ 来抑制噪声边缘附近需要小 λ 来保留对比度一个全局值无法同时满足。解决办法是给 TV 项加一个预计算的权重矩阵让每一点的正则强度不同。权重来自双三次插值结果的边缘检测图先算出一幅 0 到 1 的边缘响应图再映射成惩罚权重。x_init imresize(I_lr, s, bicubic); edge_map abs(imfilter(x_init, fspecial(sobel), replicate)); edge_map edge_map(:) / max(edge_map(:) eps); % 边缘响应大的地方权重小保留边缘平坦区权重大加强平滑 w 0.3 0.7 * (edge_map 0.35); W spdiags(w, 0, M*N, M*N); cvx_begin quiet variable x(M*N, 1) minimize( 0.5 * sum_square(D * x - y) ... lambda * norm(W * Tv * x, 1) ... lambda * norm(W * Th * x, 1) ) subject to x 0; x 1; cvx_end注意权重矩阵 W 必须是在进入 CVX 之前就固定好的常数矩阵不能用优化变量参与计算。如果你把 w 写成 x 的函数模型立刻违反 DCP 规则。我在第一次尝试时就是想让权重跟着迭代变化结果 CVX 报错这才改成用插值结果做预估计效果也不差。这个技巧让同一张图上的平坦区域和纹理区域不用再手工分成两个模型综合效果比固定 λ 好一个档。更重要的是这套框架可以整体迁移。把 D 换成单位矩阵低分辨率观测换成噪声图目标函数就变成了一个带 TV 正则的去噪模型。把 B 换成运动模糊的卷积矩阵它就变成了去模糊模型。我在做某跨平台系统的图像预处理模块时把同一个 CVX 求解函数复制了三份只改了 D 和 B 的构造部分就同时交付了上采样、去噪、去模糊三个功能。数学核心完全一致工程量只差在退化矩阵的构造上。最后说一个血泪教训我前几次跑出块状伪影时一直以为是自己 λ 没调好反复扫参扫了一个下午。后来才发现问题根本不在 λ而在 D 矩阵的构造和观测图像的降采样方式不一致。从那以后我再做任何上采样实验都坚持先造一对“已知真值”的合成数据做闭环验证再碰真实图像。希望这篇文章能帮你跳过这一段弯路直接看到凸优化上采样真正有效的样子。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/13 17:21:45 · 阅读 1.2k 图像上采样别只靠插值：用CVX凸优化重建高清图像
简介一份基于MATLAB与CVX库的图像上采样凸优化实现面向具备基础图像处理知识、希望用优化方法替代传统插值提升画质的研究者或工程师。资源围绕凸优化建模与PSNR评估两条主线展开清晰演示如何将插值权重求解转化为凸问题借助CVX自动求解从而抑制最近邻、双线性等方法易产生的伪影与边缘模糊。压缩包大小约22KB平台未给出文件总数与类型明细但从标题与描述看内容以MATLAB脚本.m为主代码结构紧凑适合直接阅读、调试和二次改造。目前已有797人学习浏览。通过该资源可掌握CVX定义、求解和验证优化问题的基本流程理解如何设计目标函数与约束条件来获得更优插值权重并学会用PSNR量化比较上采样结果相关思路还能迁移到图像去噪、超分辨率重建等任务中为后续研究提供可复用的建模参考与实验起点。1. 图像上采样别只想到插值用 CVX 做凸优化结果更稳我第一次把放大倍数调到 4 的时候双三次插值出来的边缘全是锯齿换了好几种插值核也压不下去。后来我把图像上采样当成一个“由低分辨率观测反推高分辨率原图”的反问题用 Matlab 里的 CVX 工具箱直接求解带正则项的凸优化目标反而拿到了稳定、可控的结果。这套流程不需要深度学习框架不依赖训练数据换一台机器结果完全一致特别适合做算法基线、对比实验和工业可复现场景。下面我从退化模型怎么建、目标函数怎么写、CVX 代码怎么落地三个角度展开最后再把实际跑图时踩过的坑和调参经验一并交代清楚。2. 建立模型之前为什么图像上采样是反问题又为什么凸优化能解2.1 上采样的数学本质是把丢失信息当成约束估计大多数成像过程都可以写成退化模型y D B x n其中 x 是待求的高分辨率图像B 表示光学模糊核D 表示下采样矩阵n 是噪声y 是我们手里拿到的低分辨率图像。上采样的任务不是把像素填出来而是根据 y 反推 x。这个反问题天然不适定下采样已经丢失了高频信息满足“下采样后等于 y”的高分辨率图像有无数张。插值只是从中选了一张看起来平滑的没有任何物理约束。凸优化的思路是加先验。我们把问题改写成带正则项的约束优化minimize 0.5 * ||D B x - y||² λ * R(x)前半部分是数据保真项强制重建结果在退化后与观测一致后半部分 R(x) 是正则项把“边缘清晰、局部平滑”这类先验写进去。这个形式最大的好处是当 R(x) 选 L1 范数或总变分范数时整个问题保持凸能保证收敛到全局最优不依赖初始化。从这个角度看上采样和去噪、去模糊是同一个框架下的三个实例。区别只是 D 是不是下采样矩阵、B 是不是已知的模糊核。后面第 6 章我会专门说这个迁移先记住这个模型整篇文章都是围绕它在展开。2.2 双三次插值在放大倍数变大时的边界三个具体表现双三次插值本质是局部核方法用邻域像素的加权组合来估计未知像素。它在 1.5 到 2 倍放大时表现尚可一旦超过 3 倍问题就很明显第一是边缘出现锯齿和振铃。因为高频信息已经物理性丢失局部核只能猜而核的形状固定遇到陡峭边缘就会过冲。第二是纹理区域糊成一片。头发丝、草地、栅格这类周期性纹理双三次核会把相邻像素平均掉细节全部抹平。第三是误差是结构性的不是随机噪声。这个很关键插值误差集中在边缘附近且方向一致后续做配准、测量、检测时会产生系统性偏差不是加一点平滑就能补救的。我并不是说双三次插值没有用。在做初值、显示预览、放大倍数小于 2 的场景下它依然是最快的选择。但当你需要拿上采样结果去做定量分析或者放大倍数逼近 3 以上就必须换一个能约束全局一致性的方法。凸优化正好补上这个空档。2.3 凸优化和 CVX 的适用边界不是所有能量函数都能直接丢进去CVX 是一个基于 Matlab 的建模框架它要求你的目标函数符合 DCP 规则也就是只能由凸函数、凹函数和仿射函数按固定规则组合。很多东西随手写出来是不符合规则的比如我自己就犯过把二次型展开成 (Ax-y)(A*x-y) 直接写进 CVX 的错结果直接报 DCP 错误。图像反问题里常见的几类算子都能在 CVX 里安全表达数据保真项L2 用 sum_square(Ax-b)L1 用 norm(Ax-b,1)都是标准写法正则项L1 范数、各向同性或各向异性 TV都可以拆成若干范数的组合约束条件像素范围 0x255、能量守恒、线性等式约束全部支持因此这套方法很适合做原型验证和中小尺寸图像的正式实验。代价是 CVX 内部会做大幅度的问题转换变量一多内存和求解时间会急剧上升。这决定了它的应用边界整图做适合几百像素见方的小图大图必须拆 patch。第 4 章我会给出具体做法。2.4 什么场景值得选 CVX 上采样我的判断标准有三条第一你手里没有现成的深度学习环境或者你想先跑一个不吃数据的强基线第二你的图像尺寸可控或者可以接受分块重建第三你需要在论文或项目文档里写清楚“我的模型是凸的解是全局最优的”而不是“网络训练后效果变好了”。凡是满足其中两条的都值得把 CVX 方案作为首选或对照方法。反过来如果你的低分辨率图是 1000 万像素、要求 4 倍放大且必须在秒级出结果那凸优化不是好选择patch 拼接带来的缝合问题会让你耗尽耐心。这种情况我一般直接劝退去做深度学习或传统插值加锐化的组合别硬用 CVX。3. 搭建一个能复现的上采样目标退化矩阵、正则项和 CVX 最小实现3.1 准备一对有真值可对比的测试图像做算法实验不能上来就处理真实低分辨率图因为没有高分辨率真值你根本分不清结果是变好了还是变难看了。我习惯先用一张高分辨率灰度图人为加上模糊、下采样和噪声生成一张低分辨率观测图。这样后面可以用真实高分辨率图算 PSNR。clear; clc; s 2; % 放大倍数先不要超过 2 I_hr im2double(imread(cameraman.tif)); % 模拟光学模糊高斯核尺寸和方差按相机 PSF 估计 h fspecial(gaussian, 5, 1.2); I_blur imfilter(I_hr, h, conv, replicate); % 块平均降采样等价于后面 D 矩阵定义的采样过程 I_lr imresize(I_blur, 1/s, box); % 加一点高斯噪声模拟传感器噪声 I_lr I_lr 0.003 * randn(size(I_lr)); imwrite(I_lr, obs.png); [m, n] size(I_lr); M m * s; N n * s; % 高分辨率尺寸这里有一个关键选择低分辨率观测的生成方式必须和你后面构造的下采样矩阵 D 保持一致。我用了 box 降采样它在每个 s×s 块内取平均。如果你在这里用双三次降采样而后面 D 矩阵还是块平均那模型和观测就对不上结果一定出问题。这一点在真实图像上尤其重要后面避坑章节还会展开。3.2 构造下采样矩阵 D块平均模型与向量化写法D 的尺寸是 mn 行MN 列。每一行对应低分辨率上的一个像素行内非零位置是该像素对应的高分辨率 s×s 块数值取 1/s²也就是做块内平均。这个矩阵在数学上定义得很清楚但直接写双循环在中等尺寸下也够慢所以我把索引部分改成了向量化写法。D sparse(m*n, M*N); row 0; for i 1:m for j 1:n row row 1; base ((j-1)*s) * M (i-1)*s; % 按列主序生成块内线性偏移顺序为 (行偏移, 列偏移) offset kron((0:s-1), ones(s,1)) * M ... kron(ones(s,1), (0:s-1)); D(row, base offset 1) 1 / (s*s); end end B speye(M*N); % 模糊矩阵先设为单位阵跑通后再换成真实 PSF代码里最需要注意的是 MATLAB 的线性索引按列主序排列。高分辨率图像第 i 行、第 j 列的像素在线性向量中的位置是 (j-1)*M i。所以一个块的基准点不是 (i-1)*s (j-1)*s 这种行优先写法必须先乘以 M。我在这边上过一次当构造出的 D 矩阵转置后怎么都对不上最后把第 20 行的非零位置打印出来才看出索引错位。B 矩阵当前是单位阵不代表后面不用换。它存在的意义是让模型保留“先模糊后降采样”的表达能力。等你在模拟数据上把流程验证通了再把 B 替换成真实成像系统的 PSF 对应的循环矩阵。实际情况中块平均的下采样模型本身已经隐含了一点模糊效果所以很多人直接让 B 保持单位阵也能得到可接受结果。3.3 构造差分矩阵并选择正则项为什么先用各向异性 TVTV 正则的解释是图像在局部应该是分片平滑的但边缘处允许突变。实现上需要两个稀疏差分矩阵分别计算水平方向和垂直方向的相邻像素差。Tv sparse(M*N, M*N); % 垂直方向差分 Th sparse(M*N, M*N); % 水平方向差分 for ii 1:M for jj 1:N idx (jj-1)*M ii; if ii M Tv(idx, idx) 1; Tv(idx, idx1) -1; end if jj N Th(idx, idx) 1; Th(idx, idxM) -1; end end end这里用的是各向异性 TV也就是垂直和水平方向分开算绝对值再相加|Tvx| |Thx|。另一个常见选择是各向同性 TV形式是 sqrt((Tvx)²(Thx)²)边缘更平滑但表达式稍复杂。我推荐先用各向异性版本理由很实际它在 CVX 里的表达简单不增加求解器负担而且第一次写代码时不容易出错。等你把保真项、约束和 λ 都调顺了再换各向同性版本对比效果不迟。还有一个很容易混淆的点为什么要用 L1 而不是 L2 来做正则项L2 正则也叫 Tikhonov它惩罚的是像素差平方大梯度会被压得很狠结果就是边缘变软。L1 范数对小幅值惩罚相对温和对大幅值的惩罚与幅值成正比能保留陡峭边缘。图像上采样最需要保的就是边缘所以默认 TV 而不是 Tikhonov。3.4 核心CVX 最小上采样代码下面这段就是整个流程的主干。把观测图像拉成向量声明高分辨率变量写出保真项和正则项的凸组合。y I_lr(:); lambda 0.02; % 正则化权重经验区间 0.005~0.05 cvx_begin quiet variable x(M*N, 1) minimize( 0.5 * sum_square(D * B * x - y) ... lambda * (norm(Tv * x, 1) norm(Th * x, 1)) ) subject to x 0; x 1; cvx_end I_rec reshape(x, M, N);这个目标函数里的三项各有分工。0.5 * sum_square(DBx - y) 是 L2 保真项假设噪声是高斯分布它的梯度是平滑的不会引入椒盐噪声。norm(Tvx,1) 和 norm(Thx,1) 是 TV 正则控制重建图的分片平滑程度。λ 是两者的权重平衡“贴合观测”和“符合先验”。边界约束 x 0 和 x 1 是很多人会忘记的一步。图像是 0 到 1 的归一化量如果不加约束优化器可能为了极小化保真项把某些像素推到负值或超过 1结果图虽然 PSNR 看起来还行但显示出来全是过曝或死黑的点。加了约束后问题依然是凸的CVX 处理起来没有任何障碍。我特别强调一下如果你的输入图是 0 到 255 的 uint8请在进入 CVX 之前统一转成 double 并归一化到 0 到 1。否则 λ 的取法、约束边界、TV 项的尺度全都要跟着变调参就是一场灾难。这也是我最早把代码写得像玄学的原因之一。3.5 目标函数的三个细节能量守恒、初始化与求解器第一是能量守恒。构造 D 时行和要严格等于 1代表块内求平均后能量不变。验证方法是在 Matlab 里执行 full(sum(D,2))如果某个块因为索引越界少写了像素行和会小于 1重建结果就会出现整体变暗或局部发灰。第二是初始化。CVX 本质上不使用我们的初始化做热启动它的求解器从零向量开始。所以你在外面给 x 赋任何初值都不会影响求解结果唯一能做的是通过 precompute 减少建模时间。第三是求解器选择。CVX 默认的 SDPT3 在中小规模问题上比较稳如果碰到 L1 项导致模型转换后求解过慢可以换成 SeDuMi 试试同一个目标函数有时能快一倍但精度略低。4. 工程落地参数设置、patch 拼接和客观评估4.1 关键参数一览放大倍数、模糊核、λ 的经验区间参数作用经验取值s放大倍数2 最稳妥3 以上需配合更真实的 PSF模糊核 B光学退化建模先设单位阵跑通有 PSF 测量值再替换λ正则化权重0.005~0.05图像归一化到 0~1 时有效图像尺寸决定变量个数整图建议 M*N 不超过 30000大图拆 patch求解精度cvx_precision原型验证用默认批量实验用 low 加速λ 这个参数没有万能值。我的习惯是先固定一张小图把 λ 按 0.005、0.01、0.02、0.05、0.1 扫一遍画出 PSNR 曲线。正常情况下曲线是单峰的峰值对应的 λ 就是当前场景下的最优值。如果曲线一直单调上升或下降说明保真项和正则项的量纲不一致先检查 D 是否能量守恒再检查图像是否归一化而不是去继续调 λ。4.2 大图跑不动时的 patch 拼接方案CVX 的变量规模上去后问题转换阶段就会耗尽内存。我常用的上限是 M*N 不超过 30000 个变量也就是大约 173×173 的高分辨率图。实际工程中哪有这么小的图所以必须拆块。把低分辨率图切成若干小块每个小块独立做上采样再拼回全图。拼接时用重叠和加权平均消除接缝。function I_rec upsample_cvx_patch(I_lr_patch, s, lambda) % 输入一个低分辨率小块返回对应高分辨率块 [pm, pn] size(I_lr_patch); pM pm * s; pN pn * s; % 这里省略 D、Tv、Th 的构造直接复用第 3 章的代码 y I_lr_patch(:); cvx_begin quiet variable x(pM*pN, 1) minimize( 0.5*sum_square(D*x - y) ... lambda*(norm(Tv*x,1) norm(Th*x,1)) ) subject to x 0; x 1; cvx_end I_rec reshape(x, pM, pN); end调用时按重叠步长滑动窗口重建结果先累加最后除以权重图。patch_size 48; % 低分辨率 patch 尺寸 overlap 8; % 重叠像素 step patch_size - overlap; acc zeros(M, N); cnt zeros(M, N); lambda 0.02; for i 1:step:(m-patch_size1) for j 1:step:(n-patch_size1) sub I_lr(i:ipatch_size-1, j:jpatch_size-1); rec upsample_cvx_patch(sub, s, lambda); rows (i-1)*s (1:size(rec,1)); cols (j-1)*s (1:size(rec,2)); acc(rows, cols) acc(rows, cols) rec; cnt(rows, cols) cnt(rows, cols) 1; end end % 处理不够一个 patch 的右缘和下缘 if i patch_size - 1 m for j 1:step:(n-patch_size1) sub I_lr(m-patch_size1:m, j:jpatch_size-1); rec upsample_cvx_patch(sub, s, lambda); rows (m-patch_size)*s (1:size(rec,1)); cols (j-1)*s (1:size(rec,2)); acc(rows, cols) acc(rows, cols) rec; cnt(rows, cols) cnt(rows, cols) 1; end end I_rec_patch acc ./ max(cnt, 1);重叠区域被多个 patch 重建了多次取平均可以抑制接缝。如果接缝仍然可见把重叠加大到 16 像素或者对每个 patch 的贡献加一个以中心为高值的窗函数再加权融合。代价是计算量会上升所以我在大多数场景下用重叠平均就够了。4.3 用 PSNR 和边缘响应做客观评估PSNR 是图像重建领域最常用的指标但它对边缘的刻画很弱两个 PSNR 相同的重建图可能在视觉上差别巨大。我建议同时计算 PSNR 和边缘响应差值。function p psnr_db(x, ref) % 输入图像范围 0~1 x double(x); ref double(ref); mse_val mean((x(:) - ref(:)).^2); p 10 * log10(1 / (mse_val eps)); end I_bic imresize(I_lr, s, bicubic); fprintf(bicubic PSNR: %.2f dB\n, psnr_db(I_bic, I_hr)); fprintf(cvx PSNR: %.2f dB\n, psnr_db(I_rec, I_hr)); gx fspecial(sobel); gy fspecial(sobel); edge_cvx abs(imfilter(I_rec, gx, replicate)) abs(imfilter(I_rec, gy, replicate)); edge_bic abs(imfilter(I_bic, gx, replicate)) abs(imfilter(I_bic, gy, replicate)); % 边缘强度均值差值数值越小代表边缘形态越接近真值 delta_edge mean(abs(edge_cvx(:) - edge_bic(:)));拿到这几个数值后判断标准是PSNR 比双三次高至少 0.5 dB同时边缘差值显著下降。如果你跑出来 PSNR 高了但边缘差值反而变大多半是 λ 太小放大噪声。如果 PSNR 和双三次持平但视觉上更干净也说明方法本身有效只是指标没体现出来。这种情况我会再截一块纹理区放大对比主观判断在工程报告里同样要有据可依。5. 五个避坑记录从翻车里总结的 CVX 上采样经验5.1 报错“Disciplined convex programming error”不见得是你目标函数写错现象CVX 在解析阶段直接报错提示大意是“Cannot perform the operation {convex} ./ {convex}”。第一反应是目标函数选错了但问题是这段代码明明在别人那里能跑。原因CVX 只认它自己的原子函数组合。你把二次型写成了(D*x-y)*(D*x-y)或者把 L2 范数平方写成了norm(D*x-y,2)^2这两种写法在普通 Matlab 里合法在 CVX 里不被识别为凸表达。DCP 规则不允许直接对凸函数做任意幂次运算。解决二次型一律用sum_square(D*x-y)L1 项一律用norm(...,1)不要自己展开矩阵乘积。我还在目标里写过(D*x-y).^2这种逐元素平方同样报错换成sum_square立刻通过。记住一条原则CVX 的表达式只能用原子函数搭积木不能按数值计算的习惯展开。5.2 小图能跑大图直接内存爆炸现象64×64 的低分辨率图跑得挺顺利换到 256×256 的图CVX 在建模阶段就卡住后来直接提示内存不足。原因变量数等于高分辨率像素数。256×256 的低分辨率图放大 2 倍后有 512×512262144 个变量。CVX 处理 L1 项时会把问题转成一个大规模 SOCP内部变量还会成倍扩张内存消耗远不是稀疏矩阵表面看起来那么小。解决不要在整图上直接跑把输入切成 48×48 或 64×64 的低分辨率 patch逐块求解再拼接。如果业务不允许拼接就接受现实换成其他方法。CVX 适合做中小尺寸和原型验证不适合做大图整张重建这是方法论边界不是代码优化能解决的。5.3 重建图整体发灰亮度对不上现象求解完成后图像确实比双三次插值更清楚但整体发灰白色区域被压到了 0.8 左右黑色区域抬到 0.2 左右。原因最常见的是 D 矩阵能量不守恒。比如你用双三次下采样生成低分辨率观测但 D 矩阵却按块平均构造这会导致 D*x 的值域和 y 的值域有系统偏差。另一个原因是输入是 uint8 的 0 到 255直接放进了约束 0x255但 λ 还是按 0~1 范围设的正则项相对过强把像素往中间值拉。解决先执行full(sum(D,2))检查每一行确认非零和为 1。再统一数据范围输入图在 im2double 转成 0~1 后约束写 0x1λ 按 0.005~0.05 扫描。如果问题依旧把生成低分辨率观测的方式改成和 D 相同的 block 平均模型先做到自洽再谈真实退化。5.4 输出图像像马赛克拼接痕迹和网格状纹理明显现象重建结果边缘锐利但能隐约看到低分辨率像素块的网格像用积木拼出来的画。放大倍数越大网格越明显。原因这是 D 矩阵块平均模型带来的伪影。块平均相当于用一个矩形的低通滤波器它和现实中的光学 PSF 差别很大。当块边界和图像边缘恰好重合时重建结果会强化这种格子感。另一个诱因是 λ 过大平滑项把块内部的差异压平块之间又靠边缘检测拉回形成方块轮廓。解决把模糊矩阵 B 从单位阵换成高斯核对应的循环矩阵让退化模型变成“先模糊、后块平均”这样网格伪影会明显减弱。具体做法是用fspecial(gaussian,[5 5],1.0)生成核再用稀疏 Toeplitz 结构构造 B。如果暂时不想引入 B可以把 λ 调小到 0.01 量级观察网格是否缓解但这治标不治本。5.5 跑完比双三次还差PSNR 掉 0.5 dB 以上现象CVX 老老实实收敛了可拿到 PSNR 一看还不如imresize(...,bicubic)。求解时间花了十几分钟结果出来是负优化。原因大概率是正则项权重 λ 和噪声水平不匹配。λ 太小时重建会把下采样引入的噪声放大λ 太大时边缘和纹理被过度平滑丢失双三次还能保留的局部反差。另一个原因是观测图像本身已经被处理过比如你拿来的是一个已经增强或锐化过的低分辨率图它和 y DBxn 的假设不符凸优化没有任何高频信息可以恢复。解决先把 λ 在 0.005 到 0.05 之间做一次小扫描固定一张小图画出 PSNR 曲线。如果曲线最高点仍低于双三次那就检查退化模型。用合成数据进行闭环验证从一张清晰图生成 y再上采样确认流程能超过双三次然后再应用到真实图上。真实图若效果始终不理想说明该退化场景下凸优化的假设不成立放弃这条路不要硬调参。这个判断标准也回答了很多人的疑问凸优化上采样不是万能的但它值得做因为当你匹配好退化模型后效果是稳定且可解释的。6. 进阶用边缘自适应权重把 TV 正则局部化再把同一套框架迁移到去噪固定 λ 的 TV 正则有一个尴尬平滑区需要大 λ 来抑制噪声边缘附近需要小 λ 来保留对比度一个全局值无法同时满足。解决办法是给 TV 项加一个预计算的权重矩阵让每一点的正则强度不同。权重来自双三次插值结果的边缘检测图先算出一幅 0 到 1 的边缘响应图再映射成惩罚权重。x_init imresize(I_lr, s, bicubic); edge_map abs(imfilter(x_init, fspecial(sobel), replicate)); edge_map edge_map(:) / max(edge_map(:) eps); % 边缘响应大的地方权重小保留边缘平坦区权重大加强平滑 w 0.3 0.7 * (edge_map 0.35); W spdiags(w, 0, M*N, M*N); cvx_begin quiet variable x(M*N, 1) minimize( 0.5 * sum_square(D * x - y) ... lambda * norm(W * Tv * x, 1) ... lambda * norm(W * Th * x, 1) ) subject to x 0; x 1; cvx_end注意权重矩阵 W 必须是在进入 CVX 之前就固定好的常数矩阵不能用优化变量参与计算。如果你把 w 写成 x 的函数模型立刻违反 DCP 规则。我在第一次尝试时就是想让权重跟着迭代变化结果 CVX 报错这才改成用插值结果做预估计效果也不差。这个技巧让同一张图上的平坦区域和纹理区域不用再手工分成两个模型综合效果比固定 λ 好一个档。更重要的是这套框架可以整体迁移。把 D 换成单位矩阵低分辨率观测换成噪声图目标函数就变成了一个带 TV 正则的去噪模型。把 B 换成运动模糊的卷积矩阵它就变成了去模糊模型。我在做某跨平台系统的图像预处理模块时把同一个 CVX 求解函数复制了三份只改了 D 和 B 的构造部分就同时交付了上采样、去噪、去模糊三个功能。数学核心完全一致工程量只差在退化矩阵的构造上。最后说一个血泪教训我前几次跑出块状伪影时一直以为是自己 λ 没调好反复扫参扫了一个下午。后来才发现问题根本不在 λ而在 D 矩阵的构造和观测图像的降采样方式不一致。从那以后我再做任何上采样实验都坚持先造一对“已知真值”的合成数据做闭环验证再碰真实图像。希望这篇文章能帮你跳过这一段弯路直接看到凸优化上采样真正有效的样子。希望帮到你。本文还有配套的精品资源点击获取