陈 涛 周世健 陶 欢 侯艺璇

(1. 东华理工大学 测绘工程学院, 江西 南昌 330013; 2. 南昌航空大学, 南昌 330063;3. 中国科学院地理科学与资源研究所, 北京 100101)

0 引言

油菜是我国最主要的油料作物之一,其种植面积占总油料作物种植面积的50%,产量大约占总油料产量的40%以上[1-3]。目前我国关于油菜种植分布与面积的数据是基于各地区县域统计上报的结果,更新周期较长。遥感技术的特点市快速、无破坏和大面积监测优势,对油菜种植面积监测提供了新的方法。油菜不同生长阶段对应的光谱反射率不同,表现在遥感影像上呈现出光谱信息的差异,为基于多时相的油菜遥感监测提供理论基础。目前已有一些研究结合农作物物候创建时间序列数据,以及植被指数,并将两者应用到农作物提取。在大尺度制图研究选取的影像大部分为中分辨率成像光谱仪(Moderate-resolution Imaging Spectroradiometer,MODIS)数据,王凯等[4]利用2008—2013年75个时相的MODIS-NDVI(归一化植被指数，Normalized Difference Vegetation Index)时序数据,将实地样本数据与油菜物候信息相结合,建立油菜种植面积提取模型,通过多次阈值设定提取了2009—2013年湖北省油菜种植面积。尤慧等[5]基于250 m空间分辨率的MODIS-EVI遥感数据,TM数据作为野外采样数据与MODIS EVI数据之间的过渡数据,最终,经过阈值设定提取2014—2015年间江汉平原油菜种植面积分布。在小尺度制图方面,基于高空间分辨率影像数据,从单一影像向时间序列影像数据发展。柴振刚等[6]为提高Sentinel-1 SAR(合成孔径雷达,Synthetic Aperture Radar)数据作物种植分布提取精度,以湖北省江陵县为研究区域,基于资源三号卫星CCD(电荷耦合元件，Charge-Coupled Device)融合数据(空间分辨率为2 m),对田间边界进行提取,采用SAR数据进行阈值分类、通过SAR不同时期的影像数据进行田间信息提取的比较,最终提取油菜种植分布信息。肖善才等[7]基于Landsat-8和资源3号影像,利用决策树方法,提取南京市2016年油菜种植分布信息,并进行验证。韩涛等[8]利用Sentinel-2A、Landsat-8影像以及野外调查数据,基于影像的光谱特征与植被指数信息利用不同分类方法提取南京市高淳区的油菜种植面积。

由上述研究可知,在遥感数据的选择上,从单一影像到时间序列影像的发展。基于单一影像农作物提取的方法,该方法在确定关键物候期时选取监督分类的方法对农作物进行提取;基于时间序列影像数据农作物提取方法包括:有基于单一特征参量的方法(阈值法-NDVI或EVI)、基于多特征参量的方法,也常用阈值法,除对农作物NDVI或EVI时间序列曲线构建之外,还应该考虑各个波段的影响[9-10]。因油菜种植结构较为复杂,必须对生长期的状况与NDVI进行实时监测,所以本文选取时间分辨较高,覆盖范围较大的MODIS影像。

2015年,湖南省油菜种植统计面积为1.31×106hm2、总产量210.8×104t,分别占全国比重为17.4%、14.12%,居全国前列[11]。本研究以湖南省为研究区,研究的内容包括:(1)利用MODIS-NDVI时序数据和油菜实地采样调查数据,获取油菜像元MODIS-NDVI的时序图,构建湖南省油菜生长期的NDVI时序标准曲线;(2)采用最小二乘法计算湖南省时序影像的每个像元与标准曲线的相似度,相似性越高则认为该像元是油菜的概率越高;(3)基于湖南省各区县统计上报的油菜种植面积结果,对油菜种植面积分布的相似度阈值进行优化,获得湖南省油菜种植的空间分布,并估算在该最优阈值下油菜种植的空间分布误差。

1 研究区域概况

湖南省(24°38′～30°08′N, 108°47′～114°15′E)地处云贵高原向江南丘陵和南岭山脉向江汉平原过渡的地带,地势呈现三面环山、包括平原、河湖等地貌,主要以丘陵为主,是我国油菜种植的优势产区。全省年平均气温20℃左右,冬季最冷月(1月)平均温度在4℃以上。春、秋两季平均气温大多在16～19℃之间。夏季平均气温比春、秋气温高10℃左右,年降水量大约在1 500 mm,属亚热带季风气候。油菜播种气温与之相差不大,因此研究区非常适合油菜的种植。

2 数据来源与研究方法

2.1 数据获取与预处理

本文采用的数据来自美国NASA官网提供的2017年MOD13A1 16天合成的空间分辨率为500 m的产品。将所获得的30景影像利用MRT软件进行拼接、选取,把所有影像的数学基础都统一为GCS_WGS_1984的地理坐标系,投影到WGS_1984_UTM_Zone_50N投影坐标系上,最终进行重采样成500 m。土地利用类型数据使用的清华大学宫鹏团队发布的2017年全球10 m的土地类型数据(http:∥data.ess.tsinghua.edu.cn/)[12]。在ArcGIS中,将13景影像拼接、裁剪,并重采样至500 m分辨率,最终提取出耕地类型图,而湖南省县域油菜种植面积通过统计年鉴获得。

油菜NDVI标准物候曲线的获取通过湖南省长沙市、常德市、衡阳市综合实验站获取GPS经纬度坐标,调查地点的主要分布在长沙市、衡阳市、常德市、怀化市。样点选取原则是尽量保证油菜端元是纯净的,对油菜品种具有代表性的样点进行选取[13]。在经过Google Earth确定农田信息侯,通过ArcGIS得到矢量数据并投影成与MODIS数据坐标相一致,最终从遥感数据中提取出对应像元的NDVI作为油菜的端元,确定湖南省油菜NDVI标准曲线。

2.2 研究方法

2.2.1植被指数NDVI

构建NDVI[14],它在遥感影像进行植被覆盖与植被物候研究得广泛应用,对植被生长状态、植被覆盖度都有着较高的监测效果,并在消除大气、地形以及其他辐射干扰具有很好的作用。

2.2.2 Savitzky-Golay滤波

在光谱分析中平滑滤波是常见的预处理方法之一。利用Savitzky Golay方法的优势在于可以降低噪音的干扰、提高光谱的平滑性[15-16]。而且S-G滤波最大的特点在去除噪声之后还能保证信号的形状和宽度不变。

2.2.3与标准曲线判定相似性

通过综合实验站的实地数据提取相应遥感像元,形成标准的油菜端元NDVI时序曲线图,随后利用ENVI(完整的遥感图像处理平台,The Environment for Visualizing Images)中的bandmath进行相似性计算,一次来判断每个像元的NDVI时序曲线与标准曲线的相似性,确定与标准像元相似的程度。如计算公式(1):

(1)

式中:xi为10个时相对应的NDVI值;yi为标准曲线中NDVI的值;min的值表明像元NDVI曲线与标准NDVI时序曲线在一定值域范围内的拟合程度。

2.2.4相关系数及指标

为评估最终提取结果的精度与效果,本文利用皮尔逊相关系数(r)、决定系数(R2)和均方根误差(RMSE)进行评估,具体公式参考Carlos Antonio da Silva Junior[17]。

本文的研究思路:利用清华大学宫鹏团队发布的2017年全球10 m的土地类型数据[12]从中提取出湖南省耕地类型数据与湖南省的时序,MODIS-NDVI数据进行掩膜提取,得到耕地上的NDVI;再与实地的采样点数据(滤波之后)进行匹配;最终将每个像元的曲线与标准曲线之间进行最小二乘计算,寻找最优相似度阈值的范围,从而得到相应的油菜面积。具体的技术路线如图1所示。

3 结果分析与讨论

3.1 遥感提取面积计算及分析

通过实地采样数据的油菜端元进行取平均值,得到油菜主要生长期内的NDVI时间序列曲线(图2)。NDVI等于0.4一般作为植被生长期的开始[18]。图2中可以看出油菜的NDVI的开始值在0.43,冬油菜NDVI时序曲线整体趋势,1～2月处于油菜的蕾台期,油菜将会停止生长,NDVI变化很小;3～4月正值油菜开花期叶绿素增加,NDVI曲线快速上升;5月处于成熟收获期,NDVI曲线开始下降,而在儒略日081(4月7日)时NDVI最高,是油菜提取最佳时相。根据Savitzky-Golay滤波后的图像,剔除异常值使曲线变得更加的平滑。

图2 油菜NDVI时间序列曲线

经过与标准曲线进行最小二乘计算,图3是湖南省耕地中与标准油菜曲线的相似性分布图,每取一次概率阈值就会得到相应的油菜分布,将其与统计年鉴数据进行比较,确定最优的取值范围为0.2～0.7对应的油菜分布如图4所示,此时的相关性最高。研究区采用的是空间分辨率为500 m的MODIS数据,因此一个油菜像元的面积对应的是地面上的500 m×500 m区域,只需要统计湖南省油菜像元数量即可计算油菜的种植面积,计算公式如下:

(2)

式中:S为提取面积,单位hm2;N为研究区作物提取像元数量;R为MODIS数据空间分辨率(500 m)。

利用ENVI中的分区统计工具对提取结果进行计算,得出共有137 732个油菜像元,根据公式(2)计算,遥感提取湖南省油菜种植面积为3.87×106hm2。

3.2 遥感提取面积与统计面积相关性分析

为了确定各个市级县区油菜提取面积对应相关性,分别将每个市级区县矢量图与湖南省提取结果做掩膜处理并进行分区统计。然后将分区统计的各市级县区的结果与统计年鉴数据做相关性分析。经Pearson相关系数计算,14个市中3个市的县级相关系数在0.95以上。3个市的县级相关系数在0.90～0.95之间,3个市的县级相关系数在0.75～0.85之间,5个市的县级相关系数在0.35～0.7之间。可得出,大部分市级油菜提取面积与统计年鉴数据的相关性较高,表明各市能达到一定精度要求。

图4 最优参数下的油菜种植分布注：审图号GS(2019)3266号。

3.3 精度验证及误差分析

本文利用遥感方法提取的湖南省2017年油菜面积与统计面积作验证(数据来源于湖南统计年鉴),如图4所示,根据相似性阈值范围进行取值,然后利用湖南省统计年鉴数据(县级)与遥感提取面积数据的相关性分析结果可知,当阈值取0.7时相关性最高,以此进一步确定阈值范围在0.2～0.7之间。遥感提取面积的相关系数(r)达到0.81,其值域等级为强相关[19],表明提取结果有较高的可靠性。从湖南省县级油菜遥感提取结果与统计年鉴数据的二维散点图所示(图5),R2为0.66,RMSE为31.4。精度分析表明本文方法在单一农作物遥感提取中有良好的精度。

图5 MODIS-NDVI提取结果与统计年鉴数据二维散点图

表1 湖南省县级统计年鉴数据与MODIS遥感提取面积的差值 单位：hm2

从单一影像发展到基于时间序列影像,利用油菜主要的物候期与MODIS时间序列数据相结合[9-10,20-21],通过遥感影像分析获得的主观数据和客观数据之间的比较可能会导致两种方法之间的差异,造成这些差异的可能是南省油菜种植期多云多雨,并插花种植严重,加之地块破碎,极易产生混合像元的现象[17,22-25]。为此对遥感结果与统计结果差异较大区域进行分析,如表1所示,其中宁乡市的油菜面积被高估了1.04×105hm2;而在怀化市、长沙市、张家界市、郴州市四个地区中有十个县区出现低估的现象,特别是在安化县表现出低估值最大最明显。

4 结束语

本文利用MODIS-NDVI数据描述了我国冬油菜在湖南省的空间格局分布,结合研究区的物候历,构建最小二乘的方法,分析了冬油菜在湖南省的空间格局,并提取冬油菜面积,同时利用统计年鉴数据对MODIS遥感提取结果进行精度验证与相关系数检验。得出以下结论:MODIS数据具有高时间分辨率、光谱范围广、更新频率高,覆盖了农作物的生长期,为大尺度制图提供了新的数据源。研究结果显示,与统计年鉴数据相比,县级相关性系数(r)达到0.81,R2达到0.66。此方法操作流程简单,因此选择合适的MODIS时间序列影像数据,为油菜面积遥感提取提供了可行性。