一种大豆种植区识别和面积测算方法

本技术涉及一种大豆种植区识别和面积测算方法，主要适用于基于sentinel1/2遥感时间序列影像的大豆种植区识别和面积测算。

背景技术：

1、遥感技术是观测地球表面和获取信息的先进技术，遥感卫星可以宏观、动态、持续地获取地表影像，进而实现地物的识别提取和时空分析。其中哨兵系列卫星具有较高的空间分辨率和时间分辨率，可实现大范围、长时序和多模态的地表信息获取，满足大豆种植区精确识别的要求。

2、传统的大豆种植区制图方法主要是通过人工田间调查和上报，这种方法耗费人力、物力和财力，且效率低下，难以应用在大范围的制图任务中。遥感技术凭借精细、宏观、高效的优点，能够以低成本、高效率的方式完成大范围的大豆种植区分布制图，并保证制图的客观性和准确性。

3、鉴于大豆在生长发育过程中光谱反射率会出现明显的季相变化特征，利用遥感影像时间序列识别大豆可以较为完整地刻画出大豆的生长周期变化，从而准确提取出其关键生长期特征和物候特征，为大豆种植区的精确识别提供重要信息。除光谱反射率外，遥感影像中还可以提取出丰富的植被指数，植被指数可以有效指示作物的生长状况和叶面积指数、水分含量等特征，已广泛应用在农作物识别中。

4、虽然利用光学遥感影像识别大豆的技术已相对成熟，但光学影像易受云雨影响，在阴雨天气往往无法获得有效的地面观测数据。而大豆精确识别依赖于密集、连续的影像序列，若其关键生长期的影像受到严重影响，其识别效果则会大大降低，因此需要降低云雨等不利天气因素所造成的影响。

技术实现思路

1、本技术解决的技术问题是克服现有技术中存在的上述不足，而提供一种使用方便，天气因素影响小，识别精度高的大豆种植区识别和面积测算方法。

2、本技术解决上述技术问题所采用的技术方案是：一种大豆种植区识别和面积测算方法，其特征是包括以下步骤：

3、步骤s1、通过地面调查获取目标区域内包括大豆在内的生长期重叠且相近作物的样点数据集，分析目标区域内所述作物的种植日历，确定其主要生长期的时间跨度(确定目标区域范围、作物种类、样点数据集及作物生长期)；

4、步骤s2、获取目标区域的sentinel2影像并进行预处理，选择七个原始波段并计算五个植被指数，得到sentinel2影像时间序列；获取目标区域的sentinel1影像并进行预处理，选择两个极化通道后向散射系数并计算三个极化通道组合，得到sentinel1影像时间序列。

5、步骤s3、提取sentinel2影像的七个原始波段和五个植被指数的统计特征、生长高峰期特征、谐波拟合特征以及物候特征作为大豆在sentinel2影像中的时序特征；提取sentinel1影像的两个极化通道后向散射系数和三个极化通道组合的统计特征和主成分特征作为大豆在sentinel1影像中的时序特征；在gee中获取esa 10m地表覆盖数据集，去除步骤s2、步骤s3中所有不在地表覆盖数据集范围内的数据(取地表覆盖数据与sentinel1、sentinel2数据集的交集)，得到掩膜处理后的目标区域数据；

6、步骤s4、将样点数据集上传到gee平台中，并将其划分为训练集和测试集；将步骤s3中特征参数提取到训练集样点中，构建随机森林分类模型，并将带有步骤s3所述特征参数的训练集输入到随机森林分类模型中完成模型训练；

7、步骤s5、将带有步骤3所述特征参数的目标区域数据输入到训练好的随机森林分类模型中，得到作物分类结果并从中提取出大豆种植区，计算大豆种植面积。

8、步骤s2中，获取样点数据集的sentinel2影像并进行预处理包含以下步骤：

9、在gee平台中获取大豆生长期内的所有sentinel2卫星地表反射率数据影像，将所有波段数值除以10000得到真实反射率数值；通过波段筛选选择sentinel2影像的七个原始波段，分别是蓝、红、红边1、红边2、红边3、近红外和短波红外1，在sentinel2影像的名称中分别为b2、b4、b5、b6、b7、b8和b11；计算sentinel2影像的五个植被指数，分别是增强型植被指数evi、地表水分指数lswi、归一化植被衰老指数ndsvi、红边归一化植被指数rendvi和红边位置rep；利用云得分算法识别每幅影像中的含云像素并将之去除；然后采用10天中值合成的方式等间隔生成目标区域sentinel2影像时间序列，使用线性插值的方式对10天中值数据集中的空缺像素进行填补，应用sg滤波处理sentinel2影像时间序列。

10、其中sentinel2影像的植被指数计算公式如下：

11、

12、其中，red、blue、nir分别代表sentinel2影像的红光波段、蓝光波段和近红外波段的反射率。

13、

14、其中nir和swir1分别代表sentinel2影像的近红外波段和短波红外1波段的反射率。

15、

16、其中nir和re2分别代表sentinel2影像的近红外波段和红边2波段的反射率。

17、

18、其中red和swir1分别代表sentinel2影像的红光波段和短波红外1波段的反射率。

19、

20、其中re1、re2、re3和red分别代表sentinel2影像的红边1波段、红边2波段、红边3波段和红光波段的反射率。

21、步骤s2中，获取样点数据集的sentinel1影像并进行预处理包含以下步骤：

22、在gee平台中获取大豆生长期内的所有sentinel1的grd地距影像，成像方式为iw干涉宽带扫描；通过波段筛选选择sentinel1影像的两个极化通道后向散射系数，分别是垂直同向极化和垂直/水平交叉极化，即vv和vh；计算sentinel1影像的三个极化通道组合，分别是交叉共极化之和、交叉共极化比率和雷达植被指数rvi。利用局部入射角对原始后向散射系数值进行入射角归一化，减弱入射角的影响；应用7×7精炼lee滤波对sar影像中的噪点进行滤除；然后采用10天中值合成的方式等间隔生成目标区域sentinel1影像时间序列，使用线性插值的方式对10天中值数据集中的空缺像素进行填补，最后应用sg滤波处理sentinel1影像时间序列。

23、其中sentinel1影像的极化通道组合方式如下：

24、交叉共极化之和σvh0+σvv0

25、交叉共极化比率σvh0/σvv0

26、雷达植被指数

27、其中σvh0和σvv0分别为sentinel1影像的vh和vv通道的后向散射系数值。

28、步骤s3中，所述sentinel2的统计特征包括特征参数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数；所述sentinel2生长高峰期特征为大豆生长高峰期时的特征参数值，将evi或lswi达到峰值的时间作为大豆生长高峰期；所述sentinel2谐波拟合特征是通过对evi时序曲线做二阶谐波拟合得到的拟合系数；所述sentinel2物候特征采用阈值法从evi时序曲线中估算得到的生长季起始日sos、生长季结束日eos和生长季长度los。

29、所述二阶谐波拟合模型计算公式如下：

30、

31、其中f为拟合值，a0、ak和bk分别为谐波余项、余弦系数和正弦系数，α0实际上为时序的平均值。n为谐波拟合的阶次，w为频率取1.5，t为当前日期在一年中的位置，在这里的数值为0到1之间的小数。

32、所述阈值法为一种作物物候期估计方法，为现有技术，通过对evi时序曲线使用阈值法得到：将evi时间序列的值按大小排序，计算出50％和95％分位数作为阈值，然后找到第一个超过50％分位数的evi值所对应的日期作为sos时间点，找到最后一个超过evi95％分位数的日期作为eos时间点。该方法对于大豆、玉米、水稻这种一年一熟制的作物会有较好的估计效果，因为其evi时间序列存在明显的增长期和衰退期。

33、步骤s3中，所述sentinel1统计特征包括特征参数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数；所述sentinel1主成分特征为每个特征参数在大豆生长期内进行主成分分析后的前三个主成分分量。

34、所述主成分分析步骤如下：

35、①构建时序影像数据矩阵x

36、

37、其中m为影像时间序列中的影像数量，n为每幅影像中的像元数。矩阵中第i行向量对应影像序列的第i幅影像，由影像的上下相邻两行数据首尾依次相接得到；

38、②根据原始图像数据矩阵x，计算协方差矩阵s

39、

40、其中为第i幅影像的均值。

41、③计算s矩阵的特征值和特征向量(均为m×1的列向量)；

42、④构建变换矩阵t

43、将特征值按其大小排列，排序后特征值依次为λ1≥λ2≥…λm，求得各特征值对应的单位特征向量：uj＝[u1j，u2j，…，umj]t。然后以各特征向量为列构成矩阵u，u即为[u1，u2，...，um]＝[uij]m×m，u的转置矩阵即为所求的pca变换的变换矩阵t。

44、⑤计算主分量

45、利用主成分变换矩阵计算y

46、

47、其中矩阵y的每一行即为一个主成分分量，yj＝[yj1，yj2，...，yjn]为第j主成分分量。

48、步骤s4中，所述构建随机森林分类模型通过调用gee云平台的随机森林函数完成，设置的参数包括：决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数，其中将决策树的数量设置为100，叶子节点最少样本数设为10，其他参数保持默认值。

49、所述作物为大豆、玉米和水稻，所述步骤s5包括：将带有步骤s3中的所有特征参数的目标区域数据输入到训练好的随机森林分类模型中，输出得到大豆、玉米和水稻的分类图，大豆、玉米、水稻的像素值分别是2、1、0，利用eq()函数提取出像素值为2的目标区域，得到大豆种植区，通过对所有大豆种植区面积的累加计算出大豆种植区总面积。

50、本技术与现有技术相比，具有以下优点和效果：天气因素影响小，识别精度高，使用方便。