简介这份资源面向遥感影像处理与变化检测方向的开发者、测绘及地理信息专业学生提供一套基于Python的PCA变化检测算法实现。算法结合sklearn与opencv对两期同尺寸遥感影像做主成分分析提取变化区域并支持大影像分块处理最终可将变化图斑转为矢量输出。代码中还加入基于图像处理的图斑过滤逻辑可按面积过小或长宽比过大等条件自定义剔除减少碎斑干扰。压缩包共3个文件均为py脚本整体约4KB分别承担主流程调度、PCA核心检测与矢量文件读写等职责结构精简便于二次修改与集成。目前已有1292人学习下载适合希望快速上手遥感变化检测、理解PCA降维思路并落地矢量成果的读者参考也可作为相关项目或课程实验的基础代码。1. 从两期遥感影像到变化图斑PCA 变化检测到底在做什么手头有两期同一区域的遥感影像一期是去年 6 月一期是今年 6 月领导要你圈出这一年里新增的建设用地、被采伐的林地、扩张的水体。人工比对两幅大图眼睛看花也容易漏。这时候很多人会搜「Python PCA遥感影像变化检测算法代码」想找一段能直接跑的脚本。PCA 变化检测的核心思路其实不复杂把两期影像对应波段的差值堆成一个多波段差值立方体用主成分分析PCA把这个立方体压缩到少数几个主成分上变化信息通常集中在前一两个主成分里再对主成分做阈值分割或聚类就得到变化图斑。它不需要训练样本属于无监督变化检测适合没有标注数据、只想快速拿到变化候选区的场景。这篇文章面向会用 Python 处理栅格数据、但没系统做过变化检测的从业者从数据准备、差值立方体构建、PCA 实现、阈值选取一路写到踩坑排查代码可以直接抄。2. 数据准备与差值立方体把两期影像对齐成可比较的输入2.1 为什么必须先做辐射归一化和几何配准PCA 变化检测对输入极其敏感。如果两期影像的成像季节、太阳高度角、大气条件差异大那么差值影像里大部分「变化」其实是辐射差异不是地物变化。常见做法是先做相对辐射归一化选一期作为参考用另一期的伪不变特征PIF比如深水体、裸岩、大片成熟林地做线性回归把两期拉到同一辐射尺度。几何配准同样关键配准误差超过一个像元边缘就会产生大量假变化。我一般要求两期影像的均方根误差控制在 0.5 像元以内配准用 ENVI 或 GDAL 的自动配准都行但配准后一定要目视检查几个明显地物点。提示如果两期影像来自不同传感器比如 Landsat 8 和 Sentinel-2波段设置和空间分辨率都不同需要先重采样到同一分辨率并选取对应的重叠波段否则差值没有物理意义。2.2 用 rasterio 读取并对齐两期影像下面这段代码完成三件事读取两期影像、检查波段数和尺寸是否一致、把两期影像裁剪到同一范围。这里用 rasterio 而不是 gdal是因为 rasterio 的 Python 接口更干净读进来的就是 numpy 数组方便后续做 PCA。import rasterio import numpy as np from rasterio.enums import Resampling def read_and_align(path_t1, path_t2, bands_t1, bands_t2, out_shapeNone): 读取两期影像并做基本对齐检查。 bands_t1/bands_t2: 要使用的波段索引列表从1开始 out_shape: 可选统一重采样到的 (height, width) with rasterio.open(path_t1) as src1: # 按指定波段读取得到 (bands, H, W) arr1 src1.read(bands_t1) profile src1.profile.copy() transform1 src1.transform crs1 src1.crs with rasterio.open(path_t2) as src2: arr2 src2.read(bands_t2) transform2 src2.transform crs2 src2.crs # 检查坐标系是否一致 if crs1 ! crs2: raise ValueError(两期影像坐标系不一致请先统一投影) # 检查尺寸不一致则重采样到同一尺寸 if arr1.shape ! arr2.shape: if out_shape is None: out_shape arr1.shape[1:] # 以T1的尺寸为准 with rasterio.open(path_t2) as src2: arr2 src2.read( bands_t2, out_shape(len(bands_t2), out_shape[0], out_shape[1]), resamplingResampling.bilinear ) # 转为 float32避免后续差值溢出 arr1 arr1.astype(np.float32) arr2 arr2.astype(np.float32) return arr1, arr2, profile, transform1逻辑说明src.read(bands_t1)返回的是(波段数, 高, 宽)的三维数组这是 rasterio 的默认顺序和 GDAL 一致。重采样用双线性插值对连续型光谱波段合适如果是分类图应该用最近邻。参数bands_t1和bands_t2要选物理意义对应的波段比如 Landsat 8 的蓝、绿、红、近红外对应 Sentinel-2 的 B2、B3、B4、B8不能随便凑数。2.3 构建差值立方体与归一化差值立方体就是把两期影像逐波段相减得到(波段数, H, W)的数组。但直接相减有个问题不同波段的数值范围差异大近红外波段反射率高差值绝对值也大PCA 会被大数值波段主导。所以差值后要做逐波段标准化让每个波段均值为 0、标准差为 1。def build_diff_cube(arr1, arr2): 构建差值立方体并逐波段标准化 # 逐波段差值 diff arr2 - arr1 # shape: (bands, H, W) # 逐波段标准化减去均值除以标准差 diff_norm np.zeros_like(diff, dtypenp.float32) for b in range(diff.shape[0]): band diff[b] # 只统计有效像元排除 NoData这里假设 NoData 为 0 或 NaN valid np.isfinite(band) (band ! 0) if valid.sum() 100: raise ValueError(f第 {b1} 波段有效像元过少检查 NoData 设置) mu band[valid].mean() sigma band[valid].std() diff_norm[b] (band - mu) / (sigma 1e-8) return diff_norm参数说明valid掩膜排除了 0 值和 NaN如果你的影像 NoData 是 -9999要把条件改成band ! -9999。1e-8是防止标准差为 0 的兜底。标准化之后每个波段的差值都变成无量纲的 z-scorePCA 才能公平地对待每个波段。3. PCA 降维与变化分量提取从多波段差值到变化强度图3.1 PCA 在差值立方体上的数学过程把标准化后的差值立方体(bands, H, W)重塑成(bands, N)N 是像元总数。然后计算波段间的协方差矩阵C (1/(N-1)) * X X.T形状是(bands, bands)。对 C 做特征分解得到特征值从大到小排列的特征向量。把原始数据投影到前 k 个特征向量上就得到 k 个主成分。变化信息通常集中在第一主成分PC1因为变化像元在所有波段上都有同向的差值这种共同模式方差最大。第二、第三主成分可能对应不同地物的光谱变化方向也可以纳入分析。这里有个容易翻车的地方N 通常是几百万甚至上千万直接算X X.T是(bands, N) (N, bands)结果只有(bands, bands)计算量可控。但如果你反过来算X.T X那就是(N, N)的矩阵内存直接爆掉。我见过有人这么写然后 32G 内存的机器卡死。3.2 用 numpy 实现 PCA 并提取前三个主成分def pca_on_diff(diff_norm, n_components3): 对差值立方体做PCA。 diff_norm: (bands, H, W) 标准化后的差值 返回: components (n_components, H, W), eigenvalues, eigenvectors bands, H, W diff_norm.shape # 重塑为 (bands, N) X diff_norm.reshape(bands, -1) # 计算协方差矩阵 (bands, bands) # 注意X已经逐波段标准化均值接近0 C np.cov(X) # 特征分解 eigenvalues, eigenvectors np.linalg.eigh(C) # eigh返回升序反转成降序 idx np.argsort(eigenvalues)[::-1] eigenvalues eigenvalues[idx] eigenvectors eigenvectors[:, idx] # 投影到前n_components个主成分 # eigenvectors[:, :k] shape: (bands, k) # X.T shape: (N, bands) proj eigenvectors[:, :n_components].T X # (k, N) components proj.reshape(n_components, H, W) return components, eigenvalues, eigenvectors逻辑说明np.cov(X)默认按行计算协方差X 的每一行是一个波段所以得到的是波段间协方差矩阵这正是我们需要的。np.linalg.eigh用于对称矩阵比eig更稳定更快。投影时用eigenvectors[:, :k].T X得到(k, N)再 reshape 回影像尺寸。eigenvalues可以用来判断前几个主成分解释了多少方差一般 PC1 能解释 70% 以上PC1PC2 能到 85% 以上如果 PC1 占比不到 50%说明差值立方体里噪声太大要回去检查辐射归一化。3.3 变化强度图的生成与可视化PC1 的绝对值越大说明该像元在两期之间的光谱变化越剧烈。但 PC1 有正有负正负代表变化方向不同比如从植被变裸土和从裸土变植被在 PC1 上符号相反。做变化检测时通常取 PC1 的绝对值作为变化强度或者对 PC1 做平方和。下面代码生成变化强度图并保存为 GeoTIFF。def save_change_intensity(components, profile, out_path): 把PC1的绝对值保存为变化强度图 pc1 components[0] intensity np.abs(pc1) # 归一化到 0-255 便于可视化 intensity_norm (intensity - intensity.min()) / (intensity.max() - intensity.min() 1e-8) intensity_uint8 (intensity_norm * 255).astype(np.uint8) profile.update( dtyperasterio.uint8, count1, compresslzw ) with rasterio.open(out_path, w, **profile) as dst: dst.write(intensity_uint8, 1) return intensity参数说明profile来自前面读取的影像包含了 transform、crs、width、height 等地理信息。compresslzw是无损压缩减小文件体积。保存成 uint8 是为了在 QGIS 或 ArcGIS 里直接看如果要做后续阈值分割应该保留 float32 的intensity数组。4. 阈值分割与变化图斑后处理从强度图到可用矢量4.1 阈值选取Otsu 与分位数法的取舍变化强度图是一张连续灰度图要变成二值变化/未变化图必须选阈值。最常用的是 Otsu 大津法它假设图像由前景和背景两类组成自动找类间方差最大的阈值。但遥感变化强度图的直方图往往不是双峰而是长尾分布Otsu 容易把阈值选偏低导致大量假变化。我的经验是先用 Otsu 跑一版看效果如果假变化太多改用分位数法比如取 95% 或 97% 分位数作为阈值只保留变化最剧烈的像元。分位数的选择取决于你研究区实际变化比例城市扩张区可能 5% 到 10%自然保护区可能不到 1%。from skimage.filters import threshold_otsu def threshold_change(intensity, methodotsu, percentile95): 对变化强度图做阈值分割。 method: otsu 或 percentile valid np.isfinite(intensity) vals intensity[valid] if method otsu: thresh threshold_otsu(vals) elif method percentile: thresh np.percentile(vals, percentile) else: raise ValueError(method 只支持 otsu 或 percentile) binary (intensity thresh).astype(np.uint8) return binary, thresh逻辑说明threshold_otsu来自 scikit-image输入是一维有效值数组。分位数法用np.percentilepercentile95 表示只保留强度最高的 5% 像元。返回的binary是 0/1 图1 代表变化。4.2 形态学去噪与最小图斑过滤二值图里会有大量孤立的单像元噪声以及变化区域内部的空洞。常见做法是开运算先腐蚀后膨胀去掉小噪点闭运算先膨胀后腐蚀填补空洞。然后用连通域分析去掉面积小于最小图斑阈值的斑块。最小图斑阈值根据你的制图规范来比如 0.5 公顷对应多少像元要按分辨率换算。from scipy import ndimage def clean_binary(binary, min_pixels9, open_size3, close_size3): 形态学去噪 最小图斑过滤 # 开运算去噪 opened ndimage.binary_opening(binary, structurenp.ones((open_size, open_size))) # 闭运算填洞 closed ndimage.binary_closing(opened, structurenp.ones((close_size, close_size))) # 连通域标记 labeled, num ndimage.label(closed) # 统计每个连通域面积 sizes ndimage.sum(closed, labeled, range(1, num 1)) # 保留面积大于阈值的 keep np.zeros_like(closed) for i, s in enumerate(sizes): if s min_pixels: keep[labeled i 1] 1 return keep.astype(np.uint8)参数说明open_size和close_size一般取 3对应 3x3 结构元。min_pixels9表示至少 9 个像元如果分辨率是 10 米9 个像元约 0.09 公顷偏小实际制图可能要调到 50 或 100。ndimage.label默认用 4 连通如果变化区域是斜向条带可以改成 8 连通。4.3 矢量化输出与属性统计最后把二值栅格转成矢量多边形方便在 GIS 里编辑和统计。用 rasterio.features.shapes 可以做到。from rasterio.features import shapes import geopandas as gpd from shapely.geometry import shape def binary_to_vector(binary, transform, crs, out_shp): 二值图转矢量多边形 mask binary.astype(np.uint8) results ( {properties: {value: v}, geometry: shape(s)} for s, v in shapes(mask, maskmask, transformtransform) ) gdf gpd.GeoDataFrame.from_features(results, crscrs) # 只保留 value1 的变化图斑 gdf gdf[gdf[value] 1] gdf.to_file(out_shp, encodingutf-8) return gdf逻辑说明shapes生成器逐个返回几何和值maskmask确保只处理非零区域。gpd.GeoDataFrame.from_features把结果转成 GeoDataFrame指定 crs 保证坐标系正确。输出 shapefile 时用 utf-8 编码避免中文属性乱码。5. 避坑与排查PCA 变化检测最常见的五个翻车点5.1 现象变化图斑沿影像边缘和道路大量出现原因两期影像几何配准误差大或者重采样时边缘像元被拉伸。道路、田埂这类线状地物对配准误差最敏感一个像元的偏移就能产生整条假变化带。解决回到配准步骤用至少 20 个均匀分布的控制点重新配准检查 RMSE 是否小于 0.5 像元。如果配准没问题检查重采样方法连续波段用双线性分类图用最近邻。还可以在差值前对两期影像做 3x3 均值滤波牺牲一点空间细节换取配准鲁棒性。5.2 现象PC1 解释方差不到 40%变化强度图一片模糊原因差值立方体里噪声占主导可能是辐射归一化没做或者两期影像季节差异太大。比如一期是雨季、一期是旱季水体面积变化剧烈PCA 会把水体变化当成主要模式掩盖真正的地物变化。解决先做相对辐射归一化用 PIF 做线性回归。如果季节差异无法避免考虑只用对季节不敏感的波段比如短波红外和近红外或者改用 CVA变化向量分析代替 PCACVA 对辐射差异的鲁棒性稍好。5.3 现象Otsu 阈值分割后变化像元占比超过 30%原因变化强度图直方图不是双峰Otsu 假设失效。或者影像里有大片云、云影、山体阴影这些在差值里表现为极端值拉偏了阈值。解决先做云掩膜用 QA 波段或 Fmask 算法去掉云和云影。然后改用分位数阈值从 95% 开始试逐步调到变化比例符合实际。也可以对强度图做对数变换压缩长尾后再用 Otsu。5.4 现象内存溢出程序在计算协方差矩阵时崩溃原因把差值立方体重塑成(N, bands)后误用了X.T X计算(N, N)矩阵。N 是像元总数一幅 10000x10000 的影像 N1 亿(1亿, 1亿)的矩阵需要 4e16 字节任何机器都扛不住。解决始终用np.cov(X)它内部计算的是(bands, bands)矩阵计算量只和波段数有关。如果波段数很多比如高光谱几百个波段可以先用随机采样取一部分像元估计协方差再投影全图。5.5 现象变化图斑破碎大量 1-2 像元的碎斑原因PCA 对噪声敏感差值立方体里的随机噪声在 PC1 上表现为孤立高值。阈值分割后这些噪声变成碎斑。解决在阈值分割前对 PC1 做高斯滤波或中值滤波平滑掉高频噪声。阈值分割后做形态学开运算和最小图斑过滤。如果碎斑仍然多考虑在 PCA 之前对差值立方体做 5x5 中值滤波但要注意这会模糊小面积变化。6. 进阶技巧用滑动窗口 PCA 捕捉局部变化并验证精度全局 PCA 有一个固有缺陷它提取的是整幅影像的全局变化模式如果研究区里同时存在城市扩张和森林砍伐两种变化它们的差值方向可能不同全局 PC1 只能捕捉其中方差最大的那个另一种变化会被压制。我一般会在大区域上改用滑动窗口 PCA把影像切成 512x512 的窗口每个窗口独立做 PCA取窗口内 PC1 的绝对值作为局部变化强度。这样不同区域的变化模式不会被互相掩盖。窗口之间要有 50% 重叠避免边界效应最后用加权平均融合重叠区。def sliding_window_pca(arr1, arr2, window512, stride256): 滑动窗口PCA返回全局变化强度图 bands, H, W arr1.shape intensity np.zeros((H, W), dtypenp.float32) weight np.zeros((H, W), dtypenp.float32) for y in range(0, H - window 1, stride): for x in range(0, W - window 1, stride): sub1 arr1[:, y:ywindow, x:xwindow] sub2 arr2[:, y:ywindow, x:xwindow] diff build_diff_cube(sub1, sub2) comps, _, _ pca_on_diff(diff, n_components1) local_intensity np.abs(comps[0]) # 汉宁窗加权减少边界突变 wy np.hanning(window)[:, None] wx np.hanning(window)[None, :] w wy * wx intensity[y:ywindow, x:xwindow] local_intensity * w weight[y:ywindow, x:xwindow] w intensity intensity / (weight 1e-8) return intensity参数说明window512是经验值太小则协方差估计不稳定太大则失去局部性。stride256是 50% 重叠。汉宁窗让窗口中心权重高、边缘权重低融合后过渡自然。这个方法的代价是计算量成倍增加一幅 10000x10000 的影像大约要跑 1500 个窗口每个窗口做一次 PCA用多进程可以加速。精度验证方面如果没有地面真值我一般用两种方式交叉验证一是用高分辨率影像比如 Google Earth 历史影像目视抽查 100 个随机点统计漏检和误检二是用变化前后的 NDVI 差值做独立参考看 PCA 变化图斑和 NDVI 显著下降区域的重合度。如果重合度低于 70%说明 PCA 结果不可靠要回去检查辐射归一化和阈值。我自己的习惯是任何无监督变化检测结果在交付前必须做至少 50 个点的目视抽查宁可多花半天也不要让假图斑流到下游。希望帮到你。本文还有配套的精品资源点击获取