基于Sentinel-2 时间序列数据及物候特征的棉花种植区提取

美合日阿依·莫一丁， 买买提·沙吾提,2,3, 李金朝

（1.新疆大学地理与遥感科学学院，新疆 乌鲁木齐 830046；2.新疆绿洲生态重点实验室，新疆 乌鲁木齐840046；3.智慧城市与环境建模自治区普通高校重点实验室，新疆 乌鲁木齐 840046）

棉花作为重要的经济作物，在新疆地区广泛种植。2014年以来，我国在新疆实施棉花目标价格改革，对疆内各地区棉花播种面积、结构、分布产生较大影响［1］。由于新疆环塔里木林果业的迅速发展，果树和棉花等作物种植模式越来越复杂，传统的监测方法无法满足大尺度、复杂地理环境下快速、高精度的棉花识别与面积监测的迫切需求［2］。因此，继续开展棉花高精度识别和信息提取等研究仍有重要意义。

植被物候是指植被发芽、叶片蔓延、抽穗、成熟和脱落的周期性变化，能够反映生长规律，不同作物具有不同的物候特征［3］。因此，可以将这些物候特征进行量化分类。遥感技术的发展为大面积提取农作物物候特征提供了有效的技术途径。国内外很多学者利用多时相Landsat TM［4］、MODIS［5］等中、低分辨率数据的归一化植被指数（Normalize difference vegetation index，NDVI）或增强植被指数（Enhanced vegetation index，EVI）获取作物的生长发育期变化，提取物候信息，进行农作物识别与分类［6-7］。这些影像数据具有高时间分辨率、获取方便与低成本等优势应用于地表植被的生长监测和识别，对种植结构单一的作物具有很大的优势［8-9］。但是这些影像数据混合像元量众多，光谱辨识率有限，而且使用的特征变量往往比较单一，因此，在地块破碎、作物生长周期相近、光谱差异小的情况下很难对其进行有效的区分，无法满足大范围、复杂地形条件下的棉花空间分布监测精度要求［10］。

欧洲航天局Sentinel-2卫星，共有13个波段，空间分辨可达10 m，比MODIS、Landsat-8相比在705～865 nm 特有3 个红边波段［11］，为红边波段作物遥感监测提供了数据保障，解决了棉花种植面积监测中精度与成本相互制约的问题。国内外已有很多学者通过Sentinel-2数据探讨了农业监测［12］、土地利用分类［13］、作物分类及面积提取［14］中的应用并证实了红边参数波段对监测植被健康和长势的有效性［15］。

渭干河—库车河三角洲绿洲（简称渭—库绿洲）是塔里木河流域开垦较早并在西北干旱区具有代表性的绿洲。以灌溉农业为主，耕地不断向绿洲外围扩张，种植地块较为细碎。棉花、玉米、果树等生长周期相似的秋季作物套种模式复杂。因此需要开展基于高分辨率数据的棉花识别特征和分类方法的探索。因此，本研究利用36 期多时相Sentinel-2数据，在已有NDVI 时序数据的遥感监测基础上结合改进的红边归一化植被指数（Red edge normalize difference vegetation index，RENDVI783）时序数据，利用Savitzky-Golay（S-G）滤波算法对其进行平滑和重构，提取11个关键性物候特征并其进行特征重要性评价，构建不同的分类特征数据集，用随机森林分类（RFC）方法进行分类和棉花种植面积提取。为农业部门对棉花种植面积监测提供重要的参考依据。

1 研究区概况

渭—库绿洲（82°20′～83°20′E，41°00′～41°40′N）地处天山南麓，塔里木盆地中北部，属于新疆阿克苏地区，包括库车市、沙雅县和新和县［16］。研究区平均海拔971 m，地势北高南低，自西北向东南倾斜，是干旱与极端干旱区典型的扇形平原绿洲。年平均气温10.3 ℃，年降水量50.0～66.5 mm，年均蒸发量2000～2092 mm。作物类型有棉花、玉米、小麦、核桃、红枣等，其中棉花、果树、玉米等秋收作物生长期为4—10 月。据统计，渭—库绿洲的棉花面积分别占全疆和阿克苏地区棉花面积的8.56%和38.20%，产量分别占8.41%及40.34%，是新疆棉花主产区［17］。研究区及样点分布如图1所示。

图1 研究区和样点分布Fig.1 Schematic diagram of the study area and sample distribution

2 数据与方法

2.1 数据来源与处理

野外实地考察是遥感监测研究的重要组成部分，基于遥感影像分析和解译验证依赖于实地实测数据。实验团队分别于2018 年7 月、2019 年2 月、2021 年5 月进行了野外实地考察，获取了778 个典型地物的样点，其中有532 个训练样本和246 个验证样本，并利用GPS 采集每一个样点的经纬度、周围地物分布情况并记录了其他属性，目的是为了考察影像分类所需要的各类地物的样本和不同作物的分布特征。据实地考察和对研究区的概况分类类型最终确定为棉花、玉米、果树、居民地、水体和其他6种地物。

因此，利用2018—2020 年作物连续生长季（1景/1 月）的36 景Sentinel-2 数据来构建植被指数时序数据，得到每个像元每个时相的植被指数变化信息（表1）。遥感影像云量均在5%以下，部分影像云量超过10%，但这些云主要分布在北部山区，并不影响作物分类。利用SNAP 和ENVI 等影像处理软件进行预处理，包括大气校正、几何校正、重采样、裁剪等［18］。

表1 遥感影像基本参数Tab.1 Basic parameters of remote sensing data

2.2 研究内容

对Sentinel-2 时间序列影像进行预处理，构建NDVI 和RENDVI783时序数据。通过S-G 滤波对NDVI、RENDVI783时序数据进行平滑、重构并提取11个物候特征，利用袋外误差法（OOB）对物候特征进行特征重要性评价，选取特征重要性得分较高的特征，构建物候特征优选组合。在此基础上利用重构后的时序数据（NDVI Fit 和RENDVI783Fit）、物候特征（RENDVI783Ph）、物候特征优选组合等数据集作为分类特征数据集，构建6 种不同的分类特征数据集并采用RFC 方法进行作物分类及提取。为了验证RFC 方法的分类精度，采用最大似然分类（MLC）方法与支持向量机（SVM）分类方法分别对上述6种分类特征数据集进行作物分类及提取。根据分类结果最佳特征组合的分类影像提取研究区棉花及其他作物的面积。

2.3 研究方法

2.3.1 植被指数时序数据构建为了探究棉花在生长过程中的变化，构建两种监测作物生育期和长势规律的指数（NDVI和RENDVI）。将Sentinel-2卫星红波段B4（665 nm）、红边1波段B5（705 nm）、红边2波段B6（740 nm）、红边3 波段B7（783 nm）、近红外波段B8（843 nm），按表2公式［19］计算。

表2 植被指数及计算公式Tab.2 Indices and calculating formulas of vegetation

2.3.2 植被指数时序重构及物候特征参数提取本研究对时序植被指数进行拟合、平滑、重构，并提取11 个物候特征参数，其含义［20］如表3 所示。S-G 拟合法是由Savitzky 和Golay 提出的基于最小二乘法的卷积算法，计算公式如下：

式中：Yj*为平滑后NDVI 值；j为原始NDVI 时间序列的序数；i为任意一个NDVI 时序数据；Yj+1为平滑后的NDVI 的序数；Ci为从滤波窗口首部开始第i个NDVI值的权值；N为滤波窗口大小；m为N/2滤波窗口长度［21］。

2.3.3 特征优选方法RFC方法不仅可以对遥感影像进行分类，在特征集的选择与优化中发挥着重要作用。从训练样本中随机抽取70%的样本进行分类，假设最初的训练样本集总有N个样本，每个样本都有M个特征。当样本数量较大时，约有30%的样本未被抽取，这一部分被称为OOB，可以估计分类精度，也可以计算参与分类的不同特征的重要性，特征变量重要性评估模型如下［22］：

2.3.4 分类方法与分类方案MLC是根据贝叶斯判决准则，在多类判决中用统计的方法建立判别函数集，即将遥感影像多波段数据的分布作为多维正态分布来构造判别分类函数［23］。

SVM 基于通过最小化经验风险和置信区间来最小化结构风险的统计理论。将神经网络结构选择问题转化为比较容易的核函数选择问题，有很好的优势［24］。

式中：f(x)为目标函数；x为向量之间的内积；w为权向量；b为分类阈值，决定超平面距离原点的距离。

RFC 方法的主要思路是基于遥感数据及特征变量基础上，首先对数据进行随机重采样，然后构建多个决策树，最后采用多决策树投票的方式确定数据类别归属［25］。

式中：Gini(A)为A集合的基尼指数；B为训练样本中样本种类数；b为随机选中的样本类别；pb为集合A中随机选中的样本类别b的概率；（1-pb）为样本被分错的概率。

利用重构后的NDVI时序数据（NDVI Fit）、重构后的RENDVI783时序数据（RENDVI783Fit）、物候特征（RENDVI783Ph）、物候特征优选组合4种数据构建6种不同的分类特征数据集（表4）。

表4 用于作物分类的数据组合类型特征Tab.4 Combination characteristics of data types used in crop classification

3 结果与分析

3.1 植被指数随时间变化特征

作物的生长规律是播种、出苗、抽雄、成熟和采摘过程［26］，图2 表示研究区作物生长周期中植被指数的变化，总体上看4 种植被指数呈先上升达到峰值后下降的趋势。从图2a 可以看出，3 种作物的NDVI 曲线较为相似，生长季同时开始，在生长季长度和生长季结束期有较好的区分性。在NDVI时序曲线中，3 种作物的峰值均低于0.8，红边波段加入之后，所有作物的RENDVI 曲线的峰值均高于0.8，反映出较好的区分效果（图2b～d）。这是因为红边波段（705～783 nm）更深入于作物冠层和叶片，对冠层的微小变化和衰老更加敏感。其中，在RENDVI783时间序列中，作物的波峰值较高，约达到0.9。作物曲线之间的分离较好，每种作物的物候规律更加清晰和直观。这是因为构成RENDVI783的B7（783 nm）波段与B8波段更接近，从而能够更好地表达植物的物候信息。进一步分析图2d可以看出，棉花的平均生长期约在第180 d 左右，NDVI 值从第120～140 d 苗期开始升高，在第230～250 d 花铃期达到峰值，在第270～280 d棉花停止生长时下降，变化特征及其明显。而玉米和果树的生育期为130～150 d左右，第210～240 d 达到最大值，易与棉花区分。因此，最终确定利用NDVI和RENDVI783时间序列数据进行下一步平滑、重构研究。

图2 不同数据典型地物曲线Fig.2 Curves of typical ground object

3.2 时序数据重构

在提取物候特征之前，需要用滤波函数拟合重构植被指数时序影像，最大程度地减少噪声的影响，为后期分类提供准确的数据。图3 是纯像元滤波处理前后NDVI 和RENDVI783时间序列对比曲线。结果表明，经过S-G滤波拟合后，两种时序曲线的噪音和异常值明显降低，曲线更加平滑。平滑后的曲线与研究区3 种作物生长周期基本一致，峰值和低谷期与作物的生长盛期及收割期完全吻合。

图3 S-G滤波重构前后作物NDVI和RENDVI783时序曲线Fig.3 Crop NDVI and RENDVI783 time curves before and after S-G filtering

3.3 物候特征提取

作物的物候特征指标是分类的关键，NDVI 和RENDVI783拟合中每种作物在生长季的物候指标曲线如图4 所示。从图4a 与图4d 生长季开始和生长季结束曲线来看，NDVI 比RENDVI783曲线起伏程度较小，相对平缓。在生长季开始不同作物在RENDVI783上体现的更分散，有区分性。生长季结束指标显示RENDVI783出现不同作物集聚为2 组。从图4b和4e生长季基值、拟合函数最大值以及生长季振幅3 个物候指标值来看，RENDVI783比NDVI 在不同作物的物候指标有较大差异。拟合函数最大值物候指标显示，RENDVI783的波动幅度以及差别要大于NDVI。RENDVI783对作物更敏感，棉花与其他作物的差别比较大，最容易识别；从图4c 和4f 生长季长度、大积分、小积分物候指标中可以看出RENDVI783比NDVI 波动幅度大，在指标值方面RENDVI783比NDVI 整体高。生长季长度指标处NDVI 出现轻微的集聚，RENDVI783相对来说更分散，说明不同作物在生长季长度比较容易区分。因此最终提取重构后RENDVI783时序数据的11个物候特征影像并参与到下一步的分类实验。

图4 NDVI和RENDVI783同一作物的不同物候特征值Fig.4 Different phenological eigenvalues of the same crop in NDVI and RENDVI783

基于S-G 滤波拟合重建后的RENDVI783时间序列数据集，选取了11个关键季节性物候特征的指标影像（图5）。可以看出生长季结束期的0 值区域波动较少，在拟合函数最大值和生长季振幅的特征值主要集中在研究区北部和渭干河—库车河两岸，在生长季左导数和右导数的0值区域分布较广。

图5 物候特征影像Fig.5 Phenological feature image

3.4 特征重要性评价

利用RFC 方法的OOB 法，对上述提取的11 个物候特征进行重要性评估，特征的重要性得分如图6。生长季拟合函数最大值的重要性得分最高为1.43，对接下来的分类及信息提取的贡献率最大；其次为生长季长度（1.40）；第三为生长季振幅（1.23）；生长季结束为1.16；生长季大积分、生长季小积分的重要性得分分别为1.02、1.01；其余物候特征的重要性得分均小于1。因此，选取前6个重要性得分较高的特征组成一个优选特征数据集，参与到下一步的分类实验。

图6 RENDVI783时序物候特征变量的重要性排序Fig.6 Importance ordering of feature variables of RENDVI783 time-series phenology

3.5 分类结果及精度验证

3 种分类方法对6 种分类特征方案的分类总体精度从高到低依次为：F方案＞C方案＞E方案＞B方案＞D方案＞A方案。MLC、SVM和RFC 3种方法中分类效果最好的特征组合均为F方案（RENDVI783Fit+物候特征优选组合），总体精度和Kappa系数分别达到了82.60%、88.40%、92.20%和0.82、0.87、0.92（图7）。

图7 3种分类方法不同分类方案精度雷达图Fig.7 Precision radar map of three classification methods and different classification plans

RFC 方法中，对比A 方案(NDVI Fit)和B 方案（RENDVI783Fit），B 方案分类精度提高了1.30%，Kappa 系数提高了0.02，说明红边波段的参与提高了分类精度；对比C 方案（RENDVI783Ph）与E 方案（RENDVI783Fit+ RENDVI783Ph），E 方案的总体精度和Kappa系数提高了3.40%和0.04，说明物候特征明显提高了分类精度；对比C方案和D方案（物候特征优选组合），D方案的总体精度和Kappa系数提高了2.10%和0.03，说明分类精度不会随着特征波段的增加而增加；对比E 方案和F 方案，F 方案总体精度和Kappa系数提高了3.60%和0.05，关键的特征信息避免交叉冗余信息对棉花提取的干扰，使所有类型的生产者精度和用户精度都有较大程度的改善。

分别对6种特征数据集，3种不同的分类方法所获取的棉花分类结果计算混淆矩阵，具体结果见表5。

表5 3种分类方法不同分类方案棉花分类精度比较Tab.5 Comparison of cotton classification accuracy of three classification methods and different plans /%

MLC方法分类结果从图8a～f中可以看出，图像分类效果整体上比较细碎，棉花的分类效果较差，绿洲边缘地区、东部和西部的棉花种植区域明显稀少。玉米几乎没有被识别出来，大部分玉米误分为居民地。果树沿着渭干河—库车河流域与居民地交错分布，分类面积过多，有误分的情况。

图8 不同数据集分类结果对比Fig.8 Comparison of different dataset classification results

SVM 分类结果，从图8g～l 中可以看出，没有太过细碎的部分，棉花成片状分布在绿洲，“椒盐”现象明显减少。玉米在研究区中部交错分布，果树沿着渭干河—库车河流域与居民地交错分布。研究区的建筑分布较稀疏且散落，东北部分为较大建筑集中区，但是与MLC方法分类相比分类结果比较精确。

RFC 分类结果，从图8m～r 中可以看出，与上述两种分类结果相比效果比较理想。棉花在研究区由内而外大面积种植，主要在南部和西南部以及塔里木河两岸呈片状分布，因此棉花的识别效果是最好的。有些玉米分布在研究区东部，大部分交错分布在研究区中部。果树大部分沿着2条河流域与居民地交错分布，研究区北部也有大量的种植。由于多种作物套种等复杂的种植方式，使得分类仍然存在错分和漏分现象。该分类结果中，居民地的识别也有显著的提高，主要分布在库车市、新和县和沙雅县人口集中分布的区域，误分情况明显减少。

为了突出不同分类方法的分类效果，在3 种分类方法中分别选取F方案对应的3个子区域分类结果与原始影像进行对比分析（图9）。样地1 的3 种方法分类结果局部细节图（图9b～d）可以看出，MLC方法把裸地错分为居民地，SVM分类方法中裸地错分为水体；样地2 的3 种方法分类结果局部细节图（图9f～h）可以看出，MLC 方法裸地错分为玉米，SVM分类方法中裸地错分为居民地，棉花分类面积明显比RFC 方法少；样地3 的3 种方法分类结果局部细节图（图9j～l）可以看出，MLC 方法把裸地错分为居民地，SVM分类方法把裸地错分为点状分散地水体，RFC方法效果最好。因此总结出，RFC方法最适合种植模式复杂的研究区，能够有效地识别不同的作物，并得到高精度的提取。

图9 3种分类方法的分类结果局部图Fig.9 Partial graph of classification results of three different classification methods

根据RFC 方法F 方案分类结果，提取研究区棉花、果树、玉米3种作物的面积。研究区主要作物为棉花，种植面积约为3424 km2，占研究区总面积的24.67%；其次是果园，种植面积约为1216 km2，占研究区总面积的8.76%；玉米的种植面积最少，约为733 km2，仅占研究区总面积5.28%。

4 讨论

本研究分析了Sentinel-2 红边归一化植被指数时序数据和物候特征对作物分类精度的影响及特征的重要性。利用RFC 方法对棉花分布识别与种植面积提取，研究结果在一定程度上解决了以往研究［27］中大尺度遥感影像作物识别过程中混合像元众多和特征变量不足而造成的分类精度较低的问题。同时为时序Sentinel-2数据在干旱区作物的识别及面积提取研究提供了一定的参考。

本研究主要利用新数据源、特征变量的引入和特征优选等方法验证了RFC 方法对作物分类和信息提取的可行性。近年来基于物候特征对多时相遥感影像在农作物的信息提取中得到了广泛应用，如MODIS、Landsat 等数据构建的NDVI、EVI 时序曲线结合作物的季节性特征，较好地提取了甘蔗、水稻等作物的分布和面积信息，精度均比较理想［28-29］。但是对于地块破碎，种植类型比较复杂的区域，Landsat和MODIS数据难以解决混合像元的问题，而且仅靠单一的植被指数难以达到精确分类的目的。本文利用长时间序列Sentinel-2影像不仅提高遥感影像的分辨率，同时为棉花识别提供了敏感波段和有效特征，在时间和空间以及特征量方面扩充作物遥感监测技术与方法体系［30］。另外，本文利用S-G 滤波法重建的植被指数时序数据消除原数据存在的噪音和异常值的干扰，提高数据质量的同时保持农作物物候信息的真实性［31］。基于OOB 法的物候特征优选避免冗余信息对棉花识别与提取的干扰，使分类精度得到一定程度的提升［32］。最终基于重构后的红边植被指数时序数据和物候特征优选组合的RFC 方法精度达到92.20%。研究区棉花种植面积向外延伸，在南部、西南部和塔里木河两岸呈片状分布；但一些套种模式比较复杂的区域存在漏分、误分等现象，作物识别精度有待提高。

虽然本文的研究方法和思路，为克服大范围的干旱区农作物种植面积提取所面临的空间分辨率和参与特征种类不足等问题提供一定的参考，但也存在进一步探讨的问题。Sentinel系列卫星的Sentinel-1雷达数据有超强的穿透性、不受云量和天气的影响，与Sentinel-2 光学数据具有更好地融合优势。已有学者利用基于多时相双极化SAR数据，根据研究区作物的物候特征和不同作物的生长期后向散射系数，提取了研究区棉花种植面积，分类精度达到了90.88%［33］。鉴于此，以后研究的重点在于利用Sentinel 多源遥感数据的综合应用、特征优选方法构建和智能化分类方法建立，从而进一步提高分类精度。

5 结论

（1）从Sentinel-2 NDVI 和RENDVI783植被指数的棉花生育期变化规律可以看出，棉花NDVI 指数与RENDVI783变化趋势较为一致，在5 月（苗期）到8月初（开花盛期）有明显的上升趋势，在8 月末至9月（花铃期）达到峰值。相比NDVI 时序曲线，加入红边波段后的RENDVI783时序曲线峰值从0.7 提高到0.9，易与相互区分，体现出更佳的区分效果。

（2）利用OOB法对RENDVI783的11个物候特征进行重要性评价，生长季拟合函数最大值的重要性最高，为1.43，对分类及信息提取的贡献率最大。其次为生长季长度（1.40）；第三为生长季振幅（1.23）；生长季结束、生长季大积分、生长季小积分的重要性分别为1.16、1.02、1.01。

（3）RFC 方法对F 方案的分类精度最佳。总体精度和Kappa系数分别为92.20%和0.92。

（4）根据RFC 方法F 方案的分类结果，提取了研究区内棉花的种植面积，约为3424 km2，占研究区总面积的24.67%。