CN116824384A - 一种基于标准曲线的大豆识别方法 - Google Patents
一种基于标准曲线的大豆识别方法 Download PDFInfo
- Publication number
- CN116824384A CN116824384A CN202310197975.7A CN202310197975A CN116824384A CN 116824384 A CN116824384 A CN 116824384A CN 202310197975 A CN202310197975 A CN 202310197975A CN 116824384 A CN116824384 A CN 116824384A
- Authority
- CN
- China
- Prior art keywords
- soybean
- sample
- standard curve
- sample point
- vegetation index
- 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 166
- 235000010469 Glycine max Nutrition 0.000 title claims abstract description 166
- 238000000034 method Methods 0.000 title claims abstract description 51
- 239000011159 matrix material Substances 0.000 claims abstract description 16
- 238000012545 processing Methods 0.000 claims abstract description 14
- 238000011156 evaluation Methods 0.000 claims abstract description 5
- 238000011160 research Methods 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 14
- 230000032683 aging Effects 0.000 claims description 11
- 238000007637 random forest analysis Methods 0.000 claims description 7
- 230000002194 synthesizing effect Effects 0.000 claims description 6
- 238000010586 diagram Methods 0.000 claims description 3
- 230000005059 dormancy Effects 0.000 claims description 3
- 230000035800 maturation Effects 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims 1
- 239000000523 sample Substances 0.000 description 127
- 240000008042 Zea mays Species 0.000 description 10
- 235000005824 Zea mays ssp. parviglumis Nutrition 0.000 description 10
- 235000002017 Zea mays subsp mays Nutrition 0.000 description 10
- 235000005822 corn Nutrition 0.000 description 10
- 230000008859 change Effects 0.000 description 7
- 238000002310 reflectometry Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 5
- 239000013074 reference sample Substances 0.000 description 5
- 238000012544 monitoring process Methods 0.000 description 4
- 238000001556 precipitation Methods 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 238000002834 transmittance Methods 0.000 description 4
- 230000002159 abnormal effect Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 235000013305 food Nutrition 0.000 description 2
- 239000002689 soil Substances 0.000 description 2
- 238000009331 sowing Methods 0.000 description 2
- 208000003643 Callosities Diseases 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 206010020649 Hyperkeratosis Diseases 0.000 description 1
- 235000010716 Vigna mungo Nutrition 0.000 description 1
- 240000001417 Vigna umbellata Species 0.000 description 1
- 235000011453 Vigna umbellata Nutrition 0.000 description 1
- 241000607479 Yersinia pestis Species 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 239000008157 edible vegetable oil Substances 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 244000144972 livestock Species 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 102000004169 proteins and genes Human genes 0.000 description 1
- 108090000623 proteins and genes Proteins 0.000 description 1
- 238000012950 reanalysis Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 230000007958 sleep Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 239000010981 turquoise Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Landscapes
- Image Processing (AREA)
Abstract
本申请涉及一种基于标准曲线的大豆识别方法，包括：步骤S1、卫星数据处理；步骤S2、构建时间序列归一化物候指数特征：计算植被指数获得植被指数图像，对植被指数图像使用时间序列构建时间序列植被指数特征；步骤S3、构建样本点标准曲线：依据样本点时间序列植被指数特征、位置、气象数据因子形成样本点标准曲线；步骤S4、样本点标准曲线参数信息提取与物候特征计算：样本点标准曲线中主要物候信息提取；对其它像元计算生成物候参数特征影像；步骤S5、农作物品种识别与精度评估：将样本信息输入到分类器中得到样本分类器，其它像元数据输入到样本分类器中进行识别，使用分类混淆矩阵计算识别精度。本申请自动化程度高且精度高。
Description
技术领域
本申请涉及一种基于标准曲线的大豆识别方法，主要适用于大豆种植面积的评估和监测。
背景技术
大豆是高蛋白食品、牲畜饲料的主要原料以及食用油的重要来源，在世界粮食生产中占有重要地位。传统的农业调查方式估算大豆种植面积相对比较费时费力，易受主观因素影响，结果数据亦无法提供具体空间分布信息。遥感技术可以更及时、高效和客观地实现大规模的农作物种植面积监测，且成本低廉。目前主流的遥感技术识别大豆的方法有四类：
(1)使用光谱特征识别大豆，即使用卫星数据的各个波段作为其分类特征识别大豆，也有人发现卫星数据如Sentinel-2数据的红边波段和近红外波段可以很好地区分大豆和玉米。该方法使用的是单时相影像，当研究区域较大(省级尺度以上)时，由于云雨等污染Sentinel-2单时相影像缺失不能完全覆盖整个研究区，则无法开展后续研究，单时相影像也不能体现作物的变化过程，尤其是在同一时期具有相似反射率值的作物容易混淆导致无法识别，并且由于卫星传感器、镜头异常等波段原始反射值可能存在一定程度的噪声导致产生误差。
(2)使用植被指数特征识别大豆，即使用卫星数据的各个波段计算植被指数识别大豆，如增强型植被指数EVI、归一化差异植被指数NDVI等。该方法使用的是单时相影像，当研究区域较大(省级尺度以上)时，由于云雨等污染Sentinel-2单时相影像缺失不能完全覆盖整个研究区，则无法开展后续研究。单时相影像也不能体现作物的变化过程，尤其是在同一时期具有相似反射率值的作物容易混淆导致无法识别。
(3)使用时间序列特征识别大豆，如使用时间序列卫星数据光谱波段，时间序列植被指数等，时间序列方法有计算百分位数合成一幅特征影像等，即将一段时间内的多幅连续卫星影像数据处理成一幅时间序列影像数据，能表现随时间变化特征。该方法仅处理卫星数据和植被指数，对大豆的生长特点不能可视化分析，无法找到区分大豆的关键生育期，不具备可解释性，不能解决样本获取困难问题，相关参数少，大豆识别精确度不高。
(4)使用物候特征识别大豆，通过对时间序列植被指数进行处理获得关键参数即为物候特征。该方法只用了几个关键物候点，不具备普适性，不能解决样本获取困难问题，不能适应于不同气候区，不能在没有样本的区域识别大豆，采用的技术参数较少，大豆识别精确度不高。
发明内容
本申请解决的技术问题是克服现有技术中存在的上述不足，而提供一种能够精确识别大豆且具有一定预测功能的基于标准曲线的大豆识别方法。
本申请解决上述技术问题所采用的技术方案是：一种基于标准曲线的大豆识别方法，其特征是包括以下步骤
步骤S1、卫星数据处理：
所述卫星数据处理步骤包括获取研究区内大豆生长期所有相关卫星数据，对卫星数据进行云量控制和去云处理。
步骤S2、构建时间序列归一化物候指数特征：
使用处理好的卫星数据计算研究区内所有像元的植被指数获得所有像元的植被指数图像，对每一个像元使用该像元的时间序列方法处理该像元的植被指数图像构建该像元的时间序列植被指数特征。本步骤处理的基本单位为卫星数据的像元，所述计算处理方法是逐像元计算。
步骤S3、构建样本点标准曲线：
构建样本点标准曲线有两种方法，一种是直接采用上一年度该样本点标准曲线或历年该样本点平均值；另一种方法是：获取包括大豆在内的农作物典型代表样本点(简称样本点)位置，使用S2步骤获得的样本点时间序列植被指数特征与样本点位置计算样本点时间序列曲线，使用样本点时间序列曲线并添加气象数据因子构建样本点标准曲线。
步骤S4、样本点标准曲线参数信息提取与物候特征计算：
将样本点标准曲线中的大豆样本点关键生长期信息提取出来，得到样本点的主要物候信息(样本信息)；对研究区范围内除样本点以外的所有像元(简称其它像元)计算大豆样本点关键生长期信息(指9个主要物候特征参数)，合成物候参数特征影像。
步骤S5、农作物品种识别与精度评估：
将样本信息输入到分类器中得到样本分类器，将样本点以外所有像元的物候参数特征影像输入到样本分类器中进行识别并输出识别结果，使用分类混淆矩阵计算识别精度。
所述步骤S1中对卫星数据进行云量控制是指仅采用(筛选)整体云量小于30％的影像；所述去云处理是指去除异常值，例如使用NA赋值给云、雪、阴影中的像元，达到去除异常值的效果。
所述步骤S2中，采用的植被指数为归一化物候植被指数，因其对土壤不敏感以及包含在归一化物候植被指数中的短波红外波段的反射率可以捕捉到树冠含水量的变化，使得归一化物候植被指数在大豆识别中效果较好，所述的归一化物候植被指数计算公式如下所示：
式中，NIR为近红外波段反射率，RED为红色波段反射率，SWIR为短波红外波段反射率。
所述步骤S2中，时间序列方法为谐波拟合，线性谐波模型使用一系列正弦和余弦波拟合复杂信号，已广泛用于植被指数时间序列的拟合，所述的线性谐波模型计算公式如下所示：
其中，f(t)是第t时刻拟合的植被指数值；a是常数项；b是一阶项的系数；M是谐波组合的数量，在本申请中，M取2，即二项线性谐波模型；C和D是余弦函数和正弦函数的系数；ω是一年中天数的倒数(1/365)；t是DOY(一年中的某一天)。
所述步骤S3中，使用S2步骤获得的时间序列植被指数特征与农作物典型代表样本点位置计算农作物典型代表样本点的时间序列曲线，即样本点位置与样本点时间序列植被指数特征对应起来统一在一起，可视化该位置的时间序列植被指数像元的变化情况，形成一条样本点时间序列曲线，横坐标为儒略日DOY(天)，纵坐标为归一化物候植被指数的值。
所述步骤S3中，使用样本点时间序列曲线并添加气象数据因子构建样本点标准曲线，即使用气象数据添加到样本点时间序列曲线中，得到该像元的(有可能是大豆，也有可能是玉米或其它)标准曲线，样本点标准曲线公式为：
其中，部分为谐波拟合公式，Temp为温度变量，hour为日照时长变量，rain为降雨变量，c0，c1，c2为气象权重因子。
所述步骤S4中，对样本点标准曲线中的大豆样本点关键生长期信息进行提取，得到样本信息，大豆样本点关键生长期信息采用主要物候点参数信息(9个参数)，9个参数分别为大豆生长期的三个节点的儒略日DOY值，衰老期的三个节点儒略日DOY值，以及对应三个大豆生长阶段的长度。
主要物候点(时间点)通过下述公式来确定：
u＝(Vmin-Vmax)×p+Vmin
其中，u为参考阈值，用于确定大豆样本点物候开始时间(SOS)和大豆样本点物候结束时间(EOS)，Vmin和Vmax分别为大豆样本点时间序列中的最小值和最大值，p为振幅比例，p有三个值，第一个值对应大豆样本点生长开始节点、大豆样本点休眠开始节点，第二个值对应大豆样本点生长中期节点、大豆样本点衰老中期节点，第三个值对应大豆样本点成熟期节点、大豆样本点衰老开始节点。p在本申请中分别取0.15、0.5和0.90，于是得到6个主要物候点信息，分别为大豆样本点生长期的三个节点的儒略日DOY值，衰老期的三个节点儒略日DOY值。
本申请将大豆样本点物候开始时间(SOS)和大豆样本点物候结束时间(EOS)分别指定为超过阈值u的第一天和最后一天，u为一个动态值，该值取决于大豆样本点时间序列的年内振幅。
SOS＝first(F(x)>u)
EOS＝last(f(x)>u)
其中，f(t)为步骤S3中大豆样本点标准曲线的公式，first()为大豆样本点超过阈值的第一个值所对应的儒略日DOY，last()为大豆样本点超过阈值的最后一个值所对应的儒略日DOY。
并且计算p分别取0.15、0.5和0.90时，大豆样本点物候开始时间(SOS)和大豆样本点物候结束时间(EOS)的长度length，长度计算公式为：
length＝EOS-SOS
其中，EOS为大豆样本点物候开始的儒略日DOY，SOS为大豆样本点物候结束的儒略日DOY，length为大豆样本点生长阶段时间长度。
所述步骤S4中，对研究区范围除样本点以外所有像元计算关键生长期信息，并合成物候参数特征影像的方法为：计算研究区除样本点以外所有像元在生长条件最接近的大豆样本点的主要物候点时的物候参数特征影像，共9幅特征影像作为物候参数特征。
所述步骤S5中，某像元的物候参数特征影像输入到样本分类器中进行识别时，采用生长条件最接近的样本点数据进行比对识别(例如先与生长条件最接近的大豆样本点进行比对，该像元数据与该大豆样本点数据相似度达标则标记为大豆否则如果需要再与大豆参照物进行比对如果不需要则直接标记为其它)，输出分类结果，提取出大豆并使用分类混淆矩阵评估精度，评估指标为总体精度(0A)，制图精度(PA)，用户精度(UA)。
本申请与现有技术相比，具有以下优点和效果：使用方便，自动化程度高，能够精确识别大豆。
附图说明
图1是本申请实施例的流程示意图。
图2是大豆2020年标准曲线(预测曲线)与实际的2020年时序曲线比较示意图。
图3是大豆归一化物候植被指数灾害曲线与大豆标准曲线对比图。
图4是大豆归一化物候植被指数及关键物候点设计原理示意图。
具体实施方式
下面结合附图并通过实施例对本申请作进一步的详细说明，以下实施例是对本申请的解释而本申请并不局限于以下实施例。
参见图1～图4，本申请主要适用于大豆种植面积的评估和监测，通过充分利用资源卫星遥感数据结合本申请的大豆识别方法来判断某区块(由众多的像元构成该区块的影像，如哨兵2号影像的像元空间分辨率是10m，每一个像元都是10m*10m大小，每个像元都有相对独立的参数，本申请涉及到的计算都是逐像元计算，即对组成影像的每个像元都单独计算)农产品种植种类是否为大豆，大豆种植面积的多少以及分布情况，为科学预测和监测大豆种植面积、种类、生长状况、产量大小、病虫害情况等提供基础信息。下面以黑龙江省大豆识别为例进行说明。
黑龙江省大豆播种面积和产量分别占全国大豆播种面积和产量的48.90％与46.95％，是我国大豆主产区，将黑龙江的大豆作为大豆识别方法的对象具有较强的代表性。
步骤S1、卫星数据处理；
所述卫星数据为哨兵2号卫星数据(Sentine1-2 MSI：MultiSpectralInstrument，Level-2A)，该数据来源于Google Earth Engine云平台，该数据是已经完成大气校正的地表反射率产品，无需进行预处理步骤，对2019年到2020年所有可获得影像使用云量筛选函数筛选云量小于30％的影像，并使用哨兵2号卫星数据特有的去云函数去掉异常值，最后得到一共7026幅卫星影像，完整覆盖整个黑龙江省大豆生长季节全过程。
步骤S2、构建时间序列归一化物候指数(简称时序归一化物候植被指数)；
首先使用处理好的卫星数据计算归一化物候植被指数合成为植被指数图像，归一化物候植被指数计算公式如下：
式中，NIR为近红外波段反射率，RED为红色波段反射率，SWIR短波红外波段反射率。
然后使用谐波拟合方法计算时序归一化物候植被指数(指带有时间序列特征的归一化物候植被指数)，谐波拟合公式如下：
其中，f(t)是第t时刻拟合的植被指数值；a是常数项；b是一阶项的系数；M是谐波组合的数量，在本申请中，M取2，即二项线性谐波模型；C和D是余弦函数和正弦函数的系数；ω是一年中天数的倒数(1/365)；t是DOY(一年中的某一天)。
步骤S3、构建样本点标准曲线：
首先是获取大豆及作为大豆参照物的农作物(生长季节与大豆相近的农作物，例如玉米和/或大米，大豆参照物可以有但非必须)典型代表样本点，每种农作物的样本点均由若干个像元组成且仅种植一种农作物。本实例使用的实地采集样本点数据主要为黑龙江省内的2019年315个大豆样本点和498个玉米样本，362个其它作物样本点，以及2020年1000个大豆样本点。本申请样本点某个生长季完整数据获得后可以长年使用，无需重复采样，也可以调整为最近一年的数据或若干年的平均数据；特殊情况下(例如新建大块的大豆田或者某块样本点改为其它用途)也可相应调整样本点。然后是在GEE云平台上使用S2步骤的时间序列归一化物候植被指数特征与大豆等(这里的等指大豆参照物，下同)作物典型代表样本点位置计算样本点的时间序列曲线，绘图使用了ui.Chart.image.doySeriesByRegion()函数。最后是使用时间序列曲线并添加气象数据因子构建样本点标准曲线，气象数据使用的是气象数据ERA5再分析数据(ERA5-Land Hourly-ECMWF Climate Reanalysis)中的气温和降水数据，对每一天的降水数据累加24小时得到每日降水量，对每一天的气温数据取平均值得到每日气温，根据太阳高度角、方位角和大气透过率等数据计算日照时长，将气温、降水量与日照时长三个变量加入到时间序列处理谐波拟合中，得到样本点标准曲线(简称标准曲线，包括大豆样本点标准曲线、玉米样本点标准曲线，还可以包含其它作物)。
此外，本实例为了证明获得的标准曲线具有鲁棒性，从时间尺度和灾害区域两个方面对曲线进行比较验证，首先是时间尺度上的验证，以2019年为参考年，2020年为目标年，使用2019年的归一化物候植被指数拟合时间序列曲线，在此基础上使用2020年ERA5气象数据添加气候因子，最后得到预测的2020年标准曲线。然后使用实际2020年的归一化物候植被指数拟合时间序列曲线，记为实际的2020年曲线。标准曲线与实际曲线对比图如图2所示，二者曲线十分接近，相比之下，实际曲线在某些时间段有一定的上下抖动趋势。该图表示，标准曲线能通过添加目标年的气候数据得到较好的预测效果，与目标年份的实际中值拟合曲线十分相似。然后是灾害区域的验证，将2020年实际曲线添加2019年洪涝灾害气象因子的标准曲线与2019年处于洪涝灾害的曲线进行对比，对比图如图3所示，二者在变化趋势上基本一致，而标准曲线因受到变化较大气象数据影响，在数值方面高于灾害曲线，但二者在曲线上升阶段和下降阶段十分一致，即在关键大豆生长阶段二者是十分接近的，说明标准曲线方法具有一定的鲁棒性。
所述标准曲线具体公式如下：
其中，部分为谐波拟合公式，Temp为温度变量，hour为日照时长变量，rain为降雨变量，c0，c1，c2为气象权重因子。
所述日照时长计算公式如下：
h＝∑(τi＞τmin)×Δti
其中，∑表示对每个时刻的值进行求和；τi表示第i个时刻的大气透过率；τmin表示阈值大气透过率，即太阳能够被视为“直射”到地表的最小透过率，通常为0.1～0.3；Δti表示第i个时刻的时间间隔。
步骤S4、样本点标准曲线参数信息提取与物候特征计算；对样本点标准曲线中的关键生长期信息(本实施例取大豆的6个主要物候点信息及生长期长度、生长中期长度、成熟期长度)进行提取，得到样本点的大豆或大豆参照物主要物候信息(简称样本信息或样本)；对研究区范围除样本点以外所有像元计算关键生长期信息，并合成物候参数特征影像。
对样本点标准曲线中的关键生长期信息进行提取，得到样本点的大豆或大豆参照物作物信息的方法如下：将样本点标准曲线中的大豆关键生长期主要物候点参数信息提取出来作为样本，所述主要物候点分别为大豆生长期的三个节点的儒略日DOY值，大豆衰老期的三个节点儒略日DOY值，共计6个主要物候点以及对应的大豆三个生长阶段的长度，共计9个参数。大豆样本点数据(大豆样本)作为下一步骤某像元分类为大豆的依据，所有大豆参照物(例如玉米)样本点数据(大豆参照物样本)作为下一步骤某像元分类为大豆参照物(例如玉米)的依据，以相似度是否达标为准，既与大豆样本不相似又与大豆参照物样本不相似的农作物分类为其它。
主要物候点确定的公式为：
u＝(Vmim-Vmax)×p+Vmin
其中，u为参考阈值，用于确定大豆样本点物候开始时间(SOS，如图4所示生长开始时间)和物候结束时间(EOS，如图4中休眠期开始时间)，Vmin和Vmax分别为大豆样本点时间序列中的最小值和最大值(这两个值为像元的时间序列中的最小值和最大值，在谷歌地球引擎(Google Earth Engine)云平台中使用.min()和.max()两个函数得到)，p为振幅比例，p在本申请中分别取0.15(对应大豆样本点生长开始节点、大豆样本点休眠开始节点)、0.5(对应大豆样本点生长中期节点、大豆样本点衰老中期节点)和0.90(对应大豆样本点成熟期节点、大豆样本点衰老开始节点)，参见图4，于是得到6个主要物候点信息，分别为生长期的三个节点的儒略日DOY值，衰老期的三个节点儒略日DOY值。
本申请将大豆样本点物候开始时间(SOS)和大豆样本点物候结束时间(EOS)分别指定为阈值u的第一天和最后一天，u为一个动态值，该值取决于大豆样本点时间序列的年内振幅。
SOS＝first(f(x)＞u)
EOS＝last(F(x)＞u)
其中，f(t)为步骤S3中大豆样本点标准曲线的公式，first()为大豆样本点超过阈值的第一个值所对应的儒略日DOY，last()为大豆样本点超过阈值的最后一个值所对应的儒略日DOY。
并且计算p分别取0.15、0.5和0.90时，大豆样本点物候开始时间(SOS)和大豆样本点物候结束时间(EOS)的长度length，所述长度计算公式为：
length＝EOS-SOS
其中，EOS为大豆样本点物候开始的儒略日DOY，SOS为大豆样本点物候结束的儒略日DOY，length为该大豆样本点生长阶段时间长度。
对研究区范围所有像元计算关键生长期信息并合成物候参数特征影像的方法如下：该步骤的方法与步骤S4中样本点参数信息提取方法相同，但应用对象不同，即计算研究区所有像元相对于生长条件最接近的样本点(通常采用距离最接近的样本点，一方面方便计算，另一方面相对而言土壤、气象条件比较接近，也可以采用某乡镇或某区块设立一个大豆样本点、一个大豆参照物样本点，某乡镇或某区块中所有像元均与这两个样本点数据进行比对)三个节点的儒略日DOY值，衰老期的三个节点儒略日DOY值，以及对应三个大豆样本点生长阶段的长度，每个像元的所有(9)参数合成一幅特征影像，即9参数的物候参数特征影像，该物候参数特征影像反应了该像元完整的作物信息，当其输入到样本分类器中与距离最接近的大豆样本点作物信息、大豆参照物样本点作物信息进行比对时，即可依据相似度最高原则识别出该像元种植的是大豆还是大豆参照物或者其它植物。
步骤S5、大豆提取(识别)与精度评估：采用随机森林分类器，树的棵数为500棵，其他参数使用默认值，使用步骤S4的样本输入到随机森林分类器中进行训练，得到具有所有样本信息的随机森林分类器(样本分类器)，将所有像元的大豆关键生长期主要物候特征输入到样本分类器进行识别，识别时采用生长条件最接近的样本点进行比对，若某像元数据与生长条件最接近的大豆样本点数据相似度达标(例如90％)则识别为大豆，若某像元数据与生长条件最接近的玉米样本点数据相似度达标(例如90％)则识别为玉米…识别结果通过颜色(本例仅大豆标注为绿色，其它植物均标注为土黄色，有需要的话也可以将玉米或其它农作物标注出来)标注在当年的黑龙江省大豆种植分布图中，同时计算随机森林分类混淆矩阵，得到总体精度(0A)、制图精度(PA)、用户精度(UA)和F1-Score四个指标，随机森林分类混淆矩阵定义如下所示：
公式(6)中，n表示类别数，mij表示实际属于i类的像元在分类结果图中被分到j类的像元数量，对角线元素的值即是各类别被正确分类的像元数，因此分类混淆矩阵中对角线元素的值越大，表示被正确分类的像元数越多，分类结果越可靠。利用分类混淆矩阵计算本申请总体精度(OA)、制图精度(PA)、用户精度(UA)和F1-Score。其计算公式分别为：
其中，mi+为分类混淆矩阵中的行总和，m+i为分类混淆矩阵中的列总和。
按照上述方法提取出的2020年黑龙江省大豆总体精度为86.95％，用户精度为90.91％，制图精度为86.14％，F1-Score为0.8846，并与统计数据相比，面积精度达到了95.94％，在大豆识别研究领域中识别效果较好。
本申请提出了一种基于标准曲线的大豆识别方法，该方法专注于植被指数时间序列曲线本身，通过可视化大豆与玉米等植被指数时间序列曲线，分析总结其曲线特征变化规律与差异，找到识别大豆的关键生长阶段时间序列曲线特征，提高大豆识别精度，并且通过添加气象权重因子，使标准曲线具有更强的鲁棒性，具有预测其他年份的标准曲线、解决样本获取困难问题以及在灾害区域更高容错率等优越性。
Claims (10)
1.一种基于标准曲线的大豆识别方法，其特征是包括以下步骤
步骤S1、卫星数据处理：
所述卫星数据处理步骤获取研究区内大豆生长期所有相关卫星数据，对数据进行云量控制和去云处理；
步骤S2、构建时间序列归一化物候指数特征：
使用处理好的卫星数据计算研究区内所有像元的植被指数获得所有像元的植被指数图像，对每一个像元使用该像元的时间序列方法处理该像元的植被指数图像构建该像元的时间序列植被指数特征；
步骤S3、构建样本点标准曲线；
步骤S4、样本点标准曲线参数信息提取与物候特征计算：
将样本点标准曲线中的大豆样本点关键生长期信息提取出来，得到样本点主要物候信息；对研究区范围内除样本点以外所有像元计算大豆样本点关键生长期信息，合成物候参数特征影像；
步骤S5、农作物品种识别与精度评估：
将样本信息输入到分类器中得到样本分类器，将样本点以外所有像元的物候参数特征影像输入到样本分类器中进行识别并输出识别结果，使用分类混淆矩阵计算识别精度。
2.根据权利要求1所述基于标准曲线的大豆识别方法，其特征是：所述步骤S2中的植被指数为归一化物候植被指数，归一化物候植被指数计算公式如下：
式中，NIR为近红外波段反射率，RED为红色波段反射率，SWIR为短波红外波段反射率。
3.根据权利要求1所述基于标准曲线的大豆识别方法，其特征是：所述步骤S2中的时间序列方法采用线性谐波模型，所述线性谐波模型计算公式如下：
其中，f(t)是第t时刻拟合的植被指数值；a是常数项；b是一阶项的系数；M是谐波组合的数量，在本申请中，M取2，即二项线性谐波模型；C和D是余弦函数和正弦函数的系数；ω是一年中天数的倒数(1/365)；t是DOY(一年中的某一天)。
4.根据权利要求1所述基于标准曲线的大豆识别方法，其特征是：所述步骤S3样本点标准曲线公式为：
其中，部分为谐波拟合公式，Temp为温度变量，hour为日照时长变量，rain为降雨变量，c0，c1，c2为气象权重因子。
5.根据权利要求1所述基于标准曲线的大豆识别方法，其特征是：大豆样本点关键生长期信息采用主要物候点参数信息，主要物候点通过下述公式来确定：
u＝(Vmin-Vmax)×p+Vmin
其中，u为参考阈值，用于确定大豆样本点物候开始时间和大豆样本点物候结束时间，Vmin和Vmax分别为大豆样本点时间序列中的最小值和最大值，p为振幅比例，p有三个值，第一个值对应大豆样本点生长开始节点、大豆样本点休眠开始节点，第二个值对应大豆样本点生长中期节点、大豆样本点衰老中期节点，第三个值对应大豆样本点成熟期节点、大豆样本点衰老开始节点。
6.根据权利要求1所述基于标准曲线的大豆识别方法，其特征是：步骤5中所述将样本点以外所有像元的物候参数特征影像输入到样本分类器中进行识别时，样本点以外所有像元的物候参数特征影像与样本分类器中生长条件最接近的样本点数据进行比对识别。
7.根据权利要求1至6任一权利要求所述基于标准曲线的大豆识别方法，其特征是：所述样本分类器采用随机森林分类器，分类混淆矩阵采用随机森林分类混淆矩阵，随机森林分类混淆矩阵定义如下：
公式(6)中，n表示类别数，mij表示实际属于i类的像元在分类结果图中被分到j类的像元数量，对角线元素的值即是各类别被正确分类的像元数，利用分类混淆矩阵总体精度(OA)、制图精度(PA)、用户精度(UA)和F1-Score，计算公式分别为：
其中，mi+为分类混淆矩阵中的行总和，m+i为分类混淆矩阵中的列总和。
8.根据权利要求1所述基于标准曲线的大豆识别方法，其特征是：所述构建样本点标准曲线直接采用上一年度该样本点标准曲线或者用历年该样本点平均值构成该样本点标准曲线。
9.根据权利要求1所述基于标准曲线的大豆识别方法，其特征是：所述构建样本点标准曲线包括：获取包括大豆在内的农作物典型代表样本点位置，使用S2步骤获得的样本点时间序列植被指数特征与样本点位置计算样本点时间序列曲线，使用样本点时间序列曲线并添加气象数据因子构建样本点标准曲线。
10.根据权利要求5所述基于标准曲线的大豆识别方法，其特征是：振幅比例p的三个值分别取0.15、0.5和0.90。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310197975.7A CN116824384A (zh) | 2023-03-03 | 2023-03-03 | 一种基于标准曲线的大豆识别方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310197975.7A CN116824384A (zh) | 2023-03-03 | 2023-03-03 | 一种基于标准曲线的大豆识别方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116824384A true CN116824384A (zh) | 2023-09-29 |
Family
ID=88115559
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310197975.7A Pending CN116824384A (zh) | 2023-03-03 | 2023-03-03 | 一种基于标准曲线的大豆识别方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116824384A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117111536A (zh) * | 2023-10-23 | 2023-11-24 | 上海永大菌业有限公司 | 一种菇棚环境远程控制系统及方法 |
| CN117151477A (zh) * | 2023-10-31 | 2023-12-01 | 中化现代农业有限公司 | 作物异常识别方法、装置、电子设备和存储介质 |
| CN117689959A (zh) * | 2024-01-30 | 2024-03-12 | 中交第二公路勘察设计研究院有限公司 | 一种融合植被生命周期特征的遥感分类方法 |
| CN117726902A (zh) * | 2023-12-19 | 2024-03-19 | 中国科学院空天信息创新研究院 | 一种基于Sentinel2时序曲线分布的大豆样本生成方法和系统 |
-
2023
- 2023-03-03 CN CN202310197975.7A patent/CN116824384A/zh active Pending
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117111536A (zh) * | 2023-10-23 | 2023-11-24 | 上海永大菌业有限公司 | 一种菇棚环境远程控制系统及方法 |
| CN117151477A (zh) * | 2023-10-31 | 2023-12-01 | 中化现代农业有限公司 | 作物异常识别方法、装置、电子设备和存储介质 |
| CN117151477B (zh) * | 2023-10-31 | 2024-01-23 | 中化现代农业有限公司 | 作物异常识别方法、装置、电子设备和存储介质 |
| CN117726902A (zh) * | 2023-12-19 | 2024-03-19 | 中国科学院空天信息创新研究院 | 一种基于Sentinel2时序曲线分布的大豆样本生成方法和系统 |
| CN117689959A (zh) * | 2024-01-30 | 2024-03-12 | 中交第二公路勘察设计研究院有限公司 | 一种融合植被生命周期特征的遥感分类方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110287944B (zh) | 基于深度学习的多光谱遥感影像的农作物虫害监测方法 | |
| CN116824384A (zh) | 一种基于标准曲线的大豆识别方法 | |
| CN111832518B (zh) | 基于时空融合的tsa遥感影像土地利用方法 | |
| Fraga et al. | Examining the relationship between the Enhanced Vegetation Index and grapevine phenology | |
| CN108458978B (zh) | 基于敏感波段和波段组合最优的树种多光谱遥感识别方法 | |
| Kordi et al. | Crop classification based on phenology information by using time series of optical and synthetic-aperture radar images | |
| CN113221765B (zh) | 一种基于数字相机影像有效像元的植被物候期提取方法 | |
| CN111598045A (zh) | 一种基于对象图谱和混合光谱的遥感耕地变化检测方法 | |
| CN114387516B (zh) | 一种针对复杂地形环境下中小田块的单季稻sar识别方法 | |
| Chen et al. | Investigating rice cropping practices and growing areas from MODIS data using empirical mode decomposition and support vector machines | |
| CN112434569B (zh) | 一种无人机热成像系统 | |
| CN112836725A (zh) | 基于时序遥感数据的弱监督lstm循环神经网络稻田识别方法 | |
| CN115810155B (zh) | 一种潮汐湿地的分类方法 | |
| CN114708490A (zh) | 水稻种植提取及复种指数监测方法、系统、终端及存储介质 | |
| CN112766036A (zh) | 一种倒伏玉米遥感提取方法与装置 | |
| CN117197668A (zh) | 基于深度学习的作物倒伏级别的预测方法及系统 | |
| CN115861629A (zh) | 一种高分耕地影像提取方法 | |
| Sabzchi-Dehkharghani et al. | Recognition of different yield potentials among rain-fed wheat fields before harvest using remote sensing | |
| CN114049564A (zh) | 基于高光谱遥感影像的松材线虫病等级预测模型构建方法 | |
| CN110929222A (zh) | 基于遥感植被冠层水分指数的灌溉农田识别方法 | |
| AU2021101780A4 (en) | Aboveground Biomass Estimation and Scale Conversion for Mean Regional Spectral Units | |
| CN113887493A (zh) | 一种基于id3算法的黑臭水体遥感影像识别方法 | |
| CN114140695B (zh) | 一种基于无人机多光谱遥感的茶树氮素诊断及品质指标测定的预测方法和系统 | |
| CN116310864A (zh) | 一种作物倒伏自动识别方法、系统、电子设备及介质 | |
| CN115424006A (zh) | 应用于作物表型参数反演的多源多层次数据融合方法 |
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 |