大家对大蒜应该不陌生,近几年也经常以"蒜你狠"出现在大众视野。我国是世界大蒜的主要生产国、消费国和出口国,从事大蒜生产的蒜农达500万之多,大蒜产品也远销东南亚、东亚、中东、美洲、 欧洲等地区。大蒜的种植面积是大蒜市场行情的重要影响因素之一,利用遥感手段可以快速、直观的监测大蒜种植面积,本文介绍利用遥感影像量测大蒜种植面积的技术流程。
一、技术流程
基于遥感技术监测大蒜种植面积的技术流程如下图所示:
图1.1基于遥感技术监测大蒜种植面积的技术流程图
(一)收集数据
数据主要是遥感影像数据,其他的数据还可以包括矢量数据,土地利用分类数据、耕地数据等。
遥感影像数据需要确定成像时间和空间分辨率,根据信息提取方法还需要考虑光谱分辨率。
-
成像时间:根据大蒜生长周期,北方地区的蒜苗返青是在1月~4月份,选择这个时间段成像的影像数据比较合适。
-
空间分辨率:由于大蒜的种植面积和范围不是很大,影像的空间分辨率至少15米。分辨率越高,越能识别大蒜。
-
如果使用计算机自动分类方法,至少需要3个波段,并且包括红色和近红外波段。
(二)数据处理
不同遥感图像需要采用不同处理流程,根据数据情况可以包括正射校正、图像配准、大气校正、图像融合、镶嵌与裁剪等。
(三)耕种信息提取
耕种信息提取包括人工目视解译、计算机自动提取(监督分类、决策树、面向对象)等方法,如果有土地覆盖数据、耕地信息、实地观测等辅助信息可以提高提取精度。
(四)结果分析
可以通过目视判读,或者利用验证样本对结果进行精度验证。通过统计结果得到种植面积等信息,如果有单产数据,可以粗略估算当年产量。
二、详细操作步骤
第一步:收集遥感数据
山东省是我国最大的大蒜生产区,种植面积在 13. 3万公顷以上。以山东济宁市金乡县为中心,辐射成武、巨野、定陶、 单县、嘉祥、 鱼台、 微山等周边地区,形成我国最大的大蒜产区。
这里选择选择免费的Landsat8遥感数据,分辨率为全色影像15米,覆盖金乡县及周边,包括2014年、2015年和2016年的影像,如下表所示:
数据 | 成像时间 |
LC81220362014080LGN00 | 2014年3月21日 |
LC81220362015019LGN00 | 2015年1月19日 |
LC81220362016070LGN00 | 2016年3月10日 |
注:同时也下载了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 三年大蒜种植区遥感量测结果
对栅格分类图进行统计,得到三年的种植面积量测结果如下表所示。
年份 | 影像数据 | 遥感量测结果 | 省农业厅网站公布结果(来源:金乡大蒜专业批发市场信息中心) |
2014 | 2014年3月21日 | 75.2万亩 | 62.30万亩 |
2015 | 2015年1月19日 | 69.28万亩 | 53.6万亩 |
2016 | 2016年3月10日 | 75.77万亩 |
遥感量测的结果跟网络公布的结果有较大的差别,但是从其他参考信息:"中华财险山东分公司相关负责人表示,该公司引入卫星遥感,全程监测金乡大蒜长势动态情况。遥感显示,2015年金乡大蒜种植面积69.01万亩"。公布的结果与本次遥感量测的结果相接近。