简介基于Python的哨兵二号Sentinel-2卫星数据像元三分法模型面向遥感、地理信息相关专业的学生与研究人员。包内实现最大噪声比变换MNF与像元纯度指数PPI两种光谱处理算法可对多波段影像进行降维压缩与纯像元筛选为后续像元三分法分解提供可靠输入。资源共33个文件以8个Python脚本为主体覆盖影像波段提取、归一化植被指数NDVI与差分特征指数DFI计算、像元纯度指数PPI计算、多波段融合及TIF格式输出另配22张PNG过程与结果图便于逐模块对照理解。压缩包约4.01MB结构紧凑适合直接运行与二次开发。目前已有440人浏览学习借助这套代码可快速搭建从原始影像到三分法分类的完整流程弄清MNF连续主成分降维和PPI阈值判定的实现细节为课程设计、算法复现或学位论文实验提供实用参考。1. 像元三分法模型把 Sentinel-2 影像拆成植被、非植被和水体占比做遥感课程设计或地表覆盖分析的同学大多会遇到同一个尴尬NDVI 只能告诉你“绿不绿”却说不清一个像元里到底有多少是植被、多少是裸土、多少是水体。城市里一个 10m 像元往往同时混着树冠、水泥地和阴影单一指数算出来的值既不像 A 也不像 B用来分类等于盲猜。这套基于 Python 的 Sentinel-2 像元三分法模型解决的正是“混合像元”问题——它通过 MNF 降维、PPI 提取纯像元再结合 NDVI 与 DFI 特征把每个像元分解成植被、非植被、水体三个端元的比例输出结果为后续分类或变化检测直接可用的丰度图。适合正在做课程设计、毕业设计或想从“会算指数”进阶到“会做混合像元分解”的从业者。2. 数据准备从 Sentinel-2 L2A 产品到统一分辨率波段栈2.1 为什么选 L2A 而不是 L1C拿到 Sentinel-2 数据的第一件事不是写代码是先确认产品级别。L1C 是大气顶部反射率没有做大气校正L2A 是地表反射率已经通过 Sen2Cor 处理过。像元三分法做的是地表组分分解端元光谱必须来自地表反射率否则大气散射会把植被和非植被的边界彻底搅浑。项目里 get_band.py 默认处理的就是 L2A 的 SAFE 目录结构这一点从文件命名和图层组织方式上能直接看出来。L2A 目录下 GRANULE 里每个 10m、20m、60m 波段都以独立 JP2 文件存放文件名里有 R60m、R20m、R10m 的标识。读数据时不能只用 gdal.Open 硬编码路径要用通配符匹配波段编号因为不同批次数据的分辨率目录名是一致的但绝对路径会因为解压位置不同而变化。import os import glob from osgeo import gdal import numpy as np def get_band(safe_path, band_id, target_res10): 从 L2A SAFE 目录中提取指定波段并重采样到目标分辨率 :param safe_path: SAFE 目录路径 :param band_id: 波段编号如 B4 :param target_res: 目标分辨率默认 10m :return: (波段数组, 仿射变换参数, 投影) # 匹配 GRANULE 下的 JP2 文件文件名形如 ..._B04_10m.jp2 pattern os.path.join(safe_path, GRANULE, *, IMG_DATA, *, f*{band_id}_{target_res}m.jp2) matches glob.glob(pattern) if not matches: raise FileNotFoundError(f未找到波段 {band_id}检查路径或分辨率参数) ds gdal.Open(matches[0]) arr ds.ReadAsArray().astype(np.float32) geo ds.GetGeoTransform() proj ds.GetProjection() ds None # L2A 地表反射率存储为整数除以 10000 得到 0~1 的反射率 arr arr / 10000.0 # 无效值掩膜L2A 中 0 通常表示无数据 arr[arr 0] np.nan return arr, geo, proj # 示例读 B4红和 B8近红外 red, geo, proj get_band(S2A_L2A_20240101T002001_XXXX, B04) nir, _, _ get_band(S2A_L2A_20240101T002001_XXXX, B08)这段代码做了三件事用 glob 匹配波段文件、除以 10000 转反射率、把无效值置为 NaN。除以 10000 这一步是关键直接拿原始整数算 NDVI 会得到一堆看似合理实则错误的数值后面我会在避坑章节专门展开。如果你需要 20m 波段参与计算把 target_res 改成 20 即可但后续必须统一到同一个网格上否则数组形状对不上。2.2 波段选型与统一网格的重采样策略三分法模型在 Sentinel-2 上常用的波段组合是 B2蓝、B3绿、B4红、B8近红外这四个 10m 波段加上 B11、B12 两个 20m 的短波红外波段。B8 和 B4 算 NDVIB11、B12 算 DFIRGB 三波段用于后续的可视化检查和 GF-2 融合对比。20m 波段不能直接参与逐像元计算需要先重采样到 10m 网格。常见做法是用 gdal.Warp 做最邻近重采样。最邻近不会改变原始波段的数值分布适合后续做比值指数双线性内插虽然平滑但会把 20m 的短波红外信息“涂抹”到相邻像元上对端元纯度分析有干扰。重采样之后必须检查数组形状是否与 B4、B8 完全一致一个简单的 assert 就能省掉后面所有因 shape 不匹配导致的报错。3. 特征构建从 NDVI 和 DFI 到二维特征空间3.1 NDVI 计算中的数组运算陷阱NDVI 的公式是 (B8 - B4) / (B8 B4)看起来简单但 Sentinel-2 的 B8 和 B4 都是 10m 分辨率直接做减法和除法即可。需要留意的是数据类型如果读出来是 int16除法和减法都会向下取整NDVI 值会变成一堆 0 和 1完全没有区分度。get_ndvi.py 里的实现先把数组转成 float32再对分母为 0 的位置做掩膜处理。3.2 DFI 指数的取值逻辑与物理含义DFI 在项目里用的是短波红外波段的归一化差值B111610nm和 B122190nm对水分和矿物成分敏感。裸土和干燥植被在 B11、B12 上的表现差异明显而水体在两个波段上几乎全吸收。用 DFI 和 NDVI 组成二维特征空间后植被、非植被、水体三个端元会分别落在特征空间的不同角落这是三分法能成立的基础。实际计算时要注意 B11、B12 是 20m 分辨率必须先重采样到 10m 再参与运算。DFI 的定义在公开文献里有多种写法有的取 (B11-B12)/(B11B12)有的取比值后做对数变换。项目里 get_dfi.py 用的是前者代码层面保持与 NDVI 一致的归一化差值形式便于后续在同一个特征空间里做投影分析。import numpy as np def get_dfi(b11, b12): 计算 DFI归一化差值短波红外指数 :param b11: 重采样后的 B11 波段数组 (float32) :param b12: 重采样后的 B12 波段数组 (float32) :return: DFI 指数数组值域约 [-1, 1] # 先做掩膜任一波段为 NaN 的位置直接置 NaN mask np.isnan(b11) | np.isnan(b12) dfi np.full(b11.shape, np.nan, dtypenp.float32) # 对有效像元计算归一化差值 valid ~mask numerator b11[valid] - b12[valid] denominator b11[valid] b12[valid] # 分母接近 0 时几乎不可能出现在真实地表结果置 0 防止除零 denominator np.where(np.abs(denominator) 1e-6, 1e-6, denominator) dfi[valid] numerator / denominator return dfiDFI 的计算逻辑和 NDVI 完全对称唯一多出来的是分母保护。真实遥感数据里 B11 和 B12 同时为 0 的情况极少但做批量处理时难免遇到坏像元加一个 epsilon 保护能避免运行时弹出 RuntimeWarning 导致后续端元选取时出现 NaN 传播。3.3 calculate_f.py 做的事特征融合与投影矩阵get_ndvi.py 和 get_dfi.py 分别输出单波段指数图calculate_f.py 的作用是把两个指数在像元维度上叠成一个双通道特征张量并计算每个像元在特征空间中的坐标。纯植被像元在 NDVI 上接近 0.8 以上、DFI 上接近 0水体像元在 NDVI 上为负值、DFI 上也是明显的负值裸土则落在 NDVI 低、DFI 略负或接近 0 的区域。这三个端元的坐标就是后面 PPI 要寻找的目标。这一章的程序输出是 F 矩阵形状为 (2, 行数, 列数)。注意这里通道在前是为了方便后续用 numpy 的 einsum 做批量投影比 for 循环遍历像元快两个数量级。4. 核心算法MNF 降维与 PPI 纯像元提取4.1 为什么三分法前要先做 MNF很多人第一次看到这个项目的文件列表会疑惑为什么要用 MNF 和 PPI 两个看起来像是高光谱处理的算法来处理 Sentinel-2 多光谱数据原因是波段之间相关性太强。B2、B3、B4 在植被覆盖区域的变化几乎同步直接对原始波段做端元提取端元光谱会高度共线后续的最小二乘解不稳定。MNF最大噪声比变换通过两次主成分分析把波段重组成按信噪比排序的新分量前几个分量集中了地表真实变化信息后面的分量以噪声为主。MNF 和 PCA 的核心区别是PCA 按方差排序而方差不等于信号——高相关波段里的噪声方差也可能被放大MNF 先估计噪声协方差矩阵把数据白化后再做 PCA这样排序依据的是信噪比而不是总方差。import numpy as np def mnf_transform(data, n_components10): 最大噪声比变换先差分估计噪声协方差再做白化 PCA :param data: 输入数据形状 (波段数, 像元数)已去除无效值 :param n_components: 保留分量数 :return: 变换矩阵和降维后的数据 bands, pixels data.shape # 步骤1用差分法估计噪声协方差 # 相邻像元之差主要包含噪声因为真实地表信号在空间上是平滑的 diff np.diff(data, axis1) noise_cov np.cov(diff) / 2.0 # 对噪声协方差做特征分解得到白化矩阵 eigval, eigvec np.linalg.eigh(noise_cov) eps 1e-8 whitening eigvec np.diag(1.0 / np.sqrt(eigval eps)) eigvec.T # 步骤2对白化后的数据做 PCA white_data whitening data cov np.cov(white_data) _, proj_vec np.linalg.eigh(cov) # 特征值升序排列取后 n_components 个向量 proj_vec proj_vec[:, ::-1][:, :n_components] # 组合变换矩阵 transform proj_vec.T whitening transformed transform data return transformed, transform # 实际调用时把 F 矩阵 reshape 成 (波段数, 像元数) 的二维数组 # f_stack 形状为 (2, rows, cols)去掉 NaN 像元后参与计算MNF 实现里的关键参数是噪声协方差的估计方式。差分法是最基础的方案假设相邻像元的信号几乎一致差分结果以噪声为主。实际使用时如果影像噪声较小差分法会低估噪声水平导致白化矩阵不稳定。我一般在做城市区域时会改用 3x3 窗口的局部标准差来估计噪声代价是计算时间长一点但结果更稳。4.2 PPI 纯像元提取的随机投影机制PPI像元纯度指数的目标是在 MNF 降维后的特征空间里找“最极端”的像元。算法原理是生成大组随机单位向量把每个像元投影到所有随机向量上记录该像元出现在投影极值位置的次数。出现在极值位置的次数越多说明它越靠近特征空间的凸包顶点纯度越高。import numpy as np from numpy.linalg import norm def calculate_ppi(data_mnf, n_vectors10000): 像元纯度指数计算 :param data_mnf: MNF 降维后的数据形状 (n_components, n_pixels) :param n_vectors: 随机投影向量数量常用 10000~50000 :return: ppi 分数数组长度等于像元数以及端元候选索引 n_comp, n_pix data_mnf.shape # 生成随机投影向量并归一化 rng np.random.default_rng(42) vectors rng.standard_normal((n_vectors, n_comp)) vectors vectors / norm(vectors, axis1, keepdimsTrue) # 投影计算所有像元在每条随机向量上的投影值 # 用矩阵乘法一次完成避免 for 循环 proj vectors data_mnf # 形状 (n_vectors, n_pix) # 统计每个像元出现在极值位置的次数 ppi_scores np.zeros(n_pix, dtypenp.int32) # 逐条向量找最大值和最小值位置 for i in range(n_vectors): row proj[i] ppi_scores[np.argmax(row)] 1 ppi_scores[np.argmin(row)] 1 # 按分数降序排列取前若干作为端元候选 candidate_idx np.argsort(ppi_scores)[::-1] return ppi_scores, candidate_idxPPI 的参数里最敏感的是随机向量数量。向量太少极值统计不稳定同一个端元可能被漏掉向量太多计算时间线性增长。经验值是 10000 起步如果端元在特征空间里分布分散可以加到 50000。另一个隐蔽参数是随机数种子项目里固定为 42这能保证结果可复现——做课程设计交报告时两次运行结果不一致会很尴尬。4.3 三分法端元选取与最小二乘解算PPI 选出候选像元后在特征空间里做 K-Means 聚类或直接取极值簇中心作为三个端元的特征坐标。得到端元坐标后每个像元的组分比例就是线性光谱混合模型的解像元反射率 植被比例 × 植被端元 非植被比例 × 非植被端元 水体比例 × 水体端元加上一个误差项。约束条件是三个比例之和为 1 且每个比例在 0~1 之间这是一个带约束的最小二乘问题。项目中这部分逻辑没有单独成文件而是放在 test.py 里做了完整的流程串联。求解时用 scipy.optimize.lsq_linear 或手动实现投影到单纯形上的算法都可以数据量不大时用 SLSQP 也能跑。5. 避坑排查像元三分法跑不出结果的五个真实原因5.1 NDVI 全图异常偏小或出现大量负值现象NDVI 最大值不超过 0.3水域像元甚至出现 -1 附近的极值植被像元数值和裸土难以区分。原因直接用了 JP2 原始整数 DN 值计算没有除以 10000。L2A 产品的存储值是地表反射率乘以 10000 后的整数DN 值范围 0~10000不缩放直接代入公式NDVI 依然会落在 [-1,1]但植被区域的 NDVI 会被压缩到 0.1 以下整个直方图畸形。解决在 get_band.py 里强制除以 10000.0并且检查输出数组的最大值是否小于等于 1。如果最大值还是几百几千说明缩放步骤没生效。5.2 20m 波段重采样后影像出现黑边或错位现象DFI 计算完成后影像边缘出现一条不规则的黑色区域或者地物在 B11 和 B8 之间有半个像元的位移。原因用 gdal.Warp 重采样时没有指定输出范围GDAL 默认按目标分辨率的整数网格对齐导致 20m 波段和 10m 波段的像元网格存在半个像元的偏移边缘区域出现空值。解决重采样时显式指定输出范围和输出网格用 B4 波段的 GetGeoTransform 作为基准。from osgeo import gdal def resample_to_match(src_path, ref_ds_path, out_path): 将 src 波段重采样到与参考波段完全相同的网格 ref gdal.Open(ref_ds_path) geo ref.GetGeoTransform() width ref.RasterXSize height ref.RasterYSize ds gdal.Warp( out_path, src_path, formatGTiff, xResgeo[1], yResabs(geo[5]), outputBounds(geo[0], geo[3] height * geo[5], geo[0] width * geo[1], geo[3]), resampleAlgnear, dstSRSref.GetProjection() ) ds None这里的 outputBounds 计算是四则运算但很容易写错符号——注意 GeoTransform 的第五个参数是负值北在上所以 y 方向的下边界是 geo[3] height × geo[5]而不是减号。5.3 MNF 变换后第一分量仍是噪声条纹现象MNF 第一分量图看起来不是地物信息而是横向或纵向的条纹像是传感器本身的条带噪声。原因差分法估计噪声协方差时如果影像存在系统性条带噪声相邻像元的差分结果会把条带当作“信号”纳入噪声统计导致白化矩阵失真。解决改用局部均值残差法估计噪声协方差即用 3x3 均值滤波后的影像减去原始影像残差作为噪声样本。代价是计算量大约增加两倍但对条带噪声的鲁棒性明显更好。5.4 PPI 提取的端元全是云和雪现象PPI 分数最高的前几十个像元集中在云的亮白区域或积雪覆盖区而不是目标地物。原因云和雪在短波红外波段的反射率远高于普通地表在特征空间里天然处于凸包顶点PPI 的极值统计对它们最敏感。没有做云掩膜就提端元等于把噪声当信号。解决在进入 PPI 前先用 B2蓝波段的阈值结合 NDVI 负值区域做云和雪的粗略掩膜把掩膜像元排除在 PPI 计算之外。精确做法是用 Scene Classification Layer但课程设计场景下阈值掩膜足够用。5.5 to_tif 输出的文件在 QGIS 中打不开现象to_tif.py 生成的结果 TIFF 在 GDAL 命令行里能读取但 QGIS 打开后显示全黑或报错“未找到地理参考信息”。原因输出时用了老的 GTiff 驱动默认参数没有保留投影信息或者把 NaN 值直接写入 float32 波段导致显示范围异常。解决输出时显式设置投影、仿射变换和 NoData 值并在写出前把 NaN 替换为设置的 NoData 值。QGIS 全黑通常只是显示范围问题设置“自动重新调整显示范围”即可。6. 进阶技巧GF-2 融合对比与输出 TIFF 的质量控制6.1 fuse.py 的融合策略与用途项目中 fuse.py 处理的是 GF-2 数据的融合场景。GF-2 全色影像分辨率 1m多光谱影像分辨率 4m融合后得到的 1m 多光谱影像可以用来验证 Sentinel-2 三分法结果在空间细节上的可靠性。常见做法是 Gram-Schmidt 融合把 Sentinel-2 的 10m NDVI 特征图当成低分辨率模拟影像与 GF-2 全色波段做锐化得到 1m 分辨率的特征图再和三分法输出的丰度图对比。融合的目的一是生成一张适合制图输出的高分辨率成果图二是通过两期数据的丰度变化检测验证三分法模型在城市绿地监测上的可用性。如果你手里的数据不是 GF-2而是 SPOT-6 或高分一号fuse.py 的接口同样适用只要把全色和多光谱文件路径传进去即可。6.2 to_tif.py 输出的硬性检查清单每次跑完 to_tif.py我会按固定顺序做三件事这三件事已经写进我自己的项目习惯里。第一是检查投影——用 gdal.Info 看输出文件是否带完整 WKT 投影字符串如果只有带编号的投影参数说明写入时丢了基准面信息第二是检查 NoData——用 gdal_translate -stats 看统计信息里 NoData 占比是否异常高于预期第三是检查值域——三分法的比例输出应该在 0~1 之间最小值出现负数说明端元坐标选取有误。# 1. 查看投影和空间参考完整性 gdalinfo output_abundance.tif | grep -A 3 Coordinate System # 2. 查看 NoData 值是否写入 gdalinfo output_abundance.tif | grep NoData # 3. 用 Python 检查值域分布 python -c from osgeo import gdal import numpy as np ds gdal.Open(output_abundance.tif) arr ds.ReadAsArray() print(min:, np.nanmin(arr), max:, np.nanmax(arr)) print(NaN 占比:, np.isnan(arr).mean()) 这三个命令能发现 90% 以上的输出问题。特别是最后一条我遇到过输出的植被丰度图最小值是 -0.012原因是端元选取时把两个非植被端元的光谱区分度拉得太大最小二乘解出现轻微负值此时用 clip 到 [0,1] 处理只是面子工程根因要从端元纯度和特征空间分布去找。从那以后我每次跑完一套 Sentinel-2 三分法流程都会强制多走一步把三个端元在 NDVI-DFI 特征空间里的位置画出来看一眼确认它们确实分布在三个角落而不是挤在一起。这个步骤不花一分钟却能提前拦住一半的端元选型问题。三分法模型的变量多从波段选择、缩放、重采样到 MNF 参数、PPI 随机种子任何一个环节松懈都会在结果图上以隐蔽的方式体现出来养成检查特征空间的习惯之后翻车概率会小很多。希望帮到你。本文还有配套的精品资源点击获取