简介遥感影像变化检测经典算法资源包聚焦 IR-MAD、MAD、CVA、PCA 四种常用方法面向遥感、测绘、地理信息领域的科研人员和工程师帮助快速建立从多时相影像预处理到变化图生成的完整流程。压缩包共 97 个文件、约 10.4MB以 Matlab 脚本为主附有 bmp/tif 示例影像、ENVI hdr 头文件、fig 结果图其中 IRMAD_Update、MADGet、CVADemo、PCADemo 等算法脚本均可直接运行便于对比各方法在光照、大气条件差异下的表现。已有 2008 人学习下载适合系统了解经典算法原理、在真实数据上对比检测效果或为具体课题筛选合适方法的研究者。借助泰州 TM 两期影像示例可看到多时相异常检测、变化向量分析、主成分降维和比值均值差分各自生成的中间变量与变化图并可直接替换自己的遥感影像执行实验节省从零实现的时间。1. 两期影像摆在你面前变化检测为什么绕不开这组老算法两期同一区域的影像摆在桌上人工对比通常要小半天而 IR-MAD 这类经典算法只需几分钟就能把变化斑块圈出来。遥感影像变化检测回答的问题很直接这两个时相之间地表到底哪里变了、变成了什么。标题里的四个算法——IR-MAD、MAD、CVA、PCA——是这一领域最经典的组合它们不需要任何标注样本只靠两期影像本身就能输出差异图也正因为这个特性它们至今仍被当作深度学习变化检测方法的对照基线。适合做耕地监测、违建排查、灾后评估的从业者。一个反直觉结论是PCA 这种看似做降维的算法在变化检测里恰恰被用来找最有价值的差异方向而不是压缩掉它们。2. 算法谱系先立住CVA、PCA 与 MAD 各自在解决什么这组算法经常被放在一起比较但它们的定位完全不同。CVA 是最朴素的逐波段差值思路PCA 提供了一种线性变换视角而 MAD 则是为“不变关系”专门设计的统计建模。把它们拆开看清各自的前提假设才不会在实跑时用错方向。2.1 CVA逐波段做差后变化强度和方向怎么读CVAChange Vector Analysis的操作非常简单把两期影像对应波段的灰度值相减得到一组差值向量。差值向量的模长代表变化强度向量夹角代表变化类型。比如红光波段从 120 变到 60、近红外从 90 变到 150这个差值向量指向的方向就暗示植被可能发生了变绿或变枯。实跑时有两个问题绕不开。第一是阈值怎么定常见做法是算所有像元差值模长的均值与标准差取“均值加 k 倍标准差”作为变化阈值。k 通常在 1.5 到 3 之间具体取多少得看影像方差。第二是辐射差异的干扰如果两期影像分别来自不同传感器或不同季节大气条件、太阳高度角不一样差值向量里会混入大量非地表变化成分。CVA 对这种全局性偏移几乎没有抵抗能力这是它最明显的边界。不过 CVA 的价值在于解释性。它能把“变了”细化成“从什么方向变到什么方向”在土地覆盖转移分析里很受欢迎。我一般会用它做第一步粗筛后续再用其他算法复核。2.2 PCA主成分分析不是用来降维而是用来分离差异信息PCA 在变化检测里有两条常见路线。一条是先把两期影像逐波段做差再对差值影像做主成分分析取前几个主成分作为主要变化分量另一条是把两期影像的所有波段堆叠成一个多波段影像一次性做主成分分析认为前几个主成分捕获了两期共有的稳定背景信息排在后面的主成分则更多地体现时相差异。这里要理解一个关键点PCA 的各个主成分是沿着方差最大方向排列的。两期影像里不变的地物通常灰度分布高度相关它们会集中体现在前几个主成分里变化像元在统计上属于少数派方差贡献小反而落到后面的成分中。所以做变化检测时常见的做法是丢掉前几个主成分拿后面的成分来合成差异图。这个思路和人脸识别里“特征脸”的典故一致——PCA 的基向量是被数据驱动出来的特征方向只是变化检测里我们关心的是那些方差小却语义明确的尾巴。但 PCA 有一个结构性弱点它是全局变换。整幅影像共享同一组特征向量局部区域的差异会被全局统计平均掉。如果变化区很小PCA 的效果就会变差。2.3 MAD用典型相关分析给“不变关系”建模MADMultivariate Alteration Detection的思路比 CVA 更进一层。它不是直接在原始波段空间做差而是先对两期影像分别做线性组合得到两个“典型变量” U aᵀX、V bᵀY然后计算它们的差值作为变化分量。为什么要绕这一圈因为 CVA 假设两个时相同一波段的数值可以直接相减。但现实中两期影像之间存在传感器定标差异、大气路径辐射差异直接相减会把系统性偏差当变化。MAD 用典型相关分析CCA去寻找两期影像之间最相关的线性组合找到之后再做差这样得到的差值在统计意义上更接近真实的异常扰动。数学上MAD 分量的方差等于 2(1−ρₖ)其中 ρₖ 是第 k 对典型变量之间的相关系数。相关系数越低说明这组线性组合在时序上越不一致也就是变化信息越突出。因此 MAD 分量通常按方差从大到小排列排在前面的分量包含最有辨识力的变化信号。相比 CVAMAD 对辐射偏移和波段间线性关系的干扰明显更耐受这是它在实操中最受欢迎的原因。2.4 四个算法的适用边界与选型对比算法输入形式核心原理输出结果适用场景CVA两期多波段影像逐波段差值向量变化强度图 变化方向角土地覆盖转移分析、快速粗筛PCA差值影像或堆叠影像协方差矩阵特征分解若干主成分中的差异分量全局变化模式探索、波段压缩MAD两期多波段影像典型相关分析 线性组合差若干 MAD 分量图多时相辐射不一致时的稳健检测IR-MAD两期多波段影像MAD 迭代加权逼近不变像元加权的 MAD 分量与变化概率图高精度无监督变化检测、辐射归一化选型建议很直接如果两期影像来自同一传感器、同季相、经过严格辐射定标CVA 就够用如果影像来源复杂、辐射差异明显直接用 MAD如果追求更高精度且有耐心调参数上 IR-MAD。PCA 更多是作为预处理或辅助手段出现很少单独承担最终判定。3. 跑通 MAD 最小实现数据准备、核心推导与差异图生成MAD 的原理并不复杂但真正跑通需要处理影像读写、矩阵展平、奇异值分解、差异图重排这些环节。这一章给出一套可以用 GDAL NumPy SciPy 完整复现的最小实现不依赖 ArcGIS 或 ENVI 的现成工具箱。3.1 工具选型为什么用 GDAL NumPy 而不是 ArcGIS 工具箱ArcGIS 和 ENVI 里都有变化检测工具但对生产流程有四个不友好之处一是批处理需要写 Model Builder调试不直观二是中间矩阵不透明出问题难以定位三是许可证环境对自动化部署不友好四是很难把算法嵌入到自定义的后处理管线里。用 GDAL 负责影像 IONumPy 做数组运算SciPy 做统计分布计算整个逻辑全都摊在代码里每一行都能被检查。此外公开数据集比如 Onera Satellite Change Detection 数据集都以 GeoTIFF 格式提供GDAL 是读取这些数据最通用的方式。下面的示例默认已安装 GDAL、NumPy、SciPy 且两期影像已配准到同一网格。3.2 数据准备对齐波段、掩膜无效值、展平样本矩阵MAD 的输入是两期影像各自的像元矩阵。预处理的目标有二一是把影像数组展平成“样本×波段”的矩阵方便后续协方差计算二是把无值区、云区、边界区排除掉避免污染统计量。from osgeo import gdal import numpy as np def read_to_matrix(path, valid_min0, valid_max65535): ds gdal.Open(path) arr ds.ReadAsArray() # shape: (band, row, col) ds None if arr.ndim ! 3: raise ValueError(需要多波段影像) arr arr.astype(np.float32) b, h, w arr.shape # 构建有效像元掩膜所有波段都在有效灰度范围且不是 NaN mask np.isfinite(arr[0]) for band in arr: mask (band valid_min) (band valid_max) # 展平成 (N, B) 矩阵N 为有效像元数 samples arr[:, mask].T # shape: (N, B) return samples, mask, (h, w) X, mask1, _ read_to_matrix(t1.tif) Y, mask2, _ read_to_matrix(t2.tif) # 两期影像有效像元交集 mask mask1 mask2 X X[mask] # 严格对齐后行列号一致直接用布尔掩膜筛选 Y Y[mask] print(有效像元数:, X.shape[0], 波段数:, X.shape[1])这段代码里有三个关键点需要说明。ReadAsArray()返回的通道顺序是 (band, row, col)不要和 OpenCV 的 HWC 顺序搞混。astype(np.float32)是必需的原始影像常以 UInt16 存储直接用整型求协方差时精度损失很大尤其遇到灰度值 0 和 1 附近的小数值时会严重失真。最后用有效像元交集统一两期矩阵是为了保证下一步协方差计算基于同一批像元位置。3.3 核心推导标准化、SVD 与典型变量的关系MAD 的数学核心在于典型相关分析。这里用一个实现技巧先把两期矩阵分别标准化为零均值、单位方差这样自协方差矩阵变成单位矩阵交叉协方差的 SVD 结果就直接给出两组投影方向。def mad_components(X, Y): # X, Y: (N, B)列对应波段 n, b X.shape # 对每一波段做 z-score 标准化 x_mean X.mean(axis0) y_mean Y.mean(axis0) x_std X.std(axis0) y_std Y.std(axis0) Xn (X - x_mean) / (x_std 1e-8) # 避免除零 Yn (Y - y_mean) / (y_std 1e-8) # 交叉协方差矩阵 C (Xn.T Yn) / (n - 1) # SVD左右奇异向量就是 CCA 的投影方向 U, s, Vt np.linalg.svd(C) # 投影得到典型变量 u Xn A.T, v Yn B.T A U.T # shape (b, b) B Vt # 注意 Vt 已经是转置后的 V u Xn A.T # (n, b) 每列是一个典型变量 v Yn B.T # MAD 分量对应列的差 mad u - v # 归一化每个 MAD 分量的理论方差是 2(1 - rho_k) for k in range(b): rho s[k] mad[:, k] mad[:, k] / np.sqrt(max(2 * (1 - rho), 1e-8)) return mad, s mad, rho mad_components(X, Y)这段代码是把 CCA 求解简化为一次 SVD 的关键所在标准化让两个自协方差矩阵变成单位阵CCA 的广义特征分解退化为普通 SVD。U的每一列是 X 侧投影方向的转置Vt的每一行是 Y 侧投影方向。奇异值s[k]就是第 k 对典型变量的相关系数它越接近 1说明这对变量在两期影像之间越一致对应 MAD 分量的方差越小、变化信息越弱。最后一步对每个 MAD 分量除以sqrt(2(1 - rho))很微妙。这一步的作用是把分量方差统一归一化到 1为下一步用卡方分布做统计检验做准备。如果你不做显著性检验而只是想肉眼看图这步可以省略但后续要算变化概率时必须保留。3.4 差异图生成MAD 分量平方和与阈值选取MAD 生成的是若干个差异分量。可以把多个分量平方求和得到一个近似服从卡方分布的总统计量再用生存函数映射成每个像元的“变化概率”。这是 MAD 能直接产出概率图的关键一步。from scipy import stats def mad_change_probability(mad, k_comp3): # 取前 k_comp 个方差最大的 MAD 分量 m mad[:, :k_comp] chi2_stat np.sum(m ** 2, axis1) # 卡方统计量 # 越小表示变化越显著转成“变化概率”便于可视化 p_nochange stats.chi2.sf(chi2_stat, dfk_comp) p_change 1.0 - p_nochange return p_change p_change mad_change_probability(mad) # 把概率向量摆回二维图 h, w mask.shape change_map np.full((h, w), np.nan) change_map[mask] p_changek_comp取 2 或 3 是常见经验值MAD 分量按方差降序排列前几个分量集中了主要变化信息太靠后的分量基本是噪声。卡方检验的自由度等于分量数直接决定概率分布的形态。拿到了概率图阈值怎么定先用一个简单试验在图上叠加直方图观察概率分布是否呈双峰。如果像元的概率集中在 0 附近和 1 附近阈值取两峰之间的谷底即可如果分布连续那就得人为接受一个代价权衡比如把前 5% 概率最高的像元判为变化。这块属于经验区不同影像的阈值差异很大后面避坑章会展开讲。4. 从 MAD 到 IR-MAD迭代加权的实现与参数经验MAD 一次性算完得到的结果往往不够干净那些真正的大面积变化区会像杠杆点一样拉扯协方差矩阵导致投影方向被“带偏”。IR-MADIteratively Reweighted MAD的补救思路很优雅——每次迭代都重新估计哪些像元更可能没有变化给这些像元更高权重再重新计算投影方向。4.1 为什么 MAD 会被强变化污染IR-MAD 怎么补救设想一个场景某块区域发生了大范围森林砍伐两期影像在同一位置的光谱响应差异非常大。MAD 在计算协方差时会给这些像元同等的统计地位于是协方差矩阵会被这批高强度差异牵着走最终算出来的投影方向不再是“不变背景下的最大差异”而是“背景差异与砍伐差异的混合体”。IR-MAD 的做法是给每个像元一个权重 wᵢ权重代表这个像元属于“未变化”的概率。第一轮所有像元等权等价于普通 MAD算出 MAD 分量后用卡方分布把每个像元的差异平方和映射成一个概率值差异越小的像元概率越高下一轮带着这些概率重新计算加权协方差和投影方向反复迭代直到权重分布不再明显变化或达到最大轮数。这个过程本质上是把“不变像元”的辨识从一个一次性的全局假设变成一个逐步逼近的迭代估计。每一轮协方差矩阵都更偏向于不变背景投影方向也越来越精准。这是 IR-MAD 比 MAD 精度高的根本原因。4.2 权重计算卡方分布尾概率与 mad 归一化权重计算的输入是当前轮的 MAD 分量。每个像元有 B 个分量对分量平方求和得到一个标量 qᵢ。如果该像元真的未变化qᵢ 应服从自由度为 B 的卡方分布如果变化强烈qᵢ 会落在分布的极端尾部。基于这个假设用卡方分布的生存函数把 qᵢ 映射成概率w_i chi2.sf(q_i, dfB)qᵢ 越小生存函数值越接近 1该像元被视为不变像元的可信度越高qᵢ 极大时生存函数值趋近 0其统计权重也趋近 0。这轮概率就是下一轮的权重。“mad 归一化”这个词在 IR-MAD 的实现里通常指两层操作一层是对输入矩阵做加权 z-score 标准化让协方差计算在统一尺度下进行另一层是对每个 MAD 分量除以理论标准差使分量平方和方差匹配卡方分布的自由度。这两层少了一层卡方检验的统计性质都会崩掉。当年我第一次实现时漏掉分量的方差归一化概率图几乎是一片纯黑或纯白就是这里出了问题。4.3 IR-MAD 完整迭代代码与收敛判断下面给出一个带完整迭代逻辑的 IR-MAD 实现注意看权重如何在协方差计算和标准化两个位置同时生效。def ir_mad(X, Y, n_iter10, tol1e-3): n, b X.shape w np.ones(n) / n # 初始等权 mad_final None for it in range(n_iter): # 1. 用当前权重做加权标准化mad归一化 sum_w w.sum() x_mean (X * w[:, None]).sum(axis0) / sum_w y_mean (Y * w[:, None]).sum(axis0) / sum_w Xc X - x_mean Yc Y - y_mean # 2. 加权协方差 Cxx (Xc * w[:, None]).T Xc / sum_w Cyy (Yc * w[:, None]).T Yc / sum_w Cxy (Xc * w[:, None]).T Yc / sum_w # 3. 加权 CCACholesky 白化后做 SVD Lx np.linalg.cholesky(Cxx 1e-6 * np.eye(b)) Ly np.linalg.cholesky(Cyy 1e-6 * np.eye(b)) Lx_inv np.linalg.inv(Lx) Ly_inv np.linalg.inv(Ly) K Lx_inv Cxy Ly_inv.T U, s, Vt np.linalg.svd(K) # 4. 投影方向变换回原坐标 A Lx_inv.T U # 列向量为 X 侧投影方向 V Vt.T B Ly_inv.T V # 列向量为 Y 侧投影方向 # 5. 计算 MAD 分量并做方差归一化 u Xc A v Yc B mad u - v for k in range(b): mad[:, k] mad[:, k] / np.sqrt(max(2 * (1 - s[k]), 1e-8)) mad_final mad # 6. 更新权重卡方分布生存函数 chi2_stat np.sum(mad[:, :b] ** 2, axis1) w_new stats.chi2.sf(chi2_stat, dfb) # 截断极小权重防止协方差计算时出现病态 w_new np.clip(w_new, 1e-6, 1.0) # 7. 收敛判断权重总体变化小于阈值 diff np.abs(w_new - w).sum() / n w w_new if diff tol: print(f第 {it1} 轮收敛, 权重变化量 {diff:.6f}) break return mad_final, w mad_ir, weight ir_mad(X, Y)收敛逻辑说明第 1 步的加权标准化确保均值估计不被强变化像元拉偏这是许多简化实现遗漏的地方。第 3 步加1e-6的正则项是为了防止当某个波段在掩膜后样本方差趋近 0 时 Cholesky 分解崩溃。第 6 步用当前轮 MAD 分量直接计算权重实现的是典型的 EM 式迭代先固定权重估计投影方向再固定投影方向更新权重。参数方面tol1e-3表示平均每个像元的权重相对变化小于千分之一即停止n_iter10是上限保护防止异常影像导致不收敛死循环。实际经验中多数影像在第 4 到第 6 轮收敛。4.4 三个必调参数迭代次数、收敛阈值与权重截断迭代次数最常见取 5 到 10。超过 10 轮后权重分布通常会进入“不动点”继续迭代对结果几乎没影响。但如果影像中有大量强变化区前几轮权重会把变化像元压得极低导致协方差矩阵只反映纯背景反而可能丢失真实而又微弱的渐变信号。因此我一般不建议一次拉满 20 轮而是先用 5 轮看中间结果必要时再续跑。收敛阈值 tol 默认1e-3即可。调大到这个值的 10 倍会加快速度但权重可能没坐实调小可能出现轮次耗尽也不收敛的情况。权重截断1e-6的作用是防止某些像元权重被更新成绝对 0从而彻底退出后续统计留下数值隐患。更保守的做法是截断到1e-4这时边缘像元对协方差仍有微小贡献适合变化区特别破碎的影像。还有个容易被忽略的参数是卡方自由度。理论上自由度等于参与计算的波段数 b但有些实现会只取前几个独立 MAD 分量进入卡方统计此时自由度要改成实际取用的分量数。两者混用会直接导致权重整体偏移概率图阈值完全失真。5. 遥感变化检测算法避坑五个让结果翻车的常见问题这一章的价值直接来自踩坑现场。以下五个问题是我在多个项目里反复遇到的情况几乎涵盖了 MAD/IR-MAD 类算法最容易翻车的位置。5.1 两期影像辐射不一致大面积假变化怎么排查现象结果图上某块连片区域被整体判为变化但实地核查发现地物根本没变。原因两期影像大气条件、太阳高度角或传感器定标参数不同导致同一地物在两个时相的辐射值系统性偏移。IR-MAD 虽然对线性辐射偏移有抵抗力但如果偏移是非线性的比如大气水汽对不同波段影响不同MAD 分量仍会残留大量假信号。解决先做相对辐射归一化。常见做法是选择研究区内的伪不变特征比如深水湖泊的深水区、大片裸岩、稳定不透水面用它们拟合两期影像的线性回归关系再对第二期影像做校正。IR-MAD 迭代结束后的权重图里权重接近 1 的像元恰恰可以作为伪不变特征的自动候选这是一条闭环路线。5.2 协方差矩阵奇异波段数多于有效样本数时怎么处理现象程序在np.linalg.cholesky或np.linalg.svd处报错提示矩阵不是正定。原因掩膜后有效像元数大于波段数很多时一般不会出错真正出错通常是某个波段动态范围几乎为 0或者使用了超高光谱数据时波段间高度共线导致自协方差矩阵接近奇异。解决在 Cholesky 分解前给协方差矩阵加上一个小的正则项代码里的Cxx 1e-6 * np.eye(b)就是干这个的。如果加了正则仍然报错优先检查波段选择——把相关性极高的邻接波段删掉或用 PCA 预压缩是更治本的做法。5.3 椒盐噪声严重变化图里的孤立点怎么清理现象结果图上散布大量孤立像元单看每一处都像变化但周围背景完全稳定。原因传感器噪声、配准亚像元误差、影像拉伸导致的小幅灰度抖动都会在 MAD 分量的高维空间里表现为统计显著的变化。IR-MAD 的权重迭代对这种高频噪声并不敏感因为它们的统计特征是空间不相关。解决在概率图上做空间后处理最常见的手段是中值滤波、众数滤波配合最小图斑面积约束。给一个直观经验变化检测成果图斑面积小于 3×3 像元的图斑在网格尺度下往往没有制图意义应当合并或剔除。5.4 阈值凭经验拍脑袋精度虚高的原因与矫正现象取某个阈值后目视效果和验证点精确率的数字都很漂亮但模型换到相邻区域后立刻失效。原因阈值是在单一影像上优化的天然过拟合。MAD 输出的概率分布形态在不同影像之间差异很大有的呈 U 形、有的呈单峰长尾不存在一个通用阈值。解决用 Otsu 方法在总概率分布上自动寻找分割点或使用两成分高斯混合模型拟合概率分布以两个高斯成分的交点为阈值。更稳妥的做法是用已有 Ground Truth 上的变化比例反推阈值让阈值对应的变化面积与先验比例一致。这一条是精度评估里最后一道保险也最常被忽略。5.5 配准误差导致边缘条带后处理掩膜怎么加现象房屋、道路、山脊线的边缘出现沿地物轮廓的平行条带看起来像地物“重影”。原因两期影像配准存在一两个像元的平移误差导致地物边界两侧的灰度差被算法判定为变化。MAD 对这种高频空间位移非常敏感这是统计方法普遍绕不开的结构性缺陷。解决先查看配准残余误差报告大于 0.5 像元的必须重做配准后处理阶段可以对变化概率图进行形态学开运算去除细长条带再用边缘掩膜把高梯度区域剔除掉防止线性地物边缘被误报。代价是真实沿边界发生的小幅扩展变化也会被误删这个取舍要和业务方确认清楚。6. 结果验证与进阶技巧用混淆矩阵和 Kappa 评估变化检测无监督变化检测算法最容易出现的幻觉是“看起来不错但精度不可证”。IR-MAD 输出的概率图必须经过量化评估才能进入生产流程。下面给出一个最小验证脚本搭配两个我在实践中沉淀下来的技巧。from sklearn.metrics import confusion_matrix, cohen_kappa_score # change_truth: 0/1 验证标签; p_change: 算法输出的变化概率 threshold 0.7 change_pred (p_change threshold).astype(np.int16) # 全图二值化 valid (change_truth 0) # 只评估有标签的像元 cm confusion_matrix(change_truth[valid], change_pred[valid], labels[0, 1]) tn, fp, fn, tp cm.ravel() oa (tn tp) / cm.sum() kappa cohen_kappa_score(change_truth[valid], change_pred[valid]) f1 2 * tp / (2 * tp fp fn) print(fOA{oa:.3f} Kappa{kappa:.3f} F1(change){f1:.3f})阈值threshold不是拍脑袋定的。我会先直方图画一遍p_change的分布如果双峰明显阈值放在两峰谷底如果单峰就取排序后的 95 分位点再微调。Kappa 系数比总体精度 OA 更值得盯因为在变化检测里“无变化”通常占绝大多数一个永远全判无变化的模型也能拿到 90% 的 OAKappa 能把这个幻觉直接打到 0 附近。F1 则代表对变化类别的查准与查全的折中业务上更符合实际诉求。另一个进阶技巧是把 IR-MAD 迭代结束时权重图的倒数当作“变化敏感度”的先验把它与概率图做加权融合。权重越低的像元在迭代中被证明越是偏离不变背景的强变化点直接给这些位置的概率加一个小增量可以缓解阈值分割时低对比度变化被漏判的问题。我最早吃过一次亏某次水域淹没了大片农田由于水体信号强迭代权重把水边渐变区的像元压得很低概率阈值一卡整片渐淹区被漏判。后来我就养成了一个习惯——只看最终概率图之前先单独检查 IR-MAD 每一轮权重图的变化趋势确认是否出现局部区域权重陡降。实操上还有一条值得养成惯性在跑 IR-MAD 之前把两期影像的直方图和均值差打印出来差距超过一个标准差时先做相对辐射归一化不要硬跑。这套做法陪我扛过了多次城市扩张、森林扰动和灾后评估项目。希望帮到你。本文还有配套的精品资源点击获取