6.4 图像去噪的扩散方程
6.4 图像去噪的扩散方程
巡回的最后一站最近的日常：手机夜拍的照片。把灰度图当成一张温度场——像素值就是温度——去噪与传热在数学上变成了同一件事：热传导抹平温度尖峰，也抹平噪声；但代价同样是抹平"棱角"，而图像的棱角（边缘）恰恰是要保住的信息。本节从"热核即模糊"讲起，看各向异性扩散如何在抹平噪声与保住边缘之间做取舍，再把 5.3 的正则化请进来收拢全局，顺路交代分数阶扩散的一笔应用账（旧版 5.3 节的知识在此安家）。
把照片当成温度场
灰度图像是函数 u(x, y)：定义在像素网格上的"温度分布"。噪声是叠加在它上面的高频随机扰动——按 3.1 的谱语言，噪声的能量集中在高波数段。各向同性热扩散方程 u_t = a²·Δu 恰好按指数衰减高波数（衰减率 a²k²，波数越高衰减越快），演化若干步的效果就是高斯模糊——热核卷积（2.4 的实验三）与图像处理的高斯滤波是同一条数学。去噪的机理与代价在此看得一清二楚：高频噪声被衰减的同时，图像自身的高频成分——边缘、纹理——被同等衰减。各向同性扩散是一把不分敌我的钝刀。
钝刀的困境用一个极端场景说透：一张黑白文档照片，噪声颗粒与文字笔画同处高频段，任何"频率一刀切"的滤波注定顾此失彼。出路是让扩散知道边缘在哪：边缘两侧是不同的物体，扩散不该跨过边缘；边缘内部与平坦区，扩散照常进行。这就是"各向异性"的确切含义——扩散强度不再处处相同，而是逐点、逐方向地受局部结构调控。
佩罗纳–马利克：让边缘指挥扩散
1990 年 Perona 与 Malik 提出的修正在扩散系数上做文章：
u_t = div( c(‖∇u‖)·∇u )，c(s) = 1/(1 + (s/λ)²)。
扩散系数 c 不再是常数，而是梯度的递减函数：梯度小（平坦区、噪声）时 c 接近 1，照常扩散；梯度大（跨越边缘）时 c 趋于零，扩散自动关闭。参数 λ 是"门槛"：梯度超过它的结构被当作边缘保护，低于它的被当作噪声抹平。换个读法更妙——把 div(c·∇u) 展开为 c·Δu 加 ∇c·∇u，第二项在边缘两侧方向相反，效果是锐化：大梯度处方程从"抹平"切换成"增强"，边缘不但不糊反而更清晰。
这个设计有惊无险之处也该如实交代：c 递减的方程不适定——2.3 的老朋友——数值上会放大噪声的 staircasing 伪影；工程解法是把梯度换成先高斯磨光后的梯度（正则化后向前向后扩散），或在若干步后停止迭代。"扩散多少步"本身成了正则化参数：步数少去噪不足，步数多边缘被啃——与 5.3 反演的 λ 旋钮完全同构。
图解：各向同性与各向异性的分野
两种扩散去噪对比示意
数值实验：保边去噪的全流程验收
在一幅合成图上走全流程：生成带边缘的干净图、加高斯噪声、跑各向同性与各向异性两版扩散、用两个指标验收——峰值信噪比（整体保真）与边缘保持度（边缘位置的梯度还剩多少）。
import numpy as np rng = np.random.default_rng(1) n = 64 clean = np.where((np.indices((n, n))[1] > n//2), 1.0, 0.0) # 半黑半白竖边 noisy = clean + 0.3*rng.standard_normal((n, n)) def lap(u): return (np.roll(u,1,0) + np.roll(u,-1,0) + np.roll(u,1,1) + np.roll(u,-1,1) - 4*u) u = noisy.copy() dt, steps, lam = 0.1, 40, 0.15 for _ in range(steps): gx = np.roll(u,-1,1) - u; gy = np.roll(u,-1,0) - u # 前向差分梯度 c = 1.0/(1.0 + (np.abs(gx) + np.abs(gy))/lam)**2 # 扩散系数随梯度递减 u = u + dt*(c*lap(u) + 0.2*(np.roll(c,1,1)-c)*gx) # 简化各项异性步 def psnr(a, b): return 10*np.log10(1.0/np.mean((a-b)**2)) blur = noisy + 0.1*lap(noisy) # 各向同性对照（一步热核） edge_true = np.abs(np.diff(clean[:, n//2-3:n//2+3], axis=1)).max() edge_pm = np.abs(np.diff(u[:, n//2-3:n//2+3], axis=1)).max() edge_blur = np.abs(np.diff(blur[:, n//2-3:n//2+3], axis=1)).max() print("原始含噪 PSNR:", round(psnr(noisy, clean), 2)) print("高斯型模糊 PSNR:", round(psnr(blur, clean), 2), " 边缘强度保留:", round(edge_blur/edge_true, 2)) print("各项异性 PSNR:", round(psnr(u, clean), 2), " 边缘强度保留:", round(edge_pm/edge_true, 2))
典型读数：各向同性对照的 PSNR 略升但边缘强度只剩六七成——噪声小降、边缘大伤；各向异性扩散的 PSNR 明显更高且边缘保留九成上下——两个指标同时占优。调节 λ 再跑一遍能看到经典的取舍曲线：λ 过小什么都当边缘（噪声留下），过大什么都当噪声（边缘啃掉）。把 steps 从 40 加到 200，边缘保留度缓慢下滑——迭代步数这把"正则化旋钮"的代价肉眼可见。
从扩散到正则化的收拢
把本章与 5.3 的线索拧成一股：图像去噪的三大门派可以写成同一个目标泛函家族——min ‖u − d‖² + λ·J(u)，惩罚项 J 取二范数平方（对应热扩散，即 Tikhonov）、取梯度的单调加权范数（对应各向异性扩散）、取全变分 ∫‖∇u‖（鲁道夫–奥斯古德等人的模型，保边最硬、枯造假平坦块）。扩散迭代是这一族泛函的梯度下降流——"扩散去噪"与"正则化反演"是同一枚硬币的动力学与变分两面。这一视角还顺手接住了旧版教材里分数阶偏微分方程那节的遗产：把二阶算子换成分数阶拉普拉斯，谱衰减率从 k² 变成 k 的低次幂，高频（纹理）衰减放缓——对重尾噪声与纹理密集的图像（皮肤、织物），分数阶扩散在保纹理上优于高斯热核；它的数值实现难点在非局部性（每个点受全域影响，矩阵稠密），5.2 与第四章的蒙特卡洛与预条件技术在那里再次会师。
💡 关键直觉：去噪算法的优劣不在"多平滑"，而在"对哪类结构平滑"——把先验写成惩罚范数或扩散系数，是把"我认为图像长什么样"数学化的两条等价通道。
本节要点
- 图像即温度场：噪声与边缘同居高频，各向同性热扩散模糊掉一切高波数；
- 佩罗纳–马利克：扩散系数随梯度递减，边缘处自动关闸，机制里还藏着轻度锐化；
- 不适定与旋钮：递减系数方程天然病态，门槛 λ 与迭代步数都是正则化参数；
- 双指标验收：PSNR 管整体保真、边缘保持度管关键结构，缺一即盲；
- 泛函统一：热扩散、各向异性、全变分同属一个惩罚家族，分数阶扩散是它的非局部延伸。
四张工单做完，方程的迁移能力已经验讫。最后一章把散落全册的手艺收拢成制度：用什么软件、走什么流程、上什么算力、凭什么签字放行。