简介该资源为遥感变化检测方向的实用代码包面向地理信息、计算机视觉及遥感初学者聚焦PCA主成分分析与K-means聚类两类方法的结合应用。压缩包共14个文件、大小约59KB其中algorithm.py、k_means.py、util.py与main.py构成完整执行流程方便直接运行或二次修改PCAKmeans_burn.png、PCAKmeans_forest.png、PCAKmeans_conifer.png可直观查看不同地物的变化检测结果。资源已在CSDN被254人浏览学习适合希望快速理解PCA降维、K-means聚类及变化检测实施步骤的人群。通过阅读源码与示例图像可掌握从数据预处理、特征提取到聚类分析、变化区域识别的整体思路也能了解初始聚类中心选择对结果的影响为后续环境监测、植被变化分析等项目打下基础。1. 变化检测不止目视判读这套 PCA-Kmeans 脚本为什么值得下手里两期同一地区的遥感影像想找植被破坏或火烧迹地靠人工目视判读两张图来回切圈完一个区域至少半小时还容易漏掉藏在阴影里的小块变化。用程序跑一遍变化检测几分钟出一张变化分布图再回头复核重点区域效率完全不是一个量级。这份 PCAKmeans 资源就是干这个的核心是主成分分析加 K-means 聚类用 PCA 把多波段影像压成少数几个信息最集中的主成分再交给 K-means 聚成不同的类通过比较两时相的聚类结果定位变化区域。压缩包里有算法实现、主流程脚本和三张样例结果图。适合刚开始做遥感变化检测的从业者也适合想把手写 PCA 和 K-means 流程跑通的算法同学。2. 方法原理PCA 把噪声滤掉K-means 把变化找出来2.1 PCA 在变化检测里到底做了什么降维与噪声抑制PCA 在遥感图像处理里不是新鲜事人脸识别里的特征脸就是 PCA 降维后的基向量遥感里的主成分图也是类似思路只是基向量从人脸像素变成了波段组合方向。拿哨兵二号或者 Landsat 的多光谱影像来说十几个波段之间高度相关——近红外和红边波段对植被的响应是同步变化的直接把十几张波段图逐像素做差噪声会淹没真实变化信号。PCA 做的事情是把这些相关波段重新组合成一组互不相关的通道第一主成分集中了绝大部分方差对应地表反射的整体亮度结构第二、第三主成分往往对应植被、水体、裸土这类地物特征的差异。噪声因为方差小被自然排挤到后面的主成分里丢弃后不会影响变化信息的提取。实现步骤上常见做法是先做标准化把每个波段的均值归零、方差归一到单位方差否则量纲差异会让高值波段主导整个协方差矩阵。然后算协方差矩阵做特征值分解特征值大小代表每个主成分保留的信息量特征向量就是投影方向。做主成分选择时我一般看累计方差贡献率阈值取 95% 或 99%参数常用取值建议累计方差贡献率0.95 ~ 0.99数据质量好取 0.95噪声大取 0.99主成分个数3 ~ 6超过 6 个通常已经进入噪声带是否标准化是多光谱影像必须标准化2.2 K-means 的价值无监督聚类识别变化区域K-means 在这个流程里扮演的角色是把变化信息切成可解释的区域。对两时相图像做主成分分析之后把同一位置的主成分值相减得到的差值特征就代表了该位置在两个时间点之间的光谱变化。理论上未变化的区域差值接近零变化区域差值明显偏离零。但明显到什么程度靠人工设阈值容易翻车因为不同地物类型的背景噪声水平不一样。K-means 的聪明之处在于用一种无监督的方式自动划分先在差值特征空间里随机放 k 个中心点每个像素归到距离最近的中心然后重新计算中心位置迭代直到中心不再明显移动。最终聚出来的簇中一类通常对应未变化一类或几类对应不同类型的变化。这个过程完全不需要标注数据这对遥感应用很友好——大部分历史影像根本没有对应的真值标签。实际使用中k 值选 2 还是选 3 以上取决于你要不要区分变化方向。如果只是找植被变没了这种单一变化k2 就够如果想同时分出植被变成裸土和裸土变成植被k3 或 4 更合适。但要注意K-means 聚类出的簇不一定严格对应语义类别簇的物理意义要靠后续的地物属性分析来判断。2.3 为什么不直接用差值法或监督分类三个方案的对比做变化检测最朴素的做法是波段差值——直接把两时相的对应波段相减超过阈值判定为变化。这个思路简单但阈值选择很看经验而且对于不同波段量纲不一致的问题缺少应对手段。变化向量分析 CVA 是差值法的升级版把多个波段的变化量合成一个向量用向量模长判断变化强度但同样依赖阈值。监督分类则是另一条路训练一个分类器识别地物类型再比较两时相的分类结果。这在有高质量训练样本的时候非常准但样本标注成本极高而且不同季节影像的分类器通用性差。PCA-Kmeans 的组合恰好卡在中间——不需要标注不需要调阈值PCA 做好了噪声抑制K-means 自动划分变化簇计算量也远小于深度学习方法。下表把几个方案的特点摆在了一起方法是否需要标注核心难点适用场景波段差值法否阈值难以确定单波段、变化明显且均匀变化向量分析 CVA否多维变化阈值难标定多波段变化强度评估监督分类后比较是样本标注成本高、跨期鲁棒性差有长期真值数据的区域PCA K-means否k 值选择与初始中心影响快速普查、无真值区域这套资源采用的就是最后一种方案。main.py 里把两时相影像拼接起来做 PCA再对主成分差值做 K-means属于遥感变化检测里验证过的经典流程尤其适用于火烧迹地、砍伐、城市扩张这类有明确光谱响应的地表变化。3. 代码拆解从 util.py 到 main.py 的完整流转3.1 文件清点哪几个文件是核心哪几个可以删解压之后先别急着跑把文件角色认清楚。algorithm.py 和 k_means.py 是算法的实现核心util.py 是辅助函数main.py 是流程入口。三张结果图 PCAKmeans_burn.png、PCAKmeans_forest.png、PCAKmeans_conifer.png 分别是火烧迹地、森林、针叶林场景下的输出示例跑代码前先看这几张图能建立直观印象。剩下的 .idea 目录和pycache目录是开发环境生成的临时文件。.idea 里是 PyCharm 的项目配置pycache里是 Python 编译产生的字节码缓存与算法逻辑无关删掉不影响运行。文件角色核心内容util.py工具层图像读取、数据归一化等辅助处理algorithm.py算法层PCA 主成分分析实现k_means.py算法层K-means 聚类实现main.py入口层串联完整变化检测流程PCAKmeans_burn.png示例输出火烧迹地检测结果图PCAKmeans_forest.png示例输出森林区域变化检测结果图PCAKmeans_conifer.png示例输出针叶林区域变化检测结果图3.2 algorithm.pyPCA 主成分计算的实现逻辑从文件划分推断algorithm.py 里大概率封装了一个 PCA 类或者函数输入是二维数组形状是像素数乘波段数输出是降维后的主成分得分图。手写 PCA 的实现思路一般如下核心步骤是标准化、协方差矩阵、特征值分解和投影import numpy as np def pca(data, n_componentsNone, variance_ratio0.95): # data: (n_samples, n_features)每行一个像素每列一个波段 # 1. 标准化避免量纲差异影响特征值分解 mean np.mean(data, axis0) std np.std(data, axis0) std[std 0] 1e-12 data_norm (data - mean) / std # 2. 协方差矩阵与特征值分解 cov np.cov(data_norm.T) eig_vals, eig_vecs np.linalg.eigh(cov) # 3. 特征值升序排列翻转成降序 idx np.argsort(eig_vals)[::-1] eig_vals eig_vals[idx] eig_vecs eig_vecs[:, idx] # 4. 按累计方差贡献率自动确定主成分个数 if n_components is None: total np.sum(eig_vals) cum_ratio np.cumsum(eig_vals) / total n_components int(np.searchsorted(cum_ratio, variance_ratio) 1) # 5. 投影到主成分空间 proj data_norm eig_vecs[:, :n_components] return proj, eig_vals[:n_components], eig_vecs[:, :n_components]逻辑说明第 2 步用np.linalg.eigh而不是eig因为协方差矩阵是实对称矩阵eigh专门处理对称矩阵数值稳定性更好第 3 步翻转特征值顺序是因为eigh默认按升序返回而我们想要信息量最大的主成分排在前面第 4 步的searchsorted找的是累计方差贡献率第一次超过 0.95 的位置。参数说明variance_ratio控制 PCA 保留多少信息影像噪声大时建议调到 0.99但主成分个数会增多后续聚类计算量也变大需要做取舍。n_components传固定值可以跳过自动计算比如直接传 4前 4 个主成分在多数遥感场景下已经能覆盖绝大部分地表差异。3.3 k_means.py聚类核心循环怎么写k_means.py 里应该是一个标准的 K-means 实现输入主成分差值特征输出每个像素的簇标签。手写版核心代码通常长这样import numpy as np def kmeans(X, k, max_iter100, tol1e-4, random_state42): # X: (n_samples, n_features)主成分差值特征 rng np.random.default_rng(random_state) # 从样本里随机挑 k 个像素作为初始中心远优于全零初始化 centers X[rng.choice(X.shape[0], k, replaceFalse)] for epoch in range(max_iter): # 计算每个样本到各中心的欧氏距离 # X[:, None, :] 扩展维度后与 centers 做广播运算 dist np.linalg.norm(X[:, None, :] - centers[None, :, :], axis2) labels np.argmin(dist, axis1) # 更新中心每个簇的均值作为新中心 new_centers np.array([X[labels j].mean(axis0) for j in range(k)]) # 中心点位移小于阈值则判定收敛 shift np.linalg.norm(new_centers - centers) centers new_centers if shift tol: break return labels, centers逻辑说明距离矩阵的计算是 K-means 最耗时的部分X[:, None, :] - centers[None, :, :]的广播写法把像素数和中心数两个维度同时展开一行代码算出所有距离。迭代收敛条件是中心点移动距离小于tol这是最常用的停止标准。注意如果某个簇的样本数为 0X[labels j].mean(axis0)会产生空簇实际工程里要加一个空簇重初始化逻辑。参数说明k是目标簇数变化检测中通常取 2 或 3取 2 表示只分变化与未变化取 3 及以上表示进一步区分变化类型tol设太小会多跑几次迭代设太大容易在中心点还没稳定时提前收尾一般取 1e-4random_state固定之后结果可复现这篇资源里不同随机种子跑出的变化图会不一样这一点在下一章会展开说。3.4 main.py把两张时相串成一条变化检测流水线main.py 的角色是把 util、algorithm、k_means 串起来。整体流程是读取两时相影像、逐像素拉平成二维数组、拼接后统一做 PCA、计算主成分差值、K-means 聚类、重排成图像并输出。关键在 PCA 环节import numpy as np from util import read_image, save_image from algorithm import pca from k_means import kmeans # 读取两时相影像假设通道顺序一致、尺寸一致 t1 read_image(scene_2019.tif) # (h, w, bands) t2 read_image(scene_2021.tif) h, w, num_bands t1.shape X1 t1.reshape((-1, num_bands)) X2 t2.reshape((-1, num_bands)) # 关键两时相拼接后再做 PCA确保投影方向完全一致 X_all np.vstack([X1, X2]) P_all, _, _ pca(X_all, n_components4) # 切回两时相的主成分得分并恢复图像形状 P1 P_all[:X1.shape[0]].reshape((h, w, -1)) P2 P_all[X1.shape[0]:].reshape((h, w, -1)) # 逐通道做差得到主成分差值立方体 diff P1 - P2 feat diff.reshape((-1, diff.shape[2])) # 聚成变化与未变化两类 labels, centers kmeans(feat, k2, random_state2024) change_map labels.reshape((h, w)) save_image(change_map.png, change_map)逻辑说明这里最值得注意的写法是先把X1和X2垂直堆叠成X_all再统一做 PCA。两时相影像的波段本身就来自不同时刻的观测各自做 PCA 得到的主成分方向很可能不一致——第一主成分在前一时期对应亮度在后一时期可能对应湿度——这样直接做差毫无意义。拼接统一做 PCA 之后两边被投影到同一个坐标空间逐通道相减才有物理含义。聚类输入用的是差值特征而不是单时相的主成分图这样聚出来的簇直接对应变化模式。参数说明n_components4是经验值覆盖植被、水体、裸土和整体亮度四个维度k2对应变化与未变化二分类。跑通之后可以改成 k3 到 4在变化区域内部再做细分。random_state2024固定了随机初始中心保证每次运行结果一致避免复现时出现结果对不上的尴尬。4. 避坑记录辐射归一化、K 值与初始中心的五个坑4.1 两时相各自做 PCA差值图全是噪点现象按照分别对 t1 和 t2 做 PCA再对主成分图做差的思路跑出来的结果变化图上一片花白看不出集中连片的变化区域和资源附带的 PCAKmeans_burn.png 那种清晰的火烧迹地形态完全不一样。原因PCA 是一种无监督变换它找到的主成分方向是数据本身方差最大化的方向。两时相影像的亮度分布、大气条件不同各自求出来的特征向量正交基可能差别很大甚至某个主成分的方向符号是反的。这个时候逐通道相减得到的是两个不同坐标空间之间的投影残差而不是真实的光谱变化。解决把所有像素拼起来统一做 PCA。np.vstack([X1, X2])之后再做主成分分析保证两时相共用同一组特征向量。这个操作在第三章的 main.py 里已经是标准写法但如果自行改编代码务必检查 PCA 的输入是不是拼接后的数据。4.2 辐射归一化不做大气条件差异全被当成变化现象结果图里大片区域显示变化连稳定的裸山、水库周边都被标记为变化区域虚警率高到没法看甚至整个图像一半以上像素都处于变化簇里。原因遥感影像受大气散射、太阳高度角、传感器定标增益等因素影响两时相影像即使是同一传感器辐射值也有系统差异。这种差异在光谱空间里表现为整体偏移PCA 无法区分大气引起的辐射偏移和地表真实变化聚类时自然会把偏移当成一种变化模式。解决在进入 PCA 之前做辐射归一化。快速做法是选取两时相中稳定地物区域比如裸地、大型水体、道路用这些像素做线性回归拟合把 t2 的辐射值校正到 t1 的量纲上也可以用直方图匹配让 t2 的波段直方图逼近 t1。正规场景下推荐 IR-MAD 方法做相对辐射归一化但线性回归加直方图匹配在很多项目中已经够用。4.3 K 值拍脑袋决定分割粒度忽粗忽细现象K 取 2 的时候变化检测结果把大片火烧迹地和周边新修道路混在一起没法区分不同变化类型改成 K 取 5变化区域被打碎成许多小斑块中间还夹着零散噪声簇。原因K 值直接决定聚类粒度。取 2 时所有变化被迫归入一个簇不管光谱形态差异多大取太大时噪声本身也会聚成独立的簇因为噪声在特征空间里确实有聚集倾向。遥感变化区域在光谱特征空间里往往不是均匀分布靠猜容易翻车。解决跑肘部法则。对同一份差值特征依次取 k2 到 k8记录每个 k 对应的簇内距离总和 SSESSE 下降趋势出现明显拐点的位置就是合适的 k 值。另外一个直观办法是画两时相主成分前两维的散点图直接看差值分布在特征空间里形成了几团。4.4 初始中心敏感同一份数据两次跑出不同结果现象main.py 里不设random_state连跑三次变化检测得到三张略有差异的变化图有些边缘像素在变化和未变化之间反复跳动面积统计结果也不稳定。原因K-means 的初始中心是随机从样本里挑的不同的初值会把迭代过程导向不同的局部最优解。这个现象在特征维度高、簇之间有重叠时格外明显变化检测的场景恰恰如此——变化与未变化的边界在光谱空间里是渐变的而不是清晰间断的。解决一是固定随机种子用random_state或np.random.seed锁定二是换用 k-means 初始化策略让初始中心尽量分散减少随机性带来的影响三是多次聚类取 SSE 最小结果的方案。本文的手写版 kmeans 函数支持random_state参数直接固定一个值是最省事的做法。4.5 Python 版本兼容性旧代码跑到新环境直接报错现象import 时就报ModuleNotFoundError或者运行到 PCA 步骤时报AttributeError: module numpy has no attribute float项目目录下还有algorithm.cpython-36.pyc这样的字节码文件。原因从__pycache__里的文件后缀看这套代码最初是在 Python 3.6 环境下开发的。而 2024 年之后的 NumPy 版本已经移除了np.float、np.int这类旧类型别名很多早期代码跑在新环境里第一行就翻车。pyc字节码文件也不能跨 Python 版本使用3.6 编译的字节码在 3.12 下无法加载。解决把源码里的np.float改成float或np.float64np.int改成int或np.int64。直接删除__pycache__目录让 Python 在当前环境下重新编译不要依赖旧字节码。如果还报别的兼容性错误优先检查是否有np.bool、np.object这类旧别名统一替换即可。养成这个习惯之后这套代码在 Python 3.8 到 3.12 的环境里都能跑。5. 进阶用法从变了没有到怎么变的的验证与后处理5.1 后处理形态学滤波与连通域筛选去掉椒盐噪声K-means 直接输出的变化图通常是逐像素分类结果边缘区域免不了出现孤立小斑块和椒盐状噪声。这些零散像素在视觉上干扰判读在面积统计上也会造成较大误差。常见做法是先做形态学开运算消除小噪点再做连通域标记过滤掉面积低于阈值的碎块。只有面积大到有物理意义的区域才值得保留from scipy import ndimage # change_map: 布尔类型True 表示变化区域 opened ndimage.binary_opening(change_map, iterations1) # 标记连通域统计每个连通域的面积 label_im, num ndimage.label(opened) sizes ndimage.sum(opened, label_im, range(1, num 1)) # 面积小于 50 像素的连通域视为噪声置为 False clean opened.copy() for lbl, size in zip(range(1, num 1), sizes): if size 50: clean[label_im lbl] False # 对比去噪前后的变化区域面积 before change_map.sum() after clean.sum() print(f原始变化像素: {before}, 去噪后: {after})逻辑说明binary_opening先腐蚀后膨胀能把孤立的单点噪声抹掉同时保持大块区域形状基本不变ndimage.label给每个独立区域分配一个编号ndimage.sum统计各区域的像素数量面积阈值 50 是我在一般分辨率影像上的起点如果是高分影像阈值建议按实际像元尺寸换算成最小图斑面积。参数说明iterations1控制形态学运算的强度开一次就能去掉大部分椒盐噪声开两次以上会让变化区域边界明显收缩不建议过度使用。面积阈值 50 需要根据分辨率调整——0.5 米分辨率的无人机影像上50 像素可能只有 12.5 平方米Landsat 的 30 米分辨率影像上50 像素就是 4.5 公顷量级完全不同。5.2 验证精度混淆矩阵与 Kappa 系数是必须交的作业跑出变化图后还有一个容易偷懒但绕不开的环节精度验证。没有验证的变化检测结果只能算初步产品提交给别人之前至少要做一个混淆矩阵和 Kappa 系数。做法是先在图上随机撒点然后回到两时相原始影像的对应位置目视判断这个点到底变没变得到一组真值标签再和聚类结果对比from sklearn.metrics import confusion_matrix, cohen_kappa_score import numpy as np # y_true: 目视判读真值0 未变化 / 1 变化 # y_pred: 聚类结果0 未变化 / 1 变化 cm confusion_matrix(y_true, y_pred) kappa cohen_kappa_score(y_true, y_pred) # 总体精度取混淆矩阵对角线之和占总样本的比例 total_acc np.trace(cm) / np.sum(cm) print(f混淆矩阵:\n{cm}) print(f总体精度: {total_acc:.3f}, Kappa: {kappa:.3f})逻辑说明混淆矩阵能暴露系统偏差——比如聚类结果把大量真实变化区域漏检了还是把大片未变化区域误报成了变化从矩阵的列和行分布一眼就能看出来。Kappa 系数则校正了随机一致性的影响通常高于 0.8 认为变化检测结果可用0.6 到 0.8 之间需要回头检查分类器和特征提取的合理性。参数说明采样点数根据研究区大小决定一般至少 300 个样本其中变化与未变化区域按面积比例分配每类最少要有 50 个样本支撑统计意义。目视判读时最好使用多波段合成图比如近红外、红、绿波段合成比真彩色合成对植被变化的判读更敏感。5.3 延伸一点从两时相到时间序列以及新的方向PCAKmeans 这套流程做到位之后往前的路大致有两条。一是把两时相扩展成多时相时间序列比如用 Sentinel-2 每个月一期影像做长时间序列分析PCA 主成分得分随时间的变化曲线比单次差值更能反映植被恢复、城市扩张的连续过程。二是留意一下当前 CV 领域红火的开放词汇变化检测思路——用视觉语言模型把道路扩展裸土增加这类自然语言描述直接变成检测目标。这套基于 PCA 和 K-means 的传统基线虽然朴实但数据分布摸得透后续换新方法时也清楚要对比的参照物是什么。说到底变化检测的本质不是魔法而是把数据维度压低、把变化信号和噪声分开、再把结果解释成人类可理解的语言。从 PCA 和 K-means 入手能少走很多弯路。做那一次火烧迹地实验时我因为偷懒跳过辐射归一化直接跑 PCA结果把整片山体阴影都当成变化区域标记了出来白白浪费了一个晚上复查。从那以后我每次做变化检测都强制先出一张差值统计图扫一眼数据分布确认没有明显辐射偏移再决定要不要上 PCA 和 K-means。这个习惯帮我过滤掉了至少一半的无效实验。希望帮到你。本文还有配套的精品资源点击获取