简介面向遥感与地信方向学习者这套基于Python的Sentinel-2卫星数据像元三分法模型资源以课程设计形式完整呈现高光谱遥感影像的处理流程。实现中引入最大噪声比变换对影像进行连续主成分分析式压缩与降维并借助像元纯度指数衡量每个像元的纯净程度为混合像元的三分法分解提供关键输入。资源共33个文件以8个Python脚本为核心覆盖波段读取、NDVI与DFI计算、数据融合、像元纯度指数计算等关键模块22张PNG图片直观展示中间结果与最终模型输出另附README说明文档与License文件目录组织清晰便于按需查阅和本地调试。压缩包整体仅4.01MB轻量易得适合遥感课程设计、毕业设计或初入像元分解领域的读者上手。目前已有440人学习下载借助源码与结果图对照学习可快速掌握Sentinel-2数据上MNF降维、PPI计算到三分法建模的完整链路并可直接复用或二次开发。1. 像元三分法模型是什么用 Python 啃下 Sentinel-2 混合像元这块硬骨头基于 Python 使用 Sentinel-2 卫星数据的像元三分法模型解决的是遥感里最日常的一个问题一个 10 米像元里同时压着稀疏灌丛、裸土和一点阴影硬分类器只能给它贴一个“灌丛”标签而三分法给出的是“灌丛 0.45、裸土 0.52、阴影 0.03”这样的连续占比。它的数学内核是线性光谱混合模型像元反射率约等于三类端元光谱按丰度加权求和。这篇文章用一套可直接改参数跑的代码把数据预处理、端元提取、非负最小二乘求解、GeoTIFF 输出和精度验证全流程串起来。适合正在做荒漠化监测、秸秆覆盖、水体与沉水植被识别的从业者哪怕你刚跟着 Python 零基础教程装好环境先拿一小块样区跑通再铺整景影像也完全能跟上。2. 数据准备Sentinel-2 波段选择、预处理与端元光谱提取三分法模型的效果上限在数据下锅之前就已经定了一大半。这一章不讲空泛的遥感原理只讲我用 Sentinel-2 做像元三分法时波段怎么挑、预处理脚本怎么写、端元光谱从哪儿来。2.1 波段选哪几个、要不要大气校正先把输入定死三分法的输入是每个像元的一条光谱曲线所以波段不是越多越好而是越“可区分”越好。我一般从 Sentinel-2 Level-2A 影像里取 6 个波段B2(蓝光 490nm)、B3(绿光 560nm)、B4(红光 665nm)、B8(近红外 842nm)、B11(短波红外 1610nm)、B12(短波红外 2190nm)。不取满 13 个波段的原因很直接B1 海岸气溶胶和 B9 卷云波段主要服务大气参数反演对地物丰度贡献小B5、B6、B7 三个红边波段在植被精细分类里有用但在三端元框架里很容易和 B8A、B11 强相关端元矩阵接近病态丰度解会被噪声放大。新手一上来把全部波段塞进最小二乘结果往往是一张雪花噪声丰度图——波段多不等于精度高这是我反复看到的第一类翻车。大气校正是另一个前置决定。如果下载的是 Level-1C 的 TOA 反射率先把影像过一遍 Sen2Cor命令行执行L2A_Process --res10得到 Level-2A 地表反射率产品。直接拿 L1C 做像元三分法蓝光和近红外的比值会整体偏移水、植被、裸土三个端元在特征空间里被压缩成一条线丰度系统性出错。Sen2Cor 的默认参数对中低纬度多数场景够用处理完检查新增的 SCL 场景分类波段是否完整覆盖研究区再往下走。分辨率上B2/B3/B4/B8 是 10mB11/B12 是 20m。常见做法是全部重采样到 10m20m 波段用双线性插值上采样云掩膜波段用最近邻上采样避免类别边缘被磨糊。投影统一到研究区所在 UTM 分带中国东部常用 EPSG:32650西部用 32645/32646跨带研究区全图统一到一个分带即可。2.2 用 rasterio 写一个预处理脚本重投影、重采样与云掩膜波段命名各家平台略有差异我在 ESA Copernicus 数据空间下载的 L2A 产品命名形如T50TMK_20230815_L2A_B04.tif。下面这段脚本把 6 个波段统一重采样到 B04 的 10m 网格并保存成一个多波段堆叠文件。import numpy as np import rasterio from rasterio.warp import reproject, Resampling # 以 B04 的 10 m 网格作为统一网格 REF T50TMK_20230815_L2A_B04.tif BAND_NAMES [B02, B03, B04, B08, B11, B12] OUT_STACK stack_6band.tif with rasterio.open(REF) as ref: profile ref.profile.copy() profile.update(count6, dtypefloat32, nodata-1) height, width ref.height, ref.width ref_transform, ref_crs ref.transform, ref.crs bands [] for name in BAND_NAMES: with rasterio.open(fT50TMK_20230815_L2A_{name}.tif) as src: arr np.zeros((height, width), dtypefloat32) reproject( sourcerasterio.band(src, 1), destinationarr, src_transformsrc.transform, src_crssrc.crs, dst_transformref_transform, dst_crsref_crs, resamplingResampling.bilinear) bands.append(arr) with rasterio.open(OUT_STACK, w, **profile) as dst: for i, arr in enumerate(bands, 1): dst.write(arr, i)这段脚本有点绕但对新手友好reproject同时完成投影转换和重采样src 是 20m 波段dst 是 10m 网格双线性插值会补充过渡值nodata-1是必要的否则无数据区会以 0 参与计算把水体端元的光谱均值拉低。如果影像本身已经是同一 UTM 投影只是分辨率不一致用arr src.read(1, out_shape(height,width), resamplingResampling.bilinear)会更快但跨投影场景老实走reproject。云和云影不处理掉三分法会把云边缘像元错分成水体。Sentinel-2 L2A 自带 SCL 波段类别码含义是 0 无数据、1 饱和、3 云影、4 植被、5 裸土、6 水体、7 未分类、8 中概率云、9 高概率云、10 薄卷云、11 雪。掩膜脚本如下with rasterio.open(T50TMK_20230815_L2A_SCL.tif) as src: scl src.read(1, out_shape(height, width), resamplingResampling.nearest) cloud_mask np.isin(scl, [0, 1, 3, 7, 8, 9, 10, 11]) valid_mask ~cloud_maskSCL 本身是 20m重采样到 10m 必须用最近邻保持类别边界不产生过渡值。把 7 归入无效是我比较保守的习惯如果确认研究区里大范围未分类是干净裸地可以把 7 从列表里摘掉。2.3 端元光谱提取NDVI 极值像元与端元三角形端元光谱矩阵 E 是三分法的心脏。端元有三个来源影像内部纯净像元、野外光谱仪实测、公共光谱库。对大多数监测项目我推荐第一种——影像内端元与传感器定标、大气状态自洽省去光谱库与影像之间的匹配误差这也是没有地面光谱仪时最可靠的方案。纯净像元怎么找三分法里常见做法是用 NDVI 和短波红外阈值卡出三类地物的“纯像元”再取光谱均值作为端元。stack np.stack(bands, axis-1) # (height, width, 6) red stack[..., 2] nir stack[..., 3] swir1 stack[..., 4] ndvi (nir - red) / (nir red 1e-10) # 植被端元NDVI 高短波红外明显低于近红外 veg_p (ndvi 0.6) (swir1 0.3 * nir) # 裸土端元NDVI 低短波红外通道较亮 soil_p (ndvi 0.2) (swir1 0.15) (swir1 0.45) # 水体端元近红外反射率很低 water_p (nir 0.08) (ndvi 0.05) E np.array([ stack[veg_p valid_mask].mean(axis0), stack[soil_p valid_mask].mean(axis0), stack[water_p valid_mask].mean(axis0) ]).T # 形状 (6, 3)这里的 0.6、0.2、0.08 是经验参考值不同生态区必须重调。调法不是拍脑袋画 B8 和 B11 的二维散点图看三团像元是否分得开阈值卡在密度谷上。如果某个端元的纯净像元数少于 500说明阈值太紧放宽后重取若阈值太松混入过渡像元解出的丰度会出现成片低估。valid_mask在这里是第二次保险把云影像元排除在端元统计之外。3. 核心求解用 numpy/scipy 把三分法写成逐像元丰度分解端元矩阵就绪后核心工作是把每个像元的六波段光谱拆成三个端元的丰度组合。这一章的公式和代码是全文基础后面所有输出、验证、避坑都挂在它上面。3.1 线性光谱混合模型的数学形式先写对公式再写代码像元三分法模型的数学形式是一个带约束的线性方程组。对第 i 个像元ρ_b Σ_{j1..3} f_j · ρ_{j,b} ε_b b 1..6ρ_b 是该像元在第 b 波段的反射率ρ_{j,b} 是第 j 个端元在第 b 波段的端元光谱值f_j 是第 j 个端元的丰度ε_b 是残差。加上两个物理约束f_j ≥ 0且 Σf_j 1。矩阵形式是ρ E · f εE 是 6×3 端元矩阵f 是 3×1 丰度向量。几何上理解更直观三个端元在特征空间里围成一个三角形任何一个混合像元只要落在三角形内部它的丰度就是相对端元的重心坐标。这个直觉后面很有用——如果你画出解出来的丰度大量在 0 到 1 之外说明像元跑到了三角形外要么端元选得不对要么影像里还有未清除的云影。反演前先想通这一点排起错来会快很多。3.2 用 scipy.optimize.nnls 求非负解再归一化求解方法的选择决定了你会不会在丰度图上看到负数。np.linalg.lstsq是无约束最小二乘解出的 f 经常出现负值因为它在数学上不保证物理意义。遥感界更稳妥的常规做法是带非负约束的最小二乘——scipy 的nnlsfrom scipy.optimize import nnls def unmix_pixel(pixel, E): # pixel: (6,) 一个像元的六波段反射率 # E: (6, 3) 三列分别是植被/裸土/水体端元 frac, _ nnls(E, pixel) total frac.sum() if total 0: frac frac / total return fracnnls只保证非负不保证和为 1所以归一化是必须的。为什么不用同时约束“非负且和为 1”的优化scipy 里没有直接可用的单函数接口常见做法就是nnls加归一化或者用scipy.optimize.minimize带线性约束去解后者慢一个数量级。nnls加归一化的代价是如果某个像元和三个端元都不像归一化会把误差摊给三个端元但后续残差图会把这个像元暴露出来这正是我们想要的诊断信息。3.3 逐像元求解的三种写法从死循环到分块并行最直觉的写法是逐像元循环但我得先劝退你一个大区域动辄上千万像元Python 层循环调用nnls能跑两小时。实际工程里我用三种写法应对不同阶段。第一种纯循环只用于几万个像元的调试def unmix_slow(img_2d, E): out np.zeros((img_2d.shape[0], E.shape[1]), dtypefloat32) for i in range(img_2d.shape[0]): out[i] unmix_pixel(img_2d[i], E) return out第二种伪逆加截断速度最快但精度略差适合快速预览def unmix_fast(pixel, E): frac np.linalg.pinv(E) pixel # 最小二乘闭式解 frac np.clip(frac, 0, None) # 截断负值 total frac.sum() return frac / total if total 0 else fracpinv在 6×3 维矩阵上是毫秒级运算全图一次矩阵乘就完成比逐像元nnls快几十倍。但截断负值等于人为塞回了信息丰度会整体高估一点所以它只能当“预览模式”先把端元调顺、把大问题暴露出来最后正式出图再换nnls。第三种分块加多进程这是正式出图的推荐写法from concurrent.futures import ProcessPoolExecutor def unmix_block(block, E): out np.zeros((block.shape[0], E.shape[1]), dtypefloat32) for i, p in enumerate(block): f, _ nnls(E, p) s f.sum() out[i] f / s if s 0 else f return out def unmix_raster(img, E, chunks2048): h, w, b img.shape flat img.reshape(-1, b) results [] with ProcessPoolExecutor(max_workers4) as ex: jobs [ex.submit(unmix_block, flat[i:i chunks], E) for i in range(0, flat.shape[0], chunks)] for job in jobs: results.append(job.result()) return np.vstack(results).reshape(h, w, E.shape[1])chunks2048表示每个子任务处理 2048 个像元一个块大约 48KB 内存进程之间通过返回结果汇总内存占用可控。max_workers4建议按 CPU 核心数减一设避免把机器拖死。在 Windows 上要把调用入口包进if __name__ __main__:否则多进程会重复加载主模块而报错。三种写法的取舍我做成了表写法速度约束适用场景逐像元循环慢非负归一化小样区调试pinv截断快无严格约束快速预览、调端元nnls分块多进程中等非负归一化整景正式出图4. 从丰度矩阵到 GeoTIFF结果输出、可视化与精度验证模型跑完只得到三个 numpy 数组离“能交付的成果图”还差两步写成带坐标的 GeoTIFF再做残差和回归验证。这一章把这两步讲透。4.1 把三张丰度图写成 GeoTIFF保留原始坐标信息输出 GeoTIFF 最省事的办法是复用参考波段的 profile坐标系、分辨率、数据范围全都不用重写with rasterio.open(REF) as ref: profile ref.profile.copy() profile.update(count3, dtypefloat32, nodata-1) with rasterio.open(abundance_fraction.tif, w, **profile) as dst: dst.write(frac[..., 0], 1) # 植被丰度 dst.write(frac[..., 1], 2) # 裸土丰度 dst.write(frac[..., 2], 3) # 水体丰度这里有两个细节丰度是 0 到 1 的连续值必须存 float32存成 uint8 会直接把小数抹掉nodata-1是给自己留的后路云掩膜区域保持无值将来统计面积时可以按有效像元过滤。写入后可以用 QGIS 或 GDAL 快速叠加原始真彩色影像抽查。4.2 精度验证先算模型残差再做回归精度验证分两层。第一层是模型内符合度看每个像元的拟合 RMSEh, w, _ frac.shape flat_frac frac.reshape(-1, 3) flat_img img.reshape(-1, 6) pred flat_frac E.T # 用丰度反算六波段反射率 rmse np.sqrt(np.mean((flat_img - pred) ** 2, axis1)) rmse_map rmse.reshape(h, w)反射率的量级在 0 到 0.5 之间RMSE 小于 0.02 说明三个端元把像元解释得很好如果 RMSE 在某个区域成片偏高先怀疑云影没清干净或者端元矩阵漏掉了不透水面而不是急着改算法。这一步能拦下 70% 的交付事故我每次跑完都先看残差图再看丰度图。第二层是外部验证。常见做法是拿无人机正射影像或 0.3m 高分影像在研究区随机布 30 到 50 个样点人眼目视解译每个样点的三类地物占比与模型丰度做线性回归。R² 大于 0.7 算可接受大于 0.8 算良好。没有高分辨率影像时至少要安排两个人独立解译同一批样点算一致性系数把主观误差量出来。4.3 可视化三色合成、直方图与像元散点图丰度图最终要给人看可视化不是简单打印而是检验结果的一种手段import matplotlib.pyplot as plt rgb np.stack([frac[..., 0], frac[..., 1], frac[..., 2]], axis-1) plt.imshow(rgb, vmin0, vmax1) plt.axis(off) plt.savefig(abundance_rgb.png, dpi300, bbox_inchestight) plt.figure() for i, name in enumerate([植被, 裸土, 水体]): plt.hist(frac[..., i].ravel(), bins100, alpha0.5, labelname) plt.legend() plt.savefig(abundance_hist.png, dpi300)三色合成把植被通道映射到红色、裸土到绿色、水体到蓝色空间格局一眼可见河流应该变成连续稳定的一条纯蓝带农田里应该是红绿过渡而不是马赛克。直方图用于检查丰度分布有没有大量像元顶到 1 或堆在 0 附近——单峰堆在两端说明端元阈值卡得太死把过渡像元全推成了纯像元。5. 像元三分法避坑五个高频翻车现场与排查思路三分法的代码量不大但坑都藏在数据和参数里。我把过去几年带新人时最常见的五个翻车现场整理成“现象—原因—解决”三段式每一条都是真实改过代码的教训。5.1 丰度出现大量负数或超过 1现象输出丰度图里随机出现 -0.3、1.4 这类数字。原因用了np.linalg.lstsq或np.linalg.pinv直接解线无约束最小二乘在数学上会允许任何实数解还有可能是两个端元光谱太接近端元矩阵接近奇异解对噪声高度敏感。解决正式结果一律换成 3.2 节的nnls加归一化同时检查端元矩阵条件数如果两个端元光谱的欧氏距离小于 0.02说明端元几乎平行结果不可信。这时要么换波段组合要么把两个近似端元合并后重新建模。5.2 模型残差集中在近红外和短波红外现象全局 RMSE 不高但单独看 B8 和 B11 两个波段的残差比其他波段高出一倍。原因影像还是 Level-1C 的 TOA 反射率没有做大气校正。近红外受气溶胶影响最大B11 又对水汽敏感两个波段同时偏高几乎是大气校正缺失的指纹。解决确认输入数据是 L2A不是 L1C如果只有 L1C先跑 Sen2Cor。处理完后重新提取端元并再次检查残差B8 和 B11 的残差会明显回落。5.3 端元阈值在另一个区域完全不适用现象同一套 NDVI 阈值代码在 A 区效果正常换到 B 区水体丰度几乎全为 0植被丰度普遍偏高。原因0.6、0.2、0.08 是按 A 区冬季影像调的B 区夏季植被更密、裸土更亮端元像元落在阈值区间之外。解决每换一个区域先画 B8 对 B11 的二维散点图看三个端元聚类的密度谷在哪再把阈值卡到谷上。我的习惯是做一个交互脚本在散点图上手动点选三个三角形的角点作为端元比反复调阈值快得多也少很多玄学。5.4 云掩膜没做膨胀云边像元丰度全线异常现象云边界往外一个像元出现水体丰度骤增或裸土丰度骤减的亮环。原因SCL 波段是 20m 分辨率重采样到 10m 后云边缘的混合像元仍保留在有效区SCL 的 7 未分类码里还可能混着薄云。解决对云掩膜做一次膨胀把云边界往外推一格用scipy.ndimage.binary_dilation(cloud_mask, iterations1)然后在端元提取和模型求解里统一用膨胀后的valid_mask。这条操作花十秒钟能省掉半天手动掩膜工作。5.5 全图逐像元循环跑了两个小时没出结果现象小样区秒完整景影像跑了几小时内存占用涨到几十 GB。原因一次性把全图读成 float64 数组再用单进程 for 循环调nnls。nnls本身是 C 实现但 Python 层逐像元调用和内存拷贝的开销是主要瓶颈。解决读入时指定 dtypefloat32数组体积直接减半用 3.3 节的分块多进程写法20km×20km 的范围在四进程下通常能压到十几分钟。如果急着看结果先用 pinv 截断法出预览图调好端元后再用nnls出正式图。6. 进阶方向端元数量、时间序列与验证方法的取舍三分法模型不是只能写死三个端元。当研究区里出现大面积不透水面或山体阴影时把阴影硬塞给水体或裸土会带来系统性误差。常规做法是把端元数量扩到四个或五个比如植被、裸土、水体、阴影、不透水面端元矩阵变成(6, 5)公式和求解代码不用动丰度图多出两个波段。代价是端元之间的共线性会变强结果对噪声更敏感残差图要逐个波段检查。经验法则是端元数量不要超过波段数的一半六个波段撑死塞四个端元多于这个比例就属于过拟合。另一个值得投入的方向是时间序列三分法。用同一季节的逐年 Sentinel-2 影像生成丰度序列然后做 Theil-Sen 斜率和 Mann-Kendall 趋势检验。这比单期丰度绝对值可靠得多——单期影像可能残留大气和土壤水分干扰而连续几年的相对变化能把这些噪声抵消掉。荒漠化监测项目里趋势图比丰度图更受业务方认可。验证方法上我的固定顺序是先看全局 RMSE 和残差空间分布残差图干净了再布 30 到 50 个目视解译样点做回归最后把区域统计量和统计年鉴或已有专题数据对一下量级。这三件事都闭环结果才敢交付。我自己踩过最重的坑是端元光谱总想一次定死后来养成“先 pinv 快速出预览图、再调端元、看残差图、最后 nnls 出正式图”的习惯后返工率低了很多。这套流程同样适用于你手里的研究区希望帮到你。本文还有配套的精品资源点击获取