CN114724024B - 基于云计算平台和生命周期的双季作物种植界线自动化提取方法 - Google Patents

Abstract

本发明属于遥感目标识别技术领域，公开一种基于云计算平台和生命周期的双季作物种植界线自动化提取方法，该方法首先利用云计算平台收集所有可用的遥感数据并进行预处理，获取研究区多年时间序列数据集；然后根据研究区植被的全生命周期曲线，估算作物种植或收割次数，制作双季作物种植图；最后利用核密度函数与等值线法提取不同作物种植强度变化界线。本发明能够充分利用植被全生命周期的变化量化植被年内生长周期数，表征不同土地集约化利用强度，该发明还可以用于类似作物种植强度变化界线的提取，如牧草收割次数、放牧强度以及人工林地砍伐与再种植强度，揭示土地集约化利用程度，为农牧业发展规划提供科学依据。

Description

基于云计算平台和生命周期的双季作物种植界线自动化提取方法

技术领域

本发明属于遥感目标识别技术领域，尤其涉及一种基于云计算平台和生命周期的双季作物种植界线自动化提取方法。

背景技术

在中国，一半以上的耕地实施双季种植。双季种植作为中国重要的种植制度，其面积在全球气候变化、农业政策和耕作方式变化的影响下发生了重大变化。双季作物种植北界是双季作物种植制度的空间分布北部界限，同时也是双季作物种植时空格局动态变化的重要表征。明晰双季作物种植北界的时空变迁特征，对于评估粮食生产安全，缩小粮食需求缺口和改善生态系统和人类健康具有重大意义。

双季作物种植北界可为以气象观测数据为基础的潜在北界和以遥感数据为基础的实际北界。目前对于双季作物种植北界的研究主要集中在潜在北界上，而缺少实际北界及其多年变动的研究。但由于潜在北界是基于当地的实际气候条件呈现出的最佳种植模式，然而农民更倾向于按照传统模式(过去的气候条件和经验)去进行种植。相较于实际北界，潜在北界并不能反映区域真实的耕作制度。实际北界的提取是研究其多年变动的前提。

发明内容

本发明针对现有的遥感识别技术中难以准确量化不同作物种植强度界线的技术问题，提出一种基于云计算平台和生命周期的双季作物种植界线自动化提取方法，该方法基于时空大数据平台识别作物全生命周期，计算作物种植次数，获取双季作物种植范围，耦合核密度函数和等值线法，提取作物种植北界。

为了实现上述目的，本发明采用以下技术方案：

一种基于云计算平台和生命周期的双季作物种植界线自动化提取方法，包括：

步骤1、基于Google Earth Engine云计算平台获取研究期内研究区域中所有MODIS卫星影像，并对MODIS卫星影像进行预处理，得到预处理后的卫星影像数据集；

步骤2、基于步骤1所述的卫星影像数据集，获取归一化植被指数NDVI时间序列和地表水分指数LSWI时间序列；

步骤3、对步骤2得到的NDVI时间序列进行处理，得到最大NDVI时间序列数据集；

步骤4、基于步骤3中得到的NDVI时间序列，识别NDVI时间序列曲线的波峰和波谷，得到所有波峰值及波峰总个数，并得到所有波谷值及对应的日期DOY；

步骤5、对步骤2得到的LSWI时间序列进行处理，得到平均LSWI时间序列数据集；

步骤6、根据步骤5中得到的LSWI时间序列逐像素确定LSWI阈值，得到每个像素的LSWI阈值；

步骤7、根据步骤5中得到的LSWI时间序列以及步骤6中得到的LSWI阈值识别裸土，得到裸土分布图；

步骤8、根据步骤4中得到的波峰值、波峰总个数和步骤7中得到裸土分布图确定有效峰值个数，得到峰值个数图；

步骤9、根据步骤8中得到峰值个数图确定双季作物种植像素，得到双季作物种植分布图；

步骤10、根据步骤9中得到的双季作物种植分布图，利用核密度估计法获取双季作物种植北界范围，得到双季作物种植北界范围；

步骤11、根据步骤10中得到的双季作物种植北界范围，利用等值线提取法提取双季作物种植北界线，得到双季作物种植北界图。

进一步地，所述步骤1中，对MODIS卫星影像进行预处理：

利用FMask算法对所述卫星影像进行观测值提取，去除云、云阴影、卷云和冰/雪覆盖的观测值。

进一步地，所述步骤3中，按照如下方式对NDVI时间序列进行处理：

计算每8天NDVI的最大值作为综合NDVI观测值，获得等时间间隔时间序列NDVI数据集；在缺失8天的综合NDVI观测值的地区，根据8天前、后的综合NDVI观测值进行线性插值；使用Savitzky-Golay滤波器对NDVI数据集进行平滑处理。

进一步地，所述步骤4中，按照如下方式识别NDVI时间序列曲线的波峰和波谷：

识别NDVI时间序列中的局部最大值为波峰，识别NDVI时间序列中的局部最小值为波谷，如果某一时刻的NDVI值高于该时刻前后的NDVI值，则定义为一个波峰，如果某一时刻的NDVI值低于该时刻前后的NDVI值，则定义为一个波谷。

进一步地，所述步骤5中，按照如下方式对LSWI时间序列进行处理：

计算每8天LSWI的平均值作为综合LSWI观测值，获得等时间间隔时间序列LSWI数据集；在缺失8天的综合LSWI观测值的地区，根据8天前、后的综合LSWI观测值进行线性插值。

进一步地，所述步骤6中，按照如下方式确定LSWI阈值：

T_P＝LSWI_min+(LSWI_max-LSWI_min)×0.15

其中T_P是潜在LSWI阈值，T_LSWI是最终LSWI阈值，LSWI_min和LSWI_max分别是研究期内LSWI的最小值和最大值。

进一步地，所述步骤7中，按照如下方式识别裸土：

比较LSWI时间序列中各LSWI值与LSWI阈值的大小，若LSWI值小于LSWI阈值，则被分类为裸土。

进一步地，所述步骤8中，按照如下方式确定有效峰值个数：

对于NDVI时间序列的波峰，如果波峰值大于0.5，且波峰相邻的两个波谷时期均检测出裸土，则该峰值被识别为有效峰值。

进一步地，所述步骤10中，按照如下方式获取双季作物种植北界范围：

将所获的双季作物种植分布图由栅格转换为矢量点，计算概率密度的估计量：

其中，f_h(x)是概率密度的估计量，n是双季作物种植点的数量，h是用户定义的平滑参数或带宽，x是估计点的矢量坐标，x_i是采样点的矢量坐标；K是用户定义的非负数内核函数，它被视为二次Epanechnikov内核，其定义如下：

进一步地，所述步骤11中，按照如下方式提取双季作物种植北界线：

按照选定的最佳带宽，使用百分位数方法确定具有95％百分位数的核密度估计阈值，并提取等于该阈值的等值线，选取位于最北的等值线作为双季作物种植北界。

与现有技术相比，本发明具有的有益效果：

(1)本发明充分利用了植被生命周期的变化特征来量化植被年内生长周期数，表征了不同土地种植作物强度，进而提取不同作物利用强度的界线，本发明提取了双季作物种植北界，生成的双季作物种植北界地图对评估中国北方农业生产的潜力以及农业发展规划具有重大意义；

(2)本发明可以识别其他年份和其他地区的双季作物种植农田并提取其种植北界，为种植强度以及作物种植界线的研究提供了一种新的实证研究方法，该方法也可用于其他类似作物种植季节的情景，如牧草收割次数与放牧强度以及人工林砍伐与再种植强度；

(3)本发明首次提出了一种双季作物种植北界划分方法，利用该方法实现了基于云计算平台和生命周期的双季作物种植界线自动化提取。

附图说明

图1为本发明实施例一种基于云计算平台和生命周期的双季作物种植界线自动化提取方法的流程图；

图2为本发明实施例生成的双季作物种植分布图；

图3为本发明实施例生成的双季作物种植北界图。

具体实施方式

下面结合附图和具体的实施例对本发明做进一步的解释说明：

如图1所示，一种基于云计算平台和生命周期的双季作物种植界线自动化提取方法，包括：

S1、基于Google Earth Engine云平台获取研究期内研究区中所有MODIS卫星影像，并对MODIS卫星影像进行预处理得到高质量的卫星影像数据集；作为一种可实施方式，以2000-2020年为研究期，以北京，天津，河北，陕西，山西，山东，河南和湖北在内的八个省份为研究区。

进一步地，对MODIS卫星影像进行预处理包括：

利用FMask算法对所述卫星影像进行观测值提取，去除云、云阴影、卷云和冰/雪覆盖的观测值。

S2、基于步骤S1所述的卫星影像数据集，获取归一化植被指数NDVI时间序列和地表水分指数LSWI时间序列；

具体地，NDVI和LSWI计算公式分别为：

其中ρ_NIR、ρ_RED和ρ_SWIR分别代表卫星影像数据中的近红外波段、红色波段和短波红外波段。

S3、对步骤S2得到的NDVI时间序列数据集进行处理，得到最大NDVI时间序列数据集；

进一步地，按照如下方式对NDVI时间序列进行处理：

首先计算每8天内所有NDVI的最大值作为综合NDVI观测值，获得等时间间隔时间序列NDVI数据集；其次基于8天前、后的综合NDVI观测值对空隙进行线性插值，最后使用Savitzky-Golay滤波器(S-G滤波器)平滑NDVI时间序列数据集。

S4、根据步骤S3中得到的NDVI时间序列，识别NDVI时间序列曲线的波峰和波谷，得到所有波峰值及波峰总个数，并得到所有波谷值及对应的日期DOY；

具体地，识别NDVI时间序列中的局部最大值为波峰，识别NDVI时间序列中的局部最小值为波谷，如果某一时刻的NDVI值高于该时刻前后的NDVI值，则定义为一个波峰，如果某一时刻的NDVI值低于该时刻前后的NDVI值，则定义为一个波谷。

S5、对步骤S2得到的LSWI时间序列进行处理，得到平均LSWI时间序列数据集；

进一步地，按照如下方式对LSWI时间序列进行处理：

计算每8天LSWI的平均值作为综合LSWI观测值，获得等时间间隔时间序列LSWI数据集；在缺失8天的综合LSWI观测值的地区，根据8天前、后的综合LSWI观测值进行线性插值。由于LSWI对水分状况较为敏感，因此并未对LSWI时间序列进行平滑。

S6、根据步骤S5中得到的LSWI时间序列逐像素确定LSWI阈值，得到每个像素的LSWI阈值；

进一步地，确定LSWI阈值的方法为：

T_P＝LSWI_min+(LSWI_max-LSWI_min)×0.15

其中T_P是潜在LSWI阈值，T_LSWI是最终LSWI阈值，LSWI_min和LSWI_max分别是研究期内LSWI的最小值和最大值。我们用最终的LSWI阈值(LSWI<T_LSWI)来确定裸土的日期。

S7、根据步骤S5中得到的LSWI时间序列以及步骤S6中得到的LSWI阈值识别裸土，得到裸土分布图；

具体地，比较LSWI时间序列中各LSWI值与LSWI阈值的大小，如果某一时刻的LSWI值小于LSWI阈值，则被分类为裸土。

S8、根据步骤S4中得到的波峰值及总个数和步骤S7中得到裸土分布图确定有效峰值个数，得到峰值个数图；

在作物种植前和收获后，耕地像素通常被裸土或作物秸秆覆盖，叶绿素含量急剧降低，因此NDVI值通常较低。然而，一些特殊作物(如越冬作物)在一个生长期可能出现多个波峰，且峰值均较高。NDVI曲线的噪声通常也会导致多个波峰的出现。因此，为了准确识别作物的生长周期，我们使用LSWI指数来判断生长周期的开始和结束，因为裸土或作物秸秆的LSWI值远小于作物。

具体地，对于NDVI时间序列的波峰，如果波峰值大于0.5，且波峰相邻的两个波谷时期均检测出裸土，则该峰值被识别为有效峰值。

S9、根据步骤S8中得到峰值个数图确定双季作物种植像素，得到双季作物种植分布图；

具体地，有效峰值个数等于2的像素被分类为双季作物种植。

S10、根据步骤S9中得到的双季作物种植分布图，利用核密度估计法获取双季作物种植北界范围，得到双季作物种植北界范围；

具体地，将双季作物种植强度栅格图转为矢量点，作为核密度估计法的输入参数。其次，利用核密度估计法，按照选定的最佳带宽进行分析处理，得到双季作物种植密度空间分布图。其中，在选择最佳带宽时，我们对比了3种不同的带宽(5km，10km，15km)。5km的带宽曲线过于破碎，难以反映双季种植核心区情况，而15km带宽曲线过于平滑，忽视了双季种植稀疏区和非双季种植区的分布。因此，本发明确定双季作物种植北界提取的最佳带宽为10km。核密度估计法的计算公式为：

其中，f_h(x)是概率密度的估计量，n是双季作物种植点的数量，h(h>0)是用户定义的平滑参数或带宽，x是估计点的矢量坐标，x_i是采样点的矢量坐标。K是用户定义的非负数内核函数，它被视为二次Epanechnikov内核，其定义如下：

S11、根据步骤S10中得到的双季作物种植北界范围，利用等值线提取法提取双季作物种植北界线，得到双季作物种植北界图。

具体地，按照选定的最佳带宽，使用百分位数方法确定具有95％百分位数的核密度估计阈值，并提取等于该阈值的等值线，选取位于最北的等值线作为双季作物种植北界。

为验证本发明效果，通过本发明方法，我们生成了2001、2005、2010、2015和2020年的双季作物种植分布地图，如图2所示。在2001，2005，2010，2015和2020年分别有35.74％，39.16％，35.31％，37.85％和31.99％的农田为双季作物种植区。双季作物种植的空间分布与地形特征表现出较强的一致性。其中，河北省，河南省和山东省地势平坦，是双季作物种植的主要分布省份，占研究区双季作物种植区一半以上。通过本发明方法，我们生成了2001、2005、2010、2015和2020年的双季作物种植北界地图，如图3所示。各时期双季作物种植北界主要在北京市，天津市，山西省，河北省，陕西省和河南省区域内发生变动，整体上呈现东部高纬度向西部低纬度阶级降低的态势，在空间格局变化上呈现波动趋势。

具体地，通过实地调查和谷歌影像目视解译来获取地面参考数据作为验证样本。2020年3月至2020年6月期间，在研究区进行了6次实地考察并收集了不同种植强度的地理照片。这些图片包括双季作物种植，单季种植和两年三季等种植模式。其次，结合获取的实地照片，在Google Earth上标记与实地照片位置影像颜色，纹理相似的地块，并添加相应属性(单季，双季，两年三季等)。第三，在实地考察期间，小型无人机系统被用来获取不同类型地块的多光谱影像，以供目视解译。根据获得的地面参考数据，将已知的地块类型数字化为多边形。最终分别收集了284个(993像素)双季样本和73个(1964像素)非双季样本。进一步，我们使用上述样本计算了混淆矩阵，结果如表1所示，2020年双季作物种植图的总体精度，用户精度，生产者精度分别为95.97％，96.58％和92.21％，显示出较高的分类精度。0.91的Kappa系数表明分类结果与地面参考数据具有较强的一致性。总体而言，本次研究的精度和Kappa系数较高，证明了本发明的有效性、可靠性和科学性。

表1

综上，本发明以2000-2020年为研究期，以北京，天津，河北，陕西，山西，山东，河南和湖北在内的八个省份为研究区，提供了一种可用于提取其他地区，其他年份，其他地物种类种植界限的研究思路，作物生长周期的确定是基于以下简单的原则：在播种前和收获后必须存在裸露的土壤或作物秸秆，这可以通过LSWI值来确定。在其他地区，由于不同的环境因素，例如土壤湿度，播种方式，以及作物品种的不同，双季作物种植的物候特征可能存在差异，因此可以通过修正阈值来使模型更贴切研究区域的真实情况。在连续的年份里，区域气候条件以及农民的耕种习惯变动极小，因此使用本发明可以提取不同年份的双季作物种植北界。这种双季作物种植北界的提取方法也可用于其他作物种植界限的提取，例如水稻、小麦、和玉米等。

此外，生成的多年双季作物种植北界地图对于指导中国北方粮食生产具有指导意义。通过对比多年的双季作物种植北界地图能够发现种植强度发生变化的地区，这对于了解影响种植强度变化的因素，科学有效地提高耕地集约化利用程度以及实现粮食增产具有积极作用。

目前有关于双季作物种植北界提取的研究较为缺乏，主要存在两方面的缺陷。一方面，较粗分辨率的数据(8km的NOAA、1km SPOT和MODIS)存在大量混合像元，基于这些数据的产品的准确性容易受到高的亚像素异质性限制，尤其是在以小农农场为主要特征的中国农田。另一方面，诸如邻域比较法，二次差分法和滑动分割法等传统的方法没有考虑到作物与自然植被明显不同的物候信息以及种植强度的原位信息，且主要集中在对种植强度总体分布的监测上，对种植北界及其变化的定量监测十分有限。综合考虑各种数据的时间分辨率，空间分辨率以及可用数据的年限，本发明使用了MODIS 500m数据，并充分利用了植被生命周期的变化量化植被年内和年际生长周期数，能够表征不同土地集约化利用强度，如作物种植强度、牧草收割次数与放牧强度以及人工林地砍伐与再种植强度，揭示土地集约化利用程度，为农牧业发展规划提供科学依据。这相较于其他算法，本发明建立的规则更符合作物的真实生长规律，更能准确识别双季作物种植像素。

以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (7)

1.一种基于云计算平台和生命周期的双季作物种植界线自动化提取方法，其特征在于，包括：

步骤1、基于Google Earth Engine云计算平台获取研究期内研究区域中所有MODIS卫星影像，并对MODIS卫星影像进行预处理，得到预处理后的卫星影像数据集；

步骤2、基于步骤1所述的卫星影像数据集，获取归一化植被指数NDVI时间序列和地表水分指数LSWI时间序列；

步骤3、对步骤2得到的NDVI时间序列进行处理，得到最大NDVI时间序列数据集；

步骤4、基于步骤3中得到的NDVI时间序列，识别NDVI时间序列曲线的波峰和波谷，得到所有波峰值及波峰总个数，并得到所有波谷值及对应的日期DOY；

步骤5、对步骤2得到的LSWI时间序列进行处理，得到平均LSWI时间序列数据集；

步骤6、根据步骤5中得到的LSWI时间序列逐像素确定LSWI阈值，得到每个像素的LSWI阈值；

步骤7、根据步骤5中得到的LSWI时间序列以及步骤6中得到的LSWI阈值识别裸土，得到裸土分布图；

步骤8、根据步骤4中得到的波峰值、波峰总个数和步骤7中得到裸土分布图确定有效峰值个数，得到峰值个数图；所述步骤8中，按照如下方式确定有效峰值个数：

对于NDVI时间序列的波峰，如果波峰值大于0.5，且波峰相邻的两个波谷时期均检测出裸土，则该峰值被识别为有效峰值；

步骤9、根据步骤8中得到峰值个数图确定双季作物种植像素，得到双季作物种植分布图；

步骤10、根据步骤9中得到的双季作物种植分布图，利用核密度估计法获取双季作物种植北界范围，得到双季作物种植北界范围；所述步骤10中，按照如下方式获取双季作物种植北界范围：

将所获的双季作物种植分布图由栅格转换为矢量点，计算概率密度的估计量：

其中，f_h(x)是概率密度的估计量，n是双季作物种植点的数量，h是用户定义的平滑参数或带宽，x是估计点的矢量坐标，x_i是采样点的矢量坐标；K是用户定义的非负数内核函数，它被视为二次Epanechnikov内核，其定义如下：

步骤11、根据步骤10中得到的双季作物种植北界范围，利用等值线提取法提取双季作物种植北界线，得到双季作物种植北界图；

所述步骤11中，按照如下方式提取双季作物种植北界线：

按照选定的最佳带宽，使用百分位数方法确定具有95％百分位数的核密度估计阈值，并提取等于该阈值的等值线，选取位于最北的等值线作为双季作物种植北界。

2.根据权利要求1所述的基于云计算平台和生命周期的双季作物种植界线自动化提取方法，其特征在于，所述步骤1中，对MODIS卫星影像进行预处理：

利用FMask算法对所述卫星影像进行观测值提取，去除云、云阴影、卷云和冰/雪覆盖的观测值。

3.根据权利要求1所述的基于云计算平台和生命周期的双季作物种植界线自动化提取方法，其特征在于，所述步骤3中，按照如下方式对NDVI时间序列进行处理：

计算每8天NDVI的最大值作为综合NDVI观测值，获得等时间间隔时间序列NDVI数据集；在缺失8天的综合NDVI观测值的地区，根据8天前、后的综合NDVI观测值进行线性插值；使用Savitzky-Golay滤波器对NDVI数据集进行平滑处理。

4.根据权利要求1所述的基于云计算平台和生命周期的双季作物种植界线自动化提取方法，其特征在于，所述步骤4中，按照如下方式识别NDVI时间序列曲线的波峰和波谷：

识别NDVI时间序列中的局部最大值为波峰，识别NDVI时间序列中的局部最小值为波谷，如果某一时刻的NDVI值高于该时刻前后的NDVI值，则定义为一个波峰，如果某一时刻的NDVI值低于该时刻前后的NDVI值，则定义为一个波谷。

5.根据权利要求1所述的基于云计算平台和生命周期的双季作物种植界线自动化提取方法，其特征在于，所述步骤5中，按照如下方式对LSWI时间序列进行处理：

计算每8天LSWI的平均值作为综合LSWI观测值，获得等时间间隔时间序列LSWI数据集；在缺失8天的综合LSWI观测值的地区，根据8天前、后的综合LSWI观测值进行线性插值。

6.根据权利要求1所述的基于云计算平台和生命周期的双季作物种植界线自动化提取方法，其特征在于，所述步骤6中，按照如下方式确定LSWI阈值：

T_P＝LSWI_min+(LSWI_max-LSWI_min)×0.15

其中T_P是潜在LSWI阈值，T_LSWI是最终LSWI阈值，LSWI_min和LSWI_max分别是研究期内LSWI的最小值和最大值。

7.根据权利要求6所述的基于云计算平台和生命周期的双季作物种植界线自动化提取方法，其特征在于，所述步骤7中，按照如下方式识别裸土：

比较LSWI时间序列中各LSWI值与LSWI阈值的大小，若LSWI值小于LSWI阈值，则被分类为裸土。

Images