简介Split Bregman算法是求解稀疏正则化反问题的经典迭代优化框架在图像去噪、压缩感知、医学成像等领域应用广泛。这份MATLAB源码实现了从初始化、数据拟合、正则化更新到Bregman重投影的完整流程适合算法学习者、图像处理研究者以及需要快速解决问题的工程师参考。压缩包共包含11个文件核心为.m脚本覆盖2D/3D全变分去噪、ART快速重建等模块并配有README说明、示例效果图jpg、实测数据mat以及License文件整体大小约20.4MB。源码严格遵循Split Bregman迭代思想主程序与子函数分离注释清晰并提供了可直接运行的Demo脚本运行后可得到重建效果对比图直观展示算法收敛行为和去噪/复原结果。使用者可替换输入数据或修改正则化参数快速扩展到超分辨率重建、MRI与CT图像恢复等任务也可参考源码结构将算法嵌入到自己的研究项目中。目前已有917人学习下载是一份适合教学演示与二次开发的高质量参考实现。 第一次在论文里看到 Split Bregman 这个名词时我下意识以为又是一套“数学上很优美、工程上很难落地”的方法。直到有次做图像去噪项目需要在各种 L1 正则化求解方案里挑一个能快速跑通、效果又稳的思路我才认认真真把 MATLAB 代码自己写了一遍。写完才发现核心实现加起来不过几十行比想象中简单太多。这篇文章要说的就是 Split Bregman 算法的 MATLAB 源代码实现。我会先讲清楚算法到底“拆”了什么再给出可以直接运行的 TV 图像去噪代码然后把我调参数、判断收敛和踩坑的实测经验一并放出来。不管你是做图像去噪、去模糊、压缩感知重建还是纯粹被 L1 范数优化问题的求解速度折磨过这篇内容应该都能帮上忙。1. Split Bregman 到底“拆”了什么从 TV 去噪的数学形式说起1.1 L1 范数为什么让传统优化方法头疼图像去噪里最经典的做法是求解一个 TV总变差正则化问题min_u ‖∇u‖₁ (μ/2)‖u − f‖₂²其中 f 是带噪观测图像u 是我们要恢复的干净图像∇u 是图像梯度‖·‖₁ 是 L1 范数。第一项鼓励重建结果在梯度域上稀疏也就是让图像尽量分片光滑第二项保证结果不能离观测太远。这个形式看着简单真正算起来却有个绕不开的麻烦L1 范数在零点不可导。你用梯度下降去最小化它迭代点一旦落到梯度零点附近次梯度方向会来回摇摆收敛非常慢。更麻烦的是梯度和 L1 还嵌套在一起直接求近端算子也得先处理这个复合结构。当年我一开始试着用最朴素的次梯度法跑这张图步长调小了吧几十次迭代看不到变化调大了吧能量直接震荡。后来接触到 Split Bregman 的核心思路我才意识到问题的解法其实很“拆”字诀别跟复合的 L1 硬碰硬先把里面的变量拆出来再逐个击破。1.2 变量分裂把一个难问题变成两个易问题Split Bregman 的第一板斧是变量分裂。引入两个辅助变量 dx 和 dy分别代表水平梯度和垂直梯度令 dx Dx·udy Dy·u。于是原始问题变成约束优化形式min_{u,dx,dy} ‖dx‖₁ ‖dy‖₁ (μ/2)‖u − f‖₂² s.t. dx Dx·u, dy Dy·u接下来把约束放进目标函数变成增广拉格朗日形式。这里引入罚参数 λ同时为了让约束严格成立再用 Bregman 变量 bx、by 做“记账”min_{u,dx,dy} ‖dx‖₁ ‖dy‖₁ (μ/2)‖u − f‖₂² (λ/2)(‖dx − Dx·u − bx‖₂² ‖dy − Dy·u − by‖₂²)为什么这么拆有用因为现在目标函数里L1 只跟 dx、dy 有关二次项只跟 u 的梯度有关两者不再嵌套。于是可以交替优化固定 dx、dy求解 u 的子问题。这个子问题是一个纯二次优化等价于解一个带拉普拉斯算子的线性方程组我在后面会讲为什么它能用 FFT 秒解。固定 u求解 dx、dy 的子问题。这个子问题退化成逐像素的软阈值收缩不需要解任何方程组一行代码搞定。原本一个啃不动的复合优化问题就这么被拆成了两块各自都有解析解的“软柿子”。1.3 Bregman 迭代误差是怎么“记账”并回填的有了变量分裂还不够如果直接拿罚函数方法迭代λ 有限时最终解和真实约束之间会有一个偏差。Bregman 迭代的精髓就是每轮把这种偏差“记下来”下一轮求解时再补回去。具体更新顺序是这样的用当前 dx、dy、bx、by 求解 u用更新后的 u 做软阈值更新 dx、dy更新 Bregman 变量 bx bx (Dx·u − dx) by by (Dy·u − dy)第 3 步相当于把“这一步梯度值和辅助变量之间没配平的差值”累积到 Bregman 变量里。下一次求解 u 时这个历史偏差会被带入右侧项迫使 u 逐步满足 dx Dx·u。你可以把它理解成罚函数法是每次只盯着当前的作业差Bregman 迭代则是把之前欠的账记在本子上下一次一起催收。正因为有这个“记账”机制即使 λ 取值有限算法迭代下去也能收敛到约束严格满足的最优解而不是一个被罚参数“惯坏”的近似解。到这里算法骨架已经很清楚了。接下来进入 MATLAB 代码层面看看这三板斧具体怎么落地。2. MATLAB 代码里的三块基石梯度算子、FFT 求解、软阈值2.1 离散梯度算子与 Neumann 边界条件写 Split Bregman 的 MATLAB 代码第一块基石是梯度算子。这里我选择用稀疏矩阵构造前向差分算子并用 Neumann反射边界条件也就是图像边界处梯度不跨出图像范围。对一幅 m×n 的图像 u 而言水平梯度作用于每一行所以水平差分矩阵 Gx 大小是 n×n垂直梯度作用于每一列所以垂直差分矩阵 Gy 大小是 m×m。用 MATLAB 构造这两块矩阵非常简洁Gx spdiags([-ones(n,1), ones(n,1)], [0 1], n, n); Gx(n,:) 0; % 最后一行置零边界不跨出图像 Gy spdiags([-ones(m,1), ones(m,1)], [0 1], m, m); Gy(m,:) 0;Gx 的每一行本质上是在做 “下一点减当前点” 的前向差分。最后一行置零对应的是图像最右侧像素没有右邻居梯度直接取 0。这样做的好处有两个一是矩阵天然稀疏乘法和转置乘法的计算量都很小二是 Neumann 边界条件比周期边界更贴近自然图像的实际情况不容易在图像边缘产生震荡伪影。有了这两个矩阵水平梯度就是 u * Gx垂直梯度就是 Gy * u。转置算子分别是 u * Gx 和 Gy * u。这里要注意转置方向不要搞反我见过不少人在这上面栽跟头后面踩坑章节会再提。2.2 u 子问题为什么 FFT 一次除法就搞定固定 dx、dy 后u 子问题是极小化(μ/2)‖u − f‖₂² (λ/2)(‖dx − Dx·u − bx‖₂² ‖dy − Dy·u − by‖₂²)对 u 求导并令其为零整理后得到线性方程(μI − λ(DxᵀDx DyᵀDy)) · u μf λ·Dxᵀ(dx − bx) λ·Dyᵀ(dy − by)这里 DxᵀDx DyᵀDy 是离散拉普拉斯算子。关键洞察是离散拉普拉斯在傅里叶基底下是对角的。也就是说这个看似复杂的线性系统在频域里变成了逐频点的除法。因此我们可以预计算分母kx (0:m-1) * pi / m; ky (0:n-1) * pi / n; LapEig (2 - 2*cos(kx)) (2 - 2*cos(ky)); denom mu lambda * LapEig;这里 LapEig 就是向量化的拉普拉斯特征值对应 Neumann 边界条件。有了分母之后每一轮迭代只需要两次 FFT、一次除法和一次逆 FFT复杂度是 O(N log N)大尺寸图像也毫无压力。2.3 d 子问题软阈值收缩的向量化写法固定 u 之后dx、dy 的更新是标准的 L1 近端算子问题。目标函数里关于 dx 的部分是‖dx‖₁ (λ/2)‖dx − (Dx·u bx)‖₂²这个问题的闭式解就是逐像素软阈值收缩dx max(|Dx·u bx| − 1/λ, 0) · sign(Dx·u bx)为什么阈值是 1/λ因为你把目标函数对每个像素单独拆开对一维问题 t (λ/2)(t − g)² 求极小值一阶条件会自然给出这个收缩形式。代码里用向量化写法gx u * Gx bx; thr 1 / lambda; dx max(abs(gx) - thr, 0) .* sign(gx);这里绝对不要用 for 循环加 if 逐像素判断既慢又容易出错MATLAB 的矩阵运算天然就是为这种操作设计的。3. 一份可直接复用的 Split Bregman TV 去噪代码3.1 主函数完整代码把前面三块基石拼起来就是完整的 Split Bregman TV 去噪函数。这段代码我实测过直接复制到脚本里就能跑兼容绝大多数 MATLAB 版本。function [u, energy] split_bregman_tv(f, mu, lambda, niter) % SPLIT_BREGMAN_TV Split Bregman 算法求解 TV 图像去噪 % 输入: % f - m×n 灰度图像double 类型建议归一化到 [0,1] % mu - 数据保真项权重控制去噪强度与细节保持的平衡 % lambda - 罚参数与 Bregman 迭代行为相关 % niter - 迭代次数 % 输出: % u - 去噪结果 % energy - 目标函数值历史可用来画收敛曲线 f double(f); [m, n] size(f); % ---------- 1. 构造离散梯度算子 (Neumann 边界) ---------- Gx spdiags([-ones(n,1), ones(n,1)], [0 1], n, n); Gx(n,:) 0; Gy spdiags([-ones(m,1), ones(m,1)], [0 1], m, m); Gy(m,:) 0; % ---------- 2. 初始化变量 ---------- u f; dx zeros(m, n); dy zeros(m, n); bx zeros(m, n); by zeros(m, n); % ---------- 3. 预计算频域分母 ---------- kx (0:m-1) * pi / m; ky (0:n-1) * pi / n; LapEig (2 - 2*cos(kx)) (2 - 2*cos(ky)); denom mu lambda * LapEig; Ff fft2(f); if nargout 1 energy zeros(1, niter); end % ---------- 4. 主迭代 ---------- for iter 1:niter % u 子问题频域除法求解 rhs mu * Ff lambda * fft2((dx - bx) * Gx Gy * (dy - by)); u real(ifft2(rhs ./ denom)); % d 子问题软阈值收缩 gx u * Gx bx; gy Gy * u by; thr 1 / lambda; dx max(abs(gx) - thr, 0) .* sign(gx); dy max(abs(gy) - thr, 0) .* sign(gy); % Bregman 变量更新 bx bx (u * Gx - dx); by by (Gy * u - dy); % 记录目标函数值 if nargout 1 energy(iter) sum(abs(dx(:)) abs(dy(:))) ... 0.5 * mu * sum((u(:) - f(:)).^2); end end end代码结构很简单前面是算子和分母的预计算中间是主迭代核心循环就十行左右。这也正是 Split Bregman 魅力所在把算法思想翻译成 MATLAB 代码时几乎没有多余的中间层。3.2 怎么调用它一段可跑通的演示脚本写代码容易跑通才是硬道理。下面这段演示脚本可以直接用 MATLAB 自带的 cameraman.tif 测试%% 加载图像并加噪声 img im2double(imread(cameraman.tif)); rng(42); noisy img 0.05 * randn(size(img)); noisy max(min(noisy, 1), 0); %% 调用 Split Bregman TV 去噪 [denoised, energy] split_bregman_tv(noisy, 20, 40, 50); %% 显示结果 figure; subplot(1,3,1); imshow(img); title(原图); subplot(1,3,2); imshow(noisy); title(带噪图); subplot(1,3,3); imshow(denoised); title(去噪结果); %% 画能量下降曲线 figure; semilogy(max(energy - energy(end), eps)); xlabel(迭代次数); ylabel(目标函数值); title(能量下降曲线);如果手头没有 cameraman.tif用 phantom、coins 这类内置图像也可以关键是把图像归一化到 [0,1] 之后mu 和 lambda 的经验值才更可移植。3.3 代码中几个容易写错的地方主函数虽然短但我在调试过程中发现几个特别容易出错的点单独拿出来说一下。第一个是算子转置方向。(dx - bx) * Gx是 Dxᵀ(dx − bx)u * Gx才是 Dx·u这两个方向完全不一样。我最初写的时候把某个方向搞反了结果去噪效果一直不对劲图像看起来像是被“推”着往右下角移动。第二个是 FFT 的坑。求解 u 子问题时右侧项里有 fft2左侧频域除法之后要记得用 real 取实部。理论上 ifft2 的结果虚部应该接近 0但数值误差会让虚部残留一点噪声不取 real 直接显示图像会出现奇怪的“水波纹”。第三个是软阈值不要写成gx(gx thr) 0这种形式。软阈值是先减去阈值再置零不是直接小于阈值就扔掉。写成max(abs(gx) - thr, 0) .* sign(gx)是标准做法既保证符号正确又保持向量化。注意输入图像一定要先转成 double。uint8 下的减法会直接截断成负数结果会非常离谱。这是我第一次跑这段代码时踩的坑当时输出图像全是黑的排查了半天才发现是类型问题。4. mu、lambda 与迭代次数我的调参与收敛经验4.1 从能量函数理解 mu 与 lambda 的分工参数怎么调是几乎所有接触 Split Bregman 的人都会问的问题。我的经验是先理解两个参数的分工mu 是数据保真项权重决定重建结果“贴近观测图像”和“平滑去噪”之间的最终平衡。mu 越大结果越接近带噪原图噪声残留越多mu 越小图像越平滑但细节也越容易丢失。lambda 是分裂罚参数主要影响 Bregman 迭代的收敛行为和中间路径。它不像 mu 那样直接决定最终结果的平滑度但 lambda 太小会让软阈值收缩过强每轮迭代步子迈得太大lambda 太大会让收缩接近恒等映射收敛变慢。关键体会是mu 和 lambda 的绝对大小不如它们的相对关系重要。实际使用中lambda 通常取 mu 的 2 到 5 倍。这样设置可以保证 u 子问题和 d 子问题的迭代交替时两边的“权重”不会一极独大。用一张表来概括参数偏大偏小的影响参数取值偏小取值偏大mu去噪过度边缘与细节模糊噪声残留明显去噪不充分lambda软阈值收缩过强图像发虚Bregman 迭代收敛变慢模板效果不明显niter未充分收敛能量曲线仍在下降浪费算力对结果几乎无改变4.2 我实测的参数区间与迭代次数对于灰度图像归一化到 [0,1]、高斯噪声标准差在 0.05 到 0.1 之间的场景我实测下来比较好的参数范围是mu 取 10 到 30lambda 取 30 到 60niter 取 30 到 60。如果用 PSNR 做评价指标mu 调到区间中间值附近通常就能拿到不错的结果。下面是几个我实际试过的组合噪声标准差约为 0.05 时的去噪后 PSNR% 测试不同 mu 的影响 mu_list [5, 10, 20, 40, 80]; for i 1:length(mu_list) u split_bregman_tv(noisy, mu_list(i), 40, 50); psnr_val(i) 10*log10(1 / mean((u(:)-img(:)).^2)); end实测结果大致是 mu5 时 PSNR 较低图像偏平滑mu20 附近 PSNR 最高mu80 时 PSNR 又下降因为噪声没有被有效去除。这个规律在大多数测试图上都成立。4.3 用能量曲线判断收敛而不是拍脑袋定次数迭代次数不是越大越好也不是固定的。判断收敛最靠谱的方式是看能量函数曲线的下降情况。我在函数里返回了 energy就是为了方便你画这条曲线。实践中的做法是跑 60 次迭代画semilogy(energy)或者semilogy(abs(diff(energy)) eps)看后半段的斜率。如果相邻迭代的能量差已经降到初始差值的千分之一以下基本可以认为收敛了。对于大多数图像30 到 50 次迭代已经足够继续跑更多轮也只是微调。如果发现能量曲线在最后还在明显下降说明 niter 不够如果曲线早就平了但结果还不满意那问题大概率不是迭代次数而是 mu 没有调对。这条区分逻辑能帮你少走很多弯路。5. 踩坑记录与扩展思路从去噪走向其他逆问题5.1 我在复现过程中踩过的坑除了前面提到的类型转换和算子方向还有一些更深层的坑值得记录。第一个坑是边界条件不一致导致的伪影。如果梯度算子用了 Neumann 边界但 FFT 求解时却用周期边界公式计算特征值图像边缘就会出现明显的明暗条纹。解决方法是确保两处边界条件一致。我代码里 Neum本文还有配套的精品资源点击获取 
长河君 · 2026/9/8 22:22:22 · 阅读 1.2k Split Bregman算法MATLAB实现：图像去噪实战与代码解析
简介Split Bregman算法是求解稀疏正则化反问题的经典迭代优化框架在图像去噪、压缩感知、医学成像等领域应用广泛。这份MATLAB源码实现了从初始化、数据拟合、正则化更新到Bregman重投影的完整流程适合算法学习者、图像处理研究者以及需要快速解决问题的工程师参考。压缩包共包含11个文件核心为.m脚本覆盖2D/3D全变分去噪、ART快速重建等模块并配有README说明、示例效果图jpg、实测数据mat以及License文件整体大小约20.4MB。源码严格遵循Split Bregman迭代思想主程序与子函数分离注释清晰并提供了可直接运行的Demo脚本运行后可得到重建效果对比图直观展示算法收敛行为和去噪/复原结果。使用者可替换输入数据或修改正则化参数快速扩展到超分辨率重建、MRI与CT图像恢复等任务也可参考源码结构将算法嵌入到自己的研究项目中。目前已有917人学习下载是一份适合教学演示与二次开发的高质量参考实现。 第一次在论文里看到 Split Bregman 这个名词时我下意识以为又是一套“数学上很优美、工程上很难落地”的方法。直到有次做图像去噪项目需要在各种 L1 正则化求解方案里挑一个能快速跑通、效果又稳的思路我才认认真真把 MATLAB 代码自己写了一遍。写完才发现核心实现加起来不过几十行比想象中简单太多。这篇文章要说的就是 Split Bregman 算法的 MATLAB 源代码实现。我会先讲清楚算法到底“拆”了什么再给出可以直接运行的 TV 图像去噪代码然后把我调参数、判断收敛和踩坑的实测经验一并放出来。不管你是做图像去噪、去模糊、压缩感知重建还是纯粹被 L1 范数优化问题的求解速度折磨过这篇内容应该都能帮上忙。1. Split Bregman 到底“拆”了什么从 TV 去噪的数学形式说起1.1 L1 范数为什么让传统优化方法头疼图像去噪里最经典的做法是求解一个 TV总变差正则化问题min_u ‖∇u‖₁ (μ/2)‖u − f‖₂²其中 f 是带噪观测图像u 是我们要恢复的干净图像∇u 是图像梯度‖·‖₁ 是 L1 范数。第一项鼓励重建结果在梯度域上稀疏也就是让图像尽量分片光滑第二项保证结果不能离观测太远。这个形式看着简单真正算起来却有个绕不开的麻烦L1 范数在零点不可导。你用梯度下降去最小化它迭代点一旦落到梯度零点附近次梯度方向会来回摇摆收敛非常慢。更麻烦的是梯度和 L1 还嵌套在一起直接求近端算子也得先处理这个复合结构。当年我一开始试着用最朴素的次梯度法跑这张图步长调小了吧几十次迭代看不到变化调大了吧能量直接震荡。后来接触到 Split Bregman 的核心思路我才意识到问题的解法其实很“拆”字诀别跟复合的 L1 硬碰硬先把里面的变量拆出来再逐个击破。1.2 变量分裂把一个难问题变成两个易问题Split Bregman 的第一板斧是变量分裂。引入两个辅助变量 dx 和 dy分别代表水平梯度和垂直梯度令 dx Dx·udy Dy·u。于是原始问题变成约束优化形式min_{u,dx,dy} ‖dx‖₁ ‖dy‖₁ (μ/2)‖u − f‖₂² s.t. dx Dx·u, dy Dy·u接下来把约束放进目标函数变成增广拉格朗日形式。这里引入罚参数 λ同时为了让约束严格成立再用 Bregman 变量 bx、by 做“记账”min_{u,dx,dy} ‖dx‖₁ ‖dy‖₁ (μ/2)‖u − f‖₂² (λ/2)(‖dx − Dx·u − bx‖₂² ‖dy − Dy·u − by‖₂²)为什么这么拆有用因为现在目标函数里L1 只跟 dx、dy 有关二次项只跟 u 的梯度有关两者不再嵌套。于是可以交替优化固定 dx、dy求解 u 的子问题。这个子问题是一个纯二次优化等价于解一个带拉普拉斯算子的线性方程组我在后面会讲为什么它能用 FFT 秒解。固定 u求解 dx、dy 的子问题。这个子问题退化成逐像素的软阈值收缩不需要解任何方程组一行代码搞定。原本一个啃不动的复合优化问题就这么被拆成了两块各自都有解析解的“软柿子”。1.3 Bregman 迭代误差是怎么“记账”并回填的有了变量分裂还不够如果直接拿罚函数方法迭代λ 有限时最终解和真实约束之间会有一个偏差。Bregman 迭代的精髓就是每轮把这种偏差“记下来”下一轮求解时再补回去。具体更新顺序是这样的用当前 dx、dy、bx、by 求解 u用更新后的 u 做软阈值更新 dx、dy更新 Bregman 变量 bx bx (Dx·u − dx) by by (Dy·u − dy)第 3 步相当于把“这一步梯度值和辅助变量之间没配平的差值”累积到 Bregman 变量里。下一次求解 u 时这个历史偏差会被带入右侧项迫使 u 逐步满足 dx Dx·u。你可以把它理解成罚函数法是每次只盯着当前的作业差Bregman 迭代则是把之前欠的账记在本子上下一次一起催收。正因为有这个“记账”机制即使 λ 取值有限算法迭代下去也能收敛到约束严格满足的最优解而不是一个被罚参数“惯坏”的近似解。到这里算法骨架已经很清楚了。接下来进入 MATLAB 代码层面看看这三板斧具体怎么落地。2. MATLAB 代码里的三块基石梯度算子、FFT 求解、软阈值2.1 离散梯度算子与 Neumann 边界条件写 Split Bregman 的 MATLAB 代码第一块基石是梯度算子。这里我选择用稀疏矩阵构造前向差分算子并用 Neumann反射边界条件也就是图像边界处梯度不跨出图像范围。对一幅 m×n 的图像 u 而言水平梯度作用于每一行所以水平差分矩阵 Gx 大小是 n×n垂直梯度作用于每一列所以垂直差分矩阵 Gy 大小是 m×m。用 MATLAB 构造这两块矩阵非常简洁Gx spdiags([-ones(n,1), ones(n,1)], [0 1], n, n); Gx(n,:) 0; % 最后一行置零边界不跨出图像 Gy spdiags([-ones(m,1), ones(m,1)], [0 1], m, m); Gy(m,:) 0;Gx 的每一行本质上是在做 “下一点减当前点” 的前向差分。最后一行置零对应的是图像最右侧像素没有右邻居梯度直接取 0。这样做的好处有两个一是矩阵天然稀疏乘法和转置乘法的计算量都很小二是 Neumann 边界条件比周期边界更贴近自然图像的实际情况不容易在图像边缘产生震荡伪影。有了这两个矩阵水平梯度就是 u * Gx垂直梯度就是 Gy * u。转置算子分别是 u * Gx 和 Gy * u。这里要注意转置方向不要搞反我见过不少人在这上面栽跟头后面踩坑章节会再提。2.2 u 子问题为什么 FFT 一次除法就搞定固定 dx、dy 后u 子问题是极小化(μ/2)‖u − f‖₂² (λ/2)(‖dx − Dx·u − bx‖₂² ‖dy − Dy·u − by‖₂²)对 u 求导并令其为零整理后得到线性方程(μI − λ(DxᵀDx DyᵀDy)) · u μf λ·Dxᵀ(dx − bx) λ·Dyᵀ(dy − by)这里 DxᵀDx DyᵀDy 是离散拉普拉斯算子。关键洞察是离散拉普拉斯在傅里叶基底下是对角的。也就是说这个看似复杂的线性系统在频域里变成了逐频点的除法。因此我们可以预计算分母kx (0:m-1) * pi / m; ky (0:n-1) * pi / n; LapEig (2 - 2*cos(kx)) (2 - 2*cos(ky)); denom mu lambda * LapEig;这里 LapEig 就是向量化的拉普拉斯特征值对应 Neumann 边界条件。有了分母之后每一轮迭代只需要两次 FFT、一次除法和一次逆 FFT复杂度是 O(N log N)大尺寸图像也毫无压力。2.3 d 子问题软阈值收缩的向量化写法固定 u 之后dx、dy 的更新是标准的 L1 近端算子问题。目标函数里关于 dx 的部分是‖dx‖₁ (λ/2)‖dx − (Dx·u bx)‖₂²这个问题的闭式解就是逐像素软阈值收缩dx max(|Dx·u bx| − 1/λ, 0) · sign(Dx·u bx)为什么阈值是 1/λ因为你把目标函数对每个像素单独拆开对一维问题 t (λ/2)(t − g)² 求极小值一阶条件会自然给出这个收缩形式。代码里用向量化写法gx u * Gx bx; thr 1 / lambda; dx max(abs(gx) - thr, 0) .* sign(gx);这里绝对不要用 for 循环加 if 逐像素判断既慢又容易出错MATLAB 的矩阵运算天然就是为这种操作设计的。3. 一份可直接复用的 Split Bregman TV 去噪代码3.1 主函数完整代码把前面三块基石拼起来就是完整的 Split Bregman TV 去噪函数。这段代码我实测过直接复制到脚本里就能跑兼容绝大多数 MATLAB 版本。function [u, energy] split_bregman_tv(f, mu, lambda, niter) % SPLIT_BREGMAN_TV Split Bregman 算法求解 TV 图像去噪 % 输入: % f - m×n 灰度图像double 类型建议归一化到 [0,1] % mu - 数据保真项权重控制去噪强度与细节保持的平衡 % lambda - 罚参数与 Bregman 迭代行为相关 % niter - 迭代次数 % 输出: % u - 去噪结果 % energy - 目标函数值历史可用来画收敛曲线 f double(f); [m, n] size(f); % ---------- 1. 构造离散梯度算子 (Neumann 边界) ---------- Gx spdiags([-ones(n,1), ones(n,1)], [0 1], n, n); Gx(n,:) 0; Gy spdiags([-ones(m,1), ones(m,1)], [0 1], m, m); Gy(m,:) 0; % ---------- 2. 初始化变量 ---------- u f; dx zeros(m, n); dy zeros(m, n); bx zeros(m, n); by zeros(m, n); % ---------- 3. 预计算频域分母 ---------- kx (0:m-1) * pi / m; ky (0:n-1) * pi / n; LapEig (2 - 2*cos(kx)) (2 - 2*cos(ky)); denom mu lambda * LapEig; Ff fft2(f); if nargout 1 energy zeros(1, niter); end % ---------- 4. 主迭代 ---------- for iter 1:niter % u 子问题频域除法求解 rhs mu * Ff lambda * fft2((dx - bx) * Gx Gy * (dy - by)); u real(ifft2(rhs ./ denom)); % d 子问题软阈值收缩 gx u * Gx bx; gy Gy * u by; thr 1 / lambda; dx max(abs(gx) - thr, 0) .* sign(gx); dy max(abs(gy) - thr, 0) .* sign(gy); % Bregman 变量更新 bx bx (u * Gx - dx); by by (Gy * u - dy); % 记录目标函数值 if nargout 1 energy(iter) sum(abs(dx(:)) abs(dy(:))) ... 0.5 * mu * sum((u(:) - f(:)).^2); end end end代码结构很简单前面是算子和分母的预计算中间是主迭代核心循环就十行左右。这也正是 Split Bregman 魅力所在把算法思想翻译成 MATLAB 代码时几乎没有多余的中间层。3.2 怎么调用它一段可跑通的演示脚本写代码容易跑通才是硬道理。下面这段演示脚本可以直接用 MATLAB 自带的 cameraman.tif 测试%% 加载图像并加噪声 img im2double(imread(cameraman.tif)); rng(42); noisy img 0.05 * randn(size(img)); noisy max(min(noisy, 1), 0); %% 调用 Split Bregman TV 去噪 [denoised, energy] split_bregman_tv(noisy, 20, 40, 50); %% 显示结果 figure; subplot(1,3,1); imshow(img); title(原图); subplot(1,3,2); imshow(noisy); title(带噪图); subplot(1,3,3); imshow(denoised); title(去噪结果); %% 画能量下降曲线 figure; semilogy(max(energy - energy(end), eps)); xlabel(迭代次数); ylabel(目标函数值); title(能量下降曲线);如果手头没有 cameraman.tif用 phantom、coins 这类内置图像也可以关键是把图像归一化到 [0,1] 之后mu 和 lambda 的经验值才更可移植。3.3 代码中几个容易写错的地方主函数虽然短但我在调试过程中发现几个特别容易出错的点单独拿出来说一下。第一个是算子转置方向。(dx - bx) * Gx是 Dxᵀ(dx − bx)u * Gx才是 Dx·u这两个方向完全不一样。我最初写的时候把某个方向搞反了结果去噪效果一直不对劲图像看起来像是被“推”着往右下角移动。第二个是 FFT 的坑。求解 u 子问题时右侧项里有 fft2左侧频域除法之后要记得用 real 取实部。理论上 ifft2 的结果虚部应该接近 0但数值误差会让虚部残留一点噪声不取 real 直接显示图像会出现奇怪的“水波纹”。第三个是软阈值不要写成gx(gx thr) 0这种形式。软阈值是先减去阈值再置零不是直接小于阈值就扔掉。写成max(abs(gx) - thr, 0) .* sign(gx)是标准做法既保证符号正确又保持向量化。注意输入图像一定要先转成 double。uint8 下的减法会直接截断成负数结果会非常离谱。这是我第一次跑这段代码时踩的坑当时输出图像全是黑的排查了半天才发现是类型问题。4. mu、lambda 与迭代次数我的调参与收敛经验4.1 从能量函数理解 mu 与 lambda 的分工参数怎么调是几乎所有接触 Split Bregman 的人都会问的问题。我的经验是先理解两个参数的分工mu 是数据保真项权重决定重建结果“贴近观测图像”和“平滑去噪”之间的最终平衡。mu 越大结果越接近带噪原图噪声残留越多mu 越小图像越平滑但细节也越容易丢失。lambda 是分裂罚参数主要影响 Bregman 迭代的收敛行为和中间路径。它不像 mu 那样直接决定最终结果的平滑度但 lambda 太小会让软阈值收缩过强每轮迭代步子迈得太大lambda 太大会让收缩接近恒等映射收敛变慢。关键体会是mu 和 lambda 的绝对大小不如它们的相对关系重要。实际使用中lambda 通常取 mu 的 2 到 5 倍。这样设置可以保证 u 子问题和 d 子问题的迭代交替时两边的“权重”不会一极独大。用一张表来概括参数偏大偏小的影响参数取值偏小取值偏大mu去噪过度边缘与细节模糊噪声残留明显去噪不充分lambda软阈值收缩过强图像发虚Bregman 迭代收敛变慢模板效果不明显niter未充分收敛能量曲线仍在下降浪费算力对结果几乎无改变4.2 我实测的参数区间与迭代次数对于灰度图像归一化到 [0,1]、高斯噪声标准差在 0.05 到 0.1 之间的场景我实测下来比较好的参数范围是mu 取 10 到 30lambda 取 30 到 60niter 取 30 到 60。如果用 PSNR 做评价指标mu 调到区间中间值附近通常就能拿到不错的结果。下面是几个我实际试过的组合噪声标准差约为 0.05 时的去噪后 PSNR% 测试不同 mu 的影响 mu_list [5, 10, 20, 40, 80]; for i 1:length(mu_list) u split_bregman_tv(noisy, mu_list(i), 40, 50); psnr_val(i) 10*log10(1 / mean((u(:)-img(:)).^2)); end实测结果大致是 mu5 时 PSNR 较低图像偏平滑mu20 附近 PSNR 最高mu80 时 PSNR 又下降因为噪声没有被有效去除。这个规律在大多数测试图上都成立。4.3 用能量曲线判断收敛而不是拍脑袋定次数迭代次数不是越大越好也不是固定的。判断收敛最靠谱的方式是看能量函数曲线的下降情况。我在函数里返回了 energy就是为了方便你画这条曲线。实践中的做法是跑 60 次迭代画semilogy(energy)或者semilogy(abs(diff(energy)) eps)看后半段的斜率。如果相邻迭代的能量差已经降到初始差值的千分之一以下基本可以认为收敛了。对于大多数图像30 到 50 次迭代已经足够继续跑更多轮也只是微调。如果发现能量曲线在最后还在明显下降说明 niter 不够如果曲线早就平了但结果还不满意那问题大概率不是迭代次数而是 mu 没有调对。这条区分逻辑能帮你少走很多弯路。5. 踩坑记录与扩展思路从去噪走向其他逆问题5.1 我在复现过程中踩过的坑除了前面提到的类型转换和算子方向还有一些更深层的坑值得记录。第一个坑是边界条件不一致导致的伪影。如果梯度算子用了 Neumann 边界但 FFT 求解时却用周期边界公式计算特征值图像边缘就会出现明显的明暗条纹。解决方法是确保两处边界条件一致。我代码里 Neum本文还有配套的精品资源点击获取