利用遥感影像计算大蒜种植面积

大家对大蒜应该不陌生，近几年也经常以"蒜你狠"出现在大众视野。我国是世界大蒜的主要生产国、消费国和出口国，从事大蒜生产的蒜农达500万之多,大蒜产品也远销东南亚、东亚、中东、美洲、 欧洲等地区。大蒜的种植面积是大蒜市场行情的重要影响因素之一，利用遥感手段可以快速、直观的监测大蒜种植面积，本文介绍利用遥感影像量测大蒜种植面积的技术流程。

一、技术流程

             基于遥感技术监测大蒜种植面积的技术流程如下图所示：

图1.1基于遥感技术监测大蒜种植面积的技术流程图

（一）收集数据

      数据主要是遥感影像数据，其他的数据还可以包括矢量数据，土地利用分类数据、耕地数据等。

     遥感影像数据需要确定成像时间和空间分辨率，根据信息提取方法还需要考虑光谱分辨率。

• 成像时间：根据大蒜生长周期，北方地区的蒜苗返青是在1月~4月份，选择这个时间段成像的影像数据比较合适。

• 空间分辨率：由于大蒜的种植面积和范围不是很大，影像的空间分辨率至少15米。分辨率越高，越能识别大蒜。

• 如果使用计算机自动分类方法，至少需要3个波段，并且包括红色和近红外波段。

 （二）数据处理

不同遥感图像需要采用不同处理流程，根据数据情况可以包括正射校正、图像配准、大气校正、图像融合、镶嵌与裁剪等。

（三）耕种信息提取

      耕种信息提取包括人工目视解译、计算机自动提取（监督分类、决策树、面向对象）等方法，如果有土地覆盖数据、耕地信息、实地观测等辅助信息可以提高提取精度。

（四）结果分析

      可以通过目视判读，或者利用验证样本对结果进行精度验证。通过统计结果得到种植面积等信息，如果有单产数据，可以粗略估算当年产量。

二、详细操作步骤

第一步：收集遥感数据

      山东省是我国最大的大蒜生产区,种植面积在 13. 3万公顷以上。以山东济宁市金乡县为中心,辐射成武、巨野、定陶、 单县、嘉祥、 鱼台、 微山等周边地区,形成我国最大的大蒜产区。

      这里选择选择免费的Landsat8遥感数据，分辨率为全色影像15米，覆盖金乡县及周边，包括2014年、2015年和2016年的影像，如下表所示：

注：同时也下载了LC81220362016038LGN00影像，通过目视对比，没有LC81220362016070LGN00影像更好区分大蒜。

       大蒜种植面积量测区域为山东济宁市金乡县。

第二步：Landsat8数据处理

         如下图为landsat8的处理步骤，由于三景landsat8成像时天气晴朗，后面的信息提取采用的是基于规则的面向对象提取方法，这里没有进行大气校正。如果成像天气有雾霾等因素影响，可以使用快速大气校正工具（QUick Atmospheric Correction——QUAC）进行大气校正，这样不同地物的NDVI差异较大，能提高分类精度。

图2.2 Landsat8数据处理

       下面以LC81220362016070LGN00数据处理为例、ENVI5.3.1中介绍整个处理步骤，其他版本类似，其他两个数据处理类似。

（1）打开LC81220362016070LGN00_MTL.txt和金乡县矢量文件。

（2）Landat8以RGB：swir1、NIR、Red组合显示。

（3）图层管理中，金乡县矢量图层上右键选择Zoom to Layer Extent，显示金乡县范围内的影像。

图2.3 显示金乡县范围内的影像。

（4）在Toolbox中，打开/Image Sharpening/NNDiffuse Pan Sharpening工具，

• Input Low Resolution Raster：选择多光谱图像，在选择文件对话框中，单击Spatial Subset，在打开的缩微图上，单击第二个按钮"Use View Extent"，以当前视图范围作为融合的空间子区。

• Input Higth Resolution Raster：选择全色波段图像。

• 其他参数默认。

     （5）选择输出路径和文件名。

（6）单击OK执行处理。

注：由于整景Landsat8范围较大，金乡县只占其中一小部分，使用Spatial Subset选择空间子区，可以减少运算量和文件大小。

图2.4 使用Spatial Subset选择空间子区

（7）在Toolbox中，打开/Regions of Interest/Subset Data from ROIs工具。

（8）在文件对话框中选择前面融合结果。在Spatial Subset via ROI Parameters面板中，设置Mask pixels outside of ROI？：Yes。

（9）设置输出路径和文件名。

（10）单击OK执行处理。

图2.5 Spatial Subset via ROI Parameters面板

最后得到金乡县15米的Landat8影像数据。

第三步：耕种信息提取

       如下图为大蒜耕种信息提取流程，使用面向对象信息提取方法。

（1）以RGB：swir1、NIR、Red组合显示上一步裁剪结果。该地区的冬天农作物主要为冬小麦和大蒜，在影像上的颜色比较好区分，如下图所示。

图2.6 面向对象提取大蒜种植面积流程图

（2）在Toolbox中，打开/Feature Extraction/Rule Based Feature Extraction Workflow工具。打开工作流的面板，选择待分类的影像（上一步裁剪结果）。

（3）切换到Custom Bands面板，有两个自定义波段，包括归一化植被指数或者波段比值，勾选Normalized Difference，自动选择band4和band5计算NDVI。单击Next。

图2.7 选择对话框

（4）在Object Creation面板中，设置分割阈值和合并阈值，勾选Preview预览分割和合并效果，通过目视判断，设置：

• 分割阈值（Scale Level）：20

• 合并阈值（Merge Level）：0

注：由于分割边缘是绿色，可在图层列表上右键选择Change RGB Bands…设置RGB：NIR、Red、Gree

（5）其他参数默认，单击Next。

（6）在Rule-based Classification面板中，点击

按钮，新建一个类别，在右侧Class properties下修改好类别的相应属性。在默认的属性Spectral Mean上单击，激活属性，右边出现属性选择面板，如下图所示。选择Spectral，Band下面选择Normalized Difference。在第一步自定义波段中选择的波段是红色和近红外波段，所以在此计算的是NDVI。通过拖动滑条或者手动输入确定阈值。勾选Preview可以预览不同阈值的提取结果。下面我们通过选择样本统计的方法获取分类阈值。

图2.8 Rule-based Classification面板

A、在图层管理器上，Region Mean右键打开View Metadata，复制Dataset的文件路径。

B、重新启动一个ENVI，打开前面复制的路径，同时打开金乡县15米的Landat8融合影像数据。

C、参考landsat15米融合结果影像，目视选择一部分大蒜区域的ROI，在图层管理器的ROI上右键选择Statistics，统计ROI对应的Normalized Difference值范围（Band8）。

注：ROI可以任意加载在多光谱影像或者Region Mean数据上。

（7）将上一步统计得到的阈值输入Normalized Difference属性中，勾选Preview预览提取结果，如果不合适，选择更多的ROI进行统计。

（8）单击Next按钮，在Export面板中，只选择Export Raster栅格输出。

同样的方法提取2015和2014年大蒜种植面积。

三、结果分析

       如下为提取三年的大蒜种植空间分布图。

图3.1 三年大蒜种植区遥感量测结果

       对栅格分类图进行统计，得到三年的种植面积量测结果如下表所示。

      遥感量测的结果跟网络公布的结果有较大的差别，但是从其他参考信息："中华财险山东分公司相关负责人表示，该公司引入卫星遥感，全程监测金乡大蒜长势动态情况。遥感显示，2015年金乡大蒜种植面积69.01万亩"。公布的结果与本次遥感量测的结果相接近。