基于改进型果蝇算法的无性系种子园设计

齐建东 刘春霞 崔晓晖 李 伟

(1.北京林业大学信息学院， 北京 100083； 2.北京林业大学生物科学与技术学院， 北京 100083)

0 引言

当前我国林业发展中，营林造林主要依靠林木育种。其中种子园是林业生产中使用的主要手段，种子园是由优树无性系或家系按营建设计要求，实现集约经营，以生产遗传品质和播种品质优良的林木良种为目的的特种人工林[1]。在种子园设计中，为了避免近交衰退的不利影响，需要尽量减少无性系之间的近交繁殖，以达到种子园遗传效益最大化的目标。

现今，我国大部分种子园已经结束了初级改良工作，正在向高世代育种发展，主要考虑的问题有种子园亲本材料选择、种子园设计、种子园植株管理、育种群体组建及长期育种计划制定等[2]。其中，种子园设计的常用方法有：系统设计(System design，SD)(又称为顺序错位排列)、完全随机设计(Complete random design，CR)、随机完全区组设计(Randomized complete block，RCB)等[3]。国内学者在种子园设计领域的研究较少，代表性的工作有：许鲁平[4]、申文辉等[5]和袁虎威等[6]将无性系材料分为不同区组，分别采用约束变换区组设计、约束分组结合随机完全分组设计和不平衡、不完全固定区组设计的配置方式；程祥等[7]采用无性系顺序错位排列的固定配置法进行设计，以上研究工作均未充分考虑亲本之间的不同亲缘关系以及最大化近交距离。国外学者对种子园设计研究更为深入，BELL等[8]采用重置近邻设计，LSTIBUREK等[9]提出了最小化近交(Minimum-inbreeding, MI)设计方案，EI-KASSABY等[10]提出了随机、重复、交错的无性系行(Randomized, replicated, staggered clonal-row, R2SCR)种子园设计方案，LSTIBUREK等[11]通过遗传禁忌算法以及合并独立的MI算法改进了原有的MI方案，CHALOUPKOVA等[12]提出最优近邻算法(ONA)，以上学者在种子园设计方法上更注重考虑不同的亲缘关系和亲缘比例，研究工作的共同点是：无性系的分株数量以及非平衡设计比例均为已知条件；在固定的无性系种类和无性系分株比例中完成种子园设计。

在现有研究工作中，种子园无性系的亲缘关系、数量和比例都是研究者实验假设的前提条件，如：LSTIBUREK等[9]在实验中假设种子园的无性系只存在半同胞亲缘关系或无亲缘关系等不同亲缘关系的场景，但是在现实种子园中亲缘关系和遗传多样性往往会比假设条件更复杂、更混乱。因此，明确种子园中的亲缘关系，可以有效控制近交，提高子代杂合率。

本文仅基于无性系种类和无性系之间的遗传距离，利用改进型果蝇算法来使种子园中的无性系在满足最小近亲繁殖的同时得到较高的遗传效益，以及种子园的较优配置方案和合理的亲本分株比例。最后将改进的果蝇算法与传统的种子园配置方法(CR、RCB)、遗传算法(GA)以及其他两个改进果蝇算法进行比较，以评估该设计方案的优越性，并为高世代种子园的设计提供参考。

1 种子园设计

1.1 问题描述

假设存在一个规模为M×N(M行、N列)的种子园，需要从已知遗传距离的T株亲本中挑选出合适的亲本比例进行栽种。在栽种过程中，除种子园的四周位置以外，每个中心位置由8个位置包围，这8个位置作为中心位置的近邻位置，分为两种情况：正对近邻和斜角近邻，易知种子园四角位置只存在3个近邻位置，剩下的周边位置都存在5个近邻位置。其中近邻位置作为中心位置计算近邻距离的重点考虑对象，因为亲本之间的遗传距离越小，则说明亲本之间的亲缘关系越近，所以应该选择遗传距离较大的植株作为近邻，更利于提高种子园的遗传效益。除此之外还需要考虑同一亲本的不同分株之间的距离影响，即同一亲本的不同分株应该栽种在适当远的距离，以避免自交现象以及优先交配亲本组合不断循环的现象。本文实验中，不考虑种子园地形、气候等客观因素。

1.2 目标函数

该问题的目标函数为

(i≠j,k,t且j≠k)

(1)

式中fmin——所有植株的近邻距离与同一无性系所有分株距离之和的最小值

Gij——第i株无性系和第j株无性系之间的遗传距离, 第i株无性系和第j株无性系为正对近邻

Gik——第i株无性系和第k株无性系之间的遗传距离，第i株无性系和第k株无性系为斜角近邻

dit——第i株无性系和第t株无性系之间的物理距离，第i株无性系和第t株无性系为同一无性系亲本的分株

1.3 不同目标函数对比

将本文设计的目标函数分别与现有MI算法[9]中的目标函数

(2)

以及改进MI算法[11]中存在半同胞(half-siblings, h-s)关系时的目标函数

(3)

式中dil,jl——第l株无性系的第i株和第j株分株之间的距离

Nc——无性系的数量

NR——同一无性系的分株数量

dmin——所有种类的无性系中同一无性系的任何两个分株的所有可能距离的平方倒数和

dmin(h-s)——所有具有半同胞关系的成对分株距离之和

进行对比。

本文提出的算法以及两个对比算法的最终目标都是为了最小化近亲繁殖，其中，式(2)只将同一无性系的不同分株之间的距离作为近亲繁殖的衡量标准，而同一无性系的不同分株间授粉称为自交，所以式(2)其实只考虑了自交繁殖现象，并没有考虑到不同无性系之间的遗传关系以及它们之间的交配现象。式(3)加入了半同胞情况，但是在现实的种子园中特别是高世代种子园中，无性系之间具有的亲缘关系往往更复杂，算法设计不能单纯只考虑同一无性系之间的自交和半同胞现象。

本文提出的目标函数不仅考虑了同一无性系的不同分株之间的影响，还引入了遗传距离作为不同无性系之间的亲缘关系的衡量标准，利用遗传距离计算种子园中所有不同无性系之间的近交繁殖距离。

2 材料与方法

2.1 材料与数据来源

实验材料是从内蒙古红花尔基樟子松(Pinussylvestrisvar.mongolica)国家良种基地中的初级种子园、1.5代种子园、2代种子园采集的当年生针叶3～5针的樟子松无性系材料，由于SNP分子标记技术成本较高、分型技术不太成熟，因此本实验材料是通过简单、成熟、成本低的SSR分子标记法提取的樟子松基因组DNA，从15对SSR引物中选择其中多态性良好、稳定、清晰的11对SSR引物作为实验所用引物，利用GeneMarker V2.2软件对条带信息进行比对后，再基于等位基因频率的Nei 1983距离计算得到樟子松无性系材料中不同无性系之间的遗传距离[13]，作为本文的实验数据。

2.2 研究方法

种子园设计本质上是一个二次分配问题(Quadratic assignment problem，QAP)，作为一个NP难问题，QAP的解决方法有穷举法和智能算法等。受计算能力和时间制约，穷举法只适用于解决小规模的QAP案例，智能算法是解决大规模QAP案例的可行、有效手段。

较为成熟的智能算法有遗传算法(Genetic algorithms，GA)、禁忌搜索(Tabu search，TS)、粒子群算法(Particle swarm optimization，PSO)等，其中，LSTIBUREK等[9]的MI算法利用了禁忌搜索算法，且LSTIBUREK等[11]在改进的MI算法中采用了遗传禁忌搜索算法，进一步优化了原有的MI算法，使该方案适用于更大更复杂的高世代种子园。虽然经典的智能优化算法已经得到了广泛的应用，但是各自都存在一些明显的不足，这与待解决领域问题相关。近些年新的仿生智能算法也不断出现，例如果蝇优化算法(Fruit fly optimization algorithm，FOA)、蚁群算法(Ant colony optimization，ACO)、鱼群算法(Fish swarm algorithm，FSA)、免疫算法(Immune algorithm，IA)等，相比其他优化算法，FOA具有简单、参数少、易调节、计算量小、寻优精度较高等优点，易于实现，更适合解决结构复杂的现实问题[14-15]。

2.2.1果蝇优化算法

果蝇优化算法主要分为嗅觉觅食环节和视觉觅食环节，然后对这两个觅食环节进行不断的迭代最终实现该果蝇种群的进化，获得该算法对问题的最优解。但是，到目前为止，果蝇算法多用于解决连续型目标优化问题，对于离散型目标函数研究还较少，用该算法解决的离散型目标函数最常见的有旅行商问题(Traveling salesman problem, TSP)[16-18]和调度问题[19-20]，在解决TSP问题中王克甫等[16]引入局部最优概念并采用自适应步长策略抑制了早熟现象并提高了搜索效率，段艳明等[18]将果蝇算法和遗传算法以及C2Opt算子相结合加快了解决问题的局部搜索能力和收敛速度。基于不同学者对果蝇算法的改进方法，本文针对无性系种子园配置问题对传统果蝇优化算法的不足进行了改进，提出改进型果蝇算法(Improved fruit fly optimization algorithm，IFOA)。

2.2.2改进果蝇优化算法

无性系种子园设计问题不同于连续型优化问题，对T株无性系按随机比例栽种到规模为M×N的种子园的配置方案为一只果蝇，配置方案中的任意一个位置为一个基因位，其中，改变果蝇的任意一个基因位，称为一次觅食；在觅食过程中对基因位的选择行为称为觅食的方向选择，每只果蝇的浓度判定函数为式(1)。

为了使FOA能够更好地解决本文的无性系种子园配置问题，对FOA做了以下改进：

(1)轮盘赌法初始化种群。果蝇算法的初始值会在一定程度上影响果蝇算法的收敛速度，本实验采用轮盘赌法初始化果蝇种群，将所选的T株亲本的遗传距离的倒数作为轮盘赌法中的概率p(i)，即

(4)

式中 sum(G)——候选亲本中所有遗传距离的和

G(i,j)——第i株无性系到第j株无性系的遗传距离

(2)过滤无效树种选择。在每次觅食过程中，获取种子园当前基因位的正对近邻的树种，为节约觅食搜索时间，加快搜索速度，过滤与正对近邻相同无性系的选择，避免无效的觅食。

(3)步长选择。FOA的步长会影响算法的效率和精度，当步长较大时有利于全局搜索，加快搜索速度，但精度较低；步长较小时有利于局部搜索，提高算法精度，但搜索速度会下降。传统的FOA是采用固定步长进行觅食，因此很难平衡算法的效率和精度。为兼顾两者平衡，本文引入遗传算法中的基因变异对步长进行调整，在迭代前期可设置较大的变异因子，有利于加快算法的搜索和收敛速度；在迭代后期减小变异因子，在保证算法精度的同时也有助于算法跳出局部最优。

(4)觅食方向。在迭代前期为了加快搜索速度，将觅食方向定为种子园中浓度最差的基因位，在迭代后期为了避免算法陷入局部最优，将觅食方向改为随机方向。

在求解无性系种子园配置问题时需要把FOA的连续解空间对应到种子园配置方案，其主要步骤见图1。

图1 IFOA算法流程图Fig.1 Flow chart of IFOA algorithm

3 实验与结果分析

3.1 实验准备

3.1.1数据材料

从内蒙古红花尔基樟子松国家良种基地中的1代种子园、1.5代种子园、2代种子园遗传距离的Excel文件中，随机选取其中13个无性系和它们的遗传距离作为实验数据。

3.1.2参数设置

本实验中，种子园的规模假设为9×9(即9行、9列)，初始种群大小为25，迭代次数为800，迭代前期变异因子为0.4，迭代后期变异因子为0.02，在此基础上将本实验与传统的种子园配置方法(CR、RCB)、遗传算法(GA)以及其他两个改进果蝇算法：MFFA[16]、ASFOA[18]进行了比较，以上所有算法均

采用式(1)作为目标函数。

3.1.3实验说明

该算法并不只针对数据材料中选取的13个无性系或参数设置中的9×9固定大小的种子园，它适用于任意无性系的种类和数量以及任意规模的种子园。

3.2 算法比较

3.2.1综合对比

通常全局优化算法结果并不是问题的最小值，而是最小值的近似解，即最终的结果存在一定的误差。为了避免偶然事件出现影响对比结果，将各算法在以上参数设置条件下分别执行200次后得到对比结果，见表1。

从表1可以发现，本文提出的IFOA的平均值、最大值、最小值分别为1 169.780、1 198.153、1 137.774，都明显低于其他算法，其中IFOA的最大值也明显低于其他算法的最小值，其次结果较好的为GA。相对于传统算法而言，智能优化算法具有明显优势。

表1 不同算法分别执行200次后适应度对比Tab.1 Comparison of fitness of each algorithm after 200 times respectively

3.2.2收敛速度比较

从200次试验中随机选择6次运行结果，作为各算法的收敛速度，见图2，IFOA的收敛速度以及效果明显优于其他算法。

3.2.3不同方法设计对比

将本文的IFOA算法得出的设计方案，分别与效果较好的智能算法GA设计方案以及传统算法CR设计方案进行对比，如图3所示。

图3a存在多处同一无性系的分株作为直接邻居，如第1行第3个位置和第4个位置都是无性系9的分株；图3b、3c没有同一无性系的分株作为直接邻居出现；图3c无性系3和无性系7作为直接邻居出现的概率明显高于图3b，原因是无性系3和无性系7之间的遗传距离是实验数据中遗传距离最大的组合，因此会优先交配，有利于杂交育种；图3c没有无性系2，因为无性系2与其他无性系之间的遗传距离都较小，即无性系2与其他无性系都具有较近的遗传距离，因此不建议栽种。

图2 各算法收敛情况对比Fig.2 Comparison of convergence of different algorithms

图3 13个无性系在9×9种子园的设计方案Fig.3 Design of 13 clones in 9×9 orchard

4 结束语

将无性系的分株数量作为未知条件，基于无性系之间的遗传距离，不仅考虑了同一亲本的分株影响，也考虑了近邻之间的影响，使得种子园最终能够达到较优的遗传效益，通过改进型果蝇算法实现该目标，并得到较优的种子园设计方案和亲本分株比例。计算结果表明IFOA的收敛速度和效果优于其他算法，并且IFOA的设计方案中同一无性系分株不会作为近邻出现，当存在亲本组合的遗传距离明显大于其他组合时，即亲缘关系最远组合，该组合通常会作为优先交配组合。