CN109800921B - 基于遥感物候同化和粒子群优化的区域冬小麦估产方法 - Google Patents
基于遥感物候同化和粒子群优化的区域冬小麦估产方法 Download PDFInfo
- Publication number
- CN109800921B CN109800921B CN201910094036.3A CN201910094036A CN109800921B CN 109800921 B CN109800921 B CN 109800921B CN 201910094036 A CN201910094036 A CN 201910094036A CN 109800921 B CN109800921 B CN 109800921B
- Authority
- CN
- China
- Prior art keywords
- period
- winter wheat
- data
- yield
- crop growth
- 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 66
- 239000002245 particle Substances 0.000 title claims abstract description 58
- 241000209140 Triticum Species 0.000 title claims abstract description 53
- 235000021307 Triticum Nutrition 0.000 title claims abstract description 53
- 238000005457 optimization Methods 0.000 title claims abstract description 37
- 230000012010 growth Effects 0.000 claims abstract description 93
- 230000008569 process Effects 0.000 claims abstract description 31
- 238000011161 development Methods 0.000 claims abstract description 27
- 238000011160 research Methods 0.000 claims abstract description 20
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 18
- 230000002194 synthesizing effect Effects 0.000 claims abstract description 5
- 238000004088 simulation Methods 0.000 claims description 34
- 239000002689 soil Substances 0.000 claims description 23
- 238000005070 sampling Methods 0.000 claims description 9
- 230000000630 rising effect Effects 0.000 claims description 8
- 230000009105 vegetative growth Effects 0.000 claims description 6
- 238000011156 evaluation Methods 0.000 claims description 5
- 238000000605 extraction Methods 0.000 claims description 5
- 238000001556 precipitation Methods 0.000 claims description 5
- 230000005855 radiation Effects 0.000 claims description 5
- 230000002159 abnormal effect Effects 0.000 claims description 4
- 230000015572 biosynthetic process Effects 0.000 claims description 4
- 230000035699 permeability Effects 0.000 claims description 4
- 230000029553 photosynthesis Effects 0.000 claims description 4
- 238000010672 photosynthesis Methods 0.000 claims description 4
- 230000029058 respiratory gaseous exchange Effects 0.000 claims description 4
- 230000001850 reproductive effect Effects 0.000 claims description 3
- 230000004807 localization Effects 0.000 abstract description 14
- 230000018109 developmental process Effects 0.000 description 19
- 235000013339 cereals Nutrition 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000003306 harvesting Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000013508 migration Methods 0.000 description 2
- 230000005012 migration Effects 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000002360 preparation method Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 208000005156 Dehydration Diseases 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000003967 crop rotation Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000009313 farming Methods 0.000 description 1
- 230000004720 fertilization Effects 0.000 description 1
- 238000003973 irrigation Methods 0.000 description 1
- 230000002262 irrigation Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000007261 regionalization Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 239000011163 secondary particle Substances 0.000 description 1
- 230000005068 transpiration Effects 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Cultivation Of Plants (AREA)
Abstract
本发明公开了一种基于遥感物候同化和粒子群优化算法的区域冬小麦估产方法，包括如下步骤：获取运行作物生长模型所需要的数据；获取研究区域全生育期内的遥感LAI数据，并按时间序列合成，对每个格点单元的LAI节点数据使用三次样条插值得到连续的LAI曲线；识别LAI曲线的第一拐点和第二拐点，将其分别作为冬小麦返青期和抽穗期的标志；将能够同时识别出返青期和抽穗期的格点作为冬小麦种植区域；使用二分法逐格点逐年优化作物生长发育参数；使用粒子群算法优化生物物理过程参数，以完成参数本地化；将本地化后的作物生长发育参数和生物物理过程参数带入作物生长模型运行，获得各个格点产量，汇总获得冬小麦种植区域的总产量。
Description
技术领域
本发明涉及农业遥感技术领域，更具体涉及一种基于遥感物候同化和粒子群优化算法的区域冬小麦估产方法。
背景技术
作物生长模型以光、温、水、土壤等条件为环境驱动变量，运用数学物理方法和计算机技术，对作物生育期内光合、呼吸、蒸腾等重要生理生态过程及其与气象、土壤等环境条件以及耕作、灌溉、施肥等技术条件的关系进行定量描述和预测，再现农作物生长发育及产量形成过程。由于作物生长模型机理性较强，包含参数较多，因此参数的可靠性是影响模型应用的关键问题。大多数模型是在田间尺度上建立和检验的，将其扩展到区域尺度乃至全球尺度存在着困难，不同的区域气候、土壤、作物和水文要素等不同，模型预测的准确性也不同，因此如何高效地实现作物生长模型参数本地化，获得可靠稳定的参数仍是当前研究的热点问题。
目前，利用数据同化技术把遥感反演参数融合到作物机理过程模型，是改进区域作物生长模拟精度和提高作物估产精度的重要途径。在遥感数据选择方面，随着研究范围的扩大，以 MODIS为代表的中分辨率遥感数据应用更为广泛。然而，MODIS 叶面积指数(LAI)在冬小麦区域存在的严重低估现象使得同化标准MODIS LAI后产量模拟严重偏低。同时，作物生长模型仅依靠自身模拟能力或同化关键状态变量对实际物候进行估算时，会不可避免地出现物候偏差，最终直接影响产量模拟结果。因此，如何充分挖掘现有成熟遥感产品(例如中分辨率的MODIS LAI及其衍生产品)显著提升作物生长模型的区域估产能力，满足区域农业估产的业务化需求仍是当前亟需解决的问题。
另外，由于作物生长模型的设计最初是基于田间站点尺度的，模型的参数本地化需要大量精细的田间作物生长参数和田间管理措施数据，因此当将模型推广到区域尺度时存在资料获取困难的问题。现有一些研究将田间、站点尺度的模型参数进行空间插值，由此实现参数区域化。但是插值参数受插值方法，插值区域气候、土壤、作物和水文等的影响存在很大的不确定性，常常导致区域产量模拟精度不高。
为了解决基于田间站点的作物生长模型本地化参数难以适用于区域尺度的问题，遥感数据同化作物生长模型技术被广泛应用到区域作物估产中。现有的很多遥感数据同化方案中，为了尽可能避免因遥感数据误差所带来的同化精度影响，对作物生长参数进行高密度的实地采样观测从而确保高分辨率遥感数据向状态变量(例如LAI)的反演精度，辐射模型的校验精度或实现对现有遥感产品的二次修正。然而，当研究范围进一步扩大时，大范围内高密度、高频率的地表观测需求将在时间、人力、金钱成本方面都面临巨大的压力，区域农业估产难以实现业务化应用。
因此，需要新的技术以至少部分解决现有技术中存在的局限。
发明内容
为了解决上述技术问题，本发明提供了一种基于遥感物候信息，有效实现作物生长模型参数本地化的区域作物估产方法。
根据本发明的一方面，提供一种基于遥感物候同化和粒子群优化算法的区域冬小麦估产方法，包括如下步骤：
S1:获取运行作物生长模型所需要的数据，所述作物生长模型为MCWLA-Wheat，所述数据包括气象数据和土壤属性数据；
S2:获取研究区域全生育期内的遥感LAI数据，并按时间序列合成，对每个格点单元的LAI节点数据使用三次样条插值得到连续的LAI曲线；
S3:识别S2步骤所得到的LAI曲线的第一拐点和第二拐点 (峰值)，将其分别作为冬小麦返青期和抽穗期的标志，各个格点单元的成熟期由农气站的记录进行反距离加权插值获得；
S4:将S3步骤中能够同时识别出返青期和抽穗期的格点作为冬小麦种植区域，用于实施后续的步骤；
S5：基于S3和S4步骤获得的冬小麦种植格点的物候信息，使用二分法逐格点逐年优化作物生长发育参数，所述生物发育参数包括营养生长期前期生长速率参数、营养生长期后期生长速率参数和生殖生长期生长速率参数；
S6:在S5步骤获得优化作物生长发育参数之后，再使用粒子群算法优化生物物理过程参数，以完成参数本地化，其中所述生物物理过程参数包括与土壤水分平衡、光合、呼吸作用以及产量形成模块相关的参数，
其中粒子群算法优化包括：
(1)根据生物物理过程参数的先验取值区间进行随机抽样，使用拉丁超立方体抽样法(Latin hypercube sampling)抽取N套初始化参数，称为一个初始输入粒子群；重复进行M次粒子群抽取，共形成M个初始输入粒子群共M*N套初始化参数，其中，M和N为大于1的自然数；
(2)将其中一个初始粒子群带入作物生长模型进行运行，调用粒子群优化算法模块，设定适宜度为模拟产量和实际记录产量之间的误差，使用的适宜度评价指标是标准均方根误差, 公式如下：
式中，RMSE为均方根误差，NRMSE为标准均方根误差，n是数据个数，S和O分别表示模拟和观测产量，是观测产量平均值；
当满足以下迭代终止条件之一时，则结束当前粒子群优化：
a)全局最优适宜度小于10％；
b)全局最优适宜度在连续5次循环后无法提高0.0001％；
c)粒子群连续更新超过100次。
(3)重复过程(2)直到m个初始粒子群均到达迭代终点，然后将m个粒子群的整体最优解作为全局最优参数组，完成参数本地化工作；
S7：将本地化后的作物生长发育参数和生物物理过程参数带入作物生长模型运行，获得各个格点产量，汇总获得冬小麦种植区域的总产量。
根据本发明的一个实施方案，步骤S1中所述作物生长模型为MCWLA-Wheat。
根据本发明的一个实施方案，其中步骤S1中所述气象数据包括每日最高气温、最低气温、太阳辐射、相对湿度、降水和风速，所述土壤属性数据包括基于土壤质地的渗透率和有效田间持水量。
根据本发明的一个实施方案，其中步骤S2中，所述遥感LAI 数据来自北京师范大学反演制作的GLASS LAI(Global LAnd Surface Satellites(GLASS)Leaf Area Index(LAI)) 产品。
根据本发明的一个实施方案，其中在S3步骤中，根据LAI 曲线识别返青期和抽穗期包括：对研究区域冬小麦返青期和和抽穗期的可能时期进行划定以避免将其他地物误识别为冬小麦，在划定的可能时期内进行识别，其中对返青期的识别满足如下条件：
(1)返青期节点(第一拐点)之前的遥感反演数据维持在较低值的水平，不存在异常的较高值；
(2)返青期节点(第一拐点)之后的遥感反演数据需要识别出连续上升趋势；
其中对抽穗期的识别满足如下条件：
(1)峰值节点(第二拐点)前至少有连续三期的遥感数据反映出上升趋势；
(2)峰值节点(第二拐点)后至少有连续三期的遥感数据呈现下降趋势。
根据本发明的一个实施方案，其中步骤S5中，所述作物生长发育参数来自MCWLA-Wheat模型。
与现有技术相比，本发明具有诸多的优点：
1.采用的作物生长模型是为满足区域尺度作物生长模拟需求所开发的MCWLA(Model to capture the Crop-Weather relationship over a Large Area)模型,该模型基于大区域尺度进行开发，数据准备资料获取容易。
2.使用的遥感产品包括GLASS LAI(Global LAnd Surface Satellites(GLASS)Leaf Area Index(LAI))产品,该产品是对MODIS LAI产品的改进，在生产过程中已经包括了S-G滤波预处理，在数据质量上有一定提升且具有更好的时间完整性和空间连续性，节省大量遥感数据预处理的时间成本。
3.本发明中基于遥感数据的物候反演方法可以快速准确地得到大区域的精细物候信息。
4.本发明充分发挥了遥感数据和作物生长模型的优势，把遥感观测的LAI时序曲线所存在的物候信息有效利用到作物生长模型参数本地化工作中，避免了直接使用绝对LAI数据存在产量模拟偏低的问题，同时也解决了作物生长模型在模拟过程中存在的物候偏移问题。
5.本发明在县级尺度上的粒子群参数优化方案可以快速高效地获得全局最优解，最终基于最优参数组模拟的产量精度远远高于当前同类发明。
附图说明
附图中相同的附图标记标示了相同或类似的部件或部分。本发明的目标及特征考虑到如下结合附图的描述将更加明显，附图中：
图1是根据本发明一个实施方案的基于遥感物候同化和粒子群优化算法的区域冬小麦估产方法的流程示意图。
图2是根据本发明一个实施方案的应用本发明方法模拟区域冬小麦产量的结果图。
具体实施方式
为清楚的说明本发明中的方案，下面给出优选的实施例并结合附图详细说明。以下的说明本质上仅仅是示例性的而并不是为了限制本公开的应用或用途。
应该理解的是，本发明所引用的作物生长模型和遥感模型本身是已知的，例如模型的各个子模块、各种参数、运行机制等等，因此本发明重点阐述优化参数的过程。
图1是根据本发明一个实施方案的基于遥感物候同化和粒子群优化算法的区域冬小麦估产方法的流程示意图。如图1所示，本发明提供了一种基于遥感物候信息，有效实现作物生长模型参数本地化的区域作物估产方法，具体步骤如下：
S1:获取运行作物生长模型所需要的数据；例如其中所采用的作物生长模型可以为MCWLA(Model to capture the Crop-Weather relationship over a Large Area)系列模型中的MCWLA-Wheat。当然也可以采用其他合适的作物生长模型。 MCWLA模型是一系列为满足区域尺度作物生长模拟需求所开发的作物生长模型,已经应用在了国内外多个作物生长模拟研究中，表现稳健。S1步骤所需的数据可以包括气象数据和土壤属性数据。气象数据可以包括每日最高气温、最低气温、太阳辐射、相对湿度、降水和风速。土壤属性数据可以包括基于土壤质地的渗透率和有效田间持水量。这些数据可以通过现场监测或者现有的记录数据来获取。
S2:对研究区域全生育期内的遥感LAI数据按时间序列合成，对每个格点单元的LAI节点数据使用三次样条插值得到连续的LAI曲线。使用的遥感数据可以是北京师范大学反演制作的 GLASS LAI(Global LAnd Surface Satellites(GLASS)Leaf Area Index(LAI))产品,，该产品可从北京师范大学全球变化数据处理与分析中心(http://www.bnu- datacenter.com/) 和马里兰大学全球土地利用数据库 (http://glcf.umd.edu/data/ helpme.shtml)获得。该产品是对MODIS LAI产品的改进，其在生产过程中已经包括了S-G滤波预处理，在数据质量上有一定提升且具有更好的时间完整性和空间连续性。应该理解，也可使用其他合适的遥感LAI数据。
S3:识别S2得到的LAI曲线的第一拐点和第二拐点，第二拐点即为峰值，将其分别作为冬小麦返青期和抽穗期的标志，各个格点单元的成熟期由农气站的记录进行反距离加权插值获得。
其中，S3根据LAI曲线识别返青期和抽穗期具体为：首先，对研究区域冬小麦返青期和和抽穗期的可能时期进行划定以避免将其他地物误识别为冬小麦。然后，在划定的可能时期内进行识别，其中对返青期的识别需要满足如下条件：
(1)返青期节点之前的遥感反演数据维持在较低值的水平，不存在异常的较高值；可以通过与LAI最高值相比较来判断所谓的较高值和较低值。例如可以设定大约在LAI最高值的15％以下是判断为较低值的水平，超过最高LAI最高值的50％可以判断为较高值的水平；
(2)返青期节点之后的遥感反演数据需要识别出连续上升趋势；
对抽穗期的识别需要满足如下条件：
峰值节点前至少有连续三期的遥感数据反映出上升趋势；
峰值节点后至少有连续三期的遥感数据呈现下降趋势。
影响成熟期的遥感反演的因素较多，例如天气状况，农民的种植收获习惯，格点内可能存在的其他植被的生长以及背景噪音的影响等等都会影响对遥感反演的LAI曲线末端的识别。因此，为了提高效率，各个格点单元的成熟期由农气站的记录进行反距离加权插值获得。反距离加权插值以样本点与待插值点间的距离为权重，待插值点与样本点越近，则该点权重越大，反之越小。反距离加权插值法为本领域普通技术人员所述熟知。例如其可以采用如下一般公式：
其中，为S0处的待插值点值；N为已知点数目；λi为在第i预测点位置上的权重，样本点与待插值点之间的距离越大，该值越小；Z(Si)是在Si处的测量值。
权重的确定计算公式为：
其中，di0是待插值点S0与样本点Si之间的距离；p是指数，默认值是2。
S4:将S3步骤中能够同时识别出返青期和抽穗期的格点作为冬小麦种植区域，用于后续的研究步骤。
S5：基于S3和S4步骤获得的冬小麦种植格点的物候信息，使用二分法逐格点逐年优化作物生长发育参数。
其中，S5步骤中所述的作物生长发育参数可以为 MCWLA-Wheat模型中的营养生长期前期生长速率参数、营养生长期后期生长速率参数和生殖生长期生长速率参数。针对作物生长模型中作物生长发育的模拟方程是一个时间上的单调函数，且每个发育阶段仅设置了一个待优化参数的特点，S5步骤选择了二分法实现对作物生长发育参数的优化。
S6:在S5获得最优作物生长发育参数之后再使用粒子群算法优化其他生物物理过程参数，以完成参数本地化过程。
其中，S6步骤中所述的生物物理过程参数包括涉及土壤水分平衡、光合、呼吸作用、产量形成等模块的参数，例如产量差异系数、土壤水胁迫因子对作物生长产生影响的阈值等。由于生物物理过程涉及参数较多，因此采用粒子群优化算法进行参数优化，具体如下：
(1)根据参数的先验取值区间进行随机抽样，使用拉丁超立方体抽样法(Latinhypercube sampling)抽取N套(例如10 套，当然也可以根据需要选择其他的数量)初始化参数，称为一个初始输入粒子群。为了有效避免粒子群算法在处理复杂多参数优化问题时可能陷入局部最优的情况，本研究重复进行了M 次(例如5次，当然也可以根据需要选择其他的数量)粒子群抽取，共形成M个初始输入粒子群共M×N套初始化参数，其中，M 和N为大于1的自然数。
(2)将其中一个初始粒子群带入作物生长模型例如 MCWLA-Wheat模型进行运行，调用粒子群优化算法模块，设定适宜度为模拟产量和实际记录产量之间的误差，使用的适宜度评价指标是标准均方根误差,公式如下：
式中，RMSE为均方根误差，NRMSE为标准均方根误差，n是数据个数，S和O分别表示模拟和观测产量，是观测产量平均值。
当满足以下迭代终止条件之一时，则结束当前粒子群优化：
a)全局最优适宜度(即最小NRMSE)小于10％；
b)全局最优适宜度在连续5次循环后无法提高0.0001％；
c)粒子群连续更新超过100次。
(3)重复过程(2)直到5个初始粒子群均到达迭代终点，然后将M个粒子群的整体最优解作为全局最优参数组，完成参数本地化工作。
S7：将本地化后的作物生长发育参数和生物物理过程参数带入作物生长模型运行，获得各个格点产量，汇总获得冬小麦种植区域的总产量。
逐格点运行作物生长模型模拟产量，然后汇总到县域行政单元，基于县级冬小麦单产记录，以县级模拟产量的标准均方根误差最小作为最优适宜度判断标准，获得每个县的全局最优参数组，完成模型本地化校准工作。用本地化参数运行模型，输出县级单产模拟结果。
实施例
下面以河北、山东、河南和安徽四省的冬小麦种植区域作为研究区域，示例性说明本发明的方法的具体应用。
步骤S1，选择河北、山东、河南和安徽四省的冬小麦种植区域作为研究区域，共涉及159个县。研究区域地处 111.5°E–118.2°E，32.1°N–38.8°N，地形以平原为主，占据了大部分华北平原的冬小麦主产区。气候属温带季风性气候，平均温度约15.6℃，年均降水量约834.5mm。获取以下数据：气象数据采用了拥有足够插值精度(R2>0.8)，由Yuan等人[1](Yuan W,Xu B,Chen Z,Xia J,Xu W,Chen Y,Wu X,Fu Y. (2015).Validation of China-wideinterpolated daily climate variables from 1960to 2011.Theor Appl Climatol119:689-700.)制作的全国0.5°×0.5°格网气象数据集中 2000-2008年的部分，该数据集包含了MCWLA-Wheat模型模拟所需的天步长气象数据，包括最高气温(℃)、最低气温(℃)、太阳辐射(MJ/m2)、相对湿度(％)、日总降水量(mm)、风速(m/s) 等数据。模型模拟中的土壤属性数据来源于FAO土壤数据库，包括基于土壤质地的渗透率和有效田间持水量，分辨率为5’。土壤属性数据库主要应用于土壤水平衡部分的计算。2001-2008 年县级冬小麦单产记录从中国种植业信息网和县级统计年鉴上整理得出，主要用于作物生长模型的参数本地化与模拟效果评估。农业气象站的物候记录数据从中国气象数据网获得。遥感数据使用的是2001-2008年的GLASS LAI产品，该产品的空间分辨率为1km，时间分辨率为8天，满足本发明的需求。
步骤S2，对研究区域全生育期内的GLASS LAI数据按时间序列合成，对每个格点单元的LAI节点数据使用三次样条插值得到连续的LAI曲线。
步骤S3,识别S2得到的LAI曲线的第一拐点和第二拐点，将其分别作为冬小麦返青期和抽穗期的标志。首先，对研究区域冬小麦返青期和和抽穗期的可能时期进行划定以避免将其他地物误识别为冬小麦。其中，根据研究区内的气候特征，返青期的时间设定在每年的2-3月，抽穗期设定在4月左右。然后，在划定的可能时期内进行识别，其中对返青期的识别需要满足如下条件：
(1)返青期节点之前的遥感反演数据维持在较低值的水平，本案例中选择大约在LAI最高值的15％以下，不存在异常的高值，选择超过LAI最高值的50％；
(2)返青期节点之后的遥感反演数据需要识别出连续上升趋势；
对抽穗期的识别需要满足如下条件：
(1)峰值节点前至少有连续三期的遥感数据反映出上升趋势；
(2)峰值节点后至少有连续三期的遥感数据呈现下降趋势；
影响成熟期的遥感反演因素较多，例如天气状况，农民的种植收获习惯，格点内可能存在的其他植被的生长以及背景噪音的影响等等都会影响对遥感反演的LAI曲线末端的识别。因此，为了提高效率，各个格点单元的成熟期由农气站的记录进行反距离加权插值获得。
步骤S4，基于步骤S3识别返青期和抽穗期的原则，对初步划定的约65万个1km×1km格点2001-2008年的GLASS LAI数据进行处理，如果一个格点在某一年的LAI反演数据完全符合上述提取规则，则认为当地在当年种植了冬小麦，反之则认为该格点在当前年份没有种植冬小麦。最终，保留在2001-2008 年期间至少有一年被识别为种植了冬小麦的格点，作为研究区域内冬小麦的种植区以及潜在的研究区。考虑大量格点的计算需要花费较高的时间成本，因此将格点数据分辨率重采样到8km 以提高计算效率，也为大区域业务化应用做准备。
步骤S5，针对MCWLA-Wheat模型中作物生长发育的模拟方程是一个时间上的单调函数，且每个发育阶段仅设置了一个待优化参数的特点，选择了二分法逐格点逐年实现对作物生长发育参数的优化。对返青期和抽穗期相关参数的校准是基于格点尺度以及年步长进行的，每个格点在每一年都会对相关的作物生长发育参数进行重新校准以获得与遥感反演数据相匹配的物候模拟结果。对于成熟期相关参数的校准采用了参数先验的办法，即在2001-2008年对每个格点的成熟期进行参数校准，在每个格点获得一个使模拟和农气站插值成熟期年际误差最小的成熟期作物生长发育参数，然后将该参数应用于当前格点 2001-2008年的模型模拟。
步骤S6，由于生物物理过程涉及参数较多，因此采用粒子群优化算法以县为单位进行参数优化，即同一县内格点共用一套生物物理过程参数。逐年运行作物生长模型模拟各格点产量，然后汇总到县域行政单元，基于县级冬小麦单产记录，以 2001-2008年模型模拟的县级产量的标准均方根误差最小作为最优适宜度判断标准，当满足迭代终止条件时，当前粒子群优化结束。粒子群参数优化方案具体如下：
(1)根据参数的先验取值区间进行随机抽样，使用拉丁超立方体抽样法(Latinhypercube sampling)抽取10套初始化参数，称为一个初始输入粒子群。为了有效避免粒子群算法在处理复杂多参数优化问题时可能陷入局部最优的情况，本发明重复进行了5次粒子群抽取，共形成5个初始输入粒子群共50套初始化参数。
(2)将一个初始粒子群带入MCWLA-Wheat模型进行运行，调用粒子群优化算法模块，设定适宜度为模拟产量和实际记录产量之间的误差，使用的适宜度评价指标是标准均方根误差, 公式如下：
式中n是总年份，S和O分别表示模拟和观测产量，是观测产量平均值。
当满足以下迭代终止条件之一时，则结束当前粒子群优化：
a)全局最优适宜度小于10％；
b)全局最优适宜度在连续5次循环后无法提高0.0001％；
c)粒子群连续更新超过100次。
(3)重复过程(2)直到5个初始粒子群均到达迭代终点，然后将5个粒子群的整体最优解作为全局最优参数组，完成参数本地化。
将本地化后的作物生长发育参数和生物物理过程参数带入MCWLA-Wheat模型运行，获得各个格点产量，再结合各格点的冬小麦种植比例，按县域行政边界汇总，获得区域单产结果，图2为模型模拟的2001-2008年159个县的冬小麦单产情况图。
本发明实施例所述的一种基于遥感物候同化和粒子群优化算法的区域冬小麦估产方法，充分发挥了遥感数据和作物生长模型的优势，把遥感观测的LAI时序曲线所存在的物候信息充分利用到作物生长模型参数本地化工作中，避免了直接使用绝对LAI数据存在产量模拟偏低的问题，同时也解决了作物生长模型在模拟过程中存在的物候偏移问题。本发明针对不同县的气候、土壤、水文条件以及管理措施的差异，利用粒子群优化算法为每个县因地制宜地设置参数组，极大地提高了产量模拟精度。其中，2001-2008年159个县的冬小麦产量模拟均方根误差 RMSE为465kg ha-1，决定系数R2为0.75，是目前国内同类型发明中适用研究区域最大，时间跨度最长，模拟精度最高的发明。
粮食的生产者与消费者都需要及时准确地了解粮食产量信息，根据本方法可以花费最少的时间和人力在冬小麦成熟期大面积地获得产量数据，为国家有关部门进行粮情判断、粮食调控等科学决策等提供重要的科学依据，并且可以作为粮食贸易的重要依据。
本发明充分利用了现有的成熟遥感产品在格点和县级尺度上实现作物生长模型参数本地化，节省了大量获取地面精细实测数据的时间、人力、金钱成本，避免了直接使用绝对LAI 数据存在产量模拟偏低的问题，同时也解决了作物生长模型在模拟过程中存在的物候偏移问题。大大提高了区域农业估产的精度和效率，为区域快速农业估产等业务化应用提供可能。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的装置及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。
Claims (5)
1.一种基于遥感物候同化和粒子群优化算法的区域冬小麦估产方法，包括如下步骤：
S1:获取运行作物生长模型所需要的数据，所述数据包括气象数据和土壤属性数据；
S2:获取研究区域全生育期内的遥感LAI数据，并按时间序列合成，对每个格点单元的LAI节点数据使用三次样条插值得到连续的LAI曲线；
S3:识别S2步骤所得到的LAI曲线的第一拐点和第二拐点，将其分别作为冬小麦返青期和抽穗期的标志，各个格点单元的成熟期由农气站的记录进行反距离加权插值获得；
S4:将S3步骤中能够同时识别出返青期和抽穗期的格点作为冬小麦种植区域，用于实施后续的步骤；
S5：基于S3和S4步骤获得的冬小麦种植格点的物候信息，使用二分法逐格点逐年优化作物生长发育参数，所述作物生长发育参数包括营养生长期前期生长速率参数、营养生长期后期生长速率参数和生殖生长期生长速率参数；
S6:在S5步骤获得优化作物生长发育参数之后，再使用粒子群算法优化生物物理过程参数，其中所述生物物理过程参数包括与土壤水分平衡、光合、呼吸作用以及产量形成模块相关的参数，
其中粒子群算法优化包括：
(1)根据生物物理过程参数的先验取值区间进行随机抽样，使用拉丁超立方体抽样法(Latin hypercube sampling)抽取N套初始化参数，称为一个初始输入粒子群；重复进行M次粒子群抽取，共形成M个初始输入粒子群共M*N套初始化参数，其中，M和N为大于1的自然数；
(2)将其中一个初始粒子群带入作物生长模型进行运行，调用粒子群优化算法模块，设定适宜度为模拟产量和实际记录产量之间的误差，使用的适宜度评价指标是标准均方根误差,公式如下：
式中，RMSE为均方根误差，NRMSE为标准均方根误差，n是数据个数，S和O分别表示模拟和观测产量，是观测产量平均值；
当满足以下迭代终止条件之一时，则结束当前粒子群优化：
a)全局最优适宜度小于10％；
b)全局最优适宜度在连续5次循环后无法提高0.0001％；
c)粒子群连续更新超过100次，
(3)重复过程(2)直到m个初始粒子群均到达迭代终点，然后将m个粒子群的整体最优解作为全局最优参数组，完成作物生长发育参数和生物物理过程参数本地化工作；
S7：将本地化后的作物生长发育参数和生物物理过程参数带入作物生长模型运行，获得各个格点产量，汇总获得冬小麦种植区域的总产量。
2.根据权利要求1所述的区域冬小麦估产方法，其特征在于，步骤S1中所述作物生长模型为MCWLA-Wheat。
3.根据权利要求1所述的区域冬小麦估产方法，其特征在于，步骤S1中所述气象数据包括每日最高气温、最低气温、太阳辐射、相对湿度、降水和风速，所述土壤属性数据包括基于土壤质地的渗透率和有效田间持水量。
4.根据权利要求1所述的区域冬小麦估产方法，其特征在于，在S3步骤中，根据LAI曲线识别返青期和抽穗期包括：对研究区域冬小麦返青期和和抽穗期的可能时期进行划定以避免将其他地物误识别为冬小麦，在划定的可能时期内进行识别，其中对返青期的识别满足如下条件：
(1)返青期拐点之前的遥感反演数据维持在较低值的水平，不存在异常的较高值；
(2)返青期拐点之后的遥感反演数据需要识别出连续上升趋势；
其中对抽穗期的识别满足如下条件：
(1)抽穗期拐点前至少有连续三期的遥感数据反映出上升趋势；
(2)抽穗期拐点后至少有连续三期的遥感数据呈现下降趋势。
5.根据权利要求1所述的区域冬小麦估产方法，其特征在于，步骤S5中，所述作物生长发育参数来自MCWLA-Wheat模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910094036.3A CN109800921B (zh) | 2019-01-30 | 2019-01-30 | 基于遥感物候同化和粒子群优化的区域冬小麦估产方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910094036.3A CN109800921B (zh) | 2019-01-30 | 2019-01-30 | 基于遥感物候同化和粒子群优化的区域冬小麦估产方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109800921A CN109800921A (zh) | 2019-05-24 |
| CN109800921B true CN109800921B (zh) | 2019-12-20 |
Family
ID=66560564
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910094036.3A Active CN109800921B (zh) | 2019-01-30 | 2019-01-30 | 基于遥感物候同化和粒子群优化的区域冬小麦估产方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109800921B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110426755B (zh) * | 2019-07-26 | 2021-09-14 | 南京信息工程大学 | 一种基于改进反距离加权的地面气温资料质量控制方法 |
| CN110633841B (zh) * | 2019-08-13 | 2022-04-01 | 中国农业大学 | 基于集合采样的省级范围地块尺度数据同化产量预测方法 |
| CN110766308B (zh) * | 2019-10-17 | 2020-07-10 | 中国科学院地理科学与资源研究所 | 一种基于集合同化策略的区域农作物估产方法 |
| CN112579981B (zh) * | 2020-12-24 | 2022-05-31 | 青岛大学 | 一种耦合灌溉的区域尺度冬小麦产量遥感模拟方法 |
| CN113077077B (zh) * | 2021-03-18 | 2022-05-27 | 四川农业大学 | 带状复合作物种植产量评估方法、装置、设备及存储介质 |
| CN116756691B (zh) * | 2023-06-25 | 2024-01-30 | 国家海洋环境预报中心 | 一种海洋数据同化方法、系统、电子设备及介质 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8346711B2 (en) * | 2009-11-24 | 2013-01-01 | King Fahd University Of Petroleum And Minerals | Method for identifying multi-input multi-output Hammerstein models |
| CN102651096B (zh) * | 2012-04-28 | 2015-06-03 | 中国农业大学 | 同化叶面积指数时序曲线特征的冬小麦估产方法 |
| CN104134095B (zh) * | 2014-04-17 | 2017-02-15 | 中国农业大学 | 一种基于尺度转换和数据同化的农作物产量估测方法 |
| CN109146111A (zh) * | 2017-06-27 | 2019-01-04 | 中国农业大学 | 一种基于arima-lssvm组合模型预测粮食产量的方法 |
| CN108509836B (zh) * | 2018-01-29 | 2021-10-08 | 中国农业大学 | 双极化合成孔径雷达与作物模型数据同化的作物估产方法 |
| CN108984803B (zh) * | 2018-10-22 | 2020-09-22 | 北京师范大学 | 一种农作物产量空间化的方法及系统 |
-
2019
- 2019-01-30 CN CN201910094036.3A patent/CN109800921B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109800921A (zh) | 2019-05-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109800921B (zh) | 基于遥感物候同化和粒子群优化的区域冬小麦估产方法 | |
| CN110751094B (zh) | 一种基于gee综合遥感影像和深度学习方法的作物估产方法 | |
| Huang et al. | The improved winter wheat yield estimation by assimilating GLASS LAI into a crop growth model with the proposed Bayesian posterior-based ensemble Kalman filter | |
| Gaál et al. | Modelling the impact of climate change on the Hungarian wine regions using random forest | |
| CN110909933B (zh) | 一种耦合作物模型与机器学习语言的农业干旱快速诊断和评估方法 | |
| CN109711102B (zh) | 一种作物灾害损失快速评估方法 | |
| Wang et al. | Summer maize growth under different precipitation years in the Huang-Huai-Hai Plain of China | |
| Sirsat et al. | Machine Learning predictive model of grapevine yield based on agroclimatic patterns | |
| CN105912836A (zh) | 一种纯遥感数据驱动的流域水循环模拟方法 | |
| CN110705182B (zh) | 耦合作物模型和机器学习的作物育种适应时间预测方法 | |
| US12112105B1 (en) | Soil-climate intelligent type determining method for rice target yield and nitrogen fertilizer amount | |
| CN118072178B (zh) | 基于分类后百分比数据同化的玉米产量估算方法及系统 | |
| CN114186423A (zh) | 一种雪茄烟品的适宜种植区预测评估方法及系统 | |
| CN113011372A (zh) | 一种盐碱化土地自动监测和识别方法 | |
| CN110533287A (zh) | 夏玉米干旱灾损定量测评模型的构建与应用 | |
| CN109272161A (zh) | 基于动态hi的水稻估产方法 | |
| Borus et al. | Improving the prediction of potato productivity: APSIM-Potato model parameterization and evaluation in Tasmania, Australia | |
| CN116861298A (zh) | 一种无资料地区的流域水文模型参数估计方法 | |
| CN117763970A (zh) | 基于图像深度学习的稀缺资料地区水文模型参数重建方法 | |
| Cheng et al. | Improving soil available nutrient estimation by integrating modified WOFOST model and time-series earth observations | |
| Wang et al. | Exploring the feasibility of winter wheat freeze injury by integrating grey system model with RS and GIS | |
| CN117010546A (zh) | 一种云南省次季节尺度温度异常的预测方法和装置 | |
| CN113762768B (zh) | 基于天气发生器和作物模型的农业旱灾动态风险评估方法 | |
| Sehgal | Remote sensing for crop growth and crop simulation modelling | |
| CN117744898B (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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |