基于改進型果蠅算法的無性系種子園設(shè)計

齊建東 劉春霞 崔曉暉 李 偉

(1.北京林業(yè)大學(xué)信息學(xué)院， 北京 100083； 2.北京林業(yè)大學(xué)生物科學(xué)與技術(shù)學(xué)院， 北京 100083)

0 引言

當(dāng)前我國林業(yè)發(fā)展中，營林造林主要依靠林木育種。其中種子園是林業(yè)生產(chǎn)中使用的主要手段，種子園是由優(yōu)樹無性系或家系按營建設(shè)計要求，實現(xiàn)集約經(jīng)營，以生產(chǎn)遺傳品質(zhì)和播種品質(zhì)優(yōu)良的林木良種為目的的特種人工林[1]。在種子園設(shè)計中，為了避免近交衰退的不利影響，需要盡量減少無性系之間的近交繁殖，以達到種子園遺傳效益最大化的目標。

現(xiàn)今，我國大部分種子園已經(jīng)結(jié)束了初級改良工作，正在向高世代育種發(fā)展，主要考慮的問題有種子園親本材料選擇、種子園設(shè)計、種子園植株管理、育種群體組建及長期育種計劃制定等[2]。其中，種子園設(shè)計的常用方法有：系統(tǒng)設(shè)計(System design，SD)(又稱為順序錯位排列)、完全隨機設(shè)計(Complete random design，CR)、隨機完全區(qū)組設(shè)計(Randomized complete block，RCB)等[3]。國內(nèi)學(xué)者在種子園設(shè)計領(lǐng)域的研究較少，代表性的工作有：許魯平[4]、申文輝等[5]和袁虎威等[6]將無性系材料分為不同區(qū)組，分別采用約束變換區(qū)組設(shè)計、約束分組結(jié)合隨機完全分組設(shè)計和不平衡、不完全固定區(qū)組設(shè)計的配置方式；程祥等[7]采用無性系順序錯位排列的固定配置法進行設(shè)計，以上研究工作均未充分考慮親本之間的不同親緣關(guān)系以及最大化近交距離。國外學(xué)者對種子園設(shè)計研究更為深入，BELL等[8]采用重置近鄰設(shè)計，LSTIBUREK等[9]提出了最小化近交(Minimum-inbreeding, MI)設(shè)計方案，EI-KASSABY等[10]提出了隨機、重復(fù)、交錯的無性系行(Randomized, replicated, staggered clonal-row, R2SCR)種子園設(shè)計方案，LSTIBUREK等[11]通過遺傳禁忌算法以及合并獨立的MI算法改進了原有的MI方案，CHALOUPKOVA等[12]提出最優(yōu)近鄰算法(ONA)，以上學(xué)者在種子園設(shè)計方法上更注重考慮不同的親緣關(guān)系和親緣比例，研究工作的共同點是：無性系的分株數(shù)量以及非平衡設(shè)計比例均為已知條件；在固定的無性系種類和無性系分株比例中完成種子園設(shè)計。

在現(xiàn)有研究工作中，種子園無性系的親緣關(guān)系、數(shù)量和比例都是研究者實驗假設(shè)的前提條件，如：LSTIBUREK等[9]在實驗中假設(shè)種子園的無性系只存在半同胞親緣關(guān)系或無親緣關(guān)系等不同親緣關(guān)系的場景，但是在現(xiàn)實種子園中親緣關(guān)系和遺傳多樣性往往會比假設(shè)條件更復(fù)雜、更混亂。因此，明確種子園中的親緣關(guān)系，可以有效控制近交，提高子代雜合率。

本文僅基于無性系種類和無性系之間的遺傳距離，利用改進型果蠅算法來使種子園中的無性系在滿足最小近親繁殖的同時得到較高的遺傳效益，以及種子園的較優(yōu)配置方案和合理的親本分株比例。最后將改進的果蠅算法與傳統(tǒng)的種子園配置方法(CR、RCB)、遺傳算法(GA)以及其他兩個改進果蠅算法進行比較，以評估該設(shè)計方案的優(yōu)越性，并為高世代種子園的設(shè)計提供參考。

1 種子園設(shè)計

1.1 問題描述

假設(shè)存在一個規(guī)模為M×N(M行、N列)的種子園，需要從已知遺傳距離的T株親本中挑選出合適的親本比例進行栽種。在栽種過程中，除種子園的四周位置以外，每個中心位置由8個位置包圍，這8個位置作為中心位置的近鄰位置，分為兩種情況：正對近鄰和斜角近鄰，易知種子園四角位置只存在3個近鄰位置，剩下的周邊位置都存在5個近鄰位置。其中近鄰位置作為中心位置計算近鄰距離的重點考慮對象，因為親本之間的遺傳距離越小，則說明親本之間的親緣關(guān)系越近，所以應(yīng)該選擇遺傳距離較大的植株作為近鄰，更利于提高種子園的遺傳效益。除此之外還需要考慮同一親本的不同分株之間的距離影響，即同一親本的不同分株應(yīng)該栽種在適當(dāng)遠的距離，以避免自交現(xiàn)象以及優(yōu)先交配親本組合不斷循環(huán)的現(xiàn)象。本文實驗中，不考慮種子園地形、氣候等客觀因素。

1.2 目標函數(shù)

該問題的目標函數(shù)為

(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 不同目標函數(shù)對比

將本文設(shè)計的目標函數(shù)分別與現(xiàn)有MI算法[9]中的目標函數(shù)

(2)

以及改進MI算法[11]中存在半同胞(half-siblings, h-s)關(guān)系時的目標函數(shù)

(3)

式中dil,jl——第l株無性系的第i株和第j株分株之間的距離

Nc——無性系的數(shù)量

NR——同一無性系的分株數(shù)量

dmin——所有種類的無性系中同一無性系的任何兩個分株的所有可能距離的平方倒數(shù)和

dmin(h-s)——所有具有半同胞關(guān)系的成對分株距離之和

進行對比。

本文提出的算法以及兩個對比算法的最終目標都是為了最小化近親繁殖，其中，式(2)只將同一無性系的不同分株之間的距離作為近親繁殖的衡量標準，而同一無性系的不同分株間授粉稱為自交，所以式(2)其實只考慮了自交繁殖現(xiàn)象，并沒有考慮到不同無性系之間的遺傳關(guān)系以及它們之間的交配現(xiàn)象。式(3)加入了半同胞情況，但是在現(xiàn)實的種子園中特別是高世代種子園中，無性系之間具有的親緣關(guān)系往往更復(fù)雜，算法設(shè)計不能單純只考慮同一無性系之間的自交和半同胞現(xiàn)象。

本文提出的目標函數(shù)不僅考慮了同一無性系的不同分株之間的影響，還引入了遺傳距離作為不同無性系之間的親緣關(guān)系的衡量標準，利用遺傳距離計算種子園中所有不同無性系之間的近交繁殖距離。

2 材料與方法

2.1 材料與數(shù)據(jù)來源

實驗材料是從內(nèi)蒙古紅花爾基樟子松(Pinussylvestrisvar.mongolica)國家良種基地中的初級種子園、1.5代種子園、2代種子園采集的當(dāng)年生針葉3～5針的樟子松無性系材料，由于SNP分子標記技術(shù)成本較高、分型技術(shù)不太成熟，因此本實驗材料是通過簡單、成熟、成本低的SSR分子標記法提取的樟子松基因組DNA，從15對SSR引物中選擇其中多態(tài)性良好、穩(wěn)定、清晰的11對SSR引物作為實驗所用引物，利用GeneMarker V2.2軟件對條帶信息進行比對后，再基于等位基因頻率的Nei 1983距離計算得到樟子松無性系材料中不同無性系之間的遺傳距離[13]，作為本文的實驗數(shù)據(jù)。

2.2 研究方法

種子園設(shè)計本質(zhì)上是一個二次分配問題(Quadratic assignment problem，QAP)，作為一個NP難問題，QAP的解決方法有窮舉法和智能算法等。受計算能力和時間制約，窮舉法只適用于解決小規(guī)模的QAP案例，智能算法是解決大規(guī)模QAP案例的可行、有效手段。

較為成熟的智能算法有遺傳算法(Genetic algorithms，GA)、禁忌搜索(Tabu search，TS)、粒子群算法(Particle swarm optimization，PSO)等，其中，LSTIBUREK等[9]的MI算法利用了禁忌搜索算法，且LSTIBUREK等[11]在改進的MI算法中采用了遺傳禁忌搜索算法，進一步優(yōu)化了原有的MI算法，使該方案適用于更大更復(fù)雜的高世代種子園。雖然經(jīng)典的智能優(yōu)化算法已經(jīng)得到了廣泛的應(yīng)用，但是各自都存在一些明顯的不足，這與待解決領(lǐng)域問題相關(guān)。近些年新的仿生智能算法也不斷出現(xiàn)，例如果蠅優(yōu)化算法(Fruit fly optimization algorithm，F(xiàn)OA)、蟻群算法(Ant colony optimization，ACO)、魚群算法(Fish swarm algorithm，F(xiàn)SA)、免疫算法(Immune algorithm，IA)等，相比其他優(yōu)化算法，F(xiàn)OA具有簡單、參數(shù)少、易調(diào)節(jié)、計算量小、尋優(yōu)精度較高等優(yōu)點，易于實現(xiàn)，更適合解決結(jié)構(gòu)復(fù)雜的現(xiàn)實問題[14-15]。

2.2.1果蠅優(yōu)化算法

果蠅優(yōu)化算法主要分為嗅覺覓食環(huán)節(jié)和視覺覓食環(huán)節(jié)，然后對這兩個覓食環(huán)節(jié)進行不斷的迭代最終實現(xiàn)該果蠅種群的進化，獲得該算法對問題的最優(yōu)解。但是，到目前為止，果蠅算法多用于解決連續(xù)型目標優(yōu)化問題，對于離散型目標函數(shù)研究還較少，用該算法解決的離散型目標函數(shù)最常見的有旅行商問題(Traveling salesman problem, TSP)[16-18]和調(diào)度問題[19-20]，在解決TSP問題中王克甫等[16]引入局部最優(yōu)概念并采用自適應(yīng)步長策略抑制了早熟現(xiàn)象并提高了搜索效率，段艷明等[18]將果蠅算法和遺傳算法以及C2Opt算子相結(jié)合加快了解決問題的局部搜索能力和收斂速度。基于不同學(xué)者對果蠅算法的改進方法，本文針對無性系種子園配置問題對傳統(tǒng)果蠅優(yōu)化算法的不足進行了改進，提出改進型果蠅算法(Improved fruit fly optimization algorithm，IFOA)。

2.2.2改進果蠅優(yōu)化算法

無性系種子園設(shè)計問題不同于連續(xù)型優(yōu)化問題，對T株無性系按隨機比例栽種到規(guī)模為M×N的種子園的配置方案為一只果蠅，配置方案中的任意一個位置為一個基因位，其中，改變果蠅的任意一個基因位，稱為一次覓食；在覓食過程中對基因位的選擇行為稱為覓食的方向選擇，每只果蠅的濃度判定函數(shù)為式(1)。

為了使FOA能夠更好地解決本文的無性系種子園配置問題，對FOA做了以下改進：

(1)輪盤賭法初始化種群。果蠅算法的初始值會在一定程度上影響果蠅算法的收斂速度，本實驗采用輪盤賭法初始化果蠅種群，將所選的T株親本的遺傳距離的倒數(shù)作為輪盤賭法中的概率p(i)，即

(4)

式中 sum(G)——候選親本中所有遺傳距離的和

G(i,j)——第i株無性系到第j株無性系的遺傳距離

(2)過濾無效樹種選擇。在每次覓食過程中，獲取種子園當(dāng)前基因位的正對近鄰的樹種，為節(jié)約覓食搜索時間，加快搜索速度，過濾與正對近鄰相同無性系的選擇，避免無效的覓食。

(3)步長選擇。FOA的步長會影響算法的效率和精度，當(dāng)步長較大時有利于全局搜索，加快搜索速度，但精度較低；步長較小時有利于局部搜索，提高算法精度，但搜索速度會下降。傳統(tǒng)的FOA是采用固定步長進行覓食，因此很難平衡算法的效率和精度。為兼顧兩者平衡，本文引入遺傳算法中的基因變異對步長進行調(diào)整，在迭代前期可設(shè)置較大的變異因子，有利于加快算法的搜索和收斂速度；在迭代后期減小變異因子，在保證算法精度的同時也有助于算法跳出局部最優(yōu)。

(4)覓食方向。在迭代前期為了加快搜索速度，將覓食方向定為種子園中濃度最差的基因位，在迭代后期為了避免算法陷入局部最優(yōu)，將覓食方向改為隨機方向。

在求解無性系種子園配置問題時需要把FOA的連續(xù)解空間對應(yīng)到種子園配置方案，其主要步驟見圖1。

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

3 實驗與結(jié)果分析

3.1 實驗準備

3.1.1數(shù)據(jù)材料

從內(nèi)蒙古紅花爾基樟子松國家良種基地中的1代種子園、1.5代種子園、2代種子園遺傳距離的Excel文件中，隨機選取其中13個無性系和它們的遺傳距離作為實驗數(shù)據(jù)。

3.1.2參數(shù)設(shè)置

本實驗中，種子園的規(guī)模假設(shè)為9×9(即9行、9列)，初始種群大小為25，迭代次數(shù)為800，迭代前期變異因子為0.4，迭代后期變異因子為0.02，在此基礎(chǔ)上將本實驗與傳統(tǒng)的種子園配置方法(CR、RCB)、遺傳算法(GA)以及其他兩個改進果蠅算法：MFFA[16]、ASFOA[18]進行了比較，以上所有算法均

采用式(1)作為目標函數(shù)。

3.1.3實驗說明

該算法并不只針對數(shù)據(jù)材料中選取的13個無性系或參數(shù)設(shè)置中的9×9固定大小的種子園，它適用于任意無性系的種類和數(shù)量以及任意規(guī)模的種子園。

3.2 算法比較

3.2.1綜合對比

通常全局優(yōu)化算法結(jié)果并不是問題的最小值，而是最小值的近似解，即最終的結(jié)果存在一定的誤差。為了避免偶然事件出現(xiàn)影響對比結(jié)果，將各算法在以上參數(shù)設(shè)置條件下分別執(zhí)行200次后得到對比結(jié)果，見表1。

從表1可以發(fā)現(xiàn)，本文提出的IFOA的平均值、最大值、最小值分別為1 169.780、1 198.153、1 137.774，都明顯低于其他算法，其中IFOA的最大值也明顯低于其他算法的最小值，其次結(jié)果較好的為GA。相對于傳統(tǒng)算法而言，智能優(yōu)化算法具有明顯優(yōu)勢。

表1 不同算法分別執(zhí)行200次后適應(yīng)度對比Tab.1 Comparison of fitness of each algorithm after 200 times respectively

3.2.2收斂速度比較

從200次試驗中隨機選擇6次運行結(jié)果，作為各算法的收斂速度，見圖2，IFOA的收斂速度以及效果明顯優(yōu)于其他算法。

3.2.3不同方法設(shè)計對比

將本文的IFOA算法得出的設(shè)計方案，分別與效果較好的智能算法GA設(shè)計方案以及傳統(tǒng)算法CR設(shè)計方案進行對比，如圖3所示。

圖3a存在多處同一無性系的分株作為直接鄰居，如第1行第3個位置和第4個位置都是無性系9的分株；圖3b、3c沒有同一無性系的分株作為直接鄰居出現(xiàn)；圖3c無性系3和無性系7作為直接鄰居出現(xiàn)的概率明顯高于圖3b，原因是無性系3和無性系7之間的遺傳距離是實驗數(shù)據(jù)中遺傳距離最大的組合，因此會優(yōu)先交配，有利于雜交育種；圖3c沒有無性系2，因為無性系2與其他無性系之間的遺傳距離都較小，即無性系2與其他無性系都具有較近的遺傳距離，因此不建議栽種。

圖2 各算法收斂情況對比Fig.2 Comparison of convergence of different algorithms

圖3 13個無性系在9×9種子園的設(shè)計方案Fig.3 Design of 13 clones in 9×9 orchard

4 結(jié)束語

將無性系的分株數(shù)量作為未知條件，基于無性系之間的遺傳距離，不僅考慮了同一親本的分株影響，也考慮了近鄰之間的影響，使得種子園最終能夠達到較優(yōu)的遺傳效益，通過改進型果蠅算法實現(xiàn)該目標，并得到較優(yōu)的種子園設(shè)計方案和親本分株比例。計算結(jié)果表明IFOA的收斂速度和效果優(yōu)于其他算法，并且IFOA的設(shè)計方案中同一無性系分株不會作為近鄰出現(xiàn)，當(dāng)存在親本組合的遺傳距離明顯大于其他組合時，即親緣關(guān)系最遠組合，該組合通常會作為優(yōu)先交配組合。