企业官网建设流程全解析

简介：这份资源面向遥感影像处理与变化检测方向的开发者、测绘及地理信息专业学生，提供一套基于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_shape=None): """ 读取两期影像并做基本对齐检查。 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]), resampling=Resampling.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, dtype=np.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"第 {b+1} 波段有效像元过少，检查 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-score，PCA 才能公平地对待每个波段。

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_components=3): """ 对差值立方体做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% 以上，PC1+PC2 能到 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( dtype=rasterio.uint8, count=1, compress='lzw' ) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(intensity_uint8, 1) return intensity

参数说明：profile来自前面读取的影像，包含了 transform、crs、width、height 等地理信息。compress='lzw'是无损压缩，减小文件体积。保存成 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, method='otsu', percentile=95): """ 对变化强度图做阈值分割。 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.percentile，percentile=95 表示只保留强度最高的 5% 像元。返回的binary是 0/1 图，1 代表变化。

4.2 形态学去噪与最小图斑过滤

二值图里会有大量孤立的单像元噪声，以及变化区域内部的空洞。常见做法是开运算（先腐蚀后膨胀）去掉小噪点，闭运算（先膨胀后腐蚀）填补空洞。然后用连通域分析，去掉面积小于最小图斑阈值的斑块。最小图斑阈值根据你的制图规范来，比如 0.5 公顷对应多少像元，要按分辨率换算。

from scipy import ndimage def clean_binary(binary, min_pixels=9, open_size=3, close_size=3): """形态学去噪 + 最小图斑过滤""" # 开运算去噪 opened = ndimage.binary_opening(binary, structure=np.ones((open_size, open_size))) # 闭运算填洞 closed = ndimage.binary_closing(opened, structure=np.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_pixels=9表示至少 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, mask=mask, transform=transform) ) gdf = gpd.GeoDataFrame.from_features(results, crs=crs) # 只保留 value=1 的变化图斑 gdf = gdf[gdf['value'] == 1] gdf.to_file(out_shp, encoding='utf-8') return gdf

逻辑说明：shapes生成器逐个返回几何和值，mask=mask确保只处理非零区域。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（变化向量分析）代替 PCA，CVA 对辐射差异的鲁棒性稍好。

5.3 现象：Otsu 阈值分割后变化像元占比超过 30%

原因：变化强度图直方图不是双峰，Otsu 假设失效。或者影像里有大片云、云影、山体阴影，这些在差值里表现为极端值，拉偏了阈值。

解决：先做云掩膜，用 QA 波段或 Fmask 算法去掉云和云影。然后改用分位数阈值，从 95% 开始试，逐步调到变化比例符合实际。也可以对强度图做对数变换，压缩长尾后再用 Otsu。

5.4 现象：内存溢出，程序在计算协方差矩阵时崩溃

原因：把差值立方体重塑成(N, bands)后，误用了X.T @ X计算(N, N)矩阵。N 是像元总数，一幅 10000x10000 的影像 N=1 亿，(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, window=512, stride=256): """滑动窗口PCA，返回全局变化强度图""" bands, H, W = arr1.shape intensity = np.zeros((H, W), dtype=np.float32) weight = np.zeros((H, W), dtype=np.float32) for y in range(0, H - window + 1, stride): for x in range(0, W - window + 1, stride): sub1 = arr1[:, y:y+window, x:x+window] sub2 = arr2[:, y:y+window, x:x+window] diff = build_diff_cube(sub1, sub2) comps, _, _ = pca_on_diff(diff, n_components=1) local_intensity = np.abs(comps[0]) # 汉宁窗加权，减少边界突变 wy = np.hanning(window)[:, None] wx = np.hanning(window)[None, :] w = wy * wx intensity[y:y+window, x:x+window] += local_intensity * w weight[y:y+window, x:x+window] += w intensity = intensity / (weight + 1e-8) return intensity

参数说明：window=512是经验值，太小则协方差估计不稳定，太大则失去局部性。stride=256是 50% 重叠。汉宁窗让窗口中心权重高、边缘权重低，融合后过渡自然。这个方法的代价是计算量成倍增加，一幅 10000x10000 的影像大约要跑 1500 个窗口，每个窗口做一次 PCA，用多进程可以加速。

精度验证方面，如果没有地面真值，我一般用两种方式交叉验证：一是用高分辨率影像（比如 Google Earth 历史影像）目视抽查 100 个随机点，统计漏检和误检；二是用变化前后的 NDVI 差值做独立参考，看 PCA 变化图斑和 NDVI 显著下降区域的重合度。如果重合度低于 70%，说明 PCA 结果不可靠，要回去检查辐射归一化和阈值。我自己的习惯是，任何无监督变化检测结果，在交付前必须做至少 50 个点的目视抽查，宁可多花半天，也不要让假图斑流到下游。希望帮到你。

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