简介：基于Python的哨兵二号（Sentinel-2）卫星影像像元三分法模型资源包，面向遥感科学与技术、地理信息科学等专业的课程设计和科研入门者，重点解决中等分辨率影像中混合像元分解的实现问题。资源围绕最大噪声比变换（MNF）和像元纯度指数（PPI）两个核心环节，提供八个Python脚本，完整覆盖波段提取、归一化植被指数与特征指数计算、特征融合、像元纯度分析以及三分模型构建等步骤，并配有二十二张过程结果图片和一份说明文档，便于随时对照中间结果理解算法细节。整个压缩包共三十三个文件，整体大小约四点零一兆字节，结构清晰、代码组织紧凑，适合作为课程设计参考或算法复现的蓝本。目前已有四百三十九人下载学习，对希望快速上手哨兵二号数据预处理与像元分解算法的读者具有直接参考价值。
1. 像元三分法的思路与 Sentinel-2 的自然契合点
一个 10m 分辨率的 Sentinel-2 像元，落在农田边缘时通常同时包含植被、裸土和阴影。你问 NDVI 多少，它给你一个 0.4 的中间值，但说不出这 0.4 到底是“植被稀疏的裸坡”还是“长势中等且部分被遮挡的作物”。像元三分法直接把这个混合像元拆成三个端元的比例丰度：植被占多少、土壤占多少、阴影占多少，对应输出三张连续数值图。对做植被覆盖度估算、退耕还林监测、撂荒地识别的人来说，这套方法的可解释性远超单一指数。用 Python 实现时，整个流程只需要 rasterio 读波段、NumPy 做矩阵运算、Matplotlib 出图，端元和求解过程都能在本地完整复现，下面从模型原理一直走到验证技巧。
2. 线性混合模型与 Sentinel-2 端元选取方案
2.1 像元三分法为什么用线性混合模型
像元三分法本质上是线性光谱混合模型把端元数限定为 3。它假设像元反射率是各端元反射率的加权和，权重就是丰度。公式写出来是：
R = E × F + ε
其中 R 是像元的反射率向量，E 是端元光谱矩阵，F 是丰度向量，ε 是残差。植被冠层、裸土、阴影在亚像元尺度上基本不透明，光子打到哪个端元就反射哪个端元的信息，二次反射的贡献在 10m 分辨率下大多低于传感器噪声，所以线性假设在这个场景成立，不需要引入非线性混合模型。
选择三个端元而不是两个或四个，背后是自由度问题。两个端元只能表达一条光谱线上的比例，无法处理阴影这个无处不在的暗端元；四个以上端元会导致方程组病态，丰度解在不同波段组合间剧烈跳变。三个端元的典型组合是植被-土壤-阴影，在湿地或水体密集区可以把阴影替换成阴影/水体端元。Sentinel-2 的可见光-近红外波段恰好能把这三个端元分开：植被在红波段强吸收、近红外高反射，土壤在红和近红外都相对平稳，阴影在所有波段都接近低值。
用 3 个波段解 3 个端元，方程组是方阵，能唯一求解但无法给残差。加入第 4 个波段后变成超定方程组，多出来的自由度可以计算每个像元的分解误差，这是验证模型适用性的关键工具。所以实践中的经典配置是 B3、B4、B8 三个 10m 波段做主解，B11 做超定验证。
2.2 Sentinel-2 波段选择与端元对
不是波段越多越好，关键是端元之间的光谱可分性。B3、B4、B8 三个波段覆盖了绿色植物反射峰两侧和红边起点，植被和裸土在这三个波段上的差异最明显。B11 虽然对土壤含水量敏感，但它是 20m 分辨率，参与计算前必须重采样到 10m，会引入相邻像元混合的额外误差。
| 波段 | 中心波长 | 原始分辨率 | 在三分类中的角色 |
|---|---|---|---|
| B3 绿 | 560nm | 10m | 植被中等反射、土壤中等反射，区分阴影与亮度地物 |
| B4 红 | 665nm | 10m | 植被强吸收，分离植被与非植被的关键波段 |
| B8 近红外 | 842nm | 10m | 植被高反射、土壤中等反射，阴影端元与植被端元差异最大 |
| B11 短波红外 | 1610nm | 20m | 对土壤湿度敏感，用于构建超定方程组和残差验证 |
选波段时要注意端元矩阵的条件数。B2 蓝光和 B4 红光在大气瑞利散射后的变化模式高度相关，放进同一个方程组容易让矩阵接近奇异，求解时丰度值会出现正负大幅摆动。B3、B4、B8 之间相关性相对低，端元矩阵的典型条件数在 10 到 50 之间，求解稳定。检查方法很简单，在 Python 里对端元矩阵调用np.linalg.cond(E)，条件数超过 100 就说明波段组合需要调整。
2.3 从影像本身提取端元光谱的 Python 实现
端元光谱有两个来源：光谱库和影像本身。光谱库里的纯植被光谱是实验室或机载传感器测的，波段设置和观测几何与 Sentinel-2 不一致，直接用会让丰度结果带系统性偏差。从影像本身提取端元更常见，因为大气残余、观测几何、物候状态都包含在影像里，端元和实际场景同时刻匹配。常见做法是先算 NDVI 找纯植被像元，按亮度找裸土像元，取全影像辐亮度最低的暗像元作阴影参考。
import numpy as np def compute_endmembers(nir, red, green): nir = nir.astype("float32") red = red.astype("float32") green = green.astype("float32") ndvi = (nir - red) / (nir + red + 1e-6) # 纯植被：NDVI 最高的 0.1% 像元 v_quant = np.nanpercentile(ndvi, 99.9) veg_mask = ndvi > v_quant veg = np.array([ np.nanmean(green[veg_mask]), np.nanmean(red[veg_mask]), np.nanmean(nir[veg_mask]) ]) # 裸土：NDVI 接近 0 且平均亮度最高的 0.1% 像元 brightness = (green + red + nir) / 3.0 soil_mask = (np.abs(ndvi) < 0.05) & (brightness > np.nanpercentile(brightness, 99.9)) soil = np.array([ np.nanmean(green[soil_mask]), np.nanmean(red[soil_mask]), np.nanmean(nir[soil_mask]) ]) # 阴影/水体：三波段总和最低的 0.1% 像元 total = green + red + nir dark_mask = total < np.nanpercentile(total, 0.1) shadow = np.array([ np.nanmean(green[dark_mask]), np.nanmean(red[dark_mask]), np.nanmean(nir[dark_mask]) ]) return np.vstack([veg, soil, shadow]) # shape (3, 3)，行顺序：植被、土壤、阴影代码里的百分位阈值 99.9 和 0.1 依赖影像像元总数，一景完整的 Sentinel-2 切块有千万级像元，取 0.1% 足够稳定；如果用的是几百米的小试验区域，建议放宽到 99.5 和 0.5，避免端元被个别异常像元主导。1e-6是防止 NDVI 分母为零的平滑项，np.nanmean保证云掩膜后的 NaN 不会污染端元均值。土壤掩膜里的abs(ndvi) < 0.05在植被盖度很高的影像里可能筛不出像元，备选做法是去掉 NDVI 限制、直接取亮度最高的 0.1% 但再排除掉 NDVI 高于 0.3 的像元。
3. 用 Python 搭建 Sentinel-2 三分法数据管线
3.1 预处理与数据标准化
Sentinel-2 数据有 L1C 和 L2A 两种常用级别。L1C 是大气表观反射率，使用前需要做大气校正；L2A 经 Sen2Cor 处理后已经是地表反射率，可以直接用。收到 L2A 时要注意无效值：无数据区填 0，云掩膜区在某些产品中填 0，这些 0 值如果不处理，会让端元提取里的暗像元全部落到无数据区上。下面的函数把 B3、B4、B8 读进一个三维数组，顺便把非正值替换成 NaN。
import rasterio as rio import numpy as np def read_s2_stack(band_paths): """读取三个 Sentinel-2 波段，返回 (rows, cols, 3) 的 float32 数组""" bands = [] for path in band_paths: with rio.open(path) as src: arr = src.read(1).astype("float32") arr[arr <= 0] = np.nan bands.append(arr) stack = np.stack(bands, axis=-1) return stack # 通道顺序与 band_paths 保持一致逻辑说明：band_paths按顺序传 B3、B4、B8 的文件路径，函数内逐个读取，最终堆叠成通道在最后一维的数组，正好对应端元矩阵的列顺序。arr <= 0替换为 NaN 是因为 L2A 产品中 0 通常代表无数据或无效观测，保留的话会让 percentile 统计失真。读取 20m 的 B11 时，需要先用rio.warp.reproject或 GDAL 重采样到 10m 网格，再做波段对齐，这一步放在读取之后、端元提取之前。
预处理里另一个常见问题是单波段数据的 CRS 和 transform 不完全一致。正式处理时建议先在 GIS 软件里把各波段统一到同一网格，或者用rasterio.merge做一次强制对齐，否则后面按像素运算时会出现半像素错位，导致 NDVI 计算出现条带状噪声。
3.2 云和云影掩膜：端元提取前必须做
一景影像里只要有一片薄云，云顶的高反射率就会被当成裸土端元，整个像元三分法输出都会偏移。L2A 产品附带的 SCL 波段按像元类别编码，可以直接当掩膜用。SCL 的类别定义是：0 无数据、1 缺陷像元、2 暗区、3 云影、4 植被、5 裸土、6 水体、7 低概率云、8 中概率云、9 高概率云、10 卷云、11 雪。
def mask_cloud_scl(scl_path, stack): """用 SCL 波段剔除云和云影，保留的类别用于端元提取和丰度求解""" with rio.open(scl_path) as src: scl = src.read(1) # 保留类别：暗区、云影、植被、裸土、水体 good_codes = [2, 3, 4, 5, 6] good_mask = np.isin(scl, good_codes) # 每个像元三个波段要么全保留，要么全置 NaN valid = np.broadcast_to(good_mask[..., None], stack.shape) masked = np.where(valid, stack, np.nan) return masked, good_mask这里的掩膜逻辑是整像元左右：云、云影、雪都不参与后续计算，而被保留的暗区、云影、水体对阴影端元同样有贡献。暗区在 SCL 里对应地形阴影，云影和阴影光谱接近，水体也呈暗色，把这三类都保留下来能增加暗端元的样本量。掩膜做完后，应该打印一下good_mask.mean()，如果有效像元占比低于 70%，说明影像质量有问题，后续丰度图的空洞会很大，要考虑换一景时相。
3.3 环境与包管理
这一步专门说运行环境，因为 rasterio 和 GDAL 的版本冲突是新手最常卡住的地方。创建独立 conda 环境是稳妥做法，不碰系统的 Python：
conda create -n s2env python=3.10 -y conda activate s2env conda install -c conda-forge rasterio numpy scipy matplotlib -ypip 方式也可以，但建议加国内镜像地址，rasterio 的安装包体积大，直接走默认源容易超时：
pip install rasterio numpy scipy matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple如果你在 vscode python 环境配置里发现import rasterio报 ModuleNotFoundError，先确认左下角解释器选中的是 s2env 而不是全局环境。一个容易忽略的细节是 conda 安装 rasterio 时会自动解析 GDAL 依赖，而 pip 安装的 rasterio 依赖系统里已有的 GDAL，版本不匹配就会出现典型的“请安装缺失的包以使用此工作流”类提示或 OSError。遇到这类问题，最直接的办法是把当前环境的 rasterio 和 GDAL 一起卸载，再用 conda-forge 统一安装。
4. 像元三分法主程序：求解、约束与堵住常见错误
4.1 向量化的无约束最小二乘与后处理
核心求解可以一次矩阵乘法完成。端元矩阵 E 的形状是 (3, 3)，对全影像所有像元的反射率矩阵 R（形状 N×3）求丰度矩阵 F（N×3），无约束最小二乘解是 F = R × pinv(E)。np.linalg.pinv是对 E 求伪逆，比直接求逆更稳，即使 E 的条件数稍高也能给出可用结果。物理上丰度必须满足两个约束：每个端元丰度在 0 到 1 之间、三个丰度之和等于 1。下面用后处理方式施加约束：先截断负值，再归一化到和为 1。
def unmix_vectorized(stack, endmembers): rows, cols, bands = stack.shape flat = stack.reshape(-1, bands).astype("float32") valid = np.isfinite(flat).all(axis=1) # 有效像元才参与运算，无效像元保持 NaN f_un = np.full((flat.shape[0], 3), np.nan, dtype="float32") f_un[valid] = flat[valid] @ np.linalg.pinv(endmembers).T # 后处理约束：非负 + 和为一 f_pos = np.maximum(f_un, 0.0) f_pos[~valid] = np.nan f_sum = np.nansum(f_pos, axis=1, keepdims=True) f_norm = f_pos / np.maximum(f_sum, 1e-6) # 残差 RMSE：用归一化后的丰度重建反射率，衡量分解质量 rec = np.full_like(flat, np.nan) rec[valid] = f_norm[valid] @ endmembers rmse = np.sqrt(np.nanmean((flat - rec) ** 2, axis=1)) f_3d = f_norm.reshape(rows, cols, 3) rmse_2d = rmse.reshape(rows, cols) return f_3d, rmse_2d参数和逻辑说明：endmembers的行顺序必须与stack的波段顺序对应，代码里是绿色、红色、近红外。pinv(endmembers).T的转置是因为求解公式里端元矩阵按列排列，用行向量形式的反射率做右乘时要转回去。np.maximum(f_sum, 1e-6)防止归一化时除零，这个平滑项只影响那些丰度和接近零的退化像元。RMSE 计算用的是重建反射率与实际反射率的逐波段均方根，值越接近 0 说明三端元线性组合越能解释该像元，RMSE 超过 0.05 的像元通常落在云边缘或水体波浪区这类模型不适用的位置。
4.2 严格约束版本与场景选择
后处理本质上是在无约束解上做投影，当像元光谱落在端元三角形外部时，截断加归一化的结果和真正约束最小二乘解有偏差。对精度要求高的研究，或者像元落在端元三角形外的比例超过 10% 时，用 scipy 的nnls更严谨。
from scipy.optimize import nnls def unmix_nnls(flat, endmembers, valid): """严格非负最小二乘 + 和为一归一化，适合小样本验证""" A = endmembers.T # (3, 3)，nnls 要求列是端元 out = np.zeros_like(flat) for i in range(flat.shape[0]): if not valid[i]: out[i] = np.nan continue x, _ = nnls(A, flat[i]) s = x.sum() out[i] = x / s if s > 1e-6 else x return out这个版本逐像元循环，速度慢，不做全量部署，只推荐用在验证集、小范围试验或融合时序分析的抽样点上。两种版本在大多数正常像元上给出的丰度值差异在 0.02 以内，所以全图处理用向量化版、重点区域验证用 nnls 版，是性价比最高的组合。
4.3 三个常见错误及其判断依据
| 错误现象 | 可能原因 | 检查方法 |
|---|---|---|
| 丰度图出现整片高值条纹 | 端元矩阵条件数过高，多波段相关性强 | 打印np.linalg.cond(endmembers)，超过 100 则减少波段数 |
| 阴影丰度在水体区域接近 1 | 阴影端元和水体端元混淆 | 对比 SCL 类别 6 的水体像元，若阴影丰度均大于 0.8 属正常 |
| RMSE 总体偏高且多出现在边缘 | 影像波段间未严格配准 | 检查各波段的 transform 和边界是否完全一致 |
NaN 传播是最隐蔽的问题。stack里一个波段是 NaN，np.isfinite(flat).all(axis=1)已经把该像元标记为无效，但f_norm在 reshape 回三维后可能会和原始掩膜错位，所以掩膜数组valid也要 reshape 回二维，用它统一控制后续所有分析。还有一点：不要对整幅影像直接算 percentile 来提取端元，要先按 3.2 节的云掩膜过滤，否则云边缘的混合像元会进入端元候选集。
5. 成果图输出与像元三分法的模型验证
5.1 丰度图组合输出
三分法的输出是三个通道，直接画成 RGB 假彩色图最容易看出空间格局。把植被丰度放绿色通道、土壤丰度放红色通道、阴影丰度放蓝色通道，混合色块的分布能直观显示地表覆盖的空间异质性。用 Matplotlib 输出时，注意裁剪显示范围，避免水体和云掩膜的 NaN 影响色彩拉伸。
import matplotlib.pyplot as plt from matplotlib.colors import TwoSlopeNorm fig, axes = plt.subplots(1, 3, figsize=(15, 5)) titles = ["vegetation", "soil", "shadow"] for ax, title, arr in zip(axes, titles, [f_3d[..., 0], f_3d[..., 1], f_3d[..., 2]]): im = ax.imshow(arr, cmap="YlGn", vmin=0, vmax=1) ax.set_title(title) ax.axis("off") plt.colorbar(im, ax=ax, fraction=0.046) plt.tight_layout() plt.savefig("abundance.png", dpi=150) plt.show()vmin=0, vmax=1强制映射到丰度物理范围，对比多个时相时不会因为色彩拉伸不一致而产生视觉误导。
5.2 剖面线与端元残差验证
验证三分法是否真的有效，画一条穿过明显地表边界的剖面线，读取沿线像元的丰度值看变化是否与地表对应。植被和土壤的丰度应该在边界处陡变而不是渐变拖尾，阴影丰度应在地形阴影区域升高。
from matplotlib.ticker import MaxNLocator line_y = 300 # 剖面线所在的行号 profile = f_3d[line_y, :, :] # 取一行像元 fig, ax = plt.subplots(figsize=(10, 4)) ax.plot(profile[:, 0], label="vegetation") ax.plot(profile[:, 1], label="soil") ax.plot(profile[:, 2], label="shadow") ax.xaxis.set_major_locator(MaxNLocator(6)) # 控制横坐标标签数量 ax.set_xlabel("pixel along line") ax.set_ylabel("fraction") ax.legend() plt.tight_layout() plt.show()画这类剖面图时横坐标标签很容易挤在一起，MaxNLocator(6)把标签限制在 6 个刻度内，配合tight_layout能直接解决标签重叠问题。最后对照真实地表检查：耕地区域植被丰度应高于 0.7，裸土道路交叉口土壤丰度应接近 1，云影覆盖的林地阴影丰度应显著抬升而植被丰度下降但不过度归零。RMSE 图也是一份验证材料，RMSE 高值区不该集中在场景中央的地物交界处，而是集中在云边缘、水体波浪和阴影过渡带上，分布合理就说明像元三分法在这景影像上成立。
本文还有配套的精品资源，点击获取
