EN中文

Journal of Geo-information Science >

2022 , Vol. 24 >Issue 5: 996 - 1008

DOI: https://doi.org/10.12082/dqxxkx.2022.210527

Identifying Crop Types in Early Growing Season in the Middle Reaches of Heihe River by Fusing Multi-Source Remote Sensing Data

  • YANG Zehang ,
  • WANG Wen , * ,
  • BAO Jianxiong
Expand
  • State Key Laboratory of Hydrology- water Resources and Hydraulic Engineering science, Hohai University, Nanjing 210098, China
* WANG Wen, E-mail:

Received date: 2021-08-31

  Request revised date: 2021-11-08

  Online published: 2022-07-25

Supported by

International (Regional) Cooperation and Exchange Program of National Natural Science Foundation of China(41961134003)

Copyright

Copyright reserved © 2022.
Fold

Abstract

Obtaining the condition of crops in the early growing season is important for the agricultural water resources management, especially for water allocation in water shortage areas. In this paper, the Enhanced Spatiotemporal Adaptive Fusion Model (ESTARFM) was used to fuse Sentinel-2 NDVI data with NDVI data calculated from MOD09GQ in the early growing season in the middle reaches of the Heihe River basin from March to June to construct NDVI time series, and then the random forest classification method was used to identify the crop structure in 2019. By using the Sentinel-2 NDVI time series and spatiotemporal fused NDVI time series jointly in the early growing season (i.e., from March to June), the crop classification accuracy reached 91.42 % with a kappa coefficient of 0.85. Compared with the classification results obtained using only Sentinel-2 NDVI time series, the accuracy improved by 1.01 % and the kappa coefficient improved by 0.02. The accuracy was only 1.53 % lower and the kappa coefficient was only 0.03 lower than the results obtained using Sentinel-2 NDVI time series in the whole crop growth period (from March to October). The Gini coefficient was used to evaluate the feature importance of the time series constructed by combining the Sentinel-2 NDVI time series and fused NDVI series. It was found that six of the 10 NDVI images with Gini coefficient scores higher than the average were spatiotemporally fused images, indicating that the NDVI data obtained by spatiotemporal fusion could improve the classification accuracy. The NDVI time series of different lengths were also established separately for early identification of crop structures. Results show that early identification of alfalfa and corn can be achieved as early as mid-April and late April, respectively. In addition, the classification accuracy of corn was strongly influenced by the length of NDVI time series, and early identification of corn can be achieved in late May.

Key words： crops identification; remote sensing; spatial-temporal data fusion; NDVI time series; the middle reaches of the Heihe River basin

Cite this article

YANG Zehang , WANG Wen , BAO Jianxiong . Identifying Crop Types in Early Growing Season in the Middle Reaches of Heihe River by Fusing Multi-Source Remote Sensing Data[J]. Journal of Geo-information Science, 2022 , 24(5) : 996 -1008 . DOI: 10.12082/dqxxkx.2022.210527

1 引言

农作物种植结构监测对于农业水资源管理、作物估产以及农业政策的制定等有着重大意义[1]。农业遥感中,基于一期或少数几期多光谱遥感影像进行作物种植结构提取是较为快速简单的方法。如曹卫彬等[2]基于2003年9月14日Landsat TM影像数据,通过分析棉花与其它主要作物在各波段上的光谱差异,研究开发了基于光谱特征的棉花识别模型。Mathur等[3]利用2003年9月22日IRS-1D影像数据的光谱单波段特征,采用支持向量机分类方法获得印度旁遮普区域棉花和水稻两大农作物的空间分布,精度达到91%。上述方法虽然快速简单,但也存在明显不足,如难以获得最佳物候期的影像数据,当区域作物种植结构复杂时不易区分光谱相近的不同作物。
针对单期遥感数据在作物种植结构提取中的不足,很多学者又提出了第2种方法,即基于时间序列影像的农作物种植结构提取方法。此方法可以充分利用不同作物的物候变化规律,更好地提取不同作物种植结构[4]。如边增淦等[5]利用18景GF-1 NDVI数据建立时间序列,利用多分类器组合的方法,成功对黑河中游地区作物种植结构进行提取,精度达到94.19%。Zhong等[6]引入GDD（累积积温）,并采用随机森林分类器对光谱波段、GDD、EVI、NDSVI（归一化衰败植被指数）等多个特征量以及它们之间的组合进行测试,发现四类特征量共同参与分类能获得最高的整体分类精度。
由于大多数农作物的分类产品都需要基于全年的多时相遥感数据,因此往往只能在年底或者第二年才能获取作物的种植结构,对于实际应用存在滞后性。如果可以在作物生长早期获取作物的种植情况,对于农业水资源管理,尤其是缺水地区的水量分配等具有更为重大的意义。关于作物早期识别主要有2种思路。
（1）根据历史作物生长早期遥感数据建立时间序列作为参考,从而进行当年作物的早期识别。如Cai等[7]基于2000—2015年1322幅Landsat影像数据,将2000—2014年Landsat数据建立时间序列用于训练深度神经网络模型,使用2015年的时序数据预测当年的作物类型,结果表明使用7月前的时序数据解译伊利诺斯州尚佩恩县的玉米和大豆准确率达96%。但此方法只针对特定作物提出,在作物种植结构较为复杂地区实用性还有待检验。郝鹏宇等[8]基于2006-2014年MODIS EVI时间序列数据和作物数据层(CDL)数据,分别建立美国堪萨斯州2006—2014年主要作物EVI时间序列和确定2014年训练样本,使用2014年Landsat NDVI月合成时间序列进行作物识别,结果表明当时间序列长度为4—8月时作物分类总体精度接近使用完整时间序列数据的作物分类总体精度。但是MODIS空间分辨率较低,每个像元中存在混合像元,因此从中获取的样本精度还有待检验。
（2）仅使用当年作物生长早期影像,通过引入多种遥感指数或使用多源遥感数据融合的方法,增加作物生长相关信息量,从而进行作物的早期识别。如Zhou等[9]基于4—10月无云250m MODIS 10 d合成数据,对加拿大萨斯喀彻温省南部农业区进行作物识别,结果表明在较短的时间序列数据周期下土地覆盖分类可能达到最高的精度。但是,MODIS数据像元中往往含有混合像元,而混合像元时间序列曲线较为混乱,因此,此方法在田块破碎度较高地区会对分类结构造成较大影响。针对MODIS数据空间分辨率较低的问题,Villa等[10]通过整合2013年和2014年4—7月30 m空间分辨率的Landsat 8 3个光谱指数和雷达x波段数据获得天气季节特征,基于决策树分类方法获得意大利伦巴第地区主要作物类型的作物图,作物分类总体精度大于86%。Azar等[11]基于2013年26景Landsat 8数据构建EVI、NDFI、RGRI时间序列,对意大利北部伦巴第地区的7种作物进行早期识别研究,结果表明当使用7月中旬前的时间序列分类时作物总体精度超过85%。但是,Landsat 8数据时间分辨率较低,加上云量影响,可能会丢失作物生长关键信息,在实际使用中,受到的限制较大。
本文利用遥感影像时空融合方法,融合作物生长季早期的哨兵二号(Sentinel-2)影像与MODIS 影像数据,从而解决单卫星可用遥感数据较少以及不同卫星遥感数据之间存在偏差的问题,建立高时间分辨率和高空间分辨率的时间序列,更好地实现作物的早期识别;并通过对比使用不同长度NDVI时间序列的作物识别结果,探究可以有效识别作物种植结构的最早时间。

2 研究区概况及数据来源

2.1 研究区概况与实地调查数据

黑河流域是我国西北地区第二大内陆流域,属于甘肃省,中游地区东西走向位于97.2 °E—102.1°E之间,南北走向位于38.1 °N—39.6 °N之间。本文的研究区为黑河流域中游地区的三个县级行政区,分别为甘州区、临泽县和高台县。研究区平均海拔约1422.9 m,研究区内年平均气温6~8 ℃,年平均降水量在110~370 mm之间。主要种植作物为玉米、春小麦、苜蓿和蔬菜等。在2019年5月和7月分别在研究区展开实地调查,获取一系列作物种植情况样本。其中苜蓿样本82个,玉米样本1470个,小麦样本208个,蔬菜及其他样本321个,按照7:3的比例,随机分为训练样本和检验样本。研究区和采样情况如图1所示。主要作物物候期见表1。
图1 黑河流域中游研究区与采样情况示意

Fig. 1 Schematic diagram of study area and sampling in the middle reaches of the Heihe River Basin

表1 黑河流域中游地区主要作物物候期

Tab. 1 Phenology of main crops in the middle reaches of Heihe River Basin

4月 5月 6月 7月 8月 9月 10月
玉米 播种 抽穗 收获
小麦 播种 抽穗 收获
苜蓿 播种 收割 收割 收割

2.2 遥感数据

本次研究使用数据包括哨兵2号（Sentinel-2）卫星影像数据与MOD09GQ产品数据。
Sentinel-2卫星是高分辨率多光谱成像卫星,携带一个多光谱成像仪(MSI),由2A和2B两颗卫星构成,组网之后重访周期为5 d。本次研究所用数据为Sentinel-2 L1-C级数据,从欧空局哥白尼数据中心 (https://scihub.copernicus.eu/) 下载获取,已经过正射校正和几何精校正,未进行大气校正。因此,使用欧洲航天局(ESA)发布的对Sentinel-2数据进行辐射定标和大气校正的插件Sen2Cor(http://step.esa.int/main/third-party-plugins-2/sen2cor/)对数据进行预处理,获取最终可用的影像数据,数据的投影坐标系为WGS_1984_UTM_zone_47N。本文使用了Sentinel-2的红波段 (0.665 μm)与近红外波段(0.842 μm),其空间分辨率皆为10 m。
MOD09GQ产品数据属于陆地2级标准产品,拥有250 m分辨率的MODIS 1、2波段逐日表面反射率,已进行了辐射定标和大气校正处理。利用MODIS投影转换工具MRT提取MOD09GQ数据的红波段、近红波段反射率数据层,并对其进行重投影,使其投影坐标系与Sentinel-2 MSI数据保持一致。由于原始的Sentinel-2卫星数据与MOD09GQ产品数据中有部分数据本身质量较低,且存在云量较多,无法进行有效除云的问题,故对所有数据进行筛选,挑选出无云影像以及云覆盖量少于10%的影像,研究最终所用数据如表2所示。
表2 黑河流域中游地区2019年Sentinel-2与MOD19GQ影像数据

Tab. 2 Sentinel-2 and MOD19GQ image data from the middle reaches of the Heihe River Basin in 2019

数据类型 成像时间
Sentinel-2 上半年 3月2日 3月12日 3月17日 4月6日 4月16日 4月21日
5月1日 5月11日 5月16日 6月15日 6月30日
下半年 7月5日 7月15日 7月25日 7月30日 8月4日 8月14日
9月3日 9月23日 9月28日 10月18日 10月28日
MOD09GQ 上半年 3月2日 3月4日 3月10日 3月11日 3月12日 3月17日
3月18日 3月24日 4月6日 4月14日 4月16日 4月20日
4月21日 4月22日 4月30日 5月1日 5月3日 5月11日
5月12日 5月13日 5月16日 5月22日 5月23日 5月25日
6月3日 6月7日 6月15日 6月16日 6月17日 6月30日

3 研究方法

本文通过建立黑河中游地区作物生长早期Sentinel-2 NDVI时间序列,对作物种植情况进行早期识别。考虑到云量对作物种植结构提取的影响,利用MNSPI[12]去云算法对影像进行去云处理。此外针对作物生长早期可用影像数据较少的问题,利用ESTARFM[13]时空融合算法,将Sentinel-2 NDVI影像数据与通过MOD09GQ红波段和近红外波段计算获取的NDVI影像数据进行时空融合,构建时间分辨率更高的NDVI时间序列,更好的获取作物种植情况。并使用整个作物生长期的Sentinel-2 NDVI时间序列对作物进行提取,对比不同时间序列下作物分类精度,具体研究路线如图2所示。
图2 多数据融合的作物早期识别技术路线

Fig. 2 Multi data fusion for early crop identification technology route

3.1 NDVI时间序列重构

3.1.1 MNSPI影像去云算法
邻域相似像素插值方法（MNSPI）是Zhu等基于邻域相似像元插值方法（NSPI）[14]提出一种改进的除云方法。其假设云覆盖像元周边的相似像元在时间与空间上存在连续性,并根据不同云覆盖像元的空间位置差异对光谱空间信息预测值与时间信息预测值进行加权运算,最终目标像元的预测值为：
L x , y , t 2 , b = L 1 x , y , t 2 , b / r 1 + L 2 x , y , t 2 , b / r 2 / 1 / r 1 + 1 / r 2
式中： L 1与 L 2分别为基于光谱空间信息和时间信息的目标像元预测值; r 1表示目标像元与其相似像元的平均空间距离; r 2表示目标像元与云中心的空间距离。
L 1与 L 2分别由以下公式计算得到：
L 1 ( x , y , t 2 , b ) = ∑ i = 1 N W i × L ( x i , y i , t 2 , b )
L 2 ( x , y , t 2 , b ) = L ( x , y , t 1 , b ) + ∑ i = 1 N W i × ( L ( x i , y i , t 2 , b ) - L ( x i , y i , t 1 , b ) )
式中： L ( x i , y i , t 1 , b )与 L ( x i , y i , t 2 , b )分别为 t 1与 t 2时刻第i个相似像元b波段的光谱值,N为相似像元的数量, W i为第i个相似像元的权重,由第i个相似像元与目标像元之间的欧几里得距离 D i和光谱距离 RMS D i计算而来的, RMS D i是由式（4）计算。
RMS D i = ∑ b = 1 K ( L ( x i , y i , t 1 , b ) - L ( x , y , t 1 , b ) ) 2 K
式中： K是波段的个数。
3.1.2 ESTARFM时空融合算法
有研究表明,时空数据融合方法直接应用于植被指数,由于误差传播较小,比融合各波段反射率后再计算植被指数具有更高的精度[15]。因此,本文直接利用Sentinel-2 NDVI与MOD09GQ计算得到的NDVI数据进行时空融合处理。为防止混淆进行以下说明：将Sentinel-2数据、MOD09GQ数据、Sentinel-2与MOD09GQ时空融合数据的归一化植被指数分别简记为NDVISen、NDVIMOD、NDVIFUS。
ESTAREM时空融合模型需要输入二组相同时间的高空间分辨率和低空间分辨率数据以及一期预测时间的低空间分辨率数据。因此在本文中,在进行时空融合处理时需要使用二组相同时间的NDVISen与NDVIMOD数据以及一期预测时间的NDVIMOD数据。而在时空融合数据数据筛选中,结合4.2中的结果可以发现,参与时空融合的数据与预测时间越接近,预测结果与真实结果的相似度更高,因此在具体的ESTARFM时空融合实施过程中,选用与预测时间最为接近的NDVISen与NDVIMOD数据用于进行时空融合处理。
3.1.3 Savitzky-Golay时间序列滤波
原始遥感数据集受云、气溶胶、太阳高度角等因素的影响使数据存在很多的噪声,影响了数据分析和应用的效果[16]。因此,需要对时间序列进行去噪处理。Savitzky-Golay(S-G)滤波是一种通过局部多项式回归模型实现平滑时序数据的时域低通滤波方法,其基本思想是基于多项式在滤波窗口内利用最小二乘法对数据进行最佳拟合[17],如式（5）所示。
Y j ' = ∑ i = - n n C i Y j + 1 N
式中：Y为原始数值;Y′为拟合值;Ci为第i个点的权重;N=2n+1为滤波窗口的大小。

3.2 作物分类方法

3.2.1 随机森林分类
本文采用随机森林分类法对NDVI时间序列进行分类,达到对作物的有效提取。随机森林是 Breiman等[18]提出来的一种基于多决策树投票决策的分类方法,通过在数据及其对应的特征变量上的采取随机策略Laura随机重采样来构建若干CART类型决策树(不剪枝),然后通过投票的方式对多决策树分类结果进行后处理进而确定数据的类别归属。随机森林算法调整参数少,对于训练数据的减少和噪声具有很强的鲁棒性且会对泛化误差进行内部无偏估计,具有较高的高分类精度,能够在高维数据集上高效地运行[19]。而随机森林的不足则在于解决噪声较大的分类或拟合问题时,容易产生过拟合现象[20]。随机森林分类过程中需要设定2个重要参数,一个是随机森林树的数量,树的数量越多,构建所耗时间越长,在此次研究中统一设置为100。另一个为每棵决策树在分类中使用的特征数量,数量值越高,每棵决策树稳定性也越高。本次研究中,每棵决策树分类使用的特征数量取总特征数量的平方根。
3.2.2 评估方法
本文评估主要包括去云与时空融合效果评估、作物分类精度评估以及作物识别特征可分性与重要性评估3个方面。去云与时空融合效果评估中,使用均方根误差（Root Mean Square Error, RMSE）来测量预测值和真实影像之间的差异;使用决定系数（Determination Coefficient, R2）表示预测影像和真实影像的回归关系。
在作物分类精度评价中利用总体分类精度、Kappa系数、制图精度和用户精度等对不同作物的分类精度进行分析,所有定量指标均可由混淆矩阵计算获得到。而作物识别特征重要性和可分性则分别利用Gini系数与JM(Jeffries-Matusita)距离进行评估。Gini系数是Python Sklearn库(http://www.lfd.uci.edu/~gohlke/pythonlibs/#)随机森林包中用于评价作物识别特征重要性的指数,做了归一化处理。Gini系数得分越高,对分类结果的稳定性影响越大,作物分类越重要。JM距离则可用于作物可分性判断,其计算公式为：
JM = 2 × 1 - e - B
式（7）中 B为巴氏距离,计算公式为：
B = 1 8 D 2 + 1 2 ln ∑ j + ∑ k 2 ∑ j ∑ k
其中 D 2计算公式为：
D 2 = μ j - μ k T ∑ j + ∑ k 2 - 1 μ j - μ k
式中： μ j和 μ k为参与计算的地物类型的预期值, ∑ j和 ∑ k为对应的地物类型协方差矩阵的无偏估计,JM距离取值范围为0~2,值越大,代表可分性越高。

4 结果分析

4.1 遥感数据处理结果

4.1.1 遥感影像去云结果
在本次研究中,将筛选后的云覆盖量少于10%的Sentinel-2影像和MOD09GQ数据,利用目视解译的方法,对研究区内的云进行识别提取。需要进行去云处理的数据如表3所示。
表3 需进行去云处理的数据

Tab. 3 Data to be de-clouded

数据类型 有云日期
Sentinel-2 上半年 3月17日 4月6日 5月1日 5月11日 5月16日 6月15日 6月30日
下半年 7月15日 8月4日 9月3日 9月23日 10月18日
MOD19GQ 上半年 3月17日 4月6日 4月20日 4月22日 4月30日 5月1日 5月3日
5月11日 5月12日 5月16日 5月23日 6月7日 6月15日 6月16日
6月17日 6月30日
应用MNSPI去云算法前对该算法在研究区的精度和适用性进行评价。选用同等大小的5月16日的Sentinel-2 MSI影像和5月11日的Sentinel-2 MSI影像,在其中绘制面积大小约为4.1 km2的耕地范围,用于模拟有云区域,进行模拟去云处理。同样的,使用相同大小范围的5月13日MOD09GQ影像对5月11日MOD09GQ影像进行模拟去云实验,结果如图3所示。
图3 Sentinel-2和MOD09GQ去云模拟实验真实反射率与模拟反射率散点图

Fig. 3 The scatterplot of real and simulated reflectance in Sentinel-2 and MOD09GQ cloud removal simulations

从图3中可以发现,在Sentinel-2和MOD09GQ的红波段与近红外波段,真实反射率与模拟反射率散点图的分布较为均匀,说明除云模拟结果与真实图像之间具有较高相关性。且从中可以看出Sentinel-2和MOD09GQ红波段和近红波段的均方根误差均较低,且决定系数都达到了0.85以上,相似度高。由此说明,MNSPI去云算法模拟去云影像与真实影像的总体光谱信息十分接近,除云效果良好。4.1.2 遥感影像时空融合结果将利用MNSPI模型去云后的影像数据与原始无云数据应用于时空融合处理。应用ESTARFM进行数据融合前,应先对该算法的精度和适用性进行评价。分别利用以下5组数据（① 3月2日、3月17日NDVISen和3月2日、3月12日、3月17日NDVIMOD预测3月12日NDVIFUS;② 4月6日、4月21日NDVISen和4月6日、4月16日、4月21日NDVIMOD预测4月16日NDVIFUS;③ 5月1日、5月16日NDVISen和5月1日、5月11日、5月16日预测5月11日NDVIFUS; ④ 4月21日、6月15日NDVISen和4月21日、5月11日、6月15日预测5月11日NDVIFUS;⑤ 5月16日、6月30日NDVISen和5月16日、6月15日、6月30日日预测6月15日NDVIFUS。）获取其时空融合结果,对比精度差异,从而分析ESTARFM模型的适用性,以及影响时空融合精度的因素。结果如图4所示。
图4 不同时间ESTARFM时空融合NDVI影像与真实NDVI影像散点图

Fig. 4 12 March ESTARFM fusion image effect simulation

从图4中可以发现NDVI时空融合预测值与真实值的散点图分布均匀,说明时空融合生成影像与真实影像之间的具有一定相似性。图中所有时空融合预测结果的RMSE均小于0.10,且决定系数都超过0.80,说明真实的NDVI值与预测的NDVI值之间的误差较小,拟合性较高。此外,随着时间的增长,由于作物生长速度增快,导致地物NDVI值变化速度增快,时空融合数据的精度有所降低。对比2组5月11日的ESTARFM融合结果,可以发现使用与预测时间更近的数据所得到的结果,其误差更小,相关系数更高,更能反映真实的地物情况。综上所述,ESTARFM时空融合模型能在一定程度上对缺失的NDVI值进行补充,可用于后续研究。4.1.3 NDVI时间序列去噪处理将使用ESTARFM时空融合后获取的NDVIFUS数据与原始NDVISen数据按照时间顺序进行排列,构成新的NDVI时间序列,利用S-G滤波对其进行去噪处理。分别对每种作物选取3个样本,绘制其NDVI变化曲线,结果如图5所示。
图5 各作物S-G滤波前后NDVI时序曲线比较

Fig. 5 Comparison of NDVI time series curves of different crops before and after S-G filtering

从图5中可以发现,滤波之后各类型作物生长曲线产生了明显的平滑效果,生长曲线特征更为突出,说明S-G滤波在一定程度上对数据存在的噪声实现了过滤。玉米的生长曲线在6月3日之前位于一个较低的水平,而在6月3日出现急速攀升。小麦的生长曲线在5月11日与6月16日左右出现2个波峰。苜蓿由于其生长过程中会存在收割现象,因此其生长曲线起伏较大。蔬菜的生长曲线并没有展现出明显的规律性。利用滤波后不同作物在生长早期的NDVI变化曲线差异,可以实现对地表作物种植结构的提取,从而达到对作物类型生长季早期识别的目的。

4.2 分类结果与精度评价

4.2.1 分类结果与精度评价
图6为分别利用不同时间序列数据（① 2019年3—6月11期Sentinel-2 NDVI数据;② S-G滤波后2019年3—6月11期Sentinel-2 NDVI和19期时空融合数据;③ 2019年3—10月整个作物生长期的22景Sentinel-2 NDVI数据）,采样随机森林的分类方法对作物种植结构进行提取,并对分类后结果进行主要分析的最终结果。
图6 作物分类结果

Fig. 6 Diagram of crop classification results

观察图6中的3幅分类图,发现不同土地覆被类型的空间分布大致相同。进一步针对基于3—6月11期NDVI时间序列、3—6月30期NDVI时间序列以及全年22期NDVI时间序列获得的研究区土地覆被分类结果,使用验证样本分别对其分类精度进行评价并比较,结果如表4所示。
表4 黑河中游地区2019年基于不同NDVI时间序列作物分类精度比较

Tab. 4 Comparison of crop classification accuracy based on different NDVI time series in the middle reaches of Heihe River in 2019

作物类别 3—6月11期NDVISen 3—6月11期NDVISen+8期NDVIFUS 3—10月22期NDVISen
用户精度/% 制图精度/% 用户精度/% 制图精度/% 用户精度/% 制图精度/%
苜蓿 98.84 91.20 99.48 91.26 98.99 94.65
小麦 76.67 86.07 80.17 90.44 80.01 88.23
玉米 88.50 97.20 88.67 98.44 91.80 98.47
蔬菜 70.61 43.03 80.99 43.33 83.79 49.85
总体精度/% 90.37 91.42 92.95
kappa系数 0.83 0.85 0.88
观察表4可以发现,增加19期NDVIFUS后,作物总体分类精度提高了1.05%,kappa系数提高了0.02,说明获得更细致的3—6月NDVI时间序列,增加对地“观测”频率,有利于提高黑河中游研究区作物分类精度。且增加19期NDVIFUS后的3-6月NDVI时间序列,相对利用整个作物生长期的22期NDVISen进行种植结构提取的结果精度仅相差1.53%,kappa系数仅差0.03,说明利用时空融合增加对地“观测”频率的方法,能够获取足够的作物分类精度,达到对作物类型生长季早期识别的目的。
4.2.2 特征可分性与重要性评估
利用不同作物的训练样本,计算其在不同时间序列数据（① 2019年3—6月11期Sentinel-2 NDVI数据;② S-G滤波后2019年3—6月11期Sentinel-2 NDVI和19期时空融合数据;③ 2019年3—10月整个作物生长期的22景Sentinel-2 NDVI数据）中的JM距离,从而对各作物的特征可分性进行评估,结果如图7所示。
图7 不同NDVI数据组合时间序列中每种作物与其他3种作物JM距离平均值

Fig. 7 Average JM distance of each crop from the other three crops in the time series of different NDVI data combinations

由图7可知,相比仅使用11期NDVISen,增加19景NDVIFUS数据后,各作物可分性有了明显提高。对比表4中各作物分类精度可以发现,小麦的分类精度受NDVI期数影响较为明显,在增加19期NDVIFUS后,小麦的分类精度与JM距离一样,都有了明显提升。且此时的小麦的可分性比在3—10月整个作物生长周期时间序列中的可分性更高,而实际分类精度也有略微领先,造成此现象的原因可能是由于小麦将在7月进行收获,在增加下半年NDVISen影像数据后,对其整个生长期而言仅7月数据为有效数据,剩下的8、9和10月NDVISen数据意义不大,甚至可能会增加冗余度。相反,增加19期NDVIFUS数据,可以在小麦收获之前补充其在生长过程中的数据缺失,增加有效信息量,因此分类精度反而更高。此外,苜蓿、玉米的可分性虽大幅度上升,但实际分类精度提升并不明显,说明仅使用11期NDVISen时,苜蓿与玉米的可分性就趋于饱和,在此基础上再增加数据可在一定程度上提高其分类精度,但效果不明显。蔬菜由于其种类多样,分类精度受影响因素较多,故不在此进行讨论。综上所述,时空融合增加可利用数据的方法在一定程度上可以提高作物可分性,对提高作物分类精度有一定效果。
在作物可分性分析基础上,利用Gini系数对S-G滤波后2019年3—6月11期Sentinel-2 NDVI和19期时空融合数据构成的时间序列进行作物识别特征重要性分析。结果如图8所示。
图8 2019年3—6月30期NDVI图像Gini系数得分情况

Fig. 8 Gini coefficient scores of 30 NDVI images from March to June in 2019

在图8中,高于平均得分线的共有10期NDVI数据,其中有6期为NDVIFUS数据,说明利用时空融合方法获取的NDVI数据在作物分类过程中占有重要的作用,对作物分类精度的提高起到了促进作用。

5 NDVI时间序列长度对作物分类精度影响

5.1 不同长度时间序列特征可分性分析

考虑到3月黑河中游地区主要作物还未种植,仅有部分蔬菜,因此舍弃单独对3月NDVI时间序列的作物可分性研究。将S-G滤波后的2019年3—6月11期NDVISen和19期NDVIFUS数据构成的时间序列进行重组,生成新的时间序列为： ① 3—4月6期NDVISen与9期NDVIFUS数据构成的时间序列;② 3—5月9期NDVISen与15期NDVIFUS数据构成的时间序列;③ 3—6月11期NDVISen与19期NDVIFUS数据构成的时间序列。利用各作物的训练样本,计算其在不同时间序列数据中的JM距离,从而对各作物的特征可分性进行评估。结果如图9所示。
图9 不同长度NDVI时间序列中每种作物与其他3种作物JM距离平均值

Fig. 9 Average JM distance between each crop and other three crops in NDVI time series of different length

如图9所示,在增加5月NDVI数据后,各作物JM距离明显增大,可分性提升显著。对比3—5月24期NDVI时间序列与3—6月30期NDVI时间序列中JM距离,可以发现,当增加6月NDVI数据后,各作物JM距离变化不大,说明仅使用3—5月的24期NDVI数据,各作物的可分性就已经趋向饱和,再增加数据对提高作物可分性作用不大。

5.2 基于不同长度时间序列的分类精度对比

分别利用上述不同NDVI时间序列,采用随机森林的分类方法对作物种植结构进行提取,研究不同时间长度的NDVI时间序列对作物分类的影响。结果如表5所示。
表5 不同长度NDVI时间序列作物分类精度对比

Tab. 5 Comparison of crop classification accuracy in NDVI time series of different lengths

作物类别 3—4月15期NDVISen 3—-5月24期NDVISen 3—6月30期NDVISen
用户精度/% 制图精度/% 用户精度/% 制图精度/% 用户精度/% 制图精度/%
苜蓿 97.49 91.40 99.45 90.15 99.48 91.26
小麦 72.57 88.27 80.10 87.93 80.97 90.46
玉米 88.49 97.56 87.78 98.63 88.67 98.44
蔬菜 66.91 28.65 82.99 41.49 80.99 43.33
总体精度/% 89.22 90.92 91.42
kappa系数 0.81 0.82 0.85
从表5可知,利用不同月份长度的NDVI图像对作物进行提取,作物总体分类精度均在89%以上,说明仅利用3—4月的15期NDVI图像便可获取精度较高的作物总体分类结果。但是此时小麦的用户精度仅为72.57%,与3—5月NDVI时间序列以及3—6月NDVI时间序列中提取出来的小麦精度存在较大差异,说明仅使用3—4月15期NDVI数据,无法获取较高精度的小麦提取结果。而使用3—5月24期NDVI影像时,总体分类精度与3—6月30期NDVI时仅相差0.5%,kappa系数仅差0.02,且各作物的分类精度差距很小,因此,综合以上信息,可以发现,当仅使用3—5月24景NDVI数据时,就可获取较高精度的作物分类结果。结合图5中可以看出,在仅使用3—4月15期NDVI影像数据对作物进行分类时,玉米的生长曲线与其他作物生长曲线相比,NDVI值明显偏低,全线小于0.2,差异明显。苜蓿则由于其NDVI值相较于小麦和蔬菜偏大,尤其是在4月30日,其NDVI值普遍大于0.5,相较于小麦和蔬菜表现出了一定的差异性,因此,也可获得较好的分类结果。而小麦的生长曲线与部分蔬菜相似,因此其分类结果精度受到了较大影响。在使用3—5月24期NDVI影像数据时,玉米生长曲线仅为单独的上升,与其余作物差异较明显。而小麦与苜蓿则都存在一个波峰,但是由于苜蓿是从4月6日开始,NDVI值出现快速增长,而小麦是从4月20日开始的,且小麦的生长曲线更为陡峭,因此可以利用这些特性,实现二者的较好区分。在3—6月30期NDVI时间序列中,玉米和小麦的生长曲线都表现出了一定的规律性,因此分类精度较高,而苜蓿生长曲线由于存在多个起伏,生长曲线特征较为明显,因此也可获取较高精度的分类结果。蔬菜生长曲线在 3个时间序列中表现都较为混乱,无明显特征,故精度较低。
在此基础上,以一旬的NDVI数据为“步长”,构建新的时间序列,对作物种植结构进行提取。由于蔬菜生长曲线无特定规律,增加或减少NDVI数据对其提取精度影响较小,故在此不予讨论。结果如表6所示。
表6 不同长度NDVI时间序列的作物精度提取结果对比

Tab. 6 Comparison of crop precision extraction results with different length NDVI time series

时间序列 玉米 苜蓿 小麦 总体精度/% Kappa
系数
用户精度/% 制图精度/% 用户精度/% 制图精度/% 用户精度/% 制图精度/%
3月—4月6日 83.07 79.08 85.78 77.91 65.80 58.20 78.72 0.72
3月—4月20日 85.98 86.89 93.42 89.68 71.23 76.13 86.12 0.79
3月—4月30日 88.49 97.56 97.49 91.40 72.57 88.27 89.22 0.81
3月—5月3日 88.33 98.62 97.29 91.39 73.97 87.39 89.72 0.82
3月—5月16日 88.52 98.66 97.13 91.79 78.67 87.33 89.94 0.82
3月—5月25日 87.78 98.63 99.45 90.15 80.10 87.93 90.92 0.82
如表6所示,当时间序列为3—4月20日苜蓿的分类精度达到一个较高水平,说明最早可在4月中旬实现对苜蓿的早期识别。而当时间序列为3—4月30日玉米的用户精度以及制图精度都达到一个较高水准,说明最早可在4月下旬实现对玉米的早期识别。小麦的分类精度受时间序列长度影响较大,当时间序列为3—5月25日小麦的分类精度均达到80%以上,说明直到5月下旬才可实现对小麦的早期识别。

6 结论

目前已有研究中,针对作物类型生长季早期识别的方法大部分为根据作物光谱特性,设置相关提取条件的“经验法”,在不同地区以及不同作物的识别上适用性较低。利用时空融合的方法,并不需要单独针对特定作物进行光谱信息的总结,可以通过生成高时间分辨率和高空间分辨率的时间序列来增加作物的可用生长信息,从而更好地实现对生长季早期的作物类型识别,且此方法在不同地区作物的早期识别中也拥有更强的适用性,因此更具实际使用价值。本文利用ESTARFM时空融合的方法,将利用Sentinel-2 计算得到的(NDVISen)数据与利用MOD09GQ产品计算得到的(NDVIMOD)数据进行融合,建立高时间分辨率和高空间分辨率的NDVI时间序列,利用随机森林分类方法进行黑河流域中游地区生长季早期的作物种植结构识别,结果表明：
（1）使用ESTARFM时空融合模型,将3-6月的NDVISen数据与NDVIMOD数据进行时空融合,获取由30景NDVI组成的时间序列,作物识别精度达到91.42%,比仅使用3-6月NDVISen数据的作物识别精度提升了1.05%,kappa系数提高0.02。与利用 3—10月的整个作物生长周期NDVISen数据的分类结果对比,发现精度仅相差1.53%,kappa系数仅差0.03,说明利用ESTARFM实现NDVISen数据与NDVIMOD数据的时空融合,可以提高时间序列的时间分辨率,有效提高作物分类精度。
（2）利用Gini系数对3—6月30期NDVI数据在分类中的重要程度从高到底排序,共有10期NDVI数据Gini系数得分高于平均值,其中有6景NDVI是时空融合得到的数据,说明时空融合获取的NDVI数据对作物提取精度起到促进作用。
（3）利用不同月份的NDVI时间序列进行作物分类发现,分类精度随着NDVI时间序列的增长而提高。其中仅使用3—4月的NDVI时间序列,可以较好地提取出苜蓿与玉米,但是小麦的提取精度较低。当使用3—5月的NDVI时间序列时,苜蓿、玉米和小麦的提取都达到较高精度,且与3—6月及3—10月NDVI时间序列的作物提取精度差异较小。在此基础上,以旬为单位,对作物进行提取,发现最早可在4月中旬和4月下旬分别实现对苜蓿和玉米的早期识别,而小麦的分类精度受时间序列长度影响较大,最早可在5月下旬实现对其的早期识别。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
欧阳玲, 毛德华, 王宗明, 等. 基于GF-1与Landsat8 OLI影像的作物种植结构与产量分析[J]. 农业工程学报, 2017, 33(11):147-156.

[ Ouyang L, Mao D H, Wang Z M. Analysis crops planting structure and yield based on GF-1 and Landsat8 OLI images[J]. Transactions of the Chinese Society of Agricultural Engineering, 2017, 33(11):147-156. ] DOI: 10.11975/j.issn.1002-6819.2017.11.019

DOI

[2]
曹卫彬, 杨邦杰, 宋金鹏. TM影像中基于光谱特征的棉花识别模型[J]. 农业工程学报, 2004, 20(4):112-116.

[ Cao W B, Yang B J, Song J P. Cotton recognition model based on spectral features in TM image[J]. Transactions of the Chinese Society of Agricultural Engineering, 2004, 20(4):112-116. ]

[3]
Ajay M, Foody G. Crop classification by support vector machine with intelligently selected training data for an operational application[J]. International Journal of Remote Sensing, 2008, 29(8):2227-2240. DOI: 10.1080/01431160701395203

DOI

[4]
胡琼, 吴文斌, 宋茜, 等. 农作物种植结构遥感提取研究进展[J]. 中国农业科学, 2015, 48(10):1900-1914.

[ Hu Q, Wu W B, Song Q, et al. Recent progresses in research of crop patterns mapping by using remote sensing[J]. Scientia Agricultura Sinica, 2015, 48(10):1900-1914. ] DOI: 10.3864/j.issn.0578-1752.2015.10.004

DOI

[5]
边增淦, 王文, 江渊. 黑河流域中游地区作物种植结构的遥感提取[J]. 地球信息科学学报, 2019, 21(10):1629-1641.

DOI

[ Bian Z G, Wang W, Jiang Y. Remote sensing of cropping structure in the middle reaches of the Heihe River Basin[J]. Journal of Geo-information Science, 2019, 21(10):1629-1641. ] DOI: 10.12082/dqxxkx.2019.190183

DOI

[6]
Zhong L H, Gong P, Biging Gregory-S. Efficient corn and soybean mapping with temporal extendability: A multi-year experiment using Landsat imagery[J]. Remote Sensing of Environment, 2013, 140:1-13. DOI: 10.1016/j.rse.2013.08.023

DOI

[7]
Cai Y P, Guan K Y, Peng J, et al. A high-performance and in-season classification system of field-level crop types using time-series Landsat data and a machine learning approach[J]. Remote Sensing of Environment, 2018, 210:35-47. DOI: 10.1016/j.rse.2018.02.045

DOI

[8]
郝鹏宇, 唐华俊, 陈仲新, 等. 基于历史增强型植被指数时序的农作物类型早期识别[J]. 农业工程学报, 2018, 34(13):179-186.

[ Hao P Y, Tang H J, Chen Z X, et al. Early season crop type recognition based on historical EVI time series[J]. Transactions of the Chinese Society of Agricultural Engineering, 2018, 34(13):179-186. ] DOI: 10.11975/j.issn.1002-6819.2018.13.021

DOI

[9]
Zhou F C, Zhang A N, Townley-Smith Lawrence. A data mining approach for evaluation of optimal time-series of MODIS data for land cover mapping at a regional level[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 84:114-129. DOI: 10.1016/j.isprsjprs.2013.07.008

DOI

[10]
Paolo Villa, Stroppiana Daniela, Fontanelli Giacomo, et al. In-Season mapping of crop type with optical and x-band SAR data: A classification tree approach using synoptic seasonal features[J]. Remote Sensing, 2015, 7(10):12859-12886. DOI: 10.3390/rs71012859

DOI

[11]
Ramin Azar, Villa Paolo, Stroppiana Daniela, et al. Assessing in-season crop classification performance using satellite data: a test case in northern Italy[J]. European Journal of Remote Sensing, 2016, 49(1):361-380. DOI: 10.5721/EuJRS20164920

DOI

[12]
Zhu X L, Gao F, Liu D S, et al. A modified neighborhood similar pixel interpolator approach for removing thick clouds in Landsat images[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(3):521-525. DOI: 10.1109/LGRS.2011.2173290

DOI

[13]
Zhu X L, Chen J, Feng G, et al. An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions[J]. Remote Sensing of Environment, 2010, 114(11):2610-2623. DOI: 10.1016/j.rse.2010.05.032

DOI

[14]
Chen J, Zhu X L, James E, Vogelmann. et al. A simple and effective method for filling gaps in Landsat ETM+ SLC-off images[J]. Remote Sensing of Environment, 2010, 115(4). DOI: 10.1016/j.rse.2010.12.010

DOI

[15]
Feng T, Wang Y J, Fensholt Rasmus, et al. Mapping and evaluation of NDVI trends from synthetic time series obtained by blending Landsat and MODIS data around a coalfield on the Loess Plateau[J]. Remote Sensing, 2013, 5(9):4255-4279. DOI: 10.3390/rs5094255

DOI

[16]
孙华生, 徐爱功, 林卉, 等. 时间序列植被指数频域滤波去噪算法的优化研究[J]. 遥感信息, 2013, 28(1):24-28.

[ Sun H S, Xu A G, Lin H, et al. Optimization of frequency domain denoising algorithms for time-series vegetation index[J]. Remote sensing Information, 2013, 28(1):24-28 ] DOI: 10.3969/j.issn.1000-3177.2013.01.006

DOI

[17]
程琳琳, 李玉虎, 孙海元, 等. 京津冀MODIS长时序增强型植被指数拟合重建方法适用性研究[J]. 农业工程学报, 2019, 35(11):148-158.

[ Cheng L L, Li Y H, Sun H Y, et al. Applicability of fitting and reconstruction method of MODIS long-time enhanced vegetation index in Beijing-Tianjin-Hebei[J]. Transactions of the Chinese Society of Agricultural Engineering, 2019, 35(11):148-158. ] DOI: 10.11975/j.issn.1002-6819.2019.11.017

DOI

[18]
栗云峰. 基于时间序列高分一号影像的南京市农业用地提取方法研究[D]. 南京:南京大学, 2019.

[ Li Y F. Agricultural land extraction in Nanjing based on time series GF-1 WFV images[D]. Nanjing: Nanjing University, 2019. ]

[19]
V-F Rodriguez-Galiano, Ghimire B, Rogan J, et al. An assessment of the effectiveness of a random forest classifier for land-cover classification[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 67:93-104. DOI: 10.1016/j.isprsjprs.2011.11.002

DOI

[20]
杨振忠. 基于遥感数据的水稻生育期识别研究[D]. 武汉:武汉大学, 2019.

[ Yang Z Z. Identification of growth stage in rice using remote sensing date[D]. Wuhan: Wuhan University, 2019. ]

Options
Outlines

/

〈 〉