基于多时相哨兵影像和随机森林算法的冬小麦提取

一、简介

小麦是世界三大粮食作物之一，也是我国主要粮食之一。通过对冬小麦面积的提取，可以知道冬小麦的分布区域，面积大小以及为冬小麦的产量估测做准备。本文将基于GeoScene Pro软件通过传统机器学习算法结合多实相NDVI数据来实现冬小麦面积的精准提取。

二、影像数据介绍

影像数据使用的是哨兵2经过大气校正之后的影像产品，该数据包含空间分辨率10、20和60米的数据。本课程用到的是10米的数据，其中10米的数据主要用到蓝、绿、红、近红外以及RGB图像数据。研究区域选择的是安徽省宿州市丰县。数据的时间段为2019年10月-2020年6月，即冬小麦整个生育期。其中受云层遮挡影响，期间可用影像共11景，且分布较为均匀。

三、冬小麦提取依据

基于随机森林算法和多时相影像数据提取冬小麦，最重要的依据就是，冬小麦多时相NDVI的特征曲线的独特性。如图3，蓝色是冬小麦的整个生育期的特征曲线，红色为冬小麦的整个生育期其他地表覆盖物的NDVI特征曲线。可以看出冬小麦在整个时期NDVI的值变化趋势和变化范围与其他地物相比，差异很大。我们将利用这种差异来分类冬小麦和其他地物。

查看特征曲线的流程：

• 右击多通道的NDVI图层（该图层下文会介绍如何得到），选择创建图表，点击光谱图如图1操作所示。

  

图1 创建光谱图

• 在打开的图标属性中选择点，如图2。

  

图2 光谱图

• 加载RGB模式的图层作为参考，利用上一步选择的点工具，点击图层上各类地物像素点，即可可视化显示出各类地物的特征曲线。如图3。

图3 多时相NDVI特征曲线，蓝色为冬小麦，红色为其他

四、技术路线

这个小结先总结一下冬小麦整个提取流程。

• 下载整个冬小麦生育期的可用影像数据。
• 计算多有影像的NDVI。
• 合并所有单通道的NDVI影像到多通道的NDVI影像，便于样本点向量的提取。
• 船舰矢量文件，选取冬小麦和其他地物两类样本。
• 利用随机森林模型训练和推理冬小麦。
• 计算模型的Kappa系数，确保模型可用。
• 推理结果和官网公布结果对比。

图4 基于随机森林算法提取冬小麦面积

五、具体流程

1)   GeoScene Pro的开启

• 启动 GeoScene Pro
• 如果出现提示，请使用您获得许可的帐户进行登录，如果没有许可，可以联系客服购买。
• GeoScene Pro 随即打开启动画面。选择地图来新建一个工程，输入名称和路径之后点击确定。

图5 新工程的创建

2)   影像的加载和NDVI的计算

• 打开哨兵影像所在的文件夹，如图6，选中b04和b08两个波段直接鼠标拖入图7中的任意一个红色框中即可加载。

图6 哨兵2影像10m分辨率影像文件夹

图7 影像的拖入

• 依次选择影像-栅格函数-数学分析-计算器，来打开栅格计算器。如图8。

图8 打开计算器的操作指示

图9 NDVI计算

• 设置参数：1，在栅格变量左边输入red变量名，右边输入对应的波段号为b04的影像。2，再次以同样的方式输入nir和对应的波段为b08影像。3，在表达式一栏中输入计算NDVI的表达式(nir -red) / (nir + red)。4，点击创建图层来生成NDVI。之后依次计算所有时期的NDVI。

3)   合并所有的NDVI波段并裁剪出预测的区域

多波段的合成是为了便于软件直接提取每个像素点对应的所有波段的NDVI的值。裁剪是为了便于研究某个区域的冬小麦种植分布情况。

• 波段合成，打开波段合成工具。如图10。

图10 工具所在位置

• 参数设置：1，输入栅格中，输入所有的NDVI单波段影像。2，输出栅格中，设置输出路径和名称。3，点击运行，输出合成后的多NDVI波段影像。如图11。

图11 NDVI波段合成

• 区域裁剪，拖拽的方式加载丰县的矢量行政边界。打开按掩膜提取工具。如图12。

图12 工具所在位置

• 参数设置：1，输入栅格，输入待裁剪的影像（本课程指多通道的NDVI图）。2，输入栅格或要素掩膜数据，输入丰县的矢量边界文件。3，输出栅格，输入裁剪结果的路径和文件名，后缀为.tif。如图13。

图13 按掩膜提取工具

4)   选取样本

• 打开哨兵2自带的RGB影像，作为参考底图。如图14。

    

图14 哨兵10m影像文件夹

• 创建两个矢量面文件，一个作为训练样本，一个作为测试样本。其中训练样本是为了训练分类，测试样本为了测试模型的准确性。通过“视图-目录窗格-找到存放文件夹-右击新建shape file文件”打开创建要素类，如图15。

图15 创建矢量文件指示图

• 参数设置：1，要素类位置，输入具体存放路径。2，要素类名称输入名称加.shp文件后缀。3，几何类型选polygon。4，坐标系选择和图像对应的坐标系。如图16。

图16 创建要素

• 对两个矢量（训练和测试）分别开始选取正负样本，即小麦和其他两个类别。以训练集为例，通过选中“创建的矢量-编辑-创建要素-选中待编辑的要素-选择任意多边形”来打开选取样本工具。如图17。

图17 样本的选取

• 样本选取和保存：1，鼠标左击图层，确认一点，之后继续左击，完成时连续左击两次，画完一个面。2，之后分布均匀的画上百个。3，画完之后，点击保存即可。最终程序会在这些图形中选取一个或多个点作为样本集。如图17。

• 增加字段：样本选好之后，右击矢量样本打开属性表来增加字段。目的是为了符合该工具的输入格式。

图18 矢量样本属性表

• 具体操作：

1. 点击属性表的添加，添加“classname”和“classvalue”字段。类型分别是文本和长整型。随后保存。如图19。

2. 选中classname，点击计算，输入‘dxm’，将选中的冬小麦样本命名为‘dxm’，同样的方式，将classvalue改为1。如图20。
3. 如上同样的操作选取负样本（除冬小麦之外任意地物），在属性表中选中所有的负样本（shift+左键），对选中的部分（青色高亮区域）classname设置为‘other’，class value设置为0，最终结果如图21。

图19 添加字段和保存设置

图20 添加字段、值和保存设置

图21 添加字段和字段值计算

5)   随机森林分类提取

• 影像分类，选中要分类的多通道NDVI图（丰县区域），选择影像-分类工具-分类。如图22。

图22 传统机器学习分类工具

• 参数设置：1，分类器，选择随机森林算法（当然也可尝试其他算法，操作相同）。2，训练样本，输入上文选好的样本矢量文件。3，每个类的最大样本数本课程中选了1000。其余参数默认。如图23。

图23 随机森林分类

• 分类结果展示：

图24 左为多时相NDVI图，右为分类结果

6)   计算模型准确性

这一步是为了确保训练的模型是否以达到可以接受的精度。评价指标是Kappa系数。Kappa系数用于一致性检验，也可以用于衡量分类精度，Kappa系数的计算是基于混淆矩阵的。我们将通过创建精度评估点和更新精度评估点来计算混淆矩阵和Kappa系数。

• 打开创建精度评估点工具，如图25。

图25 工具所在位置

• 参数设置：1，输入栅格数据或要素类数据，输入测试样本的矢量文件。2，输出精度评估点，输入保存的路径、名称和‘.shp’后缀。3，目标字段选择实际地表。4，随机点数量选500，也可默认或自己定。其他参数默认即可。如图26。

图26 创建精度评估点

• 更新精度评估点，打开更新精度评估点工具，图27。目的是将识别的结果和对应的自己的标注的结果放在一起，便于计算混淆矩阵。

图27 工具所在位置

• 参数设置：1，输入栅格数据或要素类数据，输入分类结果文件。2，输入精度评估点中，输入精度评估点矢量文件。3，输出精度评估点中，输入更新后的精度评估点路径、名称和后缀‘.shp’。4，目标字段，选择分类。如图28。

图28 更新精度评估点工具

• 打开计算混淆矩阵工具，图29。

图29 工具位置

• 参数设置，1，输入精度评估点中，输入更新精度评估点矢量文件。2，输出混淆矩阵：输入计算结果的输出路径、名称以及后缀‘.txt’。如图30。

图30 混淆矩阵计算结果

• 结果展示：结果是个包含混淆矩阵和Kappa系数的表。本课程计算的Kappa系数为1。模型准确性达标。如图31。

图31 混淆矩阵计算结果

7)   分析结果

• 打开构建栅格属性表工具，图32。对分类的结果创建属性表，图33。目的是为了统计分类结果中冬小麦的像素点个数，进而计算冬小麦的面积。

图32 工具位置

图33 构建栅格属性表工具

• 查看结果，右击分类结果图层，打开属性表，图34，查看值为1的个数，即冬小麦像素点的个数为4256985，其中每个像素点的面积是100 m²。换算单位后约63.9万亩，而官方公布2020年冬小麦面积的64.7万亩，图35。误差百分比0.012。

图34 统计冬小麦像素点

图35 2020年丰县冬小麦种植面积官方数据

文章来源于微信公众号【易智瑞】，原文链接：基于多时相哨兵影像和随机森林算法的冬小麦提取