基于Sentinel-2时序谐波特征的县域农作物分类

孙庆松， 张晓楠， 陈利东， 王 坤， 宋宏利

(1.河北工程大学地球科学与工程学院，河北 邯郸 056038； 2.河北工程大学矿业与测绘工程学院，河北 邯郸 056038； 3.河北省地矿局第六地质大队，河北 石家庄 050000)

农作物精细化分类是农作物长势监测、产量预测及灾害评估的重要前提和基础。及时准确地获取农作物类型、面积和空间分布信息，对于优化农作物种植结构、实现农业生产精准化管理及确保国家和地区粮食安全具有重要意义[1-3]。传统农作物类型信息的获取主要通过实地走访调查，然而进行实地精细调查费时费力，且容易受到人为因素干扰，因此更新不够及时。随着对地观测技术的发展，利用遥感数据提取农作物信息已经成为重要手段。

目前，越来越多的研究者利用植被光谱特征的季节性变化来提高作物分类的准确性，已有学者利用多时相光学遥感数据提取增强型植被指数(EVI)、归一化植被指数(NDVI)、土壤调节植被指数(SAVI)等时序指数集开展农作物分类研究，并取得了较好的分类效果[4-8]。目前国内农作物信息提取主要是依据指数特征进行分类，利用谐波特征提取农作物信息的研究较少。谐波特征分类优势在于并不仅仅依靠某时相内指数值的差异，同时利用不同作物时序曲线中唯一对应的相位值和谐波余项值进行分类[9]。谐波分析[10-12]技术最早应用于电力信号分解和水文气候预测。在作物信息提取方面最初由Jakubauskas等[13]提出，他们利用谐波分析技术对26个年内NOAA气象卫星传感器(AVHRR)获取的NDVI时间序列进行分析，并提出了一种基于NDVI值随时间变化的作物识别方法。

目前常用的遥感数据源主要包括Landsat[14]、GF系列[15]、MODIS[16]及多源融合数据[17]。与Landsat相比，Sentinel-2卫星具有高时间分辨率和高光谱分辨率，同时具有3个红边波段，近年来被广泛应用于农作物分类研究[18-19]，并展现出其在农作物分类以及农情监测应用中的优势[20-21],为县域尺度范围农作物信息提取提供了新的选择。

以往研究主要集中于影响因素可控的小区域[22]和地势平坦的西北及东北区域[23-24]，而在地块破碎、种植结构复杂的华北县域范围内，关于农作物精细化分类的研究相对较少。因此本研究结合中国探索实行耕地轮作休耕制度的[25]背景，以中国地下水漏斗区黑龙港流域南宫市为研究对象，针对县域尺度范围内农作物种植信息提取精度低的问题，利用多时相Sentinel-2数据，通过构建指数时序集及其谐波特征，组合多种分类方案，对研究区内农作物信息进行提取，以期为县域尺度农作物分类研究提供参考。

1 材料与方法

1.1 研究区概况

南宫市位于河北省东南部(37°05′～37°27′N，115°08′～115°45′E)，黑龙港流域南部(图1)，属黄河冲积平原，地面平均海拔28.6 m，地势东南稍高西北低，坡降为1/7 000，地形相对平坦。属暖温带亚湿润季风型气候，光照充足，雨热同季，年均降水量476 mm。全市境内共有3条主干河道，24条分支，均属于黑龙港水系，总长367 km。全市面积863.3 km2，总耕地面积597.3 km2，占全市面积的69%，主要农作物有冬小麦、棉花、夏玉米、谷子、红薯、辣椒等，是华北地区重要产粮基地之一。

图1 研究区及样点分布Fig.1 Distribution of the research area and sample points

南宫市主要农作物是冬小麦、棉花、夏玉米、谷子、辣椒、红薯、花生，其中，冬小麦-夏玉米是一年两季传统轮作模式，其他农作物均属于一年一季。根据实地走访调查及邢台市统计局提供的农耕数据，总结出当地各类农作物物候历(表1)。

表1 南宫市主要农作物物候期

1.2 数据源及预处理

本研究使用的遥感数据是欧洲航天局(ESA)(https：//scihub. copernicus. eu/dhus/#/home)免费提供的Sentinel-2多光谱数据。Sentinel-2包含A和B两颗卫星，单颗卫星重访周期为10 d，两颗卫星同时运行重访周期为5 d，Sentinel-2卫星最高空间分辨率为10 m，具有13个波段数据。为保证能够完全覆盖研究区内各类作物生长周期，选取2019年10月至2020年10月共63幅云量低于20%的光学影像，影像获取时间信息如表2所示，其中影像最大时间间隔为30 d，最小为2 d。最大时间间隔在1月份，处于冬小麦的生长周期内，由于冬小麦在2月至6月期间，其植被指数值与其他作物有明显差异，因此在1月份时，影像间隔对分类结果的影响较小。

表2 遥感影像数据

采用经几何精校正的L1C级产品数据，根据植被指数计算公式，研究使用蓝、绿、红、2个红边波段及近红外波段数据(其中2个红边波段空间分辨率为20 m，可见光及近红外波段空间分辨率为10 m)。利用ESA提供的Sen2cor插件对L1C级数据进行大气校正，并将所需波段的空间分辨率和反射率利用三次卷积法重采样为10 m的分辨率和反射率。为便于后续调用各波段数据计算植被指数，在SNAP软件中将同一影像中的波段数据合成为一个文件，并将文件转换为ENVI可读取的TIFF格式，之后在ENVI软件中进行裁剪、镶嵌、植被指数计算、最大值合成等操作，进而构建植被指数时间序列。利用MATLAB及IDL工具包提取时序数据中的谐波特征(对应各阶振幅值图像)，之后利用随机森林分类器进行图像解译。

1.3 实地调查数据

2019年10月至2020年9月在南宫市进行多次实地调查。为避免手持全球定位系统(GPS)定位精度问题导致的作物种类边缘误差，在样本选择时，尽可能选择地块面积在30 m×30 m以上区域，在地块中心用手持GPS接收机标记地块经纬度坐标及地块类型，样本点均匀分布于全市(图1)。共采集实地样本1 090个，其中冬小麦547个，棉花282个，谷子118个，花生36个，辣椒12个，红薯10个，建筑用地及裸地、水体和其他地类共85个。由于研究区内花生、辣椒、红薯种植数量较少，为减少样本数量悬殊导致的分类误差，故合并为一类“其他农作物”。将2/3的实地采集数据作为训练样本，1/3的实地采集数据作为验证样本。

1.4 构建植被指数时序曲线

根据作物光谱反射特性，结合Sentinel-2卫星中心波长及波段范围，为获取最佳组合分类方案，构建5种植被指数，分别为EVI、NDVI、NDVI705[26]、EVI+NDVI、EVI+NDVI+NDVI705，在构建2种组合指数方案时，采用指数值算数相加原理，其计算公式及构建规则见表3。计算每幅影像的植被指数后，采用最大值合成法(从相邻影像中提取最大像元值重构一幅新的影像)，将全年63幅影像合成为30幅，确保每月有2～3幅影像，在保证时序数据集质量的同时，能最大程度减少云雾干扰，从而构建作物生长周期变化规律曲线。

3种植被指数(EVI、NDVI、NDVI705)被广泛应用于作物估产、环境监测、火灾评估等领域。NDVI指数对于波段范围限制较小，应用广泛，能够将植被与其他地物明显区分，由于其本质是近红外波段和红外波段的非线性拉伸，数值容易饱和，导致在高密度植被区灵敏度降低。与传统NDVI不同，NDVI705指数计算选用了叶绿素吸收特征较窄的波段，更容易受到植被内叶绿素含量的影响，在植被生长变化过程中，数值变化更加敏感[26]。EVI指数的提出是为了改进NDVI在植被高密度区灵敏度低的问题，加入蓝光波段，以减少大气折射造成的影响，其时序数据被广泛应用于林地变化监测研究。

表3 植被指数计算公式及构建规则

1.5 谐波特征

谐波分析又称傅里叶分析，是指将一个随时间变化的周期性函数通过傅里叶变换分解为无穷多个谐波分量，每个分量由频率相同的正弦波和余弦波组成。在农业遥感领域，谐波分析的本质是将单维空间域影像转换为多维频率域数据，对频率域内的波段数据进行级数分解，将不同植被物候变化相关的信号分隔开，最后再将频域影像转换为空间域影像。

谐波主要特征在于前三阶谐波的振幅值和谐波余项，三阶之后的谐波特征多是各类噪音误差[10](图2)。图2a中关于时间的任意曲线(时间函数)都可以分解为图2b中的正弦和余弦波，图2c中冬小麦-夏玉米在NDVI时序中的函数可以分解为图2d中的谐波特征分量，图2b和图2d谐波曲线的区别在于振幅值和谐波余项的不同，即谐波特征不同。Jakubauskas等[13]在2001年将谐波分析法应用于NDVI时间序列分解，并根据每个谐波分量在原始数据中占总方差的百分比进行农作物分类。谐波分析数学分解公式为：

(1)

式中，f(x)为时间序列函数，a0是谐波余项，L是周期函数的时间长度，x表示时间，an是第n次谐波余弦函数的振幅，bn是第n次谐波正弦函数的振幅。各类作物在NDVI时间序列中的谐波余项和振幅值如表4所示。

图2 谐波特征示意Fig.2 Schematic diagram of harmonic characteristics

表4 各类作物谐波特征值

1.6 分类方法

2 结果与分析

2.1 不同作物时序曲线变化分析

为探究各类农作物植被指数值在生长周期内的变化趋势，根据实地调查样本数据，多次提取各类农作物植被指数值并取平均值，构建农作物生长变化曲线。

图3表明，冬小麦-夏玉米的时序曲线在3个指数波段范围内与其他地类均有明显差异。冬小麦在10月上旬播种后，育苗期一直持续到12 月上旬，同时植被指数值升高；12月中旬之后进入越冬期，气温降低，光合作用缓慢从而导致植被指数值降低；次年1月中旬至3月中旬冬小麦进入返青发育期，植被指数值迅速增长；在3月中下旬时期，该地区受到北方寒潮影响，温度骤降，植被指数值有小幅度下降；在4月至5月期间，冬小麦进入孕穗抽穗期，植被指数值迅速升高，达到生长发育期内的植被指数最大值，该时期冬小麦与其他作物物候差异最大，可以有效将冬小麦与其他作物区分开。夏玉米在6月份播种，之后进入育苗期，各类指数值呈上升趋势，7月至8月中旬进入抽穗灌浆期，指数值到达第2个峰值，9月份以后进入成熟收获期，指数值缓慢下降。

a:2019-10-05; b:2019-10-15; c:2019-10-30; d:2019-11-04; e:2019-11-14; f:2019-11-19; g:2019-12-04; h:2019-12-19; i:2020-01-03; j:2020-01-30; k:2020-02-17; l:2020-02-22; m:2020-03-10; n:2020-03-13; o:2020-03-18; p:2020-03-23; q:2020-04-12; r:2020-04-17; s:2020-04-22; t:2020-04-27; u:2020-05-02; v:2020-05-22; w:2020-06-03; x:2020-06-21; y:2020-07-06; z:2020-07-16; A:2020-08-22; B:2020-08-30; C:2020-09-04; D:2020-09-19。图3 各类农作物时序曲线Fig.3 Time series curves of different crops

棉花在4月中下旬播种，5月中下旬开始出苗，在此期间，植被指数值没有明显变化。之后经历现蕾、开花阶段，植被指数值急速上升。7月份开花以后，其植被指数值基本维持稳定，与其他非轮作农作物时序曲线在3种指数波段范围内均有明显差异，开花之后的指数值可以作为识别棉花的重要特征。9月份进入成熟期，植被指数值下降。谷子在5月中下旬播种，在出苗、拔节、抽穗时期，植被指数值呈上升趋势，8月份成熟，此时植被指数值达到峰值。

辣椒和红薯均在4月下旬播种，到6月上旬，2种作物植被指数值无明显变化；6月下旬辣椒进入分蘖期，红薯进入现蕾期，植被指数值上升；到8月份，植被指数到达峰值，成熟以后其植被指数值缓慢下降。花生在5月中旬播种，进入出苗阶段，植被指数值开始增长；6月份至7月份为幼苗期，植被指数值迅速升高；到8月份结荚时期，植被指数值保持相对稳定；9月份成熟时收获，植被指数值缓慢下降。由于这3类作物种植面积较小，因此在分类时，将其合并为一类。由图3可以看出，依靠单时相影像数据，并不能将生长周期相似的农作物明显区分，而时间序列数据可以凸显物候差异特征，从而提高分类准确率。

2.2 不同指数组合的农作物分类方案精度对比

为探讨不同指数组合及谐波特征对农作物分类精度的影响，选用EVI、NDVI、NDVI705、EVI+NDVI、EVI+NDVI+NDVI7055种植被指数，构建指数时序曲线，并提取作物生长曲线谐波特征。将分类特征和样本数据输入到随机森林分类器中，结合实地采集的验证样本，采用混淆矩阵验证分类结果，评价指标包括总体精度(Overall accuracy, OA)、用户精度(User accuracy, UA)、生产者精度(Producer accuracy, PA)、Kappa系数，结果如表5、表6、表7所示。

表5 各种分类方案总体分类精度及Kappa 系数

从原始分类结果中可以看出，当指数特征和谐波特征共同作为分类依据时，相比于仅利用指数特征进行分类，5种时序集的总体分类精度均有明显提高(表5)。其中EVI+NDVI时序集提升最为明显，提高了9.21个百分点，其次是EVI+NDVI+NDVI705时序集，提高了8.75个百分点，NDVI705时序集提升最不明显，提高了8.14个百分点，表明谐波特征的加入可以有效提高分类精度。

为进一步验证不同植被指数对于农作物分类精度的影响，选取指数特征和谐波特征共同作为分类依据的5种指数组合方案，并对其精度评价指标进行详细分析，结果见表6、表7，5种原始分类结果如图4所示。

表6 3种单一指数的农作物分类方案精度评价结果

表7 2种组合指数的农作物分类方案精度评价结果

表6表明，以3种单一植被指数作为特征波段时，采用EVI指数特征作为分类依据的总体精度最低，为83.82%，利用NDVI705指数分类时总体精度最高，其值为90.00%，NDVI分类精度比EVI分类精度提高了约5个百分点。NDVI705分类结果说明，实际验证样本与模型分类结果具有高度一致性，Kappa系数为0.88；NDVI其次，Kappa系数为0.86；EVI最低，Kappa系数为0.80。在3种分类方案中，农作物中冬小麦-夏玉米的用户精度和生产者精度均为最高，从种植模式上分析，研究区内在1月份至6月份期间只有冬小麦种植(图3)，因此物候特征明显，分类精度较高；棉花和谷子的分类精度略低，主要是由于其生长周期相似，且两类作物种植结构呈嵌套式分布(图4)，容易造成错分情况；建筑用地、裸地和水体的光谱反射率并无明显季节变换规律，与农作物物候特征有明显差异，因此分类精度较高。NDVI705指数比NDVI指数分类精度高，从算法原理上分析，较窄的红边波段相比于红光波段对于植被生长变化程度更加敏感，说明Sentinel-2特有的红边波段数据在农作物精细化分类上具有较大的应用潜力。

2种组合指数分类方案中，EVI+NDVI+NDVI705作为特征组合时，总体分类精度最高,为94.95%，比EVI+NDVI分类精度提高了约2个百分点(表7)，说明红边波段指数NDVI705的加入能够明显提高分类精度；EVI+NDVI组合指数比NDVI指数分类精度提高了约4个百分点，表明具有蓝光波段的EVI指数的加入能够提高分类精度。

EVI:增强型植被指数；NDVI：归一化植被指数；NDVI705：红边归一化植被指数。图4 不同指数组合的农作物分类效果Fig.4 Crop classification effects of different index combinations

2.3 南宫市作物种植空间分布

南宫市农作物分类总体精度最高的指数组合方案是EVI+NDVI+NDVI705，其分类效果见图4。南宫市2020年冬小麦-夏玉米、棉花、谷子、其他农作物的种植面积分别为21 580.77 hm2、16 968.16 hm2、3 000.81 hm2、403.36 hm2。种植面积最大的作物是冬小麦-夏玉米，主要分布在西北部和东南部；其次是棉花，主要位于中部和南部。冬小麦-夏玉米和棉花是当地重要的经济作物，因此呈现大面积块状分布。谷子的种植分布比较零散，且地块面积较小，多聚集在居民地附近，主要与当地的饮食习惯有关。

3 结论与讨论

本研究基于Sentinel-2光谱数据，通过构建多种植被指数及其谐波特征，采用随机森林算法，对南宫市内的主要农作物进行分类提取，结果表明：1)当指数特征和谐波特征共同作为分类依据时，相比于利用单一指数特征进行分类，5种时序集的总体分类精度均有明显提高，EVI+NDVI时序集提升最为明显，提高了9.21个百分点，NDVI705时序集提升最不明显，提高了8.14个百分点，表明时序集中谐波特征的加入可以有效提高分类精度。2)在指数特征和谐波特征共同作为识别依据的分类方案中，基于EVI+NDVI+NDVI705组合指数的分类方案精度最高，比EVI+NDVI分类精度提高了2.57个百分点，而EVI+NDVI组合分类精度比NDVI指数分类精度提高了3.53个百分点。以上结果表明，红边NDVI705和EVI的加入能够提高作物识别精度。3)通过对3种单一指数时序曲线特征的分析，可以明显看出冬小麦-夏玉米与其他农作物生长周期不一致，而其他农作物物候周期相似，且指数值峰值无明显差异，同时分类结果表明，利用单一指数进行分类并不能获得最佳的分类效果，验证了多特征组合分类方法在农作物精细化分类上的可行性。

尽管本研究取得了较好的分类效果，但仍存在一些问题。在作物生长关键期，单一数据源的影像获取频率较低，本研究没有采用多源遥感来提升遥感数据的获取概率。且本研究仅考虑了EVI和NDVI等植被指数，其他植被指数没有纳入考虑范围。因此在今后的工作中可以结合高分系列遥感数据，融合其他植被指数特征开展研究。