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

齊建東 劉春霞 崔曉暉 李 偉

(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的設計方案中同一無性系分株不會作為近鄰出現，當存在親本組合的遺傳距離明顯大于其他組合時，即親緣關系最遠組合，該組合通常會作為優先交配組合。