简介针对遥感图像变化检测任务基于主成分分析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*W1: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*W1: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);参数说明DistancesqEuclidean是默认距离但对于稀疏数据可以考虑cosine不过在主成分空间里欧氏距离更符合方差含义。Replicates5表示从5个随机初始点出发返回误差最小的结果避免局部最优。MaxIter200是最大迭代次数太大的值无意义通常100~200足够收敛。Csum是最终的聚类中心通过计算与零点的欧氏距离找变化簇比直接比较标签更可靠。如果聚类中心有正有负变化簇可能不止一个比如正方向和负方向各一个那就把与零点距离超过阈值的多个簇都合并成变化区域。3.4 聚类后的变化判定规则比较簇标签还是比较距离判定方案适合场景需要处理的坑分别对两个时相聚类比较簇标签两个时相的类别语义词义明确时簇编号顺序不一致需要匈牙利算法匹配簇中心对差异特征直接聚类只关心变化强度不要求语义正向变化和负向变化会被分成多个簇需要按中心距离合并对特征差绝对值聚类只关心是否变化不关心方向丢失改造相反的信息多方向变化被混在一起工程里我优先选第二和第三种混合。先对featDiff做聚类观察簇中心的空间分布。如果发现有两个簇中心距离零点都很远且方向相反就把它们同时归入变化类如果只有一个远说明该区域内变化单向。注意不能直接把两次分别聚类的结果做差比较因为两次聚类中心顺序可能整体轮转比较之前必须按距离最近原则重排簇号否则会得到大量虚假变化。4. k值怎么定、后处理怎么做实战调参与问题排查4.1 用轮廓系数和肘部法则反复试探kK-Means要求提前给定k这个值直接影响变化图。k2基本等价于变/不变二分类k3会在正负两个方向各分一个变化簇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与变化类别的对应关系k2适合只有变与不变的地表类型比如洪水淹没范围提取。k3适合既存在地物消失又存在新增的复杂城区正向变化和负向变化各占一个簇。k4以上除非有明确的多种变化方向需求否则不建议。判断当前数据适合几个k可以先看featDiff的直方图histogram(featDiff(:,1), 256); % 观察第一主成分的差如果直方图在0附近只有一个高峰说明大部分像素没变化k取2即可如果两侧出现不对称的拖尾取3能把正负变化分开如果两侧都出现多个峰再考虑k4。这个步骤不能省很多人上来就把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以上可以认为适合业务使用。注意计算之前要确保两幅影像地理范围完全一致并且变化区域边界的配准误差不要超过一个像素。本文还有配套的精品资源点击获取
