上个月帮朋友处理了一景哨兵二号影像目标是给某研究区做土地覆盖分类再顺带用机器学习回归反演一下植被叶绿素相对含量。整个过程从原始影像下载、波段处理、样本标注到随机森林分类、精度评估再到定量反演全部用Python走通。这条链路其实是目前遥感数据应用里最典型的一种需求组合先做图像分类再做定量评估中间穿插各种机器学习方法。这篇文章就把这条链路完整复盘一遍把每一步的关键参数、代码逻辑、踩坑记录都留下来给正在做多光谱遥感数据处理、图像分类或者想入门机器学习遥感应用的朋友做个参考。1. 案例立项多光谱遥感到底在解决什么问题1.1 多光谱和“肉眼可见光”差在哪里很多人第一次接触多光谱数据第一反应是“不就是彩色照片加几个通道吗”。如果拿普通数码相机和哨兵二号做对比人眼和相机只有红绿蓝三个波段而哨兵二号有十三个波段覆盖了可见光、红边、近红外和短波红外。关键不在于波段多而在于波段的“位置”。近红外波段对植被非常敏感因为健康植物叶片的细胞结构会强烈反射近红外而红光波段会被叶绿素强烈吸收。这两个波段一组合能算出归一化植被指数NDVI直接反映植被长势。这就是多光谱数据最核心的价值——它能看见肉眼看不见的光谱信息而光谱差异正是区分地物类别的物理基础。在做土地覆盖分类时耕地的玉米、林地、草地、水体、建筑在不同波段上的反射率曲线差异非常明显。水体在近红外和短波红外几乎不反射建筑刚好相反植被在红光低、近红外高的“红边”特征突出。这些光谱形态叠加纹理特征就是机器学习分类器判断地物类型的依据。1.2 为什么这个案例要选Python做遥感数据处理可选的工具链其实不少。传统上很多人用ENVI、ERDAS或者ArcGIS的栅格计算器这些软件在交互式操作上确实方便但遇到批量处理、参数调整、自动化流程效率就明显不够了。Python在这个领域的优势很直接rasterio和GDAL负责读写栅格numpy负责数组计算scikit-learn负责机器学习模型matplotlib画图geopandas处理矢量边界。整条数据管线都能在一个脚本里串起来参数随时改结果可复现。而且遥感影像本质就是带地理坐标系的多维数组用数组思维去处理比鼠标点来点去高效得多。对于后续要扩展到深度学习、批量处理大量影像的场景Python几乎是唯一合理的选择。1.3 案例整体流程设计这个案例的完整流程我拆成了七个环节环境准备、数据获取与预处理、样本标注与特征工程、分类模型训练、分类后处理、定量评估与精度验证、结果分析与制图。其中分类和评估是两条技术线并行。图像分类解决的是“这是什么地物”的离散问题用监督分类模型完成定量评估解决的是“植被长势怎么样”的连续问题用回归型机器学习模型完成。这两条线共用同一套多光谱数据但在样本形式、模型选择、评价指标上差别很大。这也是很多初学者容易混的地方先把这两条线的区别想清楚后面的流程就顺了。2. 环境搭建与核心工具选型2.1 用conda而不是直接pip遥感库的依赖关系比普通web开发复杂得多rasterio要依赖libgdalGDAL本身又是出了名的“安装杀手”。直接用pip install经常出现版本冲突、编译失败的报错。我的建议是先装Miniconda然后新建独立的虚拟环境让Python版本固定在3.9或3.10。为什么不用最新版本因为很多深度学习或者地理空间库对Python 3.12以上的支持在早期并不完善遇到“No matching distribution found”的错误很折腾。conda环境的好处是隔离不同项目用不同环境装坏了删掉重来不会把系统的Python搞坏。提示如果电脑里已经装了ArcGIS或其他GIS软件千万不要随便动系统级Python很可能导致软件自带的Python环境崩溃。用conda环境隔离是最安全的做法。2.2 必须装的库和版本建议我这次环境里装的库直接列出来供参考python3.9rasterio1.3.x读写GeoTIFF、JP2影像numpy1.24.x数组运算pandas样本数据管理scikit-learn1.3.x分类与回归模型scikit-imageGLCM纹理特征计算matplotlib结果出图geopandas0.14.x矢量边界读取、裁剪joblib模型保存与加载tqdm进度条处理大影像时很有用安装语句很简单conda create -n rsenv python3.9创建好环境后conda activate rsenv再逐个conda install -c conda-forge 包名。conda-forge这个频道对地理空间库维护得很及时优先用它。2.3 开发环境配置经验有些朋友会卡在“Anaconda装好了但IDE里import rasterio失败”这个问题上。核心是解释器路径选错了。在PyCharm里需要到Settings - Project - Python Interpreter选择Existing environment路径定位到conda环境里的python.exe比如C:\Users\用户名\anaconda3\envs\rsenv\python.exe。VSCode里则是在命令面板敲Python: Select Interpreter选择对应的conda环境。选对之后终端里敲python再import rasterio如果能正常通过说明环境通了。还有一个很实用的调试技巧——在脚本开头加上这几行能把环境信息和当前工作路径先打出来避免后续因为路径不对而报“file not found”import sys, os sys.path.append(os.path.dirname(__file__)) print(Python版本:, sys.version) print(当前目录:, os.getcwd())3. 数据获取与预处理决定分类上限的第一步3.1 数据源选择与波段对应关系哨兵二号L2A产品已经做过大气校正反射率数据可以直接用省去很多大气校正的麻烦。如果用的是L1C数据还需要做大气校正这里不展开但记住一点分类模型训练时用的是地表反射率还是大气顶反射率会直接影响结果的可迁移性L2A是更稳妥的选择。哨兵二号的波段分辨率不统一四个可见光波段和近红外是10米红边和短波红外是20米气溶胶和水汽波段是60米。实际使用时我一般保留十个核心波段B2、B3、B4、B8、B5、B6、B7、B11、B12再加上一个可能用到的B1气溶胶波段。把他们全部重采样到10米统一坐标系这样每个像元都对应一组完整的光谱特征向量。3.2 波段读取、裁剪与重采样先解决最基础的问题怎么把影像波段读进来。哨兵二号L2A数据下载下来是.SAFE文件夹格式里面是JP2格式的影像。用rasterio读取的时候需要定位到GRANULE子目录下的IMG_DATA文件夹。我习惯把需要的波段全部读取成一个多维数组然后保存成自己拼接的多波段GeoTIFF。这样后续每次只要读取一个文件不用反复定位原始目录。用rasterio实现波段拼接和重采样代码大概是这样import rasterio from rasterio.warp import calculate_default_transform, reproject from rasterio.enums import Resampling band_paths [ IMG_DATA/B02_10m.jp2, IMG_DATA/B03_10m.jp2, IMG_DATA/B04_10m.jp2, IMG_DATA/B08_10m.jp2, # 其他波段... ] # 读取每个波段并堆叠成数组 le_bands [] with rasterio.open(band_paths[0]) as src: profile src.profile for bp in band_paths: with rasterio.open(bp) as src_band: le_bands.append(src_band.read(1)) # 转成 (波段, 行, 列) 数组 stack np.stack(le_bands, axis0) # 写入多波段GeoTIFF with rasterio.open(subset_stack.tif, w, **profile) as dst: dst.count stack.shape[0] dst.write(stack)如果是做研究区裁剪先读入shapefile边界再用rasterio.mask.mask函数裁剪。注意要处理坐标系一致性如果shapefile是WGS84经纬度影像又是UTM投影需要先to_crs统一坐标。这一步忘了做mask出来的结果常常是空的。3.3 光谱指数计算给机器学习加“先验知识”机器学习模型虽然能自己从原始波段里学习特征但光谱指数相当于告诉模型“我已知哪个波段组合对植被敏感”这是把遥感领域的先验知识注入特征空间。实际测试下来加入NDVI、NDWI、SAVI之后分类精度通常会有几个百分点的提升。我这次计算了一组指数NDVI用红光和近红外NDWI用绿光和近红外SAVI加了土壤调节因子EVI增强植被指数用蓝光做气溶胶校正。以NDVI为例计算逻辑就是波段数组的逐像元运算import numpy as np def ndvi(nir, red): # 避免分母为0 return (nir - red) / (nir red 1e-10) # 假设 B8 是近红外B4 是红光 ndvi_arr ndvi(stack[3].astype(float), stack[2].astype(float))需要注意数据类型精度。JP2原始数据是uint16如果直接用整数类型做差值计算负数就会被截断算出来全错。必须先把波段转为float32再做运算这是多光谱数据处理里最容易踩的坑。预处理做完整套每个像元就有一个高维特征向量十个波段反射率、五个光谱指数再叠加后面要提到的纹理特征。这一步是整个流程的地基数据质量决定了机器学习分类精度的上限后面模型再怎么优化都是在这个上限之内逼近。4. 样本标注与特征工程4.1 样本怎么做最省力训练样本的质量比数量重要。我在QGIS里把真彩色影像加载进去叠加研究区的实地调查点然后手动绘制多边形ROI分别标上类别标签林地、耕地、草地、水体、建设用地。这是图像分类任务里最经典的样本生产方式也叫感兴趣区标注。样本量参考经验值每类至少五百到一千个像元。这里的“像元”不是样本点而是每个标注多边形内部的所有像元。比如画一个五十像元的耕地多边形这五十个像元都是耕地样本。标注时要注意均匀覆盖不同光谱形态不要只在影像某个角落画样本否则模型学到的只是局部光谱特征换一个区域就得重新标注。制作样本的时候直接在标注多边形范围内提取每个像元的特征向量并把对应的类别标签存成一张表。这一步的输出是一份带标签的DataFrame每一行是一个像元样本每一列是一个特征。后续训练模型吃的就是这份数据。4.2 纹理特征怎么加光谱特征只能区分不同地物在谱面上的差异但建筑和裸地在某些波段上可能长得差不多这时候纹理信息就派上用场了。比如建设用地通常是粗糙的纹理水面是平滑的纹理。我用的方法是灰度共生矩阵GLCM在scikit-image库里有现成实现。计算时先选择灰度化比较合适的波段一般是近红外波段因为它对地物结构敏感然后在一个窗口内统计像元灰度值的空间共现关系。from skimage.feature import graycomatrix, gray_levels import numpy as np # 取近红外波段 nir_band stack[3] # 转换为灰度值并归一化到8bit nir_u8 (nir_band / nir_band.max() * 255).astype(np.uint8) # 计算GLCM距离1方向0度这里用Cython底层并行会较快 def glcm_features(win_size5): from skimage.feature import local_binary_pattern # 备选特征 glcm graycomatrix(nir_u8, distances[1], angles[0], levels256, symmetricTrue, normedTrue) contrast np.sum(glcm * (np.arange(256)[:, None, None, None] - np.arange(256)[None, :, None, None])**2) return contrastGLCM的窗口大小需要试验。窗口太小纹理估不准窗口太大又会把边界模糊掉。我试过3x3、5x5、7x7最后5x5效果比较均衡。计算时注意用反射率归一化后的灰度图避免原始数值范围影响。4.3 训练集验证集划分必须分层样本划分有个容易被忽略的细节——遥感样本空间自相关性很强。同一块农田里的像元高度相似如果随机划分训练集和验证集验证集里可能全是训练集邻居模型评估结果虚高。正确的做法是分层抽样并且最好按多边形分组。也就是说保证训练集和验证集来自不同的样本多边形区域这样评估出的精度更接近真实业务效果。scikit-learn里的train_test_split默认是逐样本随机用它做遥感分类评估会乐观得离谱这个坑我踩过后面单独说。数据标准化也叫特征缩放也是不能漏的。SVM和KNN这类模型对特征尺度极其敏感十波段反射率值域可能是0到几千光谱指数是-1到1如果不做标准化模型会被数值大的特征主导小特征的有效信息被淹没。随机森林这类树模型对尺度不敏感但如果后续要对比SVM、KNN还是先统一做一遍StandardScaler更省心。5. 图像分类从随机森林到深度学习的对比5.1 为什么这个案例先选随机森林遥感图像分类里随机森林几乎是最稳妥的“起步模型”。它对高维特征不敏感不用做太多特征筛选能处理非线性关系不容易过拟合训练速度快还能量化每个特征的重要性方便做特征筛选。算法原理上随机森林就是构建一大批决策树每棵树用不同的随机子样本和随机特征子集训练最终分类结果由所有树投票决定。它吸收了一百多年前“三个臭皮匠顶个诸葛亮”这个朴素思想在多数遥感分类任务里随机森林的精度都能达到甚至超过SVM。用scikit-learn实现非常直接from sklearn.ensemble import RandomForestClassifier rf_model RandomForestClassifier( n_estimators300, max_depthNone, min_samples_split5, min_samples_leaf2, class_weightbalanced, n_jobs-1 ) rf_model.fit(X_train, y_train)参数解释一下n_estimators是树的数量越大越稳但超过一定阈值收益递减我一般取300到500min_samples_leaf控制叶子节点最小样本数设置成2可以有效抑制过拟合class_weightbalanced让每类的权重按样本量反比分配解决样本不平衡问题。5.2 不同模型横向对比光用一个模型说服力不够我这次一并跑了SVM、KNN和XGBoost做对比。SVM在小样本上精度不错但影像分类会把全研究区的每个像元都预测一遍像元数量动辄上千万SVM推理起来非常慢时间成本高。KNN实现简单但对特征尺度敏感预测时计算量也大。XGBoost梯度提升树精度高传统机器学习里最强的一类但参数多调参成本高。这几种模型在训练集的精度都接近完美但真正拉开差距的是验证集精度和推理速度。整体下来随机森林是最平衡的选择这也是遥感业务里它长期占据主流位置的原因。模型验证集精度训练速度推理速度调参难度结论随机森林高快中低首选SVM高中慢中小样本可用KNN中快很慢低不推荐大影像XGBoost更高快中高精度瓶颈期再上5.3 要不要上深度学习、Transformer类模型现在深度学习在遥感领域很火热词里也有transformer图像分类但实际落地要看样本量。CNN和Transformer类模型需要大量训练数据一套完整标注好的遥感分类数据集通常要几十万以上的样本。普通研究区人工标注几千个样本训练深度学习模型很容易过拟合效果反而不如随机森林。我的经验判断是先跑通传统机器学习流程建立基线精度。如果样本量足够大比如自动化样本生成或者有开源数据集再尝试深度学习模型做对比。图像分类模型更新迭代很快但业务场景中首先要求稳定、可控、可解释传统机器学习在这些方面优势明显。顺便回应一下热词里那个“计算机视觉和机器学习区别”的问题——在这个案例里用GLCM纹理特征加随机森林就是机器学习方法而端到端的CNN自动提取特征再分类就偏向计算机视觉了。前者特征工程靠人来设计后者特征自动学习。两者不是对立关系而是同一问题的不同路径。5.4 分类后处理椒盐噪声怎么去像素级分类输出最典型的问题是椒盐噪声——分类结果图上一颗颗孤立的小图斑看起来像撒了盐和胡椒。产生原因是单个像元光谱异常导致误分类比如阴影里的草地被分成林地。后处理我用的是众数滤波也叫majority filter。原理很简单对每个像元统计周围窗口内所有像元的类别频次把当前像元替换成窗口内的众数类别。窗口大小一般用3x3或5x5。虽然只动一个窗口但在去掉小图斑方面效果立竿见影。还可以用连通域分析把小于指定面积比如十像元的图斑合并到周围占比最大的类别里。这就更像正规业务逻辑了处理完的分类图会干净很多。后处理做完再去做精度评估通常精度值也会有轻微提升。6. 定量评估精度验证与回归反演6.1 混淆矩阵、总体精度、Kappa系数怎么算分类做完得用独立验证样本算精度。这里用到的评估指标包括总体精度、Kappa系数、各类别的生产者精度和用户精度它们都从混淆矩阵派生出来。sklearn计算非常方便from sklearn.metrics import accuracy_score, classification_report, cohen_kappa_score preds_train rf_model.predict(X_train) preds_test rf_model.predict(X_test) print(训练集精度:, accuracy_score(y_train, preds_train)) print(验证集精度:, accuracy_score(y_test, preds_test)) print(Kappa:, cohen_kappa_score(y_test, preds_test)) print(classification_report(y_test, preds_test, target_namesclass_names))classification_report里每个类别都有三个指标precision是预测为该类别的像元里有多少真属于该类recall是该类别真实像元里有多少被正确识别F1是两个的平衡。在多类别分类里还有一个宏平均和加权平均的区别宏平均不考虑样本量加权平均会按各真实类别样本量加权。遥感文献里一般同时报告总体精度和KappaKappa值大于0.8说明一致性很好。一个重要的检查项训练集精度和验证集精度之差。如果差值超过五个百分点说明过拟合了。这时候需要增加样本量、降低模型复杂度、加大正则化参数或者检查特征中是否混入了泄漏数据。6.2 像元尺度上的面积统计分类得到的是每个像元的类别标签号做定量评估时还要统计各类别面积。方法并不复杂统计每个类别的像元个数乘以单像元面积得到各自面积再算面积占比。# 每个像元的面积哨兵二号10米分辨率 pixel_area 10 * 10 # 100平方米 unique, counts np.unique(class_map, return_countsTrue) class_area counts * pixel_area total_area class_area.sum() for class_id, area_m2 in zip(unique, class_area): print(f类别{class_id}: {area_m2:.2f}平方米, 占比{area_m2/total_area*100:.2f}%)在业务中还会用面积变化做多期对比比如同一区域的六月份和九月份分别分类定量评估耕地面积、水体面积的变化量。这就是遥感图像分类在农林监测、生态评估中的常见应用方式。6.3 回归型定量评估机器学习反演植被参数除了分类还有一类定量评估任务不是判断地物是什么而是估算“多少”。比如叶绿素含量、叶面积指数LAI、土壤含水量。这类问题用回归模型解决输出的是连续数值而不是类别标签评价指标变成R平方、均方根误差、平均绝对误差。我这次用回归方式反演了研究区植被的相对叶绿素含量。用实验仪器在实地测了一批叶绿素值然后与影像上同一像元的光谱特征做关联训练一个回归模型。这里的数据组织方式和分类任务类似但标签是连续值。from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error reg_model RandomForestRegressor(n_estimators300, n_jobs-1) reg_model.fit(X_train, y_train_reg) preds_reg reg_model.predict(X_test) print(R2:, r2_score(y_test_reg, preds_reg)) print(RMSE:, np.sqrt(mean_squared_error(y_test_reg, preds_reg))) print(MAE:, mean_absolute_error(y_test_reg, preds_reg))R平方解释模型解释的方差比例0.8以上基本可以接受RMSE有量纲比如叶绿素是0到60的量程RMSE为8说明平均偏差约8个单位这个绝对数值能帮你判断模型误差是否在业务容忍范围内。反演模型训练好之后把研究区所有像元的特征灌进去输出一幅连续数值栅格图就能直观看到研究区叶绿素含量的空间分布。这整个流程就是标题里的“定量评估”和“机器学习方法应用”结合最紧密的部分。7. 常见报错与排查经验7.1 读遥感影像的经典报错处理遥感数据最常见的报错是路径和格式问题。.SAFE目录下的JP2文件扩展名是.jp2如果不指定驱动rasterio有时会报“Unrecognized dataset driver”。解决办法是显式注册驱动或者在读取前确认文件路径指向的是具体影像文件而非.SAFE文件夹。还有一类报错是投影坐标系不一致。比如用rasterio.mask.mask裁剪时掩膜矢量与影像投影不同结果可能全黑。解决办法是先用shapefile.to_crs(src.crs)转换坐标。这个检查看起来简单但实际中经常会漏掉尤其在不同来源数据叠加时。7.2 内存溢出怎么办哨兵二号研究区稍微大一点波段数组就可能到几亿像元。把所有波段一次性用np.stack堆叠可能直接内存溢出。我的做法是分块处理把影像切割成小块patch逐块计算特征并保存到磁盘最后再合并。rasterio支持窗口读取可以只读某个行列范围内的数据。这样内存占用固定在一个小patch的规模不会随着影像大小线性增长。实测下来处理一景完整的哨兵二号大区域分块处理能省掉大半内存峰值。对于整图预测分类结果也是按块预测后拼接不会爆内存。7.3 模型训练不收敛或者过拟合很多人在遥感分类里碰到的不是不收敛而是过拟合。现象是训练集精度接近百分之百验证集精度突然掉到七十多甚至更低。最常见的原因是前面提到的样本空间自相关。解决方法是按多边形分组划分训练集而不是逐像元随机划分。另一个原因是特征过多且相关性太强可以看随机森林的特征重要性把重要性极低的特征删除或者做主成分分析降维。还有一个原因是模型过于复杂比如随机森林的树深度不限制叶子节点样本数太少这时候调整max_depth和min_samples_leaf能明显改善过拟合。7.4 处理效率太慢的经验整景影像预测类别时把所有像元特征堆成一个大矩阵输给模型推理阶段可能要跑很久。sklearn的预测是单线程的随机森林的n_jobs-1在训练时有效预测时不一定全核并行。加速的办法是减少特征维度、使用批处理预测或者如果容器内存足够直接利用scikit-learn的check_inputFalse参数减少预处理开销。对于特别大的影像更务实的方法是先用小范围测试流程全通再全图预测。千万不要一上来就对整景影像跑报错排查成本太高。8. 实操心得这套流程还能怎么扩展写到这里把整个流程串起来看核心是“特征构建传统机器学习独立验证”这个三角关系。图像分类和定量评估两类任务都建立在这套框架上区别只在标签的类型和评价指标。这个框架本身非常稳定只要数据预处理做好后续扩展深度学习方法也是顺理成章的。我再分享两个实际操作中很有用的习惯。第一每个环节都保留中间结果磁盘缓存尤其是拼接好的多波段影像、特征矩阵和训练好的模型。这样每次调试都不用从原始数据重新跑一遍能省掉大量时间。第二把整套流程写成函数封装比如prepare_features()、train_classifier()、evaluate_model()方便换数据源快速复用。这套流程在做完土地覆盖分类之后下一步很自然地就能扩展到时序分析。比如拿连续几个月或者几年的同区域影像用同一套样本和模型分别分类就能定量评估水体面积变化、耕地作物长势变化甚至结合气象数据做综合分析。做遥感数据处理的初期最大的坑不是模型不会写而是数据管线跑不通。把预处理步骤反复打磨把每个环节的输入输出都验证清楚后面加多少算法都是顺理成章的事。