基于植被指数时序谱特征的大豆种植区提取方法和系统

1.本发明属于遥感影像处理技术领域，尤其涉及基于植被指数时序谱特征的大豆种植区提取方法和系统。


背景技术：

2.目前大豆种植区提取方法主要基于中高空间分辨率的遥感图像分类方法和基于中低空间分辨率的时序数据提取方法，其中基于中高空间分辨率的遥感图像分类方法一般选择几景影像，会因为缺少对时序特征的充分利用而限制分类精度。而基于中低空间分辨率时序数据的提取方法又因较低空间分辨率产生的混合像元问题而影响提取精度。随着大量中高分辨率的遥感卫星发射升空，使得利用中高空间分辨率高频度时序数据成为可能。例如，sentinel-2a/b卫星的空间分辨率为10米(蓝色、绿色、红色和近红外波段)和20米(红边1、红边2、红边3、红色边缘4、短波红外1和短波红外2波段)，它的一颗卫星的重访周期为10天，两颗互补，重访周期为5天，相对较短的重访周期可以提供作物更详细物候信息。此外，红边范围含有三个波段，这可能有助于区分形态相似作物类型之间的细微差异。因此，亟需充分挖掘中高空间分辨率高频度时序数据优势，构建大豆种植区高精度提取方法。
3.专利：一种大豆生长季空间分布图的生成方法和系统，大豆生长季的影像的提取和预处理；随机森林分类模型的构建与训练；时间窗口的设置；特征子集的选择；大豆生长季空间分布图的获得。本技术构造大豆生长季内sentinel-2光谱波段的时间序列合成影像，然后结合随机森林分类模型探究大豆最早识别的时间窗口，其次通过评估时间窗口内所有特征的重要性进一步筛选特征子集，最终绘制出大豆空间分布图。


技术实现要素：

4.为解决上述技术问题，本发明提出基于植被指数时序谱特征的大豆种植区提取方法和系统的技术方案，以解决上述技术问题。
5.本发明第一方面公开了一种基于植被指数时序谱特征的大豆种植区提取方法，所述方法包括：
6.步骤s1、遥感数据预处理：对sentinel-2卫星影像进行预处理，并进行耕地掩膜；
7.步骤s2、特征提取：提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征，简称时序谱特征；提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征，简称光谱特征；提取海拔及坡度特征，简称地形特征，并利用耕地数据对所述时序谱特征、光谱特征和地形特征分别进行掩膜；
8.步骤s3、分类模型的构建与训练：构建分类模型并通过训练设置参数；
9.步骤s4、以所述时序谱特征为主，将所述时序谱特征、光谱特征和地形特征进行组合得到待选特征组合，得到待选特征组合；将所述待选特征组合输入分类模型，计算待选特征组合的分类精度，根据分类精度得到最优特征组合；
10.步骤s5、大豆种植区提取及精度评价：将最优特征组合输入随机森林分类模型提
取大豆种植区分布图并进行精度评价。
11.根据本发明第一方面的方法在所述步骤s2中，所述提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征的方法包括：
12.利用线性谐波模型提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征；
13.所述线性谐波模型的公式为：
[0014][0015]
其中，f(t)是第t时刻的拟合的增强型植被指数的值，ak是余弦系数，bk是正弦系数，c是截距项，n是谐波级数的阶数，w是频率等于1.5；自变量t是一年中的某一天；n分别设置为1至5，每次增加1，依据大豆的时间序列中原始值与拟合值之间均方误差最小确定n的值，此处n取2；取相位和振幅作为时序特征输入，振幅定义为二维矢量[ak,bk]的长度，相位定义为二维矢量[ak,bk]形成的角度。
[0016]
根据本发明第一方面的方法在所述步骤s2中，所述提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征的方法包括：
[0017]
利用百分位数、最大值、最小值、平均值和标准差方法提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征。
[0018]
根据本发明第一方面的方法在所述步骤s2中，大豆关键生育期内掩膜后的sentinel-2影像波段包含，红边1、红边2、红边3、短波红外1和短波红外2；
[0019]
所述关键生育期为1-2月，或者为8-9月，具体包含大豆开花期、结荚期和灌浆期三个生长阶段。
[0020]
根据本发明第一方面的方法在所述步骤s2中，所述提取海拔及坡度特征的方法包括：
[0021]
利用航天飞机雷达地形测绘任务数字高程数据提取海拔及坡度特征。
[0022]
根据本发明第一方面的方法在所述步骤s4中，所述待选特征组合的具体组合方式包括：
[0023]
时序谱特征、光谱特征、地形特征、时序谱特征及光谱特征、时序谱特征及地形特征、时序谱特征和光谱特征以及地形特征，共6组待选特征组合。
[0024]
根据本发明第一方面的方法在所述步骤s4中，计算待选特征组合的分类精度的方法包括：
[0025]
采用基于混淆矩阵的随机森林分类精度评价方法计算待选特征组合的分类精度。
[0026]
本发明第二方面公开了一种基于植被指数时序谱特征的大豆种植区提取系统，所述系统包括：
[0027]
第一处理模块，被配置为，遥感数据预处理：对sentinel-2卫星影像进行预处理，并进行耕地掩膜；
[0028]
第二处理模块，被配置为，特征提取：提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征，简称时序谱特征；提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征，简称光谱特征；提取海拔及坡度特征，简称地形特征，并利用耕地数据对所述时序谱特征、光谱特征和地形特征分别进行掩膜；
[0029]
第三处理模块，被配置为，分类模型的构建与训练：构建分类模型并通过训练设置
参数；
[0030]
第四处理模块，被配置为，以所述时序谱特征为主，将所述时序谱特征、光谱特征和地形特征进行组合得到待选特征组合，得到待选特征组合；将所述待选特征组合输入分类模型，计算待选特征组合的分类精度，根据分类精度得到最优特征组合；
[0031]
第五处理模块，被配置为，大豆种植区提取及精度评价：将最优特征组合输入随机森林分类模型提取大豆种植区分布图并进行精度评价。
[0032]
根据本发明第二方面的系统，第二处理模块，被配置为，所述提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征包括：
[0033]
利用线性谐波模型提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征；
[0034]
所述线性谐波模型的公式为：
[0035][0036]
其中，f(t)是第t时刻的拟合的增强型植被指数的值，ak是余弦系数，bk是正弦系数，c是截距项，n是谐波级数的阶数，w是频率等于1.5；自变量t是一年中的某一天；n分别设置为1至5，每次增加1，依据大豆的时间序列中原始值与拟合值之间均方误差最小确定n的值，此处n取2；取相位和振幅作为时序特征输入，振幅定义为二维矢量[ak,bk]的长度，相位定义为二维矢量[ak,bk]形成的角度。
[0037]
根据本发明第二方面的系统，第二处理模块，被配置为，所述提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征包括：
[0038]
利用百分位数、最大值、最小值、平均值和标准差方法提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征。
[0039]
根据本发明第一方面的方法在所述步骤s2中，大豆关键生育期内掩膜后的sentinel-2影像波段包含，红边1、红边2、红边3、短波红外1和短波红外2；
[0040]
所述关键生育期为1-2月，或者为8-9月，具体包含大豆开花期、结荚期和灌浆期三个生长阶段。
[0041]
根据本发明第二方面的系统，第二处理模块，被配置为，所述提取海拔及坡度特征包括：
[0042]
利用航天飞机雷达地形测绘任务数字高程数据提取海拔及坡度特征。
[0043]
根据本发明第二方面的系统，第四处理模块，被配置为，所述待选特征组合的具体组合方式包括：
[0044]
时序谱特征、光谱特征、地形特征、时序谱特征及光谱特征、时序谱特征及地形特征、时序谱特征和光谱特征以及地形特征，共6组待选特征组合。
[0045]
根据本发明第二方面的系统，第四处理模块，被配置为，计算待选特征组合的分类精度包括：
[0046]
采用基于混淆矩阵的随机森林分类精度评价方法计算待选特征组合的分类精度。
[0047]
本发明第三方面公开了一种电子设备。电子设备包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时，实现本公开第一方面中任一项的一种基于植被指数时序谱特征的大豆种植区提取方法中的步骤。
[0048]
本发明第四方面公开了一种计算机可读存储介质。计算机可读存储介质上存储有
计算机程序，计算机程序被处理器执行时，实现本公开第一方面中任一项的一种基于植被指数时序谱特征的大豆种植区提取方法中的步骤。
[0049]
本发明提出的方案，提供了大区域尺度大豆种植区快速准确的识别方法，为大豆面积、产量、病虫害监测等提供科学依据。
附图说明
[0050]
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0051]
图1为根据本发明实施例的一种基于植被指数时序谱特征的大豆种植区提取方法的流程图；
[0052]
图2为根据本发明实施例的2020/2021生长季巴西大豆种植区分布图；
[0053]
图3为根据本发明实施例的一种基于植被指数时序谱特征的大豆种植区提取系统的结构图；
[0054]
图4为根据本发明实施例的一种电子设备的结构图。
具体实施方式
[0055]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例只是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0056]
本发明第一方面公开了一种基于植被指数时序谱特征的大豆种植区提取方法。图1为根据本发明实施例的一种基于植被指数时序谱特征的大豆种植区提取方法的流程图，如图1所示，所述方法包括：
[0057]
步骤s1、遥感数据预处理：对sentinel-2卫星影像进行预处理，并进行耕地掩膜；
[0058]
步骤s2、特征提取：提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征，简称时序谱特征；提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征，简称光谱特征；提取海拔及坡度特征，简称地形特征，并利用耕地数据对所述时序谱特征、光谱特征和地形特征分别进行掩膜；
[0059]
步骤s3、分类模型的构建与训练：构建分类模型并通过训练设置参数；
[0060]
步骤s4、以所述时序谱特征为主，将所述时序谱特征、光谱特征和地形特征进行组合得到待选特征组合，得到待选特征组合；将所述待选特征组合输入分类模型，计算待选特征组合的分类精度，根据分类精度得到最优特征组合；
[0061]
步骤s5、大豆种植区提取及精度评价：将最优特征组合输入随机森林分类模型提取大豆种植区分布图并进行精度评价。
[0062]
在步骤s1，遥感数据预处理：对sentinel-2卫星影像进行预处理，并进行耕地掩膜。
[0063]
具体地，1)在gee云计算平台中调用大豆整个生育期内所有的sentinel-2大气底
端反射率数据，利用质量控制qa60波段将云量限制在10％以内，并将所有波段值除以10000得到各个光谱波段反射率值，获得大豆整个生育期内预处理后的sentinel-2影像；
[0064]
此处，大豆整个生育期指大豆从播种至收获期结束。
[0065]
2)在gee平台中调用欧洲航天局(esa)2020年空间分辨率为10m的全球土地覆盖数据，利用ee.filter.eq()筛选器获得耕地数据(此处耕地数据的标签为40，故筛选器参数设置为eq(40))。
[0066]
在步骤s2，特征提取：提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征，简称时序谱特征；提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征，简称光谱特征；提取海拔及坡度特征，简称地形特征，并利用耕地数据对所述时序谱特征、光谱特征和地形特征分别进行掩膜。
[0067]
在一些实施例中，在所述步骤s2中，所述提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征的方法包括：
[0068]
利用线性谐波模型提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征；
[0069]
所述线性谐波模型的公式为：
[0070][0071]
其中，f(t)是第t时刻的拟合的增强型植被指数的值，ak是余弦系数，bk是正弦系数，c是截距项，n是谐波级数的阶数，ω是频率等于1.5；自变量t是一年中的某一天；n分别设置为1至5，每次增加1，依据大豆的时间序列中原始值与拟合值之间均方误差最小确定n的值，此处n取2；取相位和振幅作为时序特征输入，振幅定义为二维矢量[ak,bk]的长度，相位定义为二维矢量[ak,bk]形成的角度。
[0072]
所述提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征的方法包括：
[0073]
利用百分位数、最大值、最小值、平均值和标准差方法提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征。
[0074]
大豆关键生育期内掩膜后的sentinel-2影像波段包含，红边1、红边2、红边3、短波红外1和短波红外2；
[0075]
所述关键生育期为1-2月，或者为8-9月，具体包含大豆开花期、结荚期和灌浆期三个生长阶段。
[0076]
所述提取海拔及坡度特征的方法包括：
[0077]
利用航天飞机雷达地形测绘任务数字高程数据提取海拔及坡度特征。
[0078]
具体地，表1为不同类别特征信息简介。
[0079]
表1
[0080][0081]
1)构建线性谐波模型提取大豆整个生育期内预处理后的sentinel-2影像的增强型植被指数(evi)时序谱特征(简称时序谱特征)，再利用mask函数对所提时序谱特征进行耕地掩膜；
[0082]
此处，大豆整个生育期内预处理后的sentinel-2影像的增强型植被指数(evi)的计算公式为：
[0083][0084]
其中，red、blue、nir分别是大豆整个生育期内预处理后的sentinel-2影像的红色、蓝色和近红外波段值，数值2.5为增益因子的值，6、7.5和1对应校正气溶胶对冠层反射率和土壤背景信号的影响的系数。
[0085]
所述线性谐波模型的公式为：
[0086][0087]
其中，f(t)是第t时刻的拟合的增强型植被指数的值，ak是余弦系数，bk是正弦系数，c是截距项，n是谐波级数的阶数，ω是频率等于1.5；自变量t是一年中的某一天；n分别设置为1至5，每次增加1，依据大豆的时间序列中原始值与拟合值之间均方误差最小确定n的值，此处n取2；取相位和振幅作为时序特征输入，振幅定义为二维矢量[ak,bk]的长度，相位定义为二维矢量[ak,bk]形成的角度。
[0088]
2)利用百分位数(15/50/90)、最大值、最小值、平均值、标准差方法提取大豆关键生育期内预处理后的sentinel-2影像的5个波段反射率特征(简称光谱特征)，再利用mask函数对所提光谱特征进行耕地掩膜；
[0089]
此处，大豆关键生育期指大豆开花期、结荚期至灌浆期三个生长阶段。在大豆整个生育期内预处理后的sentinel-2影像的基础上，利用gee中的filterdate()筛选器获得大豆关键生育期内预处理后的sentinel-2影像。
[0090]
此处，大豆关键生育期内预处理后的sentinel-2影像的5个波段，包含红边1、红边2、红边3、短波红外1、短波红外2。在大豆关键生育期内预处理后的sentinel-2影像的基础
上，利用gee中的select()筛选出这5个波段，获得大豆关键生育期内预处理后的sentinel-2影像的5个波段。
[0091]
此处，百分位数的计算方法可表述如下：
[0092][0093][0094]
其中，k表示15/50/90，n为某个像素大豆关键生育期内预处理后的sentinel-2影像清晰观察的总数，r是第k百分位的排名，pk是第k百分位的百分位合成值。具体来说，如果公式(2)获得的r是整数，则第k个百分位是第r个和第(r 1)个值的平均像素值；如果公式(2)获得的r不是整数，则将其四舍五入到最接近的整数。
[0095]
此处，最大值、最小值、平均值、标准差的计算方法分别可表述如下：
[0096][0097][0098][0099][0100]
其中，n为某个像素大豆关键生育期内预处理后的sentinel-2影像清晰观察的总数，xj表示n个观测值中第j个观测值，max表示取n个观测值中的最大值函数，min表示取n个观测值中的最小值函数，mean表示取n个观测值中的平均值函数,stddev表示取n个观测值中的标准差函数。
[0101]
3)利用航天飞机雷达地形测绘任务(srtm)数字高程数据提取海拔及坡度特征(简称地形特征)，再利用mask函数对所提特征进行耕地掩膜。
[0102]
此处，在gee中直接调用航天飞机雷达地形测绘任务(srtm)数字高程数据，其中已包含海拔波段，利用ee.algorithms.terrain()算法计算坡度数据，获得海拔及坡度特征。
[0103]
在步骤s3，分类模型的构建与训练：构建分类模型并通过训练设置参数。
[0104]
具体地，依据随机森林(rf)/支持向量机(svm)/反向传播神经网络(bpnn)三种分类精度，最终选择rf，即随机森林，随机森林分类模型采用随机重复自抽样方法有放回的在原始数据中抽取k组数据集，每组数据集所包含的数据量为原始数据总量的约三分之二。通过k组数据集选定合适的特征节点数p构建k棵决策树，集合k棵决策树对结果进行简单的投票即可获得所需分类器。其中特征数量和决策树棵树是决定决策树生长的主要参数，该参数的最优解是通过袋外误差获取。袋外误差指未参与模型构建的三分之一数据集对常规误差进行无偏估计所获得的结果。具体计算方法是将未参与模型构建的三分之一数据集，应用生成的分类器，对其进行分类获得分类结果，由于该部分的类别为已知，故将分类器生成
的分类结果与已知类别进行比对，计算分类器所得的每一类别的错误分类结果的占比即为该类别的分类误差，将所有类别的误差通过均值计算可得该分类器的平均袋外误差。袋外误差具有高效性，且与交叉验证的结果相近，故而在随机森林分类中无需进行交叉验证或采用独立数据建立误差无偏估计。
[0105]
其具体方法包括：
[0106]
1)基于原始训练集，采用随机且有放回地抽样bagging方法产生每棵决策树的训练子集；
[0107]
2)构建随机森林分类模型，设定随机森林由多棵cart决策树组成，cart决策树采用基尼gini系数选择最优特征；其具体步骤如下：
[0108]
①
计算每个特征的基尼系数，选择基尼系数最小的特征x进行节点划分，基尼系数的计算公式如下所示：
[0109][0110]
其中，m1和m2是根据特征x的某个属性值，将m分成的两部分数据集；
[0111]
②
对于特征x，根据某一属性值将其划分为两个子集m1和m2，计算该属性值进行节点划分得到的基尼系数，计算公式如下所示：
[0112][0113]
其中，m表示给定样本集中样本数量，k表示类别个数，mi表示第i个类别的样本数量；
[0114]
③
遍历特征x的所有属性值，选择基尼系数最小的属性值作为特征x的最优划分节点值；
[0115]
④
不断遍历这棵树的特征子集，重复步骤
①
和步骤
②
直至所有的特征都被选择完毕或子数据集都属于同一类；
[0116]
3)利用训练好的随机森林模型对分类样本进行预测，其中每棵树都得到一个独立的预测结果，对每棵树的预测结果进行汇总统计，按照投票的方式，将获得票数最多的类作为最终结果。
[0117]
随机森林模型中的需要设置的参数主要有决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数。
[0118]
①
决策树的数量：决策树的数量是影响分类精度的重要参数之一，若该数值设置过小，会使模型对数据的拟合能力不足，而若决策树数目设置过大，分类精度会提高，但模型的复杂程度和运行时间会大幅增加。
②
最大特征数：随机森林在构建每棵树时，不是所有的特征都参与节点分裂过程，而是随机地选择某些特征，最大特征数是指从特征空间可随机抽取的特征的最大数量。若该数值设置过小，则决策树的分类能力较弱，反之，最大特征数设置过大，单棵决策树的分类能力可能会上升，但会导致决策树之间的相关性增加，从而使得随机森林的性能下降。常用的最大特征数取值有：所有特征数、所有特征数量的开平方、所有特征数量的对数。
③
决策树深度：即随机森林模型中决策树生长的最大深度，若深度值设置过小，单棵决策树的分类能力较弱，降低其分类精度，设置过大则容易引起模型过度拟合，同时会增加模型的复杂程度及运行时间。
④
叶子节点最少样本数：该值决定决策子
树中叶子节点是否剪枝，当子树中叶子节点地样本数目小于该值，则舍弃该节点。
⑤
节点划分的最少样本数：当模型中决策子树节点的样本数小于该值时，不再选取最优特征对其进行划分。
[0119]
在gee中调用随机森林分类模型，主要设置以下参数：森林中决策树的数量。在gee中将棵数分别设置为50至400每次增加50。依据棵树大于100且分类精度第一次达到局部最大值确定树的数量。其中，由于随机森林每次取样的随机性，为了避免每次重复实验结果的细微差别，设置随机种子(seed)为999。其他参数保持默认值。
[0120]
在步骤s4，以所述时序谱特征为主，将所述时序谱特征、光谱特征和地形特征进行组合得到待选特征组合，得到待选特征组合；将所述待选特征组合输入分类模型，计算待选特征组合的分类精度，根据分类精度得到最优特征组合。
[0121]
在一些实施例中，在所述步骤s4中，所述待选特征组合的具体组合方式包括：
[0122]
时序谱特征、光谱特征、地形特征、时序谱特征及光谱特征、时序谱特征及地形特征、时序谱特征和光谱特征以及地形特征，共6组待选特征组合。
[0123]
计算待选特征组合的分类精度的方法包括：
[0124]
采用基于混淆矩阵的随机森林分类精度评价方法计算待选特征组合的分类精度。
[0125]
具体地，1)以时序谱特征为主将不同类别特征进行组合得到待选特征组合时序谱特征、光谱特征、地形特征、时序谱特征及光谱特征、时序谱特征及地形特征、时序谱特征和光谱特征以及地形特征，共6组待选特征组合。
[0126]
2)分别计算待选特征组合，利用随机森林分类精度；
[0127]
其中随机森林分类精度评价是基于混淆矩阵的统计方法，混淆矩阵定义如下所示：
[0128][0129]
式中，n表示类别数，m
ij
表示实际属于i类的像元在分类结果图中被分到j类的像元数量，对角线元素的值即是各类别被正确分类的像元数，因此混淆矩阵中对角线元素的值越大，表示被正确分类的像元数越多，分类结果越可靠。利用混淆矩阵总体精度(oa)、kappa系数、生产者精度(pa)、用户精度(ua)和f1-score。其计算公式分别为：
[0130][0131][0132][0133]
[0134][0135]
式中，m
i
为混淆矩阵中的行总和，m
i
为混淆矩阵中的列总和；
[0136]
3)以时序谱特征为主形成不同的待选特征组合，对比不同待选特征组合随机森林分类精度，选择分类精度均最高的特征组合作为最优特征组合。
[0137]
在步骤s5，大豆种植区提取及精度评价：将最优特征组合输入随机森林分类模型提取大豆种植区分布图并进行精度评价。
[0138]
具体地，将最优特征组合输入随机森林分类模型提取大豆种植区分布图，并计算随机森林分类精度及大豆种植面积提取精度，即大豆提取面积与官方农业统计数据对比的精度，其中大豆种植面积提取精度评价指标包括：相对误差、均方根误差，计算公式如下所示，
[0139][0140]
其中，re表示相对误差，s表示提取的大豆面积，s’表示大豆面积农业统计数据；
[0141][0142]
其中，rmse表示均方根误差，fi表示第i个区域大豆提取面积，yi表示第i个区域大豆农业统计面积，n表示区域个数。
[0143]
具体计算结果：
[0144]
借助google earth engine(gee)云平台，利用多种遥感数据提取不同类别特征，然后结合随机森林分类模型探究不同类别特征组合对大豆识别精度的影响，最终选择最优特征组合提取大豆种植区。结果如下：如表2所示时序谱特征、光谱特征及地形特征组合为最优组合，由此提取巴西2020/2021生长季大豆空间分布图，如图2所示。其中，oa为0.93，kappa为0.86，大豆的pa、ua和f1分别为0.94、0.93和0.94；大豆面积为36612千公顷与农业统计数据的相对误差为-5.94％；大豆面积省级数据与农业统计数据对比均方根误差为590千公顷。
[0145]
表2
[0146][0147]
综上，本发明提出的方案提供了大区域尺度大豆种植区快速准确的识别方法，为大豆面积、产量、病虫害监测等提供科学依据。
[0148]
本发明第二方面公开了一种基于植被指数时序谱特征的大豆种植区提取系统。图3为根据本发明实施例的一种基于植被指数时序谱特征的大豆种植区提取系统的结构图；如图3所示，所述系统100包括：
[0149]
第一处理模块101，被配置为，遥感数据预处理：对sentinel-2卫星影像进行预处理，并进行耕地掩膜；
[0150]
第二处理模块102，被配置为，特征提取：提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征，简称时序谱特征；提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征，简称光谱特征；提取海拔及坡度特征，简称地形特征，并利用耕地数据对所述时序谱特征、光谱特征和地形特征分别进行掩膜；
[0151]
第三处理模块103，被配置为，分类模型的构建与训练：构建分类模型并通过训练设置参数；
[0152]
第四处理模块104，被配置为，以所述时序谱特征为主，将所述时序谱特征、光谱特征和地形特征进行组合得到待选特征组合，得到待选特征组合；将所述待选特征组合输入分类模型，计算待选特征组合的分类精度，根据分类精度得到最优特征组合；
[0153]
第五处理模块105，被配置为，大豆种植区提取及精度评价：将最优特征组合输入随机森林分类模型提取大豆种植区分布图并进行精度评价。
[0154]
根据本发明第二方面的系统，第二处理模块102，被配置为，所述提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征包括：
[0155]
利用线性谐波模型提取大豆整个生育期内掩膜后的sentinel-2影像增强型植被指数时序谱特征；
[0156]
所述线性谐波模型的公式为：
[0157][0158]
其中，f(t)是第t时刻的拟合的增强型植被指数的值，ak是余弦系数，bk是正弦系数，c是截距项，n是谐波级数的阶数，w是频率等于1.5；自变量t是一年中的某一天；n分别设置为1至5，每次增加1，依据大豆的时间序列中原始值与拟合值之间均方误差最小确定n的值，此处n取2；取相位和振幅作为时序特征输入，振幅定义为二维矢量[ak,bk]的长度，相位定义为二维矢量[ak,bk]形成的角度。
[0159]
根据本发明第二方面的系统，第二处理模块102，被配置为，所述提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征包括：
[0160]
利用百分位数、最大值、最小值、平均值和标准差方法提取大豆关键生育期内掩膜后的sentinel-2影像波段反射率特征。
[0161]
根据本发明第一方面的方法在所述步骤s2中，大豆关键生育期内掩膜后的sentinel-2影像波段包含，红边1、红边2、红边3、短波红外1和短波红外2；
[0162]
所述关键生育期为1-2月，或者为8-9月，具体包含大豆开花期、结荚期和灌浆期三个生长阶段。
[0163]
根据本发明第二方面的系统，第二处理模块102，被配置为，所述提取海拔及坡度特征包括：
[0164]
利用航天飞机雷达地形测绘任务数字高程数据提取海拔及坡度特征。
[0165]
根据本发明第二方面的系统，第四处理模块104，被配置为，所述待选特征组合的
具体组合方式包括：
[0166]
时序谱特征、光谱特征、地形特征、时序谱特征及光谱特征、时序谱特征及地形特征、时序谱特征和光谱特征以及地形特征，共6组待选特征组合。
[0167]
根据本发明第二方面的系统，第四处理模块104，被配置为，计算待选特征组合的分类精度包括：
[0168]
采用基于混淆矩阵的随机森林分类精度评价方法计算待选特征组合的分类精度。
[0169]
本发明第三方面公开了一种电子设备。电子设备包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时，实现本发明公开第一方面中任一项的一种基于植被指数时序谱特征的大豆种植区提取方法中的步骤。
[0170]
图4为根据本发明实施例的一种电子设备的结构图，如图4所示，电子设备包括通过系统总线连接的处理器、存储器、通信接口、显示屏和输入装置。其中，该电子设备的处理器用于提供计算和控制能力。该电子设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该电子设备的通信接口用于与外部的终端进行有线或无线方式的通信，无线方式可通过wifi、运营商网络、近场通信(nfc)或其他技术实现。该电子设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该电子设备的输入装置可以是显示屏上覆盖的触摸层，也可以是电子设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
[0171]
本领域技术人员可以理解，图4中示出的结构，仅仅是与本公开的技术方案相关的部分的结构图，并不构成对本技术方案所应用于其上的电子设备的限定，具体的电子设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
[0172]
本发明第四方面公开了一种计算机可读存储介质。计算机可读存储介质上存储有计算机程序，计算机程序被处理器执行时，实现本发明公开第一方面中任一项的一种基于植被指数时序谱特征的大豆种植区提取方法中的步骤中的步骤。
[0173]
请注意，以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。以上实施例仅表达了本技术的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。