尧图网站建设 
简介压缩包内汇总了遥感影像变化检测中四种经典算法的 Matlab 实现面向遥感、测绘及地理信息专业的科研人员和学生可用于多时相地表变化分析。包内共有97个文件包括.m算法脚本、ENVI标准.hdr/.tif影像数据、用于结果展示的.fig图表以及若干.bmp辅助图像压缩后约10.4MB文件组织清晰便于对照学习。其中包含IR-MAD、MAD、CVA、PCA各自的Demo与核心函数并附泰州区域的真实影像及中间结果读者可直接运行观察强度图、二值图与变化检测对比效果。此外还包括enviwrite.m、covw.m等读写与统计工具方便将算法迁移到自己的数据集。目前已有约2008人学习下载适合想通过现成代码理解算法原理、开展实验验证的入门与进阶使用者。1. 为什么传统遥感变化检测算法仍是工程首选十年前做地物变化分析大家用ENVI里的工具一步一步点现在一提起变化检测很多人第一反应就是上深度学习。但你真正落过地就会发现神经网络需要标签、需要训练数据而且在跨传感器、跨时相的影像上泛化能力并不稳。IR-MAD、MAD、CVA、PCA这四种经典算法胜在不依赖样本、数学可解释、跑起来快。这套MATLAB代码以台州区域的ENVI格式影像为样例完整实现了从影像读取、预处理到变化结果输出的全链路。对刚接触遥感影像预处理的研究生来说是理解统计类变化检测最好的入门素材对已经用深度学习做过语义分割复现的工程师同样适合拿来做候选区筛选和基线对比。2. MAD与IR-MAD从统计异常到迭代加权的演进逻辑2.1 MAD的数学建模线性组合如何分离变化信息MADMultitemporal Anomaly Detection多时相异常检测的数学内核是典型相关分析。它把两期影像看作两个随机向量X和Y试图找到一组线性组合系数a和b使得组合后的变量UaᵀX和VbᵀY之间的相关性最小。当两期影像在同一位置的像元未发生变化时U−V的值应该接近零而真正变化的区域会因为光谱关系被破坏而产生显著的偏离值。这套代码里的MADGet.m实现了完整的求解过程。核心是构造两期影像的联合协方差矩阵将其分块后求解广义特征值问题。协方差矩阵分解为Σxx、Σxy、Σyx、Σyy四个子块其中Σxy和Σyx反映了两期影像之间的交叉相关性。然后求解Σxy·Σyy⁻¹·Σyx相对于Σxx的广义特征值得到MAD变换系数。function [mads, A, B] MADGet(img1, img2) % img1, img2: 两期影像维度为 [行数, 列数, 波段数] [rows, cols, bands] size(img1); X reshape(img1, rows * cols, bands); Y reshape(img2, rows * cols, bands); % 计算联合协方差矩阵使用加权版本为IR-MAD预留接口 S covw([X, Y], ones(rows * cols, 1)); Sxx S(1:bands, 1:bands); Sxy S(1:bands, bands1:end); Syy S(bands1:end, bands1:end); Syx Sxy; % 求解广义特征值问题 [A, D] eig(Sxy * inv(Syy) * Syx, Sxx); [~, idx] sort(diag(D), ascend); A A(:, idx); % 计算第二组系数并输出MAD分量 B -inv(Syy) * Syx * A; mads X * A - Y * B; mads reshape(mads, rows, cols, bands); end这段代码的关键在于广义特征值分解的顺序。eig(A, B)求解的是AvλBv的问题这里A是Sxy * inv(Syy) * SyxB是Sxx。特征值从小到大排序后对应的特征向量依次构成了MAD分量越靠前的分量对变化越不敏感越靠后的分量越能捕捉异常。实际使用中通常取最后一个或最后几个MAD分量来生成变化强度图。2.2 IR-MAD的权重迭代如何压制异常像元干扰IR-MADIteratively Reweighted MAD在MAD上的改进集中在迭代加权机制上。MAD假设两期影像的协方差关系是均匀的但变化像元的存在会污染协方差估计导致检测结果偏移。IR-MAD的做法是每轮迭代后统计每个像元所有MAD分量的标准化平方和理论上这个统计量服从自由度为波段数的卡方分布。依据卡方分布的尾部概率为每个像元赋权变化越剧烈的像元权重越小下一轮协方差估计时其影响被有效抑制。IRMAD_Update.m中的权重计算逻辑如下function [result, weights, converged] IRMAD_Update(img1, img2, init_weights, max_iter) [rows, cols, bands] size(img1); X reshape(img1, rows * cols, bands); Y reshape(img2, rows * cols, bands); N rows * cols; weights ones(N, 1); % 初始权重全为1 tol 1e-3; converged false; for iter 1:max_iter % 用当前权重计算加权协方差 S covw([X, Y], weights); Sxx S(1:bands, 1:bands); Sxy S(1:bands, bands1:end); Syy S(bands1:end, bands1:end); Syx Sxy; % 求解MAD变换系数 [A, D] eig(Sxy * inv(Syy) * Syx, Sxx); [~, idx] sort(diag(D), ascend); A A(:, idx); B -inv(Syy) * Syx * A; % 计算MAD分量和卡方统计量 mads X * A - Y * B; variances diag(A * Sxx * A B * Syy * B 2 * A * Sxy * B); chi2_stats sum(mads.^2 ./ variances, 2); % 更新权重卡方分布的生存函数 new_weights 1 - chi2cdf(chi2_stats, bands); new_weights max(new_weights, 1e-6); % 避免除零 % 检查收敛 if norm(new_weights - weights) / norm(weights) tol converged true; end weights new_weights; if converged break; end end result reshape(mads, rows, cols, bands); end这段迭代逻辑是整个IR-MAD算法的灵魂。注意chi2cdf函数需要MATLAB的Statistics and Machine Learning Toolbox如果你没有这个工具箱可以自己实现卡方分布的数值计算——参考eigen2.m中的做法用Gamma函数近似。还有一点值得注意如果两期影像之间发生了大面积变化比如超过30%的区域IR-MAD的收敛会很困难权重均值持续波动这种情况下建议退回到普通MAD或者先做一层粗略的掩膜把明显变化的区域排除掉再做IR-MAD。注意IR-MAD的迭代不是每次都收敛到更好的结果。迭代20轮左右如果权重都已经趋近于0或1说明已经达到了稳定态如果权重一直在0.5附近徘徊多半是影像配准误差太大此时应优先检查几何校正质量而不是继续调迭代参数。3. CVA与PCA从向量差到主成分空间的检测路径3.1 CVA的向量构造与强度阈值化CVAChange Vector Analysis的实现门槛最低计算两期影像每个像元的光谱向量之差取模长作为变化强度取方向作为变化类别标识。CVADemo.m以台州区域的IntensityImage强度影像为输入先逐波段做差再合成变化强度和变化方向。% CVADemo.m 核心片段 [img1, info1] freadenvi(CVAImageTaizhou_IntensityImage.hdr); [img2, info2] freadenvi(PCAImageTaizhou_IntensityImage.hdr); % 计算逐波段差值向量 diff_vec double(img2) - double(img1); % 变化强度 差值向量的欧氏模长 change_magnitude sqrt(sum(diff_vec.^2, 3)); % 变化方向 向量角度多波段时用余弦方向 direction atan2(diff_vec(:,:,2), diff_vec(:,:,1)); % 自适应阈值均值 2倍标准差 thr mean(change_magnitude(:)) 2 * std(change_magnitude(:)); binary_cva change_magnitude thr;CVA的结果直观变化强度图可以直接出灰度图变化方向图可以用来区分不同类型的转变比如植被变裸地和水体变植被方向是不同的。但CVA对预处理极其敏感太阳高度角差异、气溶胶光学厚度变化、传感器响应差异都会在差值向量中引入系统偏差。所以CVA之前辐射归一化是必须的前置步骤。3.2 PCA的堆叠变换为什么方差最大的方向就是变化方向PCA在变化检测中的应用方式和单独做数据降维不同。这套代码采用的思路是把两期影像的所有波段堆叠成一个高维向量然后做主成分变换。堆叠后的数据中未变化区域的像元在各波段上的取值高度相关它们的方差主要来自噪声而变化区域的像元会沿某些特定方向产生更大的离差这些方向通常对应主成分分析中方差贡献最大的几个主成分。这里有一个关键认知主成分分析并不直接区分变化和噪声它只是按方差大小排序。所以当传感器噪声水平较高时第一主成分可能主要捕捉噪声而不是变化。这就解释了为什么PCADemo.m中会生成PCAImageTaizhou_IntensityImage和PCAImageTaizhou_BinaryValue两个图层——前者是主成分得分图后者是进一步做阈值分割后的二值变化图。function [score, coeff, latent] PCADemo(img1, img2) [rows, cols, bands] size(img1); X1 reshape(img1, rows * cols, bands); X2 reshape(img2, rows * cols, bands); % 两期影像堆叠 X [X1, X2]; X X - mean(X, 1); % 中心化 % SVD分解比eig数值更稳定 [U, S, V] svd(X, econ); latent diag(S).^2 / (size(X, 1) - 1); % 取前k个主成分这里取一半维度 k size(X, 2) / 2; coeff V(:, 1:k); score X * coeff; % 计算变化强度前k个主成分得分的欧氏距离 change_score sqrt(sum(score.^2, 2)); score reshape(score, rows, cols, k); endeigen2.m在这套代码中的作用是提供特征值分解的封装。如果你的数据维度较小直接用eig没问题但高光谱影像波段数多协方差矩阵达到几百乘几百svd的数值稳定性远好于eig。一个从实践中得到的经验当条件数超过10⁶时eig的结果会出现虚部噪声svd则不会。3.3 四种算法的结果形式对比算法输出形式变化判定依据对预处理的敏感度计算复杂度MAD多通道MAD分量卡方检验概率中等O(N·B²)IR-MAD多通道分量权重图加权卡方概率中等O(iter·N·B²)CVA变化强度图方向图模长阈值高需辐射校正O(N·B)PCA主成分得分图得分距离阈值低自动去冗余O(N·B²)从表中能看出如果追求计算速度CVA是最快选择如果需要自动化分析和抗噪声能力PCA更合适如果两期影像来自不同传感器IR-MAD的迭代加权机制能最大程度抵消辐射差异的影响。4. MATLAB工作流从ENVI数据到变化检测结果4.1 ENVI格式解析与数据预处理要点这套代码必须依赖freadenvi.m和enviwrite.m这两个核心IO函数它们负责读取和写出ENVI标准格式数据。ENVI格式由数据文件和.hdr头文件组成头文件的关键字段包括samples列数、lines行数、bands波段数、data type数据类型、interleave存储格式BSQ/BIL/BIP。遥感影像预处理是决定变化检测成败的第一步。拿到数据后不要急着跑算法先确认三件事两期影像的行列数是否一致、投影是否一致、像元是否严格配准。建议在arcmap或arcgis pro中加载遥感影像叠加两个时相的图层检查同名地物的偏移量。如果偏移超过半个像元先做配准再继续。另一个常见的坑是数据类型uint8和uint16直接做差值运算会溢出必须先用double()做类型转换。% 读取ENVI影像并做基础预处理 [img1, info1] freadenvi(CVAImageTaizhou_IntensityImage.hdr); [img2, info2] freadenvi(PCAImageTaizhou_IntensityImage.hdr); % 转为double并统一到[0,1]范围 img1 double(img1); img2 double(img2); if max(img1(:)) 1 img1 img1 / max(img1(:)); end if max(img2(:)) 1 img2 img2 / max(img2(:)); end % 裁剪共同区域如果行数列数不一致取交集 common_rows min(size(img1, 1), size(img2, 1)); common_cols min(size(img1, 2), size(img2, 2)); img1 img1(1:common_rows, 1:common_cols, :); img2 img2(1:common_rows, 1:common_cols, :);freadenvi.m解析头文件时需要注意data type的映射关系。ENVI中类型1是uint82是int163是int324是float325是double6是complex649是complex128。读取时用MATLAB的fread的精度参数对应即可。如果读出来的图像出现隔行错位现象多半是interleave字段解析错了——BSQ是波段顺序存储先存完第一波段再存第二波段BIL是逐行交替存储波段BIP是逐像元交替存储波段三者reshape的顺序完全不同。4.2 New_main.m中四种算法的串联执行New_main.m是整套代码的入口脚本它串联了四个算法的调用。执行流程分为五步读取影像、逐算法检测、K均值聚类生成二值图、结果对比、可视化。KmeansMap.m在这里的作用不是传统意义的分类而是把连续的变化强度图聚成变化/未变化两个类别这种方式的优势是无需人工设定阈值。% New_main.m 主流程 close all; clear; clc; % Step 1: 读取两期影像 hdr1 CVAImageTaizhou_IntensityImage.hdr; hdr2 PCAImageTaizhou_IntensityImage.hdr; [img1, info1] freadenvi(hdr1); [img2, info2] freadenvi(hdr2); % Step 2: 运行MAD和IR-MAD [MAD_result, A, B] MADGet(img1, img2); [IRMAD_result, weights, flag] IRMAD_Update(img1, img2, [], 30); % Step 3: 运行CVA [CVA_intensity, CVA_direction] CVA_Compute(img1, img2); % Step 4: 运行PCA [PCA_result] PCADemo(img1, img2); % Step 5: 将各算法的变化强度图统一做KMeans二值分割 MAD_binary KmeansMap(sum(MAD_result.^2, 3), 2); IRMAD_binary KmeansMap(sum(IRMAD_result.^2, 3) .* (weights 0.5), 2); CVA_binary KmeansMap(CVA_intensity, 2); PCA_binary KmeansMap(PCA_result, 2); % Step 6: 保存结果 enviwrite(MAD_binary, MAD_Binary, float32); enviwrite(IRMAD_binary, IRMAD_Binary, float32); enviwrite(CVA_binary, CVA_Binary, float32); enviwrite(PCA_binary, PCA_Binary, float32); % Step 7: 可视化对比 createfigure(MAD_binary, IRMAD_binary, CVA_binary, PCA_binary);值得注意的是IR-MAD的二值分割用权重做了掩膜处理。权重小于0.5的像元被认为在迭代中置信度低很多是边缘混合像元或配准残差区直接丢弃可以减少椒盐噪声。createfigure.m负责绘制对比图通常输出一个2×2的子图矩阵。如果要把结果叠加到真实影像上查看推荐使用enviwrite把二值图写出再在arcmap或arcgis pro中叠加显示。4.3 KMeans参数与收敛问题的调整KmeansMap.m调用MATLAB内置的kmeans函数默认使用K-means初始化。对于变化检测的二值分割需要设置Distance为sqEuclideanReplicates设为3~5次以避免局部最优。如果聚类结果不稳定每次运行结果都不一样是因为初始中心选择不同增加Replicates次数取最小SSE的结果即可。function binary_map KmeansMap(intensity, k) % intensity: 单通道变化强度图 % k: 聚类类别数变化检测中通常为2 pix intensity(:); % 去除NaN和Inf valid isfinite(pix); [idx, centroid] kmeans(pix(valid), k, ... Distance, sqEuclidean, ... Replicates, 5, ... MaxIter, 200, ... Display, off); % 找到强度均值较大的类别作为变化类 [~, change_class] max(centroid); binary_map zeros(size(pix)); temp zeros(length(pix), 1); temp(valid) idx; binary_map reshape(temp change_class, size(intensity)); end另一个容易忽略的细节K-means对输入数据的尺度敏感如果强度图中有极大离群值比如云或阴影区域聚类中心会被拉偏。建议在输入K-means前先做百分位截断把强度图clip到1%和99%分位数之间p_low prctile(intensity(:), 1); p_high prctile(intensity(:), 99); intensity_clipped min(max(intensity, p_low), p_high);这个预处理对阈值类方法同样有效能在不引入人工判断的前提下显著提升二值分割的稳定性。5. 结果验证与精度评估混淆矩阵与Kappa系数的落地计算拿到四张变化检测二值图怎么判断哪个效果好直接目视对比在arcmap或arcgis pro中叠加遥感影像并创建新的地类图斑矢量图层作为真值然后逐算法计算混淆矩阵。Change Result Comparing.m已经实现了结果对比框架但缺少数值评估模块这里给出补充方法。function [Kappa, F1, OA] evaluate_change(predicted, ground_truth) % predicted: 算法输出的二值变化图1变化0未变化 % ground_truth: 人工标注的真值二值图 p predicted(:); g ground_truth(:); % 混淆矩阵元素 TP sum(p 1 g 1); FP sum(p 1 g 0); FN sum(p 0 g 1); TN sum(p 0 g 0); OA (TP TN) / length(p); % F1分数 precision TP / (TP FP eps); recall TP / (TP FN eps); F1 2 * precision * recall / (precision recall eps); % Kappa系数 p0 OA; pe ((TP FP) * (TP FN) (FN TN) * (FP TN)) / length(p)^2; Kappa (p0 - pe) / (1 - pe eps); end实际项目中如果验证区域是耕地变化检测数据集这类带专题标签的数据建议按地类分层抽样评估。不同算法对线状地物道路、沟渠和面状地物耕地、建筑的检测能力差异很大CVA倾向于把线状地物检出为细碎像元IR-MAD对大面积连片变化更敏感。一个工程上的技巧是把四个算法的结果做投票集成至少两个算法判定为变化的像元才标记为变化这个简单操作通常能把F1提高3~5个百分点。在动手跑这些代码之前最后提醒一点这套代码里的LastWeight.mat保存了上一次IR-MAD迭代的权重矩阵。如果你处理的影像和台州样例数据特征差异较大务必清空这个变量重新训练不要让历史权重影响新数据的协方差估计。直接删除或注释掉load(LastWeight.mat)这一行让IRMAD_Update.m从全1权重开始迭代即可。本文还有配套的精品资源点击获取 遥感变化检测经典算法：MAD、IR-MAD、CVA与PCA的MATLAB实现
简介压缩包内汇总了遥感影像变化检测中四种经典算法的 Matlab 实现面向遥感、测绘及地理信息专业的科研人员和学生可用于多时相地表变化分析。包内共有97个文件包括.m算法脚本、ENVI标准.hdr/.tif影像数据、用于结果展示的.fig图表以及若干.bmp辅助图像压缩后约10.4MB文件组织清晰便于对照学习。其中包含IR-MAD、MAD、CVA、PCA各自的Demo与核心函数并附泰州区域的真实影像及中间结果读者可直接运行观察强度图、二值图与变化检测对比效果。此外还包括enviwrite.m、covw.m等读写与统计工具方便将算法迁移到自己的数据集。目前已有约2008人学习下载适合想通过现成代码理解算法原理、开展实验验证的入门与进阶使用者。1. 为什么传统遥感变化检测算法仍是工程首选十年前做地物变化分析大家用ENVI里的工具一步一步点现在一提起变化检测很多人第一反应就是上深度学习。但你真正落过地就会发现神经网络需要标签、需要训练数据而且在跨传感器、跨时相的影像上泛化能力并不稳。IR-MAD、MAD、CVA、PCA这四种经典算法胜在不依赖样本、数学可解释、跑起来快。这套MATLAB代码以台州区域的ENVI格式影像为样例完整实现了从影像读取、预处理到变化结果输出的全链路。对刚接触遥感影像预处理的研究生来说是理解统计类变化检测最好的入门素材对已经用深度学习做过语义分割复现的工程师同样适合拿来做候选区筛选和基线对比。2. MAD与IR-MAD从统计异常到迭代加权的演进逻辑2.1 MAD的数学建模线性组合如何分离变化信息MADMultitemporal Anomaly Detection多时相异常检测的数学内核是典型相关分析。它把两期影像看作两个随机向量X和Y试图找到一组线性组合系数a和b使得组合后的变量UaᵀX和VbᵀY之间的相关性最小。当两期影像在同一位置的像元未发生变化时U−V的值应该接近零而真正变化的区域会因为光谱关系被破坏而产生显著的偏离值。这套代码里的MADGet.m实现了完整的求解过程。核心是构造两期影像的联合协方差矩阵将其分块后求解广义特征值问题。协方差矩阵分解为Σxx、Σxy、Σyx、Σyy四个子块其中Σxy和Σyx反映了两期影像之间的交叉相关性。然后求解Σxy·Σyy⁻¹·Σyx相对于Σxx的广义特征值得到MAD变换系数。function [mads, A, B] MADGet(img1, img2) % img1, img2: 两期影像维度为 [行数, 列数, 波段数] [rows, cols, bands] size(img1); X reshape(img1, rows * cols, bands); Y reshape(img2, rows * cols, bands); % 计算联合协方差矩阵使用加权版本为IR-MAD预留接口 S covw([X, Y], ones(rows * cols, 1)); Sxx S(1:bands, 1:bands); Sxy S(1:bands, bands1:end); Syy S(bands1:end, bands1:end); Syx Sxy; % 求解广义特征值问题 [A, D] eig(Sxy * inv(Syy) * Syx, Sxx); [~, idx] sort(diag(D), ascend); A A(:, idx); % 计算第二组系数并输出MAD分量 B -inv(Syy) * Syx * A; mads X * A - Y * B; mads reshape(mads, rows, cols, bands); end这段代码的关键在于广义特征值分解的顺序。eig(A, B)求解的是AvλBv的问题这里A是Sxy * inv(Syy) * SyxB是Sxx。特征值从小到大排序后对应的特征向量依次构成了MAD分量越靠前的分量对变化越不敏感越靠后的分量越能捕捉异常。实际使用中通常取最后一个或最后几个MAD分量来生成变化强度图。2.2 IR-MAD的权重迭代如何压制异常像元干扰IR-MADIteratively Reweighted MAD在MAD上的改进集中在迭代加权机制上。MAD假设两期影像的协方差关系是均匀的但变化像元的存在会污染协方差估计导致检测结果偏移。IR-MAD的做法是每轮迭代后统计每个像元所有MAD分量的标准化平方和理论上这个统计量服从自由度为波段数的卡方分布。依据卡方分布的尾部概率为每个像元赋权变化越剧烈的像元权重越小下一轮协方差估计时其影响被有效抑制。IRMAD_Update.m中的权重计算逻辑如下function [result, weights, converged] IRMAD_Update(img1, img2, init_weights, max_iter) [rows, cols, bands] size(img1); X reshape(img1, rows * cols, bands); Y reshape(img2, rows * cols, bands); N rows * cols; weights ones(N, 1); % 初始权重全为1 tol 1e-3; converged false; for iter 1:max_iter % 用当前权重计算加权协方差 S covw([X, Y], weights); Sxx S(1:bands, 1:bands); Sxy S(1:bands, bands1:end); Syy S(bands1:end, bands1:end); Syx Sxy; % 求解MAD变换系数 [A, D] eig(Sxy * inv(Syy) * Syx, Sxx); [~, idx] sort(diag(D), ascend); A A(:, idx); B -inv(Syy) * Syx * A; % 计算MAD分量和卡方统计量 mads X * A - Y * B; variances diag(A * Sxx * A B * Syy * B 2 * A * Sxy * B); chi2_stats sum(mads.^2 ./ variances, 2); % 更新权重卡方分布的生存函数 new_weights 1 - chi2cdf(chi2_stats, bands); new_weights max(new_weights, 1e-6); % 避免除零 % 检查收敛 if norm(new_weights - weights) / norm(weights) tol converged true; end weights new_weights; if converged break; end end result reshape(mads, rows, cols, bands); end这段迭代逻辑是整个IR-MAD算法的灵魂。注意chi2cdf函数需要MATLAB的Statistics and Machine Learning Toolbox如果你没有这个工具箱可以自己实现卡方分布的数值计算——参考eigen2.m中的做法用Gamma函数近似。还有一点值得注意如果两期影像之间发生了大面积变化比如超过30%的区域IR-MAD的收敛会很困难权重均值持续波动这种情况下建议退回到普通MAD或者先做一层粗略的掩膜把明显变化的区域排除掉再做IR-MAD。注意IR-MAD的迭代不是每次都收敛到更好的结果。迭代20轮左右如果权重都已经趋近于0或1说明已经达到了稳定态如果权重一直在0.5附近徘徊多半是影像配准误差太大此时应优先检查几何校正质量而不是继续调迭代参数。3. CVA与PCA从向量差到主成分空间的检测路径3.1 CVA的向量构造与强度阈值化CVAChange Vector Analysis的实现门槛最低计算两期影像每个像元的光谱向量之差取模长作为变化强度取方向作为变化类别标识。CVADemo.m以台州区域的IntensityImage强度影像为输入先逐波段做差再合成变化强度和变化方向。% CVADemo.m 核心片段 [img1, info1] freadenvi(CVAImageTaizhou_IntensityImage.hdr); [img2, info2] freadenvi(PCAImageTaizhou_IntensityImage.hdr); % 计算逐波段差值向量 diff_vec double(img2) - double(img1); % 变化强度 差值向量的欧氏模长 change_magnitude sqrt(sum(diff_vec.^2, 3)); % 变化方向 向量角度多波段时用余弦方向 direction atan2(diff_vec(:,:,2), diff_vec(:,:,1)); % 自适应阈值均值 2倍标准差 thr mean(change_magnitude(:)) 2 * std(change_magnitude(:)); binary_cva change_magnitude thr;CVA的结果直观变化强度图可以直接出灰度图变化方向图可以用来区分不同类型的转变比如植被变裸地和水体变植被方向是不同的。但CVA对预处理极其敏感太阳高度角差异、气溶胶光学厚度变化、传感器响应差异都会在差值向量中引入系统偏差。所以CVA之前辐射归一化是必须的前置步骤。3.2 PCA的堆叠变换为什么方差最大的方向就是变化方向PCA在变化检测中的应用方式和单独做数据降维不同。这套代码采用的思路是把两期影像的所有波段堆叠成一个高维向量然后做主成分变换。堆叠后的数据中未变化区域的像元在各波段上的取值高度相关它们的方差主要来自噪声而变化区域的像元会沿某些特定方向产生更大的离差这些方向通常对应主成分分析中方差贡献最大的几个主成分。这里有一个关键认知主成分分析并不直接区分变化和噪声它只是按方差大小排序。所以当传感器噪声水平较高时第一主成分可能主要捕捉噪声而不是变化。这就解释了为什么PCADemo.m中会生成PCAImageTaizhou_IntensityImage和PCAImageTaizhou_BinaryValue两个图层——前者是主成分得分图后者是进一步做阈值分割后的二值变化图。function [score, coeff, latent] PCADemo(img1, img2) [rows, cols, bands] size(img1); X1 reshape(img1, rows * cols, bands); X2 reshape(img2, rows * cols, bands); % 两期影像堆叠 X [X1, X2]; X X - mean(X, 1); % 中心化 % SVD分解比eig数值更稳定 [U, S, V] svd(X, econ); latent diag(S).^2 / (size(X, 1) - 1); % 取前k个主成分这里取一半维度 k size(X, 2) / 2; coeff V(:, 1:k); score X * coeff; % 计算变化强度前k个主成分得分的欧氏距离 change_score sqrt(sum(score.^2, 2)); score reshape(score, rows, cols, k); endeigen2.m在这套代码中的作用是提供特征值分解的封装。如果你的数据维度较小直接用eig没问题但高光谱影像波段数多协方差矩阵达到几百乘几百svd的数值稳定性远好于eig。一个从实践中得到的经验当条件数超过10⁶时eig的结果会出现虚部噪声svd则不会。3.3 四种算法的结果形式对比算法输出形式变化判定依据对预处理的敏感度计算复杂度MAD多通道MAD分量卡方检验概率中等O(N·B²)IR-MAD多通道分量权重图加权卡方概率中等O(iter·N·B²)CVA变化强度图方向图模长阈值高需辐射校正O(N·B)PCA主成分得分图得分距离阈值低自动去冗余O(N·B²)从表中能看出如果追求计算速度CVA是最快选择如果需要自动化分析和抗噪声能力PCA更合适如果两期影像来自不同传感器IR-MAD的迭代加权机制能最大程度抵消辐射差异的影响。4. MATLAB工作流从ENVI数据到变化检测结果4.1 ENVI格式解析与数据预处理要点这套代码必须依赖freadenvi.m和enviwrite.m这两个核心IO函数它们负责读取和写出ENVI标准格式数据。ENVI格式由数据文件和.hdr头文件组成头文件的关键字段包括samples列数、lines行数、bands波段数、data type数据类型、interleave存储格式BSQ/BIL/BIP。遥感影像预处理是决定变化检测成败的第一步。拿到数据后不要急着跑算法先确认三件事两期影像的行列数是否一致、投影是否一致、像元是否严格配准。建议在arcmap或arcgis pro中加载遥感影像叠加两个时相的图层检查同名地物的偏移量。如果偏移超过半个像元先做配准再继续。另一个常见的坑是数据类型uint8和uint16直接做差值运算会溢出必须先用double()做类型转换。% 读取ENVI影像并做基础预处理 [img1, info1] freadenvi(CVAImageTaizhou_IntensityImage.hdr); [img2, info2] freadenvi(PCAImageTaizhou_IntensityImage.hdr); % 转为double并统一到[0,1]范围 img1 double(img1); img2 double(img2); if max(img1(:)) 1 img1 img1 / max(img1(:)); end if max(img2(:)) 1 img2 img2 / max(img2(:)); end % 裁剪共同区域如果行数列数不一致取交集 common_rows min(size(img1, 1), size(img2, 1)); common_cols min(size(img1, 2), size(img2, 2)); img1 img1(1:common_rows, 1:common_cols, :); img2 img2(1:common_rows, 1:common_cols, :);freadenvi.m解析头文件时需要注意data type的映射关系。ENVI中类型1是uint82是int163是int324是float325是double6是complex649是complex128。读取时用MATLAB的fread的精度参数对应即可。如果读出来的图像出现隔行错位现象多半是interleave字段解析错了——BSQ是波段顺序存储先存完第一波段再存第二波段BIL是逐行交替存储波段BIP是逐像元交替存储波段三者reshape的顺序完全不同。4.2 New_main.m中四种算法的串联执行New_main.m是整套代码的入口脚本它串联了四个算法的调用。执行流程分为五步读取影像、逐算法检测、K均值聚类生成二值图、结果对比、可视化。KmeansMap.m在这里的作用不是传统意义的分类而是把连续的变化强度图聚成变化/未变化两个类别这种方式的优势是无需人工设定阈值。% New_main.m 主流程 close all; clear; clc; % Step 1: 读取两期影像 hdr1 CVAImageTaizhou_IntensityImage.hdr; hdr2 PCAImageTaizhou_IntensityImage.hdr; [img1, info1] freadenvi(hdr1); [img2, info2] freadenvi(hdr2); % Step 2: 运行MAD和IR-MAD [MAD_result, A, B] MADGet(img1, img2); [IRMAD_result, weights, flag] IRMAD_Update(img1, img2, [], 30); % Step 3: 运行CVA [CVA_intensity, CVA_direction] CVA_Compute(img1, img2); % Step 4: 运行PCA [PCA_result] PCADemo(img1, img2); % Step 5: 将各算法的变化强度图统一做KMeans二值分割 MAD_binary KmeansMap(sum(MAD_result.^2, 3), 2); IRMAD_binary KmeansMap(sum(IRMAD_result.^2, 3) .* (weights 0.5), 2); CVA_binary KmeansMap(CVA_intensity, 2); PCA_binary KmeansMap(PCA_result, 2); % Step 6: 保存结果 enviwrite(MAD_binary, MAD_Binary, float32); enviwrite(IRMAD_binary, IRMAD_Binary, float32); enviwrite(CVA_binary, CVA_Binary, float32); enviwrite(PCA_binary, PCA_Binary, float32); % Step 7: 可视化对比 createfigure(MAD_binary, IRMAD_binary, CVA_binary, PCA_binary);值得注意的是IR-MAD的二值分割用权重做了掩膜处理。权重小于0.5的像元被认为在迭代中置信度低很多是边缘混合像元或配准残差区直接丢弃可以减少椒盐噪声。createfigure.m负责绘制对比图通常输出一个2×2的子图矩阵。如果要把结果叠加到真实影像上查看推荐使用enviwrite把二值图写出再在arcmap或arcgis pro中叠加显示。4.3 KMeans参数与收敛问题的调整KmeansMap.m调用MATLAB内置的kmeans函数默认使用K-means初始化。对于变化检测的二值分割需要设置Distance为sqEuclideanReplicates设为3~5次以避免局部最优。如果聚类结果不稳定每次运行结果都不一样是因为初始中心选择不同增加Replicates次数取最小SSE的结果即可。function binary_map KmeansMap(intensity, k) % intensity: 单通道变化强度图 % k: 聚类类别数变化检测中通常为2 pix intensity(:); % 去除NaN和Inf valid isfinite(pix); [idx, centroid] kmeans(pix(valid), k, ... Distance, sqEuclidean, ... Replicates, 5, ... MaxIter, 200, ... Display, off); % 找到强度均值较大的类别作为变化类 [~, change_class] max(centroid); binary_map zeros(size(pix)); temp zeros(length(pix), 1); temp(valid) idx; binary_map reshape(temp change_class, size(intensity)); end另一个容易忽略的细节K-means对输入数据的尺度敏感如果强度图中有极大离群值比如云或阴影区域聚类中心会被拉偏。建议在输入K-means前先做百分位截断把强度图clip到1%和99%分位数之间p_low prctile(intensity(:), 1); p_high prctile(intensity(:), 99); intensity_clipped min(max(intensity, p_low), p_high);这个预处理对阈值类方法同样有效能在不引入人工判断的前提下显著提升二值分割的稳定性。5. 结果验证与精度评估混淆矩阵与Kappa系数的落地计算拿到四张变化检测二值图怎么判断哪个效果好直接目视对比在arcmap或arcgis pro中叠加遥感影像并创建新的地类图斑矢量图层作为真值然后逐算法计算混淆矩阵。Change Result Comparing.m已经实现了结果对比框架但缺少数值评估模块这里给出补充方法。function [Kappa, F1, OA] evaluate_change(predicted, ground_truth) % predicted: 算法输出的二值变化图1变化0未变化 % ground_truth: 人工标注的真值二值图 p predicted(:); g ground_truth(:); % 混淆矩阵元素 TP sum(p 1 g 1); FP sum(p 1 g 0); FN sum(p 0 g 1); TN sum(p 0 g 0); OA (TP TN) / length(p); % F1分数 precision TP / (TP FP eps); recall TP / (TP FN eps); F1 2 * precision * recall / (precision recall eps); % Kappa系数 p0 OA; pe ((TP FP) * (TP FN) (FN TN) * (FP TN)) / length(p)^2; Kappa (p0 - pe) / (1 - pe eps); end实际项目中如果验证区域是耕地变化检测数据集这类带专题标签的数据建议按地类分层抽样评估。不同算法对线状地物道路、沟渠和面状地物耕地、建筑的检测能力差异很大CVA倾向于把线状地物检出为细碎像元IR-MAD对大面积连片变化更敏感。一个工程上的技巧是把四个算法的结果做投票集成至少两个算法判定为变化的像元才标记为变化这个简单操作通常能把F1提高3~5个百分点。在动手跑这些代码之前最后提醒一点这套代码里的LastWeight.mat保存了上一次IR-MAD迭代的权重矩阵。如果你处理的影像和台州样例数据特征差异较大务必清空这个变量重新训练不要让历史权重影响新数据的协方差估计。直接删除或注释掉load(LastWeight.mat)这一行让IRMAD_Update.m从全1权重开始迭代即可。本文还有配套的精品资源点击获取