
简介这是一套面向SAR图像去噪的MATLAB实现专注处理Gamma乘性噪声适合有一定MATLAB基础的遥感图像处理学习者用于科研实验及课程设计场景。PPBPatch-based Progressive Blending方法借助非局部patch相似性在抑制噪声的同时保留边缘与纹理细节实现中先搜索相似patch集合再按Gamma分布模型估计噪声并渐进融合从而在平滑噪声时保护目标结构。整个流程涵盖图像预处理、patch提取、相似性度量、噪声估计与融合重建形成完整可复现的去噪链路便于读者从零复现算法或作为改进基线。压缩包共13个文件、681KB包含7个.m源码、Windows/Linux等多平台mex可执行文件、示例图像lena.png及README说明m文件便于阅读和二次开发mex文件可直接运行加速README可帮助快速上手。已有326人学习/下载通过运行示例脚本可对比去噪前后PSNR/SNR指标理解patch相似性度量与迭代融合过程方便迁移到其他乘性噪声去除任务尤其适合开展SAR图像定量实验的研究者。1. SAR图像去噪里的Gamma乘性噪声为什么不能套高斯NLM星载SAR成像得到的强度图噪声机制和普通光学照片完全不同。每个分辨单元内的后向散射由大量散射体相干叠加形成的相干斑噪声在强度图上近似服从Gamma乘性分布观测值等于真实反射率乘上一个均值为1的Gamma随机量亮区噪声方差随亮度平方增长。经典的非局部均值NLM和BM3D都建立在加性高斯噪声假设上拿它们的欧氏距离直接度量patch相似性会系统性误判亮区的小差异被当成结构保留暗区的大差异又被当成噪声去除边缘和点目标很难两全。PPBProbabilistic Patch-Based方法把非局部patch的相似性判断从高斯核距离换成Gamma分布下的似然比距离再用迭代估计逐步修正权重。它在数学上和Gamma乘性噪声的物理模型对齐是SAR去噪里少见的、能在理论层面说得清且MATLAB可复现的非局部算法。下面从统计模型推到MATLAB实现再把参数怎么调、验证怎么做讲透。2. PPB的patch相似性度量从Gamma似然比到权重公式2.1 SAR强度图的统计模型与等效视数SAR强度图里观测值 (z) 与真实后向散射系数 (x) 之间是纯乘性关系[ z x \cdot n ]其中 (n) 服从形状参数为 (L)、尺度参数为 (1) 的Gamma分布。(L) 在SAR领域叫等效视数ENL由多视处理决定单视强度图 (L1)四视大约 (L4)。给定真实值 (x)条件概率密度写为[ p(z|x) \frac{L^L}{\Gamma(L) x^L} z^{L-1} e^{-L z / x} ]注意这个分布有两个关键性质。首先是均值不变噪声乘性因子的均值为1所以 (E[z|x]x)这保证了相干斑不会改变区域均值。其次是方差随亮度变化(\mathrm{Var}(z|x)x^2/L)亮度翻倍噪声标准差也翻倍。后一条直接决定了用高斯距离做SAR去噪会失败高斯距离认为所有像素的噪声水平相同而实际图像里亮区和暗区的噪声尺度相差一个数量级。PPB的出发点就是构造一个能自动按局部亮度归一化的相似性度量。2.2 单像素Gamma似然比距离判断两个像素 (a)、(b) 是否来自同一个真实反射率最自然的办法是似然比检验。原假设 (H_0)两个像素共享真实值 (\theta)备择假设 (H_1)两个像素各自独立对应自己的真实值。在 (H_0) 下(\theta) 的最大似然估计恰好是两个观测值的平均[ \hat{\theta} \frac{ab}{2} ]把分子上的 (p(a|\hat{\theta})p(b|\hat{\theta})) 与分母上的 (p(a|a)p(b|b)) 相除Gamma函数、系数项和指数项全部约干净得到一个非常干净的结果[ \Lambda(a,b) \left(\frac{4ab}{(ab)^2}\right)^L ]取负对数就得到单像素距离[ d(a,b) -\log \Lambda L \cdot \log \left(\frac{(ab)^2}{4ab}\right) ]这个距离有两个值得记住的性质。第一(ab) 时距离为0像素差距越大距离越大单调性符合直觉。第二也是最关键的(a)、(b) 同时放大相同倍数时距离不变因为它依赖的是比值。比如 (a10,b12) 和 (a100,b120) 的距离完全一样这正好匹配Gamma乘性噪声的物理特性亮区允许更大的绝对波动相似性应该由相对比例决定而不是绝对差。在MATLAB里这个距离用一行就能实现function d gamma_pixel_dist(a, b, L) % 计算Gamma乘性噪声下两个像素的概率距离 % a, b: 强度值需要为正数; L: 等效视数 d L * log( (a b).^2 ./ (4 .* a .* b) eps ); end加eps是防止强度图里有零值导致除零。这个函数只用于说明距离形式实际patch级计算应该向量化见第3章的主函数。2.3 从单像素到非局部patch权重单像素距离对相干斑过于敏感PPB的思路和NLM一样用patch内的多个像素联合判断。假设patch内各像素条件独立两个patch之间的距离就是内部所有像素距离之和[ D(i,j) \sum_{k \in \text{patch}} L \cdot \log\left(\frac{(u_{i,k}u_{j,k})^2}{4 u_{i,k} u_{j,k}}\right) ]这里 (u) 不是原始含噪图而是迭代过程中的当前估计值。权重定义与NLM保持同一形式[ w(i,j) \exp\left(-\frac{D(i,j)}{h}\right) ](h) 是平滑参数控制相似度向权重的映射强度。与高斯NLM相比两者的框架完全一致唯一区别是距离函数从欧氏距离换成了Gamma似然比距离但这一个变化就把算法从“假设加性高斯噪声”搬到了“精确匹配乘性Gamma噪声”上。PPB还有一层迭代。第一次迭代时只能用含噪观测 (uz) 计算权重加权平均后得到初步去噪结果第二次迭代用这个初步结果重新计算patch距离权重会更可信。一般迭代三次左右权重质量提升明显继续迭代边际收益递减。3. MATLAB实现PPB滤波主函数与逐像素加权流程3.1 主函数代码下面给出一个完整的PPB滤波函数针对SAR强度图设计没有依赖任何工具箱图像处理大作业和科研原型都够用。function denoised ppb_gamma_sar(z, L, h, search_r, patch_r, n_iter) % PPB_GAMMA_SAR Gamma乘性噪声下的PPB非局部去噪 % denoised ppb_gamma_sar(z, L, h, search_r, patch_r, n_iter) % z : SAR强度图double类型像素值必须为正 % L : 等效视数可用estimate_enl估计 % h : 权重平滑参数典型值20~40 % search_r: 搜索窗半径典型值7~10 % patch_r : patch半宽典型值1~2实际patch尺寸为2*patch_r1 % n_iter : 迭代次数典型值2~4 z max(z, eps); % 消除零值避免除零 u z; % 初始估计用含噪观测 [rows, cols] size(z); p search_r patch_r; % 总填充宽度 zt padarray(z, [p p], symmetric); denoised zeros(rows, cols); for it 1:n_iter ut padarray(u, [p p], symmetric); for i 1:rows for j 1:cols ci i p; % 原图(i,j)在填充图中的坐标 cj j p; ref ut(ci-patch_r:cipatch_r, ... cj-patch_r:cjpatch_r); acc 0.0; % 加权累加器 wsum 0.0; for di -search_r:search_r for dj -search_r:search_r ti ci di; tj cj dj; nb ut(ti-patch_r:tipatch_r, ... tj-patch_r:tjpatch_r); if di 0 dj 0 w 1.0; % 自身patch权重固定为1 else ratio (ref nb).^2 ./ ... (4 .* ref .* nb eps); dist L * sum(log(ratio(:) eps)); w exp(-dist / h); end acc acc w * zt(ti, tj); wsum wsum w; end end denoised(i, j) acc / wsum; end end u denoised; % 更新当前估计用于下一轮权重 end end3.2 代码逻辑与调用方式整体流程分三层。外层是迭代次数n_iter每轮都用上一轮的去噪结果更新 (u)。中层是遍历所有像素对每个像素在搜索窗内找相似patch。内层是距离计算和权重累加这是整个算法的核心计算量所在。两个实现层面的选择需要说明。权重计算用的是当前估计 (u) 的patch因为 (u) 已经比含噪图更接近真实反射率距离更可靠但加权平均的对象是原始观测 (z)而不是 (u)。这样做的原因是 (z) 对真实反射率无偏直接用 (u) 做加权平均会把迭代过程中产生的偏差逐步固化。如果想把实现改成论文中严格的迭代加权最大似然形式只需把acc acc w * zt(ti, tj)改成acc acc w * ut(ti, tj)两种写法在主循环里只差一个变量名。距离计算中的ratio是一个与patch同尺寸的矩阵log(ratio(:))把所有像素的贡献相加。这里没有调用第2章的gamma_pixel_dist因为那个函数接标量、逐像素循环太慢直接矩阵运算可以利用MATLAB原生速度。eps除了防止除零也保证了零值像素不会让log产生-Inf。调用方式很简单denoised ppb_gamma_sar(z, 4, 25, 7, 1, 3);其中z是double类型的强度图。很多坑出在输入上uint8图像必须提前转double否则refnb做的是整数运算接近0的像素会直接截断成0或255强度值范围太大比如0到4095会导致距离数值整体偏大h需要同步放大。一般建议先把图像缩放到0到1区间再处理后面调h更顺。4. PPB的5个关键参数与MATLAB加速手段4.1 参数取值范围与互相影响PPB有5个参数需要设定它们不是各自独立的调参要按顺序来。参数含义典型值调小/调大的影响L等效视数1~16调大则距离变大h需同步放大h权重平滑参数10~40调大更光滑调小保留细节和噪声search_r搜索窗半径7~10调大候选patch更多计算量平方级上升patch_rpatch半宽1~2调大结构保持更强边缘容易糊n_iter迭代次数2~4超过4轮收益下降过平滑风险增加L是最先确定的参数。合成图像可以直接知道精确值真实SAR图需要估计。常用方法是找一块均匀区域水域、阴影区等用均值平方除以方差function Lest estimate_enl(img) % 估计等效视数输入应为均匀区域的强度值 % 切忌直接对整幅含纹理图像计算方差会被纹理抬高 p img(:); Lest mean(p)^2 / var(p); end真实图像上如果实在找不到均匀区域可以先做一次 (3\times3) 均值滤波再估计但这样得到的L偏大使用时要把h也相应调大。常见的错误是拿整幅图估计L纹理区的高方差会把L压到1以下算出来的权重失去鉴别力。h的数值和L直接相关。从距离公式可以看出距离本身乘以L所以L4时的距离范围是L1时的4倍h不放大就会让除自身外的所有权重都趋近于0。一个实用经验是令 (h \alpha \cdot L \cdot (2\cdot patch_r1)^2)(\alpha) 从0.5起步。如果结果太嘈杂就加大 (\alpha)太光滑就减小。相比手动试值这个公式至少让h在不同L和不同patch尺寸之间保持可比。patch_r和search_r的选择受图像内容影响。patch_r1对应 (3\times3) patch对边缘保持最好但对单个强散射点会敏感patch_r2对应 (5\times5)结构更稳定代价是边缘过渡处被拉平。search_r超过10以后搜索窗内绝大多数patch都来自不相关的区域这些低权重patch虽然对输出的贡献不大却白白消耗了好几倍的计算时间。4.2 迭代次数不是越多越好PPB的迭代机制和多数迭代算法不同第一轮用含噪图计算权重权重质量差但能去除大部分强噪声第二轮用清除了强噪声的图像重新计算距离边缘处的权重开始变得准确第三轮到第四轮patch距离趋于稳定输出基本收敛。再往后u中的真实结构也开始被平均掉点目标和精细纹理逐步消失出现类似过度平滑的分段常量效应。判断是否过平滑的方法很简单比较两次迭代输出的差值图像。迭代4次后如果差值图里出现了明显的结构轮廓说明正在磨掉真实目标。对于大多数强度图n_iter3是开始调参的基准值不要一开始就设10。4.3 用parfor和预提取改善速度前面主函数的双层循环在256x256图上跑一轮search_r7时要计算几百万次patch距离MATLAB里可能要几十秒到几分钟。第一个加速手段是把最外层像素循环改成parforparfor i 1:rows row_out zeros(1, cols); for j 1:cols % 内部逻辑与ppb_gamma_sar一致 % 只写入 row_out(j)不修改共享变量 end denoised(i, :) row_out; endparfor要求每次迭代只写自己的输出行所以要把累积结果先放到临时行向量row_out循环结束后一次性赋值。这样改动不会影响数值结果核心计算过程在多个worker上并行。第二个常见做法是预提取patch。用im2col把整幅图的patch向量化成矩阵再按列索引搜索窗内的候选patch距离计算变成一次矩阵运算。缺点是内存占用大search_r7、patch_r1时每个中心像素要保留约450个候选patch的向量内存不足时反而比循环更慢。适合中小尺寸图像大图建议保留循环结构只加parfor和外层迭代提前终止判断。5. 用合成SAR图验证PPB效果PSNR、ENL与h自适应验证PPB有没有写对最好先用合成Gamma乘性噪声图做定量测试因为真实SAR图没有干净的参考真值。生成合成图很简单x im2double(imread(cameraman.tif)); % 模拟真实反射率 L 4; % 等效视数 rng(0); z x .* randg(L, size(x)) / L; % 乘性Gamma噪声randg生成形状参数为L、尺度为1的标准Gamma随机数除以L后均值为1乘到x上正好满足Gamma乘性模型。注意不要用randn或poissrnd前者是加性高斯后者是泊松计数都不匹配SAR强度图的噪声特性。运行PPB并计算峰值信噪比den ppb_gamma_sar(z, L, 25, 7, 1, 3); psnr 10 * log10(1 / mean((den(:) - x(:)).^2));噪声极强时PSNR只有十几dB去噪后能提升4到8dB就说明权重计算和迭代逻辑基本正确。PSNR之外还要看边缘保持指数EPI它评估的是图像在去噪后保留了多少边缘强度epi sum(sum(abs(diff(den, 1, 2)))) / ... sum(sum(abs(diff(z, 1, 2))));EPI接近1说明边缘被完整保留低于0.3说明过度平滑严重。合理区间通常在0.4到0.7具体取决于噪声强度和h。EPI过高时噪声也没怎么去掉需要结合PSNR一起判断。真实SAR图上没有真值验证思路换成视觉加统计。视觉上重点看两处点目标和细线是否存活均匀区域是否变平。统计上用第4章的estimate_enl找一个均匀区域对比去噪前后的ENL。去噪前ENL是4左右去噪后应该明显提升如果某块区域的ENL提升到几百多半是过度平滑了。最后留一个实用的h自适应技巧。真实图像不知道最优h可以在第一次迭代时收集所有搜索窗内的patch距离把中位数作为距离分布的参考% 在ppb_gamma_sar内部第一次迭代时收集dist到数组D h 3 * median(D(:));这样h会自动跟随L、图像灰度和patch尺寸变化避免每次换图像都重新试参。h和距离中位数的比值在2到4之间通常是个安全起点想更平滑就加大倍数想保留细节就取小倍数。本文还有配套的精品资源点击获取