CN116645603A - 一种大豆种植区识别和面积测算方法 - Google Patents
一种大豆种植区识别和面积测算方法 Download PDFInfo
- Publication number
- CN116645603A CN116645603A CN202310522156.5A CN202310522156A CN116645603A CN 116645603 A CN116645603 A CN 116645603A CN 202310522156 A CN202310522156 A CN 202310522156A CN 116645603 A CN116645603 A CN 116645603A
- Authority
- CN
- China
- Prior art keywords
- image
- sentinel2
- soybean
- sentinel1
- red
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 244000068988 Glycine max Species 0.000 title claims abstract description 90
- 235000010469 Glycine max Nutrition 0.000 title claims abstract description 90
- 238000000691 measurement method Methods 0.000 title claims abstract description 15
- 238000000034 method Methods 0.000 claims abstract description 23
- 238000007637 random forest analysis Methods 0.000 claims abstract description 18
- 238000012549 training Methods 0.000 claims abstract description 17
- 238000007781 pre-processing Methods 0.000 claims abstract description 11
- 238000012545 processing Methods 0.000 claims abstract description 8
- 238000012360 testing method Methods 0.000 claims abstract description 3
- 239000011159 matrix material Substances 0.000 claims description 35
- 238000002310 reflectometry Methods 0.000 claims description 22
- 230000010287 polarization Effects 0.000 claims description 15
- 230000009466 transformation Effects 0.000 claims description 12
- 238000001914 filtration Methods 0.000 claims description 10
- 239000013598 vector Substances 0.000 claims description 10
- 235000007164 Oryza sativa Nutrition 0.000 claims description 8
- 240000008042 Zea mays Species 0.000 claims description 8
- 235000005824 Zea mays ssp. parviglumis Nutrition 0.000 claims description 8
- 235000002017 Zea mays subsp mays Nutrition 0.000 claims description 8
- 235000005822 corn Nutrition 0.000 claims description 8
- 238000000513 principal component analysis Methods 0.000 claims description 8
- 235000009566 rice Nutrition 0.000 claims description 8
- 230000015572 biosynthetic process Effects 0.000 claims description 6
- 238000003066 decision tree Methods 0.000 claims description 6
- 238000003786 synthesis reaction Methods 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000012216 screening Methods 0.000 claims description 5
- WPPDXAHGCGPUPK-UHFFFAOYSA-N red 2 Chemical compound C1=CC=CC=C1C(C1=CC=CC=C11)=C(C=2C=3C4=CC=C5C6=CC=C7C8=C(C=9C=CC=CC=9)C9=CC=CC=C9C(C=9C=CC=CC=9)=C8C8=CC=C(C6=C87)C(C=35)=CC=2)C4=C1C1=CC=CC=C1 WPPDXAHGCGPUPK-UHFFFAOYSA-N 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 3
- 239000002352 surface water Substances 0.000 claims description 3
- 230000032683 aging Effects 0.000 claims description 2
- 238000004422 calculation algorithm Methods 0.000 claims description 2
- 238000010276 construction Methods 0.000 claims description 2
- 238000005388 cross polarization Methods 0.000 claims description 2
- 238000003384 imaging method Methods 0.000 claims description 2
- 238000007670 refining Methods 0.000 claims description 2
- 240000007594 Oryza sativa Species 0.000 claims 3
- 230000000694 effects Effects 0.000 description 7
- 241000209094 Oryza Species 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 238000011835 investigation Methods 0.000 description 4
- 238000005070 sampling Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000000605 extraction Methods 0.000 description 3
- 230000003287 optical effect Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000002411 adverse Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000003595 spectral effect Effects 0.000 description 2
- 238000010200 validation analysis Methods 0.000 description 2
- 230000003313 weakening effect Effects 0.000 description 2
- HQHACWYQWDKBOZ-UHFFFAOYSA-N 1-nitroso-azacyclotridecane Chemical compound O=NN1CCCCCCCCCCCC1 HQHACWYQWDKBOZ-UHFFFAOYSA-N 0.000 description 1
- 208000003643 Callosities Diseases 0.000 description 1
- 206010020649 Hyperkeratosis Diseases 0.000 description 1
- 239000000443 aerosol Substances 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000009758 senescence Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/188—Vegetation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/0464—Convolutional networks [CNN, ConvNet]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/764—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/82—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using neural networks
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A40/00—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
- Y02A40/10—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production in agriculture
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Software Systems (AREA)
- Computing Systems (AREA)
- Multimedia (AREA)
- Databases & Information Systems (AREA)
- Medical Informatics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Life Sciences & Earth Sciences (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Molecular Biology (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Image Processing (AREA)
Abstract
本申请涉及一种大豆种植区识别和面积测算方法，包括以下步骤：步骤S1、确定目标区域范围、作物种类、样点数据集及作物生长期：步骤S2、获取Sentinel2影像并进行预处理；获取Sentinel1影像并预处理；步骤S3、提取Sentinel2影像时序特征；提取Sentinel1影像时序特征；掩膜处理得到目标区域数据；步骤S4、将样点数据集上传到GEE平台中并划分为训练集和测试集；将步骤S3中特征参数提取到训练集中，构建随机森林分类模型并用所述训练集完成模型训练；步骤S5、将目标区域数据输入到随机森林分类模型中得到作物分类结果,提取大豆种植区，计算大豆种植面积。本申请天气因素影响小，识别精度高。
Description
技术领域
本申请涉及一种大豆种植区识别和面积测算方法，主要适用于基于Sentinel1/2遥感时间序列影像的大豆种植区识别和面积测算。
背景技术
遥感技术是观测地球表面和获取信息的先进技术，遥感卫星可以宏观、动态、持续地获取地表影像，进而实现地物的识别提取和时空分析。其中哨兵系列卫星具有较高的空间分辨率和时间分辨率，可实现大范围、长时序和多模态的地表信息获取，满足大豆种植区精确识别的要求。
传统的大豆种植区制图方法主要是通过人工田间调查和上报，这种方法耗费人力、物力和财力，且效率低下，难以应用在大范围的制图任务中。遥感技术凭借精细、宏观、高效的优点，能够以低成本、高效率的方式完成大范围的大豆种植区分布制图，并保证制图的客观性和准确性。
鉴于大豆在生长发育过程中光谱反射率会出现明显的季相变化特征，利用遥感影像时间序列识别大豆可以较为完整地刻画出大豆的生长周期变化，从而准确提取出其关键生长期特征和物候特征，为大豆种植区的精确识别提供重要信息。除光谱反射率外，遥感影像中还可以提取出丰富的植被指数，植被指数可以有效指示作物的生长状况和叶面积指数、水分含量等特征，已广泛应用在农作物识别中。
虽然利用光学遥感影像识别大豆的技术已相对成熟，但光学影像易受云雨影响，在阴雨天气往往无法获得有效的地面观测数据。而大豆精确识别依赖于密集、连续的影像序列，若其关键生长期的影像受到严重影响，其识别效果则会大大降低，因此需要降低云雨等不利天气因素所造成的影响。
发明内容
本申请解决的技术问题是克服现有技术中存在的上述不足，而提供一种使用方便，天气因素影响小，识别精度高的大豆种植区识别和面积测算方法。
本申请解决上述技术问题所采用的技术方案是：一种大豆种植区识别和面积测算方法，其特征是包括以下步骤：
步骤S1、通过地面调查获取目标区域内包括大豆在内的生长期重叠且相近作物的样点数据集，分析目标区域内所述作物的种植日历，确定其主要生长期的时间跨度(确定目标区域范围、作物种类、样点数据集及作物生长期)；
步骤S2、获取目标区域的Sentinel2影像并进行预处理，选择七个原始波段并计算五个植被指数，得到Sentinel2影像时间序列；获取目标区域的Sentinel1影像并进行预处理，选择两个极化通道后向散射系数并计算三个极化通道组合，得到Sentinel1影像时间序列。
步骤S3、提取Sentinel2影像的七个原始波段和五个植被指数的统计特征、生长高峰期特征、谐波拟合特征以及物候特征作为大豆在Sentinel2影像中的时序特征；提取Sentinel1影像的两个极化通道后向散射系数和三个极化通道组合的统计特征和主成分特征作为大豆在Sentinel1影像中的时序特征；在GEE中获取ESA 10m地表覆盖数据集，去除步骤S2、步骤S3中所有不在地表覆盖数据集范围内的数据(取地表覆盖数据与Sentinel1、Sentinel2数据集的交集)，得到掩膜处理后的目标区域数据；
步骤S4、将样点数据集上传到GEE平台中，并将其划分为训练集和测试集；将步骤S3中特征参数提取到训练集样点中，构建随机森林分类模型，并将带有步骤S3所述特征参数的训练集输入到随机森林分类模型中完成模型训练；
步骤S5、将带有步骤3所述特征参数的目标区域数据输入到训练好的随机森林分类模型中，得到作物分类结果并从中提取出大豆种植区，计算大豆种植面积。
步骤S2中，获取样点数据集的Sentinel2影像并进行预处理包含以下步骤：
在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影像时间序列。
其中Sentinel2影像的植被指数计算公式如下：
其中，Red、Blue、NIR分别代表Sentinel2影像的红光波段、蓝光波段和近红外波段的反射率。
其中NIR和SWIR1分别代表Sentinel2影像的近红外波段和短波红外1波段的反射率。
其中NIR和RE2分别代表Sentinel2影像的近红外波段和红边2波段的反射率。
其中Red和SWIR1分别代表Sentinel2影像的红光波段和短波红外1波段的反射率。
其中RE1、RE2、RE3和Red分别代表Sentinel2影像的红边1波段、红边2波段、红边3波段和红光波段的反射率。
步骤S2中，获取样点数据集的Sentinel1影像并进行预处理包含以下步骤：
在GEE平台中获取大豆生长期内的所有Sentinel1的GRD地距影像，成像方式为IW干涉宽带扫描；通过波段筛选选择Sentinel1影像的两个极化通道后向散射系数，分别是垂直同向极化和垂直/水平交叉极化，即VV和VH；计算Sentinel1影像的三个极化通道组合，分别是交叉共极化之和、交叉共极化比率和雷达植被指数RVI。利用局部入射角对原始后向散射系数值进行入射角归一化，减弱入射角的影响；应用7×7精炼Lee滤波对SAR影像中的噪点进行滤除；然后采用10天中值合成的方式等间隔生成目标区域Sentinel1影像时间序列，使用线性插值的方式对10天中值数据集中的空缺像素进行填补，最后应用SG滤波处理Sentinel1影像时间序列。
其中Sentinel1影像的极化通道组合方式如下：
交叉共极化之和σVH 0+σVV 0
交叉共极化比率σVH 0/σVV 0
雷达植被指数
其中σVH 0和σVV 0分别为Sentinel1影像的VH和VV通道的后向散射系数值。
步骤S3中，所述Sentinel2的统计特征包括特征参数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数；所述Sentinel2生长高峰期特征为大豆生长高峰期时的特征参数值，将EVI或LSWI达到峰值的时间作为大豆生长高峰期；所述Sentinel2谐波拟合特征是通过对EVI时序曲线做二阶谐波拟合得到的拟合系数；所述Sentinel2物候特征采用阈值法从EVI时序曲线中估算得到的生长季起始日SOS、生长季结束日EOS和生长季长度LOS。
所述二阶谐波拟合模型计算公式如下：
其中f为拟合值，a0、ak和bk分别为谐波余项、余弦系数和正弦系数，α0实际上为时序的平均值。n为谐波拟合的阶次，w为频率取1.5，t为当前日期在一年中的位置，在这里的数值为0到1之间的小数。
所述阈值法为一种作物物候期估计方法，为现有技术，通过对EVI时序曲线使用阈值法得到：将EVI时间序列的值按大小排序，计算出50％和95％分位数作为阈值，然后找到第一个超过50％分位数的EVI值所对应的日期作为SOS时间点，找到最后一个超过EVI95％分位数的日期作为EOS时间点。该方法对于大豆、玉米、水稻这种一年一熟制的作物会有较好的估计效果，因为其EVI时间序列存在明显的增长期和衰退期。
步骤S3中，所述Sentinel1统计特征包括特征参数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数；所述Sentinel1主成分特征为每个特征参数在大豆生长期内进行主成分分析后的前三个主成分分量。
所述主成分分析步骤如下：
①构建时序影像数据矩阵X
其中m为影像时间序列中的影像数量，n为每幅影像中的像元数。矩阵中第i行向量对应影像序列的第i幅影像，由影像的上下相邻两行数据首尾依次相接得到；
②根据原始图像数据矩阵X，计算协方差矩阵S
其中 为第i幅影像的均值。
③计算S矩阵的特征值和特征向量(均为m×1的列向量)；
④构建变换矩阵T
将特征值按其大小排列，排序后特征值依次为λ1≥λ2≥…λm，求得各特征值对应的单位特征向量：Uj＝[u1j，u2j，…，umj]T。然后以各特征向量为列构成矩阵U，U即为[U1，U2，...，Um]＝[uij]m×m，U的转置矩阵即为所求的PCA变换的变换矩阵T。
⑤计算主分量
利用主成分变换矩阵计算Y
其中矩阵Y的每一行即为一个主成分分量，Yj＝[yj1，yj2，...，yjn]为第j主成分分量。
步骤S4中，所述构建随机森林分类模型通过调用GEE云平台的随机森林函数完成，设置的参数包括：决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数，其中将决策树的数量设置为100，叶子节点最少样本数设为10，其他参数保持默认值。
所述作物为大豆、玉米和水稻，所述步骤S5包括：将带有步骤S3中的所有特征参数的目标区域数据输入到训练好的随机森林分类模型中，输出得到大豆、玉米和水稻的分类图，大豆、玉米、水稻的像素值分别是2、1、0，利用eq()函数提取出像素值为2的目标区域，得到大豆种植区，通过对所有大豆种植区面积的累加计算出大豆种植区总面积。
本申请与现有技术相比，具有以下优点和效果：天气因素影响小，识别精度高，使用方便。
附图说明
图1是本申请实施例的主要流程示意图。
图2是本申请实施例研究区主要作物样本点分布及耕地掩膜范围示意图。
图3是本申请实施例研究区2019年大豆种植区空间分布示意图。
具体实施方式
下面结合附图并通过实施例对本申请作进一步的详细说明，以下实施例是对本申请的解释而本申请并不局限于以下实施例。
参见图1～图3，本申请公开了一种基于Sentinel1/2遥感时序影像的大豆种植区识别和面积测算方法。本实施例的局部地理区域为中国大豆主产区黑龙江省九三垦区，研究作物为大豆作物，如图2所示表示本实施例所选择九三垦区内的三个市，分别是黑河市、绥化市和齐齐哈尔市。图1为本申请实施例基于Sentinel1/2遥感时序影像的大豆种植区识别和面积测算方法的流程图，所述方法包括：
步骤S1、准备包含大豆在内的主要作物地面调查样点，获取作物的生长期时间跨度；
所述作物地面调查样点包括大豆、水稻和玉米，通过人工野外实地采样得到，采样人员通过在农田中央手持GPS仪记录样点位置和作物类型。各类型作物的样点数量应当均衡，均匀地分布在感兴趣区域；主要作物生长期时间跨度同样通过实地调查获取，用于确定影像时间序列的时间范围，确保影像时间序列能够覆盖大豆的完整生长期，同时避免影像时间序列过长包含了无用信息。
步骤S2、Sentinel1/2时序影像获取及预处理：
在GEE中获取Sentinel2 Level1C级的TOA影像数据，时间范围为2019年5月1日-10月10日，空间范围为九三垦区的矢量边界。
将所有波段数值除以10000得到真实反射率数值，利用气溶胶波段、蓝光波段、绿光波段、红光波段和卷云波段，以及归一化湿度指数NDMI计算所有影像中像素的云得分，识别出含云像素并将之去除。
通过波段筛选选择Sentinel2影像的七个原始波段，分别是蓝、红、红边1、红边2、红边3、近红外和短波红外1，在Sentinel2影像的名称中分别为B2、B4、B5、B6、B7、B8和B11；计算Sentinel2影像的五个植被指数，分别是增强型植被指数EVI、地表水分指数LSWI、归一化植被衰老指数NDSVI、红边归一化植被指数RENDVI和红边位置REP。
其中，Red、Blue、NIR分别代表Sentinel2遥感影像的红光波段、蓝光波段和近红外波段的反射率。
其中NIR和SWIR1分别代表Sentinel2遥感影像的近红外波段和短波红外1波段的反射率。
其中NIR和RE2分别代表Sentinel2遥感影像的近红外波段和红边2波段的反射率。
其中Red和SWIR1分别代表Sentinel2遥感影像的红光波段和短波红外1波段的反射率。
其中Red和SWIR1分别代表Sentinel2影像的红光波段和短波红外1波段的反射率。
其中RE1、RE2、RE3和Red分别代表Sentinel2影像的红边1波段、红边2波段、红边3波段和红光波段的反射率。
然后采用10天中值合成的方式等间隔生成目标区域高质量Sentinel2遥感影像时间序列，使用线性插值的方式对10天中值数据集中的空缺像素进行填补，最后应用SG滤波处理插值后的时序曲线得到最终的Sentinel2遥感影像时间序列。
在GEE中获取Sentinel1 GRD地距影像数据，时间范围为2019年5月1日-10月10日，空间范围为九三垦区的矢量边界。
利用局部入射角对原始后向散射系数值进行入射角归一化，减弱入射角的影响；应用7×7精炼Lee滤波对Sentinel1影像中的噪点进行滤除。
通过波段筛选选择Sentinel1影像的VV和VH两个极化通道后向散射系数；计算Sentinel1影像的三个极化通道组合，分别是交叉共极化之和、交叉共极化比率和雷达植被指数RVI。
其中Sentinel1影像的极化通道组合方式如下：
交叉共极化之和σVH 0+σVV 0
交叉共极化比率σVH 0/σVV 0
雷达植被指数
其中σVH 0和σVV 0分别为Sentinel1影像的VH和VV通道的后向散射系数值。
然后采用10天中值合成的方式等间隔生成目标区域高质量Sentinel1影像时间序列，使用线性插值的方式对10天中值数据集中的空缺像素进行填补，最后应用SG滤波处理插值后的时序曲线得到最终的Sentinel1影像时间序列。
步骤S3、提取Sentinel2影像时间序列特征；
通过计算统计值、二阶谐波拟合以及物候估计等方法，提取Sentinel2影像的七个原始波段和五个植被指数的统计特征、生长高峰期特征、谐波拟合特征以及物候特征，作为大豆在Sentinel2影像中的时序特征。所述Sentinel2统计特征包括七个原始波段和五个植被指数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数(分位数是将一组数据从小到大排序，并计算相应的累计百分位，则某一百分位所对应数据的值就称为这一百分位的百分位数。可表示为：一组n个观测值按数值大小排列。如，处于p％位置的值称第p百分位数)。所述Sentinel2生长高峰期特征为特征参数在大豆生长高峰期时的特征值，将EVI或LSWI达到峰值的时间作为大豆生长高峰期。所述Sentinel2谐波拟合特征是通过对EVI时序曲线做二阶谐波拟合得到的拟合系数。所述Sentinel2物候特征是采用阈值法从EVI时序曲线中估算得到的生长季起始日SOS、生长季结束日EOS和生长季长度LOS(生长季起始日SOS至生长季结束日EOS之间的总时长)。
通过计算统计值、主成分分析等方法，提取Sentinel1影像的两个极化通道后向散射系数和三个极化通道组合的统计特征和主成分特征，作为大豆在Sentinel1影像中的时序特征。
所述二阶谐波拟合模型计算公式如下：
其中f为拟合值，a0、ak和bk分别为谐波余项、余弦系数和正弦系数，a0实际上为时序的平均值。n为谐波拟合的阶次，w为频率取1.5，t为当前日期在一年中的位置(例如当前日期为当年的第N天，则t＝N/365)，为0到1之间的小数。
通过计算统计值、主成分分析等方法，提取Sentinel1影像的两个极化通道后向散射系数和三个极化通道组合的统计特征和主成分特征，作为大豆在Sentinel1影像中的时序特征。所述Sentinel1统计特征包括特征参数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数。所述Sentinel1主成分特征为每个特征参数在大豆生长期内进行主成分分析后的前三个主成分分量。
所述主成分分析步骤如下：
①构建时序影像数据矩阵X
其中m为影像时间序列中的影像数量，n为每幅影像中的像元数。矩阵中第i行向量对应影像序列的第i幅影像，由影像的上下相邻两行数据首尾依次相接得到。
②根据原始图像数据矩阵X，计算协方差矩阵S
其中 为第i幅影像的均值。
③计算S矩阵的特征值和特征向量(均为m×1的列向量)。
④构建变换矩阵T
将特征值按其大小排列，排序后特征值依次为λ1≥λ2≥…λm，求得各特征值对应的单位特征向量：Uj＝[u1j，u2j，...，umj]T。然后以各特征向量为列构成矩阵U，U即为[U1，U2，...，Um]＝[uij]m×m，U的转置矩阵即为所求的PCA变换的变换矩阵T。
⑤计算主分量
利用主成分变换矩阵计算Y
其中矩阵Y的每一行即为一个主成分分量，Yj＝[yj1，yj2，…，yjn]为第j主成分分量。
在GEE(Google Earth Engine)中获取ESA 10m地表覆盖数据集，从中提取出耕地像素作为掩膜处理上述特征图，排除非耕地像素的影响。
步骤S4、构建随机森林分类模型并完成模型训练
在GEE中调用ee.Classifier.smileRandomForest()函数构建随机森林分类模型，设置的参数包括：决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数，其中将决策树的数量设置为100，叶子节点最少样本数设为10，其他参数保持默认值。将实地采样的作物样点数据集上传到GEE平台中，并将其划分为训练集和验证集。将步骤S3中的所有特征参数汇总到一起，然后使用sampleRegions()函数将所有特征提取到训练样本中，将训练样本输入到构建好的随机森林分类模型中，完成模型训练。
步骤S5、完成作物分类并提取大豆种植区，计算大豆种植面积
将带有步骤S3中的所有特征参数的样点数据集输入到训练好的随机森林分类模型中，输出得到大豆、玉米和水稻的分类图，三者的像素值分别是2、1、0。利用eq()函数提取出大豆像素，得到大豆种植区分布图。利用验证集样本计算分类混淆矩阵，通过混淆矩阵计算出分类总体精度为0.9471，kappa系数为0.9326，大豆F1得分为0.9326。通过统计大豆种植区分布图中所有像素的面积，得到大豆种植面积为2574.08千公顷，面积统计精度为92.13％。2019年九三垦区大豆种植区分布图如图3所示。
本申请基于Sentinel1/2遥感影像时间序列，充分提取大豆生长期关键特征和物候特征，能够对大豆种植区进行精确识别提取，减弱云雨天气对单独使用光学影像进行提取所造成的不良影响，使大豆种植区的提取精度得到提升。本申请能够实现市级或省级尺度的大豆种植区空间分布的精确提取，并计算区域内的大豆种植面积，为大豆作物的种植管理、产量估算和市场调动提供重要科学指导依据。
凡是本申请技术特征和技术方案的简单变形或者组合，应认为落入本申请的保护范围。
Claims (8)
1.一种大豆种植区识别和面积测算方法，其特征是包括以下步骤：
步骤S1、确定目标区域范围、作物种类、样点数据集及作物生长期：
步骤S2、获取目标区域的Sentinel2影像并进行预处理，得到Sentinel2影像时间序列；获取目标区域的Sentinel1影像并进行预处理，得到Sentinel1影像时间序列；
步骤S3、提取Sentinel2影像的统计特征、生长高峰期特征、谐波拟合特征以及物候特征作为大豆在Sentinel2影像中的时序特征；提取Sentinel1影像的两个极化通道后向散射系数和三个极化通道组合的统计特征和主成分特征作为大豆在Sentinel1影像中的时序特征；在GEE中获取ESA 10m地表覆盖数据集，去除步骤S2、步骤S3中所有不在地表覆盖数据集范围内的数据，得到掩膜处理后的目标区域数据；
步骤S4、将样点数据集上传到GEE平台中，并将其划分为训练集和测试集；将步骤3中特征参数提取到训练集样点中，构建随机森林分类模型，并将带有步骤S3所述特征参数的训练集输入到随机森林分类模型中完成模型训练；
步骤S5、将目标区域数据输入到训练好的随机森林分类模型中，得到作物分类结果并从中提取出大豆种植区，计算大豆种植面积。
2.根据权利要求1所述大豆种植区识别和面积测算方法，其特征是：获取目标区域的Sentinel2影像并进行预处理包含以下步骤：
在GEE平台中获取大豆生长期内的所有Sentinel2卫星地表反射率数据影像，将所有波段数值除以10000得到真实反射率数值；通过波段筛选选择Sentinel2影像的七个原始波段，分别是蓝、红、红边1、红边2、红边3、近红外和短波红外1，计算Sentinel2影像的五个植被指数，分别是增强型植被指数EVI、地表水分指数LSWI、归一化植被衰老指数NDSVI、红边归一化植被指数RENDVI和红边位置REP；利用云得分算法识别每幅影像中的含云像素并将之去除；然后采用10天中值合成的方式等间隔生成目标区域Sentinel2影像时间序列，使用线性插值的方式对10天中值数据集中的空缺像素进行填补，应用SG滤波处理Sentinel2影像时间序列；
其中Sentinel2影像的植被指数计算公式如下：
其中，Red、Blue、NIR分别代表Sentinel2影像的红光波段、蓝光波段和近红外波段的反射率；
其中NIR和SWIR1分别代表Sentinel2影像的近红外波段和短波红外1波段的反射率；
其中NIR和RE2分别代表Sentinel2影像的近红外波段和红边2波段的反射率；
其中Red和SWIR1分别代表Sentinel2影像的红光波段和短波红外1波段的反射率；
其中RE1、RE2、RE3和Red分别代表Sentinel2影像的红边1波段、红边2波段、红边3波段和红光波段的反射率。
3.根据权利要求1或2所述大豆种植区识别和面积测算方法，其特征是：获取目标区域的Sentinel1影像并进行预处理包含以下步骤：
在GEE平台中获取大豆生长期内的所有Sentinel1的GRD地距影像，成像方式为IW干涉宽带扫描；通过波段筛选选择Sentinel1影像的两个极化通道后向散射系数，分别是垂直同向极化和垂直/水平交叉极化，即VV和VH；计算Sentinel1影像的三个极化通道组合，分别是交叉共极化之和、交叉共极化比率和雷达植被指数RVI；利用局部入射角对原始后向散射系数值进行入射角归一化；应用7×7精炼Lee滤波对SAR影像中的噪点进行滤除；然后采用10天中值合成的方式等间隔生成目标区域Sentinel1影像时间序列，使用线性插值的方式对10天中值数据集中的空缺像素进行填补，应用SG滤波处理Sentinel1影像时间序列；
其中Sentinel1影像的极化通道组合方式如下：
交叉共极化之和σVH 0+σVVV 0
交叉共极化比率σVH 0/σVV 0
雷达植被指数
其中σVH 0和σVV 0分别为Sentinel1影像的VH和VV通道的后向散射系数值。
4.根据权利要求1所述大豆种植区识别和面积测算方法，其特征是：所述Sentinel2的统计特征包括特征参数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数；所述Sentinel2生长高峰期特征为大豆生长高峰期时的特征参数值，将EVI或LSWI达到峰值的时间作为大豆生长高峰期；所述Sentinel2谐波拟合特征是通过对EVI时序曲线做二阶谐波拟合得到的拟合系数；所述Sentinel2物候特征采用阈值法从EVI时序曲线中估算得到的生长季起始SOS、生长季结束EOS和生长季长度LOS。
5.根据权利要求4所述大豆种植区识别和面积测算方法，其特征是：所述二阶谐波拟合模型计算公式如下：
其中f为拟合值，a0、ak和bk分别为谐波余项、余弦系数和正弦系数，n为谐波拟合的阶次，w为频率取1.5，t为当前日期在一年中的位置。
6.根据权利要求1所述大豆种植区识别和面积测算方法，其特征是：所述Sentinel1统计特征包括特征参数在大豆生长期中的最大值、最小值、方差、15/50/90th分位数；所述Sentinel1主成分特征为每个特征参数在大豆生长期内进行主成分分析后的前三个主成分分量；
所述主成分分析步骤如下：
①构建时序影像数据矩阵X
其中m为影像时间序列中的影像数量，n为每幅影像中的像元数；矩阵中第i行向量对应影像序列的第i幅影像，由影像的上下相邻两行数据首尾依次相接得到；
②根据原始图像数据矩阵X，计算协方差矩阵S
其中 为第i幅影像的均值；
③计算S矩阵的特征值和特征向量；
④构建变换矩阵T
将特征值按其大小排列，排序后特征值依次为λ1≥λ2≥…λm，求得各特征值对应的单位特征向量：Uj＝[u1j，u2j，...，umj]T，然后以各特征向量为列构成矩阵U，U即为[U1，U2，...，Um]＝[uij]m×m，U的转置矩阵即为所求的PCA变换的变换矩阵T；
⑤计算主分量
利用主成分变换矩阵计算Y
其中矩阵Y的每一行即为一个主成分分量，Yj＝[yj1，yj2，...，yjn]为第j主成分分量。
7.根据权利要求1所述大豆种植区识别和面积测算方法，其特征是：所述构建随机森林分类模型通过调用GEE云平台的随机森林函数完成，设置的参数包括：决策树的数量、最大特征数、决策树深度、叶子节点最少样本数、节点划分的最小样本数，其中将决策树的数量设置为100，叶子节点最少样本数设为10，其他参数保持默认值。
8.根据权利要求1所述大豆种植区识别和面积测算方法，其特征是：所述作物为大豆、玉米和水稻，所述步骤S5包括：将带有步骤S3中的所有特征参数的目标区域数据输入到训练好的随机森林分类模型中，输出得到大豆、玉米和水稻的分类图，大豆、玉米、水稻的像素值分别是2、1、0，利用eq()函数提取出像素值为2的目标区域，得到大豆种植区。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310522156.5A CN116645603A (zh) | 2023-05-10 | 2023-05-10 | 一种大豆种植区识别和面积测算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310522156.5A CN116645603A (zh) | 2023-05-10 | 2023-05-10 | 一种大豆种植区识别和面积测算方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116645603A true CN116645603A (zh) | 2023-08-25 |
Family
ID=87623893
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310522156.5A Pending CN116645603A (zh) | 2023-05-10 | 2023-05-10 | 一种大豆种植区识别和面积测算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116645603A (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117607063A (zh) * | 2024-01-24 | 2024-02-27 | 中国科学院地理科学与资源研究所 | 一种基于无人机的森林垂直结构参数测量系统和方法 |
| CN117636152A (zh) * | 2023-10-23 | 2024-03-01 | 广东省国土资源测绘院 | 一种面向多云雨区烟草遥感监测方法及系统 |
| CN117649598A (zh) * | 2023-10-23 | 2024-03-05 | 广东省国土资源测绘院 | 基于sar影像的近海养殖空间分布信息监测方法及系统 |
| CN117690017A (zh) * | 2023-11-16 | 2024-03-12 | 宁波大学 | 一种顾及物候时序特征的单、双季水稻提取方法 |
| CN117726902A (zh) * | 2023-12-19 | 2024-03-19 | 中国科学院空天信息创新研究院 | 一种基于Sentinel2时序曲线分布的大豆样本生成方法和系统 |
-
2023
- 2023-05-10 CN CN202310522156.5A patent/CN116645603A/zh active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117636152A (zh) * | 2023-10-23 | 2024-03-01 | 广东省国土资源测绘院 | 一种面向多云雨区烟草遥感监测方法及系统 |
| CN117649598A (zh) * | 2023-10-23 | 2024-03-05 | 广东省国土资源测绘院 | 基于sar影像的近海养殖空间分布信息监测方法及系统 |
| CN117690017A (zh) * | 2023-11-16 | 2024-03-12 | 宁波大学 | 一种顾及物候时序特征的单、双季水稻提取方法 |
| CN117726902A (zh) * | 2023-12-19 | 2024-03-19 | 中国科学院空天信息创新研究院 | 一种基于Sentinel2时序曲线分布的大豆样本生成方法和系统 |
| CN117607063A (zh) * | 2024-01-24 | 2024-02-27 | 中国科学院地理科学与资源研究所 | 一种基于无人机的森林垂直结构参数测量系统和方法 |
| CN117607063B (zh) * | 2024-01-24 | 2024-04-19 | 中国科学院地理科学与资源研究所 | 一种基于无人机的森林垂直结构参数测量系统和方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116645603A (zh) | 一种大豆种植区识别和面积测算方法 | |
| CN110472184B (zh) | 一种基于Landsat遥感数据的多云雨雾地区水稻识别方法 | |
| CN110796001B (zh) | 一种卫星影像覆膜农田识别与提取方法及系统 | |
| CN101692037B (zh) | 高光谱图像和独立分量分析植物叶面叶绿素分布的方法 | |
| CN112800973A (zh) | 一种基于植被物候特征决策的互花米草提取方法 | |
| CN112836725A (zh) | 基于时序遥感数据的弱监督lstm循环神经网络稻田识别方法 | |
| CN115271217A (zh) | 一种基于无人机多源遥感数据的小麦产量预测方法 | |
| CN113205014B (zh) | 一种基于图像锐化的时序数据耕地提取方法 | |
| Rodrigues et al. | Monitoring vegetation dynamics inferred by satellite data using the PhenoSat tool | |
| CN112577906A (zh) | 城市绿地土壤含水率检测方法 | |
| CN112669363A (zh) | 城市绿地三维绿量测算方法 | |
| CN116665073A (zh) | 一种基于多源数据的玉米产量遥感估测方法 | |
| CN117011694A (zh) | 一种基于级联循环网络的林木生长参数预测方法 | |
| CN115994939A (zh) | 一种基于地面激光点云的树木叶面积估算方法 | |
| CN115359365A (zh) | 一种面向滨海潮间带盐地碱蓬识别和覆盖度估算方法 | |
| Wittamperuma et al. | Remote-sensing-based biophysical models for estimating LAI of irrigated crops in Murry darling basin | |
| CN112577954A (zh) | 城市绿地地上生物量估测方法 | |
| CN115063610B (zh) | 基于Sentinel-1、2影像的大豆种植区识别方法 | |
| Medvedeva et al. | Determination of area of drought-affected crops based on satellite data (exemplified by crops in Chuvashia in 2010) | |
| CN113505635B (zh) | 基于光学和雷达的冬小麦与大蒜混种区识别方法及装置 | |
| CN115019205B (zh) | 一种基于无人机多光谱影像的油菜花期spad和lai估测方法 | |
| CN113887493B (zh) | 一种基于id3算法的黑臭水体遥感影像识别方法 | |
| Janwale et al. | Automatic estimation of nitrogen content in cotton (Gossypium hirsutum L.) plant by using image processing techniques: A review | |
| Wang et al. | Evaluation and Mapping of Rice Flood Damage Using Domestic Remotely Sensed Data in China | |
| Hachicha et al. | Forecasting of the normalized difference vegetation index time series in Jbeniana |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |