简介面向遥感图像处理与合成孔径雷达变化检测场景的基于MATLAB的实现包覆盖四种最常见的差异图生成方法分别是比值法、差值法、均值比算法和对数比法。包内共10个文件含4个MATLAB脚本对应上述算法和6张位图示例图像包含原始雷达图像及各方法生成的差异图便于直观对比整套资源仅478KB轻量精简适合变化检测入门学习者、算法研究人员及需要快速验证差异图方法的工程人员使用。目前已有791人学习下载关注度较高。通过分析脚本和结果图可深入理解各差异图对地表回波强度变化的响应特点比值法简单直观差值法突出快速变化区域均值比法削弱雷达系统参数误差对数比法则能有效压缩亮暗目标对比度、提升检测准确性。读者可据此结合实际场景选择合适算法也可借助源码进行二次开发或算法融合。1. 先搞清楚SAR差异图为什么差值在SAR图像上会翻车两期同一区域的SAR图像想快速找出哪些地方变了——新建了建筑、水面涨了、塌了一角——最朴素的做法就是把两期图逐像素相减。但这个操作在SAR图像上几乎必翻车SAR回波自带乘性斑点噪声直接相减以后满屏椒盐点真实变化区域完全淹没在噪声里。所以SAR变化检测的第一步不是「怎么减」而是「用什么算子生成差异图」。差异图把两个时相浓缩成一张能体现差异的单通道图供阈值分割、聚类、深度学习分类后续使用。这篇笔记围绕差值图、比值图、对数比值图、均值比值图这四种最常见方案讲清公式、MATLAB实现、参数选型和踩坑记录。适合做变化检测的从业者、刚接触SAR图像处理的研究生以及准备把差异图作为深度学习辅助通道的工程师。2. 四种最常见SAR差异图的MATLAB实现公式、代码与选型对比2.1 SAR数据模型为什么算子形态由噪声性质决定SAR强度图像服从近似Gamma分布的乘性噪声模型I(i, j) σ(i, j) × n(i, j)其中σ是地物真实后向散射系数n是服从某种分布的乘性斑点噪声。对同一位置的两个时相σ可能没变但n完全不同。直接用I1减I2差值里既包含「真实变化」又包含「两期噪声之差」而乘性噪声之差的方差和信号本身量级接近信噪比极差。比值算子能把乘性项约掉I1 / I2 σ1n1 / σ2n2 ≈ σ1 / σ2 × (n1/n2)。真实后向散射的比值保留下来噪声项被压缩。这就是SAR差异图领域「能用比值就不用差值」的根本原因。后面几种差异图全是围绕「如何把乘性噪声的影响降到最低」展开的。2.2 差值图公式最简单在SAR上最不稳定差值图的公式D(i, j) I1(i, j) − I2(i, j)它理解成本最低光学图像变化检测里很常用但SAR图像上直接相减会把斑点噪声和真实变化混在一起。工程上我几乎不单独用它只拿它做对照。MATLAB实现非常简单% 读取两个时相的强度图im2double 转到 [0,1] 浮点域 I1 im2double(imread(sar_t1.png)); I2 im2double(imread(sar_t2.png)); D I1 - I2; % 差值图正值为后向散射增强负值为减弱逻辑说明必须先转double再相减。uint8图像直接相减时负数会被截断成0负变化区域直接丢失。差值图保留了变化方向增强和减弱用正负区分但代价是噪声也被完整保留。如果数据已经是单视SAR强度图而不是多视平均后的图差值图基本没法直接看。2.3 比值图与对数比值图SAR差异图的标准答案比值图公式R(i, j) I1(i, j) / I2(i, j)比值图对乘性斑点噪声天然不敏感但在分母接近0时会出现极大值导致动态范围被个别像素拉爆。所以实际使用中更常见的是对数比值图LR(i, j) |log(I1(i, j) / I2(i, j))|取对数有两个作用把乘性关系变成加性关系让直方图更接近对称分布取绝对值后增强和减弱统一映射到正值后续阈值分割只需要处理单侧分布。% 比值图 R I1 ./ (I2 eps); % 对数比值图eps 防止 log(0) 出现 -inf LR abs(log((I1 eps) ./ (I2 eps)));逻辑说明eps是MATLAB内置的极小浮点数约2.2e-16加在分母上只为了保证数值稳定不影响真实像素值。实际处理SAR图像时如果在阴影区、静水面存在大量接近0的像素建议用一个更大的常数比如0.01替代eps避免比值图出现极端尖峰。参数说明对数底默认用自然对数也可以用log2或log10只是数值范围不同对后续归一化没有影响。取了绝对值以后LR的取值范围是[0, ∞)理想情况下未变化区域集中在0附近变化区域形成一个拖尾或第二个峰。2.4 均值比值图用局部窗口把斑点噪声先压一遍均值比值图是比值图的改进版。核心思想先对两个时相分别做均值滤波再用滤波结果做比值MR(i, j) mean(I1, w) / mean(I2, w)窗口内均值对斑点噪声有抑制作用因此均值比值图在低信噪比场景下比像素级比值图更稳代价是空间分辨率下降细小变化容易被窗口抹掉。w 5; % 滤波窗口尺寸5x5 h fspecial(average, [w w]); % 生成均值滤波核 M1 imfilter(I1, h, replicate); % 对时相一做均值滤波 M2 imfilter(I2, h, replicate); % 对时相二做均值滤波 MR M1 ./ (M2 eps); % 均值比值图逻辑说明imfilter的replicate参数表示边界外采用复制边缘像素的方式填充避免边界区域因为补零产生虚假的低值。窗口尺寸w是均值比值图最重要的参数w太小噪声压不住w太大变化区域被削平。对512像素左右的实验图像w取5~7对真实大场景SAR图像我一般从w9开始尝试观察直方图双峰是否变清晰。2.5 四种差异图选型对比表四种差异图的工程属性对比如下选型直接看这张表差异图类型公式噪声抑制变化方向动态范围典型场景差值图I1 − I2差保留正负适中光学图习惯迁移、对照实验比值图I1 / I2较好丢失需额外判断极易被拉爆SAR变化检测基础方案对数比值图|log(I1/I2)|好统一为正稳定后续要做阈值分割时首选均值比值图mean(I1)/mean(I2)最好丢失稳定噪声极重、变化区域较大补充说明还有一种归一化差异图ND公式为(I1−I2)/(I1I2)本质是差值算子与比值算子的混合SAR文献里偶尔出现但与本次说的四种常见类型有重叠建议先掌握上述四种再按需扩展。工程里我个人的默认选择是对数比值图为底均值比值图做交叉验证差值图只用来展示变化方向。3. 在MATLAB里跑通一套SAR差异图生成流程3.1 数据准备没有真实SAR数据时怎么验证流程很多初学者卡在第一步——没有现成的两时相SAR图像。我不建议随便抓两张光学图硬试因为乘性噪声机制完全不同跑出来的结论没有参考价值。常见做法是用MATLAB自己生成一对带乘性噪声的模拟强度图先跑通流程、建立直觉再替换成自己的真实数据。rng(2024); % 固定随机种子保证结果可复现 N 512; % 模拟图像尺寸 % 模拟真实后向散射底图低频背景 中等幅度起伏 bg 0.6 0.3 * rand(N, N); % 两时相分别叠加乘性斑点噪声0.7~1.3倍随机起伏 I1 bg .* (0.7 0.6 * rand(N, N)); I2 bg .* (0.7 0.6 * rand(N, N)); % 人为制造变化区左上区域后向散射增强模拟新建建筑 I2(180:280, 220:320) I2(180:280, 220:320) * 2.2; % 右下区域后向散射减弱模拟地面塌陷/水体变化 I2(360:420, 80:150) I2(360:420, 80:150) * 0.4; % 转成8位图保存后续走一遍真实读图流程 imwrite(uint8(I1 * 255), sar_t1.png); imwrite(uint8(I2 * 255), sar_t2.png);逻辑说明bg模拟的是地物后向散射的宏观分布乘性噪声用均匀分布模拟虽然真实SAR斑点噪声更接近Gamma分布但均匀乘性模型足够演示差异图的行为差异。两处变化区分别覆盖增强和减弱两种情形可以在后面的差异图里直接验证对数比值图应该对增强和减弱同时敏感。参数说明rng(2024)里的数字可以改成任意整数目的是让随机序列固定便于和别人讨论时复现同一组结果。噪声幅度0.7~1.3倍、后向散射增强2.2倍这些数值是我调过的一组比较直观的配置如果你想让差异图更难或更容易调大噪声幅度或调小变化倍率即可。3.2 主流程把四种差异图一次生成并归一化对比模拟数据准备好以后主流程就是读图、运算、归一化、拼接显示。注意最后的montage只能显示[0,1]范围的图像所以必须对每种差异图做mat2gray归一化否则数值范围差异太大结果显示全黑或全白。% 读图并转浮点 I1 im2double(imread(sar_t1.png)); I2 im2double(imread(sar_t2.png)); % 1) 差值图 D I1 - I2; % 2) 比值图 R I1 ./ (I2 eps); % 3) 对数比值图 LR abs(log((I1 eps) ./ (I2 eps))); % 4) 均值比值图 w 5; h fspecial(average, [w w]); M1 imfilter(I1, h, replicate); M2 imfilter(I2, h, replicate); MR M1 ./ (M2 eps); % 统一归一化并拼接显示 maps {mat2gray(D), mat2gray(R), mat2gray(LR), mat2gray(MR)}; montage(maps, Size, [1 4]);逻辑说明mat2gray做的是min-max归一化把每张差异图的最小值映射到0、最大值映射到1。注意整个流程里每个差异图是独立归一化的所以四张图的「明暗」只表示各自内部的相对差异不能跨图比较像素值大小——这是后续踩坑章节里的高频问题。参数说明w5是最常用起始值。图里能明显看到差值图背景噪声颗粒感最重比值图因为个别极值导致整体偏暗对数比值图背景最干净、两个变化区域最清楚均值比值图边缘略糊但背景噪声最少。这就是四种差异图在视觉上最直观的区别。3.3 大图内存问题single精度与分块生成真实SAR图像动辄上万像素前面这个流程如果直接套用内存会非常紧张。算一笔账10000×10000的图像用double存储单张需要10000×10000×8字节≈800MB。主流程里I1、I2、D、R、LR、MR同时存在就是6×800MB≈4.8GB还没算中间变量普通笔记本直接卡死。常见做法有两个一是全程用single精度存储内存减半代价是精度损失对差异图后续阈值分割几乎无影响二是分块计算只保留当前块的结果。% 法一全程 single 精度内存减半 I1s single(imread(sar_t1.png)) / 255; I2s single(imread(sar_t2.png)) / 255; LRs abs(log((I1s eps) ./ (I2s eps))); % 法二分块计算对数比值图逐块写结果不保留全部中间变量 LR_block zeros(size(I1s), single); bs 1024; % 块大小1024x1024 for r 1:bs:size(I1s, 1) for c 1:bs:size(I1s, 2) r2 min(r bs - 1, size(I1s, 1)); c2 min(c bs - 1, size(I1s, 2)); block1 I1s(r:r2, c:c2); block2 I2s(r:r2, c:c2); LR_block(r:r2, c:c2) abs(log((block1 eps) ./ (block2 eps))); end end逻辑说明分块循环在MATLAB里可读性最好块大小bs决定了内存峰值。bs1024时每个临时块只有1024×1024×4字节≈4MB可以忽略不计。如果你熟悉blockproc也可以用它改写但多通道输入的处理函数回调比较绕我倾向于用这个双层循环逻辑一目了然以后改成GPU并行也容易。参数说明bs不是越大越好它只影响循环次数和临时内存不影响结果。对2万像素级图像bs取2048即可没必要为了减少循环次数把块设得过大。注意这段代码里分块只针对对数比值图其他三种差异图如果需要也可以照此扩展。4. SAR差异图踩坑记录五个常见的「现象-原因-解决」4.1 数值坑uint8溢出、除零出NaN、log(0)出-inf现象一uint8图像直接相减负值全部变成0差值图上增强区域是亮点减弱区域一片死黑完全看不到负变化信息。原因MATLAB的uint8是无符号整型负数被截断到0。解决任何算术运算前先im2double转换运算完再决定要不要转回显示格式。这属于基本功但我在实际项目里见过不止一次因为读图默认就是uint8。现象二比值图出现成片的NaN或白条。原因SAR阴影区、静水面经常有后向散射接近0的像素分母为0时结果变成InfInf参与后续运算就变成NaN。解决给分母加eps还不够建议在比值前先把两期图像整体抬升一个小的偏移量比如minVal min([I1(:); I2(:)])然后用I1-minVal0.01这类方式避免零值。偏移量越大比值越钝感需要根据数据分布权衡。现象三log计算出现-inf直方图最左边出现一个诡异的尖峰。原因log(0)直接给出负无穷mat2gray时-inf会破坏整个归一化范围。解决log前给分子分母都加eps或者用max(I, eps)统一截断。如果在计算后发现差异图里有NaN先排查是不是这一步没加保护。4.2 工程坑配准残差、斑点噪声、归一化玄学现象四差异图里建筑物边缘出现一圈细线状假变化看起来像轮廓描边。原因两个时相图像存在亚像素配准残差强边缘在作差时被梯度放大。这在高分辨率SAR数据里非常典型。解决算差异图之前先对两期图像做配准常见做法是手动选控制点之后用仿射变换对齐更快的工程做法是直接检查配准偏移量如果偏移在1~2像素内对差异图做一次3×3中值滤波把细线磨掉。配准不解决的均值比值图也会带重影只是比差值图轻一点。现象五同一张差异图用mat2gray和用0.1%~99.9%分位截断后再归一化显示出的变化区域大小完全不同。原因min-max归一化对离群值极度敏感个别极亮的噪声点会把整体动态范围拉宽真实变化区被压到暗部。解决不要直接mat2gray先做裁剪。我常用的写法是percentileClip mat2gray组合function out clipNormalize(im, low, high) p_low prctile(im(:), low); % 低分位阈值默认1 p_high prctile(im(:), high); % 高分位阈值默认99 im(im p_low) p_low; im(im p_high) p_high; out mat2gray(im); % 裁剪后再归一化 end习惯上low取1、high取99噪声特别重时改成0.5和99.5。这个函数只影响显示和后续阈值分割的输入范围不影响差异图本身的信息量但能让不同数据的可视化结果可比。5. 差异图质量不看图说话直方图、分离度与可视化验证5.1 直方图双峰才是好差异图光靠肉眼判断四种差异图的优劣很容易被显示参数误导。更可靠的快速检验是看直方图。理想差异图的直方图应该呈现双峰主峰集中在0附近未变化区域次峰在变化区域一侧或一侧拖尾。用MATLAB一行就能对比figure; subplot(2,2,1); histogram(D, 64); title(差值图); subplot(2,2,2); histogram(log(R), 64); title(比值图(log)); subplot(2,2,3); histogram(LR, 64); title(对数比值图); subplot(2,2,4); histogram(MR, 64); title(均值比值图);注意比值图要加log再看直方图否则动态范围被几个极值拉扯根本看不出分布形态。对数比值图的直方图通常最对称、双峰最明显就是因为log把乘性噪声转成了加性。5.2 用Otsu阈值和类间方差给四种差异图打分直方图不够量化。工程上我常用无监督可分离性作为差异图质量的评分先用graythresh求Otsu阈值把差异图分成变化/未变化两类再算两类之间的类间方差得分越高说明差异图越容易划分。脚本形式如下function score sepScore(im) im mat2gray(im); % 先归一化 t graythresh(im); % Otsu 阈值 mask im t; mu0 mean(im(~mask)); % 未变化类均值 mu1 mean(im(mask)); % 变化类均值 var0 var(double(im(~mask))); % 未变化类方差 var1 var(double(im(mask))); % 变化类方差 score (mu1 - mu0).^2 / (var0 var1 eps); end这个函数需要保存成sepScore.m文件或者放在R2016b及以后版本的脚本末尾作局部函数。分离度越高后续阈值分割越不容易踩边界。实际跑模拟数据时对数比值图和均值比值图的得分通常明显高于差值图这就是「眼见为虚、指标为实」的典型例子。5.3 一个固定习惯把差异图作为深度模型的第三通道差异图的价值没有被深度学习替代而是变成了廉价高效的辅助输入。我现在做两时相变化检测项目固定会生成一张归一化后的对数比值图与原图一起拼成三通道输入喂给分割网络收敛速度比单纯输入两期原图更快因为网络不需要自己去隐式学习差异表征。日常出图我也建议用parula而不是jetimagesc(LR); colormap(parula); colorbar; axis image;jet的红色会让人本能地以为「强变化危险」干扰对差异图的实际判断parula和gray更客观也方便后续标注。习惯养成以后拿到任何新数据我都会先把两期图的对数比值直方图打出来看一遍再决定用哪种差异图做主通道。这套流程投入低、可解释性强值得作为你SAR变化检测项目的第一级预处理固定下来。希望帮到你。本文还有配套的精品资源点击获取