CN109472320B - 不确定条件下作物生育期模型品种参数自动校正方法 - Google Patents

Abstract

本发明公开了一种不确定条件下作物生育期模型品种参数自动校正方法，针对在不确定条件下，相似性度量算法以及自适应差分进化算法与作物生育期模型耦合过程中一些关键环节。该方法解决选取多年份温光属性特征差异度大的数据集代表当地的气候条件，作为模型品种参数校正的样本数据集，即筛选出能够约束所有参数的关键物候期变量作为适应度函数中目标变量，当变量实测数据缺失时利用栽培知识推算其缺失值；确定具有自适应特征差分的进化算法来校正模型品种参数；最后对“异参同效”的估计参数结果，采用改进K‑means聚类算法选择出一组鲁棒性较好的代表性品种参数，使得模拟误差最小并具有较好稳定性，便于用于实际生产实践。

Description

不确定条件下作物生育期模型品种参数自动校正方法

技术领域

本发明属于数字作物技术领域，是序列数据相似性分析算法、自适应差分进化算法与作物阶段发育模型的交叉领域。提出了一种可用于多参数、非线性、不连续、多峰值的作物生育期模型品种参数的校正样本集自动选择以及参数自动校正的综合方法。

背景技术

物候期预测为作物生长模型预测提供统一时间尺度的衡量标准，是分析气候变化对作物生长发育时间影响的重要工具。生育期模型模拟结果受到基因和温光环境等多因素的影响，模型的本地化应用存在着不确定性，主要表现在数据不确定性、模型结构不确定性以及模型参数不确定性三个方面。模型参数不确定性是指受到不同自然环境条件、不同栽培管理方式以及作物品种多样性的影响，参数值无法直接测量且有多种可能的组合。虽然在以往的研究实践中，作物建模专家明确了模型参数的取值范围，但是在实际应用中，还需要结合先验知识进一步确定出物候期模型品种参数的最优解。这项工作费时费力，人工操作往往会忽略最优解。因此，在不确定条件下对模型参数实现自动校准是提高作物模型在本地应用适应性与精度的重要技术手段。

作物模型参数校准主要有基于统计推断理论的方法和基于最优化理论的方法。基于统计推断理论的方法主流是贝叶斯统计估算，如马尔科夫蒙特卡罗方法(Markov ChainMonte Carlo，MCMC)、普适似然不确定估计法(General Likelihood UncertaintyEstimation， GLUE)等，已被广泛应用于WOFOST、RZWQM2、STICS等模型的参数优化。该类方法的优点是可以对模拟结果进行不确定性分析和量化，给出模型参数的后验概率分布，缺少最优解的直接推荐；基于最优化理论是指通过优化算法在模型参数空间分布上实现自动搜索，并按照一定的寻优规则识别较优或最优的参数集合。常用方法包括最小二乘法等传统优化方法以及遗传算法(Genetic Algorithm，GA)、粒子群优化算法(Particle SwarmOptimization， PSO)，差分进化算法(Differential Evolution，DE)等进化算法。以改进的非线性最小二乘法为基础的PEST(Parameter Estimation)参数自动优化工具，在作物模型和水文模型参数优化研究中得到了广泛的应用。但在参数太多或模型结构太复杂时往往无法收敛到最优值。作物生育期模型具有非线性、非凸以及参数之间相关等特点，模型参数优化问题可以转化为有约束的多模态函数优化问题，利用全局优化方法如进化算法(GA、DE、PSO)可以提供近似全局最优解。Soundharajan和Sudheer等利用标准GA对ORYZA2000模型参数实现了自动校准，校准后的参数能够适用正常灌溉和水分胁迫下的模拟。庄嘉祥等提出了一种改进型协同进化GA算法，个体优势遗传算法(Individual Advantages GeneticAlgorithm， IAGA)，突出了精英个体的引导能力，应用到RiceGrow和ORYZA2000物候期模型的参数校准中，提高了模型应用的精确程度。

等比较了差分进化算法(DE)、协方差矩阵适应进化策略(Covariance Matrix Adaptation Evolution Strategy，CMA-ES)，粒子群优化 (PSO)和人工蜂群算法(Artificial Bee Colony，ABC)应用于SUCROS生长模型的参数估计效果，发现应用DE算法的结果误差最小。综上，进化算法作为一类智能的全局最优化算法，具有自组织、自学习等特性，可以在一定迭代次数内找到最优解或近似最优解，在作物模型参数校准领域已经有了广泛应用。作物模型的结构包括多组常微分方程，具有非线性、不连续、多峰值特点，进化算法的应用效果通常取决于对特定类型问题的适应能力，因此针对作物模型特点选择合适的进化算法尤为重要。与此同时，进化算法自身参数对算法性能有较大影响，这些参数选择往往靠经验得出，增加了优化算法在模型参数校准过程中的不确定性。

此外，已有研究发现，利用不同的校正集样本数据来校正数据驱动的模型参数后，模型的模拟结果不同。例如：梁秀英等在《基于机器视觉的玉米性状参数与近红外光谱的玉米组分含量检测方法研究》(2013年)中用不同的数据集选择方法来划分玉米籽粒样本集，结果证明不同样本划分方法对PLS模型的预测结果有影响。庄嘉祥等在《基于个体优势遗传算法的水稻生育期模型参数优化》(2013年)中对RiceGrow和ORYZA2000生育期模型参数估算的研究中发现，在已有的不同年份的实测数据中，当选择隔年的数据来校正模型参数时，相比于其他校正数据集选择方案，模型能够取得较好的模拟结果；因此，选择具有代表性的校正集样本数据亦是提高模型预测性能的关键技术之一。

综上所述，在不确定环境下进化算法在应用到作物模型参数校正过程中目前普遍存在的问题是：(1)如何选择校正数据集，充分挖掘已有数据信息来一定程度上避免数据不确定性。(2)如何选择目标变量，使得能够约束所有的待调参数。进化算法在随机搜索的过程中，仅注重目标变量的拟合情况，而参数校正时数据集不完备，仅选取部分物候期作为目标变量，生育期模型各个物候期受各因子约束的程度不同，所选择的目标物候期无法很好的约束所有的待调参数。(3)进化算法参数校正的能力通常取决于特定类型的问题，不同进化算法解决问题的有效性也有具有不确定性。RiceGrow生育期模型表现出不连续、多峰值、非线性、多极值等复杂特性，在采用一些进化算法进行参数校正时，极易陷入局部最优解。因此，针对具体的模型选择合适的算法尤为重要。(4)模型校正结果会出现“异参同效”问题，即模型参数具有不确定性。如何从多组解集中选择一组鲁棒性好的代表性参数对作物模型应用于实际生产中具有重要意义。

发明内容

本发明针对在不确定条件下，相似性度量算法以及自适应差分进化算法与作物生育期模型耦合过程中一些关键环节，包括如何选择校正样本的数据集，如何确定进化算法适应度函数中的目标变量，针对模型特点选择合适的进化算法，针对“异参同效”如何选择一组有代表性的参数四个问题提出作物生育期模型品种参数自动校正方法，该方法解决选取多年份温光属性特征差异度大的数据集代表当地的气候条件，作为模型品种参数校正的样本数据集，即筛选出能够约束所有参数的关键物候期变量作为适应度函数中目标变量，当变量实测数据缺失时利用栽培知识推算其缺失值；确定具有自适应特征差分的进化算法来校正模型品种参数；最后对“异参同效”的估计参数结果，采用改进K-means聚类算法选择出一组鲁棒性较好的代表性品种参数，使得模拟误差最小并具有较好稳定性，便于用于实际生产实践。

技术方案：

本发明公开了一种不确定条件下作物生育期模型品种参数自动校正方法，其主要包括以下步骤：

1)历史气象数据准备及选择温光条件相似性年份的气象数据样本校正数据集；

2)通过参数敏感性分析确定待调参数，并确定参数范围；

3)通过约束性分析确定关键物候期，并进行缺失物候期实测数据范围估算；

4)设置自适应控制参数改进差分进化算法参数初始范围；

5)生成初始种群；

6)构造适应度函数；

7)构造差分矢量以对目标矢量

进行变异操作获得变异矢量

同时构造自适应缩放因子F_i；

8)将第j个变量的目标矢量

与变异矢量

进行交叉操作，获得试验个体

其中自适应交叉概率因子CR由自适应缩放因子F_i确定；

9)对目标矢量

和试验个体

进行选择操作，获得子代矢量

10)重复步骤5至步骤9，直至算法运行代数达到算法终止条件。

11)重复步骤4至步骤10，直至达到算法重复次数。

12)将步骤11)获得的参数组进行分析获得分类类别数目，进而用K-means算法聚类得到聚类中心作为该品种的代表性品种参数。

具体的，步骤1)中，作物生育期模型品种参数优化需准备的数据包括气象数据、不同品种在出苗、拔节、抽穗、成熟物候时期的历史实测数据以及相应的播期、播种深度实测数据，其中，采取DTW-SPXY组合算法获得对作物生育期模型品种参数优化有影响的温光条件气象数据，具体步骤为：

步骤1.1温光条件气象数据预处理

选取全年气象数据作为样本集，将其中的每日最高温、最低温进行相加除二，转换为每日平均气温构造年份-日均温表；将气象数据中的日日照时数提出，构造年份-日日照时数表；

步骤1.2DTW距离计算

对步骤1.1中建立的年份-日均温表、年份-日日照时数表，分别计算年与年之间的DTW 距离，形成年份间气温相似性度量矩阵表和年份间日照时数相似性度量矩阵表；

步骤1.3SPXY年份选择

对步骤1.2中的年份间气温相似性度量矩阵表和年份间日照时数相似性度量矩阵表各项分别除以自身矩阵中的最大项，而后将两张表按照气温表各项乘以0.7，日照表各项乘以0.3 后相加，合并为一张表，最后使用KS法选择具有代表性的年份，所述代表性的年份所对应的气象数据作为品种参数优化的数据。

具体的，步骤2)中，使用RiceGrow软件筛选出作物生育期模型敏感性大的品种参数作为待校正参数，根据品种参数的物理特性和经验确定模型品种参数初始范围。

具体的，步骤3)中确定对模型品种参数具有关键约束性的物候期，包括：

步骤3.1定义模型品种参数集合P＝{x₁,x₂,…,x_n}，模型模拟物候期集合 W＝{Y₁,Y₂,…,Y_n}；采用以下方法确定作物模型模拟物候期对品种参数的约束能力：通过公式(3)衡量物候期Y_i与品种参数x_i的约束关系以及正反比关系；若Z_i≠0，则Y_i约束x_i；若 Z_i＝0，则Y_i不约束x_i，x_i与Y_i不相关；若Z_i＞0时，x_i与Y_i成正相关；若Z_i＜0时，x_i与Y_i成负相关；

式中：X_i表示品种参数初始值，X′_i表示品种参数改变后的值，Y_i表示物候期i日期，Y_i+1表示物候期i+1的日期，Y′_i表示改变后X_i后物候期i的日期，Y′_i+1表示改变后X_i后物候期i+1 的日期，Z_i表示品种参数i的敏感度；

步骤3.2定义关键物候期集合K＝{k1，k2，…，kn}，采用以下规则确定作物模型参数优化的关键物候期，若物候期Y_i满足以下两个条件之一时认为该物候期为关键物候期，将该物候期Y_i加入关键物候期集合K中：

(1)若某一X_i仅被唯一的Ki约束时，则Ki为关键物候期；

(2)若X_i被Kj至Kk发育阶段中的物候期约束时，则选取Kj和Kk的两个物候期作为关键物候期；

步骤3.3估算关键物候期k1基础日期

首先判断实测物候期数据是否完全包含关键物候期集合K中全部元素的实测数据，若包含，则跳过该步骤；若缺失某一关键物候期k1的实测数据，则根据公式(4)推理估算该关键物候期k1的物候日期D_j；

式中：D_j表示待估算的物候日期，D_i表示已知物候期的日期，d_ij表示物候期i与物候期 j的持续天数，i，j表示物候期索引，i>j表示物候期j在物候期i之前，其中，d_ij的值由文献资料或栽培经验获取；

步骤3.4构造关键物候期实测数据的允许误差范围

构造关键目标物候期估算数据的运行误差范围[D_j-E,D_j+E]，E为估算物候期允许误差天数。

具体的，步骤6)中适应度函数的构造，主要由目标函数和惩罚函数两部分构成，包括：

步骤6.1利用模型模拟值和实际观测值之间的根均方差RMSE构造目标函数f_target，公式见(5)；当物候期实测数据是由步骤3)推算估计获得的范围时，对该范围进行一下处理，公式见(6)；

式中：P_i表示物候期模拟值，O_i表示物候期实测值，n表示物候期实测数据样本数，i表示1,2，…，n，E表示估算物候期允许误差天数，O_max表示估算范围的上限，O_min表示估算范围的下限；

步骤6.2构造惩罚函数，对关键物候期集合K中任一元素k_i的模拟效果极差的参数集进行惩罚，降低其传入下一代的概率，惩罚函数见公式(7)：

式中：f_target表示目标函数值，P_i表示物候期模拟值，O_i表示物候期实测值，ε表示误差临界值，超过ε即认为模拟效果极差；

步骤6.3首先采用公式(8)将最小化的目标函数转化为适应度函数，再对各目标的适应度函数值进行加权求和，从而计算个体的适应度函数值：

式中：

表示物候期j的目标函数值，f_fit表示适应度函数值，n表示关键物候期个数，w_j表示物候期j的权重。

具体的，步骤7)中构造差分矢量进行变异操作，同时构造自适应缩放因子，

差分进化算法DE最基本的变异成分是父代的差分矢量，差分矢量定义为

式中r₁和r₂表示种群中两个不同个体的索引号，

和

分别表示r₁和r₂对应的种群个体；将差分矢量加到另一个随机选择的个体矢量上，就生成了变异矢量；

对每一目标矢量

变异操作为：

其中，

为随机选择的三个个体矢量中适应度最优的个体，朝适应度更优的个体矢量方向搜索，会更加容易和快速找到全局最优解，将随机选择的三个个体矢量按适应度大小进行排序，适应度次优的个体记为

最差的个体记为

F为缩放因子，根据变异操作中随机选择的三个个体的相对位置来自适应的调节：如果随机选择的三个个体的位置相邻很近，则生成的变异矢量会很小，因此应该选择较大的F，否则变异量过小，起不到变异的效果；如果三个个体的位置很分散，则应该选择较小的F，否则变异量太大而使得生成的变异个体超出变量的搜索空间；

为变异矢量；

通过下式构造自适应缩放因子F_i：

式中，f_tb、f_tm、f_tw分别为

的适应度；式中，F_l＝0.1即F_i的下限， F_u＝0.9即F_i的上限，F_i在[0.1,0.9]之间随

和

在搜索空间中的相对位置而自适应的变化。

具体的，步骤8)中构造自适应交叉概率因子CR进行交叉操作，如下式所示；

其中，rand(j)∈[0,1]，为均匀分布的随机数，j表示第j个变量(基因)，CR为交叉概率因子；CR根据个体适应度与整个种群适应度的平均值的比较来确定：若个体适应度小于平均适应度，说明个体较优(最小化优化问题)，CR应该取较小的值，避免该个体基因结构遭到过多的破坏；若个体适应度大于平均适应度，说明该个体较差，CR应取较大的值，使之接受更多变异个体的基因，促使该个体的进化；

为变异矢量，

为第j个变量的目标矢量；

为交叉操作后生成的试验个体；

第i个目标矢量的交叉概率因子CR_i的取值公式为：

式中，f_i为第i个目标矢量的适应度,f_min为当前种群中适应度最优个体的适应度，f_max为当前种群中适应度最差个体的适应度，

为当前种群适应度的平均值；CR_l和CR_u分别为CR 的下限和上限。本发明实验研究表明，CR_l＝0.1，CR_u＝0.6对大多数优化问题比较适合。

具体的，步骤9)进行选择操作，差分进化算法DE采用“贪婪”的搜索策略，交叉操作后生成的试验个体

与目标矢量

进行竞争，只有当

的适应度较

更优时才被选做子代，否则，直接将

作为子代，选择操作的公式如式(13)：

其中，

表示

的适应度，

表示

的适应度，

表示子代矢量。

具体的，步骤12)使用改进的K-Means聚类算法对步骤11)获得的参数组选择一套代表性品种参数，具体的：

步骤11)获得的参数组为矩阵R(p,n)：

n表示作物生育期模型品种参数的数量，p表示步骤4至步骤10重复进行差分进化算法的次数，用x_i,j表示第i次搜索所得的第j个参数，R(p,n)中的每一行表示一组参数；

改进的K-Means聚类算法的具体步骤为：

步骤12.1：通过曲线图像观察p个参数组的参数分布决定起始质心的数量k；

步骤12.2：在p组参数组中任取k组作为起始质心；

步骤12.3：计算其余参数组到各个质心的距离，这个距离具体指的是两个参数组a(x₁₁,x₁₂,…,x_1s…,x_1n)与b(x₂₁,x₂₂,…,x_2s…,x_2n)之间的欧式距离，其公式为(14)：

其中：a表示一个起始质心参数组或一个质心参数组，b为一个非起始质心的参数组， x_1s为一个参数组中的第s个参数，x_2s为另一个参数组中的第s个参数；

步骤12.4：将步骤12.3中的非起始质心参数组，划归到与其距离最小的质心所在的类中，并重新计算该类的质心，其公式为(15)：

其中：Q为该质心所在的类包含的参数组数量，x_q1为其中一个参数组q的第1个参数， x_qn为其中一个参数组q的第n个参数，New为新的质心；

步骤12.5：若k个质心的值均不发生变化则执行步骤12.6，否则从步骤12.3开始重复执行；

步骤12.6：找出包含参数组数量最多的那个类L，计算L中各个参数组到该类的质心的距离，找出其中距离最小的一组R(x₁₁,x₁₂,…,x_1s…,x_1n)，即为最优参数。

优选的，所述作物为水稻。

本发明的有益效果

本发明的提出一种不确定条件下作物生育期模型品种参数自动校正方法，该方法一方面通过DTW-SPXY组合算法对作物全年气象数据进行选择，获得具代表性的作物生育期模型品种参数优化的样本集数据，使得校正模型的准确度提高，且受地点不同、水稻品种不同的影响较小。另一方面对参数进行敏感性分析，排除掉敏感性低的参数，通过经验或者品种特性确定其值，降低参数不确定性。最后对校正数据集做了筛选，选出有代表性的数据。选用的控制参数改进差分进化算法，一方面已有研究证明差分算法在解决此类非线性、多极值、复杂的模型表现出了较好的性能，另一方面，该算法的控制参数随进化过程中个体适应度的变化而自适应的调整，避免了原DE算法的控制参数选择对优化问题的依赖性，提高了算法的鲁棒性。最后对“异参同效”的品种参数进行聚类分析，选出一组有代表性的品种参数。

实际应用表明，本发明避免了传统作物生育期模型品种参数优化的数据准备步骤中选取的年份数据不具代表性的问题；同时解决了作物模型参数校正过程中数据集带来的不确定性以及参数带来的不确定性。并为RiceGrow这一类相同结构的作物生育期模型选择适合的进化算法，增强了模型品种参数优化的准确性，提高了作物模型品种参数的优化效率，有助于作物生长模型的发展与应用，具有科学意义和实用价值。

附图说明

图1本发明的流程示意图

图2为图1中气象基础数据准备的流程示意图

图3为不同生态点不同品种数据下基于不同进化算法的RiceGrow物候期模型参数校准的目标函数值

图4为实施例1中不同进化算法参数优化轨迹对比

图5为实施例1中宜兴-武育粳-SACPMDE-100次重复实验参数组分布图

图6为实施例2中基于多种训练集选择方法的参数优化对比流程图

具体实施方式

下面结合实施例对本发明作进一步说明，但本发明的保护范围不限于此：

具体实施方式中以水稻为例进行说明，其它作物如小麦等亦具备同样的效果。

实施例1：

以南京农业大学的RiceGrow水稻生育期模型模拟水稻品种武育粳、武育粳3号、芦丝占、雪花粘、两优培九在宜兴、兴化、高要和南京四个地点种植的物候期变化情况。水稻品种与地点对应关系为：宜兴-武育粳，兴化-武育粳3号，高要-芦丝占，高要-雪花粘，南京-两优培九。结合附图1，具体说明本实施方式：

步骤1：历史气象数据准备及选择温光条件相似性年份的气象数据样本校正数据集。

江苏兴化、南京及广东高要3个站点的逐日气象资料，来源于国家气象局气象信息中心。江苏太湖稻区(宜兴)和溧阳的逐日气象数据来源于江苏省气象局。气象资料包括日最高温度(℃)、日最低温度(℃)、日照时数(h/d)、降水量(mm)等。播种深度为 2.5cm。参照图2的流程进行物生育期模型品种参数优化的数据准备，最终获得具代表性的年份如表1所示。

表1试验站点及年份数据

南京2007—2008年的两优培九和宜兴1994—2011年的武育粳来源于国家信息农业工程技术中心水稻田间试验资料。兴化2001—2004年的武育粳3号、高要1997—2002 年的芦丝占和高要2004—2009年的雪花粘品种的实测数据来自中国气象局气象数据中心汇总的全国各水稻种植区田间试验资料。最终获得相应年份区间的栽培信息如表2所示。

表2水稻品种栽培信息

步骤2：通过参数敏感性分析确定待调参数为温度敏感因子TS、光周期敏感因子PS、最适温度To、品种基本早熟因子IE、基本灌浆因子FDF，并确定参数范围见表3；

表3 RiceGrow物候期模型的主要参数

步骤3：通过约束性分析确定关键物候期，并进行缺失物候期实测数据范围估算；

步骤3.1分析模型模拟输出物候期Y_i到Y_i+1的持续天数受到品种参数X_i变化所带来的影响。根据公式(3)计算获得Z_i，具体数据见表4。

表4RiceGrow水稻生育期模型各物候期持续天数随参数变化的情况

步骤3.2分析物候期Yi和品种参数Xi的约束关系发现，参数FDF被成熟期唯一约束， IE和PS被拔节至抽穗阶段约束，而TS则影响整个出苗至成熟阶段。因此根据调参关键物候期的确定规则，确定RiceGrow水稻生育期模型的拔节期、抽穗期和成熟期为调参关键物候期。

步骤3.3关键约束性物候期实测数据范围的估算

缺少抽穗期和拔节期的实测数据，需要利用当地栽培专家的经验知识加以推理估算数据范围。

估算关键物候期中抽穗期和拔节期的基础日期：本实施例中根据文献资料搜集和专家调研，估算缺失物候期。

步骤3.4构造关键物候期实测数据的允许误差范围

由估算得到关键约束性物候期中拔节期和抽穗期的实测数据具有一定误差。因此，需构造拔节期和抽穗期的实测数据的允许误差范围。根据当地栽培经验，抽穗期和拔节期的误差分别设置为1天和3天。

步骤4：设置自适应控制参数改进差分进化算法SCPMDEA参数初始范围；

设置种群规模BN＝100、最大迭代次数TM＝500、待优化参数个数PN＝5，算法终止条件为迭代次数达到最大迭代次数。

同时进行IAGA、M-精英、DE、PSO、ALEMSPSO算法的对比验证，各算法的具体参数如下表所示；

表5算法参数设置

步骤5-步骤11给出了本发明提出的SCPMDEA算法具体流程，其余验证算法均为现有技术，本文不再赘述。

步骤5：初始种群的生成

种群个体的编码方式采用实数编码，每一组个体对应一套作物模型品种参数，每个基因位对应一个品种参数。在该范围内随机生成BN个个体。

步骤6：构造目标函数

适应度函数是评价个体优劣的指标，主要由目标函数和惩罚函数两部分构成。

所述步骤6中包括：

步骤6.1目标函数设计

利用模型模拟值和实际观测值之间的根均方差(RMSE)构造目标函数，公式见(5)。当物候期实测数据是由步骤3推算估计获得的范围时，对该范围进行一下处理，公式见(6)。

式中：P_i表示物候期模拟值，O_i表示物候期实测值，n表示物候期实测数据样本数，i表示1,2，…，n，E表示估算物候期允许误差天数，O_max表示估算范围的上限，O_min表示估算范围的下限。

步骤6.2构造惩罚函数，

对关键物候期集合K中任一元素k_i的模拟效果极差的参数集进行惩罚，降低其传入下一代的概率，惩罚函数见公式(7)。

式中：f_target表示目标函数值，P_i表示物候期模拟值，O_i表示物候期实测值，ε表示误差临界值，本例取5。

步骤6.3首先采用公式(8)将最小化的目标函数转化为适应度函数，再对各目标的适应度函数值进行加权求和，从而计算个体的适应度函数值。

式中：

表示物候期j的目标函数值，f_fit表示适应度函数值，N表示关键物候期个数，w_j表示物候期j的权重。

步骤7：构造差分矢量以对目标矢量

进行变异操作获得变异矢量

同时构造自适应缩放因子F_i；

步骤8：将第j个变量的目标矢量

与变异矢量

进行交叉操作，获得试验个体

其中自适应交叉概率因子CR_i由自适应缩放因子F_i确定；

步骤9：对目标矢量

和试验个体

进行选择操作，获得子代矢量

步骤10：重复步骤5至步骤9，直至算法运行代数达到算法终止条件；

步骤11：重复步骤4至步骤10，直至达到算法重复次数；

最终获得不同水稻品种校准后模型参数的平均值及验证数据集下的评价指标如表6所示。

表6参数校正的平均值及测试数据集下的评价指标(括号内为标准偏差)

不同生态点不同品种数据下基于不同进化算法的RiceGrow物候期模型参数校准的目标函数值如图3所示。通过图3可以看出：差分算法系列其目标函数平均值最小，而目标函数最优值、最差值及平均值反应了算法的优化效果优劣，再次表明DE算法系列在RiceGrow 物候期模型品种参数校准上的有效性。DE算法在宜兴武育粳、兴化武育粳3号、高要雪花粘、高要芦丝占和南京两优培九数据集上模型参数校准的目标函数最优值分别为8.76、0.4、 2.6、1.52和2.51，而SACMPDE算法参数校准结果的目标函数最优值为8.65、0.4、2.6、 1.52和2.08，结果表明SACPMDE的目标函数最优值优于标准DE算法。

目标函数值的标准偏差反映了算法应用的稳定性，DE系列目标函数标准差都基本接近0，特别是SACPMDE在兴化武育粳3号和南京两优培九数据集上目标函数标准偏差都为0，说明算法的稳定性好，可以在有限时间范围内找到近似最优解。另外在达到相同精准度下，DE 和SACPMDE分别的收敛代数为421、108、347、318、358和354、106、186、200、121，SACPMDE收敛代数要比DE收敛代数分别小67、2、161、118、237，可以看出SACPMDE算法比DE算法收敛的快。其次协同进化遗传算法系列IAGA和MECA收敛代数相对较快，但其目标函数标准偏差在0.077～0.508范围内，原因可能是提早收敛到一个局部最优解，算法不够稳定。而PSO算法系列的目标函数值和标准偏差较前2类算法都大，表明在PSO系列在应用到RiceGrow物候期参数校准上模拟精度较低同时算法寻得最优解的能力相对不稳定。

以宜兴-武育粳为例，通过6类进化算法分别对RiceGrow物候期模型进行参数校准，重复试验100次的目标函数值随迭代次数增加画出其收敛曲线。不同进化算法参数优化轨迹对比如图4所示。可以看出SACPMDE算法和DE算法校准曲线尾部最窄，多次重复试验均收敛到近似同一值。但DE算法参数是由经验确定，且SACPMDE比DE收敛的速度平均快 117代。

步骤12：改进K-means算法求得“异参同效”参数组的聚类中心。改进的K-Means聚类算法的具体步骤为：

步骤12.1：通过曲线图像观察p个参数组的参数分布决定起始质心的数量k，以宜兴- 武育粳为例，曲线图像如图5所示；

步骤12.2：在p组参数组中任取k组作为起始质心；

步骤12.3：计算其余参数组到各个质心的距离，这个距离具体指的是两个参数组a(x₁₁,x₁₂,…,x_1s…,x_1n)与b(x₂₁,x₂₂,…,x_2s…,x_2n)之间的欧式距离，其公式为(14)：

其中：a表示一个起始质心参数组或一个质心参数组，b为一个非起始质心的参数组，x_1s为一个参数组中的第s个参数，x_2s为另一个参数组中的第s个参数；

步骤12.4：将步骤12.3中的非起始质心参数组，划归到与其距离最小的质心所在的类中，并重新计算该类的质心，其公式为(15)：

其中：Q为该质心所在的类包含的参数组数量，x_q1为其中一个参数组q的第1个参数， x_qn为其中一个参数组q的第n个参数，New为新的质心；

步骤12.5：若k个质心的值均不发生变化则执行步骤12.6，否则从步骤12.3开始重复执行；

步骤12.6：找出包含参数组数量最多的那个类L，计算L中各个参数组到该类的质心的距离，找出其中距离最小的一组R(x₁₁,x₁₂,…,x_1s…,x_1n)，即为最优参数。

按照上述步骤获得结果如表7所示。

表7SACPMDE参数校正“异参同效”不同方案参数选择校正评价

实施例2：

现有技术中作物生育期模型品种参数优化的历史气象样本数据选择常用的方法有：

1、气温隔年选择法

(1)气象数据预处理

对当地全年气象数据和全生育期气象数据分别实验，理由同上。将气象数据中每日最高温、最低温进行相加除二，转换为每日平均气温构造年份-日均温表。

(2)分别对全年和全生育期的年份-日均温表，计算日均温平均值，按照日均温的平均值大小排序，抽取其中第奇数个年份的数据加入调参训练集，其余年份作为验证数据集。

2、隔年选择法

直接抽取第奇数个年份的数据加入调参训练集，其余年份作为验证数据集。

为了验证本发明提出的作物生育期模型品种参数优化的气象样本数据选择方法(DTW- SPXY组合算法)的有效性，参照附图6，申请人进行基于多种训练集选择方法的参数。

实验样本来自广东省高要地区1997-2002年的气象数据(包括最低温、最高温、日照) 和1997-2002年水稻品种卢丝占在当地的种植实测数据(包括播种、移栽、拔节、抽穗、成熟的日期)，以及江苏省宜兴地区1994-2011年的气象数据(包括最低温、最高温、日照)和1994-2005以及2009-2011年水稻品种武育粳在当地种植的实测数据(包括播种、移栽、拔节、抽穗、成熟的日期)。如表8所示：

表8实验样本品种特性、种植地点气候特征以及气象数据的年份

结合三种作物生育期模型品种参数优化的数据准备方法，以及当地全年气象数据和该品种作物对应的全生育期气象数据两种数据长度，共进行5组对比实验，分别为：①生育期数据-SPXY，②全年数据-SPXY，③生育期数据-气温隔年，④全年数据-气温隔年，⑤隔年。分别按照5组方法选出其认为的具有代表性的年份如表9(高要)和表10(宜兴)所示：

表9高要1997-2002调参训练集年份选择结果

表10宜兴1994-2005、2009-2011调参训练集年份选择结果

将上述选出的训练集作为调参算法的输入，即可得到对应的优化后的模型温度敏感性因子(TS)、最适温度(To)、光周期敏感性因子(PS)、基本早熟性因子(IE)及基本灌浆期因子 (FDF)，结果如表11：

表11不同训练集调参结果

将表11中的参数优化结果带入RiceGrow模型(RiceGrow模型是南京农业大学国家信息农业工程技术中心研制的水稻生长模拟模型的简称。RiceGrow水稻生育期模型以水稻发育生理生态过程为基础，以生理发育时间(PDT)作为定量水稻发育进程的通用尺度。每日生理发育时间由每日生理效应累积形成，而每日热效应、光周期效应及品种的基本早熟性(抽穗前)或基本灌浆期因子(抽穗后)互作共同决定每日生理效应的大小。该模型具有较强的系统性、解释性与预测性，同时模型参数少、易确定，简便实用)，替换默认的温度敏感性因子、最适温度、光周期敏感性因子、基本早熟性因子及基本灌浆期因子，再用验证集年份的实测数据计算预测结果的RMSE和MAD值，得到结果如表12：

表12模型验证结果

首先从调参算法上来看，在高要卢丝占的数据上，生育期SPXY选择训练集、生育期气温选择训练集、全年SPXY选择训练集、全年气温选择训练集、隔年选择训练集的验证结果中，SACPMDE算法的RMSE值比IAGA算法的RMSE值分别小了0.5859、0.5679、0.1173、 0.0285、0.0988，即SACPMDE算法较之IAGA算法的RMSE值分别减小了7.75％、9.68％、 2.94％、0.77％、2.80％，综合5种训练集选择算法，SACPMDE算法较之IAGA算法的验证结果RMSE值平均减小了4.79％。

在宜兴武育粳的数据上，生育期SPXY选择训练集、生育期气温选择训练集、全年SPXY选择训练集、全年气温选择训练集、隔年选择训练集的验证结果中，SACPMDE算法的RMSE值比IAGA算法的RMSE值分别小了0.0988、0.2484、3.9620、0.6685、3.0095，即SACPMDE算法较之IAGA算法的RMSE值分别减小了2.51％、6.83％、58.79％、16.76％、52.20％，综合5种训练集选择算法，SACPMDE算法较之IAGA算法的验证结果RMSE值平均减小了27.42％。

从预测准确的差异大小上来看，两种调参算法在宜兴-武育粳样本集上的差异远大于两者在高要-卢丝占样本集上的差异，前者的RMSE差值为27.42％，而后者仅为4.79％。可见，在数据集数据量增大时，SACPMDE算法和IAGA算法的差异性将更加明显。

其次，从调参算法与训练集选择算法组合后的优化能力的不确定性来看。IAGA算法在高要卢丝占上的验证结果中，不同训练集选择算法的训练集的调参验证结果由好到差依次为：全年-SPXY、全年-气温、生育期-气温、隔年、生育期-SPXY。而IAGA算法在宜兴武育粳上的验证结果中，不同训练集选择算法的训练集的调参验证结果由好到差依次为：生育期- 气温、生育期-SPXY、全年-气温、隔年、全年-SPXY。两者之间差异较大，可见时间、地点、作物品种，均有可能对IAGA算法与训练集选择算法的组合效率产生较大影响。

而SACPMDE算法在高要卢丝占上的验证结果中，不同训练集选择算法的训练集的调参验证结果由好到差依次为：全年-SPXY、全年-气温、生育期-气温、隔年、生育期-SPXY。而SACPMDE算法在宜兴武育粳上的验证结果中，不同训练集选择算法的训练集的调参验证结果由好到差依次为：隔年、全年-SPXY、全年-气温、生育期-气温、生育期-SPXY。只有隔年选择法的位置发生了变动，其他训练集选择方法的调参准确度排序没有发生变化。可见SACPMDE算法和各种训练集选择算法的组合效率，受时间、地点、作物品种变化影响较小，算法与训练集选择这种数据预处理的方式结合效果较好。

由此可见，在同样的训练集的情况下，SACPMDE算法普遍优于IAGA算法。且在训练集数量增加时，SACPMDE算法调参结果更为准确的这种优势将被扩大。且SACPMDE调参算法下，对于不同的训练集选择算法，调参结果变化更为稳定、更为规律，不会因为时间地点的变化，而对不同的训练集选择算法选择代表性数据的能力，造成大的影响。因此，下面我们着重讨论SACPMDE调参算法下，训练集选择算法对模型调参结果的影响。

在SACPMDE算法下，当样本集为高要-卢丝占时，不同训练集选择算法的训练集的调参验证结果由好到差依次为：全年-SPXY、全年-气温、生育期-气温、隔年、生育期-SPXY。当样本集为宜兴-武育粳时，不同训练集选择算法的训练集的调参验证结果由好到差依次为：隔年、全年-SPXY、全年-气温、生育期-气温、生育期-SPXY。

首先从高要-卢丝占的数据来看。从数值上看，在高要-卢丝占的样本集上，从生育期- SPXY到生育期-气温，RMSE数值减少了1.7157，相当于在生育期-SPXY的基础上提上了24.45％。从生育期-气温到全年-SPXY，RMSE数值减少了1.4214，相当于在生育期-SPXY 的基础上提上了26.82％。从这里开始变化放缓，从生育期-气温到全年-SPXY，RMSE数值减少了0.1910，仅相当于在生育期-气温的RMSE值的4.92％，但仍可以看出两种选择算法之间的准确度差距。从全年-SPXY到全年-气温，RMSE值上升了0.2645，相当于全年-气温 RMSE值的7.17％。可以看出全年-SPXY算法选择的训练集，在模型预测的结果与实测数据的差异度上，优于其他几种选择算法。

对于宜兴-武育粳的样本集分析结果与高要-卢丝占基本一致，生育期-SPXY、生育期-气温、全年-SPXY、全年-气温、隔年选择的训练集的验证结果之间的RMSE差值分别为0.4532、0.0694、0.5417、0.0224，提升百分比分别为11.80％、2.05％、16.32、0.81％。

综合两个样本集，可以看出全年-SPXY选出的训练集的调参验证结果一直是较好的，且全年-SPXY、全年-气温、生育期-气温、生育期-SPXY四种选择方法的趋势基本保持一致。虽然在宜兴-武育粳的验证结果中，隔年的结果RMSE值好于全年-SPXY，但仅相差0.81％。再综合考虑隔年选择法在高要-卢丝占数据集上的结果，SPACMDE和全年-SPXY算法的组合准确度要高于其他四种选择算法。

另一方面，申请人进一步研究了当地全年气象数据和该品种作物对应的全生育期气象数据两种数据长度对最终模型精度的影响。比较全年-SPXY、全年-气温、生育期-气温、生育期-SPXY四种选择方法。在高要-卢丝占样本集上，全年-SPXY选择法的验证结果比生育期- SPXY选择法的验证结果的RMSE值小了3.5927，MAD值小了3.0476。全年-气温选择法的验证结果比生育期-气温选择法的验证结果的RMSE值小了0.1920，MAD值小了0.1905。

同样在宜兴-武育粳的样本集结果上，我们也可以发现相同的规律。全年-SPXY选择法的验证结果比生育期-SPXY选择法的验证结果的RMSE值小了1.0644，MAD值小了0.5476。全年-气温选择法的验证结果比生育期-气温选择法的验证结果的RMSE值小了0.069476022， MAD值相等。在宜兴-武育粳的样本集上两者的差距虽然有所减小，但全年样本集比生育期样本集调参的验证结果更为准确、更为稳定的趋势并没有变化。如表13所示：

表13不同样本集长度的调参结果

根据上述分析，可以发现以全年气象数据作为样本集的选择方法，较之同种算法以生育期数据为样本集的选择结果，其验证结果准确度，都是以全年数据为样本集的更为准确，同时结果数据分布也更为集中。

本文中所描述的具体实施例仅仅是对本发明精神做举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (8)

1.一种不确定条件下作物生育期模型品种参数自动校正方法，其特征在于其主要包括以下步骤：

1)历史气象数据准备及选择温光条件相似性年份的气象数据样本校正数据集；

步骤1)中，作物生育期模型品种参数优化需准备的数据包括气象数据、不同品种在出苗、拔节、抽穗、成熟物候时期的历史实测数据以及相应的播期、播种深度实测数据，其中，采取DTW-SPXY组合算法获得对作物生育期模型品种参数优化有影响的温光条件气象数据，具体步骤为：

步骤1.1温光条件气象数据预处理

选取全年气象数据作为样本集，将其中的每日最高温、最低温进行相加除二，转换为每日平均气温构造年份-日均温表；将气象数据中的日日照时数提出，构造年份-日日照时数表；

步骤1.2 DTW距离计算

对步骤1.1中建立的年份-日均温表、年份-日日照时数表，分别计算年与年之间的DTW距离，形成年份间气温相似性度量矩阵表和年份间日照时数相似性度量矩阵表；

步骤1.3 SPXY年份选择

对步骤1.2中的年份间气温相似性度量矩阵表和年份间日照时数相似性度量矩阵表各项分别除以自身矩阵中的最大项，而后将两张表按照气温表各项乘以0.7，日照表各项乘以0.3后相加，合并为一张表，最后使用KS法选择具有代表性的年份，所述代表性的年份所对应的气象数据作为品种参数优化的数据；

2)通过参数敏感性分析确定待调参数，并确定参数范围；

3)通过约束性分析确定关键物候期，并进行缺失物候期实测数据范围估算；

4)设置自适应控制参数改进差分进化算法参数初始范围；

5)生成初始种群；

6)构造适应度函数；

7)构造差分矢量以对目标矢量

进行变异操作获得变异矢量

同时构造自适应缩放因子F_i；

8)将第j个变量的目标矢量

与变异矢量

进行交叉操作，获得试验个体

其中自适应交叉概率因子CR由自适应缩放因子F_i确定；

9)对目标矢量

和试验个体

进行选择操作，获得子代矢量

10)重复步骤5至步骤9，直至算法运行代数达到算法终止条件；

11)重复步骤4至步骤10，直至达到算法重复次数；

12)将步骤11)获得的参数组进行分析获得分类类别数目，进而用K-means算法聚类得到聚类中心作为该品种的代表性品种参数。

2.根据权利要求1所述的参数自动校正方法，其特征在于步骤2)中，使用RiceGrow软件筛选出作物生育期模型敏感性大的品种参数作为待校正参数，根据品种参数的物理特性和经验确定模型品种参数初始范围。

3.根据权利要求1所述的参数自动校正方法，其特征在于步骤3)中确定对模型品种参数具有关键约束性的物候期，包括：

步骤3.1定义模型品种参数集合P＝{x₁,x₂,…,x_n}，模型模拟物候期集合W＝{Y₁,Y₂,…,Y_n}；采用以下方法确定作物模型模拟物候期对品种参数的约束能力：通过公式(3)衡量物候期Y_i与品种参数x_i的约束关系以及正反比关系；若Z_i≠0，则Y_i约束x_i；若Z_i＝0，则Y_i不约束x_i，x_i与Y_i不相关；若Z_i>0时，x_i与Y_i成正相关；若Z_i<0时，x_i与Y_i成负相关；

式中：X_i表示品种参数初始值，X′_i表示品种参数改变后的值，Y_i表示物候期i日期，Y_i+1表示物候期i+1的日期，Y′_i表示改变后X_i后物候期i的日期，Y′_i+1表示改变后X_i后物候期i+1的日期，Z_i表示品种参数i的敏感度；

步骤3.2定义关键物候期集合K＝{k1，k2，…，kn}，采用以下规则确定作物模型参数优化的关键物候期，若物候期Y_i满足以下两个条件之一时认为该物候期为关键物候期，将该物候期Y_i加入关键物候期集合K中：

(1)若某一X_i仅被唯一的Ki约束时，则Ki为关键物候期；

(2)若X_i被Kj至Kk发育阶段中的物候期约束时，则选取Kj和Kk的两个物候期作为关键物候期；

步骤3.3估算关键物候期k1基础日期

首先判断实测物候期数据是否完全包含关键物候期集合K中全部元素的实测数据，若包含，则跳过该步骤；若缺失某一关键物候期k1的实测数据，则根据公式(4)推理估算该关键物候期k1的物候日期D_j；

式中：D_j表示待估算的物候日期，D_i表示已知物候期的日期，d_ij表示物候期i与物候期j的持续天数，i，j表示物候期索引，i>j表示物候期j在物候期i之前，其中，d_ij的值由文献资料或栽培经验获取；

步骤3.4构造关键物候期实测数据的允许误差范围

构造关键目标物候期估算数据的运行误差范围[D_j-E,D_j+E]，E为估算物候期允许误差天数。

4.根据权利要求1所述的参数自动校正方法，其特征在于步骤6)中适应度函数的构造，主要由目标函数和惩罚函数两部分构成，包括：

步骤6.1利用模型模拟值和实际观测值之间的根均方差RMSE构造目标函数f_target，公式见(5)；当物候期实测数据是由步骤3)推算估计获得的范围时，对该范围进行一下处理，公式见(6)；

式中：P_i表示物候期模拟值，O_i表示物候期实测值，n表示物候期实测数据样本数，i表示1,2，…，n，E表示估算物候期允许误差天数，O_max表示估算范围的上限，O_min表示估算范围的下限；

步骤6.2构造惩罚函数，对关键物候期集合K中任一元素k_i的模拟效果极差的参数集进行惩罚，降低其传入下一代的概率，惩罚函数见公式(7)：

式中：f_target表示目标函数值，P_i表示物候期模拟值，O_i表示物候期实测值，ε表示误差临界值，超过ε即认为模拟效果极差；

步骤6.3首先采用公式(8)将最小化的目标函数转化为适应度函数，再对各目标的适应度函数值进行加权求和，从而计算个体的适应度函数值：

式中：

表示物候期j的目标函数值，f_fit表示适应度函数值，n表示关键物候期个数，w_j表示物候期j的权重。

5.根据权利要求1所述的参数自动校正方法，其特征在于步骤7)中构造差分矢量进行变异操作，同时构造自适应缩放因子，

对每一目标矢量

变异操作为：

其中，

为随机选择的三个个体矢量中适应度最优的个体，朝适应度更优的个体矢量方向搜索，会更加容易和快速找到全局最优解，将随机选择的三个个体矢量按适应度大小进行排序，适应度次优的个体记为

最差的个体记为

F为缩放因子，根据变异操作中随机选择的三个个体的相对位置来自适应的调节；

为变异矢量；

通过下式构造自适应缩放因子F_i：

式中，f_tb、f_tm、f_tw分别为

的适应度；式中，F_l＝0.1即F_i的下限，F_u＝0.9即F_i的上限，F_i在[0.1,0.9]之间随

和

在搜索空间中的相对位置而自适应的变化。

6.根据权利要求1所述的参数自动校正方法，其特征在于步骤8)中构造自适应交叉概率因子CR进行交叉操作，如下式所示；

其中，rand(j)∈[0,1]，为均匀分布的随机数，j表示第j个变量，CR为交叉概率因子；CR根据个体适应度与整个种群适应度的平均值的比较来确定：若个体适应度小于平均适应度，CR应该取较小的值，避免该个体基因结构遭到过多的破坏；若个体适应度大于平均适应度，CR应取较大的值，使之接受更多变异个体的基因，促使该个体的进化；

为变异矢量，

为第j个变量的目标矢量；

为交叉操作后生成的试验个体；

第i个目标矢量的交叉概率因子CR_i的取值公式为：

式中，f_i为第i个目标矢量的适应度,f_min为当前种群中适应度最优个体的适应度，f_max为当前种群中适应度最差个体的适应度，

为当前种群适应度的平均值；CR_l和CR_u分别为CR的下限和上限。

7.根据权利要求1所述的参数自动校正方法，其特征在于步骤9)进行选择操作，差分进化算法DE采用“贪婪”的搜索策略，交叉操作后生成的试验个体

与目标矢量

进行竞争，只有当

的适应度较

更优时才被选做子代，否则，直接将

作为子代，选择操作的公式如式(13)：

其中，

表示

的适应度，

表示

的适应度，

表示子代矢量。

8.根据权利要求1所述的参数自动校正方法，其特征在于步骤12)使用改进的K-Means聚类算法对步骤11)获得的参数组选择一套代表性品种参数，具体的：

步骤11)获得的参数组为矩阵R(p,n)：

n表示作物生育期模型品种参数的数量，p表示步骤4至步骤10重复进行差分进化算法的次数，用x_i,j表示第i次搜索所得的第j个参数，R(p,n)中的每一行表示一组参数；

改进的K-Means聚类算法的具体步骤为：

步骤12.1：通过曲线图像观察p个参数组的参数分布决定起始质心的数量k；

步骤12.2：在p组参数组中任取k组作为起始质心；

步骤12.3：计算其余参数组到各个质心的距离，这个距离具体指的是两个参数组a(x₁₁,x₁₂,…,x_1s…,x_1n)与b(x₂₁,x₂₂,…,x_2s…,x_2n)之间的欧式距离，其公式为(14)：

其中：a表示一个起始质心参数组，b为一个非起始质心的参数组，x_1s为一个参数组中的第s个参数，x_2s为另一个参数组中的第s个参数；

步骤12.4：将步骤12.3中的非起始质心参数组，划归到与其距离最小的质心所在的类中，并重新计算该类的质心，其公式为(15)：

其中：Q为该质心所在的类包含的参数组数量，x_q1为其中一个参数组q的第1个参数，x_qn为其中一个参数组q的第n个参数，New为新的质心；

步骤12.5：若k个质心的值均不发生变化则执行步骤12.6，否则从步骤12.3开始重复执行；

步骤12.6：找出包含参数组数量最多的那个类L，计算L中各个参数组到该类的质心的距离，找出其中距离最小的一组R(x₁₁,x₁₂,…,x_1s…,x_1n)，即为最优参数。

Images