简介：基于偏微分方程（PDE）的图像去噪MATLAB实现资源，系统覆盖各向异性扩散、方向扩散、高阶扩散、全变分去噪等多种经典方法，面向高校本科及硕士阶段的图像处理、信号处理等教研学习场景。压缩包共16个文件，其中包含9个.m源程序，涵盖全变分去噪、方向扩散、高阶扩散等多种核心算法，并配有边缘停止函数、自动参数计算、信噪比评估等辅助脚本，便于完整复现实验；另有2篇PDF文档，分别讲解基于偏微分方程的图像去噪原理和混合噪声去噪研究方法，可帮助读者从理论到实践建立系统认知；同时附有3张结果示意图和1个动态效果图，直观展示不同算法对噪声图像的平滑与边缘保持能力。整个资源包仅3.27MB，但代码、文档、图示配置齐全，结构简洁。目前已有187人学习下载。通过研读代码与文档，读者可掌握偏微分方程在图像去噪中的建模和离散化求解技巧，并能直接运行main.m观察去噪效果，通过修改参数对比不同扩散方式的特点，为深入研究图像修复、边缘检测等方向奠定基础。
1. 图像去噪遇到瓶颈时，PDE 是那条最值得走的路
做图像去噪做过一段时间的人都会有这种感觉：高斯滤波、中值滤波、小波阈值，调参调到吐，边缘还是被抹平。偏微分方程（PDE）方法给我的第一印象是数学味道太重，但真正把 Perona-Malik 模型跑通、看到带噪平面图上的文字边缘被保住、噪声被平滑掉之后，我意识到这其实是比一堆滤波核更值得放进工具箱的基线方案。这篇笔记不只讲概念，我会把离散格式怎么推、Matlab 代码怎么写、参数怎么调、边界条件和步长有哪些坑，一条条讲透。适合手里有 Matlab、想把去噪效果从“能用”往上拉一档的人，也适合要对比算法效果、需要可复现基线的工程场景。
2. 从热传导方程到各向异性扩散：PDE 去噪的基本原理与两种模型选型
2.1 热传导方程为什么让图像变糊：各向同性扩散的退化逻辑
PDE 去噪的起点是一维热传导方程在二维图像上的推广：把图像像素值当作温度场，噪声是高温点或低温点的突变。热传导的物理规律是温度差越大，热量流动越快，反映到数学上就是 ∂u/∂t = div(∇u)。对图像来说，这个扩散过程会快速拉平像素间的差值，噪声被抹掉，但同时也会把边缘的梯度一起抹掉。原因在于这个算子的扩散系数是常数，对每个方向一视同仁，所以叫各向同性扩散。
用 Matlab 实际跑一遍最常见的热传导迭代，你会发现第 5 次迭代时噪声已经明显减轻，但第 20 次迭代后整张图像蒙了一层雾，边缘过渡区变宽，细节损失严重。这就是线性扩散的退化逻辑：它不知道哪些梯度是噪声造成的，哪些梯度是真实边缘。噪声和边缘在局部都是高频变化，唯一的区别是尺度。所以后来做 PDE 去噪的人都在解决同一个问题：让扩散在边缘方向变慢，在边缘垂直方向变快，或者直接判断梯度大小来决定扩散强度。
从实现角度看，各向同性扩散只需要一个拉普拉斯算子，代价极低，但效果连经典的导向滤波都不如。真正让 PDE 方法翻身的思路是改变扩散方程里面的系数项，让方程本身学会“看梯度”。这就是下一节要展开的 Perona-Malik 模型。
2.2 Perona-Malik 模型与扩散系数的三种常见表达式
Perona-Malik 模型的核心是把热传导方程中的常数系数改成图像梯度模的单调减函数，扩散方程为 ∂u/∂t = div(g(|∇u|)∇u)。当梯度模很大时，g 趋向于 0，扩散基本停止，边缘保留；当梯度模很小时，g 接近 1，噪声被平滑掉。
扩散系数 g(|∇u|) 常用的有两种经典写法：指数型 g = exp(-(|∇u|/K)^2)，以及反比例型 g = 1/(1 + (|∇u|/K)^2)。指数型对梯度的响应更敏感，小梯度时扩散能力强，大梯度时几乎完全锁死；反比例型则更平稳，不容易在处理强边缘时产生阶梯效应。两者中的 K 是控制边缘保留强度的阈值，K 越大，被判定为“边缘”的梯度阈值越高，平滑力度越大；K 太小，噪声点也会被视为边缘，导致去噪不彻底。
第三种表达式在工程里也常见，是 g = 0.5 * (1 - (|∇U|/K)^2) 当梯度小于 K 时取反比例补集。它本质上是一个软阈值设计，但实际效果与反比例型差别不大，我建议新手先从反比例型入手，因为它对 K 的变化不那么敏感，不容易跑出一个满是黑斑的结果。
选哪个表达式不是数学偏好问题，而是噪声类型问题。对于高斯白噪声，反比例型配合 20 到 40 次迭代表现稳定；对于椒盐噪声，指数型的锁定效果更利于保护极细纹理，但如果 K 没调好，椒盐点周围的孤立像素反而会被当成边缘保留下来。
2.3 模型选型：P-M、曲率驱动与总变分，怎么按噪声和场景选
P-M 模型是入门首选，但要意识到它的边界：扩散系数只看像素点的梯度模，不看梯度方向，所以当边缘出现断裂时，断裂处也会被平滑。第二类常见 PDE 模型是曲率驱动扩散，比如曲率方程，让扩散沿等照度线方向发生，边缘断裂被修复的能力更强，但迭代步长不得不设置得很小，否则数值不稳定，Matlab 上跑一张 512×512 图往往要等几十秒。
第三类模型是总变分去噪，其本质是求解一个能量最小化问题，min ∫(|∇u| + λ/2(u - I)²)dxdy，迭代时会用到梯度下降或对偶方法。它的优点是去噪效果和边缘保持都相对均衡，但参数 λ 难调，且迭代过程比 P-M 更容易产生贴片状伪影。
所以我的选型经验是：如果是工程 demo、对比实验、批量图像处理，优先 P-M；如果图像里有大量长度较长但细的线条，比如血管、纹理扫描图，用曲率驱动更稳；如果不仅要视觉效果，还要让下游检测算法精度更高，总变分往往是最稳妥的候选。接下来先给出一个最小可跑的 P-M 扩散函数，这条代码我反复用了很多年。
function u = pm_diffuse(u0, K, dt, n_iter, lambda) % u0 输入灰度图像 double, 范围 0-255 % K 边缘阈值，越大平滑越强 % dt 时间步长，需要满足 dt <= 0.25 % n_iter 迭代次数 % lambda 梯度求导时是否使用中心差分，1 使用，0 使用前向差分 u = u0; for t = 1:n_iter % 四邻域梯度 gradN = u([2:end end], :) - u; gradS = u([1 1:end-1], :) - u; gradE = u(:, [2:end end]) - u; gradW = u(:, [1 1:end-1]) - u; if lambda % 中心差分梯度估计 gradMagN = sqrt(gradN.^2 + 1e-10); gradMagS = sqrt(gradS.^2 + 1e-10); gradMagE = sqrt(gradE.^2 + 1e-10); gradMagW = sqrt(gradW.^2 + 1e-10); else gradMagN = abs(gradN); gradMagS = abs(gradS); gradMagE = abs(gradE); gradMagW = abs(gradW); end % 扩散系数 cN = 1 ./ (1 + (gradMagN / K).^2); cS = 1 ./ (1 + (gradMagS / K).^2); cE = 1 ./ (1 + (gradMagE / K).^2); cW = 1 ./ (1 + (gradMagW / K).^2); % 显式更新 u = u + dt * (cN.*gradN + cS.*gradS + cE.*gradE + cW.*gradW); end end代码里需要注意几个点：边界像素通过复制延拓来处理，u([2:end end], :)表示把第二行复制到底部作为底部边界，这样不至于引入不必要的零梯度；中心差分对梯度模的估计更准确，但会放大噪声影响，所以代码留了lambda开关做对比；+1e-10是为了防止梯度模为 0 时除零。显式更新的时间步长限制是dt <= 1/4，原因是二维显式格式稳定条件为dt <= 1/4，我一般取 0.2，既不超限，又不会因为太小让收敛慢到看不见。
3. 用 Matlab 实现可复用的 P-M 去噪函数：离散格式到参数选择
3.1 显式有限差分格式推导：从连续方程到四邻域迭代
上一章的代码能跑，但如果你只拷贝不推导，后面调参会像开盲盒。这里把离散化过程说清楚。连续方程 ∂u/∂t = div(g(|∇u|)∇u) 在二维平面展开后，对每个像素 (i,j)，得到关于 x 和 y 方向的扩散项。常用的做法是把它拆成四个方向的通量：北、南、东、西，每个方向的扩散通量 = 扩散系数 × 该方向的梯度。
离散化时，梯度用最近邻差分：北方向梯度取当前像素减去上方像素，之所以取负号，是为了和图像坐标系的行方向对齐。然后每个方向的扩散系数由该方向梯度的模决定，这就是各向异性的来源。最终更新公式是 u_new = u_old + dt * (cN∇N + cS∇S + cE∇E + cW∇W)。注意这个公式里没有除以像素间距，因为在图像处理中像素间距固定为 1，所以离散格式天然把间距吸收进系数里。
稳定条件来自线性化分析：如果把扩散系数近似为常数 g，则迭代格式变成一个显式热传导方程更新，傅里叶分析表明步长必须小于等于 1/4，否则高频误差会指数放大，表现出来就是图像上出现棋盘格噪声且越迭代越重。这是 PDE 去噪里最容易踩的硬性约束，比任何超参数都重要。
中心差分与单向差分的差别也需要理解：中心差分对梯度模的估计考虑了相邻两侧的像素，边缘定位更准，但对孤立噪点的响应会比前向差分大。在含椒盐噪声的图上，如果 K 又取得小，椒盐点会被识别成边缘，迭代后保留成大黑点或大白点。所以看噪声类型选差分方式，比盲目追求“高级”格式更实际。
3.2 可直接运行的完整 Matlab 代码：读图、去噪、做对比
下面是一个完整脚本，包含加噪、去噪、输出 PSNR 三个环节。我不想只给一个函数，因为很多人拿过去跑不知道输入范围是什么，结果图全黑或全白。
% 读取灰度图并归一化到 0-255 范围 img = imread('demo.png'); if size(img, 3) == 3 img = rgb2gray(img); end u0 = double(img); % 添加高斯噪声 rng(0); uNoisy = u0 + 25 * randn(size(u0)); % 调用扩散函数 K = 20; dt = 0.2; n_iter = 30; uDenoised = pm_diffuse(uNoisy, K, dt, n_iter, 1); % 计算 PSNR psnr_noisy = psnr(uint8(uNoisy), uint8(u0)); psnr_denoised = psnr(uint8(uDenoised), uint8(u0)); fprintf('噪声图 PSNR: %.2f dB\n', psnr_noisy); fprintf('去噪图 PSNR: %.2f dB\n', psnr_denoised); % 显示 subplot(1,3,1); imshow(uint8(u0)); title('原始图'); subplot(1,3,2); imshow(uint8(uNoisy)); title('加噪图'); subplot(1,3,3); imshow(uint8(uDenoised)); title('P-M去噪');运行时请记住一个前提：pm_diffuse函数里的输入默认范围是 0 到 255。很多写过传统滤波器的人习惯把图归一化到 0 到 1，直接套进来会发现 K 和梯度的量纲完全对不上。梯度模的范围会缩小到 1/255，原来的 K=20 变得过于巨大，结果是整张图被严重平滑，边缘也没了。这是最常见的翻车现场。
关于 PSNR 的对比，我不能只告诉你“数值变高就是好”。PDE 去噪经常出现数值提升不明显、但主观边缘更干净的情况，因为 PSNR 对整体像素误差敏感，对结构边缘损失不敏感。所以脚本里我习惯同时输出去噪前后的差图，看差图中是否出现明显的边缘高亮，如果边缘被高亮，说明边缘细节被过度平滑了。
3.3 核心参数对照：迭代次数、步长 K 的选择策略
参数表比代码本身更重要，因为 P-M 模型的所有玄学都集中在三个超参数上。迭代次数决定扩散进行的长度，太少了噪声残留，太多了细节消失。K 控制边缘判断阈值，K 太小扩散系数几乎处处为 1，退化成高斯滤波；K 太大扩散系数处处接近 0，去噪无效。dt 则受稳定性限制，不能作为自由参数看待。
| 参数 | 常见范围 | 对结果的影响 | 调参建议 |
|---|---|---|---|
| n_iter | 10~80 | 迭代越多平滑越强，边缘逐渐变糊 | 先固定 30 次再微调 |
| K | 5~50 | K 越小边缘越尖锐，噪声越难去除 | 从噪声标准差的两倍开始试 |
| dt | 0.1~0.25 | 超过 0.25 必然发散 | 不要超过 0.24，留出余量 |
| lambda | 0 或 1 | 中心差分保边更好，前向差分更稳 | 噪声强时用 0，细节多时用 1 |
迭代次数和 K 存在明显的互相补偿关系：K 增大 30% 而迭代次数减少一半，可能得到相近的整体平滑度，但边缘保留趋势不同。我一般会先调 K，把边缘调到一个可接受状态，再用迭代次数去平衡噪声剩余量。这样做的好处是避免“越迭代越干净但边缘越糊”的恶性循环。
K 的一个经验参考值是噪声标准差的两倍。以 0-255 灰度范围为基准，高斯噪声 sigma=25 时 K=20 起步是合理的，但如果噪声是 sigma=10，K=20 会把许多真实细节当平滑对象。反过来，sigma=40 时 K=20 会导致椒盐状残留。另一个更稳的做法是把扩散系数的分母改成梯度的中位数，让 K 随图像内容自适应，这一步在后面的进阶技巧里会细讲。
4. 边界条件与步长的四个隐藏坑：排查 PDE 去噪的翻车时刻
4.1 去噪图像边缘发黑或发亮：边界延拓选择不当
现象：迭代十几步后，图像四条边出现一条明显的暗带或亮带，越靠边越严重。原因：边界像素的邻居不存在，很多实现里直接补零，导致边界梯度被严重高估，扩散系数变得异常，扩散行为与内部不一致。
解决：使用复制延拓，也就是把最外一圈像素向外复制。Matlab 里可以用u([1 1:end-1], :)与u([2:end end], :)配合实现。另一个常见选项是对称延拓，效果相似，但复制延拓更节能且稳定。如果你在调参时发现不同 K 下边界伪影差异很大，先检查边界处理代码。边界条件对内部去噪结果影响并不大，但对 PSNR 数值却有几 dB 的区别，评估时一定要避开边界统计。
4.2 迭代过程中出现棋盘格状噪声：步长 dt 超限或梯度过大
现象：前 10 次迭代还很正常，第 15 次开始图像出现周期为 1 像素的交替高亮点，之后全局发散。原因：显式格式的稳定性条件没满足。虽然 dt 设了 0.2，但图像灰度范围是 0-255，扩散系数最大接近 1，实际稳定性分析应基于dt * max(c) <= 0.25，如果梯度归一化没做好，等效步长会被放大。
解决：固定dt = 0.2只是一个保险值，更稳妥的做法是每个迭代步都动态计算最大扩散系数：max_c = max(max(cN(:)), max(cS(:)))，然后让dt = 0.25 / max_c - 0.01。这样即使某一步扩散系数异常增大，也不会直接爆炸。如果你用的是 0-1 归一化输入，dt=0.2 通常安全；但一旦图里有超高对比度的直线边缘，梯度模可能瞬间超过 K 数十倍，扩散系数虽然被压下去了，可梯度本身参与更新时依然可能局部超限。
4.3 迭代越多边缘越重的“描边感”：K 设置偏向低阈值
现象：去噪图边缘不但没有变柔和，反而出现像水墨画一样的黑色描边，边缘两侧各有一条晕带。原因：K 太小，边缘像素的扩散系数接近 0，而边缘紧邻的过渡带像素梯度中等，扩散系数偏大，导致过渡带被进一步平滑，视觉上边缘被“强调”了。
解决：把 K 提高到噪声标准差的三倍以上。另一个配套操作是给扩散系数的分母里加入一个小的正则项eps，让梯度模为零处的扩散系数不是 1，而是 1/(1+eps)，避免数值抖动。这件事在数学上叫正则化，在工程上就是防呆。如果提高 K 后噪声残留增多，不要急着继续加迭代，而是返回去检查差分格式是不是选了中心差分，中心差分在高噪声下反而更容易把残留噪声边缘化。
4.4 去噪后图像整体偏亮或偏暗：数据类型转换问题
现象：算法结果是 double 类型，显示前没有做数据类型转换，或者只做了uint8()截断，导致负值被截成 0，图像整体偏黑，亮度直方图变形。原因：扩散更新过程中像素值可能超过 0-255 范围，或者低于 0，uint8()只截断不缩放，丢失了动态范围。
解决：在显示前使用mat2gray做归一化，或者手动u = min(max(u, 0), 255)。对于需要保留物理意义的数据，比如医学影像灰度值，不能直接归一化到 0 到 255，而是要保存一个全局最小值最大值，去噪后按同一范围映射。这个问题在量化对比时最容易暴露，因为误用uint8()截断会凭空增加 3 到 5 dB 的 PSNR“假提升”。
5. 评估与验证：用模拟噪声图量化 PDE 去噪的效果
5.1 PSNR 与 SSIM 的完整评估脚本：加四种噪声做对照
只拿两张图看效果很难说服自己。我建议建立一个小型评估流程：合成四组带噪图，分别是高斯噪声 sigma=10、sigma=30、椒盐噪声密度 0.05、泊松噪声。每组用一个基准参数跑一遍，把去噪图的 PSNR 和 SSIM 都打出来，再和原图算一遍梯度差，作为边缘保留的代理指标。
noise_types = {'gaussian', 'gaussian', 'salt & pepper', 'poisson'}; params = struct('sigma', [10, 30], 'density', [0.05], 'lambda', [1]); K = 20; dt = 0.2; n_iter = 30; for i = 1:4 if i == 1 nimg = imnoise(uint8(u0), 'gaussian', 0, (10/255)^2); elseif i == 2 nimg = imnoise(uint8(u0), 'gaussian', 0, (30/255)^2); elseif i == 3 nimg = imnoise(uint8(u0), 'salt & pepper', 0.05); else nimg = imnoise(uint8(u0), 'poisson'); end nd = double(nimg); result = pm_diffuse(nd, K, dt, n_iter, 1); psnr_after = psnr(uint8(result), uint8(u0)); ssim_after = ssim(uint8(result), uint8(u0)); grad_before = sum(abs(gradient(double(nimg))), 'all'); grad_after = sum(abs(gradient(result)), 'all'); fprintf('噪声类型 %d: PSNR=%.2f, SSIM=%.4f, 梯度保留率=%.2f%%\n', i, psnr_after, ssim_after, 100*grad_after/grad_before); end这里需要解释几个设计决定。梯度保留率不是越高越好，因为噪声本身也有梯度；但如果你把去噪图和干净图的梯度做比较，会发现边缘区域梯度若掉到 60% 以下，主观上就明显发虚。作为一个辅助指标，我一般看的是边缘区域梯度保留率，而不是全图统计。全图梯度保留率会被大面积平坦区域的噪声掩盖，变得没有区分度。
SSIM 在这个场景下比 PSNR 更可信，因为它的计算考虑局部结构和亮度对比。PDE 去噪一个常见问题是整体噪声被平滑后 PSNR 提升很多，但 SSIM 上不去，这说明结构信息被破坏了。如果出现这种结果，优先调小 K，增加迭代次数，而不是反过来。
5.2 参数网格搜索：K 与迭代次数如何协同寻优
参数搜索不要手调，写一个两层循环，把 K 从 5 到 40、n_iter 从 10 到 80 各取若干值，跑一个 5×5 的网格。输出最优 PSNR、最优 SSIM 以及两者的妥协点。代码量很小，但得到的结果比手动试错可靠得多。
K_list = [5, 10, 15, 20, 30, 40]; iter_list = [10, 20, 30, 50, 80]; best_psnr = 0; best_ssim = 0; for K = K_list for n_iter = iter_list result = pm_diffuse(double(uNoisy), K, 0.2, n_iter, 1); p = psnr(uint8(result), uint8(u0)); s = ssim(uint8(result), uint8(u0)); if p > best_psnr best_psnr = p; best_K_psnr = K; best_iter_psnr = n_iter; end if s > best_ssim best_ssim = s; best_K_ssim = K; best_iter_ssim = n_iter; end end end网格搜索也有坑：如果噪声是椒盐噪声，PSNR 最优的 K 往往偏大，结果图边缘模糊但数值好看。所以最后选参时建议以 SSIM 为主，PSNR 为辅，并且一定要看一眼差图。很多情况下最优参数在数值上的领先非常小，可能只有 0.3 dB，但主观差的十万八千里，这时候就不要再继续加参数密度了，先检查边界处理和数据范围。
5.3 一个有效的目测验收方法：细节放大加残差图
数值指标只能过滤出绝对不合格的参数，不能告诉你真实视觉感受。我的习惯是取一块包含边缘、平坦区域、纹理的局部区域，把原始清晰图、去噪图、噪声图并排放大到 200%，观察边缘周围有没有晕染，平坦区域有没有残留颗粒。
另外一个被低估的工具是残差图，即噪声图减去去噪图。残差图里应该只包含噪声，如果残差图里能看到清晰的结构轮廓，说明算法把真实信号也抹掉了。这个判断比任何曲线都直接。残差图的显示需要做归一化，否则像素值范围很小，屏幕上看起来一片灰，什么都看不出来。
6. 从 demo 到可用：把 P-K 阈值改成自适应的进阶做法
入门代码里 K 是全局常数，这在实际图像上限制非常大。背景平坦区域梯度和边缘附近纹理区域的梯度动态范围可能相差两个数量级，一个全局 K 注定顾此失彼。我后来做的第一个有效改进是把 K 变成逐像素的：每个像素的 K(x,y) 取该像素邻域梯度模的中位数乘一个缩放因子。这样平坦区域的中位数梯度很小，局部 K 小，扩散充分；纹理区域中位数梯度大，局部 K 大，扩散减弱。
Matlab 里可以用nlfilter做，但速度很慢；我习惯先算梯度模图，再用imfilter配合一个 5×5 中值滤波器模拟局部中位数，代码简化为Kmap = img * medfilt2(gradMag, [5 5])，然后把这个 Kmap 逐像素代入扩散系数。效果是去噪图在纹理丰富区域的细节明显比全局 K 好，而且迭代次数敏感性降低。
第二个值得做的改进是加入边缘停止机制的松弛项：让扩散系数在连续迭代中缓慢恢复，即 gk_new = α*gk_old + (1-α)*g_current，α 取 0.8 到 0.9。这样做的好处是避免迭代后期扩散系数归零，算法空转。很多 P-M 实现跑 30 次迭代之后基本不再变化，就是因为扩散系数全被锁死，加了这个松弛项后，尾段迭代还能继续平滑残留下的小尺度噪声。
第三个技巧是用两步策略做复杂噪声：先用一个相对小的 K 跑 10 次迭代去掉明显离群点，再用一个较大的 K 跑 20 次迭代做整体平滑。这个思路来源于对噪声尺度的分离假设，对混合噪声效果比单次调参更稳定。我踩过最深的坑是把所有希望都寄托在 K 上，K 调小边缘好但噪声残留重，K 调大噪声净但图发虚，后来才意识到问题的本质是噪声尺度不单一。
现在我的 PDE 去噪脚本里默认不再用固定 K，而是先算梯度中位数，再用两步策略跑迭代。这个改动让原来那些容易翻车的图基本都稳住了。如果你也遇到参数老是换图就不灵的情况，建议优先把 K 改成自适应，再回头看别的问题。希望这篇笔记能帮你少走一段我走过的弯路。
本文还有配套的精品资源，点击获取
