基于Landsat-8和Sentinel-2时间序列合成影像的山区甘蔗种植区提取

康倩文　徐伟恒　王雷光　洪泽湖　刘运

关键词：甘蔗；物候；时间序列影像；影像合成；谷歌地球引擎

中图分类号：S566.1 文献标识码：A

甘蔗是全球第一大糖料作物和第二大生物能源作物，中国是世界上最主要的甘蔗种植国家之一[1]。从甘蔗中提取的蔗糖是中国食品行业中不可替代的重要原料，而云南省作为我国第二大蔗糖生产基地[2]，蔗糖產业是云南省的高原特色产业，更是扶贫产业[3]。云南的甘蔗主要种植于丘陵山区旱坡地，准确识别甘蔗种植区分布范围，对有关部门制定相应的土地政策，为甘蔗产量和效益评估提供有效参考，对当地农业可持续发展、保障食糖供给安全都具有重要现实意义。

农作物的遥感识别方法主要包括两大类：单时相遥感影像识别和多时相遥感变化监测[4]。多时相遥感影像相较单时相影像更能反映植被的季节变化以及物候特征，有效增加农作物的识别精度，目前已经被广泛应用于农作物的识别与监测中[5]。如边增淦等[6]利用2018 年18 景GF-1 WFV影像构建出的归一化植被指数时间序列并采用分层策略逐步提取出黑河流域中游地区的作物种植结构。刘吉凯等[7]选用多时相GF_1 WFV 高分辨率遥感数据，采用多时相迭代法构建甘蔗提取决策树模型提取出广西江州区的甘蔗种植区，总体精度（overall accuracy， OA）达到95.26%。大量研究表明，联合多源遥感影像可以增加影像数量，增加对地观测有效像元次数[8]，有效弥补单种数据源数据量少的局限性[9]，是实现多云雾地区、大范围农作物种植信息提取的主要手段[10]。如谭深等[11]融合Landsat-8 与Sentinel-2 光学数据、Sentinel-1 微波SAR 数据和其他辅助数据，解决了南方地区多云天气导致影像质量不佳的问题，准确、高效地提取了海南省的水稻种植范围，OA为93.2% 。ZHENG 等[12] 基于Landsat-7/8 、Sentinel-1/2 合成的时间序列影像利用TWDTW法评估作物的物候相似性，并识别出中国的甘蔗种植区，SR 数据的OA 为93.47%。

不同作物具有不同的物候特征，植被物候能够反映植被生长规律，国内外学者通过提取植被物候特征，进行甘蔗等农作物的识别与分类。ALEXANDRE 等[13]利用多时相MODIS-EVI 数据通过甘蔗与其他作物的增强型植被指数物候特征的差异实现巴西圣保罗州农作物的分类。张东东等[14]通过选择特定时相的HJ-1 CCD 影像，利用甘蔗生长特点，采用决策树分类模型区分出中国南方地区的甘蔗与其他地物，OA 达到93.2%，Kappa 系数为0.81。WANG 等[15]基于Sentinel-1、Sentinel-2 和Landsat 影像构建的时间序列数据，采用基于像元和物候学的方法绘制了2018 年广西省甘蔗分布图， 甘蔗的制图精度（produceraccuracy， PA）和用户精度（user accuracy， UA）分别为88%和96%。

综上所述，目前学者们已经利用多时相、多数据源并结合物候特征的方法对甘蔗种植区进行了提取，但是缺乏专门针对山区甘蔗提取的研究。由于山区地物类型复杂，耕地碎片化现象严重，耕地地块细小狭窄且结构复杂[16]，导致耕地信息难以快速、精准获取。为明确山区甘蔗种植区分布的情况和规律，本研究综合利用多源、多时相遥感数据与植被物候的特点，以云南省玉溪市新平县为研究区，采用谷歌地球引擎（Google EarthEngine， GEE）云计算平台，基于Landsat-8 和Sentinel-2 长时间序列合成影像数据，依据甘蔗与不同地物之间的物候特征、地形特征，借助物候参数，探究适用于山区甘蔗种植区提取的方法。旨在开发一种适合提取山区甘蔗种植区的有效算法，为当地农业部门制定合理的土地利用政策提供参考。

1 材料与方法

1.1 研究区概况

新平县位于云南省玉溪市西部，地形复杂，多山地、丘陵，耕地空间破碎化，种植结构复杂，是云南省甘蔗主要生产县之一，也是云南元江干热河谷的典型山区农业县[17] 。新平县地处23°38′15″～24°26′05″N，101°16′30″～102°16′50″E，全县面积约为4223 km²（图1）。研究区地形以山地为主，地处哀牢山中段东麓，县境山区面积达4139.6 km²，占县域面积的98%；地势西北高、东南低，海拔为418～3143 m。新平县年均降雨800 mm 左右，降雨量少，蒸发量大，气候干热，全年无霜，有利于蔗糖的积累[18]，适宜栽种甘蔗。蔗糖产业是新平县的传统产业，是当地农民增收、企业增效的主要来源之一。

1.2 甘蔗提取技术路线

首先，将Landsat-8 OLI 和Sentinel-2 MSI 影像数据经过去云、光谱波段校准后计算遥感指数得到遥感指数时间序列，并将影像进行合成、插值、平滑，构建出以10 d 为间隔的高质量无云、无空洞的时间序列合成影像；其次，利用时间序列合成影像提取耕地图层，并基于耕地图层，利用甘蔗的物候特征、地形特征提取甘蔗，并进行精度验证（图2）。

1.3 数据源及预处理

1.3.1 遥感影像及预处理 本研究基于GEE 平台调用了2019 年10 月1 日到2021 年7 月1 日的覆盖研究区所有的Landsat-8 OLI 与Sentinel-2MSI 影像共计1022 张。所用Landsat-8 OLI 影像、Sentinel-2 MSI 影像分别来自美国地质调查局提供的Landsat-8 一级表面反射率数据集和Sentinel-2 的Level-2A 级地球表观反射率数据集（表1），2 种数据集均已经过几何校正、辐射定标、大气校正， 同时利用Landsat-8 OLI 的QA_PIXEL 波段和Sentinel-2 MSI 的QA 波段分别对数据集进行云和云阴影的掩膜。研究区内像元的总观测次数的空间分布见图3A，研究区内像元的良好观测次数即掩膜云和云阴影后的像元的有效观测次数的空间分布见图3B，统计结果显示，所有像元在2019 年10 月1 日到2021 年7 月1 日期间至少有35 次良好的观测，图3A 与图3B 中右下角观测次数较高的原因是该区域恰好位于Landsat-8 OLI 行列号为129/43、129/44、130/43的影像的重叠区域以及Sentinel-2 MSI 行列号为N0208_R104 、N0209_R104 、N0300_R104 、N0208_R061、N0209_R061、N0300_R061 的影像的重叠区域。

本研究采用ZHANG 等[19]提供的校正系数对Landsat-8 和Sentinel-2 影像进行光谱波段校准，以消除不同卫星传感器之间的差异。计算每张影像的光谱指数： 归一化植被指数（ normalizeddifference vegetation index， NDVI）[20]、陆表水分指數（land surface water index， LSWI）[21]、增强型植被指数（enhanced vegetation index， EVI）[22]和修正后归一化差异水体指数（modified normalizeddifference water index， MNDWI）[23]，计算公式如式（1）～（4）所示。

1.3.2 数字高程模型 数字高程模型（digitalelevation model， DEM）可生成坡度、坡向、剖面图等，辅助地形地貌分析。本研究所采用的DEM开源数据来源于美国国家航空航天局提供的第三级产品， 空间分辨率为30 m ， 产品ID 为“USGS/SRTMGL1_003”。由于新平县多山地的地形地貌，因此本研究考虑引入地形因子来辅助分析甘蔗的地形特征。

1.3.3 野外调查数据 为确定新平县甘蔗的物候特征，本团队于2022 年4 月进行实地调研。样本集由实地调查定位出在2020 年作物的种植类型未发生改变的样本和谷歌地球结合目视解译的解译样本组成，包括常绿植被、水体、不透水层、甘蔗与其他农作物。甘蔗与其他农作物的样本组成如表2 所示。

1.4 时间序列影像合成

影像合成压缩了数据的时间维度，减少了数据冗余，合成后的时间序列影像十分规则有序[24]，有利于研究地物的物候特征。最大值合成法可避免云雨天气对遥感影像计算的NDVI 值存在低值噪声的影响[25]，因此，本研究以10 d 为影像合成的时间窗口，以最大值合成法生成NDVI 时间序列，以均值合成法生成LSWI、EVI、MNDWI 时间序列，并将合成后的影像全部重采样为10 m。采用线性插值法[25]填补影像由于去云和去云影造成的空洞，利用移动窗口大小为9，滤波器阶数为2 的Savitzky-Golay 滤波器来平滑NDVI 时间序列，以消除数据噪声影响[26]。基于此，本研究构建了新平县2019 年10 月1 日至2021 年7 月1 日的时间序列合成影像数据集。

1.5 甘蔗种植区提取算法

实地调研显示，新平县的主要土地类型大致分为4 类：常绿植被、水体、不透水层和耕地。借鉴前人的理论基础通过遥感指数时间序列对4类地物进行提取，识别算法如表3 所示。其中，常绿植被在时间序列中被定义为全年存在超过90%的观测值的LSWI 大于0 且EVI 大于0.2；水体在时间序列中被定义为全年存在至少5%的观测值的MNDWI 大于NDVI 且NDVI 小于0.1；不透水层在时间序列中被定义为全年超过90%的观测值的LSWI 小于0.2。耕地图层通过掩膜常绿植被、水体、不透水层的图层范围得到，并通过在耕地图层中研究甘蔗与其他作物的物候差异来实现对甘蔗种植区的提取。

植物物候现象是环境条件季节和年际变化最直观、最敏感的生物指示器[30]。植被指数时间序列能够较好地反映植被的代谢强度及其季节性变化及年际变化。通过研究分析甘蔗与研究区内主要农作物（水稻、玉米、香蕉）的NDVI 时间序列曲线（图4），发现新平县甘蔗大多在3 月末至4 月初播种，此时甘蔗的NDVI 值最低，9 月甘蔗的NDVI 值上升到0.7 以上，11—12 月甘蔗成熟，12 月至次年4 月为甘蔗收割期，甘蔗的NDVI 值开始下降；研究区内水稻大多为双季稻，在1 年内有2 个生长周期，与甘蔗明显不同；玉米比甘蔗的种植时间晚而收获时间早；香蕉与甘蔗的NDVI 曲线类似，但是甘蔗收割后表现出与裸土相似的NDVI 值，与香蕉明显不同。甘蔗的生长周期较长，一般在8 个月至1 年以上[31]，并且甘蔗是一种高生物量作物，因此甘蔗有较长一段时间的NDVI 值处于较高水平。

为了定量研究甘蔗与其他农作物的物候特征差异，我们引入物候参数来说明。植被指数时间序列被广泛应用于植被物候参数提取[32]。利用物候参数可以较好地进行农作物分类研究[33]。本研究根据NDVI 时间序列曲线提取的物候参数如图5 所示，物候参数的说明如表4 所示，其中基值（base of season， BOS）为生长季阶段NDVI 曲线左边最小值（a 点）与右边最小值（f 点）的均值，振幅（amplitude， AMP）为生长季阶段NDVI 曲线最大值与基值的差值，上升时间（rise time， Tr）表示生长季阶段NDVI 从振幅的10%（b 点）上升到振幅的90%（c 点）所用的时间，下降时间（fall time， Tf）表示生长季阶段NDVI 从振幅的90%（d 点）下降到振幅的10%（e 点）所用的时间，上积分（above integral of season， AIOS）表示生长季阶段NDVI 曲线与基值之间的积分，下积分（below integral of season， BIOS）表示生长季阶段基值与坐标0 轴之间的积分[34]。本研究选择物候参数上升时间、下降时间、上积分、下积分作为甘蔗提取的物候特征，并绘制甘蔗与其他农作物的训练样本在物候参数上升时间、下降时间、上积分、下积分上的直方图分布（图6）。本研究基于90%的置信区间获取甘蔗训练样本的物候参数范围，并依据此范围进行研究区内甘蔗种植区的提取。物候参数阈值分布范围为：上升时间大于30 d 且小于190 d，下降时间大于70 d 且小于200 d，上积分大于90 且小于164，下积分大于26 且小于82。

由于新平县地形以山地为主，因此本研究考虑引入海拔、坡度因子作为山区甘蔗提取的辅助特征。据研究报道，云南省甘蔗海拔分布范围为422～1300 m[35-37]，坡度达到25°以上。实地调研显示， 新平县甘蔗种植区大致分布在海拔480～1330 m，坡度30°以下。因此综合考虑研究报道以及实地调研结果调整研究区内提取甘蔗种植区的海拔及坡度因子范围为海拔为400～1400 m，坡度在30°以下。

综上，本研究提出的甘蔗种植区提取算法为：30

2 结果与分析

2.1 主要土地类型提取结果

本研究提取的新平县常绿植被、水体、不透水层和耕地如图7A 所示。图7B 为图7A 中标记为1 区域的细节放大图，绿色区域为新平县常绿植被提取结果的局部示意图；图7C 为图7A 中标记为2 区域的细节放大图，蓝色区域为新平县水体提取结果的局部示意图；图7D 为图7A 中标记为3 区域的细节放大图，粉色区域为新平县不透水层提取结果的局部示意图；图7E 为图7A 中标记为4 区域的细节放大图，黄色区域为新平县耕地提取结果的局部示意图。由图7A 可知，新平县耕地主要集中在新平县中部及北部。

2.2 甘蔗提取结果及精度评价

为了验证本研究算法对研究区甘蔗的提取精度，依据实地调研情况在Google Earth 中随机选取274 个甘蔗样本与470 个其他农作物样本作为验证样本，基于混淆矩阵构建OA、Kappa 系数、PA 和UA 进行山区甘蔗提取精度评估。新平县甘蔗提取的OA 达97.07%，Kappa 系数为0.83，从分类精度评价指标来看，提取OA 较高，整体分类精度满足空间分析与实际应用需求。甘蔗提取的PA 与UA 分别达到88.86%和80.57%，其他农作物PA 与UA 分别达到97.88%和98.89%，提取面积为7705 hm²，提取结果较好。

本研究提取的2020 年新平县的甘蔗种植区分布图见图8A，其中图8B 和图8C 分别为图8A中标记为1、2 区域的细节放大图（透明度60%的浅绿色区域为识别的甘蔗）。由图8A 可知，2020年新平县甘蔗种植区大多分布于新平县中南部及北部，集中于漠沙镇、戛洒镇、老厂乡，其他地区甘蔗种植区零星分布。从甘蔗种植区分布地形地貌来看，2020 年新平县甘蔗种植区多分布于新平县北部及中部偏南，分布于地势较平坦、海拔相对较低的地带，由于新平县地势西北高、东南低，甘蔗种植区也呈现东南多于西北的现象。进一步分析研究区内甘蔗种植区的空间分布情况可知，新平县甘蔗地块分散且不规则，地块边界模糊。该分布图可以作为开展甘蔗种植区长势、灾害等遥感监测的基础数据。从乡镇尺度来看，新平县甘蔗分布面积较多的区域是漠沙镇、戛洒镇、老厂乡，其中漠沙镇甘蔗分布最多（图9）。漠沙镇、戛洒镇、老厂乡的甘蔗种植面积分别为2786、1458、1332 hm²，分别占新平县甘蔗种植面积的36.16%、18.92%、17.29%。新平县甘蔗分布面积较少的区域分别是古城街道、桂山街道、建兴乡，古城街道甘蔗分布最少，古城街道、桂山街道、建兴乡的甘蔗种植面积分别为0.87、54.00 、13.70 hm²，分别占新平县甘蔗种植面积的0.01%、0.70%、0.18%。分析发现新平县各乡镇甘蔗提取分布趋势情况与实地调研结果一致。

2.3 物候参数敏感性分析

对本研究中用到的上升时间、下降时间、上积分、下积分4 个物候参数进行敏感性分析（图10）。可以发现，选取上升时间、下降时间、上积分、下积分这4 个物候参数作为提取山区甘蔗种植区的物候特征时，甘蔗提取算法中甘蔗的UA为88.85%达到最高，PA 为80.57%，此时甘蔗的错分误差最小，漏分误差最大。当甘蔗提取算法中物候参数分别设为下降时间、上积分、下积分（Tf+AIOS+BIOS），上升时间、上积分、下积分（Tr+AIOS+BIOS），上升时间、下降时间、上积分（Tr+Tf+AIOS），上升时间、下降时间、下积分（Tr+Tf+BIOS）这4 种情况时，甘蔗提取算法中甘蔗的UA 分别下降为88.50%、79.37%、75.35%、67.20%，PA 分别上升为82.42%、84.40%、83.52%、85.00%。物候参数敏感性分析表明，选取上升时间、下降时间、上积分、下积分4 个物候参数参与甘蔗提取算法时甘蔗提取效果最好，4个物候参数在提取甘蔗时都是重要的物候参数，都是甘蔗制图算法的重要组成部分。本研究甘蔗制图算法中所用物候参数对于耕地中甘蔗种植区的提取有较好的提取能力，有助于甘蔗与其他农作物的物候特征分析。

3 讨论

本研究提出了适用于山区复杂地形，基于Landsat-8 OLI 与Sentinel-2 MSI 合成的时间序列遥感影像数据的甘蔗种植区的提取方法。以典型山区甘蔗种植区新平县为研究区，生成一幅10 m空间分辨率的2020 年新平县甘蔗种植区分布图。本研究中甘蔗种植区提取的OA 为97.07%，Kappa系数为0.83， 提取精度较为良好， 表明基于Landsat-8 OLI 与Sentinel-2 MSI 合成的时间序列遥感影像数据结合光谱指数特征、物候特征与地形特征可以实现山区甘蔗种植区的提取并绘制出甘蔗分布图，精度满足基本需求。提取的新平县甘蔗种植区面积为7705 hm²，其中漠沙镇、戛洒镇、老厂乡的甘蔗种植最多，总共占新平县甘蔗种植区的72.37%，古城街道、桂山街道和建兴乡的甘蔗种植最少，总共占新平县甘蔗种植区的0.89%，可见乡镇间甘蔗种植区分布差异显著。本研究提出的甘蔗種植区提取算法的UA和PA 不是很高，其原因是新平县耕地种植结构复杂，地块破碎，中等空间分辨率的影像像元内可能包含多种地类；本研究选择可免费获取的Landsat-8 和Sentinel-2 影像在多云天气的新平县很难保证每个像元在10 d 内有一幅高质量的卫星影像，经过插值后的像元值在反映作物生长的真实情况时存在一定的不确定性，因此未来可利用不受天气影响的Sentinel-1 A/B 影像来辅助进行甘蔗种植区的提取。

新平县是典型的山区甘蔗种植区，在利用物候特征提取山区甘蔗时，多源、多时相遥感数据的融合可极大提高遥感影像的空间与时间分辨率，尽可能避免山区多云雨、复杂的地形地貌因素导致山区遥感影像质量不高[38]的问题，为山区甘蔗提取与精细化制图提供一种有效参考。本研究基于NDVI 时间序列数据，借助甘蔗与常绿植被、水体、不透水层、其他农作物在光谱指数特征、物候特征、地形特征上的差异，采用上升时间、下降时间、上积分、下积分4 个物候参数以及海拔、坡度因子确定提取甘蔗的最佳阈值，相较张东东等[14]学者选择用多时相影像，单一光谱特征，利用纯净甘蔗训练样本进行决策树分类的方法来进行甘蔗种植区的提取，本研究不仅增加了研究区内像元良好观测次数，克服了山区等多云雨地区影像质量不高的问题，并且可以更充分地利用甘蔗的生长特征及物候特征，对甘蔗样本长势情况的依赖相对较小。本研究利用4 个物候参数确定了甘蔗提取的最佳物候参数阈值，上升时间、下降时间、上积分、下积分4 个物候参数在提取甘蔗时都是重要的物候参数，为未来其他地区的山区甘蔗种植区提取提供重要的物候特征参考。

研究区内地貌复杂，耕地破碎度较高，由于研究区海拔及气候等因素的影响，耕地在休耕期会出现杂草丛生，作物收割残留枝叶等现象，在影像上导致光谱混淆，容易出现“异物同谱”现象。因此，使得在甘蔗种植区提取过程中难免出现错提或漏提现象。今后考虑将使用更高分辨率的遥感数据或运用更小的时间分辨率组成的时间序列数据来进一步提高提取的精度，建立普适性更高的山区甘蔗种植区提取算法。本研究的提取精度虽然满足县域的山区甘蔗提取精度，适用于作为当地农业部门的决策建议，但是甘蔗提取的用户精度和制图精度还是有待提高。希望之后考虑通过引入甘蔗的纹理特征等，提高山区甘蔗的提取精度，实现山区甘蔗的精细化提取。