简介面向遥感影像变化检测研究的一套经典算法实现包涵盖基于多变量统计的IR-MAD、MAD变化向量分析的CVA以及主成分分析的PCA四种主流算法适合测绘、遥感及相关专业学生和工程师用于地表变化监测、土地利用变化分析等场景。压缩包共九十七个文件大小约十点四兆其中既有Matlab源码与可直接运行的示例脚本也有实验影像、结果图表及ENVI头文件等配套数据目录层次清晰便于按模块调用和二次开发。目前已有两千零一十四人学习下载。资源提供统一的算法入口和多种对比脚本可一键运行各算法并输出变化强度图、二值图等中间结果帮助读者结合泰州地区两期真实影像理解不同算法在光照、大气影响下的表现差异以及参数设置对检测效果的影响。掌握这些经典算法可帮助科研与工程人员更准确地识别地表变化服务于土地利用、城市规划和环境监测等领域。1. 为什么老算法至今仍是遥感变化检测的起手式遥感影像变化检测看似已经有大量深度学习方法但真到了生产项目里面对多时相、多传感器、辐射不一致的影像绝大多数团队第一版跑通的依然是经典算法。IR-MAD、MAD、CVA、PCA这四兄弟之所以被反复提起不是因为它们论文年代久远而是因为它们不依赖标注样本、计算可控、结果可解释。尤其是IR-MAD在辐射归一化和变化检测两个任务上都是公认的基线比很多花哨的模型稳定得多。这套算法的核心思路并不复杂找两个时相影像之间的统计差异。CVA直接用变化向量长度衡量差异强度PCA先把高维波段压缩成几个主成分再求差MAD和IR-MAD则是通过典型相关分析找到“变化最不明显”的投影方向再做差分。实际落地时它们的适用场景有明确分工CVA适合多波段影像快速筛查PCA适合变化信息集中在前几个主成分的数据IR-MAD则是对辐射差异大、噪声多的影像最有鲁棒性的选择。本文按从易到难的顺序把这四个算法拆开讲清楚给出可运行的代码和参数建议最后集中说坑。2. CVA与PCA先上手变化向量分析的两条捷径2.1 CVA的原理变化强度和变化方向各是什么CVAChange Vector Analysis是变化检测领域最直观的方法。它的思路是把每个像素在不同时相的多波段光谱值看作两个向量两个向量之差就是“变化向量”。变化向量有长度和方向两个属性长度代表变化强度方向代表变化类型。比如植被变裸土和耕地变水体虽然变化强度可能接近但光谱变化的方向差异很大方向可以辅助区分变化类别。数学表达很简单。设时相1像素光谱向量为 X1(x1_1, x1_2, ..., x1_b)时相2为 X2(x2_1, x2_2, ..., x2_b)b是波段数。变化向量为D X2 - X1变化强度 ||D|| sqrt(sum((x2_i - x1_i)^2))变化方向则由D的分量比值决定。实际使用时每个像素算出一个强度值得到一张单波段灰度图再通过阈值分割变成变化/不变二值图。这里有一个容易忽略的前提CVA对辐射一致性要求极高。如果两期影像来自不同传感器、不同太阳高度角或不同大气条件未做归一化直接算差值结果会被辐射差异主导真实变化反而被淹没了。因此CVA在实操中通常是配合直方图匹配或后续要讲的IR-MAD辐射归一化一起用的。2.2 用Python实现CVA的完整流程用Python做CVA的最小流程分五步读取两期影像——波段对齐——计算变化向量模长——阈值分割——输出结果。以下是核心代码。import numpy as np from osgeo import gdal def read_image(path): ds gdal.Open(path) band_count ds.RasterCount bands [ds.GetRasterBand(i 1).ReadAsArray().astype(np.float32) for i in range(band_count)] return np.stack(bands, axis-1), ds def cva_change_detection(img1, img2, thresholdNone): # 输入img1、img2形状均为 (h, w, b)要求波段顺序一致 assert img1.shape img2.shape, 两期影像shape不一致检查波段数和行列数 change_vector img2 - img1 magnitude np.sqrt(np.sum(change_vector ** 2, axis-1)) if threshold is None: # 默认用均值 1倍标准差作为阈值 threshold magnitude.mean() magnitude.std() change_mask (magnitude threshold).astype(np.uint8) return magnitude, change_mask img1, ds1 read_image(t1.tif) img2, ds2 read_image(t2.tif) magnitude, change_mask cva_change_detection(img1, img2) # 写输出 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(cva_magnitude.tif, ds1.RasterXSize, ds1.RasterYSize, 1, gdal.GDT_Float32) out_ds.GetRasterBand(1).WriteArray(magnitude) out_ds.SetGeoTransform(ds1.GetGeoTransform()) out_ds.SetProjection(ds1.GetProjection()) out_ds.FlushCache()逻辑说明read_image把两期影像读成(h, w, b)的numpy数组cva_change_detection先做差值再沿波段轴求欧几里得模长。threshold参数如果不给代码会用均值加一倍标准差作自适应阈值。参数说明这里的自适应阈值适合变化面积占比不太高的场景如果变化区域超过30%均值加一倍标准差会漏检大量真实变化此时应该改用OTSU或人工目视选定阈值。另一个关键点是波段对齐——不同传感器的波段设置不同比如Sentinel-2的B4和Landsat-8的B4波长范围不一样简单粗暴直接相减会引入系统性误差。2.3 PCA变化检测为什么先压缩再差分更稳PCA主成分分析用于变化检测的常见做法有两种一种是对两期影像分别做PCA再对前几个主成分做差另一种是把两期影像叠在一起做PCA变化信息会在特定主成分上集中体现。业界主流用的是前者因为物理含义清晰——每个时相的主成分代表该时相的主要光谱变异方向相同序号的PC之间做差就能突出真实变化。和CVA的直接波段相减相比PCA有两个实际优势。第一它对波段间相关性的处理更好高光谱或多光谱影像波段之间往往高度相关直接相减会把冗余噪声叠加PCA把信息压缩到前几个主成分噪声被留在尾部分量里。第二PCA能天然处理波段数量不一致的情况只要两个时相分别做主成分取各自的前k个PC做差即可。但PCA不是万能的。它最典型的翻车场景是“主成分序号不对应”——两个时相的PC1可能代表完全不同的物理含义强行相减反而引入虚假变化。实操中通常先看两个时相的特征值分布和特征向量确认PC的物理对应关系后再做差。2.4 PCA的Python实现从特征值到变化图from sklearn.decomposition import PCA def pca_change_detection(img1, img2, n_components3): h, w, b img1.shape # 把三维影像展开成二维矩阵每行是一个像素 X1 img1.reshape(-1, b) X2 img2.reshape(-1, b) pca1 PCA(n_componentsn_components) pca2 PCA(n_componentsn_components) pc1 pca1.fit_transform(X1) # 形状 (h*w, n_components) pc2 pca2.fit_transform(X2) # 逐主成分做差值再计算综合变化强度 diff pc2 - pc1 magnitude np.sqrt(np.sum(diff ** 2, axis-1)).reshape(h, w) return magnitude, pca1, pca2 # 使用示例 magnitude, pca1, pca2 pca_change_detection(img1, img2, n_components3) # 查看特征值占比判断前几个主成分是否足够 print(t1 前3主成分解释方差比:, pca1.explained_variance_ratio_) print(t2 前3主成分解释方差比:, pca2.explained_variance_ratio_)逻辑说明PCA模型对每个时相独立拟合fit_transform返回降维后的主成分得分将对应序号的PC相减后求模长得到变化强度。n_components的取值需要根据特征值分布来定如果前3个主成分能解释85%以上方差取3就够如果解释率偏低要适当增加分量。参数说明这里有一个在实战中很重要的细节——PCA拟合前要不要做标准化。如果各波段量纲差异大比如有的波段反射率范围0-1有的亮度温度范围280-300K不做标准化会让高方差波段主导主成分建议先做z-score标准化。我一般会加一行X1 (X1 - X1.mean(axis0)) / X1.std(axis0)。和CVA一样PCA变化检测输出的强度图也需要阈值分割才能变成最终的变化区域。阈值选择方法和CVA通用可以用OTSU也可以用均值加标准差。但要注意PCA变化检测的强度值分布往往更集中直接套CVA的阈值经验值容易失效。一个实用的做法是先把强度图做百分位拉伸比如2%-98%再做阈值分割效果会稳定很多。3. MAD与IR-MAD从典型相关到迭代加权3.1 MAD的数学原理找到变化最不明显的方向MADMultivariate Alteration Detection的出发点比CVA和PCA都更讲究。它不直接比较原始波段而是先对两个时相的影像分别做线性组合然后找让“组合后差值方差最大化”的投影方向。用典型相关分析CCA的语言说MAD要找的是两组波段各自的线性组合U和V使得U和V之间的相关性最小——相关性最小意味着这两个组合反映的差异最大也就是“变化最显著”的方向。具体来说MAD定义第k个变化变量为MAD_k a_k^T * X1 - b_k^T * X2其中a_k和b_k是第k对典型向量通过CCA求解。求解过程等价于对两期影像的交叉协方差矩阵做广义特征分解。计算得到的所有MAD分量之间互不相关且方差从大到小排列。实际应用中通常取前几个MAD分量来代表主要变化信息最后一个MAD分量理论上方差最小。MAD的一个重要优势是它对辐射增益和偏置的变化有天然抵抗能力——因为CCA不依赖绝对辐射值而依赖协方差结构。只要两期影像的辐射关系近似线性MAD结果就不会被整体亮度差带偏。这就是为什么MAD常用作辐射归一化的预处理步骤。3.2 IR-MAD改进了什么EM迭代与权重更新MAD有个明显软肋它对噪声和高杠杆点非常敏感。如果影像里有云、阴影或大面积饱和像素协方差矩阵估计会被这些异常值带偏得到的投影方向就失真。IR-MADIteratively Reweighted MAD的解决思路是引入迭代加权——每轮MAD计算后根据当前变化变量的大小给每个像素分配权重变化剧烈的像素大概率是真实变化或噪声权重降低稳定像素权重升高然后用加权后的统计量重新计算MAD循环往复直到收敛。这个机制本质上是在做鲁棒估计。加权函数通常取卡方分布的概率密度w_i 1 - F_chi2(chi2_i)其中F_chi2是卡方分布的累积分布函数自由度等于MAD分量数。权重越接近1说明该像素在各MAD分量上的累积变化量越小越接近0说明变化越大。经过若干轮迭代权重会收敛最终的MAD分量对噪声离群点稳健得多。IR-MAD的实际效果有两个明显体现一是变化检测的虚警率大幅下降二是得到的MAD分量可以直接用来反推两期影像间的线性辐射变换参数因此它被广泛用作辐射归一化的“黄金步骤”——先做IR-MAD选出不变特征点再用这些点做回归比全局直方图匹配可靠得多。3.3 IR-MAD的Python实现核心步骤与收敛判据Python里没有开箱即用的IR-MAD库需要自己组合numpy和scipy完成。实现分六步读数据、CCA初始化、计算MAD分量、卡方权重更新、加权重估、迭代收敛判断。以下是精简的实现骨架。import numpy as np from scipy.stats import chi2 from scipy.linalg import eigh def irmad_change_detection(img1, img2, n_componentsNone, max_iter50, tol1e-3): h, w, b img1.shape X1 img1.reshape(-1, b).T # (b, n_pixels) X2 img2.reshape(-1, b).T n_pixels X1.shape[1] if n_components is None: n_components b # 去掉均值 X1_c X1 - X1.mean(axis1, keepdimsTrue) X2_c X2 - X2.mean(axis1, keepdimsTrue) # 初始权重全为1 w np.ones(n_pixels) mad_list [] for it in range(max_iter): # 按权重计算协方差矩阵 w_sqrt np.sqrt(w) C11 (X1_c * w_sqrt) (X1_c * w_sqrt).T / np.sum(w) C22 (X2_c * w_sqrt) (X2_c * w_sqrt).T / np.sum(w) C12 (X1_c * w_sqrt) (X2_c * w_sqrt).T / np.sum(w) # 求解广义特征值问题得到典型相关系数 # 构建分块矩阵特征值为典型相关系数的函数 M np.linalg.inv(C11) C12 np.linalg.inv(C22) C12.T evals, evecs eigh(M) # 按特征值降序排列特征向量就是a_k idx np.argsort(evals)[::-1] a evecs[:, idx[:n_components]] # (b, n_components) # 由a和典型相关系数求b rho np.sqrt(np.clip(evals[idx[:n_components]], 0, 1)) b np.linalg.inv(C22) C12.T a / rho.reshape(1, -1) # 计算MAD分量 mad a.T X1_c - b.T X2_c # (n_components, n_pixels) mad_list.append(mad.copy()) # 计算卡方统计量并更新权重 mad_var np.var(mad, axis1) chi2_stat np.sum(mad ** 2 / mad_var.reshape(-1, 1), axis0) new_w 1 - chi2.cdf(chi2_stat, dfn_components) new_w np.clip(new_w, 0.01, 1.0) if np.max(np.abs(new_w - w)) tol: w new_w break w new_w # 用最后的MAD分量计算变化强度 mad_final mad_list[-1] magnitude np.sqrt(np.sum(mad_final ** 2, axis0)).reshape(h, w) return magnitude, w.reshape(h, w)逻辑说明第1步先把影像展开成(b, n)矩阵并去均值第2步根据当前权重计算加权协方差矩阵收敛前权重是均匀的第一轮等价于普通MAD第3步求解广义特征值问题得到典型向量a和b这是MAD的核心。第4步计算MAD分量第5步用卡方分布更新权重——这一步是IR-MAD区别于MAD的关键权重越小代表该像素变化越大下一轮协方差估计时它的影响被抑制。第6步用权重最大变化量做收敛判据。参数有三个需要重点说明。n_components一般取b到b-2之间如果影像波段数超过8取4到6就够因为后面的MAD分量噪声主导对变化检测贡献不大。max_iter设30到50足够收敛IR-MAD通常在10轮以内就稳定了超过30轮还不收敛要检查是否出现了不收敛的震荡。tol建议设1e-3这个精度足够设太小只增加计算量不增加精度。注意这个实现做了简化完整IR-MAD还需要在每次迭代中更新X1和X2的均值因为加权后均值会变并把均值项纳入MAD计算。实际操作中均值漂移对结果影响不大但对严谨性有要求时不要省。3.4 IR-MAD的三个典型应用方式IR-MAD在实践中有三种用法。第一种是直接做变化检测把最后的MAD分量求平方和得到变化强度图再阈值分割。这是最经典的应用。第二种是辐射归一化预处理利用IR-MAD迭代收敛后的权重选取权重0.9的像素作为“不变点”用这些像素拟合两个时相之间的线性回归关系把时相2影像校正到时相1的辐射水平然后重新做CVA或PCA。这种做法比全局直方图匹配好很多因为不变点都是真正的稳定地物没有混入水体、植被这些辐射变化大的区域。第三种用法是作为深度学习训练样本的生成器。在没有人工标注的情况下用IR-MAD得到一个初步变化图置信度高的正负样本可以直接作为半监督学习的种子。这个方法在生产项目里很实用。4. 阈值分割四类算法殊途同归的最后一道关4.1 变化强度图到变化掩膜的三种阈值策略无论是CVA、PCA还是IR-MAD最后一步都是把连续的变化强度图变成一张“变/不变”的二值图。这一步的阈值选择不当前面积累的效果会被全部抵消。工程上常用三种策略。第一种是均值加k倍标准差。这个策略适合变化面积占比小且强度分布近似正态的数据。优点是计算简单缺点是当变化面积大或强度分布偏态时误差大。第二种是OTSU大津法它自动寻找使类间方差最大的阈值不依赖正态假设对双峰分布效果非常好。第三种是双阈值法先取一个高阈值确定“确定变化”区域再取低阈值确定“疑似变化”区域中间的灰度通过区域生长决定归属。双阈值法适合变化区域边界破碎的场景。实际项目中我一般先看强度图的直方图形态再选策略直方图如果呈现明显双峰直接OTSU如果是长尾分布先用OTSU再用均值加标准差做对比两套结果差异大的话要检查输入数据质量。另外要注意阈值分割前先做轻微的高斯平滑避免单像素噪声造成椒盐效果。4.2 OTSU与自适应阈值代码与适用场景from skimage.filters import threshold_otsu def threshold_with_otsu(magnitude): # 直接OTSU thresh threshold_otsu(magnitude) mask (magnitude thresh).astype(np.uint8) return mask, thresh def threshold_adaptive(magnitude, sigma_scale1.0): # 均值 sigma倍数 mean magnitude.mean() std magnitude.std() thresh mean sigma_scale * std mask (magnitude thresh).astype(np.uint8) return mask, thresh # 使用示例 magnitude np.load(magnitude.npy) mask_otsu, th_otsu threshold_with_otsu(magnitude) mask_adapt, th_adapt threshold_adaptive(magnitude, 1.0)逻辑说明threshold_otsu直接调skimage的OTSU实现返回二值掩膜和阈值threshold_adaptive实现了均值加标准差策略sigma_scale控制敏感度。参数说明sigma_scale取0.8到1.2之间比较合理小于0.8会把大量背景像素误检为变化大于1.5则严重漏检。对CVA强度图我常用的值是1.0对IR-MAD强度图因为分布更集中sigma_scale取0.6到0.8更好。4.3 分割后处理面积滤波与连通域分析阈值分割完成后不能直接交出去一个典型的中间结果包含大量独立的小图斑——它们可能是传感器噪声或单像素的配准误差。工程上必须做面积滤波。做法是对二值图做连通域标记统计每个连通域的像素数把小于min_area的连通域整个置为0。这个操作对减少虚警非常有效。from scipy.ndimage import label, sum as ndi_sum def area_filter(mask, min_area100): labeled, num_features label(mask) sizes ndi_sum(mask, labeled, range(num_features 1)) # sizes[0]是背景去掉 keep_labels np.where(sizes min_area)[0] filtered np.isin(labeled, keep_labels).astype(np.uint8) return filtered mask_filtered area_filter(mask_otsu, min_area50)逻辑说明label函数给每个独立的连通区域编号ndi_sum统计每个编号的像素数最后用np.isin把面积小于阈值的区域全体置0。min_area取值取决于影像分辨率0.5米分辨率下50像素大约对应12.5平方米可以作为最小图斑30米分辨率下建议min_area不小于10。参数说明面积滤波只处理“独立的碎图斑”如果变化的楼栋和道路是连通的大块不会被误删。更精细的做法是结合形状特征——比如用长宽比滤除道路用复杂度滤除阴影但常规项目做到面积滤波就够。5. 变化检测的五大常见坑与排查手记5.1 两期影像分辨率不一致导致大量伪变化现象变化检测结果出现大片沿道路和地物边界的线条状伪变化。原因时相1分辨率10米时相2分辨率5米两者没有重采样到同一像元大小直接做逐像素比较导致地物边界错位。这是我见过最频繁的翻车现场。解决先检查两期影像的分辨率、投影和像元偏移务必重采样到相同空间分辨率后对齐。可以用gdal.Warp统一投影再用gdal.ReprojectImage做像元对齐。分辨率不一致时优先采用较高分辨率的像元大小变化检测结果会更精细。配准误差本身也是伪变化的来源如果精度要求高还要做亚像元配准。5.2 辐射差异压过真实变化CVA亮度和光谱都超了现象整个画面几乎全被判定为“变化”尤其是大面积裸地和山体阴影区。原因两期影像成像季节不同或大气条件差异大导致无变化区域的辐射值系统性偏移。CVA和PCA直接比较波段值对此非常敏感。解决先做辐射归一化。IR-MAD本身就是最好的归一化工具——先用5.3节的IR-MAD计算权重提取不变点做线性回归校正再做CVA。这个流程用了很多次比任何直方图匹配都稳定。另外要检查输入的波段是否已经做过大气校正工程上至少要做到TOA反射率否则辐射偏移很难通过算法修正。5.3 IR-MAD迭代权重全部趋近于零结果全黑现象IR-MAD运行正常但权重图全黑变化强度图没有有效信息。原因典型相关分析中某个特征值计算为负数或零使得rho取到0b矩阵出现无穷大MAD分量方差突变。通常是波段数多于样本数或协方差矩阵奇异。另一个常见原因是在带权重的协方差计算中某些波段的标准差为0比如海岸带影像中全水体的波段。解决第一步检查数据中是否有恒定值波段直接剔除第二步在计算b时对rho加一个极小量epsilon1e-8避免除零第三步在权重更新时加上下限约束不要把任何像素的权重压到完全为0。这些修复在代码层面都是两三行的事情但能救回整个流程。5.4 子像素级变化被漏检阈值定得太高现象真实变化区域是细小地物比如单棵树木、临时建筑变化面积占全图不到1%阈值分割后全部被滤掉。原因均值加标准差和OTSU都服务于“大变化量”目标。OTSU在类不平衡严重变化像素占比极低时求出的阈值往往偏高。解决处理细碎变化时改用固定阈值策略先用目视选几个已知变化区统计它们的变化强度值取最小值作为阈值下限。或者用双阈值法低阈值提到全图前95%分位数再由人工检查高阈值与低阈值之间的模糊区域。5.5 云和阴影没有预处理产生大块伪异常现象某期影像上有一块云影检测结果中这块区域全部变成“变化”。原因云和阴影的光谱特征和真实地物变化非常相似尤其在可见光波段。没有任何算法能只靠统计模型区分云和真实变化。解决在变化检测之前必须先做云掩膜把云和云影覆盖的像素排除在统计计算之外。可以从公开的云掩膜产品如Sentinel-2的SCL查得也可以通过光谱阈值自己生成。IR-MAD的鲁棒加权能降低云的影响但不会完全消除。如果两期都有云先做云掩膜再插值填补或者把云区像素直接标记为“未知”不要强行给出变化判断。6. 进阶玩法IR-MAD辐射归一化与人工复核配合IR-MAD最大的实战价值其实不是变化检测而是辐射归一化。前面提到过用权重0.9的像素作为不变点做线性回归这里把具体做法展开。假设时相1是基准时相2需要校正到基准的辐射水平用IR-MAD筛选的不变点分别拟合每个波段的线性关系y_baseline a * y_target b用numpy的polyfit就能做。关键是拟合后要做验证把拟合参数应用到整幅时相2影像后再跑一次IR-MAD或直接算CVA看整体变化强度是否显著下降。如果下降不明显说明不变点选取有问题需要降低权重阈值重新选点。另一个进阶应用是变化结果的人工复核策略。算法输出的变化图任何时候都只能当作“预筛选结果”。在项目交付阶段我会把变化检测结果叠加到两期影像上按空间位置抽检每5000个检测图斑抽查10%统计检测准确率和漏检率。这套策略的价值在于——算法是黑匣子但复核流程能把黑匣子的风险控制住。最后说一个排查技巧当算法输出结果和目视判断严重不符时不要急着调参数先用三维散点图看两期影像在关键波段的联合分布。如果分布呈明显的非线性形状说明辐射关系不是线性的任何基于协方差的方法都受限。这时候需要先做更严格的大气校正或者改用相对辐射归一化方法再回到经典算法上来。这么些年下来我深刻体会到经典算法不是终点但永远是第一步最稳的起点。希望帮到你。本文还有配套的精品资源点击获取