企业官网建设流程全解析

简介：遥感影像变化检测经典算法（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：多时相统计异常的第一步

MAD（Multitemporal 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-value，p 值越小，说明该像元是变化的可能性越大。

% 从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：光谱空间里算出变化方向和强度

CVA（Change 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：在降维空间里浓缩变化信息

PCA（Principal 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 选型判断与场景匹配

四种算法没有绝对的优劣，只有适不适合。我用一张表给出选型建议：

如果你手头只有两期影像，不确定预处理能做到什么程度，我的习惯是先跑 IR-MAD，它抗干扰能力最强。如果后续要做变化类型细分，再叠加 CVA 的方向信息。这块选型判断也是近几年遥感变化检测研究里一直讨论的问题，即使现在出现了开放词汇变化检测这类新方向，经典算法的地位依然稳固。

3. 代码包落地复现：从 ENVI 数据到变化图

3.1 文件结构：先读懂这套代码怎么组织

打开压缩包后，看到十多个.m文件和若干.hdr、.mat、无后缀数据文件，第一反应可能有点懵。把文件按功能分组后就清楚了。

主流程是：读取两期 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 值选多少，直接决定结果图的信息粒度。K=2 时就是单纯的变化/不变二分类，K=4 时可以把变化细分成微变、轻度、中度、剧烈四档，但 K 太大也会把噪声单独聚一类。我的建议是结合地表实际情况选 K，比如城市扩张场景下，K=3 到 K=5 比较合理。

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 有严格定义：1=byte，2=integer，12=unsigned integer，4=float，等等。头文件与实际数据存储不一致时，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% 的像素，其余类别几乎看不见，或者聚类结果在多次运

本文还有配套的精品资源，点击获取