简介：Split Bregman 算法由 Goldenberg 与 Baraniuk 于 2009 年提出，是一种将复杂优化问题分解为数据拟合与正则化两个子问题的迭代优化方法，在图像去噪、图像复原、压缩感知和医学成像中应用广泛。MATLAB 源码实现共包含 11 个文件，以 6 个 .m 脚本和函数为主，覆盖二维/三维全变分去噪、ART 重建、快速算法版和演示程序，并附带 2 张算法公式示意图、1 份说明文档、1 个测试数据 .mat 文件及许可文件，压缩包整体约 20.4MB。代码结构按算法步骤组织，核心脚本涵盖 ART 重建、快速 ART 重建、二维/三维全变分去噪以及结果可视化工具，可直接运行演示程序观察 Split Bregman 迭代重建效果。已有 918 人学习下载。借助测试数据与公式示意图，读者既能对照公式推演 Split Bregman 迭代步骤，也能通过实际数据验证算法收敛性与去噪效果，适合图像处理、信号处理、压缩感知与医学成像等方向的研究生和工程师参考、修改并迁移到自己的项目中。

1. Split Bregman 算法：为什么 MATLAB 源码里十几行循环就能收敛

第一次在 MATLAB 里跑 Split Bregman 去噪时，我把主循环摊开在命令行逐行读，发现所谓新算法其实就是“交替解一个线性方程、做一次软阈值收缩、再更新一次辅助罚项”三个动作循环。真正费时间的不是理解算法，而是梯度算子、收缩算子和罚项符号这些细节——参数设不好，迭代几百次照样不收敛。

这份源码解决的是全变分（TV）正则项下的图像去噪、去模糊和压缩感知重建问题，核心价值是把原来需要外部凸优化工具箱的难问题压缩成几十行纯 MATLAB 代码。适合正在做图像去噪实验、论文里需要对比 TV 方法，或者想理解 L1 优化到底怎么落地的人。

如果你是第一次接触 bregman 迭代，建议先跑通 demo 再回头读原理；如果是老手，直接跳到最后两章看参数和残差曲线的坑。

2. 从 Bregman 迭代到变量分裂：源码的三个核心模块拆解

2.1 全变分模型为什么难解

TV 去噪的目标函数可以写成这样：最小化(λ/2)||u - f||² + ||D_x u||₁ + ||D_y u||₁。第一项要求重建图 u 不要偏离观测图 f 太远，第二、三项要求 u 的梯度尽量稀疏。稀疏梯度意味着图像内部平滑、边缘处有大的跳变，这正是自然图像的结构特征。

这个模型看着简单，真要数值求解时处处别扭。L1 范数在零点的不可导性，让普通梯度下降走到平坦区域附近时来回抖动，步长调小了像蜗牛爬，调大了直接振荡。我最早用 ISTA（迭代软阈值）硬解过，收敛速度慢得让人失去耐心；后来换成 Bregman 迭代，收敛速度明显上了一个台阶。

Split Bregman 的精髓不在“Bregman”这个名字唬人，而是把原问题里不可导的||D u||₁分裂成对辅助变量 d 的惩罚项：先解光滑的 u 子问题，再用软阈值闭式解 d 子问题，最后用 b 变量把两步之间的误差“记回来”。这套思路绕开了 L1 不可导的雷区，每一步都是现成的 MATLAB 操作。

2.2 变量分裂：把不可导项挪到辅助变量上

具体做法是引入d_x、d_y两个辅助变量，让它们去逼近D_x u和D_y u。原问题变成对 u 和 d 交替优化：固定 d 更新 u 时，目标函数是二次型，有闭式解；固定 u 更新 d 时，等价于对残差做一次软阈值收缩。麻烦在于两个子问题来回切换会累积误差，所以又引入 b 变量作为“后悔药”，把每次分裂产生的误差留到下一步补偿。

整个循环里反复出现的符号，先把它们捋清楚，看源码时不至于迷路。

u 更新时，d 和 b 固定，问题变成标准的最小二乘，可以用 FFT 一步解出；d 更新时，u 和 b 固定，变成逐点软阈值。b 的更新在每一轮最后执行，它记住的是“当前差了多少”，下一轮把差值补偿回 u 的求解里。这就是 Bregman 迭代“不走回头路但保留误差记忆”的核心。

2.3 源码里的三个算子与主循环：可以直接跑的 MATLAB 骨架

先看梯度算子和它的伴随算子。前向差分算子负责把图像变成梯度场，伴随算子负责把梯度场映回图像域，对应负散度。两者配合 FFT 才能在频域一步求解 u 子问题。

% D.m 前向差分梯度算子，周期边界 % 输入 u：M*N 灰度图；输出两个 M*N 差分矩阵 function [Dx_u, Dy_u] = D(u) Dx_u = [diff(u, 1, 2), u(:, 1) - u(:, end)]; Dy_u = [diff(u, 1, 1); u(1, :) - u(end, :)]; end % 伴随算子（周期边界下）水平分量 function Dtx_d = Dtx(d) Dtx_d = circshift(d, [0 1]) - d; end % 伴随算子（周期边界下）垂直分量 function Dty_d = Dty(d) Dty_d = circshift(d, [1 0]) - d; end

逻辑说明：diff(u,1,2)沿列方向做前向差分，得到 M×(N-1)，补上最后一列形成周期边界闭环，保证输出和 u 同尺寸。伴随算子里的circshift是实现转置差分的便捷写法，circshift(d,[0 1])把矩阵沿列方向整体右移一格，减掉自己就得到D^T的效果。

参数说明：d 必须和 u 同尺寸；周期边界意味着图像左右、上下被当成首尾相连，这为后面的 FFT 求解提供了前提。

再看软阈值收缩算子，这是整个算法里唯一出现非线性操作的地方。

% 软阈值（近端算子）：x 是残差，gamma 是阈值 function y = shrink(x, gamma) y = max(abs(x) - gamma, 0) .* sign(x); end

逻辑说明：绝对值小于 gamma 的系数直接置零，大于 gamma 的向零点方向收缩 gamma。这个操作对所有元素并行，没有任何循环，MATLAB 跑起来非常快。

参数说明：gamma 在去噪场景取1/beta。beta 越大，阈值越小，d 越接近 D(u)，残差保留得越多。

最后是主循环。建议把下面这段封装成split_bregman_denoise.m函数，后续所有实验都调用它。

% 预计算频域分母：利用周期边界 + FFT 对角化 [M, N] = size(f); [Kx, Ky] = meshgrid(0:N-1, 0:M-1); % 注意：Kx 的每行是列频率 Fx = exp(1i * 2 * pi * Kx / N) - 1; % 水平差分特征值 Fy = exp(1i * 2 * pi * Ky / M) - 1; % 垂直差分特征值 denom = lambda + beta * (abs(Fx).^2 + abs(Fy).^2); % 初始化辅助变量 dx = zeros(size(f)); dy = dx; bx = dx; by = dx; for k = 1:maxIter % 第 1 步：解 u 的线性子问题，FFT 一步完成 rhs = lambda * f + beta * (Dtx(dx - bx) + Dty(dy - by)); u = real(ifft2(fft2(rhs) ./ denom)); % 第 2 步：软阈值收缩，更新辅助变量 d vx = Dx(u) + bx; vy = Dy(u) + by; dx = shrink(vx, 1 / beta); dy = shrink(vy, 1 / beta); % 第 3 步：更新 Bregman 罚项 b bx = bx + (Dx(u) - dx); by = by + (Dy(u) - by); end

逻辑说明：第一步求解的是带罚项的二次型，denom提前算好避免每轮重复构造。第二步把D(u)+b送入软阈值，得到新的 d。第三步 b 累加的是“当前梯度与辅助变量之间的差值”，这个差值恰恰是下一轮要补偿的信息，也是 Bregman 迭代区别于普通罚方法的根本。

参数说明：lambda建议从 0.05 起步，beta从 0.5 起步，maxIter先给 80。初始化时 u 可以直接用 f，d 和 b 全零即可。如果输入是 RGB 图，先转灰度，否则尺寸不匹配会让diff作用到三维矩阵上，结果完全不是想要的。

3. 跑通去噪流程：参数设置与效果验证的五个关键点

3.1 一键式 Demo：从读图到出图的完整脚本

拿到源码包后，先别急着改算法，把下面这个 demo 完整跑通，确认环境没问题再动手。

% demo_run.m 跑通 Split Bregman TV 去噪 clear; close all; % 1. 读图并转灰度 double，范围 [0,1] f = im2double(imread('cameraman.tif')); if size(f, 3) == 3, f = rgb2gray(f); end % 2. 加高斯白噪声，固定随机种子方便复现 rng(42); g = f + 0.05 * randn(size(f)); % 3. 设置参数：lambda 控制保真，beta 控制分裂罚 lambda = 0.05; beta = 0.5; maxIter = 80; % 4. 调用封装好的主循环 u = split_bregman_denoise(g, lambda, beta, maxIter); % 5. 输出 PSNR fprintf('PSNR in = %.2f dB, out = %.2f dB\n', ... my_psnr(g, f), my_psnr(u, f)); % 6. 三张图对比 subplot(1,3,1); imshow(f); title('原图'); subplot(1,3,2); imshow(g); title('含噪'); subplot(1,3,3); imshow(u); title('重建'); % 如果 Image Processing Toolbox 没有 psnr 函数，用这个替代 function p = my_psnr(a, b) mse = mean((a(:) - b(:)).^2); p = 10 * log10(1 / max(mse, eps)); end

逻辑说明：第 1 步用im2double把 uint8 图像转成 double 且归一化到 [0,1]，这一步不做的话后面全白或 NaN。第 2 步加标准差 0.05 的高斯噪声，固定随机种子让每次结果可复现。第 6 步的psnr函数在部分旧版 MATLAB 里需要 Image Processing Toolbox，没有就用手写函数替换。

参数说明：demo 里的 lambda、beta、maxIter 是保守起点，跑完应该能看到 PSNR 提升 2~4 dB。如果输出比你预期差很多，先检查是否用了 RGB 图直接进主循环。

3.2 参数经验表：lambda、beta、maxIter 怎么配

这套算法的参数集中且敏感，先把经验范围摆在这里，再逐项解释。

lambda 的物理含义是“重建图对含噪图的信任程度”。噪声强时，观测值不可信，lambda 要调小，让正则项占据主导；如果图像本身纹理丰富，lambda 调太小会把细节当噪声抹掉，此时应该适当加大。

beta 控制 d 与 D(u) 的绑定紧密度。beta 偏小时，软阈值幅度大，每步改动 d 的步子小，收敛慢；beta 偏大时，每步 d 变动剧烈，能量曲线可能出现锯齿状震荡。我一般先固定 beta=0.5，把 lambda 调到视觉满意，再回头微调 beta。

maxIter 不是越多越好。很多实现默认给 200，但实际跑到 60 次左右曲线已经走平。后面第 6 章会讲怎么用能量曲线判断该停在哪，这里先记住一个原则：与其盲目加大迭代次数，不如先确认参数组合正确。

3.3 看结果：量化指标与视觉检查

PSNR 是去噪实验最常用的量化指标，但只靠它判断效果会翻车。PSNR 对“整体像素误差”敏感，对“边缘是否锐利”不敏感，一张被磨皮的图 PSNR 可能不低，视觉上却很差。

正确的验收姿势是量化指标加视觉检查一起看。视觉检查三件事：边缘是否保持锐利；平滑区域是否出现块状或网格；重建图像有没有横竖条纹。如果边缘发虚，多半是各向同性/各向异性的选择问题（见第 4 章）；如果有条纹，先检查有没有踩第 5 章的 meshgrid 坑。

调参到这个阶段就是玄学现场：同一组参数在不同图像上表现完全不同，没有万能组合。我的习惯是准备两张图，一张纹理丰富、一张大片平滑区，参数组合必须同时在这两张图上表现稳定才算合格。

4. 场景改造：从各向异性 TV 到压缩感知与彩色图

4.1 各向同性 TV：把两个 shrink 合成一个联合收缩

第 2 章给出的主循环是各向异性 TV，即水平方向和垂直方向独立做软阈值。实际使用中，各向同性 TV 更常见，因为它在斜向边缘上响应一致，不会出现横竖分家的伪影。改动只有一个地方：把两个独立的 shrink 替换成联合收缩。

% 各向同性 TV 的 d 更新，替换掉两个独立 shrink smag = sqrt(vx.^2 + vy.^2); t = max(smag - 1/beta, 0) ./ max(smag, eps); dx = vx .* t; dy = vy .* t;

逻辑说明：smag是两个梯度分量的模长，阈值因子 t 由模长共同决定，而不是分别对vx和vy独立收缩。这样梯度方向保持不变，只有幅值被压缩，斜向边缘不会被切成水平加垂直两段。

参数说明：eps加在分母上防止零向量除零。阈值依然是1/beta，和独立 shrink 一致。切换成各向同性后，重建图像整体会更加平滑自然，但边缘会稍微柔和一点，如果追求极致锐利可以退回各向异性版本。

4.2 压缩感知重建：换保真项，保留 Bregman 骨架

Split Bregman 之所以通用，是因为它不关心保真项长什么样。去噪时保真项是||u-f||²，压缩感知时变成||A u - y||²，其中 A 是测量算子、y 是测量数据。Bregman 迭代的骨架不变，变的只是第一步求解 u 时的线性系统。

经典的应用是部分傅里叶采样重建，常见于磁共振成像类实验。此时 A 就是“取 k 空间部分坐标”，实现上极其简洁：每轮 Bregman 迭代结束后，把测量过的频域坐标覆盖回原值。

% 数据一致性投影：在每次 Bregman 迭代结束后执行 U = fft2(u); U(mask) = Y(mask); % 只覆盖被采样位置 u = real(ifft2(U));

逻辑说明：mask 是逻辑矩阵，true 表示该频点被采集；Y 是测量得到的 k 空间数据。这个投影操作把 u 的频域强制拉回测量值，未采样位置不受影响，由 TV 正则项自行填补。

参数说明：这种写法只适用于 A 是傅里叶下采样算子的情况。如果 A 是随机高斯矩阵这类一般线性算子，不能靠投影解决，必须把欧拉方程改成(A'A + βL'L) u = A'y + β(...)的形式再求逆，那就要引入共轭梯度法了。压缩感知场景下 lambda 通常调得很小甚至趋向 0，因为数据保真由投影步骤负责，正则项负责“编造”未采样位置的合理结构。

4.3 彩色图像处理：通道循环与联合 TV

直接把灰度算法套到 RGB 三个通道上可以跑，但结果往往出现伪彩，因为三个通道独立收缩会破坏通道间相关性。常见做法是转到 YCbCr 空间，只对亮度通道 Y 做 TV 去噪，色度通道 Cb、Cr 做轻度平滑即可。

% 彩色图像：只对亮度通道做 TV，色度通道轻度高斯滤波 yuv = rgb2ycbcr(f); y = split_bregman_denoise(yuv(:,:,1), lambda, beta, maxIter); cb = imgaussfilt(yuv(:,:,2), 1.5); cr = imgaussfilt(yuv(:,:,3), 1.5); u = ycbcr2rgb(cat(3, y, cb, cr));

逻辑说明：人眼对亮度细节最敏感，对色度细节不敏感，所以色度通道交给高斯滤波完全够用，还省下了大量计算。

参数说明：imgaussfilt的 sigma 取 1.5 是个中庸值，色度噪声重就加大到 2~3。如果 MATLAB 版本没有rgb2ycbcr，可以手动用标准转换矩阵实现，效果差别不大。

5. 避坑与常见问题：五个最容易翻车的细节

5.1 输入没有转 double：全白图与 NaN 事故

第一次跑这个源码时，直接把imread读进来的 uint8 矩阵传进主循环，输出一片白，命令行里到处是 NaN。这是因为 MATLAB 的 uint8 运算会截断：图像差值一旦为负直接被置成 0，后面 FFT 和除法全乱套。

解决方式是入口处强制统一数据格式。在 demo 最前面加一行f = im2double(f)，如果输入是三维图先转灰度。这是一条铁律：所有 Bregman 变量和中间结果都必须是 double，只有最后展示或保存时才允许转回 uint8。

5.2 meshgrid 维序搞反：横竖条纹分家的「鬼图」

症状是重建图像水平方向干净、垂直方向有条纹，或者反过来，看起来像把图像网格撕裂了。原因是写频域特征值时把meshgrid两个参数的顺序搞反了。

[Kx, Ky] = meshgrid(0:N-1, 0:M-1)中，Kx 的每一行是列频率，尺寸是 M×N；Ky 的每一列是行频率。如果把两个范围对调，水平差分的特征值会错位，导致频域分母不对应，重建自然就歪了。解决方法是写完后立刻自检：size(Kx)必须等于[M, N]。把它当成断言写在预计算后面，比靠眼睛看输出靠谱得多。

5.3 d 更新里符号写反：能量曲线震荡的元凶

现象是能量曲线前几十次迭代正常下降，之后突然反弹并持续震荡，图像在锐利和虚化之间反复横跳。原因多半是 d 更新里把+bx写成了-bx。Bregman 罚项 b 的符号必须和 b 的更新公式严格配对：b 按b + (D(u)-d)更新，d 更新时就要用D(u)+b；如果一边加号一边减号，整个反馈回路变成正反馈，误差被越放越大。

我从不硬记符号，每次拿到一份新的 Bregman 代码，先用 5×5 随机小矩阵做一次目标函数数值核对，确认符号没问题再上真实图像。这个习惯帮我避免了很多无效调参。

5.4 各向同性阈值写成逐分量：边缘发虚的坑

如果你把算法从各向异性改成各向同性，但 d 更新仍然用两个独立的shrink，你会得到一张边缘发虚、像被磨皮过的图。原因前面说过：独立收缩把斜向边缘的梯度拆成水平和垂直两段，各自被不同程度压缩，斜向边缘就糊了。

解决方式就是第 4.1 节的联合收缩代码，注意分母的max(smag, eps)不能省。如果做完联合收缩后图像还是虚，检查smag是否用的sqrt(vx.^2+vy.^2)，有人会图省事写成abs(vx)+abs(vy)，那又是另一种错误。

5.5 边界处理不匹配：周期边界与图像接缝

FFT 求解隐含周期边界假设，图像左右边界被视为相邻，上下边界同理。如果读图后做过padarray复制边缘，或者对图像做了分块处理，周期边界假设就被破坏了，重建图边缘出现肉眼可见的接缝。

解决方式有两种。如果图像不大，直接全图做 FFT 求解，不 padding，这是最省事的路径。如果图像太大必须分块，块与块之间要留重叠区域，并且重叠部分最终取均值融合。另一种思路是改用 DCT 边界：把 FFT 对角化换成dct2，分母换成 DCT 特征值，这种边界假设是“对称延拓”，更适合自然图像。不少新版实现直接内置了 DCT 选项，选型时优先考虑。

6. 验证收敛性：残差曲线与迭代终点的判断技巧

6.1 把能量曲线记录下来

调参时最怕黑匣子：不知道迭代到第几步才算出个能用的结果。我在主循环里加一段能量记录，每次迭代都算一次目标函数值，结束后画成对数坐标曲线。

E = zeros(maxIter, 1); for k = 1:maxIter % ... 主循环体 ... [du_x, du_y] = D(u); E(k) = 0.5 * lambda * norm(u(:) - f(:), 2)^2 + ... sum(abs(du_x(:))) + sum(abs(du_y(:))); end semilogy(max(E, eps), 'LineWidth', 1.2); xlabel('iteration'); ylabel('energy (log)');

逻辑说明：能量由保真项和正则项两部分组成。因为能量跨度可能很大，用semilogy画对数坐标更容易看出平台期。max(E, eps)防止数据中有零导致画图跳点。

参数说明：E 的绝对值大小没有标准，重要的是形态。

6.2 看图说话：什么时候停

我第一次做压缩感知实验时，傻乎乎固定 500 次迭代，每跑一次等三分钟，改一个参数又等三分钟，后来发现很多结果其实跑到 80 次已经不动了。从此拿到任何 Bregman 实现，第一件事就是画能量曲线，确认它走平之后再谈视觉效果。参数调完之后，再花十秒钟看一眼曲线形态，能拦住九成以上的无效调参。希望帮到你。

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