基于MODIS NDVI数据的辽宁省玉米种植面积提取研究

宫诏健,田景仁,陈 杰,朱明宇,宋文锦,朱 红,黄荟羽,张运芝,王 莹,贺 鑫,刘志邦

(1.辽宁省营口经济技术开发区气象局,辽宁 营口 115007;2.辽宁省庄河市气象局,辽宁 庄河 116400;3.辽宁省气象装备保障中心,辽宁 沈阳 110166;4.辽宁省辽阳县气象局,辽宁 辽阳 111200;5.辽宁省营口市气象局,辽宁 营口 115007)

0 引言

玉米是辽宁省第一大种植作物,利用遥感资料准确、及时地提取玉米种植面积信息,对提高玉米综合管理水平、预测区域性粮食产量、优化农作物布局和确保粮食安全等具有重要意义[1-5]。

传统的农作物种植面积获取采用农业调查队统计抽样或参照行政级别逐级上报的方式,这两种传统面积获取方式工作量大,人工成本高,耗费时间长,且很难分析大尺度空间作物种植分布情况[6-7]。遥感提取作物种植面积研究发展的时间虽然不长,但其具有覆盖区域广、获取频率高、精确、快速、经济等优势而逐渐成为农作物种植面积提取研究的重要手段[8]。目前农作物种植面积遥感提取研究多集中在高光谱、高空间分辨率的数据源上[9-10],如Grace K等使用World View高分影像遥感提取了海地和危地马拉作物的种植面积[11]; Mirik M等使用TM影像提取了德克萨斯州冬小麦的种植面积[12];李丹等使用Hyperion高光谱影像遥感获取了广州市荔枝的种植信息[13];余凌翔等使用HJ1/CCD影像遥感提取了西双版纳橡胶的分布信息[14]。高分辨率遥感数据具有较高的作物分类识别精度,但对于区域性作物种植面积提取而言,存在选取范围内遥感数据覆盖不全、获取难、成本高、时效差等诸多壁垒,中、低空间分辨率的多光谱数据在遥感提取区域性作物种植面积研究方面具有明显优势[15-18]。常用于农作物监督分类、种植信息提取的植被指数有归一化植被指数(Normalized difference vegetation index, NDVI)、比值植被指数(Ratio vegetation index, RVI)、垂直植被指数(Perpendicular vegetation index, PVI)、增强植被指数(Enhanced vegetation index, EVI)等[19]。其中NDVI植被指数能够很好地反映出农作物的生长信息,体现不同农作物的物候特征,是目前卫星遥感监督分类作物应用最稳定、最普遍的参数之一[20-22]。使用拟合、滤波法重构植被指数时序曲线提取农作物种植面积是一种较为常用的方法。对于辽宁省春玉米而言,基于非对称高斯函数拟合法重构的NDVI时序曲线保真性能好,基于Savitzky-Golay滤波方法重构的NDVI时序曲线平滑性能好,使用这两种方法重构的NDVI时序曲线都与玉米实际的生长状况接近[23]。

因此,本研究选定辽宁省春玉米种植面积为研究对象,利用S-G滤波法和AG拟合法分别重构春玉米NDVI时序曲线,选用动态振幅阈值法、拐点法和最大值法提取了辽宁省春玉米的关键物候期,并联合决策树分类法对辽宁省春玉米分布进行了监督分类,以期为使用卫星遥感数据提取辽宁省春玉米种植面积提供一种有效的方法。

1 材料与方法

1.1 研究区域概况

辽宁省(东经119°～126°，北纬39°～43°)位于中国东北南部,北与吉林省、内蒙古自治区毗连,西与河北省相接,东有鸭绿江,南邻渤海、黄海。辽宁省总面积为14.8万km2,地势为北高南低,东西两侧丘陵向中部平原下降。辽宁省气候较为宜人,四季分明,属温带大陆性季风气候区,日照丰富,积温偏高,温差较大,降水不均,雨热同期。

辽宁省是中国重要的玉米种植区之一,省内绝大多数地区以种植春玉米为主,年均种植面积约为161.5万hm2,约占全国总种植面积的6.47%;年均产量约为1069万t,约占全国玉米总产量的8.61%[24]。

1.2 数据材料

1.2.1 遥感类数据来源 本研究使用2017年MODIS NDVI数据中的3级陆地标准数据产品MOD13Q1(MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid),其时空分辨率为16 d、250 m,从美国国家航空航天局(NASA)官方网站下载。

1.2.2 统计类数据来源 选用辽宁省气象局11个农业气象站(昌图、沈北、建平、新民、灯塔、建昌、绥中、岫岩、宽甸、瓦房店、庄河，其分布见图1)的观测数据作为辽宁省春玉米物候期遥感提取的验证数据,由辽宁省气象信息中心提供;选用2018年辽宁省统计年鉴提供的统计数据作为春玉米种植面积遥感提取的验证数据。

图1 辽宁省玉米实验区分布

1.3 研究方法

1.3.1 数据预处理 选用NASA研制的MODIS重投影工具(MRT)对2017年MOD13Q1数据产品进行格式转换、基准投影、波段选择和栅格拼接等预处理工作,其具体参数如表1所示。

表1 MRT投影参数设置

选用ENVI软件和辽宁省行政边界矢量对拼接后的数据进行剪裁处理,得到2017年辽宁省春玉米发育期的NDVI栅格数据(图2)。

图2 辽宁省MODIS NDVI栅格数据(2017年5～9月)

以辽宁省气象局11个农业气象站的地理坐标制作美国标准信息交换代码(ASCII)文本,使用ENVI Classic中的ROI Tool模块提取所需的2017年5～9月实验研究点的NDVI栅格信息。

1.3.2 时间序列遥感数据的重构方法 MODIS NDVI数据中的3级陆地标准数据产品的主要合成算法是角度归一化合成方法,能够有效地减少太阳高度角、暗影和残云的影响,科学地提高合成数据的一致性和质量,但该数据产品仍受水汽、臭氧、降水等因素的干扰,使其产生不可忽略的、不合理的波动[25-27]。因此,本研究使用TIMESAT软件提供的以下2种重构方法去除时序曲线的噪声。

1.3.2.1 赛维迪克·格雷滤波方法(Savitzky-Golay, S-G) S-G滤波方法是移动平均算法的一种改良方法,是一种基于最小二乘卷积的移动平均算法,其权重由滤波窗口中做最小二乘多项式拟合而得出的次数决定[28]。滤波公式如下:

(1)

式(1)中:Y为NDVI原值；Y′为拟合值；Ci是第i个NDVI的权重值；N为卷积数目,数值上为滤波窗口宽度的2n+1。滤波效果由滤波窗口大小和平滑多项式的次数控制,滤波窗口影响了重构后曲线的平滑能力及其斜率变化速率;平滑多项式次数的高低会使平滑结果出现过度平滑或平滑次数不足的现象[29]。

1.3.2.2 非对称高斯函数拟合法(Asymmetric Gaussian model, AG) AG拟合法是基于高斯函数的分段最小二乘拟合算法,在运算时函数先提取时序曲线内的极大值和极小值,使用获取的极值将曲线分为多个区间,分别在各自区间内进行局部拟合,再将各段的拟合曲线进行整体拼接。局部拟合函数公式如下:

f(t)≡f(t;c,d)=c1+c2g(t;d)

(2)

式(2)中:线性系数c=(c1,c2)控制曲线的振幅和基值;非线性系数d=(d1,d2,…,di)控制基函数的形状g(t;d)。

基函数公式如下:

(3)

式(3)中:d1为极大、极小值对应变量t的位置参数;d2,d3,d4和d5分别决定了极值左右半曲线的宽度和陡峭度。

1.3.3 参数设置 在重构时序曲线时,重构模型的参数设置对重构结果的优劣有重要影响。参照前人的研究成果,本研究选用如下参数:有效值域0.00～1.00,包络迭代次数3,拟合峰值0.5,去噪声模式2:STL replace,滑动窗口大小4,自适应强度3.0。

1.3.4 玉米关键物候期的提取方法 使用Origin软件中的4次拟合平滑功能对重构后的曲线进行平滑处理,选用动态振幅阈值法、曲线拐点法和最大值法提取辽宁省春玉米的关键物候期(出苗期、拔节期、抽雄期)。选用动态振幅阈值法提取春玉米的出苗期,以春玉米生长季开始前裸地的NDVI值为基值,以距离基值振幅10%的日序为春玉米的出苗期。选用拐点法提取春玉米的拔节期,在玉米营养生长过程中,其生长速率为“慢-快-慢”的“S”型曲线,且其拔节期生长速率最快,因此时序曲线上升阶段的拐点即为春玉米的拔节期。选用最大值法提取春玉米的抽雄期,玉米抽雄期为从营养生长向生殖生长转变的时期,此时玉米的平均叶面积最大,其NDVI时序曲线达到峰值,因此以NDVI最大值的日序为春玉米的抽雄期[6]。

1.3.5 玉米种植面积的提取方法 选用关键期时序提取法获取辽宁省春玉米的种植面积。以遥感反演的春玉米关键物候期(出苗期、拔节期、抽雄期)对应日序临近的2个NDVI数据作为遥感提取目标(B1～B4),从辽宁省11个农业气象站重构平滑后的NDVI数据集中检索提取目标的NDVI值，作为决策树值域。使用ENVI软件中的决策树分类模块对遥感目标进行分类提取,并通过比值计算获得辽宁省春玉米的种植面积。决策树监督分类流程如图3所示。

图3 实验区春玉米监督分类流程

选用极值法分别筛选玉米出苗期和拔节期的决策树值域(T1和T2～T5)。在春玉米抽雄期易发生NDVI过饱和现象,因此选用四分位法筛选玉米抽雄期决策树值域(T6、T7)。决策树监督分类阈值见表2。

表2 决策树监督分类的阈值

2 结果与分析

2.1 重构方法的分析与评价

春玉米原NDVI时序曲线异常值波动明显。本研究选用S-G和AG重构法纠正时序曲线的异常值,使用四次拟合法进行平滑处理,重构平滑后的时序曲线(图4)基本上保持原曲线的形态,且除噪声效果明显。

以2017年绥中农业气象站为例。

本研究使用重构曲线与原曲线的拟合度(R2)和均方根误差(RMSE)来评价重构方法的平滑能力和保真能力，结果如表3所示。从表3可以看出：基于S-G滤波法的各实验站点的平均拟合度为0.983，优于AG重构法的平均拟合度0.956,即S-G滤波法有更好的平滑能力;基于S-G滤波法和AG重构法的各实验站点的平均RMSE分别为0.196和0.197,即这两种重构法的保真能力基本相当。基于对上述两种指标的综合评价,用S-G滤波法重构后的时序曲线综合性能更优。

表3 两种NDVI时序曲线重构方法的效果评估

2.2 春玉米关键物候期的对比分析

2.2.1 出苗期、拔节期和抽雄期的对比分析 分别对使用两种重构曲线提取的春玉米出苗期、拔节期和抽雄期日序与实际观测物候期日序进行比较，结果如图5所示，基于S-G重构法提取的各物候期精度普遍低于5 d,精度劣于5 d的样本为6个(占18.1%),出苗期、拔节期与抽雄期的平均误差分别为3.9、2.9和2.6 d;基于AG重构法提取的各物候期精度普遍低于6 d,精度劣于6 d的样本为4个(占12.1%),出苗期、拔节期与抽雄期的平均误差分别为4.5、2.8和3.2 d。对于春玉米出苗期和抽雄期的提取效果而言,基于S-G重构的时序曲线有更好的提取精度;对于春玉米拔节期的提取效果而言,基于AG重构的时序曲线的提取精度略好。

图5 2017年辽宁省春玉米关键物候期提取值与观测值比较

2.2.2 不同地区关键物候期的对比分析 从表4可以看出，两种重构时序曲线在不同地区的提取精度有明显差异,其中基于S-G重构曲线的关键物候期提取精度在灯塔市实验区最优,其MAE为1.55 d;在宽甸县实验区最劣,其MAE为3.96 d。基于AG重构曲线的关键物候期的提取精度在新民县实验区最优,其MAE为2.31 d;在建平县实验区最劣,其MAE为5.00 d。基于各地的遥感提取效果,S-G重构曲线的平均提取精度为3.11 d，优于基于AG重构曲线的平均提取精度3.45 d。

表4 基于S-G与AG重构法的时序提取精度(MAE)对比 d

2.2.3 玉米种植面积提取分析 基于S-G重构法提取的辽宁省2018年春玉米种植面积为2.54万km2,提取精度为94.35%;基于AG重构法提取的春玉米种植面积为2.4万km2,提取精度为89.29%(图6)。整体而言,基于S-G重构曲线的春玉米种植面积提取精度优于AG重构曲线的面积提取精度。通过ENVI的三窗功能对比发现,两者提取差异主要集中在辽宁东、西部的山区,说明相较于AG重构曲线,S-G重构曲线对海拔高、坡度陡、玉米地块分散破碎、纯净度不高的像元有更好的提取效果。

a：S-G重构法；b：AG重构法。蓝色区域为玉米种植分布。

3 结论与讨论

本研究从春玉米物候信息遥感研究的特殊性出发,使用2017年4～9月辽宁省MOD13Q1数据产品中的NDVI波段构建时序曲线,采用S-G重构法、AG重构法和四次拟合平滑法除噪滤波平滑曲线,应用动态振幅阈值法、拐点法和最大值法提取春玉米关键物候期,并联合决策树监督分类法估算了辽宁省春玉米种植面积。研究结果表明:S-G重构法对时序曲线的平滑能力和综合重构性能更优,S-G重构法和AG重构法对时序曲线的保真能力基本相当。基于S-G重构法提取的出苗期和抽雄期优于基于AG重构法的提取结果,基于AG重构法提取的拔节期略优于S-G重构法的提取结果。基于S-G重构法对各实验区关键物候期的整体提取效果好于基于AG重构法的整体提取效果。基于S-G重构法的辽宁省春玉米种植面积的提取精度为94.35%,优于基于AG重构法的面积提取精度89.29%。整体而言,S-G重构时序曲线的综合性能更优,对关键物候期的提取效果更好,对春玉米种植面积的提取精度更高,因此联合S-G时序重构法和决策树监督分类法可为省级遥感估算春玉米种植面积提供一种新的途径。

本研究在以下3个方面还需继续改进:一是春玉米监督分类精度有待提高。辽宁省中部为平原地区,春玉米种植基本连片,像元的纯净度较高,提取效果好；但辽宁省东西两侧为丘陵,春玉米种植破碎,像元易混入杂物,从而降低了春玉米种植面积的提取精度。今后的研究需采用多元数据融合技术,以提高融合后的MODIS复合数据影像的空间分辨率,获取更高精度的春玉米监督分类,减少同物异谱。二是春玉米提取面积区域需要精细化划分。由于省域提取面积范围较大,受地理位置、环境气候等影响,各地区春玉米的种植品种、物候和长势存在一定的差异。本研究受数据资源的限制,没有分区域精细化遥感提取种植面积,导致部分县(区)玉米种植面积的提取精度下降。因此，今后的研究需构建省内精细化的种植面积提取区域,优化种植面积提取模型。三是数据产品的内在性能有待优化。由于NDVI数据产品易出现饱和现象,致使春玉米成熟期附近的监测效果下降,因此今后的研究需通过对比分析增强型植被指数(EVI)数据产品来克服NDVI数据产品的不足。