基于多时相Sentinel-2A的县域农作物分类

2019-10-10吕玉娜李纯斌李全红

农业机械学报 2019年9期

吴静吕玉娜李纯斌李全红

(1.甘肃农业大学资源与环境学院，兰州 730070； 2.甘肃农业大学管理学院，兰州 730070)

0 引言

农作物种植种类、面积和产量估算是国家粮食生产安全和经济运行的重要监测指标。遥感技术由于能快速获取作物类型和分布，已成为大尺度范围农作物监测的主要手段之一[1]。农作物的遥感精细分类是精确估产的前提，是目前农业遥感的研究热点之一[2-5]。

充分利用农作物的典型反射光谱特征和季相节律特征是区分不同作物的关键理论依据[6]，也是利用光学遥感手段进行农作物分类的主要着力点。国内外相关研究多采用农作物生长发育过程中多时相的光学遥感数据提取时间序列的植被指数曲线来分类。采用的遥感数据包括MODIS[7-11]、Landsat[2-3,12-14]、HJ卫星[15-17]、GF-1[18-20]等以及多元数据融合[21-24]，计算的植被指数包括NDVI、EVI等；提取的农作物信息包括单一种类农作物识别及分布信息[3-4,7,14,17]，以及多种农作物种植结构[18-27]。

图1 景泰县地理位置和行政区划Fig.1 Location and administrative zoning map of Jingtai County

近年来，国内对多种农作物分类的研究主要集中于影响因素可控的小区域(如农场)[1,15-16,21]或地势平坦的东北平原[7-10,19,22,24]，在地形较复杂的县域尺度上进行农作物分类有一定的挑战性，因为可能遇到地块较为破碎，分布比较零散、作物种植结构相对复杂等问题。

Sentinel-2A卫星能够提供10 d重访周期和最高10 m空间分辨率的多光谱数据，尤其是其红边波段数据为农作物类型识别和产量监测提供了强有力的数据支持[26]。目前对于利用Sentinel-2A数据的时空优势及其时序特征实现农作物提取的相关研究正在陆续展开[1,26-28]。然而农作物种植结构具有强烈的区域性，不同的区域地形地貌、气候、水文特征不同，农作物的种类、分布、物候特征也不同，Sentinel-2A数据的全面应用需要以实际案例研究为基础，探讨其在具体区域农作物分类方面的能力和应用潜力。本文以西北内陆温带干旱区的景泰县为研究区，在实地调查研究区主要作物及分布的基础上，结合农作物物候信息、土地变更调查数据，利用2018年多时相Sentinel-2A数据，计算时序NDVI和RENDVI指数及其组合为特征参数，采用随机森林法对研究区的主要农作物进行分类，并结合实地调查选取验证数据进行结果验证，对比分析不同特征参数提取作物信息的精度，以探讨多时相Sentinel数据的不同特征参数应用于干旱区县域农作物分类研究的可行性，以期为当地农作物种植结构信息的提取提供新思路。

1 研究区与数据源

1.1 研究区概况

景泰县位于甘肃省中部(103°33′～104°43′E，36°43′～37°38′N)，河西走廊东端，甘、宁、蒙三省(区)交界处，是黄土高原与腾格里沙漠过渡地带，如图1所示。全县面积5 483 km2，境内海拔1 276～3 321 m，地势由西南向东北倾斜，地貌类型以倾斜平原为主，地势相对平坦。景泰县属于温带干旱型大陆气候，年均降水量185 mm，年均蒸发量3 038 mm；年日照时数2 652 h，日照百分率60%，年平均太阳辐射量619 kJ/cm2，年均温度8.2℃，大于等于0℃年活动积温3 614.8℃，大于等于10℃年有效积温3 038℃，无霜期191 d，适宜一年一季的农业生产[29]。黄河流经县境110 km，境内有被誉为“中华之最”的景电一、二期高扬程提灌工程，农业生产条件较为便利，是《全国新增1 000亿斤粮食生产能力规划(2009―2020年)》中确定的甘肃省7个产粮大县之一。

本研究为减少其他地类的干扰、提高作物分类精度，首先提取出耕地信息，然后在耕地范围内提取作物信息。农作物主要有春小麦、玉米、水稻、秋油葵、胡麻、洋葱、大棚作物、马铃薯、砂田瓜果9种类型，其中，春小麦-秋油葵是景泰县传统轮作模式，除此之外，其他作物种植都属于一年一季。各种作物物候期见表1。

表1 景泰县主要农作物物候期Tab.1 Crops development period of Jingtai County

注：表中数据来源于实地走访调查及参考中国气象数据网(http://data.cma.cn/)

1.2 遥感数据及预处理

Sentinel-2A携带一台多光谱成像仪(MSI)，幅宽达290 km，重访周期10 d。从可见光和近红外到短波红外可覆盖13个光谱波段，最高空间分辨率为10 m，在光学遥感数据中，Sentinel-2A数据是唯一在红边范围含3个波段(中心波长分别为560、665、705 nm)的数据，为细致地监测植被生长情况提供了更多可能性。

通过实地走访调查和物候期分析，选取2018年4个时期的Sentinel-2A卫星遥感影像(成像时间为5月15日、7月24日、8月23日和9月22日，产品级别为Level-1C)。数据来源于Copernicus Open Access Hub(https:∥scihub.copernicus.eu)。覆盖研究区的图像包括4景，编号分别为48SUG、48SVG、48SUF、48SVF。为保证数据质量，在进行作物精细分类前，利用Sen2cor软件包对Sentinel-2A影像进行大气校正；利用SNAP软件将校正结果重采样为ENVI格式，在ENVI软件中选取蓝绿红3个可见光波段、红边波段1和近红外波段(编号分别为2、3、4、5、8波段，其中第5波段空间分辨率为20 m，其他波段空间分辨率为10 m)进行波段叠加、影像拼接及裁剪，得到包含5个波段的研究区影像图。根据土地利用变更数据，掩膜获取耕地分布信息。

1.3 农作物训练样本及验证数据集

2018年7月底对研究区开展野外实地考察，并记录所到之处的地理位置和作物种类。利用实地调查数据，结合农事历，对影像进行目视解译，选取8 760个样本点，包括水稻、玉米、春小麦、春小麦+秋油葵、胡麻、洋葱、大棚作物、马铃薯和砂田瓜果9种作物类型。样本点分布较为均匀，其中2/3作为训练样本，1/3作为验证数据。

2 研究方法

2.1 研究思路

影像预处理之后，首先计算各个时期遥感影像的植被指数，构建特征参数数据集，分析不同作物不同指数的时序变化曲线；然后采用随机森林分类法进行农作物识别得到分类结果，最后根据验证样本进行精度评价。具体流程如图2所示。

图2 技术流程图Fig.2 Flowchart of technical route

2.2 NDVI和RENDVI计算及时序变化曲线提取

NDVI和RENDVI是农作物监测中常用的光谱指数，计算式为

(1)

(2)

式中NIR、R、VRE——Sentinel-2A数据的波段8、4、5的反射率

分别计算4个时期的两种指数图像，再进行组合。

由于不同指数特征或其组合对分类结果的精度影响及贡献程度不同[30]，本文设计了3种特征指数构建方案(表2)，共提取了5种特征指数作为分类特征：NDVI、RENDVI、NDVI+RENDVI、NDVI-RENDVI和NDVI&RENDVI(即NDVI和RENDVI通过layer stack组合)。

利用不同时期指数图，统计分析不同作物的特征指数时序变化曲线，反映不同作物类型的物候差异，分析不同作物类别的可分性以及时相选择的合理性。

2.3 多种特征的随机森林分类

随机森林分类是组合多棵决策树的预算结果，然后进行投票判断，预测准确率高，是多种遥感分类方法中精度较高的一种方法[31]。

本文将5种特征指数图像与训练样本一起输入到随机森林分类器，得到5种分类结果。

2.4 精度验证

本文利用混淆矩阵进行分类结果的精度验证。基于混淆矩阵，可以计算总体分类精度、Kappa系数、每一类别的制图精度和用户精度。

表2 特征指数构建方案Tab.2 Scheme of feature parameters

总体分类精度是指所有被正确分类的像元数量与总像元数量的比值，计算式为

(3)

Kappa系数计算式为

(4)

式中k——混淆矩阵行列数，代表分类的类别数量

xii——混淆矩阵中对角线上的值，代表正确分类的像元数

N——验证的像元总数

xi+——混淆矩阵第i行元素相加之和

x+i——混淆矩阵第i列元素相加之和

3 结果与分析

3.1 不同作物指数特征的时间变化分析

基于训练样本数据，统计各种农作物的时序植被指数特征(NDVI、RENDVI)均值，时序植被指数变化曲线如图3所示。

图3 不同农作物的特征指数时序变化曲线Fig.3 VI curves of different crops

从图3中可以看出，春小麦-秋油葵的特征曲线有2个峰，为一年两熟的种植模式，而其他作物的特征曲线都只有1个峰，为一年一熟的种植模式；在各种特征曲线上，洋葱、大棚作物、砂田瓜果3种作物都表现出较低值，并且生长期内起伏不大，尤其是砂田瓜果能与其他作物明显区分。

根据图3所示特征，9种作物大致可分为3种类型：高值型(水稻、玉米、胡麻、马铃薯)，曲线形态表现为钟型，5月特征值较低且比较集中，7月达到峰值，8月略有下降，9月特征值分散下降。低值型(洋葱、大棚作物、砂田瓜果)，曲线形态表现为比较平缓，全年特征值较低，7月有一小峰值。开口型(春小麦、春小麦-秋油葵)，曲线形态表现为开口型，即5月已有较高的特征值，9月仍处于较高值。

可以看出，5月是区分开口型和其他类型的最佳时期，7、8月是区分高值型和低值型的最佳时期。开口型春小麦、春小麦-秋油葵在5月中旬表现出较高的特征值，与其他类型作物的特征参数差异较大，可以明显区分；高值型、低值型作物在5月下旬的特征参数比较接近，难以区分。从农事历来看，大部分作物在5月处于苗期，植被指数值较低；到7月下旬，两种类型的差异明显分化，并保持到8月下旬；到9月下旬，大部分作物已成熟收割，高值型和低值型作物的植被指数值降低，两种类型的分异趋缓。

高值型的4种作物在5、7月影像上指数值都比较接近，难以区分，根据农事历，利用表现出成熟期时差的8、9月影像可以区分：马铃薯与胡麻成熟得较早，而水稻和玉米还处于生育期末期，因此马铃薯与胡麻的特征值低于水稻和玉米；低值型的砂田瓜果全年的特征值都低于其他作物，洋葱和大棚作物彼此之间的差异也是显著可见，洋葱的特征值全年高于大棚作物；开口型的春小麦和春小麦-秋油葵的特征参数到8月下旬秋油葵长势旺盛期可以明显区分。

各种作物的RENDVI曲线总体形态与NDVI曲线相似，但曲线起伏趋缓，高值型的作物在生长旺盛期(7月)的特征值分异比NDVI更加明显。

3.2 不同特征指数分类的精度对比

为了评估不同指数及其组合特征对分类精度的影响，本文采用随机森林算法对NDVI、RENDVI、NDVI+RENDVI、NDVI-RENDVI、NDVI&RENDVI 5种特征指数分别进行分类，并利用混淆矩阵进行分类精度验证，结果如表3所示。

5种分类结果中，采用NDVI&RENDVI特征组合分类的总体精度最高，其次是NDVI，最后是NDVI-RENDVI，NDVI&RENDVI分类精度比NDVI分类精度提高了约4个百分点，比NDVI-RENDVI分类精度高出约10个百分点。并且NDVI&RENDVI分类结果与验证数据取得了较好的一致性，Kappa系数为0.83，NDVI次之，Kappa系数为0.78，最低的是NDVI-RENDVI，Kappa系数为0.71。

表3 各种分类结果总体精度及Kappa系数Tab.3 Overall accuracy and Kappa coefficient of different classifications

NDVI-RENDVI分类结果精度最低。从算法原理上分析，NDVI和RENDVI两值相减将共同的红波段消除而用近红外与红边波段相减，削弱了植被由吸收谷红光波段到反射峰近红外的反射陡坎效应，导致NDVI-RENDVI分类的精度最低。

单独利用RENDVI分类的精度较低，但是将RENDVI特征与NDVI特征组合到一起的精度有明显提高。

从NDVI&RENDVI分类结果来看，春小麦、春小麦-秋油葵、砂田瓜果、马铃薯的分类效果最好(用户精度大于90%)，胡麻、洋葱、大棚作物分类效果相对较差(用户精度小于80%)，被误分为马铃薯、砂田瓜果的概率较高；在NDVI特征分类结果中，玉米、水稻、大棚作物3类的精度较低，而在NDVI&RENDVI特征组合分类结果中，这3类作物的用户精度有了显著提高，从而提高了总体的分类精度(提高了4.06个百分点)。

3.3 研究区作物种植结构与分布

景泰县2018年农作物NDVI&RENDVI分类结果如图4所示。

图4 2018年景泰县作物类型空间分布图Fig.4 Spatial distribution of crops in Jingtai County in 2018

由图4可以看出，具有干旱区特色的砂田瓜果(包括小金瓜、籽瓜等)播种面积最大，主要分布在海拔较高、气温偏低的正路镇、寺滩乡和喜泉镇；大棚种植方式在景泰县比较普遍，广泛分布于全县；玉米、洋葱、马铃薯、胡麻是景泰县重要的农作物，主要分布在北部和中部比较平坦的地区，红水镇、漫水滩镇、条山镇、草窝滩镇等地；水稻面积较小，主要分布在水源条件较好的五佛乡。 2018年景泰县作物种植结构如表4所示。