哨兵2号多时相植被指数作物分类及监测

祁向前,孙德浩,贾连星

(1.龙岩学院 资源工程学院，福建 龙岩 364012;2.中国矿业大学 环境与测绘学院,江苏 徐州 221116;3. 黑龙江科技大学 矿业工程学院,哈尔滨 150022)

以往对农作物的监测(分类、面积变化等)大多靠人工现场采集信息，既耗时又费力，近年来随着现代化脚步的加快，多种类的卫星可以满足不同的监测需求，分辨率的提高也保障了监测信息的准确性，卫星用于监测农业生产，能够体现其快速、客观、准确的特点，有利于及时掌握农情基本信息，进行宏观调控。国内外学者研究表明，利用多时相遥感影像植被指数和不同作物的物候特征结合进行作物识别与监测是一种行之有效的方法。例如Jakubauskas等[1]采用谐波算法对NOAA/AVHRR-NDVI时序数据进行研究，得到作物分类信息。Moody[2]等利用归一化差值植被指数(Normailized Difference Vegetation Index,NDVI )时间序列的Fourier变换来划分土地的覆盖类型等。Lanjeri[3]等利用多时相掩膜法成功提取了葡萄园的种植面积。竞霞等[4]利用多时相Landsat TM数据，提取冬小麦不同时相的作物波谱信息，成功监测了北京地区的冬小麦种植结构调整情况。Young和[5]Running[6]通过分析研究区植被类型的年NDVI变化情况，并结合作物的物候特征对植被类型进行分类。Conese等[7]对遥感数据进行了单时相和多时相的研究，发现多时相数据能够极大改善分类精度。贾树海等[8]针对3个不同时期的遥感影像，结合不同作物的物候特征提取其NDVI植被指数，采用监督分类的方法，较好地提取了花生的种植面积并制图。

哨兵2号是由两颗极轨卫星组成的星座，具有高分辨率多光谱成像的特点，一颗卫星的重访周期为10 d，两颗互补，重访周期为5 d，在红边范围内有较多波段，有利于对植被的监测。国内外均有学者利用单时相哨兵2号数据实现了作物识别[9-10]。NDVI作为一种植被指数，被广泛应用于反映作物在不同季节、不同生育期表现出来的生理特征，也有效的应用于土地利用覆盖监测、植被覆盖密度评价、作物识别和作物产量预报等[11-13]。哨兵数据的波段最高分辨率为10 m，单纯依靠目视解译对作物进行分类较为困难，因此构建多时相NDVI植被指数来研究就能够很好地提高分类精度[14-16]。

上述文献主要是针对不同种类农作物进行分类或者是对单一农作物进行面积提取，涉及到的农作物精细分类(例如夏玉米、春玉米、山地夏玉米等)及多种农作物种植面积变化研究还很少。文中以河北省石家庄市灵寿县的多种作物及地貌类型为例，采用多时相决策树分类法，利用所建立的作物及地貌样本点不同时期内监测到的NDVI值建立时序NDVI曲线，并在具有明显区分特征的月份确定阈值，对研究区域内的多种作物及地貌类型进行分类研究。

1 研究区与数据处理

1.1 研究区

灵寿县(113°45′E～114°28′E，38°16′N～38°48′N)位于河北省西部，太行山东侧，在省会石家庄西北35 km处，东南与正定县毗连，南与鹿泉区隔滹沱河相望，西同平山县接壤，西北隅与平山县、阜平县及山西省五台县接境(图1)。全境呈西北—东南走向，宽15～20 km，长约100 km，灵寿县属北温带亚湿润气候，处于半干旱、半湿润季风区。年平均气温13.7 ℃，无霜期190 d左右，年平均降雨量417 mm。全县总面积1 066.2 km2，耕地面积521.12 km2，约占全县总面积的48.87%，主要分布在中南部区域，双熟制耕地主要有冬小麦-花生和冬小麦-玉米大豆，单熟制耕地主要有春玉米、夏玉米、山地夏玉米，县域的北边多为山区林地，其他的植被及地貌类型有杂草、河流湖泊和居民地道路裸地等。土壤有黏土、沙壤等多种类型土质，河流主要有滹沱河、慈河、松阳河、卫水河等。

图1 2020年研究区域

1.2 数据获取与预处理

根据作物的物候特征了解到作物有明显植被指数区分度的月份为4月、6月、8月和10月，后续用所监测样本点的时序植被指数变化分析也验证了该想法，并且考虑到云量会对影像分类精度造成影响，因此通过综合比较各期影像云量覆盖率，分别选取2019—2020年两年的8期影像，每年4期，分别是4月、6月、8月和10月(表1)，涵盖了全部农作物的主要生长季(其中4月、6月、10月影像无云量覆盖，8月影像云量覆盖度为0.8%)。8期影像数据均来自欧洲航天局(ESA)(https://scihub.copernicus.eu/dhus/#/home)，哨兵2号由哨兵2A和哨兵2B两颗卫星组成，重访周期为5 d，共由13个波段组成，其中较多的红外波段有利于植被的监测，分辨率有10 m、20 m、60 m，经过处理后可达到10 m。网站的影像数据出现离线状态时将其收藏并发送下载请求，过半个至一个工作日即可下载。

表1 影像清单及各农作物发育期

从欧空局下载下来的数据是L1C级，无需进行辐射定标，经过大气校正后进行重采样(10 m)，进行波段融合，剔除掉不需要的海岸/气溶胶波段(band.1)和水蒸气波段(band.9)。经过裁剪和空间配准，得到处理后的影像数据。

图2 技术流程图

2 方法研究

2.1 技术路线

研究技术路线如图2所示，在2019年3月通过外业GPS定位点来建立样本点监测NDVI值的变化(图3)，从北至南分别为杂草、林地、山地-夏玉米、河流、湖泊、春玉米、夏玉米、冬小麦-花生、居民地、道路、裸地，冬小麦-玉米大豆，监测时间持续到10月底所有作物收割结束(或正在收割)的NDVI值的变化情况，绘制NDVI时间序列曲线。了解农作物的物候特征，选取具有代表性的8期多时相遥感数据，预处理后计算每幅影像的归一化植被指数NDVI，NDVI算式如式(1)。根据农作物多时相NDVI变化情况建立分类判别规则，确定相应的阈值将各种作物及其他植被和地貌进行分类并进行分类统计，得到两年以来各种作物的种植面积的变化情况。为了验证分类结果的准确性，于2019年8月及2020年8月到现场勘察，采集多个样本点建立验证样本集，对当年的作物及其他地貌分类通过混淆矩阵进行精度评价来验证分类的准确性。

(1)

式中:BNIR为哨兵2号近红外波段(第八波段)的地表反射率;BRED为哨兵2号红光波段(第四波段)的地表反射率。

图3 样本点位置

2.2 植被指数曲线及可分性

本研究区域内典型的地表类型有双熟制耕地(冬小麦-花生、冬小麦-玉米大豆)，单熟制耕地(春玉米、夏玉米、山地-夏玉米)，林草类(林地、杂草)，非植被地表(居民地道路裸地、河流湖泊)。其中非植被地表的居民地道路裸地的植被指数常年较低，在0～0.1附近，河流湖泊的植被指数全年低于0，在-0.1～0之间。双熟制耕地NDVI曲线具有双峰，根据4月影像可以与单熟制作物区分开，10月底具有较低植被指数可与林地做区分，花生收割时间相比大豆较早，10月底花生收割完(植被指数低于0.20)，大豆正在收割(植被指数大于0.20)。单熟制耕地中，春玉米4月份播种，6月与其它两种玉米相比有较高的植被指数，夏玉米和山地-夏玉米6月播种，此时植被指数很低，三者在8月植被指数到达峰值，其中山地-夏玉米由于多种植在山地丘陵地区，自然条件恶劣且密度低，植被指数不超过0.6。林草类植被由于不是人工种植，在8月后植被指数并不会出现大降幅，其中林地常年较高，杂草多出现在道路两旁及山区，植被指数在8月到达峰值时低于作物的峰值(图4)。

图4 农作物及林草类NDVI时序曲线

在构建的NDVI时序曲线中，非植被地表由于全年低植被指数因此有着极高的识别度；双熟制耕地作物及林地特征差异明显，识别度较高；3种玉米具有相同的植被指数变化趋势，在起伏及峰值上有差异，也较容易识别；杂草种类多样，密度低，分布不集中，不能获得全部种类杂草的植被指数变化情况，只能得到一个大致变化范围，识别度一般。综上所述，利用NDVI植被指数来对作物及其他地表类型进行识别分类具有较好的可行性。

2.3 决策树方法用于作物及地貌识别

决策树作为一种分类方法，具有灵活、直观、运算效率高的优点。它是一种从无次序、无规则的样本数据集中推理出决策树表示形式的分类规则方法，采用自顶向下的递归方式，在决策树的内部节点进行属性值的比较并根据不同的属性值判断从该节点向下的分支，在决策树的叶节点得到结论。该方法在遥感影像分类和信息提取中已得到广泛应用。文中根据上述分析，设定不同的阈值来构建决策树对作物及地貌进行分类研究。

阈值的设置必须满足具有明显区分度的要求，才能保证较高精度的分类。分别将4月、6月、8月、10月影像的NDVI值设为B1、B2、B3、B4。根据参考时序曲线以及多次改进阈值可以得到，当B2大于0.6且B4大于0.5时可以提取林地区域。其余地貌类型都满足B2小于0.6，在此基础上，当B1大于0.35时可将冬小麦-花生和冬小麦-玉米大豆和其他作物和地貌类型做区分，根据B4大于0.20可以提取冬小麦-玉米大豆，B4小于0.20可以提取冬小麦-花生。其余作物和地貌类型满足B1小于0.35，在此基础上当B3大于0.65时可将夏玉米和春玉米提取出来，此时当B2大于0.45可提取春玉米，B2小于0.45可提取夏玉米。剩下两种植被类型山地-夏玉米和杂草都满足B3大于0.41且小于0.65，再根据B4大于0.20提取杂草，B4小于0.20提取山地-夏玉米。最后的地貌类型为非植被地表，在上述基础上，B1大于0可提取居民地道路裸地，B1小于0可提取河流湖泊。

3 结果分析与精度验证

3.1 分类结果及分析

根据多时相遥感数据获取的NDVI确定的最佳分类阈值，采用决策树分类方法得到2019年和2020年9种地貌类型的分类结果如表2所示。

表2 2019年和2020年各地表类型分类占比 %

决策树分类结果如图5、图6所示。

图5 2019年决策树分类结果

图6 2020年决策树分类结果

根据分类结果可知，与2019年相比，2020年林地面积增加了1.6%，造成略微增加的原因是2019年8月影像西北地区存在部分云，导致西北部分林地未分类(黑色区域)；河流湖泊下降0.2%，造成差异的原因是部分河流湖泊表面漂浮植被和藻类，使得少部分错分为植被；冬小麦-花生增加了4.2%，冬小麦-玉米大豆减少了3.8%，第一年双熟制耕地大部分种植冬小麦-玉米大豆导致土壤中的磷元素过度消耗，为了避免连作对土壤造成伤害(土壤板结、土壤团粒结构破坏、土壤疏松度下降)，因此第二年双熟制耕地大部分改种植冬小麦-花生；同理，春玉米增加了17.4%，夏玉米减少了19.6%也是为了降低对土壤的破坏；山地-夏玉米由于多种植在丘陵山区且密度较低，对土壤的伤害可忽略不计，所以种植面积基本无变化(增加了0.06%)；居民地道路裸地增加了2.8%，主要是城镇化水平的提高导致居民地面积的增加；杂草减少了3.9%，杂草种类多密度低分散广较难区分，每年的面积大致维持在10%左右。

3.2 精度验证

图7 验证样本点位置

分别在2019年8月和2020年8月(植被最茂盛期)进行实地调查，再结合哨兵影像目视判读等方式解译出地表类型，以一个像元为基本单位建立验证样本集，由于各地表类型所占面积差异过大，因此针对占比大的地貌类型选取较多验证样本点，占比少的地貌类型选取较少验证样本点(图7)，以kappa系数、制图精度和用户精度对两年的分类结果进行验证。两年的kappa系数分别为0.885 8和0.910 0，总体分类精度分别为90.76%和92.54%，9种植被类型的制图精度和用户精度如表3所示。根据kappa系数的判别原则：0～0.40(不含0.40)分类一致性差，0.40～0.80(不含0.80)分类中度一致，0.80～1分类高度一致，可知文中9种植被类型的遥感分类总体在精度上达到了高度一致。其中，各地貌类型样本点被正确分类个数及占比如表4所示。

通过精度验证可知，林地、居民地道路裸地、河流湖泊面积大、地块连片分布，有着明显且易区分的NDVI值，面积估算精度很高，大致在95%左右；冬小麦-花生和冬小麦-玉米大豆和其他作物相比有着独特的双峰特性，面积估算精度较高，大致在90%～95%之间；春玉米、夏玉米、山地-夏玉米单熟制作物种植较分散，NDVI变化曲线一致，根据峰值特性进行区分，面积估算精度一般，大致在80%左右；杂草由于其野生性，分布极其分散，很多杂草区域长宽小于10 m，小于哨兵2号最小像元面积，会与周围地物形成混合像元，因此面积估算精度较差，大致在70%左右。

表3 2019年和2020年9种植被类型分类的制图精度和用户精度 %

表4 2019年和2020年各地貌类型验证样本点个数及正确分类个数

4 结束语

文中研究了基于2019年所监测的多时相NDVI曲线来构建阈值，根据2019年和2020年4月、6月、8月、10月共8期影像，利用决策树方法将耕地类、林草类以及非植被地表类共9种植被类型进行了分类，最后根据验证样本集利用混淆矩阵对分类结果进行精度验证，结果表明：

1)该方法较好地提取了河北省灵寿县2019年和2020年两年9种不同的植被信息，其中包括3种不同品种的玉米，该技术方法不但可以对不同作物及植被类型进行分类，而且可以较好地进行精细分类(3种玉米的分类)，同时监测作物的时空变化规律，可以对来年的作物种植趋势进行预测，在一定程度上为政府部门的行政决策提供依据。

2)采用哨兵2号多波段构建NDVI值，并结合植被的多时相NDVI构建决策树进行分类，与传统的监督分类和非监督分类相比，该方法更加注重挖掘和分析不同植被的光谱信息差异，增加了植被类型可分性的信息量，分类识别的针对性更强，显著提升了分类精度。

3)基于哨兵2号时序植被指数的作物精细分类及时空监测上还有很多方面需要改进。8月份的影像大多含有云层，云层的存在极大降低了分类精度；部分河流湖泊边缘表面漂浮着植被，会将该部分水体错分为植被；玉米杂草分布不规则，混合像元较多，对面积提取的精度造成一定的影响。

针对上述问题今后可以将研究工作转向哨兵1号影像的分类，利用哨兵1号雷达数据具有穿透性的特点，可以避免云层的影响，同时构建雷达反射系数对植被类型进行分类，将两种分类结果进行比较以提高整体分类精度。其次通过上述研究可知作物的种植成数(即面积比例)和地貌的形状指数对分类结果影响较大，作物所种植的区域越集中越齐整，面积提取的精度也就越高，因此也可以将研究重心放在农作物的地块特征上。