CN110309985B - 一种农作物产量预测方法及系统 - Google Patents
一种农作物产量预测方法及系统 Download PDFInfo
- Publication number
- CN110309985B CN110309985B CN201910618693.3A CN201910618693A CN110309985B CN 110309985 B CN110309985 B CN 110309985B CN 201910618693 A CN201910618693 A CN 201910618693A CN 110309985 B CN110309985 B CN 110309985B
- Authority
- CN
- China
- Prior art keywords
- yield
- model
- evi
- comprehensive
- data
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 42
- 230000012010 growth Effects 0.000 claims abstract description 107
- 241000209140 Triticum Species 0.000 claims abstract description 53
- 235000021307 Triticum Nutrition 0.000 claims abstract description 53
- 238000004519 manufacturing process Methods 0.000 claims abstract description 24
- 238000011160 research Methods 0.000 claims description 18
- 238000001556 precipitation Methods 0.000 claims description 15
- 230000003698 anagen phase Effects 0.000 claims description 9
- 238000010276 construction Methods 0.000 claims description 9
- 238000012911 target assessment Methods 0.000 claims description 6
- 238000010219 correlation analysis Methods 0.000 claims description 4
- 238000007405 data analysis Methods 0.000 claims description 3
- 230000008859 change Effects 0.000 description 18
- 238000011161 development Methods 0.000 description 12
- 230000018109 developmental process Effects 0.000 description 12
- 230000002159 abnormal effect Effects 0.000 description 10
- 230000000694 effects Effects 0.000 description 8
- 238000009331 sowing Methods 0.000 description 8
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 6
- 230000008878 coupling Effects 0.000 description 5
- 238000010168 coupling process Methods 0.000 description 5
- 238000005859 coupling reaction Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 230000007547 defect Effects 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 239000002689 soil Substances 0.000 description 4
- 210000004291 uterus Anatomy 0.000 description 4
- 241000238631 Hexapoda Species 0.000 description 3
- 235000013339 cereals Nutrition 0.000 description 3
- 238000012795 verification Methods 0.000 description 3
- 241000196324 Embryophyta Species 0.000 description 2
- 241000607479 Yersinia pestis Species 0.000 description 2
- 239000003513 alkali Substances 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 239000008267 milk Substances 0.000 description 2
- 210000004080 milk Anatomy 0.000 description 2
- 235000013336 milk Nutrition 0.000 description 2
- 238000010899 nucleation Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 239000000443 aerosol Substances 0.000 description 1
- 230000009418 agronomic effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000008014 freezing Effects 0.000 description 1
- 238000007710 freezing Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 239000003621 irrigation water Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 239000000047 product Substances 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000007637 random forest analysis Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 229920006395 saturated elastomer Polymers 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 239000013589 supplement Substances 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/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
-
- 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
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Business, Economics & Management (AREA)
- Engineering & Computer Science (AREA)
- Strategic Management (AREA)
- Human Resources & Organizations (AREA)
- Economics (AREA)
- Theoretical Computer Science (AREA)
- Marketing (AREA)
- General Physics & Mathematics (AREA)
- General Business, Economics & Management (AREA)
- Physics & Mathematics (AREA)
- Tourism & Hospitality (AREA)
- Animal Husbandry (AREA)
- Life Sciences & Earth Sciences (AREA)
- Agronomy & Crop Science (AREA)
- Mining & Mineral Resources (AREA)
- General Health & Medical Sciences (AREA)
- Marine Sciences & Fisheries (AREA)
- Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Development Economics (AREA)
- Game Theory and Decision Science (AREA)
- Entrepreneurship & Innovation (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种农作物产量预测方法，所述预测方法，首先，建立了基于趋势产量模型、综合气候指数模型和生长关键期EVI的包含未知系数的综合估产模型，然后利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型；最后，获取目标估产年份的气象数据和EVI数据，并将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。本发明构建了基于趋势产量、综合气候指数和生长关键期EVI的综合估产模型，以从更全面的角度实现冬小麦全生育期内的产量估算，提高冬小麦的估产精度。
Description
技术领域
本发明涉及农作物管理领域，特别涉及一种农作物产量预测方法及系统。
背景技术
小麦是世界三大主要粮食作物之一，其产量直接关系到各国的粮食安全，因此，对小麦生长环境、过程的监测、管理一直以来是各国学者关注的焦点问题，有效的产量估算则是各国政府决策管理、资源调控甚至粮食贸易的重要参考依据。
目前，小麦的估产方式主要有农学估产、农业气象估产、遥感估产和混合估产四种方式。
农学估产是通过模拟不同作物品种生长过程中水分、土壤、气候等影响因子，构建作物生长的模拟模型，即作物模型，通过作物模型估算作物产量，常见的作物模型有WOFOST模型、CASA模型等。作物模型充分考虑了作物本身的生理特征和外界环境的胁迫因素，能够揭示作物生长过程和产量形成的复杂机理，是一种理想的估产模型，在实际研究中也能获得较好的估产结果。但是，由于作物模型的复杂性，在大面积、复杂生产环境的应用中，往往无法满足其对大量的调优参数的获取要求，使得该类模型难以全面推广应用。
农业气象估产主要是利用影响作物生长发育的气象数据与产量的统计关系，结合反映生产力变化的趋势产量构建气象估产模型进行估产，是典型的经验半经验模型。气象估产模型的本质是以趋势产量为基础，引入易于获取的气象数据构建。模型简单易用，数据可获取性强，已有研究表明其具有一定的估产精度，是目前气象、农业、保险等行业普遍采用的常规估产方式。然而，由于气象估产模型只考虑了气象因子对估产周期内作物产量的影响，却忽略人类活动、自然灾害等随机因素的影响，使得模型在异常年份的估产往往偏离实际产量较大。另外，也有研究表明气象要素仅能解释产量变化的1/3，甚至更少，因此需要引入额外的辅助信息来提高气象估产模型精度。
遥感估产则是利用现代空间信息技术中高时空尺度的对地探测能力，以遥感反演的反映作物长势状态的植被指数与产量之间的统计关系，建立遥感植被指数模型实现产量估算。由于遥感反演的作物植被指数能够反映作物在各种环境下的生长状态，且具有覆盖范围广、周期性强等特点，植被指数获取及模型的实现过程也较为简单，因此，在大范围作物估产中得到了广泛应用。但是，因为遥感对地观测能力本身的缺陷，使得遥感估产存在局限，即当遥感反演的植被指数无法进一步反应作物的生长状态时，模型的估产能力就会下降。例如，周磊、Wang Quan的研究表明，当作物叶面积指数(LAI)达到4.0时，其归一化植被指数(NDVI)出现饱和，基于NDVI的估产模型就无法进一步反映植被的生长状况。此外遥感数据的不确定性也使模型不稳定，例如，云、云阴影、噪声以及遥感数据的滞后性等使得遥感数据不能准确反映作物的真实生长状态，影响估产精度。
综上所述，农学估产模型因其苛刻的应用条件而难以推广，气象估产模型、遥感估产模型因为在精度和稳定性上存在较大的不确定性而难以常规应用，而已有研究表明，在应用气象估产时，引入遥感反演的光谱信息来补充气象要素提供信息的不足可提高估产的精度。因此，近年来，使用多种数据、多种模型相耦合的方法构建估产模型成为热点。例如，王恺宁等组合冬小麦灌浆期四种不同植被指数建立模型，弥补单植被指数的缺陷以提高遥感估产的精度。Umer Saeed等使用随机森林方法结合气象数据、MODIS NDVI数据建立巴基斯坦Punjab省冬小麦预测模型。Manjunath、李付琴等[33]把遥感光谱信息和气象条件结合建立统一的模型，比单纯的遥感模型和气象模型都要优越，说明外界环境对产量的影响很大，在进行遥感估产时不能忽视气象条件的影响。这些模型简单易行，也取得较好的结果，却忽略了历史趋势对产量的影响，且气象因子仅选取几个相关性较高的时期，无法全面反映冬小麦全生育期内气象要素的变化对产量的影响。
发明内容
本发明的目的是提供一种农作物产量预测方法及系统，以从更全面的角度实现冬小麦全生育期内的产量估算，提高冬小麦的估产精度；本发明的预测系统和方法不仅限于对冬小麦的产量的预测，还可以实现其它农作物产量的预测。
为实现上述目的，本发明提供了如下方案：
本发明提供一种农作物产量预测方法，所述预测方法包括如下步骤：
建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型；
获取估产区域的历史数据，所述历史数据包括作物历史产量数据、作物历史空间分布数据、作物历史EVI数据和历史气象数据；
利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型；
获取目标估产年份的气象数据和EVI数据；
将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。
可选的，所述建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型，具体包括：
建立趋势产量模型；
建立综合气候指数模型；
根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型。
可选的，所述建立趋势产量模型，具体包括：
构建表示年份X与年份X对应的趋势产量Yt关系的三次多项式：Yt＝α*X3+β*X2+γ*X+φ，其中，α为三次项系数、β为二次项系数、γ为一次项系数和φ为常数项；
根据所述作物历史产量数据采用滑动平均法获得历史年份对应的历史趋势产量数据；
根据历史年份对应的历史趋势产量数据，求解所述三次多项式中的三次项系数α、二次项系数β、一次项系数γ和常数项φ，得到趋势产量模型。
可选的，所述建立综合气候指数模型，具体包括：
建立农作物的综合气候指数模型：
其中，M表示农作物的综合气候指数，n表示生育期数；F(ci)表示第i个生育期的综合气候指数，F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)，其中，wp，wt，ws分别是降水、温度、日照适宜度指数的权重系数；F(pi)，F(ti)，F(si)分别为第i个生育期降水、温度、日照适宜度指数；wi表示第i个生育期的权重系数，
Ri为第i个生育期的综合气候指数与波动产量的相关系数。
可选的，所述根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型，具体包括：
根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型：Y＝f(Yt,M,EVIm)＝aYt+bM+cEVIm+d，Y为估测单产；Yt为趋势产量；M为综合气候指数；EVIm为生长关键期EVI指数；a，b，c和d分别为第一系数项、第二系数项、第三系数项和第四系数项。
一种农作物产量预测系统，所述预测系统包括：
包含未知系数的综合估产模型建立模块，用于建立趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型；
历史数据获取模块，用于获取估产区域的历史数据，所述历史数据包括作物历史产量数据、作物历史空间分布数据、作物历史EVI数据和历史气象数据；
未知系数求解模块，用于利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型；
目标估产年份数据获取模块，用于获取目标估产年份的气象数据和EVI数据；
产量预测模块，用于将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。
可选的，所述包含未知系数的综合估产模型建立模块包括：
趋势产量模型建立子模块，用于建立趋势产量模型；
综合气候指数模型建立子模块，用于建立综合气候指数模型；
包含未知系数的综合估产模型建立子模块，用于根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型。
可选的，所述趋势产量模型建立子模块，具体包括：
三次多项式构建单元，用于构建表示年份X与年份X对应的趋势产量Yt关系的三次多项式：Yt＝α*X3+β*X2+γ*X+φ，其中，α为三次项系数、β为二次项系数、γ为一次项系数和φ为常数项；
历史产量数据分析单元，用于根据所述作物历史产量数据采用滑动平均法获得历史年份对应的历史趋势产量数据；
三次多项式系数求解单元，用于根据历史年份对应的历史趋势产量数据，求解所述三次多项式中的三次项系数α、二次项系数β、一次项系数γ和常数项φ，得到趋势产量模型。
可选的，所述综合气候指数模型建立子模块，具体包括：
综合气候指数模型建立单元，用于建立农作物的综合气候指数模型：
其中，M表示农作物的综合气候指数，n表示生育期数；F(ci)表示第i个生育期的综合气候指数，F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)，其中，wp，wt，ws分别是降水、温度、日照适宜度指数的权重系数；F(pi)，F(ti)，F(si)分别为第i个生育期降水、温度、日照适宜度指数；wi表示第i个生育期的权重系数，
Ri为第i个生育期的综合气候指数与波动产量的相关系数。
可选的，所述包含未知系数的综合估产模型建立子模块，具体包括：
包含未知系数的综合估产模型建立单元，用于根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型：Y＝f(Yt,M,EVIm)＝aYt+bM+cEVIm+d，Y为估测单产；Yt为趋势产量；M为气候指数；EVIm为生长关键期EVI指数；a，b，c和d分别为第一系数项、第二系数项、第三系数项和第四系数项。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明提供一种农作物产量预测方法，所述预测方法，首先，建立了基于趋势产量模型、综合气候指数模型和生长关键期EVI的包含未知系数的综合估产模型，然后利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型；最后，获取目标估产年份的气象数据和EVI数据，并将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。本发明构建了基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的综合估产模型，以从更全面的角度实现冬小麦全生育期内的产量估算，提高冬小麦的估产精度，为行业部门进行冬小麦产量的业务化估算提供一种新的、简单易行的参考方案。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明提供的农作物产量预测方法的流程图；
图2为本发明提供的求解包含未知系数的综合估产模型中的未知系数的原理图；
图3为本发明提供的农作物产量预测系统的结构图；
图4为本发明提供的实验研究区(估产区)的概况图；
图5为本发明提供的不同模型的的预测结果对比图；
图6为本发明提供的不同模型的拟合精度对比图；
图7为本发明提供的波动产量异常值分布图；
图8为本发明提供的不同模型的异常年份预测结果对比图。
具体实施方式
本发明的目的是提供一种农作物产量预测方法及系统，以从更全面的角度实现冬小麦全生育期内的产量估算，提高冬小麦的估产精度。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。
本研究考虑作物最终产量受气候条件、土壤属性、作物种子遗传、田间管理以及政策导向因素等多方面综合影响，使用相对简单的统计方法，综合考虑社会因素、自然因素与冬小麦产量的关系，以趋势模型结果作为产量变化基础，引入气候指数和遥感反演的增强型植被指数(Enhanced Vegetation Index，EVI)分别反映气候变化和环境综合变化对产量造成的胁迫，构建由趋势产量、气候指数以及EVI耦合的估产模型，以期从更全面的角度实现冬小麦全生育期内的产量估算，提高冬小麦的估产精度，为行业部门进行冬小麦产量的业务化估算提供一种新的、简单易行的参考方案，具体包括：
如图1所示本发明提供一种农作物产量预测方法，所述预测方法包括如下步骤：
步骤101，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型；
具体包括：
建立趋势产量模型；构建表示年份X与年份X对应的趋势产量Yt关系的三次多项式：Yt＝α*X3+β*X2+γ*X+φ，其中，α为三次项系数、β为二次项系数、γ为一次项系数和φ为常数项；根据所述作物历史产量数据采用滑动平均法获得历史年份对应的历史趋势产量数据；根据历史年份对应的历史趋势产量数据，求解所述三次多项式中的三次项系数α、二次项系数β、一次项系数γ和常数项φ，得到趋势产量模型。
建立综合气候指数模型：建立农作物的综合气候指数模型：
其中，M表示农作物的综合气候指数，n表示生育期数；F(ci)表示第i个生育期的综合气候指数，F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)，其中，wp，wt，ws分别是降水、温度、日照适宜度指数的权重系数；F(pi)，F(ti)，F(si)分别为第i个生育期降水、温度、日照适宜度指数；wi表示第i个生育期的权重系数，
Ri为第i个生育期的综合气候指数与波动产量的相关系数。
建立包含未知系数的综合估产模型：根据所述趋势产量模型、综合气候指数模型和生长关键期EVI，建立包含未知系数的综合估产模型：Y＝f(Yt,M,EVIm)＝aYt+bM+cEVIm+d，Y为估测单产；Yt为趋势产量；M为综合气候指数；EVIm为生长关键期EVI；a，b，c和d分别为第一系数项、第二系数项、第三系数项和第四系数项。
步骤102，获取估产区域的历史数据，所述历史数据包括作物历史产量数据、作物历史空间分布数据、作物历史EVI数据和历史气象数据；
步骤103，利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型。
以冬小麦为例，如图2所示，具体包括：
1)根据作物历史产量数据，确定包含未知系数的综合估产模型的趋势产量模型的值。
2)利用公式Yf＝Y-Yt，将趋势产量模型输出的历史趋势产量Yt与历史实际产量的差，作为历史波动产量Yf；
根据历史年份的历史气象数据确定每个生育期的综合气候指数：F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)，
其中，
式中：F(pi)为i生育期降水适宜度，pi为i生育期内总降水量，pai为相应生育期内多年降水量平均值。
式中：j为日序；F(tij)为第i个生长期的第j日温度适宜度；tj为第j日的日平均温度；til，tih，tio为冬小麦第i个生长发育所需要的最低、最高和最适温度，其值见表1。当温度适宜度计算值为负值时，则赋值为0，Ji表示第i个生长期的天数。
式中F(si)为i生育期日照适宜度，si为i生育期内总日照时数，S为i生育期内总日照时数的多年平均值。
表1冬小麦各生育期最低tl、最高th、最适温度to(单位摄氏度/℃)
根据冬小麦发育阶段内各生长期的综合气候指数与波动产量之间的相关系数来确定各生育期综合气候指数的权重系数，
Ri为第i个生育期的综合气候指数与波动产量的相关系数。
将每个生育期的综合气候指数和各生育期综合气候指数的权重系数，代入包含未知系数的综合估产模型的综合气候指数模型，计算综合气候指数值。
3)首先基于历史冬小麦空间分布数据、作物历史EVI数据提取估产区冬小麦各生长期的EVI指数，并取得其空间均值：
式中：
为研究区EVI平均值，n为研究区内冬小麦纯净像元的个数(生长期的个数)，EVIi为像元值(第i个生长期的EVI指数)。
其次通过相关性分析选用估产区与最终产量相关性最高时期的EVI均值作为生长关键期EVI指数。
4)并将作物历史产量数据代入包含未知系数的综合估产模型的综合气候指数模型。
获得多组历史数据，采用步骤1)-4)方式建立方程组，通过求解方程组求解系数。
步骤104，获取目标估产年份的气象数据和EVI数据；
步骤105，将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。
如图3所示，本发明还提供了一种农作物产量预测系统，所述预测系统包括：
包含未知系数的综合估产模型建立模块301，用于建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型；
历史数据获取模块302，用于获取估产区域的历史数据，所述历史数据包括作物历史产量数据、作物历史空间分布数据、作物历史EVI数据和历史气象数据；
未知系数求解模块303，用于利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型；
目标估产年份数据获取模块304，用于获取目标估产年份的气象数据和EVI数据；
产量预测模块305，用于将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。
所述包含未知系数的综合估产模型建立模块301包括：
趋势产量模型建立子模块，用于建立趋势产量模型；
综合气候指数模型建立子模块，用于建立综合气候指数模型；
包含未知系数的综合估产模型建立子模块，用于根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型。
其中，所述趋势产量模型建立子模块，具体包括：三次多项式构建单元，用于构建表示年份X与年份X对应的趋势产量Yt关系的三次多项式：Yt＝α*X3+β*X2+γ*X+φ，其中，α为三次项系数、β为二次项系数、γ为一次项系数和φ为常数项；历史产量数据分析单元，用于根据所述作物历史产量数据采用滑动平均法获得历史年份对应的历史趋势产量数据；三次多项式系数求解单元，用于根据历史年份对应的历史趋势产量数据，求解所述三次多项式中的三次项系数α、二次项系数β、一次项系数γ和常数项φ，得到趋势产量模型。
所述综合气候指数模型建立子模块，具体包括：综合气候指数模型建立单元，用于建立农作物的综合气候指数模型：
其中，M表示农作物的综合气候指数，n表示生育期数；F(ci)表示第i个生育期的综合气候指数，F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)，其中，wp，wt，ws分别是降水、温度、日照适宜度指数的权重系数；F(pi)，F(ti)，F(si)分别为第i个生育期降水、温度、日照适宜度指数；wi表示第i个生育期的权重系数，
Ri为第i个生育期的综合气候指数与波动产量的相关系数。
所述包含未知系数的综合估产模型建立子模块，具体包括：包含未知系数的综合估产模型建立单元，用于根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型：Y＝f(Yt,M,EVIm)＝aYt+bM+cEVIm+d，Y为估测单产；Yt为趋势产量模型；M为气候指数模型；EVIm为生长关键期EVI指数；a，b，c和d分别为第一系数项、第二系数项、第三系数项和第四系数项。
为了验证本发明的预测方法及系统的有效性，本发明还提供了一种验证方法。
所示本发明选取中国河北省11个县作为验证研究区(目标估产区)。如图4所示，河北省(113.85E-119.81E，36.1N-42.58N)，位于太行山东侧、华北平原北部，有坝上高原、燕山和太行山山地和河北平原三大地貌，地形复杂多样，属于温带大陆性季风气候，是中国北方适宜冬小麦生长的主产区，其物候期一般在10月中上旬至次年6月中下旬，主要发育期如表2所示。11个县的冬小麦种植环境如表3所示，本发明选取的11个县代表不同地形、不同管理水平下冬小麦的生长发育状况且气象观测条件较好以确保气象数据的准确性。
表2冬小麦发育期
| 发育期 | 时间 | 发育期 | 时间 |
| 播种期 | 10-01～10-10 | 拔节期 | 03-21～03-31 |
| 出苗期 | 10-11～10-20 | 抽穗期 | 04-01～04-20 |
| 分蘖期 | 10-21～11-30 | 乳熟期 | 04-21～05-10 |
| 越冬期 | 12-01～02-20 | 成熟期 | 05-11～06-10 |
| 返青期 | 02-21～03-20 |
表3研究区冬小麦种植环境
本发明所使用数据如表2所示，遥感数据为2001年至2016年16年冬小麦生育期内16天合成的MOD13Q1 EVI产品数据，空间分辨率为250m，并对EVI进行SG滤波处理，以消除云、气溶胶等噪声的影响。并基于冬小麦空间分布数据提取出对应的EVI，筛选出EVI纯净像元。
气象数据从气象数据共享中心获取，包括2001年至2016年11个气象站点，如图4所示，冬小麦生育期内的最高温度、最低温度、平均温度、降水量、以及日照时数。
产量数据为河北省统计年鉴2001年至2016年县级产量数据，包括种植面积和产量，通过种植面积和县级总产计算各县的单位产量(kg/ha)如表4所示。
表4统计单产(kg/ha)
首先，利用趋势产量模型，对河北省11个县2001-2016年的冬小麦单产数据使用5年直线滑动平均分离趋势产量如表5所示，其次，利用历史数据求解三次多项式的系数建立关于时间的趋势模型，结果如表6。
表5趋势产量(kg/ha)
表6趋势产量模型
| 县名 | 模型 | R<sup>2</sup> |
| 邢台县 | Y<sub>t</sub>＝0.3414*x<sup>3</sup>-2053*x<sup>2</sup>+4.117e<sup>6</sup>*x-2.751e<sup>9</sup> | 0.960 |
| 泊头市 | Y<sub>t</sub>＝1.711*x<sup>3</sup>-1.032e<sup>4</sup>*x<sup>2</sup>+2.074e<sup>7</sup>*x-1.39e<sup>10</sup> | 0.981 |
| 饶阳县 | Y<sub>t</sub>＝-0.2615*x<sup>3</sup>+1577*x<sup>2</sup>-3.168e<sup>6</sup>*x+2.122e<sup>9</sup> | 0.922 |
| 霸州市 | Y<sub>t</sub>＝0.7914*x<sup>3</sup>-4774*x<sup>2</sup>+9.601e<sup>6</sup>*x-6.436e<sup>9</sup> | 0.966 |
| 乐亭县 | Y<sub>t</sub>＝0.3213*x<sup>3</sup>-1935*x<sup>2</sup>+3.883e<sup>6</sup>*x-2.598e<sup>9</sup> | 0.797 |
| 遵化市 | Y<sub>t</sub>＝0.6071*x<sup>3</sup>-3656*x<sup>2</sup>+7.339e<sup>6</sup>*x-4.911e<sup>9</sup> | 0.977 |
| 固安县 | Y<sub>t</sub>＝0.1109*x<sup>3</sup>-664.1*x<sup>2</sup>+1.326e<sup>6</sup>*x-8.827e<sup>9</sup> | 0.984 |
| 丰润区 | Y<sub>t</sub>＝-0.4352*x<sup>3</sup>+2626*x<sup>2</sup>-5.283e<sup>6</sup>*x+3.543e<sup>9</sup> | 0.886 |
| 南宫市 | Y<sub>t</sub>＝0.3401*x<sup>3</sup>-2051*x<sup>2</sup>+4.125e<sup>6</sup>*x-2.765e<sup>9</sup> | 0.997 |
| 徐水县 | Y<sub>t</sub>＝-0.8192*x<sup>3</sup>+4939*x<sup>2</sup>-9.926e<sup>6</sup>*x+6.649e<sup>9</sup> | 0.972 |
| 黄骅市 | Y<sub>t</sub>＝-0.9325*x<sup>3</sup>+5605*x<sup>2</sup>-1.123e<sup>7</sup>*x+7.5e<sup>9</sup> | 0.973 |
Yt为趋势产量，单位：kg/ha，x表示年份，如计算2016年趋势单产，则x＝2016。
首先以公式Yf＝Y-Yt得到11个县的波动产量，如表7所示。
表7波动产量(kg/ha)
| 年份 | 邢台县 | 泊头市 | 饶阳县 | 霸州 | 乐亭县 | 遵化市 | 固安县 | 丰润区 | 南宫市 | 徐水县 | 黄骅市 |
| 2001 | 43.989 | 69.796 | 47.504 | -191.675 | -136.684 | -216.224 | -28.618 | -178.392 | 56.296 | 112.205 | 861.109 |
| 2002 | 22.520 | -34.584 | -94.223 | -129.372 | -11.316 | -12.791 | -51.872 | 53.925 | -14.103 | 260.804 | -625.998 |
| 2003 | 114.310 | -158.347 | -173.075 | -29.145 | 22.294 | 77.913 | -8.117 | -39.925 | 3.753 | -672.159 | -480.030 |
| 2004 | 61.911 | 249.474 | 345.229 | 173.178 | 32.401 | 76.888 | 18.960 | 3.895 | 25.589 | 258.937 | 1140.688 |
| 2005 | -47.780 | 29.825 | -32.233 | 89.873 | -36.459 | -59.945 | 34.677 | -6.330 | 31.245 | 129.562 | -72.927 |
| 2006 | -89.823 | -53.033 | -37.734 | -11.974 | -49.555 | -13.299 | 25.101 | -69.882 | -40.021 | 75.218 | -986.241 |
| 2007 | 32.943 | -144.618 | -163.221 | -145.667 | 119.980 | -45.233 | -112.401 | 114.469 | -32.708 | -210.692 | 390.664 |
| 2008 | 39.456 | 109.550 | 10.994 | 2.077 | -54.088 | 31.390 | 51.571 | -33.269 | 18.361 | 86.040 | 379.163 |
| 2009 | -12.944 | 138.940 | 100.614 | 83.421 | 1.170 | 87.794 | 21.240 | 30.089 | 42.031 | 28.804 | -162.175 |
| 2010 | -39.216 | -121.249 | 27.995 | 54.722 | 36.076 | -34.312 | 20.074 | -3.617 | 1.457 | -59.762 | -10.378 |
| 2011 | 80.308 | -49.303 | -31.603 | -46.725 | -31.652 | -42.382 | -15.672 | -15.262 | -17.354 | 7.057 | -15.123 |
| 2012 | -28.790 | 28.855 | -4.760 | -34.807 | -8.113 | -32.455 | -17.604 | -4.061 | -2.681 | 30.469 | -28.753 |
| 2013 | -79.704 | -10.880 | -26.009 | -40.577 | -41.553 | 39.558 | -31.703 | -59.168 | -53.366 | -26.382 | -35.023 |
| 2014 | 66.396 | 15.910 | 34.764 | 52.664 | 57.300 | 14.799 | 41.897 | 60.528 | 58.674 | 45.030 | 89.797 |
| 2015 | 42.368 | 35.575 | -1.736 | 50.493 | 25.019 | 28.104 | 26.189 | 42.792 | 29.912 | -15.135 | -22.514 |
| 2016 | -53.652 | -35.912 | -18.841 | -51.741 | -38.817 | -52.179 | -32.177 | -53.356 | -51.678 | -44.407 | -31.434 |
其次，获取各生育期的总降水量、每日平均温度、个生育期内总日照时数，计算得到各生育期降水、温度、日照适宜度指数F(pi)，F(ti)，F(si)，采用通径分析法计算同一生育期内这三者分别对波动产量的直接影响，即直接效应，并将直接效应归一化后作为该生育期内降水、温度、日照适宜度的权重。结果如表8所示。
表8不同生长期的适宜度指数
然后，基于公式F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)求取各区县的不同生育期F(ci),并与波动产量求取相关系数Ri，结果如表9所示，最后根据公式
求取各生育期F(ci)的权重Wi，结果如表10所示。
表9各县各生育期气候指数与波动产量相关性Ri
| 发育期 | 泊头市 | 饶阳县 | 霸州市 | 遵化市 | 固安县 | 丰润区 | 南宫市 | 徐水区 |
| 播种期 | -0.384 | -0.204 | -0.248 | 0.418 | -0.155 | 0.201 | -0.301 | 0.277 |
| 播种至出苗期 | -0.433 | -0.283 | 0.540<sup>*</sup> | 0.404 | 0.637<sup>**</sup> | -0.485 | -0.330 | 0.535<sup>*</sup> |
| 播种至分蘖期 | -0.490 | 0.335 | 0.442 | 0.411 | 0.582<sup>*</sup> | -0.395 | -0.364 | 0.701<sup>**</sup> |
| 播种至越冬期 | 0.531<sup>*</sup> | -0.299 | 0.443 | 0.399 | 0.591<sup>*</sup> | -0.497<sup>*</sup> | -0.401 | 0.686<sup>**</sup> |
| 播种至返青期 | 0.678<sup>**</sup> | 0.576<sup>*</sup> | 0.606<sup>*</sup> | 0.583<sup>*</sup> | 0.593<sup>*</sup> | -0.507<sup>*</sup> | 0.765<sup>**</sup> | 0.863<sup>**</sup> |
| 播种至拔节期 | 0.644<sup>**</sup> | 0.571<sup>*</sup> | 0.635<sup>**</sup> | 0.520<sup>*</sup> | 0.599<sup>*</sup> | -0.631<sup>**</sup> | 0.816<sup>**</sup> | 0.860<sup>**</sup> |
| 播种至抽穗期 | 0.707<sup>**</sup> | 0.598<sup>*</sup> | 0.635<sup>**</sup> | 0.621<sup>*</sup> | 0.594<sup>*</sup> | -0.601<sup>*</sup> | 0.809<sup>**</sup> | 0.842<sup>**</sup> |
| 播种至乳熟期 | 0.740<sup>**</sup> | 0.672<sup>**</sup> | 0.634<sup>**</sup> | 0.646<sup>**</sup> | 0.666<sup>**</sup> | -0.624<sup>**</sup> | 0.804<sup>**</sup> | 0.823<sup>**</sup> |
| 播种至成熟期 | 0.742<sup>**</sup> | 0.721<sup>**</sup> | 0.747<sup>**</sup> | 0.747<sup>**</sup> | 0.688<sup>**</sup> | 0.775<sup>**</sup> | 0.808<sup>**</sup> | 0.826<sup>**</sup> |
*表示通过置信度0.05的显著性水平检验，**表示通过置信度0.01的显著性水平检验。
表10各生育期F(ci)的权重Wi
最后即可通过公式
计算得到M。根据式
将各县冬小麦生育期内每16天获取的县级纯净EVI均值与统计产量做相关性分析，筛选与统计产量相关性最高的时期的EVI，其结果如表11所示：
表11选取的EVI时期及相关性
| 县 | 时期 | R | 县 | 时期 | R |
| 邢台县 | 05-09 | 0.669<sup>**</sup> | 固安县 | 04-07 | 0.597<sup>*</sup> |
| 泊头市 | 05-09 | 0.723<sup>**</sup> | 丰润区 | 05-25 | 0.663<sup>**</sup> |
| 饶阳县 | 05-09 | 0.66<sup>**</sup> | 南宫市 | 05-25 | 0.839<sup>**</sup> |
| 霸州市 | 05-25 | 0.858<sup>**</sup> | 徐水县 | 12-19 | 0.612<sup>*</sup> |
| 乐亭县 | 04-07 | 0.624<sup>**</sup> | 黄骅市 | 05-09 | 0.931<sup>**</sup> |
| 遵化市 | 01-17 | 0.719<sup>**</sup> |
*表示通过置信度0.05的显著性水平检验，**表示通过置信度0.01的显著性水平检验。
使用2001～2016年的统计产量与趋势产量、气候指数、关键期EVI，通过多元回归的方法确定本研究设计的耦合估产模型参数，结果如表12所示。
表12耦合估产模型
**表示模型通过置信度0.01的显著性水平检验，R2 adj代表调节R2,Yt为趋势模型产量，M为气候指数，EVI为最佳估产期EVI。
研究比较11个研究县的不同模型预测结果与统计值，如图5所示，深色实线为统计值，浅色实线是仅使用趋势构建模型的预测结果，深色虚线是使用趋势、气候指数构建模型的预测结果，浅色虚线是使用趋势、气候指数、EVI构建模型的预测结果。结果表明在相邻年份产量变化相对平稳的地区如南宫市、霸州市，三种模型预测值与统计值相差不大，这是由于田间管理、品种改良等人工影响使得气候变化对产量的影响很小且并未发生异常天气或灾害情况，因此仅使用趋势构建模型(以下称趋势模型)就可以得到很好的结果，但是对于气候变化或其他因素导致产量突然降低或升高的年份，趋势模型存在严重的高估或低估现象如固安县、饶阳县、乐亭县等县。以趋势、气候指数构建的模型(以下称趋势气候模型)和本研究模型(以下称趋势气候EVI模型)虽然也存在高估和低估现象，却有效的改善了这种状况，这是因为趋势气候模型在趋势模型的基础上考虑了气候变化对产量丰欠的影响，而趋势气候EVI模型则在趋势气候模型基础上加入了其他因素的影响。因此趋势气候EVI模型在受病虫害、干热风等其他因素影响的年份预测结果要优于趋势气候模型，如黄骅市、徐水县因灾害产量突然降低的年份。
判定系数(R2)可以用来评价模型的拟合优度但会随自变量增加而变大，而调整的判定系数(R2 adj)同时考虑了样本量与模型中自变量的个数，不会随模型中自变量个数的增加而接近1。均方根误差(RMSE)通过计算预测值与观测值之间的偏差来估算模型的精度，其值越小模型越好。
因此实验比较各模型的R2 adj和RMSE来评价各模型的拟合精度，结果如图6所示。表明使用趋势构建的模型虽然具有一定的预测效果，但是在不同地区表现却相差很大，该模型在南宫市R2 adj达到0.95以上，在黄骅市、徐水县却不到0.7，RMSE变化区间在50～589kg/ha，说明使用趋势构建的模型在应用时不稳定，鲁棒性较差。趋势气候模型和趋势气候EVI模型的R2 adj最低分别是0.845和0.838，最高都是0.993，RMSE变化区间分别为24.4～342.1kg/ha和22.1～237.1kg/ha，表明在模型中加入气候指数不仅提高了模型的拟合优度还降低了RMSE，R2 adj从0.7提高到0.845，RMSE从589kg/ha降低到342.1kg/ha，有效地提高了估产模型的稳定性和鲁棒性。而加入了EVI的模型虽然在各县的RMSE要低于趋势气候模型，但其R2 adj除黄骅市外未有显著提升，如饶阳县、乐亭县。这是由于黄骅市遭受干热风影响，而气候指数仅考虑了温度、光照、降水对产量的影响，没有考虑干热风的因素，这些因素对产量的影响却在EVI信息上表现出来，因此趋势气候EVI模型的R2 adj在黄骅市最高，说明趋势气候EVI模型在异常年份要优于趋势气候模型。
为了进一步比较模型在异常年份的预测结果，实验使用箱式图选出各县波动产量异常的值，如图7所示，并比较模型在这些异常年份的预测结果与统计值，如图8所示。如图8所示的结果表明趋势模型在异常年份拟合优度仅达到0.855，而其他两个模型精度均在0.95以上，优于趋势模型，其中趋势气候EVI模型预测结果与统计值拟合优度最高，达到了0.98，比趋势气候模型高0.02。趋势气候EVI模型的RMSE为226kg/ha，比趋势气候模型的RMSE低140kg/ha，比趋势模型的RMSE低近500kg/ha，其拟合直线斜率和截距也是最小的，分别为0.974和114，优于趋势气候模型拟合直线斜率0.956和181，说明趋势气候EVI耦合的模型预测结果更接近于统计值。
总体上，三种模型在各县表现效果排序为趋势气候EVI模型>>趋势气候模型>>趋势模型，反映出趋势气候EVI模型预测结果更接近于实测的统计值，特别是对于异常年份产量的预测。但是对于非异常年份的预测，相比于趋势气候模型，趋势气候EVI模型预测效果提升很小，而模型的拟合优度甚至略有降低，可能是由于MODIS数据混合像元的影响，EVI不能准确反映作物生长状态。
本研究提出的趋势、气候指数、EVI耦合模型综合考虑了不同模型的特点，趋势模型可以稳定反应冬小麦单产在时间尺度上的变化，却很难体现出气候或其他因素造成的单产波动，而气候指数和EVI可以表达空间尺度上气候状况和冬小麦生长状态变化引起的单产波动。与已有的研究方法相比，本模型考虑了更多的影响因素，相比于传统的趋势模型，本模型可以反映出气候变化或其他因素导致的单产变化，而同时使用气候和光谱信息表达单产的变化更加准确，相比于遥感模型和气象模型，本模型能更好的反映出单产变化的整体趋势。因此本模型在各县都能得到不错的结果，可用于大范围区域的产量估算。
为了验证模型的可用性，本研究以2017年统计产量为真值，将趋势模型、趋势气候模型、趋势气候EVI模型预测结果进行对比，结果见表13。由于缺少2017年丰润区单产数据，本研究未做2017年丰润区的验证。
表中三种模型在10个验证区域的精度均达到0.9以上，得到了非常好的效果。其中，趋势模型在徐水区精度最高(精度0.999)，由于该区的水域面积大、水资源丰富，耕地基本都是水浇地(占耕地的90％)，气象因素中的降水的多或少对该地区冬小麦的生长发育影响很小，而降水量的多少可能会导致模型估算结果偏高或偏低。因此，仅使用趋势模型就能得到很好的估算精度，其他两种加入气候指数、EVI反而会因为自身精度的问题不能准确表达这些因素的影响从而造成模型精度降低。但趋势模型在其他区域低于其他两种模型，特别是在遵化市(精度0.922)、黄骅市(精度0.942)、乐亭县(精度0.946)明显偏低，这是由于其他地区缺水干旱，灌溉水源不足导致降水等因素对冬小麦的生长发育影响很大，其产量会随气象因素的变化而出现大幅波动。而引入气候指数的趋势气候模型弥补了趋势模型的缺点，使模型估产精度得到提升，其中在遵化县提升最大，为5个百分点，说明了模型中加入气候指数的重要性。而引入EVI的耦合模型使得其中5个县在趋势气候模型精度的基础上又得到了提升，其中在黄骅市提升最大，为2.5个百分点，这是由于黄骅市农作物主要是“靠天收”，降水量对冬小麦影响较大，受盐碱地等自然资源的限制且易受干热风、冻害、虫害等灾害影响。乐亭县所处的自然资源条件与黄骅市相似，同样也是盐碱地，受干旱影响严重，耦合模型在该县表现也较好，说明了趋势气候模型中加入EVI构建模型可以弥补趋势模型、气候指数无法估算受病虫害等其他非气象因素对冬小麦发育状况影响的不足。邢台县、霸州市、固安县、遵化市、南宫市、泊头市、饶阳县主要是受降水等气象因素影响，受其他灾害性因素影响较小，因此，在对于这些地区，趋势气候模型与耦合模型表现相似，并无太大差别，其中在饶阳县和遵化市精度相差最大均为0.005，而这两个模型精度在这些地区均高于趋势模型。
表13冬小麦不同估产模型预测结果
| 县名 | 1 | 2 | 3 | 4 | 2<sup>*</sup> | 3<sup>*</sup> | 4<sup>*</sup> |
| 邢台县 | 5238.078 | 5383.442 | 5440.720 | 5382.779 | 0.972 | 0.961 | 0.972 |
| 泊头市 | 5751.123 | 5992.192 | 6040.703 | 5978.793 | 0.958 | 0.949 | 0.960 |
| 饶阳县 | 6220.804 | 6146.639 | 6383.222 | 6181.669 | 0.988 | 0.973 | 0.993 |
| 霸州市 | 5735.274 | 5717.380 | 5953.158 | 5715.461 | 0.996 | 0.962 | 0.996 |
| 乐亭县 | 5852.689 | 6000.225 | 6167.768 | 6060.019 | 0.975 | 0.946 | 0.964 |
| 遵化市 | 5216.412 | 5338.665 | 5625.858 | 5360.559 | 0.977 | 0.922 | 0.972 |
| 固安县 | 5954.992 | 5994.816 | 6148.148 | 5993.660 | 0.993 | 0.967 | 0.993 |
| 南宫市 | 5561.002 | 5648.011 | 5644.134 | 5652.050 | 0.984 | 0.985 | 0.983 |
| 徐水区 | 6624.797 | 6801.861 | 6623.128 | 6823.519 | 0.973 | 0.999 | 0.970 |
| 黄骅市 | 3291.573 | 3297.351 | 3103.419 | 3375.590 | 0.998 | 0.942 | 0.974 |
1为真实值(kg/ha)；2为趋势气候EVI模型结果(kg/ha)；3为趋势模型结果(kg/ha)；4为趋势气候模型结果(kg/ha)；2*为趋势气候EVI模型精度；3*为趋势模型精度；4*为趋势气候模型精度。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明提供一种农作物产量预测方法及系统，本发明的趋势、气候指数、EVI耦合的估产模型在县级尺度上取得较高的精度，充分考虑了时空尺度上冬小麦单产变化及影响因素，比仅使用趋势构建的模型更加稳定，特别是在产量异常年份的模型结果优于其他没有使用遥感的两种模型，对冬小麦估产具有一定的指导意义。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
Claims (6)
1.一种农作物产量预测方法，其特征在于，所述预测方法包括如下步骤：
建立基于趋势产量、综合气候指数和生长关键期EVI指数的包含未知系数的综合估产模型；
农作物的综合气候指数模型为：
其中，M表示农作物的综合气候指数，n表示生育期数；F(ci)表示第i个生育期的综合气候指数，F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)，其中，wp，wt，ws分别是降水、温度、日照适宜度指数的权重系数；F(pi)，F(ti)，F(si)分别为第i个生育期降水、温度、日照适宜度指数；wi表示第i个生育期的权重系数，
Ri为第i个生育期的综合气候指数与波动产量的相关系数,
pi为i生育期内总降水量，pai为相应生育期内多年降水量平均值；
包含未知系数的综合估产模型：
Y＝f(Yt,M,EVIm)＝aYt+bM+cEVIm+d，Y为估测单产；Yt为趋势产量；M为综合气候指数；EVIm为生长关键期EVI指数；a，b，c和d分别为第一系数项、第二系数项、第三系数项和第四系数项；
基于历史冬小麦空间分布数据、作物历史EVI数据提取估产区冬小麦各生长期的EVI指数，并取得其空间均值：
式中：
为研究区EVI平均值，n为研究区内冬小麦纯净像元的个数，即生长期的个数，EVIi为像元值，即第i个生长期的EVI指数；
通过相关性分析选用估产区与最终产量相关性最高时期的EVI均值作为生长关键期EVI指数；
获取估产区域的历史数据，所述历史数据包括作物历史产量数据、作物历史空间分布数据、作物历史EVI数据和历史气象数据；
利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型；
获取目标估产年份的气象数据和EVI数据；
将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。
2.根据权利要求1所述的农作物产量预测方法，其特征在于，所述建立基于趋势产量、综合气候指数和生长关键期EVI指数的包含未知系数的综合估产模型，具体包括：
建立趋势产量模型；
建立综合气候指数模型；
根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型。
3.根据权利要求2所述的农作物产量预测方法，其特征在于，所述建立趋势产量模型，具体包括：
构建表示年份X与年份X对应的趋势产量Yt关系的三次多项式：Yt＝α*X3+β*X2+γ*X+φ，其中，α为三次项系数、β为二次项系数、γ为一次项系数和φ为常数项；
根据所述作物历史产量数据采用滑动平均法获得历史年份对应的历史趋势产量数据；
根据历史年份对应的历史趋势产量数据，求解所述三次多项式中的三次项系数α、二次项系数β、一次项系数γ和常数项φ，得到趋势产量模型。
4.一种农作物产量预测系统，其特征在于，所述预测系统包括：
包含未知系数的综合估产模型建立模块，用于建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型；
农作物的综合气候指数模型为：
其中，M表示农作物的综合气候指数，n表示生育期数；F(ci)表示第i个生育期的综合气候指数，F(ci)＝wp*F(pi)+wt*F(ti)+ws*F(si)，其中，wp，wt，ws分别是降水、温度、日照适宜度指数的权重系数；F(pi)，F(ti)，F(si)分别为第i个生育期降水、温度、日照适宜度指数；wi表示第i个生育期的权重系数，
Ri为第i个生育期的综合气候指数与波动产量的相关系数,
pi为i生育期内总降水量，pai为相应生育期内多年降水量平均值；
包含未知系数的综合估产模型：
Y＝f(Yt,M,EVIm)＝aYt+bM+cEVIm+d，Y为估测单产；Yt为趋势产量；M为气候指数；EVIm为生长关键期EVI指数；a，b，c和d分别为第一系数项、第二系数项、第三系数项和第四系数项；
基于历史冬小麦空间分布数据、作物历史EVI数据提取估产区冬小麦各生长期的EVI指数，并取得其空间均值：
式中：
为研究区EVI平均值，n为研究区内冬小麦纯净像元的个数，即生长期的个数，EVIi为像元值，即第i个生长期的EVI指数；
通过相关性分析选用估产区与最终产量相关性最高时期的EVI均值作为生长关键期EVI指数；
历史数据获取模块，用于获取估产区域的历史数据，所述历史数据包括作物历史产量数据、作物历史空间分布数据、作物历史EVI数据和历史气象数据；
未知系数求解模块，用于利用所述历史数据求解包含未知系数的综合估产模型中的未知系数，得到综合估产模型；
目标估产年份数据获取模块，用于获取目标估产年份的气象数据和EVI数据；
产量预测模块，用于将所述气象数据、所述EVI数据输入所述综合估产模型，预测目标估产年份的产量。
5.根据权利要求4所述的农作物产量预测系统，其特征在于，所述包含未知系数的综合估产模型建立模块包括：
趋势产量模型建立子模块，用于建立趋势产量模型；
综合气候指数模型建立子模块，用于建立综合气候指数模型；
包含未知系数的综合估产模型建立子模块，用于根据所述趋势产量模型和所述综合气候指数模型，建立基于趋势产量模型、综合气候指数模型和生长关键期EVI指数的包含未知系数的综合估产模型。
6.根据权利要求5所述的农作物产量预测系统，其特征在于，所述趋势产量模型建立子模块，具体包括：
三次多项式构建单元，用于构建表示年份X与年份X对应的趋势产量Yt关系的三次多项式：Yt＝α*X3+β*X2+γ*X+φ，其中，α为三次项系数、β为二次项系数、γ为一次项系数和φ为常数项；
历史产量数据分析单元，用于根据所述作物历史产量数据采用滑动平均法获得历史年份对应的历史趋势产量数据；
三次多项式系数求解单元，用于根据历史年份对应的历史趋势产量数据，求解所述三次多项式中的三次项系数α、二次项系数β、一次项系数γ和常数项φ，得到趋势产量模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910618693.3A CN110309985B (zh) | 2019-07-10 | 2019-07-10 | 一种农作物产量预测方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910618693.3A CN110309985B (zh) | 2019-07-10 | 2019-07-10 | 一种农作物产量预测方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110309985A CN110309985A (zh) | 2019-10-08 |
| CN110309985B true CN110309985B (zh) | 2022-05-03 |
Family
ID=68079797
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910618693.3A Active CN110309985B (zh) | 2019-07-10 | 2019-07-10 | 一种农作物产量预测方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110309985B (zh) |
Families Citing this family (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111062526B (zh) * | 2019-12-09 | 2022-03-08 | 北京师范大学 | 一种冬小麦单产预测方法及系统 |
| CN111191543A (zh) * | 2019-12-20 | 2020-05-22 | 湖南城市学院 | 油菜产量估测方法 |
| CN111260241B (zh) * | 2020-01-23 | 2021-02-02 | 北京师范大学 | 一种指数型农业保险产品设计方法及其产品 |
| CN111612236A (zh) * | 2020-05-14 | 2020-09-01 | 中电工业互联网有限公司 | 一种虫情实时分析预测方法、系统及存储介质 |
| CN111639803A (zh) * | 2020-05-29 | 2020-09-08 | 福州市规划设计研究院 | 一种应用于气候变化情景下区域未来植被指数的预估方法 |
| CN111797835B (zh) * | 2020-06-01 | 2024-02-09 | 深圳市识农智能科技有限公司 | 一种病症识别方法、病症识别装置及终端设备 |
| CN112183428A (zh) * | 2020-10-09 | 2021-01-05 | 浙江大学中原研究院 | 一种小麦种植区域分割与产量预测方法 |
| CN112840977A (zh) * | 2020-12-31 | 2021-05-28 | 航天信息股份有限公司 | 一种基于小麦关键生育时期预测小麦产量的方法及系统 |
| CN112990692B (zh) * | 2021-03-11 | 2023-05-30 | 中国水利水电科学研究院 | 一种基于分生育阶段多参数耦合的冬小麦产量估计方法 |
| CN113255999B (zh) * | 2021-05-25 | 2024-04-12 | 中国农业科学院农业信息研究所 | 基于多模型的农作物产量评估方法、装置、设备及介质 |
| CN113222294B (zh) * | 2021-06-07 | 2023-09-26 | 中国农业科学院农业信息研究所 | 一种大豆投入单产预测方法及系统 |
| CN113657781A (zh) * | 2021-08-23 | 2021-11-16 | 北京师范大学 | 一种适用于极端气候条件下的小麦估产方法及系统 |
| CN114118679B (zh) * | 2021-10-14 | 2022-09-16 | 农业农村部规划设计研究院 | 一种基于时序遥感数据的作物单产及长势评估方法 |
| CN114529091A (zh) * | 2022-02-22 | 2022-05-24 | 山东省农业科学院 | 一种融合气象数据的作物产量预测系统 |
| CN114331753B (zh) * | 2022-03-04 | 2022-06-14 | 阿里巴巴达摩院(杭州)科技有限公司 | 农事智能化的方法、装置及控制设备 |
| CN115545311A (zh) * | 2022-10-11 | 2022-12-30 | 中国科学院西北生态环境资源研究院 | 一种农作物估产方法、装置、存储介质及电子设备 |
| CN116821816B (zh) * | 2023-05-17 | 2024-05-28 | 河南农业大学 | 一种基于加权随机森林的干热风预测方法 |
| CN117634225B (zh) * | 2024-01-24 | 2024-04-30 | 中国农业科学院农业资源与农业区划研究所 | 基于作物生长通用模型和曲线拟合法的油料作物产量模拟调参方法和系统 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108154193B (zh) * | 2018-01-16 | 2021-10-08 | 黄河水利委员会黄河水利科学研究院 | 一种长时间序列降水数据降尺度方法 |
| CN108898499B (zh) * | 2018-05-31 | 2021-08-31 | 河南省气象科学研究所 | 一种冬小麦干热风保险天气指数计算方法 |
| CN108921351A (zh) * | 2018-07-06 | 2018-11-30 | 北京兴农丰华科技有限公司 | 基于趋势产量和气象产量的作物产量预测方法 |
| CN109359862B (zh) * | 2018-10-17 | 2020-09-18 | 北京师范大学 | 一种粮食作物实时估产方法及系统 |
-
2019
- 2019-07-10 CN CN201910618693.3A patent/CN110309985B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN110309985A (zh) | 2019-10-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110309985B (zh) | 一种农作物产量预测方法及系统 | |
| CN109829234B (zh) | 一种基于高分辨率遥感数据和作物模型的跨尺度高精度动态作物长势监测和估产方法 | |
| Ma et al. | Estimating regional winter wheat yield by assimilation of time series of HJ-1 CCD NDVI into WOFOST–ACRM model with Ensemble Kalman Filter | |
| Guerschman et al. | Land use impacts on the normalized difference vegetation index in temperate Argentina | |
| Zhuo et al. | Crop yield prediction using MODIS LAI, TIGGE weather forecasts and WOFOST model: A case study for winter wheat in Hebei, China during 2009–2013 | |
| CN113268923B (zh) | 一种基于模拟多光谱的夏玉米产量估算方法 | |
| Tardáguila et al. | Variations of soil properties affect the vegetative growth and yield components of “Tempranillo” grapevines | |
| CN109711102B (zh) | 一种作物灾害损失快速评估方法 | |
| CN109757175B (zh) | 一种基于无人机监测的玉米水肥一体化变量施肥方法 | |
| CN105184445A (zh) | 一种玉米干旱气象灾害多年平均损失率的计算方法 | |
| CN113673490B (zh) | 一种物候期自适应的作物生理参数遥感估测方法及系统 | |
| CN107860889A (zh) | 土壤有机质的预测方法和设备 | |
| CN110069032A (zh) | 一种基于小波神经网络的茄子温室环境智能检测系统 | |
| CN114818888B (zh) | 基于多通道卡尔曼滤波的土壤成分数据融合方法和系统 | |
| CN111044516A (zh) | 一种水稻叶绿素含量遥感估测方法 | |
| CN110569605A (zh) | 一种基于nsga2-elm的粳稻叶片氮素含量反演模型方法 | |
| CN115018105A (zh) | 一种冬小麦气象产量预测方法及系统 | |
| CN114062439B (zh) | 一种利用时间序列遥感影像联合估算土壤剖面盐分的方法 | |
| Bregaglio et al. | Improving crop yield prediction accuracy by embedding phenological heterogeneity into model parameter sets | |
| Liu et al. | Winter wheat yield estimation based on assimilated Sentinel-2 images with the CERES-Wheat model | |
| Yuan et al. | Research on rice leaf area index estimation based on fusion of texture and spectral information | |
| CN112945881B (zh) | 一种基于高光谱特征参数的马铃薯叶片含水量监测方法 | |
| Zou et al. | Combining spectral and texture feature of UAV image with plant height to improve LAI estimation of winter wheat at jointing stage | |
| CN112380497B (zh) | 用于区域玉米的秸秆系数估算方法及系统 | |
| Zhang et al. | Effects of direct heat stress on summer maize and risk assessment |
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 |