企业官网建设流程全解析

简介：基于PCA的遥感影像变化检测工具包，适合遥感、GIS方向的开发者或科研人员快速对比两期影像差异并做自动化后处理。核心代码使用sklearn完成主成分分析与差值计算，借助OpenCV实现图像处理，针对大尺寸影像设计了内存友好型方案；检测结果可输出为灰度变化图，并通过内置shape_file_io模块自动转为矢量面文件，便于在GIS平台中进一步分析。包内共7个文件，以4个Python脚本为主，分别承担核心算法、主流程控制、矢量读写等职责，同时包含依赖清单txt等配套内容，压缩包仅9KB，轻量、易部署。已有44人学习下载。对于需要批量处理或二次开发的用户，可自由调节最小面积阈值、长宽比等图斑过滤参数以剔除噪声，并借助周边子模块拓展更多变化检测流程；注意两期影像需行列数一致，存在几何差异时应预先配准或重采样。

1. 别急着上深度学习：一个Python实现的PCA遥感影像变化检测工具包，专治“只有两期影像没有样本”的落地场景

做遥感变化检测，很多从业者一上来就想堆模型。可真的到了生产项目里，样本往往只有零星几个目视解译的图斑，深度学习根本喂不饱。这时候最踏实的反而是老牌方法：对两期影像做差，再把差值波段堆叠起来做PCA，取第一主成分当变化强度图，阈值分割后矢量化、过滤碎斑，一条链路下来，成果就能直接交付。这个Python实现的PCA遥感影像变化检测工具包就是这么一套东西，核心流程包含矢量化与图斑过滤，能把“像元级检测结果”变成“可入库的矢量图斑”。文章会拆开这条链路的每一步，给出可复现的代码和参数设置。适合手里有双时相影像、想在一周内出一版成果图斑的从业者，也适合想搞懂PCA到底怎么用进遥感项目的学生。

2. PCA遥感变化检测的原理与数据准备：主成分到底在替你看什么

2.1 双时相差异堆叠：PCA能检测变化的数学逻辑

先抛一个容易被忽略的事实：PCA本身并不认识“变化”这个词，它能干的只有一件事——把一组高相关波段压缩成若干个互不相关的主成分，方差越大，代表波段间共同变化越剧烈。遥感变化检测正是利用这一点。

常见做法不是拿两期影像直接做PCA，而是先做差。设两期影像各有B个波段，逐波段相减得到B个差值波段，把这B个差值波段堆叠成一个三维数组。像元值在差值波段上的绝对值越大，说明该位置在两个时间点的光谱特征差异越明显；PCA对这个数组降维时，第一主成分会顺着“所有差值波段都在变”的方向投影，也就是变化最集中的方向。这个方向上的得分图，就是我们要的变化强度图。

为什么用PCA而不是简单地对差值波段求均值或只选某一个差值波段？因为地表变化往往不是单一的——比如植被退化在红光波段和近红外波段都有响应，但响应强度不同；水体变化可能在短波红外上最明显。单个差值波段容易漏掉只在特定波段表现突出的变化，而均衡取均值又会把“所有波段都很强”和“某一个波段特别强”混为一谈。PCA利用协方差结构重新加权，让投影方向自动适应数据的实际响应模式，这比人工拍脑袋挑波段要稳得多。

这套逻辑的一个隐含前提是：变化信号在多个波段里具有相关性，而噪声通常不具备这种相关性。PCA相当于做了一次信号增强、噪声抑制的投影。这也是为什么它对辐射归一化要求高——如果两期影像之间的系统辐射差异（比如太阳高度角不同、传感器增益不同）比真实变化还大，PCA第一主成分压出来的就是“整景辐射差异”，而不是区域性的地物变化。

2.2 工具包的整体处理链路：从原始影像到矢量图斑的五步流水线

我一般会把工具包按下面五步组织，每步的输入输出都分开，方便定位问题：

这个链条里最容易翻车的地方有两处：一处是辐射归一没做好导致PCA检测结果里全是系统性差异；另一处是矢量化后碎斑太多，面积过滤参数一紧又大面积漏检。后面第三章和第四章会对这两处重点展开。工具包如果只是把PCA封装成一个函数，那离“能用”还差得远，真正的工程量都在边界处理上。

2.3 数据对齐与辐射归一：决定检测成败的前置条件

先说配准。如果两期影像存在超过半个像元的偏移，差值图里会在地物边缘产生大量“伪变化”，PCA会把它们当成真实响应。我处理数据前必做一步：用两期影像的参考控制点或者同名地物点做一次检查，最省事的方式是直接套用已有正射影像产品，或用最近邻配准跑一遍，看叠加后地物边缘的重合情况。如果偏移超过一个像元，先做配准再进入工具包，这一步省不得。

再说辐射归一，常见做法有两种：线性拉伸和直方图匹配。线性拉伸的思路是假设两期影像在无变化区域（如裸地、稳定水体）的辐射存在一个线性关系，通过选取稳定的暗目标和亮目标，计算增益和偏移做校正。直方图匹配更粗暴：把后一期的累计直方图匹配到前一期的分布上，优点是自动，缺点是会把真实的大面积变化当成辐射差异一起抹掉——如果一个区域原本是大片农田，两期之间变成了建设用地，直方图匹配会把这个变化拉平一部分。我的习惯是先做线性拉伸，再检查归一化后的差值图像，看直方图主体是否接近零中心，如果还有明显的整体偏移，再做直方图匹配兜底。

在工具包里，归一化函数会放在PCA之前，独立的预处理脚本中。千万不要把它和PCA放进同一个函数——因为生产上经常出现只对其中一期做了重投影，另一期忘了做的情况，分开写更容易排查。

3. 用Python实现PCA核心：协方差矩阵、特征分解与变化强度图构建

3.1 从零实现PCA：为什么不用scikit-learn，以及最小可跑代码

很多初学者第一反应是调用sklearn.decomposition.PCA。但遥感影像有个现实问题：一个3000×3000像素、8波段的影像，展开后是900万×8的矩阵。sklearn的PCA默认对样本维度做中心化后计算协方差，在这个规模下矩阵本身不大（8×8），但fit过程里的中间数组会吃掉不少内存。更关键的坑在于，sklearn的PCA对输入数据的dtype和内存布局敏感，float64直接翻倍内存。影像计算本来就应该用float32，所以我在工具包里直接用NumPy实现，协方差矩阵只有B×B，几百个波段也就几十万个数，完全没有必要引入额外依赖。

import numpy as np from osgeo import gdal def read_tif(path): """读取多波段影像，返回(h, w, b)的float32数组和空间参考信息""" ds = gdal.Open(path) h = ds.RasterYSize w = ds.RasterXSize b = ds.RasterCount data = np.zeros((h, w, b), dtype=np.float32) for i in range(1, b + 1): data[:, :, i - 1] = ds.GetRasterBand(i).ReadAsArray().astype(np.float32) geo = ds.GetGeoTransform() proj = ds.GetProjection() ds = None return data, geo, proj def pca_diff_change(t1, t2, n_comp=1, standardize=True): """对两期影像差值堆叠做PCA，返回变化强度图与特征值""" assert t1.shape == t2.shape, "两期影像尺寸/波段数不一致" diff = t1 - t2 # shape: (h, w, b) h, w, b = diff.shape flat = diff.reshape(-1, b).astype(np.float32) # 剔除全零/无效像元，避免统计被边框噪声带偏 valid = np.any(flat != 0, axis=1) flat_valid = flat[valid] if standardize: mean = flat_valid.mean(axis=0) std = flat_valid.std(axis=0) std[std < 1e-6] = 1.0 # 防止常数为0的波段导致除零 norm = (flat_valid - mean) / std else: norm = flat_valid cov = np.cov(norm, rowvar=False) # 协方差矩阵，大小 b x b eigval, eigvec = np.linalg.eigh(cov) # 对称矩阵用eigh，数值更稳定 idx = np.argsort(eigval)[::-1] # 降序排列 eigval, eigvec = eigval[idx], eigvec[:, idx] score = norm @ eigvec[:, :n_comp] # 投影到前n_comp个主成分 change = np.zeros((h, w), dtype=np.float32) change[valid] = score[:, 0] return change, eigval, eigvec

代码逻辑说明：核心步骤是先做差，然后把二维影像展平成样本矩阵。standardize=True时按波段做z-score标准化，消除不同波段的辐射量纲差异——比如近红外波段的数值范围可能远大于蓝波段，不标准化的话，PCA会优先保留数值大的波段，这不符合变化检测的期望。np.linalg.eigh专门用于对称矩阵求特征值，比np.linalg.eig更快也更稳。change[valid] = score[:, 0]这一步把无效像元保留了原来的0值，避免在后续阈值分割时被误判。

参数方面，standardize推荐保持True；只有当两期影像已经做过严格的辐射归一、各波段数值范围本来一致时才可设为False。n_comp默认1，即只取第一主成分，后面第6章会讲什么时候需要取多个主成分做合成。

3.2 构建变化强度图：第一主成分的解读与可视化

拿到change数组后，要写回GeoTIFF才能进后续流程。写盘这个操作有个容易被忽视的细节：GDAL写入float32 GeoTIFF时，需要显式设置波段的数据类型，否则默认可能是Byte型，精度丢失后往往把微弱变化直接丢失。

def write_tif(path, data, geo, proj, dtype=gdal.GDT_Float32): rows, cols = data.shape driver = gdal.GetDriverByName('GTiff') ds = driver.Create(path, cols, rows, 1, dtype) ds.SetGeoTransform(geo) ds.SetProjection(proj) ds.GetRasterBand(1).WriteArray(data) ds.FlushCache() ds = None

主成分得分的符号本身没有物理意义——正分可以是变化也可以是负变化的方向选择问题。我在工具包里约定：change数组取绝对值后作为变化强度图，这样无论特征向量符号怎么翻转，下游步骤都稳定。这是踩过坑之后养成的习惯，细节在第5章避坑里展开。

可视化的参数建议：用2%～98%的百分位拉伸来做灰度显示，千万别用线性拉伸到全局最大最小值，否则个别极端像元会把整个图的对比度吃掉。变化强度图在这个阶段已经可以直接叠到原始影像上做初步目视判断了。

3.3 阈值分割：Otsu、百分比阈值与双峰校验

变化强度图是一张连续值的图，要变成“哪里变化了”，必须做阈值分割。工具包里我实现了三种方式：

import cv2 def threshold_change(change_img, method='otsu', percent=5.0): img = np.abs(change_img) # 去掉无效值后归一化到0-255 valid = img[img > 0] if len(valid) == 0: return np.zeros_like(img, dtype=np.uint8) if method == 'otsu': norm = cv2.normalize(img, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8) th, mask = cv2.threshold(norm, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU) return th, mask if method == 'percent': th = np.percentile(valid, 100 - percent) mask = (img >= th).astype(np.uint8) * 255 return th, mask

Otsu的原理是最大化类间方差，它在直方图呈双峰形态时效果最好。变化检测的特殊性在于：如果整个区域几乎没有变化，直方图只有单峰，Otsu会把靠近主峰边缘的噪声切出来，产生大量假图斑。所以我要求每次跑完阈值分割后，都按下面的逻辑做双峰校验：直方图峰值数与波谷深度的比值如果低于经验值，就自动切到百分比阈值法。

百分比阈值的含义是“假设全图有百分之几的像元发生了变化”，percent=5.0表示取第95百分位以上的像元为变化区。这个参数需要按业务场景调整，比如林地变化检测通常变化区域小于5%，我们设3左右；城区扩展检测变化面积往往更大，可以放宽到8～10。

3.4 内存不够怎么办：分块协方差与降采样策略

整景影像太大时，直接ReadAsArray()读完可能就快扛不住了。常见的绕法有两种。第一种是降采样：用GDAL的Overview或者ReadAsArray(buf_xsize, buf_ysize)读到目标分辨率，先把PCA参数（均值、标准差、特征向量）算出来，再用这些参数对原始分辨率的差值影像做投影。这个做法的前提是变化信号在降采样后依然存在，对大面积地物变化适用，对小碎斑变化不友好。

第二种是分块协方差：把影像切成若干块，每块算协方差矩阵，加权平均得到全图协方差。块与块之间的统计量可以并行计算，再汇总。这个方案更严谨，但实现起来要注意每块的有效像元数差异很大，加权时用有效像元数做权重。我在工具包里默认提供降采样路线，稳定够用，分块协方差作为可选函数留给熟练用户。

4. 栅格转矢量与图斑过滤：把像素结果变成可直接入库的成果

4.1 用GDAL Polygonize完成栅格矢量化

变化掩膜是一张0/1栅格，要变成矢量图斑，标准的做法是GDAL的Polygonize函数。它会把栅格中相邻的同一值像元合并成单个多边形。有一个绕不过去的坑：Polygonize输出的是包含外环的Polygon，但面积计算、拓扑修复（比如自相交）不做，生成的Shapefile里多边形与多边形之间经常存在微小的缝隙或重叠，这些要靠后面的过滤和修复步骤去收尾。

from osgeo import ogr, gdal def raster_to_vector(mask_path, out_shp): mask_ds = gdal.Open(mask_path) src_band = mask_ds.GetRasterBand(1) # 仅对值非0的像元（即变化区）做多边形提取 mask_band = src_band.GetMaskBand() drv = ogr.GetDriverByName('ESRI Shapefile') if drv.Exists(out_shp): drv.DeleteDataSource(out_shp) dst_ds = drv.CreateDataSource(out_shp) dst_lyr = dst_ds.CreateLayer('change', geom_type=ogr.wkbPolygon) # 添加属性字段 dst_lyr.CreateField(ogr.FieldDefn('DN', ogr.OFTInteger)) dst_lyr.CreateField(ogr.FieldDefn('AREA', ogr.OFTReal)) gdal.Polygonize(src_band, mask_band, dst_lyr, -1, [], callback=None) dst_ds = None mask_ds = None

Polygonize之后，多边形属性中DN值为1的就是变化区，DN值为0的是背景。我在工具包里默认把DN=0的多边形删掉，只保留变化图斑，否则后续过滤全是重复计算。注意gdal.Polygonize的第四个参数-1表示不进行八邻域连接，只做四邻域。遥感变化图斑绝大多数是四邻域连接就行，八邻域会把对角线接触的碎斑连成一片，容易吞掉细长的真实变化。反正后面还有图斑过滤，宁可在这一步保守，也不要在一步合并。

4.2 图斑过滤三件套：面积阈值、形状因子与边界优化

矢量化出来的初始图层几乎一定带着大量“看起来像变化、实际上只是噪声”的碎斑。工具包里我把过滤逻辑抽成了三个独立参数：最小面积、最小形状因子、最大边界像素比例。

面积阈值是硬标准——比如最小面积设为1000平方米，小于这个面积的图斑直接删除。这里有个容易搞混的地方：如果矢量图层的坐标系是地理坐标系（经纬度），直接算出来的面积单位是度，必须先把图层投影到以米为单位的投影坐标系再计算面积。我会在工具包里加一个函数，用ogr.Layer的GetSpatialRef()判断坐标系类型，非投影坐标系自动转UTM再面积过滤。

形状因子我用的是“紧凑度”：4π × 面积 / 周长²。圆的紧凑度是1，条带状图斑的紧凑度趋近于0。真实的地表变化（建筑扩展、植被退化）形状通常比较饱满，而噪声和边缘残差多为细长条。经验值：小于0.1的图斑直接删，不要心疼。

边界像素比例指图斑里落在原始掩膜边界上的像元占比——这个指标主要用于筛掉那些“像素级异常”的小图斑。因为栅格矢量化之后，一些图斑其实是椒盐噪声连成的小串，虽然面积达标，但内部孔洞极多、边界曲折。这类图斑用面积+紧凑度已经能过滤掉大半，边界比例参数作为补充，一般设置为0.3～0.5。

import geopandas as gpd import numpy as np def filter_polygons(gdf, min_area=1000, min_compactness=0.1, min_border_ratio=0.0): # 面积计算需要投影坐标系 if gdf.crs and not gdf.crs.is_projected: gdf = gdf.to_crs(gdf.estimate_utm_crs()) gdf['area_m2'] = gdf.geometry.area perimeter = gdf.geometry.length gdf['compactness'] = 4 * np.pi * gdf['area_m2'] / (perimeter ** 2 + 1e-8) # 边界比率：周长中多边形外环与背景接触的比例，这里用面积等效估算 cond_area = gdf['area_m2'] >= min_area cond_comp = gdf['compactness'] >= min_compactness filtered = gdf[cond_area & cond_comp] return filtered.reset_index(drop=True)

这段代码强依赖GeoPandas，需要读者本地装好geopandas、shapely、pyogrio这三件套。如果不想引入GeoPandas，也可以直接遍历ogr图层逐多边形做面积和周长计算，逻辑一样，只是写起来啰嗦。参数上，min_area一定要结合影像分辨率来定——2米分辨率影像上1000平方米相当于250个像元，0.5米分辨率影像上则是4000个像元，这两者的滤除效果截然不同，建议先用统计函数打印初始图斑面积分布，再定阈值。

4.3 输出与整合：GeoPackage与Shapefile的生产习惯

把过滤后的图斑写盘时，我一般直接用GeoPandas的to_file，格式选GeoPackage而不是Shapefile。原因一是GeoPackage是单文件、不用管一堆附属文件，二是属性能直接带类型信息，Shapefile的字段名长度限制和中文乱码问题则一言难尽。

def export_result(gdf, out_path): if gdf.crs is None: gdf = gdf.set_crs('EPSG:4326') # 按面积降序，方便在外业核查时优先看大面积变化 gdf = gdf.sort_values('area_m2', ascending=False) gdf.to_file(out_path, driver='GPKG')

输出前还有一件重要的事：把面积字段的精度控制一下，round到整数平方米就行，不然一个几百兆的GeoPackage里全是小数点后6位的面积属性，纯粹浪费空间。外业核查时用Google Earth或某离线地图打开GeoPackage，直接按面积排序浏览即可，这就是成果交付的常用形态。

5. PCA变化检测的五个常见坑与排查手册

5.1 主成分符号漂移：前后两次运行结果颜色反了

现象：同一组数据昨天跑出的变化强度图，高亮区域和今天的完全互补，或者同一块变化区域在两次运行中一个显示白色一个显示黑色，拿到阈值分割阶段后甚至检测区域完全不变。

原因：PCA的特征向量不是唯一的——特征向量乘以-1依旧是特征向量。当协方差矩阵的特征值非常接近（比如两个主成分方差差不多），特征向量的方向由数值计算的微小扰动决定，符号随机翻转。

解决：工具包里始终对变化强度图取绝对值再进下游，这是第一道保险。如果业务上必要区分变化方向（例如NDVI差值符号有明确意义），则在工具包外部做方向校验，比如拿一期NDVI差值与PCA得分做相关性分析，相关性为负就把特征向量乘-1。我自己的习惯是：输出图斑时只报告“哪里变了”，不区分增加或减少，方向性分析交给业务层，PCA层不背这个锅。

5.2 云和阴影被当成变化

现象：两期影像中一期有云，另一期无云，PCA的差值图中云覆盖区域出现成片的“变化”，面积大、形状完美，几乎无法过滤掉。

原因：云区的光谱特性与地面差异极大，在差值波段上的响应远强于大多数真实变化，PCA第一主成分自然优先表达这块。更阴险的是，云阴影也有类似效果，而且色调更接近真实地物的暗色。

解决：唯一可靠的办法是提前做掩膜。用两期影像各自的云掩膜做并集，把云和阴影覆盖的像元在差值堆叠前就置为无效值。工具包里我写了一个mask_cloud预处理函数，支持传入外部云掩膜文件；没有外部掩膜时，可以用亮度阈值粗筛——高云区在可见光波段的亮度明显偏高于周边地物。千万不能指望PCA自己“忽略”云区，它没这个能力。掩膜后的空洞，要在最终图斑过滤阶段用边界平滑或多边形简化把它闭合。

5.3 整景影像内存溢出

现象：一个10000×10000像素、8波段的影像，运行pca_diff_change时报MemoryError。

原因：diff.reshape(-1, b)瞬间生成了一亿行的二维数组，即使使用float32也需要约3.2GB内存，加上中间临时数组，内存直接爆掉。

解决：两个方向。方向一是降采样到1/4分辨率或更低，先把特征向量算出来，再用这个向量对全分辨率差值影像做逐块投影，投影时每块只负责一个2000×2000的窗口，逐块写盘，峰值内存能降到原来的1/25。方向二是用分块协方差，避免一次性载入全图。我在第三章给过降采样方案，注意降采样的方式要用rasterio或gdal的ReadAsArray(buf_xsize=buf_xsize)而不是先用原分辨率读完再numpy切片，原因还是内存。如果影像来自多个传感器，降采样前还需要统一重采样到相同的像元大小和网格原点，这一步在工具包外部做。

5.4 碎斑多到没法看：先形态学后矢量化才是正解

现象：Polygonize出来的矢量图层有几万个图斑，过滤参数调来调去都不理想——面积阈值调大，真实变化被一起削掉；调小，碎斑依然惨不忍睹。

原因：错失了在栅格阶段做形态学后处理的机会。阈值分割后、矢量化前，应该先对掩膜做开运算（先腐蚀后膨胀）和闭运算（先膨胀后腐蚀）。开运算的目标是去掉变化掩膜里的孤立点和小突起，闭运算用来填补变化区内部的孔洞——这些孔洞矢量化之后会变成内环，增加矢量数据的复杂度。

解决：OpenCV的两行代码，cv2.morphologyEx(mask, cv2.MORPH_OPEN, kernel)和cv2.morphologyEx(open_mask, cv2.MORPH_CLOSE, kernel)。kernel大小按影像分辨率定，2米分辨率影像用3×3或5×5足矣；0.5米分辨率影像可以用7×7，去掉的碎斑更多，但也更容易抹掉小尺寸真实变化。我的参数习惯是先开运算后闭运算，open的kernel比close小一档。形态学做完再矢量化，初始图斑数一般能降一个数量级。

5.5 直方图没有双峰：Otsu翻车时的替代方案

现象：变化区域占整景面积的不到1%，变化强度图的直方图只有一个主峰，Otsu求出的阈值落在主峰大腿上，检测结果全是主峰边缘的半噪声像元。

原因：Otsu天然假设数据是两类混合且比例接近。当两类比例悬殊时，类间方差最大化会偏向把占比大的类内部切开，而不是把两类分开。

解决：弃用Otsu，改用百分比阈值。先对全部变化强度图做百分位统计，取95分位、97分位甚至更高作为阈值，完全由业务预设变化比例决定。另一种替代是K-Means二分类，但效果和Otsu类似，同样对类别不均衡敏感。工具包里我把百分比阈值作为Otsu失败时的首选兜底方案。如果你不确定当前数据适不适合Otsu，就先跑一个快速诊断：打印直方图的前三个峰的位置和峰值，如果两个峰的高度差超过5倍，直接别用Otsu。

6. 结果可信度验证与进阶技巧：用混淆矩阵和主成分能量说话

工具包交付前，必须要回答一个问题：检测结果到底准不准。我一般会从图上随机抽100个图斑点（面积分档抽，大图斑多抽），配合原始影像人工判读每个点到底是真变化还是伪变化，然后用混淆矩阵和Kappa系数量化精度。

from sklearn.metrics import confusion_matrix, cohen_kappa_score def validate_by_samples(pred_ids, ref_ids): """ pred_ids: 图斑ID序列, 1表示检测为变化, 0表示非变化 ref_ids: 人工判读结果, 1表示确实变化, 0表示未变化 """ cm = confusion_matrix(ref_ids, pred_ids, labels=[0, 1]) kappa = cohen_kappa_score(ref_ids, pred_ids) return cm, kappa # 输出格式：cm[0,0]为真未变化, cm[0,1]为误检, cm[1,0]为漏检, cm[1,1]为正确检出

我习惯把Kappa在0.7以上定义为“结果可用”，0.4～0.7定义为“仅能做大范围排查”，低于0.4说明整个参数链有问题，回去检查辐射归一波段顺序或阈值选择。还有一个容易被忽略的参数调节思路：当漏检和误检不可兼得时，宁可在阈值分割阶段保留更多候选图斑，把误检压力留给后期图斑过滤。因为图斑过滤是有明确几何标准的，可以解释给外业人员听；而阈值切得太紧导致的漏检，可能在成果里彻底消失，无法挽救。

进阶技巧一个：不要只盯着第一主成分。工具包默认输出PC1，但可以先看特征值占比——eigval[0] / eigval.sum()，假如PC1贡献率只有40%而PC2、PC3各占25%左右，说明变化模式非常分散，光用PC1会漏掉一部分只在PC2上表现突出的变化。此时可以把PC1、PC2、PC3的能量合成一张综合变化强度图，score = sqrt(w1*pc1² + w2*pc2² + w3*pc3²)，权重用各自特征值占比。这样做的代价是噪声响应也会更显著，阈值分割时百分比参数要往高设一点。用这个技巧处理过一片混合了砍伐和建设变化的区域，漏检率降了大约三分之一。

工具包做到这个程度，基本字段、流程、坑都齐了。我自己的习惯是：每个项目跑到这一步，都会把这次的阈值参数、特征值占比、Kappa结果存成一个文本放到成果目录里——下次遇到类似的数据，先拿这些参数做初值，而不是从头调一遍。这样工具包越用越顺手，踩过的坑才真正变成了经验。希望这篇能帮你在PCA变化检测上少走几步弯路。

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