约束性知识与精英个体策略遗传算法融合的作物生育期模型品种参数优化方法

摘要

本发明提出一种约束性知识引导精英个体策略遗传算法的作物生育期模型参数优化方法，属于作物模型品种参数优化的最优化方法领域。其步骤：作物生育期模型品种参数优化的数据准备；作物生育期模型品种参数初始范围的等级划分；确定对模型品种参数具有关键约束性的物候期；关键约束性物候期实测数据范围的估算；设置改进型遗传算法控制参数的初始值；生成初始种群；构造适应度函数；构造个体的基因位方向因子；采用轮盘赌方法构造选择算子；构造方向引导交叉算子；构造方向与精英个体引导的变异算子；基于精英个体的局部搜索。本发明克服了效率低、优化结果不符合品种生物学含义等问题，增加了参数优化的准确性和效率，具有科学意义和实用价值。

说明书

一、技术领域

本发明属于农业信息技术领域，是领域知识引导进化算法与作物生长模型的交叉领域。涉及一种可用于非线性、不连续、多参数的作物生育期机理模型品种参数估算的自动优化方法，可用于小麦、水稻等禾谷类作物。具体地讲，该方法是一种约束性知识与精英个体策略遗传算法融合的作物生育期模型品种参数优化方法。

二、背景技术

气候变化问题已成为国际社会关注的热门话题，被列为全球十大环境问题之一。20世纪80年代初，科学家们开始发现全球的气候正发生着显著的变化。在2007年气候变化专门委员会(IPCC)发布的第四次评估报告指出：全球气候变暖已是不争的事实。农业生产与气候息息相关，气候变暖导致的极端气候事件趋多、趋强。这势必会对农业及生态环境等领域产生深远的影响。气候变化会导致粮食生产的不稳定。因此，气候变暖对农业生产的影响评价研究近年来受到重视，一方面积极开展作物对气候变化响应的观测试验研究，另一方面通过作物生长模型模拟评估气候变化对作物系统的可能影响，为应对气候变化条件下，农业生产策略与措施变化提供基础数据与决策支持。与观测试验研究相比，利用作物生长模型进行请模拟研究，由于它具有机理强、成本低、时效高、变量易于控制等优点，目前已成为气候变化影响评估的主要研究方法，被世界各地的研究者普遍采用。

作物生长模型是基于作物生理学和生态学原理，运用系统动力学的方法，综合考虑作物品种、栽培措施与气候环境、土壤条件之间的互作关系，量化作物生长发育、光合生产、器官建成及产量形成等生理过程的数学模型。由于作物生长模型充分考虑了作物对土壤、天气、水分、养分等多种影响因子的机理响应，能够模拟不同自然环境和管理措施下作物的生长状况，因此是有效探索气候变化对农业生产的影响，评价作物系统生产力的有力工具。作物生育期模型是作物生长模型的重要组成部分，它通过定量模拟作物对温度和光照的反应，预测和划分作物的发育进程，是作物光合作用、干物质积累与分配、养分吸收与转移、产量品质形成等方面模拟的基础，直接影响作物生长模型的模拟结果。

作物生长模型通过品种参数量化不同品种的遗传与生理特性。由于自然环境复杂性，作物生长模型应用中具有较强的不确定性，需要根据应用品种和应用区域的试验数据，对模型品种参数进行校正。目前作物生长模型参数校正方法主要分为知识驱动法以及数据驱动法。知识驱动法是指作物模型专家以作物模型建模机理、作物栽培经验、作物品种特性等专业知识为引导，采用手工试错或计算机循环迭代的方式，在分阶段、分目标等策略指导下不断试错寻找最优的参数集合。该方法以领域知识为支撑，得到的结果较为可靠，但当应用尺度增大时，为了使模型对全部试验资料均取得较好的模拟效果，调参的数据量及工作量增大，要求调参人员具有较为完备的作物模型知识和丰富的模型调参经验，因而受调参人员的主观因素影响较大；且循环试错的过程需要不断重复迭代，费事费力，参数优化效率较低。数据驱动法是将该问题转化为最优化问题，采用演化计算等全局优化算法，通过大量的随机搜索试验寻找模型参数的最优解。遗传算法是一种常用的随机化、启发式优化算法，具有鲁棒性强、搜索效率高等显著特点，适合求解作物长模型参数优化这样一类结构复杂、非线性程度高、搜索空间巨大且非单峰的问题，已经在ORYZA2000、CropSyst、RiceGrow、SUCROS等模型中取得了应用。现有研究更为关注优化算法本身的改进，如在种群划分时加入动态聚类的思想、采用两级遗传算法的方式进行参数优化等，提升了算法的优化性能。但在这类算法改进中，往往进化寻优的角度增加算法的优化策略，往往控制参数众多，算法时间和空间复杂度变高。此外，作物生育期模型的品种参数反应了作物品种的遗传特性，而基于遗传算法的作物模型参数的校正与优化仅从数值寻优的角度寻找品种参数最优解，往往忽略了作物模型品种参数所表达的生理遗传特性，应用过程中会造成优化后的品种参数出现不符合作物品种生物学意义的现象。

本发明采用两种方法有效的融合的方式，将作物品种及模型知识融合到改进型遗传算法中，设计了一种适用于作物生育期模型品种参数自动校正的方法，对快速而准确的优化作物模型参数具有重要的实际应用价值。

三、发明内容

本发明提出一种约束性知识与精英个体策略遗传算法融合的作物生育期模型参数自动优化方法(Constraints Knowledge and Elite Individual Strategy GeneticAlgorithm Fusion，简称CK-EISGA)。该方法将品种特征知识与改进型精英策略遗传算法有效融合，解决作物生育期模型品种参数自动校正时，由于模型算法多峰非线性以及进化算法优化结果的生物学意义不明确等问题，可快速准确地估算作物生育期模型品种参数。

本发明所述的约束性知识引导精英策略遗传算法的作物生育期模型品种参数优化方法，包括以下步骤：

步骤1：驱动数据准备

作物生育期模型品种参数优化需准备的基础数据包括气象数据(最低温、最高温、日照时数或辐射量)、不同品种在出苗、拔节、抽穗、成熟等物候时期的历史实测数据以及相应的播期、播种深度等实测数据。

步骤2：品种参数初始范围的等级划分

作物生育期模型的品种参数量化了作物生长发育的遗传特性，建模时确定的初始取值范围通常涵盖了某种类型作物的遗传特性，而没有针对一类具体品种遗传特性的强弱进一步细化，会造成优化算法搜索空间的扩大，参数优化后会出现模型模拟值与实测值误差虽小，但品种参数已失去生物学含义的问题。因此，品种参数的初始范围应根据具体品种的遗传特性缩小限制在某一区间内，按公式(1)(2)对品种参数范围进行等级划分。

$>R_{min}^{ij}=R_{\min }^i+j*\frac{1}{N}*(R_{max}^i-R_{\min }^i)---\left(1\right)$>

$>R_{max}^{ij}=R_{min}^i+(j+1)*\frac{1}{N}*(R_{\max }^i-R_{min}^i)---\left(2\right)$>

式中：表示品种参数i原始范围最大值，表示品种参数i原始范围最小值，表示参数i重新划分等级后，第j级特性的最小值，表示参数i重新划分等级后，第j级特性的最大值，N表示参数i品种特性的等级个数。

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

作物生育期模型模拟了多个物候阶段，每个物候阶段通常受到不同发育因子的影响且影响程度不同，造成了各物候阶段对品种参数的约束程度不同。关键物候期是指对待优化品种参数具有明显约束能力的物候期。只有当选用了关键物候期作为模拟目标，并分析与品种参数的数据相关性，才能使优化效果更为准确。

所述步骤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。

$>Z_i=\frac{(Y_{i,j+1}-Y_{i,j})-(Y_{i,j+1}^{\prime }-Y_{i,j}^{\prime })}{X_i^{\prime }-X_i}---\left(3\right)$>

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

步骤3.2定义关键物候期集合K＝{k₁，k₂，…，k_n}，采用以下规则确定作物模型参数优化的关键物候期，：

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

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

步骤3.3判定品种参数与模型模拟物候期的相关性

根据步骤3.1得到的Z的情况能够确定品种参数X_i与物候期持续天数之间存在一定的相关性，当Z＞0时，X_i与成正相关，当Z＜0时，X_i与成负相关，当Z＝0时，X_i与不相关。

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

由于作物生长需要延续若干个月，品种参数优化时所选用的物候期实测数据往往不够完备，会出现缺少步骤2确定的部分关键约束性物候期k_i的实测数据，需要利用当地栽培专家的经验知识加以推理估算数据范围。

所述步骤4中包括：

步骤4.1估算关键物候期k₁的基础日期

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

$>D_j=\left(\begin{array}{cc}{D}_i+d_{ij},& i<j\\ D_i-d_{ij},& i>j\end{array}\right)---\left(4\right)$>

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

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

由估算得到关键物候期k_i的日期具有一定误差。因此，构造关键物候期估算数据的运行误差范围[D_j-E，D_j+E]E为相差范围。

步骤5设置改进型遗传算法控制参数的初始值

主要包括：种群规模BN、最大迭代次数TM、交叉概率Pc、变异概率Pm、待优化参数个数PN、及算法终止条件等；算法终止条件可采用到达最大迭代次数或达到收敛条件，由用户自定义。

步骤6初始种群的生成

种群个体采用实数编码方式，每一个个体对应一套模型品种参数X＝{x_i，x₂，…，x_n}。在步骤2划分的品种参数范围内随机生成BN套个体。

步骤7适应度函数的构造

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

所述步骤7中包括：

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

$>f_{t\arg et}=\sqrt{{\Sigma }_{i=1}^n\frac{{(P_i-O_i)}^2}{n}}---\left(5\right)$>

$>O_i=\left(\begin{array}{cc}{O}_{max},& P_i-O_i>E\\ O_i,& |P_i-O_i|<E\\ O_{\min },& P_i-O_i<-E\end{array}\right)---\left(6\right)$>

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

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

$>f_{t\arg et}=\left(\begin{array}{cc}{f}_{t\arg et},& P_i-O_i>ϵ\\ 2*f_{t\arg et},& P_i-O_i>ϵ\end{array}\right)---\left(7\right)$>

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

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

$>f_{fit}={\Sigma }_{j=1}^N{W}_j\frac{1}{1+f_{t\arg et}^j}---\left(8\right)$>

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

步骤8构造个体的基因位方向因子

传统遗传算法的交叉和变异算子通常缺少方向性。而在交叉和变异过程中，某一基因位的交叉或变异方向与下一代期望调整方向有可能不一致，导致了作物模型参数优化背景下个体更新的无效搜索。由步骤2.3可以获得，作物生育期模型的品种参数的变化规律与各物候期模拟结果存在着正相关或负相关性，评价个体时，能够得到个体物候期模拟结果与实测结果之间的正负关系，利用该正负关系，结合模型参数与物候期间的相关性，能够对个体中每一基因位(即品种参数)下一代更新期望进化的方向加以有效引导，提高全局搜索的效率。

作物模型品种参数个体的基因位期望进化方向可由公式(9)计算。

$>d_i=\left(\begin{array}{cc}1,& (P-O)*Z_i<0\\ 0,& (P-O)*Z_i=0\\ -1,& (P-O)*Z_i>0\end{array}\right)---\left(9\right)$>

式中：d_i表示个体中第i个基因位(品种参数)的期望调整方向，P表示任意一个约束个体中第i基因位的关键物候期模拟值，O表示任意一个约束个体中第i基因位的关键物候期实测值，Z_i表示个体中第i个品种参数与物候期的正反比关系，计算公式见公式(3)。

步骤9采用轮盘赌选择算子，首先根据个体适应度值大小，计算个体遗传进入下一代的概率。根据轮盘赌选择原则，从父代中选取部分个体。

步骤10：构造方向引导交叉算子。在模拟二进制交叉算子的基础上，结合个体基因位方向引导因子，对各基因的交叉方向进行引导，如下式所示：

$>y_i^1=\left(\begin{array}{cc}{x}_i^1+0.5[(x_i^1-x_i^2)-\beta (x_i^1-x_i^2)],& (x_i^1-x_i^2)*d_i^1>0\\ x_i^1-0.5[(x_i^1-x_i^2)-\beta (x_i^1-x_i^2)],& (x_i^1-x_i^2)*d_i^1<0\end{array}\right)---\left(10\right)$>

$>y_i^2=\left(\begin{array}{cc}{x}_i^2+0.5[(x_i^2-x_i^1)+\beta (x_i^1-x_i^2)],& (x_i^2-x_i^1)*d_i^2>0\\ x_i^2-0.5[(x_i^2-x_i^1)+\beta (x_i^1-x_i^2)],& (x_i^2-x_i^1)*d_i^2<0\end{array}\right)---\left(11\right)$>

$>\beta =\left(\begin{array}{cc}{\left(2u\right)}^{\frac{1}{{\eta }_c+1}}& u\le 0.5\\ {\left(2\right(1-u\left)\right)}^{\frac{1}{{\eta }_c+1}}& u>0.5\end{array}\right)---\left(12\right)$>

式中，表示两个父代个体的第i位基因，表示两个子代个体的第i位基因，U表示区间[0，1]内的均匀分布随机数，η_c表示交叉参数，表示两个父代个体第i位基因方向因子。

步骤11：构造方向与精英个体引导的变异算子。在非均匀变异算子的基础上，利用个体基因位方向引导因子与种群精英个体综合引导变异算子，从而加快算法收敛到全局最优解的速度，公式见(13)。

$>X_k^{\prime }=X_k+d_k(U_{max}^k-X_k)(1-r^{(1-\frac{t}{T})b})+Cr(X_k^{Best}-X_k)---\left(13\right)$>

式中：X′_k表示变异后第k位基因，X_k表示变异前第k位基因，d_k表示第k位基因的方向因子，γ表示[0，1]内的均匀分布随机数，t表示当前进化代数，T表示最大进化代数，b表示非均匀度决定系数，c表示加速系数，表示种群中第k位基因的最大值，表示最优个体第k位基因。步骤12：基于精英个体的局部搜索。由于种群精英个体与全局最优解之间的亲和度要大于种群中的其他个体，利用精英个体信息产生多个与精英个体有较大亲和度的个体，可以提高种群优势个体的多样性，增加算法搜索到全局最优解的概率。因此在精英个体附近进行随机扰动，产生一定数量的新个体，提高算法后期局部搜索能力。新个体产生公式如下：

$>X_k^{\prime }=X_k^{Best}\pm (U_{max}^k-U_{min}^k)/r---\left(14\right)$>

式中：X′_k表示新个体第k位基因，表示最优个体第k位基因，表示种群中第k位基因的最大值，)表示种群中第k位基因的最小值，r表示[0，1]内的均匀分布随机数。

步骤13重复步骤7至步骤12，直至算法运行代数达到算法终止条件。

本发明与现有的作物模型参数估算方法相比，其改进主要体现如下：

(1)根据待优化品种的生物学特性对作物生育期模型的品种参数初始范围进行了进一步划分，，缩小遗传算法随机搜索空间，解决了基于遗传算法进行作物模型参数优化可能得到的不符合品种生物学特性的参数集问题。

(2)给出了关键约束性物候期的确定方法和日期范围的估算方法，为作物生育期模型参数优化的目标变量选择提供了有效方法。

(3)构造模型输出物候期与品种参数间的基因位方向引导因子，对遗传算法的交叉变异操作进行引导，提高了进化算法在该类应用中有效搜索能力。

(4)个体变异操作时加入精英个体的信息进行引导，提高了算法向最优解附近靠拢的速度。并对精英个体附近进行局部搜素，提高了种群中优势个体的多样性，增加算法搜索到全局最优解的概率。

有益效果：

本发明提出一种约束性知识与精英个体策略遗传算法融合的作物生育期模型品种参数优化方法，该方法一方面将作物品种特性的相关约束性知识引入到改进型遗传算法优化框架中，得出约束参数的初始范围和关键物候期的实测值范围，使得到的品种参数更符合作物品种的生物学特性。另一方面通过改进具有方向性的精英个体遗传算法的种群更新策略，提高了算法寻找最优解的有效性和寻优效率。

实际应用表明，本发明克服了传统参数优化方法中由于大量的迭代试错而带来的优化效率低下、调参周期长等问题，并解决了由于传统遗传算法的随机搜索机制而带来的优化结果不符合品种生物学含义的问题，增强了模型品种参数优化的准确性，提高了作物模型品种参数的优化效率，有助于作物生长模型的发展与应用，具有科学意义和实用价值。

四、说明书附图

图1本发明的流程示意图

图2基于精英个体策略遗传算法与蜜蜂进化型遗传算法的种群进化分布图

图3各项约束性知识的改进效果(图3-1为作物模型品种范围约束对参数PVT分布结果的影响，图3-2为关键约束性物候期实测数据约束对品种参数IE分布结果的影响。图3-3为方向因子对收敛速率的影响)

图4为本方法在徐州、济宁、潍坊、赣榆四个地点模拟效果的1∶1关系图；

图5为徐州、济宁、潍坊、赣榆四地的历史实测数据(图5-1为开花期和成熟期的实测数据，图5-2为四个地点的平均温度变化图，图5-3四个地点的辐射量变化图)。

五、具体实施方式

以下通过实施案例对本发明进一步描述

实施例

以南京农业大学的WheatGrow小麦生育期模型模拟济南13号小麦品种在山东省济宁市、胶州市、潍坊市和江苏省徐州市四个地点种植的物候期变化情况。结合附图1，具体说明本实施方式：步骤1：作物生育期模型品种参数优化的数据准备

选取济南13号小麦品种在山东省济宁市、胶州市、潍坊市和江苏省徐州市四个地点种植的历史数据，试验地点及年份数据见表1，各套数据均记录开花期和成熟期的实测日期，见附图5-1。气象数据均来自当地气象站，包括日最高温度(℃)、日最低温度(℃)、辐射量(MJ/m2))等，见附图5-2和5-3。播种深度为2.5cm。

表1试验选用的地点及年份数据

步骤2：WheatGrow小麦生育期模型品种参数范围划分。

小麦为冷季作物，其生育期的特性包括感温性、感光性、基本早熟型和春化效应，WheatGrow的品种参数TS、PS、IE和PVT分别表现了以上四种特性，根据公式(1)对品种参数初始范围进行划分，结果见表2。

表2 WheatGrow小麦生育期模型参数初始范围划分及生物学含义

步骤3：确定WheatGrow小麦生育期模型具有关键约束性的物候期

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

表3 WheatGrow小麦生育期模型各物候期持续天数随参数变化的情况

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

步骤3.3确定品种参数与模型模拟物候期的相关性

根据步骤2.1得到的Zi的情况能够确定品种参数X_i与物候期Y_i持续天数之间存在一定的相关性。根据表4中Zi值，确定品种参数TS、PS和PVT与物候期持续天数之间呈正比关系，品种参数IE、FD与物候期持续天数之间呈反比关系。

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

徐州等四个地点的实测数据仅包含开花期和成熟期，缺少抽穗期和拔节期的实测数据，需要利用当地栽培专家的经验知识加以推理估算数据范围。

步骤4.1估算关键物候期中抽穗期和拔节期的基础日期

根据文献资料搜集和专家调研，确定济南13号品种在江苏徐州等地种植时，拔节期距抽穗期约35天，抽穗期距开花期约5天。徐州地点实测数据中1982年的开花期为5月2日，则根据公式(4)，估算抽穗期为4月27日，拔节期为3月23日。其他地点和年份的估算方法与上例相同。

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

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

步骤5设置算法的初始控制参数

设置种群规模BN＝100、最大迭代次数TM＝100、交叉概率Pc＝0.8、变异概率Pm＝0.05、待优化参数个数PN＝5，算法终止条件为迭代次数达到最大迭代次数。

步骤6初始种群的生成

种群个体的编码方式采用实数编码，每一组个体对应一套作物模型品种参数，每个基因位对应一个品种参数。由于济南13号品种的春化特性为半冬性，因此根据表3，将PVT参数的初始范围设置为[30，40]，其他参数为初始参数。在该范围内随机生成BN个个体。

步骤7适应度函数的构造

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

所述步骤7中包括：

步骤7.1目标函数设计

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

$>f_{t\arg et}=\sqrt{{\Sigma }_{i=1}^n\frac{{(P_i-O_i)}^2}{n}}---\left(5\right)$>

$>O_i=\left(\begin{array}{cc}{O}_{max},& P_i-O_i>E\\ O_i,& |P_i-O_i|<E\\ O_{\min },& P_i-O_i<-E\end{array}\right)---\left(6\right)$>

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

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

$>f_{t\arg et}=\left(\begin{array}{cc}{f}_{t\arg et},& P_i-O_i>ϵ\\ 2*f_{t\arg et},& P_i-O_i>ϵ\end{array}\right)---\left(7\right)$>

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

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

$>f_{fit}={\Sigma }_{j=1}^N{W}_j\frac{1}{1+f_{t\arg et}^j}---\left(8\right)$>

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

步骤8构造方向因子

传统遗传算法的交叉和变异算子通常缺少方向性。而在交叉变异等算子中，某一基因位的交叉或变异方向与其期望调整方向有可能不一致，导致了无效搜索。由于作物生育期模型的品种参数的变化规律与各物候期模拟结果存在一定相关性，主要表现为正相关和负相关，由步骤3.3可以获得。在评价个体时，能够得到个体物候期模拟结果与实测结果之间的正负关系，而利用该正负关系，结合模型参数与物候期间的相关性，能够对个体中每一基因位下一步期望进化的方向进行引导。

作物模型品种参数个体中基因位期望进化的方向可由公式(9)

$>d_i=\left(\begin{array}{cc}1,& (P-O)*Z_i<0\\ 0,& (P-O)*Z_i=0\\ -1,& (P-O)*Z_i>0\end{array}\right)---\left(9\right)$>

式中：d_i表示个体中第i个品种参数的期望调整方向，P表示任意一个约束个体中第i基因位的关键物候期模拟值，O表示任意一个约束个体中第i基因位的关键物候期实测值，Z_i表示个体中第i个品种参数与物候期的正反比关系，计算公式见步骤3.3。

步骤9采用轮盘赌选择算子，首先根据个体适应度值大小，计算个体遗传进入下一代的概率。根据轮盘赌选择原则，从父代中选取部分个体。

步骤10：构造方向引导交叉算子。在模拟二进制交叉算子的基础上，结合个体基因位的方向引导因子，对各基因的交叉方向进行引导，如下式所示：

$>y_i^1=\left(\begin{array}{cc}{x}_i^1+0.5[(x_i^1-x_i^2)-\beta (x_i^1-x_i^2)],& (x_i^1-x_i^2)*d_i^1>0\\ x_i^1-0.5[(x_i^1-x_i^2)-\beta (x_i^1-x_i^2)],& (x_i^1-x_i^2)*d_i^1<0\end{array}\right)---\left(10\right)$>

$>y_i^2=\left(\begin{array}{cc}{x}_i^2+0.5[(x_i^2-x_i^1)+\beta (x_i^1-x_i^2)],& (x_i^2-x_i^1)*d_i^2>0\\ x_i^2-0.5[(x_i^2-x_i^1)+\beta (x_i^1-x_i^2)],& (x_i^2-x_i^1)*d_i^2<0\end{array}\right)---\left(11\right)$>

$>\beta =\left(\begin{array}{cc}{\left(2u\right)}^{\frac{1}{{\eta }_c+1}}& u\le 0.5\\ {\left(2\right(1-u\left)\right)}^{\frac{1}{{\eta }_c+1}}& u>0.5\end{array}\right)---\left(12\right)$>

式中，表示两个父代个体的第i位基因，表示两个子代个体的第i位基因，U表示区间[0，1]内的均匀分布随机数，η_c表示交叉参数，一般取值为2，表示两个父代个体第i位基因方向因子。

步骤11：构造方向和精英个体引导的变异算子。在非均匀变异算子的基础上，利用个体基因位的方向引导因子和种群精英个体引导变异算子，从而加快算法收敛到全局最优解的速度，公式见(13)。

$>X_k^{\prime }=X_k+d_k(U_{max}^k-X_k)(1-r^{(1-\frac{t}{T})b})+Cr(X_k^{Best}-X_k)---\left(13\right)$>

式中：X′_k表示变异后第k位基因，X_k表示变异前第k位基因，d_k表示第k位基因的方向因子，r表示[0，1]内的均匀分布随机数，t表示当前进化代数，T表示最大进化代数，b表示非均匀度决定系数，c表示加速系数，表示最优个体第k位基因，表示种群中第k位基因的最大值。

步骤12：基于精英个体的局部搜索。种群精英个体与全局最优解之间的亲和度要大于种群中的其他个体，利用精英个体信息产生多个与精英个体有较大亲和度的个体，可以提高种群优势个体的多样性，增加算法搜索到全局最优解的概率。因此在精英个体附近进行随机扰动，产生一定数量的新个体，提高算法后期局部搜索能力。新个体产生公式如下：

$>X_k^{\prime }=X_k^{Best}\pm (U_{max}^k-U_{min}^k)/r---\left(14\right)$>

式中：X′_k表示新个体第k位基因，表示最优个体第k位基因，表示种群中第k位基因的最大值，表示种群中第k位基因的最小值，r表示[0，1]内的均匀分布随机数。

步骤13重复步骤7至步骤12，直至算法运行代数达到最大迭代次数。

基于精英个体策略的改进型遗传算法(EISGA)的改进效果见表4和附图2至附图4。表4为基于精英个体策略的改进型遗传算法与其他遗传算法在f₁-f₅五个测试函数上的平均进化代数及收敛率的比较结果；附图2为EISGA与蜜蜂进化型遗传算法在有关过程中种群的分布情况。试验结果能够发现，基于精英个体策略的改进型遗传算法的收敛速度快且精度更高，表明变异算子及精英个体局部搜索等策略增强了算法的寻优能力；进化过程的种群分布情况表明了EISGA具有更强的后期在最优解処的局部搜索能力，使算法搜索到全局最优解的概率增大。

f₁＝100(x₁²-x₂)²+(1-x₁)²，x₁，x₂∈[-5，5]

f₂＝(-1)(x²+2y²-0.3cos(3x)-0.4cos(3y))+4，x，y∈[-1，1]

$>f_3=\frac{x_1^2+x_2^2}{2}-cos\left(2x_1\right)cos\left(2x_2\right),x_i\in [-10,10]$>

f₄＝|(1-x)x²sin(200x)|，x∈[0，1]

$>f_5=0.5-\frac{{\sin }^2\sqrt{x_1^2+x_2^2}-0.5}{{[1-0.001(x_1^2+x_2^2)]}^2},x_i\in [-2,2]$>

表4各测试函数的平均进化代数及其收敛率

约束性知识与精英个体策略遗传算法融合的优化方法在作物生育期模型参数优化中的各项改进如表5、附图3和附图4；表5为分别采用关键约束目性物候期的估算日期和估算误差范围构造适应度函数的优化结果对比；图3-1为作物模型品种范围约束对参数PVT分布结果的影响，图3-2为关键约束性物候期实测数据约束对品种参数IE分布结果的影响。图3-3为徐州和济宁两地试验中，方向因子对收敛速率的影响。图4为本方法在徐州、济宁、潍坊、赣榆四个地点模拟效果的1∶1关系图；试验结果能够发现，在作物模型品种参数优化中，在加入PVT品种参数范围约束后，能够有效的约束PVT的优化结果均符合半冬性的品种特性；加入关键约束性物候期的实测数据约束后，品种参数IE的分布也有明显约束，并且采用关键物候期估算范围构造适应度函数的优化效果较好，表明估算关键物候期的实测数据范围对品种参数的分布具有明显的约束作用，并且能够减小经验估算的误差；加入方向因子后，算法的寻优速度有所提高，表明了针对该应用问题方向算子是有效的。算法优化结果在四个地点均取得了良好的优化效果，证明了本方法的有效性。

表5分别采用估算的关键物候期的物候日期和物候范围构造适应度函数的对比结果