简介：针对遥感图像变化检测任务，基于主成分分析（PCA）与K-Means聚类的无监督算法无需标签数据，即可通过对比不同时相的卫星影像识别地表显著变化。该资源面向图像处理、数据挖掘和遥感应用开发者，提供了一套完整的MATLAB实现，便于学习算法原理并移植到实际项目中。压缩包共包含3个文件，分别是.m源码文件、.md说明文档和.pdf原理论文，整体仅354KB，轻量且易于下载。源码覆盖预处理、PCA降维、主成分选取、K-Means聚类、变化检测及后处理等完整链路，读者可以对照文档逐步调试参数，理解每个环节的作用；PDF论文则提供了算法背景与实验结果，便于深入掌握无监督变化检测的来龙去脉。目前已有668人学习查看，适用于环境监测、城市扩张分析、灾害评估等典型遥感场景，也适合作为课程设计或毕业设计的算法参考。
1. 为什么变化检测要选PCA和K-Means？一个不用标签的切入点
拿到两个时相的遥感影像，没有地面真值，却需要在几天内圈出新增建筑或毁林范围，这种任务不该先想着攒标签。直接对差分影像做阈值，辐射差异和配准误差会制造大量假变化。把PCA和K-Means串成无监督管道，等于先用主成分分析把影像的高频噪声和多波段相关性剥掉，再把像素在低维特征空间里按距离聚成几类，最后比较类别归属是否有跳变。整个过程不需要一个标注样本，也绕开了训练样本不均衡的问题。这个思路在MATLAB里落地很快，ChangeDetection_PCA_KMeans.m就是这个流程的骨架示例，适合遥感算法复现、无监督学习入门以及想快速出变化草图的工程场景。
2. 从影像对到特征矩阵：PCA在双时相变化检测里的数学角色
2.1 变化检测的本质：把"变了没有"变成"特征是否漂移"
变化检测不是分类，不关心像素本身是水体还是房屋，它只回答"同一位置在两个时间点是否发生显著差异"。最朴素的方法是直接计算两时相的差，再设置阈值把大于阈值的像素标为变化。这个做法在理想情况下成立，但真实遥感影像存在三个干扰源：传感器噪声、太阳高度角差异引起的辐射变化、以及影像配准后的亚像素错位。三者叠加会使得差分图上的噪声分布不是高斯，固定阈值很难同时兼顾漏检和误检。PCA从这个困境里给出另一条路：把影像从波段空间投影到一组按方差排列的正交方向，前几个方向保留了地物结构的主能量，噪声被压到后面的成分中。于是"有没有变化"被转化为"同一像素在低维特征空间里的位置是否发生漂移"，这种表示更鲁棒。
2.2 PCA如何用协方差矩阵抓住主要变化特征
设一幅影像有B个波段，重排成N行B列的矩阵X。先对每列减去均值，得到中心化矩阵。B个波段之间的协方差矩阵是一个B×B对称半正定矩阵，特征值分解后得到特征向量（主轴方向）和特征值（方差大小）。把中心化数据乘上取前n列的特征向量矩阵，就得到投影后的主成分。有一点要注意：PCA对量纲敏感，如果传感器的各波段曝光时间不同导致数值范围差异大，需要先按列除以标准差做标准化。不过对同一传感器同时获取的多光谱影像，各波段数值范围一致，做中心化就够了，强行标准化反而会放大噪声波段的贡献。
在双时相场景里，关键问题是PCA的基向量应该从哪里来。常见做法是把两个时相的像素全部堆叠成一个(2N)×B的矩阵再算PCA，好处是两个时相共享同一套特征基向量，主成分坐标可以直接做减法；如果分别对两幅影像做PCA，基向量不同，后续比较就需要先做向量空间对齐，徒增麻烦。
2.3 用MATLAB验证PCA降维效果
下面对一幅8波段的多光谱影像跑一次PCA，观察方差分布。代码用的是MATLAB统计工具箱的pca函数。
im1 = imread('time1.tif'); % 假设已经完成几何配准 im2 = imread('time2.tif'); [H, W, B] = size(im1); X = double(reshape(im1, H*W, B)); Y = double(reshape(im2, H*W, B)); C = [X; Y]; % 两时相堆叠 C = C - mean(C, 1); % 中心化 [coeff, score, latent] = pca(C); % coeff: 特征向量，latent: 特征值 explained = 100 * latent / sum(latent); cumExplained = cumsum(explained); plot(1:B, cumExplained, 'o-'); grid on; xlabel('主成分序号'); ylabel('累积方差贡献率 (%)');这里pca默认已经做了中心化，手动减去均值是为了让score解释更直观。coeff的每一列是一个主成分方向，按latent降序排列。score里行数等于C的行数，前半部分是时相1的投影坐标，后半部分对应时相2。cumExplained可以直接看出前3个主成分是否已经达到85%。如果达不到，说明波段间相关性弱，需要保留更多成分，或者检查数据是否存在坏线等离群值。
2.4 选多少个主成分？方差贡献率的工程判据
| 累计贡献率 | 主成分数量含义 | 适用场景 |
|---|---|---|
| < 85% | 信息压缩过度，细碎变化易丢失 | 大范围粗检测、快速预览 |
| 85% ~ 95% | 常规选择，保留主要地物结构 | 建筑扩张、植被覆盖变化 |
| > 95%、< 99% | 噪声开始明显，但能保留小目标 | 水体边界微变、高光谱数据 |
| > 99% | 几乎所有信息，维度仍较高 | 数值实验、严谨对比时备用 |
这个经验区间来自我处理Landsat和Sentinel-2数据的习惯，不同传感器差异很大。高光谱上百个波段用0.95截断仍然可能保留40多个维度，而多光谱的3~4个主成分通常就足够。实际调参时可以把nComp设成2到5跑一遍，用第5章里提到的kappa系数来收口，不要只凭贡献率决定。
3. K-Means聚类定位变化像素：MATLAB逐句拆解
3.1 为什么在主成分空间上聚类，而不是原始像素空间
K-Means通过迭代把样本分配到最近的聚类中心，衡量距离时用的是欧氏距离。多光谱波段之间常常存在高相关性，比如近红外和红波段的植被信息重叠，直接用原始波段距离会放大相关区域的权重；而且波段数量接近十维后，距离会逐渐趋同。PCA做完之后，各主成分正交且方差递减，正交消除了相关性，方差递减让前几个维度主导距离计算。在这个空间里做聚类，簇的形状更接近球形，K-Means的假设更容易满足。另一个实际原因是计算量，主成分从10维降到3维后，kmeans的迭代距离运算量大幅下降，对几百兆像素的大幅影像尤其明显。
3.2 从两个时相的主成分到变化特征
这里要决定把什么喂给K-Means。最直接的方案是取两时相的差异得分作为特征，featDiff = s2(:,1:nComp) - s1(:,1:nComp)，这样得到的是一个N行nComp列的矩阵，每行代表一个像素的变化向量。向量长度代表变化强度，方向代表变化性质。K-Means会把方向接近、幅度接近的分到同一簇。另一种方案是把两个时相的主成分拼接成2*nComp维特征，适合想同时保留时相绝对特征的场景，但维度翻倍，聚类稀疏性变差。变化检测任务里我倾向用差异特征，因为最终关心的就是delta，而不是绝对状态。
两种特征构造方式可以用下面代码切换：
% 假设已经通过共享PCA得到score s1 = score(1:H*W, :); s2 = score(H*W+1:end, :); nComp = 3; % 方式一：差值特征 featDiff = s2(:, 1:nComp) - s1(:, 1:nComp); % 方式二：绝对值差值，忽略变化方向 featAbs = abs(featDiff); % 方式三：拼接特征，保留时相原始信息 featConcat = [s1(:, 1:nComp), s2(:, 1:nComp)];实际使用中，方式一的聚类中心会落在正负两个方向；方式二把所有变化拉到正半轴，聚类中心个数可以少一个；方式三适合两个时相本身地物类别差别很大的场景，比如前一年全场是农田、后一年全场是裸土，这时差值特征会掩盖整体状态迁移，拼接特征反而能区分。
3.3 核心文件ChangeDetection_PCA_KMeans.m的常见流程
脚本的组织方式通常如下：读影像 → 归一化 → PCA投影 → 构造差异特征 → kmeans → reshape成变化图。我按这个顺序实现，关键代码是：
% 取前nComp个主成分 nComp = 3; s1 = score(1:H*W, 1:nComp); % 时相1的主成分得分 s2 = score(H*W+1:end, 1:nComp); % 时相2的主成分得分 % 构造差异特征 featDiff = s2 - s1; % featDiff = abs(s2 - s1); % 如果不关心变化方向，改用绝对值 % K-Means聚类 k = 3; % 簇数量，常见设置为3 rng(2024); % 固定随机种子 [idx, Csum] = kmeans(featDiff, k, ... 'Distance', 'sqEuclidean', ... 'Replicates', 5, ... 'MaxIter', 200); % 按簇中心的范数判断哪个簇是变化类 distFromZero = sqrt(sum(Csum.^2, 2)); changeCluster = find(distFromZero == max(distFromZero)); changeMask = reshape(idx == changeCluster, H, W);参数说明：Distance=sqEuclidean是默认距离，但对于稀疏数据可以考虑'cosine'，不过在主成分空间里欧氏距离更符合方差含义。Replicates=5表示从5个随机初始点出发，返回误差最小的结果，避免局部最优。MaxIter=200是最大迭代次数，太大的值无意义，通常100~200足够收敛。Csum是最终的聚类中心，通过计算与零点的欧氏距离找变化簇，比直接比较标签更可靠。如果聚类中心有正有负，变化簇可能不止一个，比如正方向和负方向各一个，那就把与零点距离超过阈值的多个簇都合并成变化区域。
3.4 聚类后的变化判定规则：比较簇标签还是比较距离？
| 判定方案 | 适合场景 | 需要处理的坑 |
|---|---|---|
| 分别对两个时相聚类，比较簇标签 | 两个时相的类别语义词义明确时 | 簇编号顺序不一致，需要匈牙利算法匹配簇中心 |
| 对差异特征直接聚类 | 只关心变化强度，不要求语义 | 正向变化和负向变化会被分成多个簇，需要按中心距离合并 |
| 对特征差绝对值聚类 | 只关心是否变化，不关心方向 | 丢失"改造相反"的信息，多方向变化被混在一起 |
工程里我优先选第二和第三种混合。先对featDiff做聚类，观察簇中心的空间分布。如果发现有两个簇中心距离零点都很远且方向相反，就把它们同时归入变化类；如果只有一个远，说明该区域内变化单向。注意不能直接把两次分别聚类的结果做差比较，因为两次聚类中心顺序可能整体轮转，比较之前必须按距离最近原则重排簇号，否则会得到大量虚假变化。
4. k值怎么定、后处理怎么做：实战调参与问题排查
4.1 用轮廓系数和肘部法则反复试探k
K-Means要求提前给定k，这个值直接影响变化图。k=2基本等价于"变/不变"二分类；k=3会在正负两个方向各分一个变化簇；k大于等于4时，多个变化簇会把同一变化按强度拆开，你很难跟业务解释"中强度变化"和"高强度变化"的区别。因此实际测试范围2到6就够。
除了轮廓系数，我常用evalclusters批量扫描：
% 采样10%像素做评估，防止数据量太大 rng(42); sampleIdx = randperm(size(featDiff,1), round(size(featDiff,1)*0.1)); featSample = featDiff(sampleIdx, :); eva = evalclusters(featSample, 'kmeans', 'Silhouette', ... 'KList', 2:6); plot(eva); kRecommended = eva.OptimalK;evalclusters的第二个参数可以传函数句柄，也可以传'kmeans'字符串。Silhouette要反复计算样本间的距离，比CalinskiHarabasz慢，但更直观。注意这里只对抽样数据求最优k，不是用全图，否则运行时间会成倍增加。采样比例10%对统计聚类结构足够，但如果你要检测极小的变化目标，采样会漏掉它们，这时应该改用分层采样或直接全图但把样本类型限制为每两百行取一个。
4.2 聚类数k与变化类别的对应关系
k=2适合只有"变"与"不变"的地表类型，比如洪水淹没范围提取。k=3适合既存在地物消失又存在新增的复杂城区，正向变化和负向变化各占一个簇。k=4以上除非有明确的多种变化方向需求，否则不建议。判断当前数据适合几个k，可以先看featDiff的直方图：
histogram(featDiff(:,1), 256); % 观察第一主成分的差如果直方图在0附近只有一个高峰，说明大部分像素没变化，k取2即可；如果两侧出现不对称的拖尾，取3能把正负变化分开；如果两侧都出现多个峰，再考虑k=4。这个步骤不能省，很多人上来就把k设成3，结果变化图把轻微的物候差异也当成独立一类，后期要花几倍时间清洗。
4.3 后处理：形态学开闭运算与连通域去噪
聚类出的变化图是逐像素标签，虽然比阈值可靠，仍然有散点噪声。原因来自两部分：PCA对小块噪声的响应，以及kmeans对边界的抖动。后处理的顺序固定为：中值滤波 → 开运算 → 闭运算 → 面积过滤。
% 得到变化掩模后 changeMask = (idx == changeCluster); changeMask = medfilt2(changeMask, [3 3]); changeMask = imopen(changeMask, strel('disk', 2)); changeMask = imclose(changeMask, strel('disk', 3)); changeMask = bwareaopen(changeMask, 50); % 像素数阈值 % 可选：填充空洞，让变化区闭合 changeMask = imfill(changeMask, 'holes');imopen是先腐蚀再膨胀，删除小于结构元素的小斑点；imclose是先膨胀再腐蚀，填补区域内的空隙。磁盘半径2和3是针对中分辨率影像（10~30m）的经验值，高分辨率影像可以适当增大到5。bwareaopen依据八连通域统计面积，把面积小于阈值的目标去掉，50个像素对于城市变化检测是合理的下限。如果最后结果里变化区域形状变得过于平滑，把bwareaopen的阈值调小，或直接用bwpropfilt按纵横比过滤细长伪变化。
4.4 遇到的坑：数据范围不一致、内存不足、随机初始化的确定性
| 现象 | 常见原因 | 对策 |
|---|---|---|
| 结果出现大面积横条纹 | 两时相辐射归一化没做 | 做直方图匹配或线性回归校准 |
| kmeans报内存不足 | featDiff全图参与迭代 | 抽样估计中心，再用knnsearch分派全图 |
| 相同代码结果每次不同 | 没有固定随机种子 | rng(固定值)，并设置Replicates |
| 变化区域完全错位 | 影像未配准或坐标系不一致 | 检查地理元数据，用imregister重配准 |
| PCA贡献率突然下降 | 影像存在云遮挡或坏像素 | 提前掩膜或插值，别让云参与特征分解 |
这些坑几乎每个遥感数据处理项目都会遇到。内存问题尤其常见，比如Landsat 8全分辨率约3000万像素，构造差异特征后是3000万×3的double矩阵，占720MB，再运行kmeans，迭代时内存会再翻几倍。所以我在大尺度任务里总会走"抽样聚类+全图归类"两步，而不是把全图喂给kmeans。步骤是用randperm取5%像素估计聚类中心，然后对全部像素调用knnsearch找最近中心，速度快且占用小很多。注意抽样前要把no data值用NaN屏蔽，否则NaN在PCA里会传播。
5. 一个提精度的小技巧：用PCA残差辅助K-Means判定变化
5.1 主成分空间上聚类可能漏掉的小变化
PCA的截断既去噪也丢信息。前几个主成分抓的是全局协方差最大的方向，一个只有几百平方米的新建围挡在整幅影像中能量很小，它的信号可能完全落在第五、第六个主成分上。如果你只拿前三个成分做聚类，这个细碎变化基本被洗掉。要补救，可以用PCA残差来生成一个补充变化层。
5.2 构造残差图像并与聚类结果融合
残差定义是原始影像与用前nComp主成分重建影像之间的差。任何一个像素，如果在低维重建后和原始值差异大，说明该像素包含低维主成分没有描述的信息；把这个残差取绝对值并求波段最大值，就能得到一幅针对"小变化"的响应图。再用Otsu对其阈值化，得到残差变化掩模，与K-Means得到的掩模合并。
% 重建两个时相 meanC = mean(C, 1); recon1 = (s1 * coeff(:,1:nComp)') + meanC; recon2 = (s2 * coeff(:,1:nComp)') + meanC; % 两时相重建误差 err1 = abs(double(im1) - reshape(recon1, H, W, B)); err2 = abs(double(im2) - reshape(recon2, H, W, B)); residual = max(err1, err2); % 取两个时相中较大的残差 residualMap = max(residual, [], 3); % 归一化并做Otsu阈值 resNorm = residualMap / max(residualMap(:)); thr = graythresh(resNorm); residMask = imbinarize(resNorm, thr); residMask = bwareaopen(residMask, 20); % 融合策略1：只要一个认为变化就标记，提高召回率 finalMask = changeMask | residMask;residual用max而不是相减，是考虑到变化可能在时相1或者时相2任一侧出现。max(residual, [], 3)取波段维最大值，保证任何一个波段有明显残差都会被保留。融合用逻辑或会提高召回率，适合先圈范围再人工核验；如果希望结果干净，可以改成changeMask & residMask，但会漏掉K-Means能识别、而PCA残差不敏感的大面积渐变变化。具体选哪种取决于后续是用变化图做统计还是做执法取证。
5.3 验证变化检测结果：kappa系数与混淆矩阵
算法说到最后要用数字证明自己，尤其当你想把这个脚本用于生产。如果有标注的真值图，可以这样算：
gt = imread('ref_change.png') > 0; finalMask = imresize(finalMask, size(gt), 'nearest'); % 对齐尺寸 cm = confusionmat(gt(:), finalMask(:)); po = trace(cm) / sum(cm(:)); pe = sum(sum(cm,1) .* sum(cm,2)) / sum(cm(:))^2; kappa = (po - pe) / (1 - pe); fprintf('Kappa: %.3f\n', kappa);confusionmat中第一列是真值背景，第二列是真值变化；对角线之和除以像素总数就是总体精度。kappa的基准是随机分类的期望一致率，超过0.6说明结果有明显一致性，0.8以上可以认为适合业务使用。注意计算之前要确保两幅影像地理范围完全一致，并且变化区域边界的配准误差不要超过一个像素。
本文还有配套的精品资源，点击获取
