简介基于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 绿560nm10m植被中等反射、土壤中等反射区分阴影与亮度地物B4 红665nm10m植被强吸收分离植被与非植被的关键波段B8 近红外842nm10m植被高反射、土壤中等反射阴影端元与植被端元差异最大B11 短波红外1610nm20m对土壤湿度敏感用于构建超定方程组和残差验证选波段时要注意端元矩阵的条件数。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 环境是稳妥做法不碰系统的 Pythonconda create -n s2env python3.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求丰度矩阵 FN×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(axis1) # 有效像元才参与运算无效像元保持 NaN f_un np.full((flat.shape[0], 3), np.nan, dtypefloat32) 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, axis1, keepdimsTrue) 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, axis1)) 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里一个波段是 NaNnp.isfinite(flat).all(axis1)已经把该像元标记为无效但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, cmapYlGn, vmin0, vmax1) ax.set_title(title) ax.axis(off) plt.colorbar(im, axax, fraction0.046) plt.tight_layout() plt.savefig(abundance.png, dpi150) plt.show()vmin0, vmax1强制映射到丰度物理范围对比多个时相时不会因为色彩拉伸不一致而产生视觉误导。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], labelvegetation) ax.plot(profile[:, 1], labelsoil) ax.plot(profile[:, 2], labelshadow) 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 高值区不该集中在场景中央的地物交界处而是集中在云边缘、水体波浪和阴影过渡带上分布合理就说明像元三分法在这景影像上成立。本文还有配套的精品资源点击获取 Python实现Sentinel-2像元三分法：从线性光谱混合模型到植被覆盖度估算
简介基于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 绿560nm10m植被中等反射、土壤中等反射区分阴影与亮度地物B4 红665nm10m植被强吸收分离植被与非植被的关键波段B8 近红外842nm10m植被高反射、土壤中等反射阴影端元与植被端元差异最大B11 短波红外1610nm20m对土壤湿度敏感用于构建超定方程组和残差验证选波段时要注意端元矩阵的条件数。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 环境是稳妥做法不碰系统的 Pythonconda create -n s2env python3.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求丰度矩阵 FN×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(axis1) # 有效像元才参与运算无效像元保持 NaN f_un np.full((flat.shape[0], 3), np.nan, dtypefloat32) 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, axis1, keepdimsTrue) 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, axis1)) 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里一个波段是 NaNnp.isfinite(flat).all(axis1)已经把该像元标记为无效但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, cmapYlGn, vmin0, vmax1) ax.set_title(title) ax.axis(off) plt.colorbar(im, axax, fraction0.046) plt.tight_layout() plt.savefig(abundance.png, dpi150) plt.show()vmin0, vmax1强制映射到丰度物理范围对比多个时相时不会因为色彩拉伸不一致而产生视觉误导。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], labelvegetation) ax.plot(profile[:, 1], labelsoil) ax.plot(profile[:, 2], labelshadow) 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 高值区不该集中在场景中央的地物交界处而是集中在云边缘、水体波浪和阴影过渡带上分布合理就说明像元三分法在这景影像上成立。本文还有配套的精品资源点击获取