遥感影像变化检测四算法与Matlab代码实践

简介遥感影像变化检测经典算法IR-MAD、MAD、CVA、PCA的Matlab实现代码包面向遥感与地信领域科研人员及学生用于多时相影像中地表变化的自动识别与分析尤其适合光照、大气影响或多传感器数据场景。压缩包共97个文件、约10.4MB以m源码脚本为主附带tif/bmp样例影像、hdr头文件、fig图形输出和mat数据文件可支撑从数据读取、算法运行到结果可视化的完整流程。已有2008人学习下载。代码包含四种经典算法的实现与演示如IRMAD_Update、MADGet、CVADemo、PCA以及envi读写、协方差计算等核心工具函数。同时提供CVA/PCA在泰州影像上的强度图、二值图与变化对比图用户可直接运行New_main.m等主程序并借助KmeansMap、Change Result Comparing等脚本做进一步比较分析便于理解算法原理、复现实验结果。1. 遥感影像变化检测四种经典算法与一套可直接跑的 Matlab 代码做双时刻变化检测时很多人第一步想到的就是两期影像直接相减。我最早在泰州一组多光谱影像上这么干结果变化图里全是条纹噪声光照稍微变了点整个图就像洒了墨水。后来换用 IR-MAD、MAD、CVA、PCA 这套经典算法组合同样的数据变化区域清晰到可以直接出报告。这套资源是一个包含完整 Matlab 实现的代码包里面既有多时相统计方法也有光谱向量差分和主成分分析路线适合刚接触变化检测的研究生、做遥感应用落地的工程师以及被预处理误差折磨到怀疑人生的同行。本文不重复教科书公式推导直接从代码落地讲清楚每个算法的输入输出、参数含义和最容易翻车的地方。2. 原理与选型四种算法的数学逻辑与适用边界2.1 MAD多时相统计异常的第一步MADMultitemporal Anomaly Detection这个名字容易让人误以为是找时间序列里的突变点但在双时刻变化检测里它实际做的是把两期影像看成两组多波段随机变量通过典型相关分析找出一组线性变换让变换后的两期分量差值方差最小化。这组差值就是 MAD 变量理论上包含的变化信息最集中而噪声被压到最低。代码包里MADGet.m做的就是这件事读入两期影像后先计算协方差矩阵和互协方差矩阵再解一个广义特征值问题。核心数学过程可以理解为寻找投影向量 a 和 b使得% MAD核心变换的简化表达实际实现见MADGet.m % X1、X2 分别为两期影像的像素矩阵行像素数列波段数 [~, D] eig(inv(S22) * S21 * inv(S11) * S12); % D的对角线元素对应典型相关系数投影后做差得到MAD变量这里的S11、S22是两期影像各自的协方差矩阵S12、S21是互协方差矩阵。MAD 变量的方差越小说明两期影像在该方向上的统计相关性越强一旦该方向出现显著差异基本可以断定是真实地表变化而不是随机噪声。实际使用时要留意MAD 假设两期影像的像元一一对应所以配准精度直接决定结果质量。泰州数据是已经做过几何配准的拿到其他数据时一定要先检查同名点误差我一般要求 RMS 控制在 0.5 个像元以内。2.2 IR-MAD迭代加权让变化检测更稳IR-MAD 的全称是 Iteratively Reweighted Multivariate Alteration Detection是 MAD 的迭代改进版。基础的 MAD 只跑一轮遇到两期影像存在部分区域明显变化时这些变化像元会干扰协方差矩阵的估计导致变换结果被污染。IR-MAD 的做法是先跑一轮 MAD算出每个像元的变化概率然后给大概率不变的像元更高的权重用加权后的样本重新估计协方差矩阵再更新 MAD 变换迭代若干轮直到权重变化收敛。代码包里IRMAD_Update.m是核心迭代函数LastWeight.mat保存的就是最后一次迭代的权重矩阵。迭代过程中变化概率通过卡方分布计算MAD 变量经过标准化后其平方和近似服从卡方分布自由度等于波段数。利用这个分布可以算每个像元的p-valuep 值越小说明该像元是变化的可能性越大。% 从IRMAD_Update.m中提取的关键步骤用于理解迭代逻辑 % chi2pdf用于计算卡方概率密度判断像元变化概率 p 1 - chi2cdf(sum(mad_squared, 2), n_bands); weight 1 - p; % 变化概率高的像元权重小反之权重大 % 用weight重新计算协方差进入下一轮迭代这里有个重要参数迭代次数。IRMAD_Update.m里通常设置最大迭代次数我一般先设 20 轮并观察权重变化的绝对值。如果连续两轮的权重变化小于 1e-5 就提前终止可以节省大量计算时间。IR-MAD 对光照差异和大气影响的鲁棒性在四种算法里最强不同传感器获取的两期影像也能处理这是它最大的价值。2.3 CVA光谱空间里算出变化方向和强度CVAChange Vector Analysis的思路最简单直观。把每个像元在两期影像上的光谱特征看成多维空间里的一个向量变化检测就是计算这两个向量的差。差向量的模长代表变化强度方向代表变化类型。% CVADemo.m 中变化向量的核心计算 % pixel1、pixel2 都是维度等于波段数的行向量 change_vector pixel2 - pixel1; change_magnitude norm(change_vector); % 变化强度 change_direction change_vector / change_magnitude; % 变化方向模长超过阈值就判定为变化方向则可以用来区分不同的变化类型比如植被变成裸土和水体变浑浊方向在光谱空间里指向不同的区域。这对土地利用变化分析特别有用CVAImageTaizhou_IntensityImage和CVAImageTaizhou_BinaryValue就是强度和二值化的结果。CVA 的优点是结果直观、计算量小但硬伤是极度依赖辐射校正。两期影像的大气条件稍微不同光谱向量的模长就会产生系统性偏差导致大量伪变化。用 CVA 前一定要做大气校正或相对辐射归一化否则就是个华丽的噪声探测器。2.4 PCA在降维空间里浓缩变化信息PCAPrincipal Component Analysis在变化检测里有两种常见用法。第一种是把两期影像的所有波段堆叠成一个高维特征空间然后做主成分分析期望把变化信息集中到前几个主成分上第二种是对两期影像的差值影像做 PCA直接提取变化的主要模式。代码包里PCA.m和eigen2.m负责特征值分解PCADemo.m展示完整流程。% PCA.m 中的核心计算 % stacked_image 是把两期影像波段堆叠后的矩阵行像素数列波段数 cov_matrix cov(stacked_image); [eigenvectors, eigenvalues] eig(cov_matrix); % 按特征值降序排列取前几个主成分用于后续分析 [~, idx] sort(diag(eigenvalues), descend); principal_components stacked_image * eigenvectors(:, idx(1:k));PCA 的优势在于去相关和降噪。多光谱波段之间存在共线性直接逐波段做差值会把噪声重复放大而 PCA 将信息浓缩到少数几个主成分信噪比更高。PCAImageTaizhou_IntensityImage和PCAImageTaizhou_BinaryValue分别对应主成分强度图和二值变化图。这里的 k 值选取是个技术活。经验做法是看特征值贡献率曲线取累计贡献率达到 90% 以上的主成分数。另外提醒一句PCA 在图像处理里应用非常广除了变化检测人脸识别里的特征脸本质也是 PCA所以在不同资料里看到 PCA 出现不必意外核心思想是一样的。2.5 选型判断与场景匹配四种算法没有绝对的优劣只有适不适合。我用一张表给出选型建议算法适用场景主要限制本代码包对应文件MAD配准质量高、双时相统计特征稳定的影像受变化像元干扰一次变换可能不稳健MADGet.mIR-MAD光照/大气差异大、多传感器数据迭代计算量大需要调收敛条件IRMAD_Update.mCVA对变化类型分类有需求、光谱信息丰富的场景对辐射校正极其敏感CVADemo.mPCA数据量大、需要快速压缩变化信息主成分物理含义不明确阈值难解释PCADemo.m如果你手头只有两期影像不确定预处理能做到什么程度我的习惯是先跑 IR-MAD它抗干扰能力最强。如果后续要做变化类型细分再叠加 CVA 的方向信息。这块选型判断也是近几年遥感变化检测研究里一直讨论的问题即使现在出现了开放词汇变化检测这类新方向经典算法的地位依然稳固。3. 代码包落地复现从 ENVI 数据到变化图3.1 文件结构先读懂这套代码怎么组织打开压缩包后看到十多个.m文件和若干.hdr、.mat、无后缀数据文件第一反应可能有点懵。把文件按功能分组后就清楚了。文件/目录作用New_main.m整合主流程按顺序调用各算法freadenvi.m/enviwrite.m/enviwrite2.mENVI 标准格式读写CVADemo.mCVA 算法完整流程入口IRMAD_Update.mIR-MAD 迭代更新核心MADGet.mMAD 变换核心计算PCADemo.m/PCA.mPCA 算法流程入口及分解实现eigen2.m/covw.m/DAcom.m特征分解、加权协方差、典型相关辅助函数KmeansMap.mK-means 聚类后处理把连续变化强度变成类别Change Result Comparing.m多种算法结果对比出图createfigure.m绘图样式辅助data/泰州影像原始数据及中间结果*Taizhou*.hdr各类结果图的 ENVI 头文件主流程是读取两期 ENVI 格式影像 → 计算变化强度 → 阈值化或聚类得到二值变化图 → 结果对比和出图。理解这个流程后再逐个脚本调试就不容易迷路。3.2 数据读取先用 freadenvi.m 打通数据入口代码包处理的泰州影像用的是 ENVI 标准格式即一个二进制数据文件加一个.hdr文本头文件。freadenvi.m负责读数据enviwrite.m和enviwrite2.m负责写结果。先把读取这块跑通后面所有算法才能正常输入。% 读取ENVI影像的标准调用方式 % 参数依次为文件名不含扩展名、文件路径 [data, info] freadenvi(CVAImageTaizhou_IntensityImage, pwd); % data是三维数组行×列×波段 % info中包含了samples、lines、bands、data type、interleave等关键元数据freadenvi.m会自动读取同名.hdr文件根据 header 里的samples、lines、bands和data type来决定二进制文件的读取方式。这里要注意data type必须与实际二进制内容匹配泰州数据一般是 uint8 或 uint16如果 header 写的是 int16 而实际是 uint16读出来的矩阵会整体偏移后续所有计算全部失真。enviwrite.m写文件的逻辑类似它会在写出二进制文件的同时生成配套.hdr。我一般建议在输出文件名里带上算法标识比如IRMAD_result_2024.hdr方便后续多个实验对比。3.3 三个演示入口看懂 CVADemo、IRMAD_Update 与 PCADemo代码包里每个算法都有独立演示脚本。先说CVADemo.m它最直观读取两期影像的光谱强度值逐像元做向量差得到强度图再根据阈值生成二值图CVAImageTaizhou_BinaryValue。% CVADemo.m 核心流程片段 img1 freadenvi(image_t1, pwd); img2 freadenvi(image_t2, pwd); % 计算变化向量模长 mag sqrt(sum((double(img2) - double(img1)).^2, 3)); % 设定阈值生成二值变化图 threshold 50; % 根据影像灰度范围调整 binary_map mag threshold; % 写结果 enviwrite(binary_map, CVAImageTaizhou_BinaryValue, pwd);这段代码里的threshold是 CVA 唯一的参数但它最敏感。泰州影像的城市区域和农田区域在光谱变化幅度上差了一个量级固定阈值必然顾此失彼。通常做法是先看变化强度图的直方图取双峰分布的谷底作为初始阈值。如果谷底不明显就用 Otsu 方法来算。再看IRMAD_Update.m的调用逻辑。IR-MAD 需要两期影像的原始多波段数据它的输入不是单一强度图而是多个波段的完整矩阵。代码里会先初始化所有像元权重为 1然后迭代更新。% New_main.m 中调用IRMAD相关函数的示意 % 初始化权重全1矩阵大小等于影像像元数 w ones(num_pixels, 1); for iter 1:max_iter % 用当前权重计算加权协方差矩阵 cov_w covw(img1_t, img2_t, w); % DAcom.m实现典型相关分析返回变换向量 [a, b] DAcom(cov_w); % 计算MAD变量 mad_variate img1_t * a - img2_t * b; % 更新权重并检查收敛条件 w_new update_weights(mad_variate); if max(abs(w_new - w)) tol break; end w w_new; end这段代码展示了 IR-MAD 的迭代逻辑。covw.m接收三个参数前两个是两期影像的像元值矩阵第三个是权重向量返回加权协方差矩阵。DAcom.m做典型相关分解相当于对加权后的两组变量求投影方向。update_weights函数就是前面提到的卡方分布概率计算代码包直接集成在IRMAD_Update.m里。PCADemo.m的路数不一样。它先把若干波段的影像向量化堆叠再做协方差分解得到的PCAImageTaizhou_BinaryValue不是某种地表覆盖的直观结果而是数据方差最大的方向。这一步需要先想清楚PCA 输出的前几个主成分里到底哪个主成分最能代表变化我一般会先画出前三个主成分的灰度图用眼睛确认变化主要集中在哪个通道再接阈值化。3.4 后处理KmeansMap 与 Change Result Comparing 让结果可读原始算法输出一般是强度图直接看也能看出大致变化区域但要出报告或者做定量分析还需要把连续值变成离散类别。KmeansMap.m做的事情就是把变化强度图聚类成几类每一类对应不同的变化程度。% KmeansMap.m 的典型调用方式 % 输入PCAImageTaizhou_IntensityImage连续强度图 % 输出聚类后的类别图每个像素值为1~K之间的整数 K 4; % 聚类类别数 [cluster_map, cluster_centers] kmeans(mag_vector, K, Distance, sqEuclidean);这里的 K 值选多少直接决定结果图的信息粒度。K2 时就是单纯的变化/不变二分类K4 时可以把变化细分成微变、轻度、中度、剧烈四档但 K 太大也会把噪声单独聚一类。我的建议是结合地表实际情况选 K比如城市扩张场景下K3 到 K5 比较合理。Change Result Comparing.m是一个集成脚本它把 MAD、IR-MAD、CVA、PCA 四种算法的结果并排显示方便对比。跑完这个脚本你能直观看到每个算法在同一数据集上的表现差异用于算法选型和写论文的对比图都很方便。4. 避坑笔记参数设置与环境依赖的五个翻车点4.1 没做辐射归一化CVA 输出全场飘红现象CVADemo 跑完变化强度图几乎全图高亮看不出具体变化区域二值化后噪声点密集成片。原因泰州影像虽然经过了基础预处理但两期影像之间的大气条件差异仍然存在直接逐像元做光谱向量差所有像元都会叠加一个系统偏差。CVA 对辐射校正的敏感度在四种算法里最致命。解决跑 CVA 之前先做相对辐射归一化。常见做法是取影像中的不变目标如深水、裸岩、大型建筑屋顶作为伪不变特征点用线性回归把两期影像的灰度值拉到同一个水平。代码包里虽然没有提供专门的归一化脚本但可以先用enviread把两期影像读进来手动选几个不变区域计算增益和偏移后校正第二期影像再跑 CVADemo。从那以后我每次拿到新数据第一步都是先看两期影像的同名地物光谱值差多少再决定用什么算法。4.2 IR-MAD 迭代不收敛权重矩阵发散现象IRMAD_Update.m 跑了很久LastWeight.mat里的权重值出现大量负数或极端值变化概率图看起来像随机噪声。原因一个典型原因是两期影像的辐射差异太大初始 MAD 变换找不准投影方向迭代过程中权重计算不稳定另一个原因是波段数太少比如只有 3 个波段卡方分布的自由度太低概率估计方差大。解决第一对输入影像做归一化把所有波段缩放到 0-1 范围提高数值稳定性第二在权重更新时加入平滑约束比如用 3×3 窗口对权重图做均值滤波抑制逐像元的剧烈跳变第三适当增大收敛阈值比如tol从 1e-5 改成 1e-4避免陷入无意义的精细迭代。实测泰州数据上把迭代上限设 20、收敛阈值设 1e-4 时结果最稳定。4.3 ENVI 头文件 data type 与二进制实际内容不匹配现象freadenvi 读出来的矩阵全是一个值或者数值范围明显不对比如影像本该是 0-255 的灰度值读出来却是几万。原因.hdr文件里的data type字段写错了。ENVI 格式的 data type 有严格定义1byte2integer12unsigned integer4float等等。头文件与实际数据存储不一致时Matlab 按头文件解析二进制自然全盘出错。解决用十六进制编辑器打开原始数据文件看前几十个字节的数值范围判断实际数据类型。如果数据文件里出现大于 255 的值基本可以断定是 uint16 或 int16。改好.hdr里的data type字段后重新读取即可。freadenvi.m本身只是个解析器不要指望它能自动纠错。4.4 两期影像行列数不一致导致矩阵维度报错现象MADGet.m 或 IRMAD_Update.m 运行到协方差矩阵计算那一步报错提示矩阵维度不一致。原因两期影像的采样范围或者裁剪范围没有对齐比如一个 512×512另一个 520×520。典型相关分析要求两个输入矩阵的行数像元数严格相等。解决读取影像后立即检查size(img1)和size(img2)不一致时以较小的一期影像为基准裁剪。我习惯写一个对齐函数% 影像对齐以较小尺寸为基准裁剪 [rows, cols, bands] size(img1); img2 img2(1:rows, 1:cols, :); % 并更新对应的坐标信息避免几何偏移注意这里的裁剪不是单纯的矩阵截取如果两期影像的投影坐标不完全重合必须先用地理坐标做重采样对齐否则只是数字上对上了空间位置上仍然错位变化检测结果没有地理意义。4.5 K-means 聚类数 K 是拍脑袋定的现象KmeansMap 聚类结果中某一类占据了 95% 的像素其余类别几乎看不见或者聚类结果在多次运本文还有配套的精品资源点击获取