简介这份《生态系统遥感解译与野外核查》资料面向从事生态监测、遥感解译与国土生态评估的技术人员及高校相关专业师生帮助解决生态状况调查中解译标准不统一、野外核查流程不规范的问题。内容围绕全国生态状况调查评估技术规范的编制说明展开涵盖项目背景与主要工作过程、标准制修订必要性分析、国内外相关标准对比、基本原则与技术路线以及适用范围、规范性引用文件、术语定义、遥感解译、野外核查与生态系统类型解译结果等主要技术内容并延伸到同类标准水平对比、实施管理措施与技术措施建议、环境效益及经济技术分析可帮助读者理解从数据获取、图像处理、解译到实地验证的完整链条。资源为1个PDF文件压缩包约265KB轻量便于随时查阅与打印。目前已有222人学习下载适合需要把握生态遥感监测标准化思路、撰写技术规范或开展野外核查工作的读者参考。1. 生态系统遥感解译与野外核查为什么必须成对做很多团队交出的第一版成果是一张生态系统类型分布图加一份精度评价表反馈回来的问题却是核查点的原始记录在哪、抽样方案怎么定的、为什么这几个类型一个样点都没有。生态系统遥感解译和野外核查在项目管理里常被拆成两件事一个靠算法一个靠出差一旦拆开做精度数字就成了没有证据链的孤证。解译是把影像的光谱、物候、纹理映射到一套类型体系上的推断过程核查则是用可重复的抽样和现场记录去约束这个推断的误差范围。前者决定你能分到多细后者决定你敢报多大的精度。做生态状况评估、自然保护地监测、碳汇本底调查、生态修复成效评估的工程师基本都会踩到这条线刚接手一份只有几百页工作方案 PDF、没写过一行分类代码的人也一样。往下走真正卡人的不是算法新旧而是分类体系、训练样本、抽样方案、验证口径这四件事能不能对齐。2. 解译底图与分类体系把口径定死再动算法2.1 生态系统类型体系怎么拆到可解译的最小单元分类体系不能直接从生态学教科书上抄。教科书按植被型、群系、群丛往下分一个群丛可能只有几公顷10 米分辨率的影像根本分不出来。落到工程里先把类型压到一级类森林、灌丛、草地、湿地、农田、建成区、裸地/稀疏植被、水体这八类在多数区域能做到 80% 以上的可分离度。二级类要按能不能被特征区分来定而不是按科学上该不该分来定。比如常绿针叶林和落叶阔叶林在生长季峰值期的单时相影像上几乎不可分但加上冬季影像或者 NDVI 年振幅就能拉开草地和农田在 7 月的华北平原光谱高度重叠得靠地块边界规整度和收割前后的物候突变来分。最小上图单元MMU必须显式写进方案。10 米影像常用 0.1 公顷换算成 10×10 像元20 米影像通常放到 0.25 公顷。不写 MMU后处理去小斑块时就没有依据做出来的图斑数量会是实际管理单元的几倍。混合像元也要提前约定地类边界内外各 1 个像元不参与训练样本选取这一点在山区和破碎地块尤其要紧。一级类主要区分特征10m 影像建议 MMU单类训练样本下限森林NDVI 高、红边斜率大、纹理粗糙0.1 hm²80灌丛NDVI 中高、纹理细碎、物候振幅小0.1 hm²60草地生长季 NDVI 峰值低、年振幅大0.1 hm²60湿地NDWI 高、季节性淹水、植被混杂0.1 hm²50农田地块规整、收割期突变、多峰物候0.1 hm²80建成区NDBI 高、纹理规则、时序稳定0.1 hm²50裸地亮度高、NDVI 低且平稳0.25 hm²50水体NDWI 高、NDVI 负值0.1 hm²402.2 Sentinel-2 与 Landsat 的时相选择与云掩膜处理时相选择比影像源选择更影响结果。北方地区类型区分主要靠 6 到 8 月的生长季峰值合成落叶与常绿、农田与草地的区分则要补一期 11 月到次年 2 月的冬季影像。南方多云雨地区一个月的无云影像可能只有两三景得把窗口放到整个生长季做中值合成。下面这段在 GEE 里完成 Sentinel-2 L2A 的生长季去云合成用的是 SCL 分类层加云概率层双重掩膜。// GEESentinel-2 L2A 生长季无云中值合成 var roi ee.Geometry.Rectangle([110.2, 30.1, 111.8, 31.5]); var s2 ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED) .filterDate(2024-06-01, 2024-09-30) .filterBounds(roi) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 40)); var cld ee.ImageCollection(COPERNICUS/S2_CLOUD_PROBABILITY) .filterDate(2024-06-01, 2024-09-30) .filterBounds(roi); // 按 system:index 把云概率影像挂到 SR 影像上 var joined ee.Join.saveFirst(cloud).apply({ primary: s2, secondary: cld, condition: ee.Filter.equals({leftField: system:index, rightField: system:index}) }); // 双重掩膜云概率 40%同时剔除 SCL 中的云、云影、卷云 var maskCloud function(img) { var prob ee.Image(img.get(cloud)).select(probability); var scl img.select(SCL); var sclMask scl.neq(3) // 云影 .and(scl.neq(8)) // 中概率云 .and(scl.neq(9)) // 高概率云 .and(scl.neq(10)); // 卷云 return img.updateMask(prob.lt(40)).updateMask(sclMask); }; var composite ee.ImageCollection(joined).map(maskCloud) .select([B2,B3,B4,B5,B6,B7,B8,B11,B12]) .median() // 中值合成抗残留云 .multiply(0.0001) // 反射率缩放 .clip(roi); Export.image.toDrive({image: composite, description: eco_comp_2024, region: roi, scale: 10, maxPixels: 1e13});filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 40))这一条把整景云量过高的影像先筛掉阈值放到 40% 而不是 10%是为了在南方保住足够的影像数量剩下的像元级云由掩膜处理。ee.Join.saveFirst把云概率集合按system:index对齐到 SR 影像两套集合的 index 是一致的这个 join 只是把两个集合合并成一条记录不改变像元。SCL 层里剔掉的四类对应云影、中概率云、高概率云、卷云水体6和雪11没有剔如果研究区有季节性积雪需要额外加条件。中值合成对残留薄云比均值合成稳健得多代价是物候信息被压平所以时序特征要单独从原始影像集合里算不能从这个合成结果里取。Landsat 8/9 的 Collection 2 Level-2 可以直接用LANDSAT/LC08/C02/T1_L2和LC09地表反射率缩放系数是 0.0000275 加 -0.2和 Sentinel-2 不一样两块数据混用时必须各自缩放后再拼接否则同一类型的光谱值会系统性偏移。2.3 训练样本数量怎么算分层抽样的最小样本量拍脑袋定每类 50 个点是最常见的做法也是最容易在评审时被问住的。简单随机抽样下比例估计的样本量公式是n Z² · p · (1-p) / e²取置信水平 95%Z1.96、预期精度 p0.5最保守、允许误差 e5%得到 n≈384。这是总样本量的下限不是每类下限。每类实际需要的样本数取决于该类在区域内的面积占比和类内方差占比小的类型湿地、水体方差大即使面积很小也要保证足够的绝对点数。import numpy as np def stratified_samples(weights, p_expected, target_se0.01, z1.96): weights: 各层面积占比和为 1 p_expected: 各层预期用户精度如 0.85 target_se: 总体精度估计的目标标准误0.01 表示 1% 返回总样本量 n 和各层分配数 n_i W np.asarray(weights, dtypefloat) p np.asarray(p_expected, dtypefloat) S np.sqrt(p * (1 - p)) # 层内标准差 n (np.sum(W * S) / target_se) ** 2 alloc n * (W * S) / np.sum(W * S) # 按 W*S 比例分配 return int(np.ceil(n)), np.ceil(alloc).astype(int) W [0.42, 0.08, 0.13, 0.05, 0.20, 0.06, 0.04, 0.02] # 各类型面积占比 p [0.90, 0.80, 0.85, 0.75, 0.90, 0.88, 0.80, 0.95] n, alloc stratified_samples(W, p, target_se0.01) print(n, alloc)target_se0.01表示希望总体精度的标准误控制在 1% 量级对应 95% 置信区间约 ±2%这是一个偏严的设定实际项目里 0.015 到 0.02 也很常见。分配用的是 Olofsson 提出的按W_i × S_i比例分配比按面积比例分配更合理面积占比大但类内一致的层比如水体少分一些面积小但异质的层比如湿地多分一些。算出来的每层分配数还要过一遍地板任何一层的样本数低于 30都要手动抬到 30 以上否则那一类的精度估计没有统计意义。2.4 底图自检指数阈值与异常像元排查合成完先别急着分类跑一遍指数分布检查。水体的 NDVI 应该整体在 -0.1 以下如果出现大片水体像元 NDVI 大于 0.2大概率是云影残留或者水体边界混合像元没有剔除。建成区的 NDBI 应该在 0.1 以上如果低于 0要么是建筑密度太低要么是大气校正有问题。import rasterio, numpy as np with rasterio.open(eco_comp_2024.tif) as src: red src.read(3).astype(float32) nir src.read(7).astype(float32) swir src.read(8).astype(float32) ndvi (nir - red) / (nir red 1e-6) ndbi (swir - nir) / (swir nir 1e-6) # 统计每个指数的分布检查是否出现异常值域 for name, arr in [(NDVI, ndvi), (NDBI, ndbi)]: valid arr[np.isfinite(arr)] print(name, np.percentile(valid, [1, 5, 50, 95, 99]).round(3)) # 计数可疑像元NDVI 大于 0.2 但被标为水体的区域 susp np.isfinite(ndvi) (ndvi 0.2) print(可疑像元数, int(susp.sum()))百分位打印出来后重点看 1% 和 99% 分位是否落在合理值域内。NDVI 的 99% 分位如果超过 0.95通常是残云1% 分位如果低于 -0.9多半是阴影或者水体边界。这一步不解决分类问题但能在训练样本选取之前发现底图缺陷比分类完再返工便宜得多。3. 特征工程与分类器参数把容易混的类型分开3.1 光谱、指数、纹理、物候四类特征的组合方式单靠原始波段绝大多数项目的结果都止步于 75% 左右的总体精度。要往上走得按四类特征堆叠。光谱特征是 B2 到 B12 的地表反射率Sentinel-2 的红边波段 B5、B6、B7 对植被类型区分贡献明显尤其是区分针叶林和阔叶林因为红边位置会随叶绿素含量和叶结构变化。指数特征是 NDVI、NDWI、NDBI、SAVI、EVI 以及差值指数比如 NDVI 减去 NDWI 能强化湿地和陆生植被的边界。纹理特征用 GLCM窗口 3×3 或 5×5取均值、方差、熵、同质性四个量。灌丛和幼龄林在光谱上接近但灌丛的纹理熵明显更高、同质性更低。纹理窗口不能太大5×5 在 10 米影像上对应 50 米地面范围再大就会把地块边界的信息抹掉。物候特征是从时序影像集合里降维出来的常用的有各月 NDVI 中值、NDVI 年振幅、NDVI 最大值出现月份、生长季长度。农田和草地的分离几乎全靠这一类特征草地 NDVI 峰值出现在 7 月下旬、单峰华北平原冬小麦是双峰5 月和 9 月各一个峰。这套特征需要单独从原始影像集合按月合成再提取不能从中值合成结果里算。特征类型具体变量典型贡献场景注意事项光谱B2-B8、B11、B12、红边 B5-B7基础区分必须做反射率缩放指数NDVI、NDWI、NDBI、SAVI、EVI水陆边界、植被活力分母加 1e-6 防除零纹理GLCM 均值/方差/熵/同质性灌丛与幼林、建成区窗口 3×3 或 5×5物候月均 NDVI、年振幅、峰值月份农田与草地、落叶与常绿需从原始集合单独计算地形高程、坡度、坡向山地垂直带、阴影区坡向做 sin/cos 分解3.2 随机森林在 GEE 和本地 Python 里的两套实现随机森林在生态系统分类里几乎是默认选择原因是它能吃高维特征、不需要正态假设、对样本不平衡比最大似然宽容。GEE 里的实现// GEE随机森林分类 var bands [B2,B3,B4,B5,B6,B7,B8,B11,B12, NDVI,NDWI,NDBI,ELEV,SLOPE,GLCM_ENT]; var samples forest_grass_water_etc.merge(grass).merge(water); // 各类型样本合并 var classifier ee.Classifier.smileRandomForest({ numberOfTrees: 500, // 树的数量500 是精度与耗时的平衡点 variablesPerSplit: 4, // 每次分裂随机选的特征数默认 sqrt(特征数) bagFraction: 0.6, // 每棵树的样本抽样比例 minLeafPopulation: 5, // 叶节点最小样本防止过拟合 seed: 42 // 固定随机种子保证可复现 }).train({ features: samples, classProperty: class_id, inputProperties: bands }); var classified composite.select(bands).classify(classifier); // 输出特征重要性用于反向删减冗余特征 print(importance, classifier.explain().get(importance));numberOfTrees从 100 加到 500总体精度通常能提升 2 到 4 个百分点再加到 1000 收益就很小了而训练时间线性增长。bagFraction设 0.6 而不是默认的 0.5是在样本总量偏少时保留更多信息如果每类样本超过 500 个回到 0.5 更稳。seed必须固定否则每次跑出来的分类图会有细微差异评审时说不清楚。本地 Python 版本适合需要精细调参和交叉验证的场景from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import confusion_matrix, cohen_kappa_score import numpy as np X np.load(features.npy) # 形状 (n_samples, n_features) y np.load(labels.npy) X_tr, X_te, y_tr, y_te train_test_split( X, y, test_size0.3, stratifyy, random_state42) # 分层切分保证每类都有测试样本 rf RandomForestClassifier( n_estimators500, max_featuressqrt, min_samples_leaf3, class_weightbalanced_subsample, # 小类样本自动加权 oob_scoreTrue, n_jobs-1, random_state42) rf.fit(X_tr, y_tr) pred rf.predict(X_te) print(OOB, round(rf.oob_score_, 4)) print(Kappa, round(cohen_kappa_score(y_te, pred), 4)) print(重要性前五, np.argsort(rf.feature_importances_)[-5:])stratifyy在样本不平衡时很关键不加的话小类可能一个测试样本都分不到。class_weightbalanced_subsample让每棵树在自助采样时按类频率反向加权比全局balanced更适应分层结构。oob_score是袋外估计能免费得到一版泛化误差参考但它和独立测试集算出来的精度通常有 1 到 3 个百分点的差距正式报告里仍然要用独立核查数据。3.3 面向对象分割的尺度参数怎么调地块破碎、类型边界清晰的区域面向对象比逐像元效果好。GEE 里能直接用 SNIC 超像素分割// GEESNIC 超像素分割再用分割单元均值做分类 var seg ee.Algorithms.Image.Segmentation.SNIC({ image: composite.select([B2,B3,B4,B8,B11,B12]), size: 32, // 种子点间距像元对应地面 320 米 compactness: 0.3, // 形状紧致度越小边界越贴合 connectivity: 8, neighborhoodSize: 256 }); // seg 带出 segments 波段和每个波段的均值波段 var segMean seg.select([B2_mean,B3_mean,B4_mean,B8_mean,B11_mean,B12_mean]);size决定分割尺度10 米影像上常用 20 到 50对应地面 200 到 500 米。设小了图斑细碎、计算量大设大了会把不同地类合并进同一单元。compactness在 0.1 到 0.9 之间调值小的时候边界贴合实际地物轮廓但形状不规则值大的时候形状规整但会切掉边角。实操里的判断标准是分割边界和影像上的地块边界目视吻合度达到八成即可不必追求完美因为后续分类器会修正一部分分割误差。3.4 分类后处理众数滤波与最小图斑合并分类图出来一定有椒盐噪声尤其是逐像元分类。两步后处理先做 3×3 众数滤波再做连通域分析去除小于 MMU 的斑块。import rasterio, numpy as np from scipy import ndimage with rasterio.open(classified.tif) as src: arr src.read(1) profile src.profile # 众数滤波用 3x3 邻域内出现次数最多的类别替换中心像元 from scipy.ndimage import generic_filter mode_filter lambda v: np.bincount(v.astype(int)).argmax() smoothed generic_filter(arr, mode_filter, size3) # 去除小于 100 个像元约 1 公顷的斑块 for cls in np.unique(smoothed): if cls 0: continue mask smoothed cls labeled, num ndimage.label(mask) sizes ndimage.sum(mask, labeled, range(1, num 1)) for i, s in enumerate(sizes, start1): if s 100: # 用邻域众数填充小斑块 region labeled i dil ndimage.binary_dilation(region, iterations2) neigh smoothed[dil ~region] if neigh.size: smoothed[region] np.bincount(neigh.astype(int)).argmax() profile.update(dtypeuint8, count1) with rasterio.open(classified_clean.tif, w, **profile) as dst: dst.write(smoothed.astype(uint8), 1)众数滤波的窗口不要超过 3×3。5×5 会把窄条状的河渠、田埂、防护林带整体抹掉这些地物在生态评估里恰恰是重点。小斑块填充用邻域众数而不是直接归并到面积最大的邻类是因为小斑块周边可能有多个类型简单归并会系统性偏向大类。填充前记录一下被填充的斑块数量和总面积占比如果超过 5%说明分类器在某一类上本身就偏弱得回头查特征而不是靠后处理掩盖。4. 野外核查的抽样设计与现场作业4.1 分层随机抽样的核查点数量计算核查点不是多跑几个地方而是按统计设计抽的。分层随机抽样里每层分配数按W_i × S_i比例确定这套方法在 Olofsson 等人的精度评价框架里有完整推导。核心是把分类结果本身当作分层依据每层按面积占比和预期类内一致性分配核查点。import numpy as np def check_points(W, p_user, target_se0.012): W: 各层面积占比来自分类结果统计 p_user: 各层预期用户精度通常取 0.7~0.95 target_se: 目标标准误0.012 对应 95% 置信区间约 ±2.4% W np.asarray(W, dtypefloat) p np.asarray(p_user, dtypefloat) S np.sqrt(p * (1 - p)) n (np.sum(W * S) / target_se) ** 2 alloc np.ceil(n * (W * S) / np.sum(W * S)).astype(int) alloc np.maximum(alloc, 30) # 每层至少 30 个点 return int(np.ceil(n)), alloc W [0.42, 0.08, 0.13, 0.05, 0.20, 0.06, 0.04, 0.02] p [0.88, 0.75, 0.82, 0.70, 0.90, 0.85, 0.78, 0.95] n, alloc check_points(W, p) print(f总点数 {n}各层分配 {alloc})面积占比从分类结果里统计出来这一步很关键如果分类把灌丛高估了按这个权重分配就会在灌丛上多布点反而浪费野外工作量。稳健的做法是用上一轮核查后校正过的面积重新分配。alloc np.maximum(alloc, 30)强制每层至少 30 点是为了任何一类都能单独算出有意义的精度区间代价是总面积小的类型会占用偏多的野外时间但这是可以接受的取舍。4.2 核查点布设规则与不可达点的替换算出各层点数后在每层内部做简单随机抽样落点然后过三道筛子。第一道剔除距离类型边界小于 1 个像元的点这些点是混合像元现场核查结果和解译类型对不上是正常的不构成误差。第二道剔除落在水域中央、悬崖、军事管理区等不可达位置的点。第三道检查空间均衡性如果某个 10 公里格网内落了超过 5 个点按层内随机重抽。替换规则必须写进方案并留痕某个点确定不可达时在同层内、同一纬度带内、距原点最近的备选点中抽取一个同时在记录表里写明原编号、替换编号、替换原因。不做留痕的替换在评审时等于没有抽样设计因为无法判断替换是否引入了空间偏倚。实际布点时我会多抽 20% 的备选点和主点一起导出成 GPX 或者 GeoPackage现场直接用手机或者手持终端打开。这一步省下来的时间远比多抽的点值钱。4.3 现场表单字段与采集工具配置表单字段设计决定了这批数据能不能用于精度评价。字段分三组定位信息、类型判定信息、佐证信息。类型判定必须给出和分类体系一一对应的编码不能让调查员自由填写文字类型名否则回填时要做一遍人工映射错误率极高。字段名类型必填说明point_id文本是与抽样方案一致的唯一编号lon / lat数值是WGS84 十进制度保留 6 位gps_accuracy数值是终端报告的定位精度米eco_type_code单选是与解译体系一致的类型编码dominant_species文本否优势种最多填三种cover_pct数值是总盖度百分比canopy数值否郁闭度森林类型必填height_m数值否平均高度森林和灌丛必填photo_id文本是照片编号与文件命名对应match_flag单选是一致 / 不一致 / 无法判定note文本否备注不一致时必须说明原因采集工具用 KoboToolbox、ODK 或者 QField 都行共同要求是能离线填写、能带 GPS 定位、能拍照并在文件里记录坐标。XLSForm 的 survey 表结构大致是这样!-- XLSForm survey 表节选用于离线野外核查 -- survey text namepoint_id label核查点编号/ geopoint namelocation label定位/ select1 nameeco_type_code label生态系统类型 itemlabel森林/labelvalue1/value/item itemlabel灌丛/labelvalue2/value/item itemlabel草地/labelvalue3/value/item itemlabel湿地/labelvalue4/value/item itemlabel农田/labelvalue5/value/item /select1 decimal namecover_pct label总盖度(%)/ select1 namematch_flag label与解译是否一致 itemlabel一致/labelvalue1/value/item itemlabel不一致/labelvalue0/value/item itemlabel无法判定/labelvalue9/value/item /select1 text namenote label备注 relevant${match_flag}0/ /surveyrelevant${match_flag}0让备注字段只在不一致时出现能显著降低调查员乱填的比例。geopoint类型会自动记录经纬度、海拔和定位精度三个值比手动填经纬度可靠。表单最后要留一个隐藏字段记录填写时间戳用于事后检查数据是否集中在某一天批量补录。4.4 GPS 精度控制与照片命名规范手持终端在树冠下定位偏差可能超过 15 米而核查点要判定的对象可能只有 10 米宽。控制办法是到达点位后站定等待 30 秒以上让终端收敛到 5 米以内再记录如果始终收不敛沿着点位前后各走 10 米取三次定位取中值并在备注里写明。森林类型的郁闭度观测要离开点位 15 米以上避免自己站在样方中心挡住卫星视线。样方大小按类型定草地 1 米×1 米灌丛 5 米×5 米森林 20 米×20 米湿地在植被均匀的区域用 2 米×2 米。这个不一致会让统计变得麻烦解决办法是在表单里加一个样方面积字段回填时统一折算成百分比。照片命名统一成point_id_方位_序号.jpg比如P0132_E_01.jpg方位用 N/E/S/W。每个点至少拍四张东、南、西、北各一张再加一张带编号牌的特写。照片不按规范命名几百个点的照片在回填阶段就是一场灾难靠时间戳和 GPS 反查的效率极低。5. 精度评价与误差溯源5.1 混淆矩阵计算与精度指标解读核查数据回来之后把它和解译结果按point_id关联构成参考-预测配对再算混淆矩阵。import numpy as np from sklearn.metrics import confusion_matrix, cohen_kappa_score labels [1, 2, 3, 4, 5, 6, 7, 8] y_true np.load(ref_label.npy) # 野外核查判定的类型 y_pred np.load(map_label.npy) # 解译图在该点上的类型 cm confusion_matrix(y_true, y_pred, labelslabels) oa np.trace(cm) / cm.sum() kappa cohen_kappa_score(y_true, y_pred, labelslabels) # 生产者精度漏分和用户精度错分 pa np.diag(cm) / cm.sum(axis1) # 行归一化 ua np.diag(cm) / cm.sum(axis0) # 列归一化 f1 2 * pa * ua / (pa ua 1e-9) print(总体精度, round(oa, 4), Kappa, round(kappa, 4)) for i, lb in enumerate(labels): print(f类{lb}: PA{pa[i]:.3f} UA{ua[i]:.3f} F1{f1[i]:.3f})cm.sum(axis1)是按参考类求和得到生产者精度也就是该类在现实中真正被判对的比例反映漏分cm.sum(axis0)按预测类求和得到用户精度反映错分。两个指标经常被搞反判断方法是问自己这个分母是我在地面上看到的那一类还是我在图上画出来的那一类。Kappa 现在争议比较大很多期刊要求同时报总体精度和各类的置信区间Kappa 作为辅助指标列出即可不要只报 Kappa。5.2 面积估计的偏差校正分类图的面积统计直接拿来用是有偏的因为错分和漏分不对称。基于分层抽样和核查数据可以做无偏的面积估计Â_i A_total × Σ_j W_j × (n_ji / n_j)其中 W_j 是第 j 层的面积占比n_j 是第 j 层的核查点数n_ji 是第 j 层中被核查判定为第 i 类的点数。这个估计量的好处是只用核查数据的判定结果不依赖分类器在该层的正确率假设。算完之后把校正面积和原始分类面积并排列一张表两者差异超过 15% 的类型基本可以确定是分类器系统性偏差需要回去查特征或者样本。5.3 误差溯源清单与迭代修正拿到混淆矩阵后不要急着调参重跑。先按误差类型归类再针对性修。误差现象常见原因修正动作草地大量错分为农田单时相特征不足补入 NDVI 年振幅和峰值月份灌丛与幼林互错纹理特征窗口不当GLCM 窗口从 3×3 改到 5×5湿地漏分严重季节性淹水未捕捉增加春季和秋季各一期影像水体错分到建成区阴影和深色建筑干扰加入 NDWI 并提升云影掩膜阈值某一类 UA 高 PA 低样本偏向典型区域在类型边缘和过渡带补样本山地区域整体偏移地形阴影未校正加入坡度和坡向特征或做地形校正迭代的顺序是先修特征再补样本最后才动分类器参数。绝大多数精度问题出在特征和样本上调numberOfTrees这类参数能带来的提升通常不超过 2 个百分点。每轮迭代后重新算一遍混淆矩阵重点看上一轮问题类的 PA 和 UA 是否同时改善——只有一个改善说明误差只是在类之间转移没有真正消除。一个实用的技巧是给每轮迭代留一份误差点集把所有核查判定和解译不一致的点单独导出成一个点图层叠上影像和特征值人工看二十个点模式通常就很清楚了。这二十个点看下来得到的信息比跑十次网格搜索参数多得多。本文还有配套的精品资源点击获取
生态系统遥感解译与野外核查：从分类体系到精度评价闭环
简介这份《生态系统遥感解译与野外核查》资料面向从事生态监测、遥感解译与国土生态评估的技术人员及高校相关专业师生帮助解决生态状况调查中解译标准不统一、野外核查流程不规范的问题。内容围绕全国生态状况调查评估技术规范的编制说明展开涵盖项目背景与主要工作过程、标准制修订必要性分析、国内外相关标准对比、基本原则与技术路线以及适用范围、规范性引用文件、术语定义、遥感解译、野外核查与生态系统类型解译结果等主要技术内容并延伸到同类标准水平对比、实施管理措施与技术措施建议、环境效益及经济技术分析可帮助读者理解从数据获取、图像处理、解译到实地验证的完整链条。资源为1个PDF文件压缩包约265KB轻量便于随时查阅与打印。目前已有222人学习下载适合需要把握生态遥感监测标准化思路、撰写技术规范或开展野外核查工作的读者参考。1. 生态系统遥感解译与野外核查为什么必须成对做很多团队交出的第一版成果是一张生态系统类型分布图加一份精度评价表反馈回来的问题却是核查点的原始记录在哪、抽样方案怎么定的、为什么这几个类型一个样点都没有。生态系统遥感解译和野外核查在项目管理里常被拆成两件事一个靠算法一个靠出差一旦拆开做精度数字就成了没有证据链的孤证。解译是把影像的光谱、物候、纹理映射到一套类型体系上的推断过程核查则是用可重复的抽样和现场记录去约束这个推断的误差范围。前者决定你能分到多细后者决定你敢报多大的精度。做生态状况评估、自然保护地监测、碳汇本底调查、生态修复成效评估的工程师基本都会踩到这条线刚接手一份只有几百页工作方案 PDF、没写过一行分类代码的人也一样。往下走真正卡人的不是算法新旧而是分类体系、训练样本、抽样方案、验证口径这四件事能不能对齐。2. 解译底图与分类体系把口径定死再动算法2.1 生态系统类型体系怎么拆到可解译的最小单元分类体系不能直接从生态学教科书上抄。教科书按植被型、群系、群丛往下分一个群丛可能只有几公顷10 米分辨率的影像根本分不出来。落到工程里先把类型压到一级类森林、灌丛、草地、湿地、农田、建成区、裸地/稀疏植被、水体这八类在多数区域能做到 80% 以上的可分离度。二级类要按能不能被特征区分来定而不是按科学上该不该分来定。比如常绿针叶林和落叶阔叶林在生长季峰值期的单时相影像上几乎不可分但加上冬季影像或者 NDVI 年振幅就能拉开草地和农田在 7 月的华北平原光谱高度重叠得靠地块边界规整度和收割前后的物候突变来分。最小上图单元MMU必须显式写进方案。10 米影像常用 0.1 公顷换算成 10×10 像元20 米影像通常放到 0.25 公顷。不写 MMU后处理去小斑块时就没有依据做出来的图斑数量会是实际管理单元的几倍。混合像元也要提前约定地类边界内外各 1 个像元不参与训练样本选取这一点在山区和破碎地块尤其要紧。一级类主要区分特征10m 影像建议 MMU单类训练样本下限森林NDVI 高、红边斜率大、纹理粗糙0.1 hm²80灌丛NDVI 中高、纹理细碎、物候振幅小0.1 hm²60草地生长季 NDVI 峰值低、年振幅大0.1 hm²60湿地NDWI 高、季节性淹水、植被混杂0.1 hm²50农田地块规整、收割期突变、多峰物候0.1 hm²80建成区NDBI 高、纹理规则、时序稳定0.1 hm²50裸地亮度高、NDVI 低且平稳0.25 hm²50水体NDWI 高、NDVI 负值0.1 hm²402.2 Sentinel-2 与 Landsat 的时相选择与云掩膜处理时相选择比影像源选择更影响结果。北方地区类型区分主要靠 6 到 8 月的生长季峰值合成落叶与常绿、农田与草地的区分则要补一期 11 月到次年 2 月的冬季影像。南方多云雨地区一个月的无云影像可能只有两三景得把窗口放到整个生长季做中值合成。下面这段在 GEE 里完成 Sentinel-2 L2A 的生长季去云合成用的是 SCL 分类层加云概率层双重掩膜。// GEESentinel-2 L2A 生长季无云中值合成 var roi ee.Geometry.Rectangle([110.2, 30.1, 111.8, 31.5]); var s2 ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED) .filterDate(2024-06-01, 2024-09-30) .filterBounds(roi) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 40)); var cld ee.ImageCollection(COPERNICUS/S2_CLOUD_PROBABILITY) .filterDate(2024-06-01, 2024-09-30) .filterBounds(roi); // 按 system:index 把云概率影像挂到 SR 影像上 var joined ee.Join.saveFirst(cloud).apply({ primary: s2, secondary: cld, condition: ee.Filter.equals({leftField: system:index, rightField: system:index}) }); // 双重掩膜云概率 40%同时剔除 SCL 中的云、云影、卷云 var maskCloud function(img) { var prob ee.Image(img.get(cloud)).select(probability); var scl img.select(SCL); var sclMask scl.neq(3) // 云影 .and(scl.neq(8)) // 中概率云 .and(scl.neq(9)) // 高概率云 .and(scl.neq(10)); // 卷云 return img.updateMask(prob.lt(40)).updateMask(sclMask); }; var composite ee.ImageCollection(joined).map(maskCloud) .select([B2,B3,B4,B5,B6,B7,B8,B11,B12]) .median() // 中值合成抗残留云 .multiply(0.0001) // 反射率缩放 .clip(roi); Export.image.toDrive({image: composite, description: eco_comp_2024, region: roi, scale: 10, maxPixels: 1e13});filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 40))这一条把整景云量过高的影像先筛掉阈值放到 40% 而不是 10%是为了在南方保住足够的影像数量剩下的像元级云由掩膜处理。ee.Join.saveFirst把云概率集合按system:index对齐到 SR 影像两套集合的 index 是一致的这个 join 只是把两个集合合并成一条记录不改变像元。SCL 层里剔掉的四类对应云影、中概率云、高概率云、卷云水体6和雪11没有剔如果研究区有季节性积雪需要额外加条件。中值合成对残留薄云比均值合成稳健得多代价是物候信息被压平所以时序特征要单独从原始影像集合里算不能从这个合成结果里取。Landsat 8/9 的 Collection 2 Level-2 可以直接用LANDSAT/LC08/C02/T1_L2和LC09地表反射率缩放系数是 0.0000275 加 -0.2和 Sentinel-2 不一样两块数据混用时必须各自缩放后再拼接否则同一类型的光谱值会系统性偏移。2.3 训练样本数量怎么算分层抽样的最小样本量拍脑袋定每类 50 个点是最常见的做法也是最容易在评审时被问住的。简单随机抽样下比例估计的样本量公式是n Z² · p · (1-p) / e²取置信水平 95%Z1.96、预期精度 p0.5最保守、允许误差 e5%得到 n≈384。这是总样本量的下限不是每类下限。每类实际需要的样本数取决于该类在区域内的面积占比和类内方差占比小的类型湿地、水体方差大即使面积很小也要保证足够的绝对点数。import numpy as np def stratified_samples(weights, p_expected, target_se0.01, z1.96): weights: 各层面积占比和为 1 p_expected: 各层预期用户精度如 0.85 target_se: 总体精度估计的目标标准误0.01 表示 1% 返回总样本量 n 和各层分配数 n_i W np.asarray(weights, dtypefloat) p np.asarray(p_expected, dtypefloat) S np.sqrt(p * (1 - p)) # 层内标准差 n (np.sum(W * S) / target_se) ** 2 alloc n * (W * S) / np.sum(W * S) # 按 W*S 比例分配 return int(np.ceil(n)), np.ceil(alloc).astype(int) W [0.42, 0.08, 0.13, 0.05, 0.20, 0.06, 0.04, 0.02] # 各类型面积占比 p [0.90, 0.80, 0.85, 0.75, 0.90, 0.88, 0.80, 0.95] n, alloc stratified_samples(W, p, target_se0.01) print(n, alloc)target_se0.01表示希望总体精度的标准误控制在 1% 量级对应 95% 置信区间约 ±2%这是一个偏严的设定实际项目里 0.015 到 0.02 也很常见。分配用的是 Olofsson 提出的按W_i × S_i比例分配比按面积比例分配更合理面积占比大但类内一致的层比如水体少分一些面积小但异质的层比如湿地多分一些。算出来的每层分配数还要过一遍地板任何一层的样本数低于 30都要手动抬到 30 以上否则那一类的精度估计没有统计意义。2.4 底图自检指数阈值与异常像元排查合成完先别急着分类跑一遍指数分布检查。水体的 NDVI 应该整体在 -0.1 以下如果出现大片水体像元 NDVI 大于 0.2大概率是云影残留或者水体边界混合像元没有剔除。建成区的 NDBI 应该在 0.1 以上如果低于 0要么是建筑密度太低要么是大气校正有问题。import rasterio, numpy as np with rasterio.open(eco_comp_2024.tif) as src: red src.read(3).astype(float32) nir src.read(7).astype(float32) swir src.read(8).astype(float32) ndvi (nir - red) / (nir red 1e-6) ndbi (swir - nir) / (swir nir 1e-6) # 统计每个指数的分布检查是否出现异常值域 for name, arr in [(NDVI, ndvi), (NDBI, ndbi)]: valid arr[np.isfinite(arr)] print(name, np.percentile(valid, [1, 5, 50, 95, 99]).round(3)) # 计数可疑像元NDVI 大于 0.2 但被标为水体的区域 susp np.isfinite(ndvi) (ndvi 0.2) print(可疑像元数, int(susp.sum()))百分位打印出来后重点看 1% 和 99% 分位是否落在合理值域内。NDVI 的 99% 分位如果超过 0.95通常是残云1% 分位如果低于 -0.9多半是阴影或者水体边界。这一步不解决分类问题但能在训练样本选取之前发现底图缺陷比分类完再返工便宜得多。3. 特征工程与分类器参数把容易混的类型分开3.1 光谱、指数、纹理、物候四类特征的组合方式单靠原始波段绝大多数项目的结果都止步于 75% 左右的总体精度。要往上走得按四类特征堆叠。光谱特征是 B2 到 B12 的地表反射率Sentinel-2 的红边波段 B5、B6、B7 对植被类型区分贡献明显尤其是区分针叶林和阔叶林因为红边位置会随叶绿素含量和叶结构变化。指数特征是 NDVI、NDWI、NDBI、SAVI、EVI 以及差值指数比如 NDVI 减去 NDWI 能强化湿地和陆生植被的边界。纹理特征用 GLCM窗口 3×3 或 5×5取均值、方差、熵、同质性四个量。灌丛和幼龄林在光谱上接近但灌丛的纹理熵明显更高、同质性更低。纹理窗口不能太大5×5 在 10 米影像上对应 50 米地面范围再大就会把地块边界的信息抹掉。物候特征是从时序影像集合里降维出来的常用的有各月 NDVI 中值、NDVI 年振幅、NDVI 最大值出现月份、生长季长度。农田和草地的分离几乎全靠这一类特征草地 NDVI 峰值出现在 7 月下旬、单峰华北平原冬小麦是双峰5 月和 9 月各一个峰。这套特征需要单独从原始影像集合按月合成再提取不能从中值合成结果里算。特征类型具体变量典型贡献场景注意事项光谱B2-B8、B11、B12、红边 B5-B7基础区分必须做反射率缩放指数NDVI、NDWI、NDBI、SAVI、EVI水陆边界、植被活力分母加 1e-6 防除零纹理GLCM 均值/方差/熵/同质性灌丛与幼林、建成区窗口 3×3 或 5×5物候月均 NDVI、年振幅、峰值月份农田与草地、落叶与常绿需从原始集合单独计算地形高程、坡度、坡向山地垂直带、阴影区坡向做 sin/cos 分解3.2 随机森林在 GEE 和本地 Python 里的两套实现随机森林在生态系统分类里几乎是默认选择原因是它能吃高维特征、不需要正态假设、对样本不平衡比最大似然宽容。GEE 里的实现// GEE随机森林分类 var bands [B2,B3,B4,B5,B6,B7,B8,B11,B12, NDVI,NDWI,NDBI,ELEV,SLOPE,GLCM_ENT]; var samples forest_grass_water_etc.merge(grass).merge(water); // 各类型样本合并 var classifier ee.Classifier.smileRandomForest({ numberOfTrees: 500, // 树的数量500 是精度与耗时的平衡点 variablesPerSplit: 4, // 每次分裂随机选的特征数默认 sqrt(特征数) bagFraction: 0.6, // 每棵树的样本抽样比例 minLeafPopulation: 5, // 叶节点最小样本防止过拟合 seed: 42 // 固定随机种子保证可复现 }).train({ features: samples, classProperty: class_id, inputProperties: bands }); var classified composite.select(bands).classify(classifier); // 输出特征重要性用于反向删减冗余特征 print(importance, classifier.explain().get(importance));numberOfTrees从 100 加到 500总体精度通常能提升 2 到 4 个百分点再加到 1000 收益就很小了而训练时间线性增长。bagFraction设 0.6 而不是默认的 0.5是在样本总量偏少时保留更多信息如果每类样本超过 500 个回到 0.5 更稳。seed必须固定否则每次跑出来的分类图会有细微差异评审时说不清楚。本地 Python 版本适合需要精细调参和交叉验证的场景from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import confusion_matrix, cohen_kappa_score import numpy as np X np.load(features.npy) # 形状 (n_samples, n_features) y np.load(labels.npy) X_tr, X_te, y_tr, y_te train_test_split( X, y, test_size0.3, stratifyy, random_state42) # 分层切分保证每类都有测试样本 rf RandomForestClassifier( n_estimators500, max_featuressqrt, min_samples_leaf3, class_weightbalanced_subsample, # 小类样本自动加权 oob_scoreTrue, n_jobs-1, random_state42) rf.fit(X_tr, y_tr) pred rf.predict(X_te) print(OOB, round(rf.oob_score_, 4)) print(Kappa, round(cohen_kappa_score(y_te, pred), 4)) print(重要性前五, np.argsort(rf.feature_importances_)[-5:])stratifyy在样本不平衡时很关键不加的话小类可能一个测试样本都分不到。class_weightbalanced_subsample让每棵树在自助采样时按类频率反向加权比全局balanced更适应分层结构。oob_score是袋外估计能免费得到一版泛化误差参考但它和独立测试集算出来的精度通常有 1 到 3 个百分点的差距正式报告里仍然要用独立核查数据。3.3 面向对象分割的尺度参数怎么调地块破碎、类型边界清晰的区域面向对象比逐像元效果好。GEE 里能直接用 SNIC 超像素分割// GEESNIC 超像素分割再用分割单元均值做分类 var seg ee.Algorithms.Image.Segmentation.SNIC({ image: composite.select([B2,B3,B4,B8,B11,B12]), size: 32, // 种子点间距像元对应地面 320 米 compactness: 0.3, // 形状紧致度越小边界越贴合 connectivity: 8, neighborhoodSize: 256 }); // seg 带出 segments 波段和每个波段的均值波段 var segMean seg.select([B2_mean,B3_mean,B4_mean,B8_mean,B11_mean,B12_mean]);size决定分割尺度10 米影像上常用 20 到 50对应地面 200 到 500 米。设小了图斑细碎、计算量大设大了会把不同地类合并进同一单元。compactness在 0.1 到 0.9 之间调值小的时候边界贴合实际地物轮廓但形状不规则值大的时候形状规整但会切掉边角。实操里的判断标准是分割边界和影像上的地块边界目视吻合度达到八成即可不必追求完美因为后续分类器会修正一部分分割误差。3.4 分类后处理众数滤波与最小图斑合并分类图出来一定有椒盐噪声尤其是逐像元分类。两步后处理先做 3×3 众数滤波再做连通域分析去除小于 MMU 的斑块。import rasterio, numpy as np from scipy import ndimage with rasterio.open(classified.tif) as src: arr src.read(1) profile src.profile # 众数滤波用 3x3 邻域内出现次数最多的类别替换中心像元 from scipy.ndimage import generic_filter mode_filter lambda v: np.bincount(v.astype(int)).argmax() smoothed generic_filter(arr, mode_filter, size3) # 去除小于 100 个像元约 1 公顷的斑块 for cls in np.unique(smoothed): if cls 0: continue mask smoothed cls labeled, num ndimage.label(mask) sizes ndimage.sum(mask, labeled, range(1, num 1)) for i, s in enumerate(sizes, start1): if s 100: # 用邻域众数填充小斑块 region labeled i dil ndimage.binary_dilation(region, iterations2) neigh smoothed[dil ~region] if neigh.size: smoothed[region] np.bincount(neigh.astype(int)).argmax() profile.update(dtypeuint8, count1) with rasterio.open(classified_clean.tif, w, **profile) as dst: dst.write(smoothed.astype(uint8), 1)众数滤波的窗口不要超过 3×3。5×5 会把窄条状的河渠、田埂、防护林带整体抹掉这些地物在生态评估里恰恰是重点。小斑块填充用邻域众数而不是直接归并到面积最大的邻类是因为小斑块周边可能有多个类型简单归并会系统性偏向大类。填充前记录一下被填充的斑块数量和总面积占比如果超过 5%说明分类器在某一类上本身就偏弱得回头查特征而不是靠后处理掩盖。4. 野外核查的抽样设计与现场作业4.1 分层随机抽样的核查点数量计算核查点不是多跑几个地方而是按统计设计抽的。分层随机抽样里每层分配数按W_i × S_i比例确定这套方法在 Olofsson 等人的精度评价框架里有完整推导。核心是把分类结果本身当作分层依据每层按面积占比和预期类内一致性分配核查点。import numpy as np def check_points(W, p_user, target_se0.012): W: 各层面积占比来自分类结果统计 p_user: 各层预期用户精度通常取 0.7~0.95 target_se: 目标标准误0.012 对应 95% 置信区间约 ±2.4% W np.asarray(W, dtypefloat) p np.asarray(p_user, dtypefloat) S np.sqrt(p * (1 - p)) n (np.sum(W * S) / target_se) ** 2 alloc np.ceil(n * (W * S) / np.sum(W * S)).astype(int) alloc np.maximum(alloc, 30) # 每层至少 30 个点 return int(np.ceil(n)), alloc W [0.42, 0.08, 0.13, 0.05, 0.20, 0.06, 0.04, 0.02] p [0.88, 0.75, 0.82, 0.70, 0.90, 0.85, 0.78, 0.95] n, alloc check_points(W, p) print(f总点数 {n}各层分配 {alloc})面积占比从分类结果里统计出来这一步很关键如果分类把灌丛高估了按这个权重分配就会在灌丛上多布点反而浪费野外工作量。稳健的做法是用上一轮核查后校正过的面积重新分配。alloc np.maximum(alloc, 30)强制每层至少 30 点是为了任何一类都能单独算出有意义的精度区间代价是总面积小的类型会占用偏多的野外时间但这是可以接受的取舍。4.2 核查点布设规则与不可达点的替换算出各层点数后在每层内部做简单随机抽样落点然后过三道筛子。第一道剔除距离类型边界小于 1 个像元的点这些点是混合像元现场核查结果和解译类型对不上是正常的不构成误差。第二道剔除落在水域中央、悬崖、军事管理区等不可达位置的点。第三道检查空间均衡性如果某个 10 公里格网内落了超过 5 个点按层内随机重抽。替换规则必须写进方案并留痕某个点确定不可达时在同层内、同一纬度带内、距原点最近的备选点中抽取一个同时在记录表里写明原编号、替换编号、替换原因。不做留痕的替换在评审时等于没有抽样设计因为无法判断替换是否引入了空间偏倚。实际布点时我会多抽 20% 的备选点和主点一起导出成 GPX 或者 GeoPackage现场直接用手机或者手持终端打开。这一步省下来的时间远比多抽的点值钱。4.3 现场表单字段与采集工具配置表单字段设计决定了这批数据能不能用于精度评价。字段分三组定位信息、类型判定信息、佐证信息。类型判定必须给出和分类体系一一对应的编码不能让调查员自由填写文字类型名否则回填时要做一遍人工映射错误率极高。字段名类型必填说明point_id文本是与抽样方案一致的唯一编号lon / lat数值是WGS84 十进制度保留 6 位gps_accuracy数值是终端报告的定位精度米eco_type_code单选是与解译体系一致的类型编码dominant_species文本否优势种最多填三种cover_pct数值是总盖度百分比canopy数值否郁闭度森林类型必填height_m数值否平均高度森林和灌丛必填photo_id文本是照片编号与文件命名对应match_flag单选是一致 / 不一致 / 无法判定note文本否备注不一致时必须说明原因采集工具用 KoboToolbox、ODK 或者 QField 都行共同要求是能离线填写、能带 GPS 定位、能拍照并在文件里记录坐标。XLSForm 的 survey 表结构大致是这样!-- XLSForm survey 表节选用于离线野外核查 -- survey text namepoint_id label核查点编号/ geopoint namelocation label定位/ select1 nameeco_type_code label生态系统类型 itemlabel森林/labelvalue1/value/item itemlabel灌丛/labelvalue2/value/item itemlabel草地/labelvalue3/value/item itemlabel湿地/labelvalue4/value/item itemlabel农田/labelvalue5/value/item /select1 decimal namecover_pct label总盖度(%)/ select1 namematch_flag label与解译是否一致 itemlabel一致/labelvalue1/value/item itemlabel不一致/labelvalue0/value/item itemlabel无法判定/labelvalue9/value/item /select1 text namenote label备注 relevant${match_flag}0/ /surveyrelevant${match_flag}0让备注字段只在不一致时出现能显著降低调查员乱填的比例。geopoint类型会自动记录经纬度、海拔和定位精度三个值比手动填经纬度可靠。表单最后要留一个隐藏字段记录填写时间戳用于事后检查数据是否集中在某一天批量补录。4.4 GPS 精度控制与照片命名规范手持终端在树冠下定位偏差可能超过 15 米而核查点要判定的对象可能只有 10 米宽。控制办法是到达点位后站定等待 30 秒以上让终端收敛到 5 米以内再记录如果始终收不敛沿着点位前后各走 10 米取三次定位取中值并在备注里写明。森林类型的郁闭度观测要离开点位 15 米以上避免自己站在样方中心挡住卫星视线。样方大小按类型定草地 1 米×1 米灌丛 5 米×5 米森林 20 米×20 米湿地在植被均匀的区域用 2 米×2 米。这个不一致会让统计变得麻烦解决办法是在表单里加一个样方面积字段回填时统一折算成百分比。照片命名统一成point_id_方位_序号.jpg比如P0132_E_01.jpg方位用 N/E/S/W。每个点至少拍四张东、南、西、北各一张再加一张带编号牌的特写。照片不按规范命名几百个点的照片在回填阶段就是一场灾难靠时间戳和 GPS 反查的效率极低。5. 精度评价与误差溯源5.1 混淆矩阵计算与精度指标解读核查数据回来之后把它和解译结果按point_id关联构成参考-预测配对再算混淆矩阵。import numpy as np from sklearn.metrics import confusion_matrix, cohen_kappa_score labels [1, 2, 3, 4, 5, 6, 7, 8] y_true np.load(ref_label.npy) # 野外核查判定的类型 y_pred np.load(map_label.npy) # 解译图在该点上的类型 cm confusion_matrix(y_true, y_pred, labelslabels) oa np.trace(cm) / cm.sum() kappa cohen_kappa_score(y_true, y_pred, labelslabels) # 生产者精度漏分和用户精度错分 pa np.diag(cm) / cm.sum(axis1) # 行归一化 ua np.diag(cm) / cm.sum(axis0) # 列归一化 f1 2 * pa * ua / (pa ua 1e-9) print(总体精度, round(oa, 4), Kappa, round(kappa, 4)) for i, lb in enumerate(labels): print(f类{lb}: PA{pa[i]:.3f} UA{ua[i]:.3f} F1{f1[i]:.3f})cm.sum(axis1)是按参考类求和得到生产者精度也就是该类在现实中真正被判对的比例反映漏分cm.sum(axis0)按预测类求和得到用户精度反映错分。两个指标经常被搞反判断方法是问自己这个分母是我在地面上看到的那一类还是我在图上画出来的那一类。Kappa 现在争议比较大很多期刊要求同时报总体精度和各类的置信区间Kappa 作为辅助指标列出即可不要只报 Kappa。5.2 面积估计的偏差校正分类图的面积统计直接拿来用是有偏的因为错分和漏分不对称。基于分层抽样和核查数据可以做无偏的面积估计Â_i A_total × Σ_j W_j × (n_ji / n_j)其中 W_j 是第 j 层的面积占比n_j 是第 j 层的核查点数n_ji 是第 j 层中被核查判定为第 i 类的点数。这个估计量的好处是只用核查数据的判定结果不依赖分类器在该层的正确率假设。算完之后把校正面积和原始分类面积并排列一张表两者差异超过 15% 的类型基本可以确定是分类器系统性偏差需要回去查特征或者样本。5.3 误差溯源清单与迭代修正拿到混淆矩阵后不要急着调参重跑。先按误差类型归类再针对性修。误差现象常见原因修正动作草地大量错分为农田单时相特征不足补入 NDVI 年振幅和峰值月份灌丛与幼林互错纹理特征窗口不当GLCM 窗口从 3×3 改到 5×5湿地漏分严重季节性淹水未捕捉增加春季和秋季各一期影像水体错分到建成区阴影和深色建筑干扰加入 NDWI 并提升云影掩膜阈值某一类 UA 高 PA 低样本偏向典型区域在类型边缘和过渡带补样本山地区域整体偏移地形阴影未校正加入坡度和坡向特征或做地形校正迭代的顺序是先修特征再补样本最后才动分类器参数。绝大多数精度问题出在特征和样本上调numberOfTrees这类参数能带来的提升通常不超过 2 个百分点。每轮迭代后重新算一遍混淆矩阵重点看上一轮问题类的 PA 和 UA 是否同时改善——只有一个改善说明误差只是在类之间转移没有真正消除。一个实用的技巧是给每轮迭代留一份误差点集把所有核查判定和解译不一致的点单独导出成一个点图层叠上影像和特征值人工看二十个点模式通常就很清楚了。这二十个点看下来得到的信息比跑十次网格搜索参数多得多。本文还有配套的精品资源点击获取