种子园是指以生产优良品质林木种子为目的，由经过筛选的优良无性系或家系按营建设计要求建立的特种人工林(杨培华等, 2005)。无性系配置是种子园营建过程中的重要内容，其目的是实现无性系间的随机交配，保证种子园子代具有较高的遗传增益，同时维持一定的多样性。

目前，对高世代种子园营建的研究分为种子园亲本及园址选择、种子园设计和种子园管理(王章荣, 2012)。其中，在种子园设计的研究中，国内学者研究不多，通常采用随机设计和系统设计(Hodge et al., 1993)，如采用不完全、不平衡固定区组设计种子园(袁虎威等, 2016)，采用无性系顺序错位排列的固定配置法进行设计(程祥等, 2016)，且并未考虑栽植时亲本间的近交对子代的影响。国外在此方面研究较为深入，将种子园的配置设计问题划分为二次分配问题(quadratic assignment problem, QAP)，提出了完全随机区组设计(White et al., 2007)、最小近交设计(minimum inbreeding design, MI)(可选型交配，减少近亲繁殖)(LstibůRek et al., 2010; 2015)、随机复制无性系行交叉法(randomized, replicated, staggered clonal-row, R2SCR)(El-Kassaby et al., 2014)及最优邻域算法(Chaloupková et al., 2016)。国外的研究主要集中于如何设计一种模型以适应各种需要的初代种子园无性系配置设计，通过评价设计方案中各亲本间最小位置距离之和来判断配置设计方案的优劣；国内在近20年内，采用计算机手段进行种子园无性系配置设计的研究工作有限(梁一池, 1991; 1995;贾乃光, 1985; 许鲁平, 2002)。综上，国内外的研究均未基于亲本间遗传距离来考虑近交对种子品质的影响。为提高目标性状的遗传增益，扩大种子园子代遗传多样性，本文依据亲本间遗传距离数据利用计算机手段对种子园无性系进行配置设计。

通常情况下，使用计算机方法基于亲本间遗传距离的种子园无性系配置设计可采用穷举法、贪心法、动态规划等传统方法进行解决，但传统方法仅适合亲本规模小的情况，当亲本规模增大时，上述方法存在成本高、代价高以及实际效用价值低等缺陷。基于遗传距离的种子园设计是QAP问题，同样也是一个NP-hard(non-deterministic polynomial hard)问题，该类问题多使用优化算法解决。常用优化算法有模拟退火算法(simulated annealing，SA)、遗传算法(genetic algorithm，GA)、粒子群算法(particle swarm optimization，PSO)等。优化算法在森林培育领域内的应用研究主要集中在造林规划方面(林晗等，2010; 吴承祯等，1997；2000)，即从林地资源优化利用的角度出发，综合各约束条件以得到最佳经济效益；在林木育种领域，仅有LstibůRek等(2015)采用遗传-禁忌搜索方法对最小近交设计中搜索过程进行研究，国内尚未见到公开发表文献。

本文将利用提出的改进型自适应并行遗传算法(improved adaptive parallel genetic algorithm，IAPGA)，基于亲本间遗传距离进行种子园无性系配置优化设计。目标是在有限时间内，得到近似最优的配置方案，以减小近交、最大化子代适应性的遗传基础，保持遗传多样性，为高世代种子园的无性系配置提供参考。

1 种子园设计问题 1.1 问题描述

为增加拓宽子代遗传基础的概率，可在规划种子园无性系配置方案时，依据遗传距离，平衡地设计栽植方案。因此，可将最大化子代适应性的遗传基础的目标转化为，求解亲本间的相似度之和的最小值。为计算需要，暂忽略样地的地形和地貌特征，假设种子园的形状为规则的。若种子园规模为M*N，即M行N列，则种子园中各栽植位置的编号如图 1所示。

1.2 问题抽象 1.2.1 编码方式

采用Grefenstette编码，将种子园的栽植设计方案转变为基于位置关系的序列，供函数计算使用。使用该编码方式可避免常规的交叉运算和变异运算产生的无意义的设计方案。

1.2.2 目标函数

该问题的目标函数如式(1)所示。

$F = \min \sum f = \min \sum\limits_{i = 1}^{M*N} {f\left({{x_i}} \right)} ; $

(1)

式(1)的约束条件如式(2)所示：

$ i = \left({1, 2, 3, \cdots M \times N} \right), \;\;\;\;{x_i} \to i, \;\;\;\;x \in {{ \mathit{ \boldsymbol{Z}}}^ + }。$

(2)

式(1)中，f为某亲本在其邻域范围内的相似度之和，F为所有亲本在其邻域范围内的相似度之和的最小值；式(2)中，i的取值范围为从1到M*N，代表各个位置，x_i→i代表各位置的亲本服从i的变化，x的取值为正整数。本文所指亲本的邻域范围为与该亲本距离最近一层的植株。如图 1，对图中任一点，需考虑其周围8个位置上的亲本对该点的影响，这8个位置分为2种：一种是其上下左右对应的位置，依次为x-N、x+N、x-1、x+1；另一种为4个斜对角的位置点，依次为x-N-1、x-N+1、x+N-1、x+N+1。若假设种子园的设计形状为规则的，可将直接相连的位置点的距离设为1，斜对角相连的位置点的距离设为$\sqrt 2$。由于距离越远的位置上的亲本对该点的影响越小，距离越近的亲本则对该点的影响越大，因此在求解4个斜角位置上亲本对该点的影响时，需除以$\sqrt 2$以消除距离的影响。因此，每个点的相似度求解如式(3)所示。

$ f\left({{x_i}} \right) = \sum\limits_j {S\left({{x_i}, {x_j}} \right) + \sum\limits_k {S\left({{x_i}, {x_k}} \right) \div \sqrt 2 } } ; $

(3)

式(3)的约束条件如式(4)所示：

$ \begin{array}{l} \;\;\;\;\;\;\;i = \left({1, 2, 3, \cdots M \times N} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_i} \to i, \;\;\;\;x \in {{ \mathit {\boldsymbol{Z}}}^ + }\\ \;\;\;\;\;\;\;\;\;\;\;j = \left\{ {i + 1, i - 1, i + N, i - N} \right\}\\ k = \left\{ {i + N + 1, i + N - 1, i - N + 1, I - n - 1} \right\}。\end{array} $

(4)

式中：S函数为相似度函数，x_j为与x_i直接相连的点，x_k为与x_i对角相连的点，j和k的取值范围如式(4)所示。

以上目标函数仅适用于非边界亲本求解相似度，为使边界点的求解方式同样适用该适应度函数，将原有M行N列改为M+2行和N+2列，即在原有的种子园的边界外增加1圈值为0的点，这些新添加的值为0的点同所有亲本的相似度值为0，这样方便边界点的计算，如图 2所示。

2 材料与方法 2.1 材料及数据来源

研究材料选自内蒙古红花尔基樟子松国家良种基地中初级种子园、1.5代种子园、2代种子园樟子松无性系材料，于5月份摘取当年生针叶3~5针，装管标记冷冻保存。对研究材料采用SSR分子标记法(Morgante et al., 1993; 方攀等，2013)，使用Gene Maker V2.2软件读取针叶的毛细管检测结果，导出引物条带大小等数据信息并整理成excel文件；之后，采用PowerMaker V3.25软件计算基于等位基因频率的Nei 1983距离(Nei et al., 1983)，得到亲本间遗传距离矩阵表，对角线左下方数据为遗传距离，即本文的试验数据。

2.2 研究方法

遗传算法(GA)(Holland, 1975)中，遗传算子中的交叉、变异概率等控制参数事先确定，与进化过程无关，这不利于个体的全局搜索和局部搜索，易早熟。本文提出的改进型自适应并行遗传算法(IAPGA)是在GA的基础上，结合文献(Srinivas et al., 1994; 吴浩扬等, 1999; 蒋建林等, 2012)的研究方法及本文的研究目标进行相应的调整。Srinivas等(1994)首次提出自适应遗传算法(adaptive genetic algorithm, AGA)，使交叉、变异概率可随种群的适应度变化而动态生成；吴浩扬等(1999)对AGA中的种群平均适应度值重新定义，将种群中个体适应度值优于种群平均适应度值的个体的平均适应度值作为新的种群平均适应度值；蒋建林等(2012)提出了启发式多点交叉，在减少重复个体的生成方面有较好表现，促进对解空间的搜索。本文首先在AGA的基础上，结合蒋建林等(2012)的方法提出了改进型自适应遗传算法(improved adaptive genetic algorithm, IAGA)，使交叉、变异概率可自适应调整；之后利用遗传算法自身的并行性，结合计算机的高速并行性，提出IAPGA。相较于GA，IAPGA扩大了对解空间搜索，减少人为控制参数的影响，可根据每代种群的自身特点自适应进行调整。本文提及到的“遗传算子”、“种群”、“子群体”等概念均属于遗传算法中的术语。

2.2.1 IAPGA算法流程

1) 初始化种群，即种子园设计方案群体, 进行相关参数设置:子种群数目、迭代次数、子群体迭代次数；

2) 划分父代种群为几个子群体；

3) 各子群体分别进行选择、交叉、变异操作产生子代，即采用遗传算子生成新的种子园设计方案；

4) 判断各子群体是否达到最大子群体迭代次数，若无执行2)3)，若达到执行5)；

5) 合并各子群体所产生的子代及其父代，计算目标函数，即适应度值，采用精英策略选择适应度值最好的部分个体以及小部分较差个体，组成新的子代，并记录保存；

6) 判断是否达到最大迭代次数，若无返回2)顺序执行，若达到则输出当前适应度值最优个体。

2.2.2 遗传算子

1) 选择算子：选择是从当前种群中选择适应度值低的个体遗传到下一代，通过选择操作可以将优良基因遗传到下一代种群中。进行选择时，采用轮盘赌算法，依据适应度值比例进行选择。

2) 交叉、变异算子：交叉操作将已选择个体配对后进行基因交换，产生新个体。变异操作可以避免种群中的个体过早收敛，维持种群多样性。

3) 个体精英保存策略：使用精英策略可以保持优势个体，保证算法的收敛性。本文对每一代个体的产生都采用精英策略，合并各子群体产生的子代及其父代，将使用精英策略保存的个体作为下一代的父代。未添加精英策略的并行遗传算法虽然也可以得到较优的解，但是每一代解的波动性较大，易丢失最优解。

2.2.3 交叉、变异概率的自适应调整

1) 交叉概率：在AGA(Srinivas et al., 1994)中，首次提出自适应交叉概率计算公式，但需要预先设定参数，为减少人为设定参数对种群进化的影响，引入蒋建林等(2012)提出的启发式多点交叉，如式(5)，提出新的交叉概率求解方法，见式(6)。该方法可根据每代种群及不同个体的适应度值动态生成交叉概率。

$ P_{\rm{c}}^{\rm{J}} = 0.5 + \frac{{{F_{\rm{m}}} - {F_{\rm{n}}}}}{{{F_{\rm{m}}}{\rm{ + }}{F_{\rm{n}}}}}; $

(5)

$ {P_{\rm{c}}} = \left\{ \begin{array}{l} P_{\rm{c}}^{\rm{J}} - \frac{{P_{\rm{c}}^{\rm{J}}\left({F' - {F_{\min }}} \right)}}{{{F_{\min }} + {F_{{\rm{avg}}}} + 1}}, \;\;\;\;F' \le {F_{{\rm{avg\_new}}}}\\ P_{\rm{c}}^{\rm{J}} + \frac{{P_{\rm{c}}^{\rm{J}}\left({F' - {F_{\min }}} \right)}}{{{F_{{\rm{avg\_new}}}} - {F_{\min }}}}, \;\;\;\;F' > {F_{{\rm{avg\_new}}}} \end{array} \right.。$

(6)

式中，F′是参与交叉操作的2个父代个体的适应度值的最大值；F_min是每代种群中个体的最小适应度值；F_avg是每代种群的平均适应度值，F_{avg_new}是每代种群中个体适应度值优于F_avg的个体的平均适应度值；F_m、F_n是进行交叉操作的父代的适应度值。

2) 变异概率：AGA(Srinivas et al., 1994)中变异概率的计算公式同样需要预先设定参数，为减少参数对种群进化的影响，提出新的变异概率求解公式，如式(7)。

$ {P_{\rm{m}}} = \left\{ \begin{array}{l} \frac{{\left({F' - {F_{\min }}} \right)}}{{{F_{\min }} + {F_{{\rm{avg}}}} + 1}}, \;\;\;\;\;\;F \le {F_{{\rm{avg\_new}}}}\\ 1, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;F > {F_{{\rm{avg\_new}}}} \end{array} \right.。$

(7)

式中，F是变异操作中被选中的个体的适应度值；F_min是每代种群中个体的最小适应度值；F_avg是每代种群的平均适应度值，F_{avg_new}是每代种群中个体适应度值优于F_avg的个体的平均适应度值。

3 结果与分析 3.1 试验准备 3.1.1 数据准备

基于已有的樟子松亲本间的遗传距离矩阵，取其中24个亲本间的遗传距离数据作为试验数据。

3.1.2 参数设置

在本试验中，将种子园规模设为4行6列，据此优化种子园的无性系配置设计方案并比较。初始种群规模设为100；在GA中，交叉概率取0.8，变异概率取0.2；IAGA中，无需设定参数；在IAPGA中，划分4个子种群，各子群体迭代代数为20；在AGA_WU(吴浩扬等, 1999)中，k₁、k₂分别取0.3、0.8；在AGA_Srivinas(Srinivas et al., 1994)中，依据参考文献k₁、k₃取1，k₂、k₄取0.5；在SA中，初始温度为1，终止温度为0.001，退火系数取0.996 6，变异概率取0.8；在PSO中，参数c₁、c₂取2。

计算环境：windows10 64bit操作系统，8G内存，3.4 GHz Core i7处理器，python2.7。

3.2 各算法综合对比 3.2.1 收敛状态下各算法对比

预先使用穷举法对种子园无性系配置设计方案的结果进行寻优，按照4行6列的规模，得到最优适应度值为65.292 6。当亲本的数量越大时，使用穷举法的时间开销会呈指数级增长，无法满足实际需求。

本文所研究问题为全局优化问题，使用优化算法进行寻优所得结果通常为最小值或最小值的近似解，即与最小值存在一定误差的解。表 1为各算法分别执行200次后，达到收敛时得到适应度值的对比结果，此处所指收敛状态分为2种，一种为算法适应度值达到最优值，另一情况为连续3 000代适应度值未发生变化。

由表 1可知，本文提出的IAPGA的最小误差值为0.007 4，最大误差值为1.067 4，平均值和最优值之间的误差为0.541 4；相较于传统的GA、AGA_Srivinas、AGA_WU、SA和PSO，IAPGA的最优适应度值的平均值更小，并且方差也最小。虽然最小值出现在Srivinas提出的算法中，但是在试验中仅出现1次，具有偶然性，且该方法的方差稍大于IAPGA的方差，IAPGA相对更稳定。表 1所展示方法的平均适应度值均优于顺序排列的方法，在试验测试过程中，将亲本随机排列栽植得到的适应度值通常大于70，采用优化算法对种子园设计方案进行寻优在很大程度上减少了近交的影响。

3.2.2 收敛速度比较

设定总的迭代次数为5 000次，每隔50次输出结果，对比展示各算法的收敛速度。图 3展示某随机6次运行结果。从图中可以看出本文提出的IAPGA的收敛态势和效果均优于GA、SA。

3.3 种子园无性系配置 3.3.1 不同方法设计对比

图 4直观地展示了在种子园无性系配置中，各位置上的亲本受近邻亲本影响的程度，颜色越浅越趋近于蓝色，代表处于该位置的亲本受近交影响越小，越有利于扩大子代的遗传基础。其中，(A)图为亲本随机排列的种子园设计方案，(B)图为亲本按原顺序排列的种子园设计方案，(C)图为使用传统计算方法，即贪心算法得到的种子园无性系配置设计方案，(D)图为使用IAPGA优化的种子园无性系配置设计方案。可观察到，使用本文方法设计的种子园，各位置的亲本受近交的影响可达到相对最小，更有利于保持遗传多样性；在随机排列、顺序排列中，单个位置亲本的最大相似度之和达到4.5，仅有位于4角的亲本所受影响最小，缺乏对内部亲本近交影响的考虑；贪心算法得到的排列方案虽优于顺序排列法，但其单点位置的相似度最大值仍大于IAPGA方法，IAPGA在单点位置上最大相似度不超过3.5；且从整体上来看，使用IAPGA时，处于非边缘位置亲本受近邻亲本的影响更均匀，比原有的设计方法更平稳。

3.3.2 不同行列数设计对比

图 5为使用IAPGA优化设计的种子园无性系配置方案，对比展示不同排列结构时，各位置亲本受近交影响程度。从图 5可以看出，按照2行12列规模得到的种子园方案中各位置亲本受近交的影响最小，3行8列其次，处于非边缘位置的亲本的数量越多，受近交影响的亲本的数量也越多；在实际应用中，最优解即最佳种子园无性系配置方案与行列数的设定有关，其具有在该行列数下解的唯一性，与采用何种方法进行种子园的配置设计无关。因此，在进行种子园设计时，预先结合基地中样地的规模进行种子园的合理规划设计，可以有效减少部分亲本间的近交。

4 结论

种子的遗传多样性受亲缘关系影响，当已知各亲本间遗传距离，可利用该数据合理规划种子园无性系配置，以减少近亲繁殖的影响，扩宽子代的遗传基础，保持遗传多样性。本文提出的改进型自适应并行遗传算法(IAPGA)可在初始种群规模较小的情况下，较原有的遗传算法和模拟退火算法得到更优解；通过利用交叉、变异概率的自适应生成，减少种群中优秀个体参与进化过程的概率，增加较差个体的进化概率，保持优势个体的稳定性，并结合计算机的并行能力，扩大对解空间的搜索。计算结果表明，在较短时间内，使用IAPGA得到的种子园无性系配置方案可使遗传距离较近的无性系在配置上保持最大距离，减少近亲交配机会，在一定程度上扩大了子代的遗传基础，最终实现基于遗传距离的种子园无性系优化配置设计，既满足了目标要求也满足了时间需求。

林木种子园是一个复杂的交配系统，本研究主要解决了种子园无性系配置设计中的近亲交配问题，改进型自适应并行遗传算法仅将亲本间的遗传距离作为唯一的约束条件。下一步的研究工作将围绕种子园不规则地块、无性系多少、无性系分株数，以及亲本数量无法按照相应规模拆分等情况开展；之后将结合不同亲本的花期、花粉、育种值等约束条件修改目标函数，进一步优化种子园设计方案，为高世代种子园建设提供参考。