本发明涉及一种基于双目标免疫粒子群算法，在粒子群算法中加入拥挤距离和免疫进行改进解决干旱地区缺水情况下农作物生育期水资源分配问题，属农业灌溉领域。
背景技术：
我国是一个干旱缺水严重的国家，由于我国人口的不断增长，工业化发展的不断扩大，加之因为人类活动造成的环境污染，水资源越来越匮乏。我国是一个农业大国，总用水量中农业用水量占了绝大部分，但是即使是这样，我国还有许多农田没有足够的水进行灌溉，水资源的匮乏使得有些地区的作物无法按照充分需水量进行灌溉，进而影响产量，因此研究非充分灌溉，实现农业水资源的优化配置至关重要。
灌溉水资源优化配置主要研究在水量有限的情况下，根据作物的生育期特点，在每个生育期进行合理的水量灌溉，使作物的产量达到最大的经济效应。静态模型描述了作物最终产量与水分的宏观关系，直接把作物产量表示为作物生育期内腾发量的函数，可通过对试验数据进行回归分析来确定。静态模型分为全生育期水分生产函数模型和分阶段水分生产函数模型，本发明中的静态模型选用jensen模型，即
式中，etmij为作物i在其第j生育阶段充分灌溉下的最大需水量，单位为mm；etij为作物i在其第j生育阶段充分灌溉下的实际需水量，单位为mm；m为作物i的剩余阶段数；λij为作物i在其第j生育阶段的缺水敏感指数；tmi为作物i的最大产量，单位为kg/hm2；yi为作物i的实际产量，单位为kg/hm2。
粒子群算法(particleswarmoptimization，pso)是一种群体智能的优化算法。该算法于1995年由kennedy和eberhart提出，pso算法源于对鸟类捕食行为的研究。当鸟类进行捕食时，会搜寻当前距离食物最近的鸟的周围区域，然后快速的靠近食物，这样能够有效简单地找到实物。pso算法就是从这种生物种群行为特征中得到启发并用于求解优化问题。算法中每个粒子都代表问题的一个可能解，每个粒子经过适应度函数的计算都有一个对应的适应度值，粒子速度更新决定了粒子移动的方向和距离，而速度本身跟自己的全局最优位置和个体最优位置有关，根据这两者进行调整，能够在不断的迭代更新中找到适应于问题的最优解。
粒子群算法通过追寻个体极值和全局极值来达到整体粒子最优，收敛速度较快，但是在迭代次数不断增加的情况下，各粒子相似度越来越高，可能陷入局部最优解而无法跳出，优化的结果并不十分理想，本章加入拥挤距离和免疫的概念，提出一种改进的粒子群算法，避免了粒子过度聚集在一起，且优化的非劣解集更好。
技术实现要素：
本发明的目的是：解决干旱地区农作物缺水条件下农业水资源的合理分配，以到达农作物产量最大而灌溉水量最小。
为了达到上述目的，本发明的技术方案是提供了一种基于双目标免疫粒子群算法的灌溉水资源优化配置方法，其特征在于，使用基于拥挤距离的免疫粒子群算法找寻使得农作物产量达到最大而灌溉水量达到最小的一组解集，包括以下步骤：
步骤1、分析农作物的生育周期，将农作物的整个生长过程划分为不同的生育阶段；
步骤2、查找当地农作物每个生育阶段所对应的降雨量和作物的最大需水量；
步骤3、将灌溉-产量静态模型描述为：
式中，yi为作物i的实际产量；wij为作物i在其第j生育阶段的灌水量，wij＝etij-peij，etij为作物i在其第j生育阶段的实际需水量，peij为作物i在其第j生育阶段的有效降雨量；etmij为作物i在其第j生育阶段充分灌溉下的最大需水量；tmi为作物i的最大产量；
步骤4、设置初始的灌溉水量，进行粒子群算法的初始化；
步骤5、粒子群算法随机生成一代粒子之后，根据拥挤距离进行粒子的排序，选择拥挤距离较大的粒子放入抗体库；
步骤6、进行粒子的速度和位置的更新操作；
步骤7、在抗体库中选择一个粒子作为抗体，对新生成的粒子进行交叉和变异；
步骤8、生成步骤3所得的灌溉-产量静态模型的作物产量最大而灌溉水量最小的非劣解集并更新抗体库；
步骤9、判断是否达到最大代数，如果没有，继续重复步骤4～7；若达到最大代数，则输出得到的作物产量最大而灌溉水量最小的非劣解集。
优选地，在所述步骤1中，若农作物在某一生育阶段的降雨量多于最大需水量，则舍弃该生育阶段。
优选地，在所述步骤4中，粒子群算法的初始化包括以下步骤：
将每个生育阶段的灌溉量作为一个5维的粒子，每一维都是实数编码，在初始化粒子的过程中，根据当前生育阶段的降雨量和作物的蒸腾量设置当前生育阶段的灌溉水量约束，并进行总体水量的约束，然后随机生成种群数大小的粒子群。
优选地，步骤5中，拥挤距离的定义为：
拥挤距离是用来计算某个体与其他个体之间的距离，用以表征个体间的拥挤程度，拥挤距离越大，表明个体间就越不拥挤，种群的多样性就越好；
拥挤距离的计算方法为：
对于多目标中的每一个目标，将所有种群个体的值先除以(种群个体中最大的值+1)的方法压缩到(-1，1)的区间上，然后进行升序排序，设定一个拥挤距离的集合，将排序后的个体中最大和最小值对应于拥挤距离集合中的位置处设为inf，其他个体对应位置的拥挤距离为该个体相邻的前后两个个体在(-1，1)区间映射后的目标函数值之差，最后的拥挤距离为所有目标的拥挤距离之和。
优选地，所述步骤6包括以下步骤：
在一个d维搜索空间中，由n个粒子组成的种群x＝(x1，x2，...，xn)，其中第i个粒子表示一个d维的向量xi＝(xi1，xi2，...，xid)t，代表第i个粒子在d维搜索空间中的位置，亦代表问题的一个潜在解；根据目标函数计算出每个粒子赌赢的适应度值，而对应的第i个粒子的速度为vi＝(vi1，vi2，...，vid)t；其个体最优为pi＝(pi1，pi2，...，pid)t；种群的全局最优为pg＝(pg1，pg2，...，pgd)t；在每次的迭代更新过程中，粒子的更新速度表示为：
ω(k)＝ωstart-(ωstart-ωend)*(k/tmax)2
vidk+1＝ωvidk+c1r1(pidk-xidk)+c2r2(pgdk-xidk)
xidk+1＝xidk+vidk+1
式中，ω(k)为第k次迭代的惯性权重；ωstart为初始惯性权重；ωend为迭代至最大次数时的惯性权重；tmax为最大迭代代数；vidk+1为第k+1次迭代的第i个粒子的速度，i＝1，2，...，n；vidk第k次迭代的第i个粒子的速度；c1和c2为加速度因子；r1和r2是分布于[0，1]区间的随机数；pidk第k次迭代的第i个粒子的个体最优；xidk第k次迭代的第i个粒子的位置；pgdk第k次迭代的种群的全局最优；xidk+1为第k+1次迭代的第i个粒子的位置。
优选地，步骤7中，对粒子群中的所有粒子进行交叉的过程为：
依照交叉的概率在种群中随机选择一个个体，然后从抗体群中拥挤距离大的前5个中选择一个个体，随机选择5位实数中的两位，对应着将抗体群中的抗体同等位上的值赋给粒子种群中选择的个体中的同一位；
进行变异的过程为：
先遍历一遍粒子，如果发现双目标都劣于抗体群中选择的抗体粒子的双目标，就对当前粒子进行全部位的变异，如果原有粒子跟抗体粒子是非劣解的关系，则依照概率进行随机位变异。
优选地，在步骤8中：
非劣解集就是所有非劣解的集合，而非劣解就是当一个粒子不受其他粒子支配时，其就为非劣解，如果个体p至少有一个目标比个体q好，而且个体p的所有目标都不比个体q的差，那么称个体p支配个体q，或者称个体q受个体p支配；
更新抗体库的方法为：
计算当前粒子群的适应度并和抗体库中的粒子进行比较，如果现有的支配抗体群中的粒子，则替换抗体群中的粒子；如果两者为非支配关系，查看此时的抗体库中粒子的个数，如果个数为r，r＜q，而新的非劣解集的个数假设为p，如果p＜q-r，则将p全部根据拥挤距离排序后放入抗体群中，如果p＞q-r，则在p根据拥挤距离排序后选取q-r个拥挤距离大的粒子放入抗体库中；另一方面，如果此时的抗体库中粒子的个数就是q，p根据拥挤距离排序后，如果其个数小于5个，则替换抗体库中后面p个个体，如果其个数大于5个，则选取p中拥挤距离大的5个去替换抗体群中最后的5个个体；
同时将当前的非劣解集替换之前的非劣解集。
本文发明提供的一种基于双目标免疫粒子群算法的灌溉水资源优化配置方法，根据拥挤距离来选择合适的灌溉水量粒子，同时向最好的灌溉水量粒子进行交叉和变异，从而在干旱地区的农作物生长中进行合理的水资源配置，以达到使得作物产量最大而灌溉水量最小的目的。
本发明具有如下优点：
粒子群算法中加入了拥挤距离，改变了粒子群算法容易陷入局部最优解的问题，同时引入了免疫中的交叉和变异，使种群的多样性提高，结果显示在干旱地区缺水条件下冬小麦5个生育期的水资源合理分配能够使得冬小麦的产量达到最大的同时灌溉水量最小。
附图说明
图1是本发明的流程图；
图2是本发明的改进的粒子群算法搜索最优解的方法流程图；
图3是本发明的农业灌溉水量的最优水配置情况下灌水总量和产量的关系结果；
图4是灌溉水量总量和单位水量的产量之间的关系图。
具体实施方式
下面结合具体实施方式，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。
采用matlab对干旱地区灌溉水资源优化配置进行仿真，仿真部分流程图如附图2，具体步骤如下：
(1)分析农作物的生育周期，将农作物的整个生长过程划分为不同的生育阶段。选取冬小麦为研究对象，则方法研究冬小麦的生长周期被划分为播种～越冬、越冬～返青、返青～拔节、拔节～抽穗、抽穗～灌浆和灌浆～成熟。由于越冬～返青季节降雨量多于作物的蒸腾量，因此舍弃该阶段。
(2)查找当地农作物每个生育阶段所对应的降雨量和作物的最大需水量。
资料选取自山西省霍泉县试验区小麦灌溉试验，其1998年小麦的作物敏感系数、作物生育阶段最大需水量和降雨量的资料如表1：
表1山西省霍泉县1998年小麦作物敏感系数、作物阶段最大需水量和降雨量
(3)选择灌溉-产量静态模型
在实际灌溉中，有很多因素影响作物的产量，包括风速、空气湿度和土壤湿度等等，但是在用投入产出关系来确定农作物的产量时，如果仅仅考虑灌溉对产量的影响，则灌溉-产量模型可以描述为
yi＝f(wij)
式中，yi为作物i的实际产量，单位为kg/hm2；wij为作物i在其第j个生育阶段的灌水量，单位为mm。
在缺水的条件下，可以忽略土壤水的水分交换、地表径流等因素，土壤水量平衡公式可简化成
wij＝etij-peij
式中，etij为作物i在其第j生育阶段的实际需水量，单位为mm；peij为作物i在其第j生育阶段的有效降雨量，在非充分灌溉条件下，假设降雨量等于有效降雨量。
作物各个生育期的灌水量都直接影响作物的整体产量，而且内部的前后影响较为复杂，但是仍然可以使用许多数学模型的结构关系表征作物不同生育阶段灌溉水量对产量的相互影响，应用比较普及的jensen乘法模型
式中etmij为作物i在其第j生育阶段充分灌溉下的最大需水量，单位为mm；m为作物i的剩余阶段数；λij为作物i在其第j生育阶段的缺水敏感指数；ymi为作物i的最大产量，kg/hm2。
综合上述表达式，得到作物灌溉-产量模型为
本实施例中，参数设置如表1，则双目标问题为
ym为冬小麦1999年的实际产量，wi为冬小麦第i阶段的灌水量，pei为冬小麦第i阶段的降雨量，m为冬小麦的5个生育期，etmi为冬小麦第i阶段的最大需水量。
(4)设置初始的灌溉水量，进行粒子群算法的初始化。
将每个阶段的灌溉量作为一个5维的粒子，每一维都是实数编码，在初始化粒子的过程中，根据该阶段的降雨量和作物的最大需水量设置该阶段的灌溉水量约束，具体方法为作物的最大需水量减去该阶段的降雨量并进行合理的放大，则每个阶段的灌溉水量约束条件为
即第一阶段的水量约束条件在0～20，第二阶段水量约束条件为0～45，后面的阶段以此类推。
在matlab设置算法的种群个数为100，迭代次数为200，c1为1.2，c2为0.8，ωstart为1.2，ωend为0.9，维数为5维，交叉概率为0.75，变异概率为0.3。然后依据上述约束条件生成100个5维的实数编码数据作为初始粒子。
(5)粒子群算法随机生成一代粒子之后，计算其适应度函数值，然后根据拥挤距离进行粒子的排序，排序的方法具体为：
对于多目标中的每一个目标，将所有种群个体的值先除以(种群个体中最大的值+1)的方法压缩到(-1，1)的区间上，如下公式：
data＝data./(1+max(abs(data)))
式中，data为其中一个目标的所有个体的值的集合。然后进行升序排序：
[sorteddata，index]＝sort(data)
式中，sorteddata为单个目标值排序之后的集合，index为单个目标值排序之后对应的索引集合。
设定一个拥挤距离的集合，将排序后的个体中最大和最小值对应于拥挤距离集合中的位置处设为inf，
crowdingdistance([index(1)，index(end)])＝inf
式中crowdingdistance为拥挤距离的集合，index(1)为单个目标值排序之后的索引集合的第一个值，index(end)为单个目标值排序之后的索引集合的最后一个值。其他个体对应位置的拥挤距离为该个体相邻的前后两个个体在(-1，1)区间映射后的目标函数值之差。
crowdingdistance(index(i))＝data(index(i+1))-data(index(i-1))；
i＝2，3，4，...，xsize，xsize为最大代数，index(i)为单个目标值排序之后的索引集合的第i个值。最后的拥挤距离为所有目标的拥挤距离之和。
选择拥挤距离较大的粒子放入抗体库。抗体库的数量设置成q个，在排序过后的非劣解集中根据轮盘赌的方法选择前5个拥挤距离较大的个体作为全局最优粒子。非劣解进行排序后，如果其个数小于抗体库的个数q，则全部放入抗体库，如果其个数大于抗体库的个数q，则选取q个拥挤距离较大的值放入抗体群。
(6)进行粒子的速度和位置的更新操作。
粒子群算法的数学表示为在一个d维搜索空间中，由n个粒子组成的种群x＝(x1，x2，...，xn)，其中第i个粒子表示一个d维的向量xi＝(xi1，xi2，...，xid)t，代表第i个粒子在d维搜索空间中的位置，亦代表问题的一个潜在解。根据目标函数可以计算出每个粒子赌赢的适应度值，而对应的第i个粒子的速度为vi＝(vi1，vi2，...，vid)t；其个体最优为pi＝(pi1，pi2，...，pid)t；种群的全局最优为pg＝(pg1，pg2，...，pgd)t；
更新的公式为：
ω(k)＝ωstart-(ωstart-ωend)*(k/tmax)2
vidk+1＝ωvidk+c1r1(pidk-xidk)+c2r2(pgdk-xidk)
xidk+1＝xidk+vidk+1
其中，ωstart为初始惯性权重，ωend为迭代至最大次数时的惯性权重，k为当前迭代代数，tmax为最大迭代代数，i＝1，2，...，n：d＝1，2，...，d；ω为惯性权重；vid为粒子的速度；c1和c2是非负的常数，称为加速度因子；r1和r2是分布于[0，1]区间的随机数。pidk为第k代的全局极值，pgdk为第k代的个体极值。
(7)在抗体群中选择一个粒子作为抗体，对新生成的粒子进行交叉和变异。
依照交叉的概率在种群中随机选择一个个体设为chrom1，从抗体库中拥挤距离大的前5个中选择一个个体设为chrom2，随机选择chrom1中的两位，对应着将chrom2的同等位上的值赋给chrom1，其中pos1和pos2为随机生成的5位中的不同的两位。即
chrom1(pos1)＝chrom2(pos1)
chrom1(pos2)＝chrom2(pos2)
变异的过程为：先遍历一遍粒子，如果发现双目标都劣于抗体群中选择的抗体粒子的双目标，就对当前粒子进行全部位的变异，如果原有粒子跟抗体粒子是非劣解的关系，则依照概率进行随机位变异。变异的公式为
u＝chrom1(pos)
u1＝chrom2(pos)
u2＝u1-u
delta＝u2*rand(1)
u＝u+delta
其中chrom1为种群中每一个粒子，chrom2为抗体库中拥挤距离前5个中的随机一个粒子，pos为随机位，u2为原有粒子和抗体粒子同位的差值。delta为u2乘以一个随机数，可以看到u是朝着u1的方向去加上或者减去delta。最后将u的值替换chrom1(pos1)的值即可。
(8)生成一组非劣解集并更新抗体库，具体方法为：
计算当前粒子群的适应度并和抗体库中的粒子进行比较，如果现有的支配抗体库中的粒子，则替换抗体库中的粒子。如果两者为非支配关系，查看此时的抗体库中粒子的个数，如果个数为r，r＜q，而新的非劣解集的个数假设为p，如果p＜q-r，则将p全部根据拥挤距离排序后放入抗体库中，如果p＞q-r，则在p根据拥挤距离排序后选取q-r个拥挤距离大的粒子放入抗体库中。另一方面，如果此时的抗体库中粒子的个数就是q，p根据拥挤距离排序后，如果其个数小于5个，则替换抗体库中后面p个个体，如果其个数大于5个，则选取p中拥挤距离大的5个去替换抗体群中最后的5个个体。同时将当前的非劣解集替换之前的非劣解集。
(9)判断是否达到最大代数，如果没有，继续重复步骤(4)～(7)；若达到最大代数，则输出结果。
当达到最大代数时，输出抗体库中的非劣解集，其能均匀的描述灌溉水量总量和冬小麦产量的关系，结果如图3所示，通过数据的调用，也可以得到每个非劣解对应的5个阶段的灌溉水量。随着灌溉水量的增加，农作物的产量逐渐增加，当灌溉水量达到农作物充分灌溉的情况时，此时再增加水量，产量反而减少，从图4可以看出，水量的增加，单位水量的产量反而在减小，证明了该方法符合冬小麦生长过程中的边际效应。