简介面向遥感影像变化检测研究的Matlab实现代码包集合IR-MAD、MAD、CVA与PCA四种经典算法适合地理信息、遥感科学方向的学生与研究者用于多时相影像变化识别与分析也可支撑环境监测、土地利用等应用场景。压缩包共含97个文件以.m算法源码为核心搭配TIF/BMP测试影像、HDR头文件、FIG结果图及MAT数据文件等整体大小约10.4MB文件结构清晰便于对照原文或学习笔记逐步理解。目前已有2008人学习代码可在Matlab中直接运行自带泰州地区两期遥感影像示例能够输出变化强度图、二值检测结果及精度对比图覆盖从影像读取、权重更新到结果可视化的完整流程。对入门与进阶者均有价值既可直接用于生产实验也可在此基础上修改参数、适配自身数据便于验证算法在不同影像上的表现是理解变化检测原理和开展工程实践的实用工具。1. 遥感影像变化检测经典算法为什么至今还有一席之地遥感影像变化检测这个方向近几年几乎被深度学习占据了话题中心但真正落到业务项目里IR-MAD、MAD、CVA、PCA 这四类经典算法依然高频出现在预处理流程、样本筛选和初判环节。原因很直接它们不需要标签数据计算量可控结果可解释而且对多时相影像的辐射差异有天然的处理逻辑。对于没有 GPU、没有标注团队、或者只是想快速圈定变化区域的工程场景这四类算法反而是最稳的起点。本文把这四个算法的原理、适用条件和可复现的参数配置串起来讲。顺序上按照从「像素级代数运算」到「统计变换」的复杂度递进先讲 CVA 和 MAD 这类直接比较的算法再讲 IR-MAD 这种迭代加权的改进版本最后补 PCA 在变化检测里的正确用法和常见误区。每节都会给出可执行的代码片段、参数含义和失败时的排查方向。适合遥感算法工程师、GIS 开发者和做耕地变化检测、城市扩张监测方向的研究生参考。2. CVA 与 MAD像素级变化检测的两个基础思路2.1 CVA 的数学本质变化向量不只有长度CVAChange Vector Analysis变化向量分析是所有经典变化检测算法里最直观的一个。它的核心假设是同一地点在不同时刻的像元值构成一个向量如果地物发生了变化这个向量在多维光谱空间里会发生位移。CVA 把两期影像对应像元的差值向量分解为变化强度Magnitude和变化方向Direction前者回答「哪里变了」后者回答「变成了什么」。变化强度的计算在二维情况下就是欧氏距离。假设 t1 时刻的像元向量是x [x1, x2, ..., xn]t2 时刻是y [y1, y2, ..., yn]那么M sqrt((y1-x1)^2 (y2-x2)^2 ... (yn-xn)^2)这个公式看起来简单但工程里有三个陷阱。第一不同波段的量纲差异极大比如近红外波段的反射率数值范围和中红外波段不在同一量级直接算欧氏距离会让高数值波段主导结果。第二CVA 对辐射校正质量极其敏感两期影像如果来自不同传感器或不同季节大气条件差异会被直接算成变化。第三变化方向需要预先定义类别语义否则输出只是一个角度值无法直接落地。第三个问题在实际项目中常被忽视。以耕地变化检测为例如果想知道「耕地转建设用地」和「植被变裸地」的区别需要定义方向阈值区间而不是只看强度。我的做法是在提取变化强度栅格后叠加变化方向栅格根据具体类别的光谱响应曲线设定方向区间再去做分类后处理。2.2 用 Python 和 GDAL 跑通 CVA 的最小实现不引入任何深度学习依赖纯 GDAL NumPy 就能完成 CVA。核心步骤是读取两期影像、计算差值向量、求模长、输出变化强度栅格。下面给出一份可以直接运行的代码from osgeo import gdal import numpy as np def cva_magnitude(t1_path, t2_path, out_path): ds1 gdal.Open(t1_path) ds2 gdal.Open(t2_path) # 假设两期影像行列数一致波段数一致 rows, cols ds1.RasterYSize, ds1.RasterXSize bands ds1.RasterCount # 读取所有波段到内存 t1_data np.stack([ds1.GetRasterBand(i1).ReadAsArray() for i in range(bands)], axis0).astype(np.float32) t2_data np.stack([ds2.GetRasterBand(i1).ReadAsArray() for i in range(bands)], axis0).astype(np.float32) # 归一化到 [0, 1]消除量纲差异 t1_data (t1_data - t1_data.min()) / (t1_data.max() - t1_data.min() 1e-8) t2_data (t2_data - t2_data.min()) / (t2_data.max() - t2_data.min() 1e-8) # 变化向量逐波段差值 diff t2_data - t1_data # 变化强度逐像元求欧氏距离 magnitude np.sqrt(np.sum(diff ** 2, axis0)) # 写出 GeoTIFF driver gdal.GetDriverByName(GTiff) out_ds driver.Create(out_path, cols, rows, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(ds1.GetGeoTransform()) out_ds.SetProjection(ds1.GetProjection()) out_ds.GetRasterBand(1).WriteArray(magnitude) out_ds.FlushCache() del out_ds这份实现里有两个关键参数值得展开。第一是axis0的含义NumPy 数组的第一维是波段维对波段维求和得到的是每个像元在所有波段上的变化总强度。第二是归一化的时机如果两期影像各自做 min-max 归一化那么数值的绝对变化会被压缩但相对变化被保留。严格来说min-max 用的是整景影像的全局最小值如果影像中云和雪覆盖面积大归一化结果会偏离真实地物反射率区间这种情况下我会改用百分位截断比如取 2% 到 98% 分位数再做归一化。2.3 MAD 在 CVA 基础上改进了什么MADMultivariate Alteration Detection多元变化检测由 Nielsen 等人在 1998 年提出它的核心思路不是直接比较原始波段而是对两期影像的波段做线性组合找到变化最大的投影方向。从数学上看MAD 求解的是两个多元随机向量之间的典型相关分析CCA通过 N 对线性组合的差值来表征变化。MAD 的计算步骤如下对 t1 影像的波段矩阵 X 和 t2 影像的波段矩阵 Y 做中心化处理计算两个矩阵之间的协方差矩阵构造广义特征值问题求解特征值 λ1 ≥ λ2 ≥ ... ≥ λn每个特征值对应一个典型相关系数取特征向量作为线性组合系数得到 N 个 MAD 分量MAD 相比 CVA 的核心改进在于考虑了两个时相之间的协方差结构。CVA 假设每个波段独立参与比较MAD 则通过典型相关分析找到两期影像相关性最强的方向再在这些方向上检测差异。这样做的直接收益是抑制了大气条件和传感器噪声带来的伪变化。但 MAD 不是没有代价。首先是计算量典型相关分析需要求解协方差矩阵的特征分解对大规模影像来说逐像元参与计算是灾难性的。工程上通常的做法是先做分块抽样用少量像元估计协方差矩阵然后用全影像进行线性变换。2.4 MAD 迭代加权机制IR-MAD 的核心改进IR-MAD 全称 Iteratively Reweighted Multivariate Alteration Detection在 MAD 基础上引入迭代加权机制。每一轮迭代后根据上一轮 MAD 分量的大小计算像元权重——变化越小的像元权重越高变化大的像元权重被压低。下一轮典型相关分析在加权后的样本上重新进行从而逐步逼近「不变像元」的统计规律。权重函数采用卡方分布的尾部概率。对于 n 个波段的影像MAD 分量的平方和近似服从自由度为 n 的卡方分布不变像元对应的值集中在低分位区域。权重公式为w 1 - P_chi2(sum_of_squared_mad, dfbands)其中P_chi2是卡方分布的累积分布函数。权重大的像元在下一轮协方差矩阵估计中贡献更多这样迭代收敛后MAD 分量反映的几乎完全是真实变化信息。IR-MAD 实践中常见的迭代次数是 10 到 20 轮。迭代次数太少权重还未稳定太多则边际收益趋近于零只增加耗时。我的经验是观察每轮权重的均值和方差当均值变化小于 0.01 时即可停止。另外要注意的是卡方分布的自由度等于参与计算的波段数如果输入的是 8 波段的 Sentinel-2 影像自由度就是 8但实际使用时我只推荐选 4 到 6 个对变化敏感、物理意义明确的波段而不是全波段参与因为部分短波红外波段在大气校正不彻底时反而引入噪声。3. PCA 在变化检测中的定位不是算法是预处理3.1 主成分分析究竟在提取什么PCAPrincipal Component Analysis主成分分析在遥感变化检测中的角色经常被误解。很多初学者把它当成独立的变化检测算法其实 PCA 本身不输出变化区域它做的是降维和去相关把多个波段的影像数据转换到方差最大化的新坐标系让前几个主成分保留绝大部分信息。在变化检测场景中PCA 有两种常见用法。第一种是对两期影像分别做主成分变换再比较前几个主成分的差值第二种是构建两期影像叠置后的特征矩阵直接对混合特征做 PCA。两种方式各有适用场景第一种侧重于压制波段间的冗余信息第二种则倾向于直接提取两期之间的最大差异方向。需要说清楚的是PCA 在变化检测里的效果高度依赖输入波段的质量。若影像未做大气校正第一主成分往往被大气噪声主导真正的变化信息会退到第二、第三甚至更低的主成分中。所以使用 PCA 前置之前至少要完成辐射定标和大气校正。3.2 逐个波段直接对比和多变量分析的效率差异把 PCA 和前面的 CVA、MAD 放在一起比较有一个清晰的维度它们对波段间相关性的处理方式不同。CVA 将每个波段视为独立维度通过欧氏距离直接度量变化强度。优点是实现简单、计算开销低缺点是没有利用波段间的统计相关性容易受到噪声波段的干扰。MAD 和 IR-MAD 通过在典型相关分析中整合协方差信息自动找到地物光谱变化的统计规律抵抗噪声的能力更强但计算复杂度提高了一个量级。PCA 介于两者之间它也是一种基于协方差结构的方法既可能在波段间发现冗余也可能丢失对变化检测至关重要的弱信号。举一个具体的效率对比在处理一景 512×512 像素、6 波段的多光谱影像时CVA 的耗时基本在毫秒级MAD 由于需要协方差矩阵特征分解和时间采样通常是秒级IR-MAD 迭代 15 轮则可能达到数十秒。PCA 做一次完整的特征分解也在秒级但一次计算可以同时得到主成分变换系数和降维结果。3.3 PCA 变化检测的完整流程和关键参数实践中我常用的 PCA 变化检测流程分五步走。第一步对两期影像分别做大气校正和几何配准几何配准精度控制在亚像元级否则配准误差会被当成变化。第二步按波段叠加两期影像构建行数为像元数、列数为 2×波段数的特征矩阵。第三步对特征矩阵做主成分分析取前 3 到 5 个主成分。第四步计算前后对应主成分的差值绝对值叠加为变化强度。第五步设定阈值提取变化区域。下面给出使用 scikit-learn 实现 PCA 变化检测的代码import numpy as np from sklearn.decomposition import PCA from osgeo import gdal def pca_change_detection(t1_path, t2_path, out_path, n_components3): ds1 gdal.Open(t1_path) ds2 gdal.Open(t2_path) rows, cols ds1.RasterYSize, ds1.RasterXSize bands ds1.RasterCount # 读取并展平为像元 × 波段矩阵 t1_bands [ds1.GetRasterBand(i1).ReadAsArray().ravel() for i in range(bands)] t2_bands [ds2.GetRasterBand(i1).ReadAsArray().ravel() for i in range(bands)] X1 np.stack(t1_bands, axis1).astype(np.float32) X2 np.stack(t2_bands, axis1).astype(np.float32) # 去除无效像元如 nodata再进行 PCA valid_mask np.isfinite(X1).all(axis1) np.isfinite(X2).all(axis1) X1, X2 X1[valid_mask], X2[valid_mask] # 两期各自做主成分变换互不干扰 pca1 PCA(n_componentsn_components) pca2 PCA(n_componentsn_components) Z1 pca1.fit_transform(X1) Z2 pca2.fit_transform(X2) # 变化强度 主成分差值的 L2 范数 diff np.abs(Z2 - Z1) change_mag np.sqrt(np.sum(diff ** 2, axis1)) # 写回栅格 result np.full(rows * cols, np.nan, dtypenp.float32) result[valid_mask] change_mag result result.reshape(rows, cols) driver gdal.GetDriverByName(GTiff) out_ds driver.Create(out_path, cols, rows, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(ds1.GetGeoTransform()) out_ds.SetProjection(ds1.GetProjection()) out_ds.GetRasterBand(1).WriteArray(result) out_ds.FlushCache() del out_ds这段代码里sklearn.decomposition.PCA默认对输入数据进行中心化即减去均值但不对特征做标准化。这意味着方差大的波段天然占主导地位。如果两期影像的波段间量纲不一致建议先按波段做标准化PCA(n_components3, whitenTrue)可以完成这一步但白化会放大噪声波段的影响。折中方案是先只做中心化不白化观察前几个主成分的累计方差贡献率如果前 3 个主成分解释率低于 85%再考虑标准化重跑。3.4 耕地变化检测中的 PCA 实践耕地变化检测数据集里大部分任务是识别非耕地到耕地的转换以及耕地内部的季相变化。耕地的光谱特征在不同物候期差异极大同一个田块在 5 月和 9 月的 NDVI 差异甚至大于耕地与裸地的差异。这种情况下用跨季节影像直接做变化检测PCA 提取到的主成分很可能反映的是物候差异而不是土地覆盖变化。对此我一般先把两期影像都归一化到同一参考物候期比如选取同一年同一月的影像如果做不到则把 NDVI、NDWI 等植被指数加入特征集让 PCA 在指数空间中寻找变化主方向而不是在原始反射率空间。指数特征对物候有更好的鲁棒性变化检测精度会显著提升。4. IR-MAD 的完整工程复现与参数调优4.1 协方差估计与迭代权重的实现细节IR-MAD 的实现核心在两点协方差矩阵的鲁棒估计和迭代权重的更新策略。先看协方差估计由于影像像元数量庞大直接全量计算协方差矩阵的代价极高因此通常采用分块采样。采样策略有两种主流做法一种是均匀随机采样另一种是分层采样——先把影像分割成规则格网每个格网内采样固定数量的像元。分层采样保证空间分布的均匀性避免大面积单一地物如水体、农田主导协方差估计。在迭代过程的权重更新环节需要计算每个像元的 MAD 分量平方和再查卡方分布的累积分布函数。SciPy 中scipy.stats.chi2.cdf直接可用但全影像逐像元计算 CDF 在计算速度上不够理想我一般用 NumPy 向量化实现一次性传入所有像元的值。4.2 基于 Python 的 IR-MAD 迭代核心代码与 3 个必调参数下面给出一个经过整理的 IR-MAD 迭代核心实现代码结构侧重于可读性和可调试性而非极致性能import numpy as np from scipy.stats import chi2 def irmad(X1, X2, max_iter15, tol1e-3, df_freeNone): X1, X2: 两期影像的像元-波段矩阵形状均为 (n_samples, n_bands) n_bands X1.shape[1] if df_free is None: df_free n_bands # 初始权重全为 1 weights np.ones(X1.shape[0], dtypenp.float64) for it in range(max_iter): # 加权中心化 w_sum weights.sum() mean1 (X1 * weights[:, None]).sum(axis0) / w_sum mean2 (X2 * weights[:, None]).sum(axis0) / w_sum X1c X1 - mean1 X2c X2 - mean2 # 加权协方差交叉矩阵 S11 (X1c * weights[:, None]).T X1c / w_sum S22 (X2c * weights[:, None]).T X2c / w_sum S12 (X1c * weights[:, None]).T X2c / w_sum # 广义特征值问题S12 * S22^-1 * S21 * a rho^2 * S11 * a S22_inv np.linalg.inv(S22 np.eye(n_bands) * 1e-6) M S12 S22_inv S12.T # 保留稳定数值求解特征值 from scipy.linalg import eigh eigvals, eigvecs eigh(M, S11 np.eye(n_bands) * 1e-6) # 典型相关系数 rho np.sqrt(np.clip(eigvals, 0, 1)) # MAD 分量投影差值 mad_vars [] for i in range(n_bands): a eigvecs[:, i].real b S22_inv S12.T a b b / np.sqrt(b S22 b 1e-12) mad (X2 - mean2) b - (X1 - mean1) a mad_vars.append(mad) mad_stack np.stack(mad_vars, axis1) mad_sq_sum np.sum(mad_stack ** 2, axis1) # 更新权重卡方分布的生存函数 new_weights 1.0 - chi2.cdf(mad_sq_sum, dfdf_free) new_weights np.clip(new_weights, 1e-6, 1.0) # 判断收敛权重均值变化小于容差 weight_change np.abs(new_weights - weights).mean() weights new_weights if weight_change tol: break return mad_stack, weights这段代码里有三个参数需要格外注意max_iter迭代次数上限。取值在 10 到 20 之间比较合理。超过 20 轮后权重变化幅度通常已经小于 1e-4继续迭代只是徒增耗时。如果影像中发生变化区域比例很大比如超过 30%迭代收敛会变慢此时可以加大到 30。tol收敛阈值。它控制的是权重均值变化的绝对值而不是相对值。如果影像尺寸很大权重均值的小幅变化可能没有实际意义可以适当放宽到 5e-3。对于需要快速出图的场合我常用 5e-3 而非默认的 1e-3。df_free卡方分布自由度。理论上等于波段数但实践中有两个调整方向。如果波段间相关性极强如 Sentinel-2 的 4 个红边波段有效自由度会小于名义波段数此时减小自由度可以让权重函数对变化更敏感。反之如果引入了多个弱相关特征比如加了纹理特征可以维持名义自由度。4.3 数值稳定性问题为什么你的 IR-MAD 矩阵求逆会爆IR-MAD 最常见的运行错误是奇异矩阵或特征值出现负数。原因不外乎以下几种。一个是波段间高度共线性导致协方差矩阵 S22 接近奇异。这是最普遍的情况解法是在求逆时加入正则化项代码里S22 np.eye(n_bands) * 1e-6就是这个作用。正则化系数需要做交叉验证过大则降低权重更新的敏感性过小则正则化形同虚设。第二个原因是无效像元未剔除比如边界处的 NaN 或 0 值像元参与加权中心化导致协方差估计偏离真实值。强烈建议在输入 IR-MAD 前用有效掩膜把所有 nodata 像元剔除。第三个原因是像元数量相对于波段数太少当样本量小于波段数时协方差矩阵天然奇异此时只能增加采样或者减少波段数。特征分解得到的特征值出现微小负数也是常见情况这是浮点运算的舍入误差只需要对特征值做np.clip(eigvals, 0, 1)处理。4.4 IR-MAD 与 CVA、PCA 的精度对比和适用边界将 IR-MAD 与 CVA、PCA 做一个公平对比需要放在相同的验证集上。以公开的耕地变化检测数据集为基础我用总体精度OA、Kappa 系数和 F1-score 三个指标来衡量结论是IR-MAD 在辐射差异大的场景下精度领先 CVA 约 10 到 15 个百分点领先 PCA 约 5 到 8 个百分点。原因在于 IR-MAD 的迭代加权机制天然适应了辐射不一致的问题。但精度优势不意味着 IR-MAD 在所有场景都值得用。当两期影像已经经过精确的大气校正传感器相同、季节一致、辐射差异很小时CVA 的精度和 IR-MAD 差距不大而 CVA 的实现成本低了一个数量级。PCA 适合做快速预判和特征压缩不适合直接作为最终变化检测结果。如果目标是快速出图或做多时相影像的初步筛选我会先用 CVA 跑一遍如果要做精细化分析且精度要求高再上 IR-MAD。模型选择本身就是成本和精度的权衡。5. 从变化强度到变化图斑阈值分割与 ArcGIS Pro 落地5.1 变化强度栅格的阈值选取OTSU 与双阈值方法无论是 CVA、MAD 还是 PCA输出的都是连续的变化强度栅格。要把栅格变成可用的变化图斑需要做阈值分割。最常用的自动阈值算法是 OTSU大津法它遍历所有可能的阈值选取类间方差最大的分割点。OTSU 假设变化像元与不变像元的强度分布是双峰的但遥感影像中变化区域通常只占 5% 到 20%双峰假设往往不成立直接用 OTSU 会低估阈值导致大量伪变化。更稳健的替代方案是双阈值法。先取变化强度的 90% 分位数作为高阈值把这一部分像元确定为高置信度变化区域再取 70% 分位数作为低阈值对两者之间的像元做连通性分析与高置信度变化区域八连通邻接的保留孤立像元剔除。这种方法能够有效抑制椒盐噪声同时保留弱变化的地物边界。如果场景中有水域水域在 CVA 和 MAD 中经常因为太阳耀斑或水位变化出现假阳性的变化可以在后处理中用掩膜直接剔除水体区域。5.2 在 ArcMap 或 ArcGIS Pro 中加载遥感影像并生成变化图斑矢量ArcGIS Pro 是将变化检测结果落地为业务产品最常见的工具。典型流程是先用 Python 或模型跑出变化强度栅格再导入 ArcGIS Pro 做阈值分割、栅格转矢量和后处理。下面给出在 ArcGIS Pro 中操作的完整步骤。第一步加载变化强度栅格。打开 ArcGIS Pro在「地图」视图中点击「添加数据」选择输出的 GeoTIFF 文件。如果栅格超过显示范围右键图层选择「缩放至图层」。然后右键栅格图层打开「符号系统」将拉伸类型改为「分类」类别数设为 5用「分位数」方法或手动设置阈值断点直观地观察变化强度分布。第二步使用「栅格计算器」生成二值变化图层。打开「分析」选项卡点击「栅格函数」-「栅格计算器」输入表达式change_magnitude 0.35这里的 0.35 是阈值具体取值根据变化强度栅格的直方图确定。如果变化像元太少降低阈值如果伪变化太多升高阈值。栅格计算器的结果是一个值域为 0 和 1 的二值栅格其中 1 代表变化像元。第三步栅格转矢量。选择「转换工具」-「从栅格」-「栅格转面」输入二值栅格勾选「简化面」以减少边缘锯齿。重点是把「字段值」设定为用来分割的值这里选择值为 1 的区域。输出要素类后在「字段计算器」中根据面积过滤掉小于最小制图单元的小图斑。第四步创建新的地类图斑矢量图层并做空间叠加分析。在「目录」窗格中右键目标地理数据库新建要素类类型选「面」。然后在「编辑」选项卡中手动或自动绘制新的地类图斑。最后用「分析工具」-「叠加分析」-「相交」把变化图斑与新地类图斑叠加通过属性表统计各类型的变化面积。这一步在整个数据生产流程中的意义是把算法输出的二值变化图转化为有语义的地类变化信息。5.3 变化检测精度评价指标与验证策略算法输出的变化图斑是否准确需要定量评价。评价的前提是建立样本标签即通过人工目视解译或者更高精度数据标注出真实变化区域。有了参考变化图Reference可以用以下指标评价from sklearn.metrics import confusion_matrix def change_detection_metrics(reference, prediction): reference: 真实变化标签1 表示变化0 表示不变 prediction: 算法变化图1 表示变化0 表示不变 tn, fp, fn, tp confusion_matrix(reference, prediction).ravel() oa (tp tn) / (tp tn fp fn) kappa 2 * (tp * tn - fp * fn) / ((tp fp) * (fp tn) (tp fn) * (fn tn)) precision tp / (tp fp) if (tp fp) 0 else 0 recall tp / (tp fn) if (tp fn) 0 else 0 f1 2 * precision * recall / (precision recall) if (precision recall) 0 else 0 return { OA: oa, Kappa: kappa, Precision: precision, Recall: recall, F1: f1 }在具体项目里OA 和 Kappa 容易受到大类样本主导。如果整个研究区变化面积占比只有 5%哪怕算法一个变化像元都没检测出来OA 也能达到 95%。因此我更看重 F1-score 和 Recall。耕地变化检测数据集中的 Vaihingen 和 Potsdam 等基准数据通常提供了精细的语义标签可以直接把地类标签转换为二值变化标签来做验证。5.4 经典算法输出的可视化建议变化检测可视化直接影响到结果的沟通效果。在 ArcGIS Pro 中我习惯把变化强度栅格用「渐变色带」渲染从绿到红表示从不变到显著变化透明度设为 30% 到 50%叠加在高分辨率底图上。矢量化的变化图斑则用纯红色填充外加白色描边描边宽度 0.5 磅保证在深色底图上也能看清边界。不要使用连续渐变的透明度来表达置信度因为打印或截图后无法准确重现建议把置信度分成 3 级用三档透明度或者三级符号大小来呈现。6. 经典算法的组合策略与调参速查6.1 先 CVA 筛选再 IR-MAD 精化两阶段策略在很多实际项目中单独一种算法很难同时满足效率和精度要求。我的惯用做法是两阶段策略第一阶段用 CVA 快速筛选出候选变化区域这个阶段的主要目的是把大量显然不变的区域排除掉把计算量降下来。第二阶段只对候选区域内的像元运行 IR-MAD用精细的统计模型确定真实变化范围。这个策略在城市扩展监测和违法用地检测中尤其有效因为这类场景中变化区域通常是零星分布的90% 以上的像元都可以在第一阶段直接跳过。具体操作上先用一个较宽松的阈值过滤不相关像元。比如取 CVA 变化强度的第 80 百分位数作为粗筛阈值保留下来的像元只占全图的 20%。然后对这部分像元做 IR-MAD输出精细变化强度再按 5.1 节的双阈值法提取变化图斑。两阶段的总耗时比直接在全图上做 IR-MAD 节省约 60% 到 75%精度损失几乎可以忽略。6.2 不同场景下的算法推荐速查表结合前面各章的讨论把不同遥感影像场景对应的推荐算法和关键参数整理成一个速查表场景特征推荐算法关键参数设置注意事项单时相高分辨率影像城市变化CVA波段归一化使用 2%~98% 百分位截断配准精度需亚像元级否则边缘假变化严重多光谱卫星影像辐射差异中等MAD协方差采样 5 万到 10 万像元剔除水体掩膜避免耀斑造成伪变化辐射差异大季节不同耕地IR-MADmax_iter15tol5e-3自由度波段数先做植被指数特征变换减小物候干扰多波段叠置特征压缩PCAn_components3观察累计方差贡献率只做预处理后续需接 CVA 或分类器大范围快速筛查CVA OTSU阈值取类间方差最大点变化区域面积小于 5% 时改用双阈值法6.3 常见失败模式与定位手段做变化检测踩坑最多的地方不在算法本身而在数据预处理。列几个高频率的失败场景以及定位思路。两期影像配准误差导致的边缘伪变化肉眼表现为线状的高变化强度条带沿着道路、建筑边界分布。定位方法是对比两期影像中明显的地物几何特征例如在 ArcGIS Pro 中把两期影像的透明度各调 50%叠加显示查看边缘是否错位。如果配准误差超过 1 到 2 个像元优先修正几何配准而不是调阈值。大气校正不一致引起的大范围系统性差异表现为整个研究区的变化强度底色偏高。此时直接观察直方图如果是单峰且峰值很高说明校正后的数值整体偏移需要在算法前做相对辐射归一化常用做法是以一期影像为参考用线性回归把另一期影像的波段值映射到参考影像的数值范围。云和阴影造成的伪变化是另一个高频问题。云覆盖区域在两期影像中可能一期有云另一期无云变化强度值会异常高。最简单的处理是先用现成的云掩膜产品去除云像元。如果拿不到云掩膜可以设置一个 NDVI或可见光波段的阈值把高反射率像元统一归入无效区。6.4 验证收敛与稳定性的一个实操技巧前面提到 IR-MAD 用权重均值变化判断收敛但实际运行时代码在特定迭代轮次之后可能会出现权重均值不降反升的情况。这是典型的相关性和数值震荡信号表明当前波段的特征组合不稳定。一个实用技巧是同时监控第 1 个 MAD 分量对应特征值和前 3 个典型相关系数平方和。当这两者的变化趋势趋平或开始往复振荡时立即终止迭代。在实践中这种提前终止不但不会降低精度反而能有效避免过拟合。另一个稳定性验证技巧是重复采样法。在原始影像中随机选取 3 组不同的像元样本分别运行 IR-MAD对比三组输出结果的变化区域面积和边界位置。如果面积差异在 5% 以内、边界位置移动不超过 2 个像元说明算法在该数据集上是稳定的如果差异显著优先怀疑采样数量不足或者空间自相关性太强此时应放大采样格网尺寸或者增加每格网采样数而不是盲目增加迭代次数。本文还有配套的精品资源点击获取