CN114118679A - 一种基于时序遥感数据的作物单产及长势评估方法 - Google Patents
一种基于时序遥感数据的作物单产及长势评估方法 Download PDFInfo
- Publication number
- CN114118679A CN114118679A CN202111200212.0A CN202111200212A CN114118679A CN 114118679 A CN114118679 A CN 114118679A CN 202111200212 A CN202111200212 A CN 202111200212A CN 114118679 A CN114118679 A CN 114118679A
- Authority
- CN
- China
- Prior art keywords
- yield
- time
- vegetation index
- crop
- remote sensing
- 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.)
- Granted
Links
- 230000012010 growth Effects 0.000 title claims abstract description 95
- 238000011156 evaluation Methods 0.000 title claims abstract description 47
- 230000001186 cumulative effect Effects 0.000 claims abstract description 57
- 238000000034 method Methods 0.000 claims abstract description 40
- 238000012544 monitoring process Methods 0.000 claims abstract description 21
- 238000000605 extraction Methods 0.000 claims abstract description 14
- 238000011160 research Methods 0.000 claims description 14
- 238000012545 processing Methods 0.000 claims description 12
- 238000009826 distribution Methods 0.000 claims description 7
- 238000004519 manufacturing process Methods 0.000 claims description 7
- 238000010899 nucleation Methods 0.000 claims description 7
- 238000007781 pre-processing Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000001914 filtration Methods 0.000 claims description 6
- 238000000611 regression analysis Methods 0.000 claims description 6
- 238000005520 cutting process Methods 0.000 claims description 4
- 238000003908 quality control method Methods 0.000 claims description 4
- 238000005516 engineering process Methods 0.000 abstract description 5
- 244000068988 Glycine max Species 0.000 description 31
- 235000010469 Glycine max Nutrition 0.000 description 31
- 230000000694 effects Effects 0.000 description 10
- 238000012360 testing method Methods 0.000 description 10
- 238000004458 analytical method Methods 0.000 description 9
- 238000011161 development Methods 0.000 description 9
- 230000018109 developmental process Effects 0.000 description 9
- 238000012417 linear regression Methods 0.000 description 6
- 238000009331 sowing Methods 0.000 description 4
- 238000009825 accumulation Methods 0.000 description 3
- 238000011835 investigation Methods 0.000 description 3
- 239000002028 Biomass Substances 0.000 description 2
- 238000010220 Pearson correlation analysis Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 235000013305 food Nutrition 0.000 description 2
- 230000009666 qualitative growth Effects 0.000 description 2
- 238000002310 reflectometry Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 241000196324 Embryophyta Species 0.000 description 1
- 244000061176 Nicotiana tabacum Species 0.000 description 1
- 235000002637 Nicotiana tabacum Nutrition 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 241000282376 Panthera tigris Species 0.000 description 1
- 240000008042 Zea mays Species 0.000 description 1
- 235000005824 Zea mays ssp. parviglumis Nutrition 0.000 description 1
- 235000002017 Zea mays subsp mays Nutrition 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 235000005822 corn Nutrition 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 239000010931 gold Substances 0.000 description 1
- 229910052737 gold Inorganic materials 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 230000029553 photosynthesis Effects 0.000 description 1
- 238000010672 photosynthesis Methods 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 230000009663 quantitative growth Effects 0.000 description 1
- 230000029058 respiratory gaseous exchange Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000012358 sourcing Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000008961 swelling Effects 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 230000017260 vegetative to reproductive phase transition of meristem Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
- G06Q10/063—Operations research, analysis or management
- G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
- G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
Landscapes
- Business, Economics & Management (AREA)
- Human Resources & Organizations (AREA)
- Engineering & Computer Science (AREA)
- Economics (AREA)
- Strategic Management (AREA)
- Physics & Mathematics (AREA)
- General Business, Economics & Management (AREA)
- Entrepreneurship & Innovation (AREA)
- Marketing (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Tourism & Hospitality (AREA)
- Development Economics (AREA)
- Educational Administration (AREA)
- Quality & Reliability (AREA)
- Operations Research (AREA)
- Game Theory and Decision Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Agronomy & Crop Science (AREA)
- Animal Husbandry (AREA)
- Marine Sciences & Fisheries (AREA)
- Mining & Mineral Resources (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本申请涉及遥感技术的应用，具体涉及一种基于时序遥感数据的作物单产及长势评估方法，其中，作物单产评估方法包括如下步骤：获取MODIS遥感数据；对MODIS遥感数据进行归一化植被指数提取，以构建NDVI时序数据集；生产累积植被指数变化率时序产品集；获取作物单产统计数据集；构建作物生长周期内不同时间节点所对应的单产预测模型；进行目标作物估产。本申请的基于时序遥感数据的作物单产及长势评估方法中，以累积植被指数变化率作为作物估产与长势监测评估的特征参量，提高了作物估产和长势监测评估的有效性、合理性，且更易于实现。
Description
技术领域
本申请涉及遥感技术的应用，具体涉及一种基于时序遥感数据的作物单产及长势评估方法。
背景技术
在全球自然环境和国际形势复杂多变的情况下，及时、准确地评估作物长势状况，预测作物产量对国家粮食安全管理至关重要。传统的农作物估产、长势评估是采用人工区域调查方法，速度慢、工作量大、成本高。
近年，随着遥感技术的发展，其在农作物估产、长势评估方面的运用已经步入了新的发展阶段。当前作物长势遥感技术总体可以分为两类，即定性长势监测方法和定量长势监测方法。其中，定性长势监测方法是利用植被光谱响应敏感波段构建能够反映作物生长状况的遥感指数例如植被指数、叶面积指数等进行监测。
但是，基于植被指数的作物估产、长势评估方法中，由于缺乏地面观测数据，导致产量过程性预测不够稳定、长势评价标准偏于主观等问题。
发明内容
为了解决现有技术中存在的至少一个技术问题，本申请提供了一种基于时序遥感数据的作物单产及长势评估方法。
第一方面，本申请公开了一种基于时序遥感数据的作物单产评估方法，包括如下步骤：
步骤1、构建作物单产预测模型
步骤101、获取研究区域内多个年份范围内的时序遥感数据(例如MODIS遥感数据)；
步骤102、对时序遥感数据进行归一化植被指数提取，以构建NDVI时序数据集；
步骤103、提取植被指数特征值及其对应的时间节点：结合历年时序遥感数据提取作物生长过程特征曲线，对播种期、成熟期最低谷值和生长期最高峰值时点进行位置识别，对应时段分别记为
并提取相应NDVImin_s、NDVImin_m、NDVImax_g值；并基于已经提取的最低谷值和最高峰值，分别用算术平均法得到生长前期中值NDVImid_gf，其对应时段记为
和生长后期中值NDVImid_gb，其对应时段记为
步骤104、计算不同时间点的累积植被指数变化率
先计算各个时间节点的植被指数变化率；然后依次计算
到
至
之间任一时间节点的累积植被指数变化率；
其中，植被指数变化率和累积植被指数变化率通过如下公式获得：
式中，CRVIi为i时段的植被指数变化率，NDVIi和NDVIi-1分别为i时段和i-1时段的植被指数，i表示时序产品对应的时间节点，i-1为i紧邻的前一个时间节点；ACRVIi为截止i时间节点的累积植被指数变化率，t0为播种期最低值出现时间，T为特定时间节点；
步骤105、获取所述多个年份所对应的作物单产统计数据，以形成作物单产统计数据集；
步骤106、将所述累积植被指数变化率时序产品集和对应的作物单产统计数据集进行回归分析，得出不同特定时间节点上针对研究区域作物单产预测模型；
步骤2、单产评估
获取待评估年份的遥感数据，根据步骤103的方法提取植被指数特征值及其对应的估产时间节点，根据步骤104的方法计算待评估年份不同时间节点的植被指数变化率及累积植被指数变化率，将其代入步骤106中得出的对应时间节点的作物单产预测模型中，经计算输出目标作物单产评估结果。
根据本申请的至少一个实施方式，如果估产时没有完整的作物生长季相关遥感数据，采用历史多年对应时间节点的有效均值特征影像作为未结束年份的准同期数据，随着作物生长期发展，进行数据的迭代更新。
需要说明的是，根据本发明的思路，采用的遥感数据源除MODIS外，还可以是其他卫星遥感数据；另外，除归一化植被指数(NDVI)外，可以采用其它反映作物生长状况的特征参量数据，例如叶面积指数(LAI)、增强植被指数(EVI)。
根据本申请的至少一个实施方式，在所述步骤102中，对MODIS遥感数据进行归一化植被指数提取之前，还包括：
对MODIS遥感数据进行预处理，所述预处理包括拼接处理、剪裁处理以及质量控制处理。
根据本申请的至少一个实施方式，在所述步骤102中，对预处理后的MODIS遥感数据进行归一化植被指数提取之前，还包括：
利用中高分辨率遥感数据进行监测区域目标作物识别，得到目标作物分布数据；
利用目标作物分布数据对预处理后的MODIS遥感数据进行掩膜处理。
根据本申请的至少一个实施方式，在所述步骤102中，对掩膜处理的MODIS遥感数据进行归一化植被指数提取之后，还包括：
进行滤波处理，本发明的实施方案中优选采用Savitzky-Golay滤波处理。
根据本申请的至少一个实施方式，数据时间步长为8d，所述多个年份指8年以上。
根据本申请的至少一个实施方式，在所述步骤2中进行目标作物估产时，目标作物的估产时间节点选自
其中，m取1-9之间自然数。
第二方面，本申请还公开了一种基于时序遥感数据的作物长势评估方法，包括如下步骤：
步骤301、采用上述第一方面中任一项所述的基于时序遥感数据的作物单产评估方法，得到目标作物的单产评估结果；
步骤302、根据得到的单产评估结果，结合单产等级划分标准进行单产等级评估，其中，单产等级包括高产、中等水平、低产；
步骤303、根据单产等级划分标准和产量预测模型测算各特定时间节点长势植被指数特征参量的等级标准后，评估待评估区域目标作物的长势。其中，单产等级的高产、中等水平、低产分别对应长势等级的好、中等、差。
根据本申请的至少一个实施方式，在所述步骤302中，单产等级以多年单产均值μ和标准差σ确定，将(μ-0.5σ)≤评估产量≤(μ+0.5σ)定级为中等水平，评估产量大于(μ+0.5σ)定级为高产，评估产量小于(μ-0.5σ)定级为低产。
本申请的创新及有益技术效果：
本发明提出了一个新的作物估产与长势监测评估的特征参量：累积植被指数变化率。应用试验显示，该参量较累积NDVI、累积植被变化速率在作物估产和长势评估方面更具优势；另外，确立了基于产量的相对客观且应用便捷的长势遥感监测评估模式。
具体创新点如下：
1)本申请的基于时序遥感数据的作物单产及长势评估方法中，以累积植被指数变化率作为作物估产与长势监测评估的特征参量，一方面利于反映作物在特定生长时段内生物量累积变化状况，另一方面可以降低作物时空异质性影响，增强不同时空范围的可比性；
2)本申请提出了累积植被指数变化率特征参量的提取方法，提出了合理的累积时段范围并明确了适宜的产量评估时间区间(比如本发明在示例性实施例中发现，合理的累积时段范围是自
到按照数据时间步长依次划分的包含
至
在内的任一时间节点；适宜的产量评估时间区间是生长前期中值对应时间节点+8d至生长前期中值对应时间节点+72d的区间)，从而提高了累积植被指数变化率特征参量用于作物估产和长势监测评估的有效性；
3)本申请中提出了基于预测产量的长势量化评估方法，基于该模式的长势监测评估更加客观合理，且易于实现。
附图说明
图1是本申请基于时序遥感数据的作物单产评估方法以及长势评估方法的总体流程图；
图2是某一具体作物(大豆)生长过程植被指数特征曲线示意图；
图3是本申请一具体研究试验区域(巴西巴拉那州西部的托莱多市)的空间位置示意图；
图4是本申请一具体研究试验中，特定时间节点累积植被指数变化率与大豆产量散点图；
图5是本申请一具体研究试验中，反映多时段累积植被指数变化率与单产相关性的示意图；
图6是本申请一具体研究试验中，反映多时段累积NDVI与单产相关性的示意图；
图7是本申请一具体研究试验中，反映多时段累积植被变化速率与单产相关性的示意图；
图8是本申请一具体研究试验中，多时段累积植被指数特征参量相对误差对比图；
图9是本申请一具体研究试验中，特定年份大豆开花结荚期长势监测图。
具体实施方式
为使本申请实施的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行更加详细的描述。
1.研究方法
1.1遥感数据处理
1、基于GIS软件或相关平台，对MODIS遥感数据进行预处理，包括拼接、剪裁、质量控制等；
2、利用中高分辨率遥感数据，结合物候信息、解译标志数据等，进行监测区域作物识别，得到作物分布数据，本发明中以大豆作物为例；
3、利用作物分布数据对MODIS时序数据进行掩膜处理，得到目标作物区域时序数据集；
4、生产归一化植被指数(NDVI)产品，进行Savitzky-Golay滤波处理，构建NDVI 时序数据集；
NDVI是由红光和近红外两个波段得到的，由于该两个波段是植物光谱中反映光合作用和呼吸作用的最重要波段，该指数能够反映植被生长差异特征，且指数计算简单，易于实现，被广泛用于作物长势监测中。其形式如下：
式中，Rnir为近红外波段反射率，Rred为红光波段反射率。
1.2累积植被指数变化率(Accumulated Changing Rate of VI，ACRVI)
以变化指数表征的植被生长过程通常包括正增长阶段和负增长阶段(即下降阶段)，利用前后植被指数差值可以表征其增长趋势，为突出表征某一时点植被变化程度和趋势，本研究提出植被指数变化率，形式见公式(2)；为表征植被生长过程变化程度，提出累积植被指数变化率，形式见公式(3)：
式中，CRVIi为i时段的植被指数变化率，NDVIi和NDVIi-1分别为i时段和i-1时段的植被指数，i表示时序产品对应的时间节点，i-1为i紧邻的前一个时间节点；
ACRVIi为截止i时间节点的累积植被指数变化率，t0为初始时间，本研究以播种期最低值出现时间为t0，T为特定时间节点。
1.3累积植被指数变化率特征值提取
结合作物生长过程曲线，提取5个植被指数特征值，包括播种期最低谷值(记为NDVImin_s)、生长期最高峰值(记为NDVImax_g)、生长前期最低-最高的中值(简称生长前期中值，记为NDVImid_gf)、生长后期最高-最低的中值(简称生长后期中值，记为NDVImid_gb) 和成熟期最低谷值(记为NDVImin_m)。
首先，根据研究区大豆种植习惯，确定监测起止时段(以月记)；
其次，采用二次差分法对播种期、成熟期最低谷值和生长期最高峰值时点进行位置识别，对应时段分别记为
再提取相应NDVI值，基于已经提取的最低谷值和最高峰值，分别用算术平均法得到生长前期中值(其对应时段记为
)和生长后期中值(其对应时段记为
)；
再次，计算
到
的CRVI值；
最后，以
为累计植被指数变化率统计起点，以
为统计终点，依次计算
到上述统计区间内按数据时间步长划分的不同时间节点的累积植被指数变化率ACRVI，从而得到不同年份截止不同时间点的累积植被指数变化率序列，数据时间步长为8d。
需要指出的是，通常作物产量预测和长势评估是在作物生长季中进行的，即可能没有完整的作物生长季相关遥感数据，因此，本研究参考Hoolst R V等[1]的方法得到多年不同时点的有效均值特征影像并以此作为未结束年份的准同期数据，随着作物生长期发展，进行数据的迭代更新。
1.4其他特征参量提取
为客观评估基于累积植被指数变化率的估产和长势评估效果，本研究同步计算得到较为常用的两个特征参量，即累积NDVI和累积植被指数变化速率。计算公式如下：
式中，ANDVIi为截止i时间节点的累积植被指数，t0为初始时间，本研究以播种期最低值出现时间为t0；SRVIi为i时段时间节点的植被指数变化速率，t为i时段时间节点和i-1时段时间节点的时间间隔；ASRVIi为截止i时间节点的累积植被指数变化速率。
1.5精度分析方法
为了评价模型的精度,采用常用的模型检验指标对其进行分析评价。通过计算模型模拟值和实测值之间的相对误差、决定系数对它们之间的符合度进行分析。计算公式如下：
式中，Re为相对误差，yi和y分别为预测值和实际值；R2为决定系数，n为样本数。R2取值范围为(0，1)，其值越接近1，表示模拟值和实测值的符合度越高，相反，越接近 0表明符合度越低。
示例性应用：
1)研究区域
本研究区为位于巴西巴拉那州西部的托莱多市，其经度范围是-54°3′～-53°32′，纬度范围是-24°27′～-24°56′，海拔在254～690m之间，平均海拔为503m，全县面积为1196.2km2(图3)。托莱多气温通常在9.4℃～30℃之间。肥沃的土壤和相对平坦的地势使该市成为农业发达区。它是该州重要粮食生产基地之一，其农业GDP在巴拉那州和南部地区排名首位，农业增加值在全国位列第三[2]。主要农作物类型为大豆、玉米、烟叶、咖啡等。其中，大豆种植生长季为9月～次年3月。
2)数据来源与处理
2.1)统计调查数据
托莱多市大豆年度单产统计数据和苗情进度与长势调查数据来源于巴西巴拉那州农村经济部(DERAL)。其中，使用的单产统计数据年限为2010/11～2020/21，苗情进度与长势调查数据年限为2019/20～2020/21。
2.2)遥感及基础地理数据
本研究使用中分辨率成像光谱仪(Moderate Resolution ImagingSpectroradiometer， MODIS)的MOD09A1(地表反射率)8d合成产品，空间分辨率为500m，数据时间为2010 年9月～2021年3月。对其进行拼接、剪裁、质量控制等预处理。由于使用农作物掩膜之后的预报结果准确率最高[3]，因此，利用Landsat8数据基于机器学习方法对研究区进行大豆种植区域提取，并以此为掩膜数据，完成掩膜处理。再进行归一化植被指数(NDVI) 提取，经过Savitzky-Golay[4]滤波处理后，生成大豆产区NDVI时序产品。
本研究中NDVI时序产品的时间表达以8d合成的截止日期为准。
2.3)提取植被指数特征值及其对应的特定时间节点
结合2010年9月～2021年3月时序遥感数据提取作物生长过程特征曲线，对各年的播种期、成熟期最低谷值和生长期最高峰值时点进行位置识别，对应时段分别记为
并提取相应NDVImin_s、NDVImin_m、NDVImax_gNDVI值；并基于已经提取的最低谷值和最高峰值，分别用算术平均法得到生长前期中值NDVImid_gf，其对应时段记为
和生长后期中值NDVImid_gb，其对应时段记为
2.4)计算不同时间节点的累积植被指数变化率
累积植被指数变化率通过如下公式获得：
式中，CRVIi为i时段的植被指数变化率，NDVIi和NDVIi-1分别为i时段和i-1时段的植被指数，i表示时序产品对应的时间节点，i-1为i紧邻的前一个时间节点；ACRVIi为截止i时间节点的累积植被指数变化率，t0为播种期最低值出现时间，T为特定时间节点。
步骤106、将所述累积植被指数变化率时序产品集和对应的作物单产统计数据集进行回归分析，得出不同特定时间点上针对研究区域作物单产预测模型。
3)研究结果与分析
3.1)累积植被指数特征参量与产量关系分析
3.11)
和
的累积植被指数特征参量与单产关系分析
对2010/11-2018/19年的
的累积植被指数变化率与大豆单产进行回归分析，如图4(a)所示，累积植被指数变化率与大豆单产呈现正向线性关系。
利用Pearson相关分析计算相关系数为0.88，决定系数为0.77，并通过显著性水平为0.01的极显著性检验；利用线性回归模型，对2019/20和2020/21两个年份单产进行预测，得到相对误差分别为-2.51％和2.14％，平均相对误差为-0.18％。
累积NDVI与大豆单产的相关系数为0.79，决定系数为0.62，通过显著性水平为0.05 的显著性检验。利用线性回归模型，得到2019/20和2020/21两个年份预测单产结果，其相对误差分别为13.97％和27.20％，平均相对误差为20.59％。
累积植被指数变化速率与大豆单产的相关系数为0.89，决定系数为0.8，通过显著性水平为0.01的极显著性检验。利用线性回归模型，得到2019/20和2020/21两个年份预测单产结果，其相对误差分别为10.74％和0.78％，平均相对误差为-5.76％。
经比较认为，在该时段，基于累积植被指数变化率进行产量预测效果最佳，累积植被指数变化速率预测效果次之，累积NDVI预测效果不佳。
对2010/11-2018/19年的
的累积植被指数变化率与大豆单产进行回归分析，如图4(b)所示，累积植被指数变化率与大豆单产呈现负向线性关系。
利用Pearson相关分析计算相关系数为0.54，决定系数为0.288，未通过显著性水平为0.05的检验。利用线性回归模型，对2019/20和2020/21两个年份单产进行预测，得到相对误差分别为-21.48％和-4.48％，平均相对误差为-12.98％。
累积NDVI与大豆单产的相关系数为0.79，决定系数为0.63，未通过显著性水平为0.05的显著性检验。利用线性回归模型，得到2019/20和2020/21两个年份预测单产结果，其相对误差分别为-8.82％和4.65％，平均相对误差为-2.09％。
累积植被指数变化速率与大豆单产的相关系数为0.77，决定系数为0.59，通过显著性水平为0.05的显著性检验。利用线性回归模型，得到2019/20和2020/21两个年份预测单产结果，其相对误差分别为-18.62％和-0.91％，平均相对误差为-9.77％。
该阶段，仅累积植被指数变化速率通过显著性检验，但平均相对误差偏高，虽然累积NDVI平均相对误差较低，好于另外两个特征参量拟合效果，但未通过显著性检验，因此认为，在该阶段三种特征参量均不适宜进行产量预测。
3.12)多时段累积植被变化特征参量与单产关系分析
以生长前期中值对应时间节点
为起始时间，分别以其顺次加8d，直至加 80d为计算累积植被变化特征参量的终止时间，共计11个监测时间模式，对 2010/11-2018/19共计9年的大豆单产与11个监测时间模式的三个累积植被指数变化特征参量进行相关性分析，结果如图5-图7所示(其中，**为极显著相关，*为显著相关)。
累积植被指数变化率和累积植被指数变化速率与单产相关系数时序尺度上分布较为相似，呈抛物线形状，在前期和后期均较低，相关关系表现为不显著，随着生长阶段发展，显著性水平表现为显著或极显著，相关关系在生长前期中值对应时间节点+40d时达到峰值，分别为0.96和0.97。累积NDVI除生长前期中值对应时间节点+24d时出现低值外，总体呈现线性递增趋势，即随生长阶段发展，相关系数逐渐上升，在生长后期达到高值，为0.89，且表现为极显著性相关。
生长前期中值主要为大豆苗期-分枝期，刚进入该阶段时，累积植被指数变化率、累积NDVI及累积植被指数变化速率三个特征参量与单产相关系数均较低，且相关关系未达到显著水平。随着大豆生长，累积植被指数变化率和累积植被指数变化速率与单产呈现显著或极显著相关关系，前者相关系数总体高于后者，可认为累积植被指数变化率在生长前期中值对应时间节点+8d至生长前期中值对应时间节点+72d期间，即TNDVImid_gf+8d～TNDVImid_gf+72d期间，对大豆单产均有比较好的反映，可以用于大豆估产，是三个特征参量中的优选指标。另外，累积NDVI与大豆单产相关系数在生长前期中值对应时间节点 +72d及以后表现好于另外两个特征参量，可以作为全生育期产量估测参量。
进一步，对2010/11-2018/19年的
的累积植被指数特征参量与大豆单产进行回归分析，建立不同时期的线性拟合模型，并基于模型对2019/20和2020/21两个年份单产进行预测，得到基于三个累积植被指数特征参量的单产预测相对误差结果，如图8所示：
基于累积NDVI的预测结果总体偏高，且相对误差较高，其中，2019/20预测相对误差在-13.53％～36.63％之间，2020/21预测相对误差在-0.07％～43.71％之间。
基于累积植被指数变化速率的预测结果总体偏低，在大豆生长中后期相对误差两年份差异较大，其中，2019/20预测相对误差在-17.76％～4.33％之间，2020/21预测相对误差在 -7.35％～5.3％之间。
基于累积植被指数变化率的预测结果在前期偏高、后期偏低，其中，2019/20预测相对误差在-13.52％～14.65％之间，2020/21预测相对误差在-7.43％～6.98％之间。
对2019/20和2020/21两年相对误差进行平均，得到三个特征参量的相对误差均值时序曲线，其中，基于累积NDVI的相对误差均值在-6.8％～40.2％之间，基于累积植被指数变化率的相对误差均值在-6.8％～6.6％之间，基于累积植被指数变化速率的相对误差均值在-11.0％～-0.6％之间。
进一步，对2019/20和2020/21两个年份基于三个累积植被指数特征参量的单产预测相对误差进行分时段平均，一个时段为
一个时段为
结果如表1所示。
表1多时段累积植被指数特征参量不同时段平均相对误差对比
可以看出：基于累积NDVI的预测相对误差总体最高，2019/20年度多时段平均相对误差为23.28％，2020/21年度平均相对误差为32.74％。
平均相对误差最小的为累积植被指数变化率，其中，2019/20年度多时段平均相对误差为0.48％，2020/21年度平均相对误差为0.71％。基于累积植被指数变化速率的预测效果次之，2019/20的为-6.96％，2020/21的为-1.6％。
由于
多为大豆的开花结荚-鼓粒期，是大豆产量形成关键阶段，因此对该阶段预测平均相对误差进行比较分析。
基于累积NDVI的预测相对误差总体最高，2019/20年度的平均相对误差为33.22％， 2020/21年度平均相对误差为40.65％。平均相对误差总体最小的为累积植被指数变化率，其中，2019/20年度的平均相对误差为0.11％，2020/21年度的平均相对误差为3.69％。基于累积植被指数变化速率的预测效果次之，2019/20的为-6.42％，2020/21的为1.17％。
3.2)基于累积植被指数变化特征参量的典型年长势分析
基于对累积植被指数变化特征参量与单产关系分析结果，认为在特定预测精度范围内，由于该累积植被指数变化特征参量可以反映作物未来产量水平和趋势，因此可以利用相关参量进行不同时段的长势监测评估。
选取2010/11～2020/21间典型年，包括2016/17(最高产年)、2019/20(上一年)、2020/21(当年)三个年份的最低谷值-最高峰值时段进行应用分析。
各特定时间节点长势等级需基于单产等级和单产预测模型进行推算，其中，单产等级以多年单产均值μ和标准差σ确定，本研究中，产量在[μ-0.5σ，μ+0.5σ]区间定级为中等水平(其对应的长势等级为中等)，大于(μ+0.5σ)定级为好(其对应的长势等级为好)，小于(μ-0.5σ)定级为差(其对应的长势等级为差)。
基于该标准得到基于累积植被指数变化特征参量的三个年份最高峰值期的长势评价结果，如图9所示。
利用DERAL2019/20和2020/21两年同期调查统计数据进行对比分析和相对误差分析，如下表2所示。
表2长势等级与精度统计
根据统计数据，2019/20年为3846kg/ha，2020/21年大豆单产为3350kg/ha。按年份，基于三个累积植被指数变化特征参量的长势监测图(图9)总体上都能体现现出2020/21年长势不及前两年份，根据表2看出，该年度苗情为好的比例低于2019/20。这是由于该年度大豆种植和出苗期出现了较重干旱，导致部分地区播种偏晚、苗情偏弱，影响了生物量积累，长势表现偏差。
为评估基于三个累积植被指数变化特征参量的长势评价效果，结合统计调查数据，进行了相对误差测算。
由表2可知，基于累积植被指数变化率的2019/20年度的相对误差分别为0.5％、1.3％和-4.4％；2020/21年度的相对误差分别为-1.9％、-0.2％和23.1％。两年份中，好和中等等级的平均相对误差分别为-0.7％和0.5％，差等级偏低，为9.4％，总体上，相对误差低于10％，尤其是好和中等的相对误差非常低，在-1.9％～1.3％之间。
基于累积NDVI的2019/20年度的相对误差分别为7.6％、4.3％和-53.4％；2020/21年度的相对误差分别为20.9％、-47.8％和-52.3％。两年份中，不同等级的相对误差差异较大，尤其是中等和差的等级，这是由于基于累积NDVI的估产在常年估计时有偏高、高产年估计时有偏低的趋势，因此导致在长势评估中也出现了较大的偏差。
基于累积植被指数变化速率的2019/20年度的相对误差分别为-3.4％、29.1％和3.0％； 2020/21年度的相对误差分别为-3.6％、9.8％和2.6％。两年份中，好和差等级的平均相对误差分别为-3.5％和2.8％，中等等级偏低，为19.5％。其评估效果介于前两特征参量之间。
综上，总体看，累积植被指数变化率是三个特征参量中最具可行性和有效性的长势评估参量。
以上所述，仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以所述权利要求的保护范围为准。
参考文献：
[1]Hoolst R V,Eerens H,Haesen D,et al.FAO's AVHRR-based AgriculturalStress Index System(ASIS)for global drought monitoring[J].InternationalJournal of Remote Sensing,2016,37(1-2):418-439.
[2]Manfrin,J.,Ceraso,A.F.F.,&
Jr,A.C.Evaluation of OdoursPerception by Residents of the Municipality of Toledo/ParanáState/Brazil[J].Journal of Scientific Research and Reports,2018,21(2),1-9.https://doi.org/10.9734/JSRR/2018/44890.
[3]李新星.基于时间序列植被指数的水稻发育期提取和估产[D].浙江大学,2017.
[4]边金虎,李爱农,宋孟强,et al.MODIS植被指数时间序列Savitzky-Golay滤波算法重构[J]. 遥感学报,2010,14(04):725-741.
Claims (9)
1.一种基于时序遥感数据的作物单产评估方法，包括如下步骤：
步骤1、构建作物单产预测模型
步骤101、获取研究区域内多个年份范围内的时序遥感数据；
步骤102、对时序遥感数据进行归一化植被指数(Normalized Difference VegetationIndex，NDVI)提取，以构建NDVI时序数据集；
步骤103、提取植被指数特征值及其对应的时间节点：结合历年时序遥感数据提取作物生长过程特征曲线，对播种期、成熟期最低谷值和生长期最高峰值时点进行位置识别，对应时段分别记为
并提取相应NDVImin_s、NDVImin_m、NDVImax_g值；基于已经提取的最低谷值和最高峰值，分别用算术平均法得到生长前期中值NDVImid_gf，其对应时段记为
和生长后期中值NDVImid_gb，其对应时段记为
步骤104、计算不同时间节点的累积植被指数变化率；
先计算各个时间节点的植被指数变化率；然后依次计算
到
至
之间任一时间节点的累积植被指数变化率；
其中，植被指数变化率和累积植被指数变化率通过如下公式获得：
式中，CRVIi为i时段的植被指数变化率，NDVIi和NDVIi-1分别为i时段和i-1时段的植被指数，i表示时序产品对应的时间节点，i-1为i紧邻的前一个时间节点；ACRVIi为截止i时间节点的累积植被指数变化率，t0为播种期最低值出现时间，T为特定时间节点；
步骤105、获取所述多个年份所对应的作物单产统计数据，以形成作物单产统计数据集；
步骤106、将所述累积植被指数变化率时序产品集和对应的作物单产统计数据集进行回归分析，得出不同时间节点上针对研究区域作物单产预测模型；
步骤2、单产评估
获取待评估年份的遥感数据，根据步骤103的方法提取植被指数特征值及其对应的估产时间节点，根据步骤104的方法计算待评估年份不同时间节点的植被指数变化率及累积植被指数变化率，将其代入步骤106中得出的对应时间节点的作物单产预测模型中，经计算输出目标作物单产评估结果。
2.根据权利要求1所述的方法，其特征在于，如果单产评估时缺乏完整的作物生长季相关遥感数据，则采用历史多年对应的时间节点的有效均值特征影像作为未结束年份的准同期数据，随着作物生长期发展，进行数据的迭代更新。
3.根据权利要求1所述的方法，其特征在于，在所述步骤102中，对MODIS遥感数据进行归一化植被指数提取之前，还包括：
对MODIS遥感数据进行预处理，所述预处理包括拼接处理、剪裁处理以及质量控制处理。
4.根据权利要求1所述的方法，其特征在于，在所述步骤102中，对预处理后的MODIS遥感数据进行归一化植被指数提取之前，还包括：
利用中高分辨率遥感数据进行监测区域目标作物识别，得到目标作物分布数据；
利用目标作物分布数据对预处理后的MODIS遥感数据进行掩膜处理。
5.根据权利要求4所述的方法，其特征在于，在所述步骤102中，对掩膜处理的MODIS遥感数据进行归一化植被指数提取之后，还包括：滤波处理。
6.根据权利要求1所述的方法，其特征在于，数据时间步长为8d，所述多个年份指8年以上。
7.根据权利要求6所述的方法，其特征在于，在所述步骤2中进行目标作物估产时，目标作物的估产时间节点选自
其中，m取1-9之间的自然数。
8.一种基于时序遥感数据的作物长势评估方法，包括如下步骤：
步骤301、采用权利要求1-7任一项所述的基于时序遥感数据的作物单产评估方法，得到待评估区域目标作物的单产评估结果；
步骤302、根据得到的单产评估结果，结合单产等级划分标准进行单产等级评估，其中，单产等级划分为高产、中等水平或低产；
步骤303、根据单产等级划分标准和产量预测模型测算各特定时间节点长势植被指数特征参量的等级标准后，评估待评估区域目标作物的长势，其中，单产等级的高产、中等水平、低产分别对应长势等级的好、中等、差。
9.根据权利要求8所述的作物长势评估方法，其特征在于，在所述步骤302中，单产等级以多年单产均值μ和标准差σ确定，(μ-0.5σ)≤评估产量≤(μ+0.5σ)为中等水平，评估产量大于(μ+0.5σ)为高产，评估产量小于(μ-0.5σ)为低产。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111200212.0A CN114118679B (zh) | 2021-10-14 | 2021-10-14 | 一种基于时序遥感数据的作物单产及长势评估方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111200212.0A CN114118679B (zh) | 2021-10-14 | 2021-10-14 | 一种基于时序遥感数据的作物单产及长势评估方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114118679A true CN114118679A (zh) | 2022-03-01 |
| CN114118679B CN114118679B (zh) | 2022-09-16 |
Family
ID=80375655
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111200212.0A Active CN114118679B (zh) | 2021-10-14 | 2021-10-14 | 一种基于时序遥感数据的作物单产及长势评估方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114118679B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115511224A (zh) * | 2022-11-11 | 2022-12-23 | 中国农业科学院农业资源与农业区划研究所 | 天地一体化的作物长势智能监测方法、装置及电子设备 |
| CN116976516A (zh) * | 2023-08-01 | 2023-10-31 | 中国科学院空天信息创新研究院 | 一种干旱地区作物单产早期预测方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106772429A (zh) * | 2016-12-24 | 2017-05-31 | 福州大学 | 基于生长盛期nmdi增减比值指数的玉米自动制图方法 |
| CN108362267A (zh) * | 2018-01-09 | 2018-08-03 | 浙江大学 | 基于卫星数据的湿渍害胁迫下油菜产量损失遥感定量评估方法 |
| CN110309985A (zh) * | 2019-07-10 | 2019-10-08 | 北京师范大学 | 一种农作物产量预测方法及*** |
| WO2020132092A1 (en) * | 2018-12-19 | 2020-06-25 | The Board Of Trustees Of The University Of Illinois | Apparatus and method for crop yield prediction |
| WO2021007352A1 (en) * | 2019-07-08 | 2021-01-14 | Indigo Ag, Inc. | Crop yield forecasting models |
-
2021
- 2021-10-14 CN CN202111200212.0A patent/CN114118679B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106772429A (zh) * | 2016-12-24 | 2017-05-31 | 福州大学 | 基于生长盛期nmdi增减比值指数的玉米自动制图方法 |
| CN108362267A (zh) * | 2018-01-09 | 2018-08-03 | 浙江大学 | 基于卫星数据的湿渍害胁迫下油菜产量损失遥感定量评估方法 |
| WO2020132092A1 (en) * | 2018-12-19 | 2020-06-25 | The Board Of Trustees Of The University Of Illinois | Apparatus and method for crop yield prediction |
| WO2021007352A1 (en) * | 2019-07-08 | 2021-01-14 | Indigo Ag, Inc. | Crop yield forecasting models |
| CN110309985A (zh) * | 2019-07-10 | 2019-10-08 | 北京师范大学 | 一种农作物产量预测方法及*** |
Non-Patent Citations (2)
| Title |
|---|
| ZHONGLIN JI 等: "Prediction of crop yield using phenological information extracted from remote sensing vegetation index", 《SENSORS》 * |
| 陈鹏飞 等: "基于环境减灾卫星时序归一化植被指数的冬小麦产量估测", 《农业工程学报》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115511224A (zh) * | 2022-11-11 | 2022-12-23 | 中国农业科学院农业资源与农业区划研究所 | 天地一体化的作物长势智能监测方法、装置及电子设备 |
| CN116976516A (zh) * | 2023-08-01 | 2023-10-31 | 中国科学院空天信息创新研究院 | 一种干旱地区作物单产早期预测方法 |
| CN116976516B (zh) * | 2023-08-01 | 2024-03-15 | 中国科学院空天信息创新研究院 | 一种干旱地区作物单产早期预测方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114118679B (zh) | 2022-09-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Fischer et al. | Global agro-ecological zones (gaez v4)-model documentation | |
| CN114118679B (zh) | 一种基于时序遥感数据的作物单产及长势评估方法 | |
| CN109711102B (zh) | 一种作物灾害损失快速评估方法 | |
| CN111598019A (zh) | 基于多源遥感数据的作物类型与种植模式识别方法 | |
| CN109800921B (zh) | 基于遥感物候同化和粒子群优化的区域冬小麦估产方法 | |
| Okada et al. | Projecting climate change impacts both on rice quality and yield in Japan | |
| Zhao et al. | Normalized NDVI valley area index (NNVAI)-based framework for quantitative and timely monitoring of winter wheat frost damage on the Huang-Huai-Hai Plain, China | |
| Joshi et al. | Winter wheat yield prediction in the conterminous United States using solar-induced chlorophyll fluorescence data and XGBoost and random forest algorithm | |
| Rao et al. | Remote sensing: A technology for assessment of sugarcane crop acreage and yield | |
| CN107437262B (zh) | 作物种植面积预警方法和*** | |
| CN111985728A (zh) | 一种建立有机高粱产量预测模型的方法 | |
| WO2024099475A1 (zh) | 基于长时序ndvi的土壤重金属胁迫甄别方法 | |
| CN101430735B (zh) | 一种保护性耕作模式选择方法 | |
| CN111223002B (zh) | 一种玉米区域干物质产量或青储产量评估方法和*** | |
| CN116151454A (zh) | 一种多光谱无人机预测矮林芳樟精油产量的方法及*** | |
| Bhimani et al. | Forecasting of groundnut yield using meteorological variables | |
| Sudianto et al. | Early warning for sugarcane growth using phenology-based remote sensing by region | |
| CN111626638A (zh) | 夏玉米倒伏气象等级评估模型的构建与应用 | |
| Ji et al. | Using NDVI time series curve change rate to estimate winter wheat yield | |
| Garde et al. | Weather based pre-harvest forecasting of wheat at Ghazipur (UP) | |
| Mech | Status of tea production in Assam: past trends and its future projections | |
| Gill et al. | Development of agrometeorological models for estimation of cotton yield | |
| CN118410896A (zh) | 一种基于作物长势的水稻产量预测方法及装置 | |
| CN117744898B (zh) | 一种大田粮食作物产量年际预测模型构建方法 | |
| Zhang | Monitoring agricultural behavior under climate change with cloud computing and satellite imagery |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20240314 Address after: 100125 Building 2, national agricultural exhibition hall, No. 16, North East Third Ring Road, Chaoyang District, Beijing Patentee after: Big data development center of the Ministry of agriculture and rural areas Country or region after: China Address before: 100125, No. 41, wheat Street, Chaoyang District, Beijing Patentee before: Ministry of Agriculture and Rural Planning and Design Institute Country or region before: China |
|
| TR01 | Transfer of patent right |