田 苗 单 捷 卢必慧 黄晓军

(江苏省农业科学院农业信息研究所， 南京 210014)

0 引言

及时准确获取水稻种植面积及空间分布，对水稻生产管理、农业政策制定及粮食安全分析等具有重要意义。目前水稻面积的遥感提取方法主要有常规分类法和基于MODIS的时间序列方法。常规分类法一般使用中、高分辨率的遥感数据[1-3]，而此种数据重访周期长，易受天气影响，不易抓住水稻生长的物候信息，因而，难以进行大范围的水稻面积业务化运行。MODIS数据具有较高的时间分辨率和较广阔的空间覆盖，而且可以免费获取，所以适合大区域水稻面积提取和推广应用[4-5]。基于MODIS时间序列的方法往往根据水稻的物候历，确定水稻识别的移栽期、生长期和收获期等关键时期，根据水稻在不同时期表现的遥感特征指数差异来识别水稻[6-11]。较少有学者考虑与水稻生产能力密切相关的季节积分和生长季振幅两个物候指标。经过众多学者对植被物候的研究，开发了许多判断植被关键生育期的方法，主要包括植被指数阈值法[12]、植被指数中点法[13]、最大变化斜率法[14]、植被指数滑动移动平均法[15]、通过NDVI时间序列累积生长天数来计算物候信息[16]、植被物候期的频率分布与遥感数据相结合法[17]等。目前图像处理已经进入人工智能时代，机器学习回归算法的高度适用性和通用性使其比传统的拟合函数能更精确地估计物候趋势。

本研究以江苏省为研究区域，应用机器学习算法对MODIS-EVI时间序列进行重构，并提取水稻物候信息，选取季节积分和生长季振幅两个指标，采用阈值法提取水稻种植面积，并应用统计数据和Landsat8 OLI影像对提取结果进行精度验证。

1 研究区域与数据获取

1.1 研究区域概况

江苏省地处亚热带向暖温带过渡的区域，其自然条件对水稻生长有利，是全国水稻主产省之一。江苏省南北纵跨5个纬度(30°45′～35°20′N)，东西横跨5个经度(116°18′～121°57′E)，其中，北部地处我国籼稻种植最北边缘，南部地处我国两熟制粳稻种植最南边缘，东部受海洋气候影响，西部受陆地气候影响，形成了境内多种稻作生态区、多种类型稻作熟制和稻作方式、籼粳交错和多种熟制品质并存的水稻生产格局，在全国稻作生产中具有特殊意义。

1.2 数据获取与预处理

本文采用2019—2020年地表反射率数据(MOD09A1)作为基础数据，数据来源于美国国家航空航天局网站(https:∥modis.gsfc.nasa.gov/)，MOD09A1数据的空间分辨率为500 m，时间分辨率为8 d。使用MODIS重投影工具(MODIS Re-projection tool, MRT)对数据进行重采样、格式转换、拼接和重投影等预处理。采用中高空间分辨率的Landsat8 OLI遥感影像数据进行水稻信息目视解译，选取水稻建模样点及精度评估样点，建立水稻面积提取模型并对提取结果进行精度验证。统计数据来自2019年和2020年江苏省农村统计年鉴，包括2019年和2020年江苏省各地级市水稻种植面积。

江苏省所在的MODIS影像序列号为h27v05和h28v05。文献[18]研究证明覆盖作物关键生育期的时间序列长度即可满足一定分类精度，因此，本研究选择的影像时期为2019年和2020年6月2日—11月9日(DOY153～DOY313)，包含水稻移栽期至成熟期。对得到的图像进行波段运算得到EVI时间序列数据集。EVI计算公式为

(1)

式中ρnir——近红外通道反射率

ρred——红色通道反射率

ρblue——蓝色通道反射率

L——土壤调节参数，取1

2 研究方法

2.1 水稻物候提取

基于遥感手段监测植被物候涉及两个重要步骤，一是植被指数的重构，二是植被物候参数提取方法的选择。目前，图像处理已进入了人工智能时代，机器学习算法已成为时间序列数据处理的有力工具[19-21]。机器学习回归算法(Machine learning regression algorithm, MLRAs)的高度适应性和通用性使其比传统的拟合函数能更精确地估计物候趋势[22]。高斯过程回归(Gaussian process regression，GPR)[23]作为一种最优机器学习回归算法在成功重建植被指数和提取可靠的物候指标方面具有潜力。本研究应用GPR算法对时间序列进行重构，在此基础上提取物候信息，提取方法为：生长季开始日期/生长季结束日期在曲线的左/右部分达到曲线上升/衰减部分的季节振幅的30%时被识别。物候指标包括：生长季振幅、最大值、最大值日期、生长季开始日期(Start of season, SOS)、生长季结束日期(End of season, EOS)、季节积分、生长季长度。

在生长发育期识别的各个指标中，季节积分和生长季振幅是两个非常重要的指标，不仅显示了水稻的生育特点，还能够反映水稻的生长状况[24]。经过反复实验，本研究选用这两个指标监测江苏省水稻种植面积。

2.2 水稻种植面积提取

本研究采用阈值法开展水稻面积提取。阈值法是图像分割中的经典方法，它利用图像中要提取的目标与背景在灰度上的差异，通过设置阈值来把像素级分成若干类，从而实现目标的提取。最优阈值选择方法通常有人工经验选择法、直方图法和最大类间方差法[25]。本研究阈值确定方法选用直方图法，结合Landsat8选取的水稻样点和统计面积数据，分别确定江苏省13个地级市季节积分和生长季振幅的阈值。确定步骤为：

(1)以高分辨率遥感影像(Landsat8 OLI)为基础，以目视解译的方式在每个地级市选择若干个水稻样点(图1)，并提取各个水稻样点对应的MODIS-EVI时间序列。

图1 2019年水稻样点分布图Fig.1 Distribution map of rice sampling points in 2019

(2)应用MODIS-EVI时间序列提取每个水稻样点的物候信息，得到每个水稻样点的季节积分和生长季振幅指标。

(3)以地级市作为分区，应用每个区域的若干个样点，做该区域水稻的季节积分频数分布图和生长季振幅频数分布图。

(4)根据频数分布图，初步确定季节积分和生长季振幅的阈值，并将该值用于提取该区域水稻种植面积。

(5)将步骤(4)提取的水稻种植面积与该区域水稻种植面积的统计数据对比分析，根据频数分布重新调整两个指标的阈值，使水稻种植面积的总体误差控制在±10%以内，并将该阈值作为该区域的水稻提取阈值。

(6)应用2019年的数据，经过步骤(1)～(5)，确定各个地级市水稻提取的阈值，应用该阈值提取2020年的水稻种植面积，并对2020年的提取结果进行精度评价。

水稻面积提取精度评价以两种方式进行：与统计面积数据进行对比分析；通过高分辨率影像选取水稻和非水稻样点，建立混合矩阵对提取结果进行定量分析。

3 结果分析

3.1 水稻物候信息提取

本研究以2019年的数据为建模数据，2020年的数据为精度验证数据，采用2019年中、高空间分辨率的Landsat8 OLI遥感影像数据进行水稻信息目视解译，在江苏省内选取了432个水稻样点(图1)，应用GPR方法将2019年水稻样点时间序列(6月2日—11月9日)进行重构，重构后的时间序列时间间隔为1 d(图2a)。图2a显示GPR模型可以很好地将时间序列中的异常值和缺失值进行重构，使曲线更加平滑合理，在此基础上提取水稻物候信息，物候提取参数生长季振幅设置为30%(图2b)。

图2 单点时间序列重构及水稻生育期提取Fig.2 Single-point time series reconstruction and extraction of rice growth period

图2b显示了水稻物候提取参数，其中横坐标为日期，第1天为时间序列开始日期，即6月2日，通过与江苏省水稻物候数据对比，该模型提取的SOS和EOS分别对应江苏水稻的分蘖期和成熟期。按以上方法提取所有水稻样点的物候信息，并分区域整理所有样点的季节积分和生长季振幅。

在水稻样点物候信息提取的基础上，逐像素提取2019年和2020年水稻物候信息，得到2019年和2020年生长季振幅和季节积分分布(图3)，从图3中可以看出，生长季振幅和季节积分的分布具有很好的区域性，水体和城市所在地区的值都较低，能够很好的与周围地物分开，植被覆盖区显示也较清晰，有利于水稻物候信息提取。

图3 江苏省2019年和2020年生长季振幅和季节积分Fig.3 Amplitude and area total of Jiangsu Province in 2019 and 2020

3.2 水稻种植面积提取

分析2019年432个水稻样点的生长季振幅和季节积分，通过直方图和水稻面积的统计数据确定13个市两个指标的阈值(图4)，图4显示，生长季振幅和季节积分阈值的分布也具有区域相似性，如南京市、苏州市、无锡市、镇江市、常州市，两个指标的阈值均相似，淮安市和宿迁市也具有相似性，连云港市、南通市和盐城市具有相似性等。

图4 江苏省13个市的生长季振幅和季节积分阈值Fig.4 Amplitude and area total thresholds of 13 cities in Jiangsu Province

将阈值相似的区域合并(图5)，图5充分体现了江苏省水稻种植的特点：东部受海洋气候的影响(连云港市、盐城市和南通市)，其水稻种植具有一定的相似性，西部受陆地气候的影响，形成了多种稻作生态区，南北纵跨5个纬度，南北水稻种植方式和格局有较大差异。徐州市作为江苏省最北端的地级市与苏南5市的差异较大，苏中水稻种植情况也具有一定差异性，徐州市、连云港市和盐城市等苏北地区也具有一定的相似性(图4)。因此，图4和图5从侧面反映了模型建模的合理性。

图5 阈值合并后的分区Fig.5 Partition map after threshold merging

根据2019年数据确定13个市的生长季振幅和季节积分的阈值，分别提取2020年13个市的水稻种植面积及空间分布，合并得到江苏省2020年水稻种植面积及空间分布(图6)。从图中可以看出，水稻的分布具有很好的区域性，能够较好地将水体、城市和山区等剔除，符合江苏省水稻种植的分布情况。

3.3 精度验证

首先，应用2020年江苏省13个地级市水稻种植面积的统计数据与MODIS提取数据进行对比分析(表1)。表1显示2020年提取的全省水稻总面积的相对误差为-6.10%。在13个市的面积统计中，连云港、淮安、宿迁、徐州、盐城、扬州、泰州、南通和苏州9个市的相对误差都在±10%以内，南京、镇江、常州、无锡4个市的相对误差较大。经分析，可能原因有：①在建模过程中，苏南4市选取的样点较少，阈值的确定不够准确。②苏南4市水稻种植面积较少，在绝对误差相同的情况下，相对误差较高。③苏南地形复杂，水稻种植情况不够稳定，每年的种植情况变化较大(从统计面积上看，2020年苏南4市的水稻种植面积较2019年减少了近10%)，且水稻种植的空间分布情况变化也较大，因此，应用2019年的数据建立的模型，不适用于2020年的水稻面积提取。④在本研究中使用的MODIS影像空间分辨率为500 m，即一个像元内存在像元“同物异谱”和“同谱异物”的现象，对面积较小、地块破碎的地区提取精度不够，从而影响水稻的提取精度。除苏南4市外，江苏省大部分区域水稻种植面积随时间变化具有一定的稳定性，可以用该模型进行水稻面积的提取。在后续研究过程中，还可以不断地加入历史数据或最新数据，不断更新调整阈值，从而更准确地把握每个区域的水稻种植特点，使参数具有稳定性、代表性和可推广性。

表1 2019年和2020年江苏省水稻面积提取结果与农业统计数据对比Tab.1 Comparison of rice area extraction results and agricultural statistics in Jiangsu Province in 2019 and 2020

其次，仅仅依靠统计面积的比较进行精度验证，在理论上存在一定程度的不可靠性，因为相同的面积可能会有不同的形式组合，用高分辨率的遥感影像验证大尺度(低分辨率)提取精度，这种方法的优点在于省时省力效率高，弥补缺乏野外验证的不足。本研究利用Landsat8和目视解译方法选取各个地类的样点共3 423个，通过混淆矩阵对提取结果进行定量分析(表2)，分析结果表明，水稻的制图精度为92.90%，用户精度为89.09%。计算得水稻提取的总体精度为92.55%，Kappa系数为0.846 3。

表2 混淆矩阵精度评价Tab.2 Confusion matrix accuracy evaluation

将非水稻样点分为水体、建筑、其他植被及其他4类，分析各个类别的提取精度(表3)。其中水体、建筑及其他类型由于与水稻差异较大，精度分别为99.52%、98.11%和92.33%，而其他植被由于与水稻具有相似性，因此错分为水稻的概率较大，其精度为87.73%。

表3 非水稻类型精度评价Tab.3 Accuracy evaluation for non-rice types

在精度评价方面，尽管已经开展了大量工作，但是对这种光谱指数时间序列水稻面积提取方法的精度评价还不够充分，本研究主要是针对统计年鉴统计数据进行的行政区尺度水稻面积评估，虽然本研究也应用高空间分辨率影像，提取样点进行了精度验证，但仍不能全面完整地体现面积提取方法的空间分布精度，对于水稻面积总量漏分和误分的空间分布的精度评价工作还有待进一步深入研究。

4 结论

(1)基于机器学习的GPR模型能够较好地重构MODIS-EVI时间序列数据集的异常点和缺失值，为后续水稻物候期的提取提供了较好的基础。根据提取的水稻物候信息，提取水稻种植面积，精度较好，与统计数据相比，大部分地级市的水稻提取精度都在90%以上。与高分辨率影像数据相比，水稻提取的总体精度为92.55%，Kappa系数为0.846 3，水稻的制图精度为92.90%，水稻的用户精度为89.09%。由于其他植被与水稻生长具有一定的相似性，其他植被类型的错分率较高，其精度为87.73%，水体、建筑和其他类型的错分率都较低，精度都在90%以上。

(2)江苏省水稻种植方式、种植格局有较大差异，本研究分13个地级市确定水稻物候指标的阈值，通过阈值相似性合并，很好地体现了江苏省水稻种植的特点。该结论一方面说明了本研究的合理性，另一方面说明对大区域水稻种植面积的分析与建模，需要分区域进行。

(3)本研究仅应用2019年的数据作为建模数据，将模型用于2020年水稻面积提取，并取得了较好的结果。