ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

西双版纳橡胶成林数据集:从遥感分类到变化检测实践

西双版纳橡胶成林数据集:从遥感分类到变化检测实践 简介西双版纳橡胶成林分布及变化数据集2014-2018是一份面向农林遥感、地理信息与区域可持续发展研究者的基础地理数据包聚焦中国最大天然橡胶产区五年间的橡胶林空间分布与动态变化可服务于农业政策评估、生态环境影响分析及教学科研。压缩包共39个文件以GeoTIFF栅格、Shapefile矢量、Excel统计表及辅助投影/元数据文件为主其中tif与shp分别刻画逐年橡胶林分布和边界范围xlsx提供种植面积、产量等关键统计指标。资源整体体积仅3.37MB轻量易用适合在ArcGIS、QGIS等平台快速加载分析。目前已有219人浏览学习。通过对比2014年与2018年数据可直接获得橡胶林扩展区、缩减区及不变区栅格图层结合研究区域边界和基础统计表格能够复现橡胶林动态变化的制图与量化分析流程为理解橡胶种植与热带生态环境之间的关联提供扎实的数据支撑。1. 从一份.rar看西双版纳橡胶林这五年发生了什么橡胶林是西双版纳最典型的土地利用类型之一它的分布直接关系到当地生态安全、水土保持和边疆经济。拿到“西双版纳橡胶成林分布及变化数据集2014-2018.rar”这份数据第一反应不该是急着解压看图而是先想清楚这里面装的是什么时相、什么分辨率、什么分类体系的栅格或矢量以及它能否支撑我后续的面积统计、变化检测和空间分析。这类数据集一般是基于 Landsat 8 OLI 影像结合物候特征分类得到的年际橡胶成林分布图2014 到 2018 五年正好覆盖了橡胶树从落叶期到展叶期的完整物候循环也跨过了 2015 年前后胶价波动对种植面积的影响窗口。对遥感、GIS 和农林从业者来说这份数据的价值不在于“有一张图”而在于能从中提取出可验证的时空变化证据对做生态模型的人来说它又是土地利用输入数据的可替代来源。本文就从解压预处理开始把它背后的分类逻辑、变化检测方法和精度验证口径逐个拆开。2. 解包与预处理先把数据从压缩包变成可分析的图层2.1 压缩包内部结构的典型形态与编码问题这类数据集的.rar内部通常会按年份组织文件常见形态是rubber_2014.tif、rubber_2015.tif这类单波段栅格取值一般为 0/1 或 1/2其中 1 表示橡胶成林0或 2表示非橡胶也有的数据集会附带一个Change_2014_2018.tif直接标记了“稳定橡胶”“新增橡胶”“流失橡胶”和“非橡胶”。少数情况下里面的主力成果是 Shapefile 格式的矢量多边形如rubber_2018.shp。矢量格式更便于做产权地块级别的统计但 2014 到 2018 的五期矢量如果拓扑不一致叠加分析时会出现碎面所以栅格反而更适合做严格的时间序列对比。解压本身不是难点坑在文件名编码。如果这份压缩包是在 Windows 环境下用中文文件名打包的在 Linux 上解压常会出现乱码# 先检查压缩包里的文件名编码 unrar l 西双版纳橡胶成林分布及变化数据集2014-2018.rar # 如果文件名乱码优先用 unar 按 GBK 解码 unar -e gbk 西双版纳橡胶成林分布及变化数据集2014-2018.rar # 或者用 7z 解压后批量改回可读文件名 7z x 西双版纳橡胶成林分布及变化数据集2014-2018.rarunar -e gbk这个参数的作用是强制把压缩包内的 GBK 编码文件名转换成 UTF-8避免后期在 Python 里读文件路径时因为编码不一致报No such file or directory。实际场景里我还会先unrar l列一遍文件清单对照查看是否存在 README、精度验证表如accuracy.csv或者 ROI 样本点文件这些往往比主分类图更值得优先阅读因为它们决定了你后续做变化检测时哪些图斑是可信的。2.2 坐标参考与栅格属性的统一五期影像分别打开没问题但一旦要做叠加分析坐标系不统一就是最大的隐患。西双版纳地区的遥感数据常见有两种坐标系WGS84 经纬度EPSG:4326和 UTM 47N 投影WGS 84 / UTM zone 47NEPSG:32647。前者适合做地图展示和经纬度采样后者适合做面积统计因为它的单位为米直接算像素个数乘分辨率就是面积。用 GDAL 检查所有期次是否在同一坐标系下for f in rubber_201*.tif; do echo $f gdalinfo $f | grep -E Origin|Pixel Size|EPSG|Data type done如果输出显示各期数据来自不同坐标系需要统一投影后再进行变化分析# 统一到 UTM 47N保留最近邻重采样以维持类别值不变 gdalwarp -t_srs EPSG:32647 -r near -of GTiff rubber_2014_wgs.tif rubber_2014_utm.tif这里有一个容易忽略的细节分类栅格的值是类别标签不是连续数值重采样时只能选nearest最近邻不能选bilinear或cubic否则会产生 0.5、0.67 这类不存在的类别值后续统计直接报错。另一点是检查每个栅格的NoData值是否统一。有的数据集把研究区外设为 0有的设为 255也有的不设置会直接影响你统计“橡胶面积”的结果。2.3 波段取值与属性类别的语义核对在动手做变化分析之前先打开属性表或者用 Python 数一遍各类别的像素量级确认分类体系是否符合预期这比画图更实在import rasterio import numpy as np for year in [2014, 2015, 2016, 2017, 2018]: with rasterio.open(frubber_{year}_utm.tif) as src: arr src.read(1) vals, counts np.unique(arr, return_countsTrue) px_area abs(src.transform.a * src.transform.e) # 单像素面积平方米 print(f{year}: , end) for v, c in zip(vals, counts): if v ! src.nodata: print(f类别{v} 面积{c * px_area / 1e6:.2f} km² | , end) print()这一步的作用是快速核验数据质量。如果 2014 年橡胶面积只有几十平方公里而 2016 年突然暴涨到几千平方公里大概率不是真实的扩张而是分类体系变了比如某一年把橡胶幼林也纳入了“成林”此时必须回到元数据或附带文档去核对分类阈值盲目拿着这个结果去写报告误差会很大。我在实际处理时还会顺手输出每期影像的覆盖范围看五期之间是否存在有效区域的偏移避免后续把“有效覆盖变化”误判成“橡胶面积变化”。3. 橡胶成林分类的技术路径这类分布图是怎么生产出来的3.1 为什么选择物候特征而非单时相光谱橡胶树和热带雨林在生长旺季的光谱特征差异很小单拿一张 12 月或 5 月的影像去做监督分类橡胶树很容易和天然林混在一起。但橡胶树有一个非常关键的物候窗口——每年 1 月到 2 月橡胶树会落叶叶片落光后树冠近乎光秃此时 NDVI 显著降低3 月到 4 月抽芽展叶NDVI 快速回升。天然常绿林没有这样强烈的冬季衰退过程所以在时间序列上橡胶林的 NDVI 曲线呈现一个明显的“V”形谷底这个特征在橡胶成林树龄 7 年以上、冠层封闭度足够上表现得尤为稳定。这份数据集的生产常见做法就是基于这个物候节奏构建时间序列特征。具体来说对一年内所有无云像元提取以下指标生长季 NDVI 最大值反映冠层茂盛程度落叶期 NDVI 最小值反映落叶强度落叶期与生长季的 NDVI 差值短波红外波段SWIR1在落叶期的反射率展叶期的物候速率NDVI 上升斜率这些特征输入随机森林分类器配合目视解译的样本点做训练可以在研究区内获得较高的分类精度。用这种方式2014 到 2018 年每年单独分类、每年独立选样本比直接把五期影像叠在一起做一次性分类更能保留年际间的独立变化信号但代价是需要为每一年都做一轮样本标注和精度评估。3.2 训练样本、分类器与关键参数样本设计是整个分类流程里最影响精度的环节。常见错误是只在路网两侧和村寨附近选样本导致模型对偏远区域的橡胶林缺乏泛化能力。更稳妥的做法是分层随机抽样先把当年影像做一次非监督聚类比如 ISODATA 分出 20 类在每个聚类簇内随机布设样本再结合高分辨率影像和人工经验判定类别这样样本能覆盖多样的光谱场景。随机森林分类器里有 3 个参数最值得调参数常见取值范围作用与风险n_estimators200 - 500太小模型不稳定太大对精度提升有限且耗时翻倍max_depth10 - 20过深会过拟合样本橡胶林和天然林的边界会被噪声打碎max_features特征数的平方根防止高相关特征主导分裂保持模型多样性from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split # X 为样本特征矩阵每行一个样本列含 NDVI/SWIR 等特征y 为类别标签 X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.3, stratifyy, random_state42 ) model RandomForestClassifier( n_estimators300, max_depth15, max_featuressqrt, n_jobs-1, random_state42 ) model.fit(X_train, y_train) print(OOB score:, model.oob_score_) print(Feature importance:, model.feature_importances_)stratifyy是保证训练集和验证集中橡胶/非橡胶的比例与总体一致的关键参数否则当橡胶面积占比只有 5% 时随机切分很容易让验证集里一个橡胶样本都没有。oob_score_是随机森林自带的袋外精度估计适合在还没做独立验证集时快速评估当前特征组合有没有救。3.3 分类后处理从像元噪声到图斑连贯分类结果直接使用时会出现两种典型噪声一是孤立像元二是椒盐效应。橡胶成林在空间上通常是连片分布的——这和幼林补植、胶园间作有本质区别因此可以用多数滤波Majority Filter去除孤立小斑块。但要提醒的是滤波窗口设太大会抹掉真实的零散胶林比如山地陡坡上的小块胶园这部分胶林虽然面积小却是生态争议的焦点。from scipy.ndimage import median_filter import rasterio with rasterio.open(rubber_2018_utm.tif) as src: arr src.read(1) profile src.profile # 3x3 中值滤波只作用于分类值nodata 保留 filtered median_filter(arr, size3, modenearest) with rasterio.open(rubber_2018_cleaned.tif, w, **profile) as dst: dst.write(filtered, 1)modenearest的作用是让研究区边缘的像元在滤波时用最邻近的有效值填充避免边缘出现一圈 0 值伪变化。这一步对后续计算面积很重要——如果边缘噪声没清理干净面积统计里会混入大量碎斑。做完滤波后我还习惯做一步连通域分析统计每个橡胶图斑的像元个数再人工查验一下面积小于 0.5 公顷约 6 个 30m 像元的小斑块在高分影像里是否真实存在。数据集如果标注了“成林”小斑块往往指向误分类而非真实胶园。4. 变化检测与时空分析从两期图层中提取可信的变化信号4.1 像元级叠加与转移矩阵拿到五期成林分布图之后最直接的变化分析方法就是逐像元叠加生成转移矩阵。具体做法是把 2014 年的类别作为行、2018 年的类别作为列统计每个像元从初始类别到最终类别的流向。以两期为例import numpy as np import rasterio with rasterio.open(rubber_2014_utm.tif) as s14, \ rasterio.open(rubber_2018_utm.tif) as s18: a14 s14.read(1).astype(np.int16) a18 s18.read(1) # 有效像元掩膜剔除 nodata mask (a14 0) (a18 0) # 计算转移矩阵行2014列2018共2类 transition np.zeros((2, 2), dtypenp.int64) for i in range(2): for j in range(2): transition[i, j] np.sum((a14[mask] i) (a18[mask] j)) print(转移矩阵2014 - 2018) print(transition)转移矩阵的四个数值分别对应稳定非橡胶、橡胶新增、橡胶流失、稳定橡胶。用面积表示为类别流向含义备注非橡胶 → 橡胶新增成林可能由幼林成熟或新开垦两种原因造成橡胶 → 非橡胶橡胶流失可能因砍伐、灾害或改种需逐斑核实橡胶 → 橡胶稳定橡胶2014 已是成林2018 仍为成林非橡胶 → 非橡胶无变化不参与变化讨论这一步的计算量不大但要特别小心坐标系和像元对齐问题。如果两期栅格的像元边界有半个像元的偏移叠加后会出现大量沿着地块边缘的“带状伪变化”面积看起来动辄几十平方公里。遇到这种情况先检查 TIFF 文件的transform是否完全一致不一致时用gdalwarp -tr 30 30 -tap强制对齐网格。4.2 年际时间序列与变化轨迹的识别五年五期数据只看首尾会丢失中间过程。比如某个像元在 2014 年是橡胶、2015 年变成非橡胶、2016 年又变回橡胶这种“抖动变化”大概率不是真实变化而是分类噪声。而 2015 年到 2017 年连续三年都是橡胶才算稳定信号。处理这类问题常见做法是引入时间一致性判别规则只有变化后状态至少连续维持两年的像元才被认定为真实变化否则标记为“不确定变化”。这可以用一组滑动窗口的多数投票来实现import numpy as np # 假设 five_years 形状为 (5, rows, cols)值为 0 或 1 five_years np.stack([a2014, a2015, a2016, a2017, a2018], axis0) # 时间维做 3 年窗口中值滤波 from scipy.ndimage import median_filter smoothed median_filter(five_years, size(3, 1, 1), modenearest) # 重新生成稳定变化掩膜任何相邻两年的变化都视为候选变化 change_candidates np.abs(np.diff(smoothed, axis0)).sum(axis0)size(3, 1, 1)表示只在时间维上做 3 年窗口滤波不改动空间结构这样能消除单一年份的孤立误分类。做完这一步再看首尾变化得到的面积会更接近真实。转移矩阵在变化检测中的意义除了算面积还能给变化图谱Change Map提供类别语义让后续的地块级核实有明确的目标位置。4.3 地块级汇总与变化热点识别像元级统计解决“变了多少”但回答不了“变在哪”。把变化结果转成矢量图斑再按乡镇或自然村做空间统计是更贴近业务需求的输出。典型做法是import geopandas as gpd import rasterio from rasterio.features import shapes with rasterio.open(change_2014_2018.tif) as src: change src.read(1) transform src.transform # 变化图斑矢量化值为1新增橡胶2流失0无变化 results [] for geom, value in shapes(change, mask(change 0), transformtransform): results.append({geometry: geom, change_type: int(value)}) gdf gpd.GeoDataFrame(results, crsEPSG:32647) gdf[area_km2] gdf.geometry.area / 1e6 # 按 change_type 聚合 summary gdf.groupby(change_type)[area_km2].sum() print(summary)这里把变化类型直接映射成整数后续用乡镇行政边界做空间连接就能输出每个乡镇的新增/流失统计表。注意矢量化之前变化栅格最好经过一次小斑块过滤如小于 3 个像元的图斑直接赋值为 0否则矢量文件会膨胀到几十万个碎面空间连接和渲染会明显卡顿。5. 精度验证与同数据源之间的交叉检验5.1 验证样本设计与混淆矩阵任何分类数据集都需要回答一个重要问题图中那些“橡胶”在图上的位置与真实地面情况有多吻合。2014 到 2018 年间没有同步的地面调查点做支撑时一个可行的替代方案是用更高分辨率的影像如 Google Earth 历史影像做分层随机抽样目视判读。每期抽 300 到 500 个样本点分层依据是类别面积比例橡胶类多抽非橡胶类按面积比例抽。验证完成后整理成验证点表格计算混淆矩阵验证点总数参考橡胶参考非橡胶用户精度分类橡胶1521889.4%分类非橡胶2320790.0%生产者精度86.9%92.0%总体精度 89.5%总体精度达到 85% 以上这个数据集用于区域尺度分析基本可靠但如果要精确到村镇尺度的面积统计用户精度低于 90% 时就要谨慎了。我判断一个数据集能不能用的标准是总体精度只是底线重点看橡胶类有没有系统性的漏分或错分因为任何误判在时间序列上会被放大成伪变化。5.2 历史高分影像复检与差异区域清单输出验证不做全图重点抽变化区域。方法是在变化矢量图斑内按面积大小分层采样——大图斑10 公顷全数检查中图斑抽取 30%小图斑抽取 10%。Google Earth 历史影像没有 OpenAPI 可以直接下载多数情况下是人工目视判读但流程可以半自动将图斑中心点导出为 KML按2014 年前/2018 年后两个时间点分别加载历史影像逐点检查。最后形成一份差异清单列出“数据集标记为新增但影像上看起来不像橡胶”的图斑编号、经纬度和面积。这份清单比一个笼统的精度表更有用——它能反过来帮你判断是数据集本身的质量问题还是 2018 年后橡胶树被砍伐了。5.3 与公开全球/区域数据产品的对比检验如果手里有同时段其他数据源做一次空间一致性对比是最快的检验手段。常见选择数据源分辨率可对比内容全球地表覆盖数据如 GlobeLand302010 与 2020 两期30m人工地表大类比对全球森林变化数据GFCHansen 等30m橡胶林区域树木覆盖度损失/增加趋势同区域其他学者的年度分类结果30m逐期面积曲线与空间重叠率比对方法不复杂把另一份数据裁剪到研究区逐像元计算类别一致率。差异大的区域导出为掩膜进一步判断是分类体系的差异比如 2018 年新植幼林是否算橡胶还是两类数据系统的系统偏差。曲线的整体走向一致、局部差异集中在边界图斑基本就可以放心使用。这份数据集中最有价值的其实是五年连续分布图本身——比起单年数据它能支撑时间一致性的判别这也是我在本地处理时最看重的能力把任意两年叠起来出变化再把变化结果和外部数据交叉核对一遍得到的结果就可以真正用于文章、报告或模型输入了。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进