無性系種子園最小近交配置研究

齊建東 買晶晶 劉春霞 李 偉

(1.北京林業大學信息學院，北京 100083；2.中國科學院軟件研究所，北京 100190；3.北京林業大學生物科學與技術學院，北京 100083)

0 引言

無性系種子園是以優良無性系個體為材料、用無性繁殖的方式建立的種子園。無性系種子園中的親本來源清楚，有利于保持樹源的優良品質，便于集約經營管理，但容易產生自交現象。因此，需要通過對種子園的無性系進行合理配置，以避免自交和近交現象。

20世紀60年代，我國試建了第1批初級無性系種子園。20世紀70年代末和80年代初，造林工作得到迅速發展。目前，很多初代林木樹種的改良工作已經完成，正進入高世代改良階段。在種子園設計方面，國內學者多采用傳統的無性系配置方法進行配置設計，如賈乃光[1]、程祥等[2]選用順序錯位排列設計法，梁一池[3]利用隨機完全區組設計方法，許魯平[4]采用約束變換區組設計，申文輝等[5]使用約束的隨機完全區組設計，鄭仁華等[6]、謝汝根[7]、蘇順德等[8]均采用完全隨機排列法進行種子園設計。現實中的種子園是一個龐大而復雜的交配系統，在進行設計工作時，以上學者沒有考慮親本的親緣關系，僅有袁虎威等[9]從分子水平上獲取優良單株之間的親緣關系，并采用不平衡、不完全固定區組方法將實驗分析得到的親緣關系結果應用于第2代種子園的無性系配置設計；王晴等[10]引用無性系親本間的遺傳距離，設計了一種改進型自適應并行遺傳算法，實現了田間設計。LSTIBREK等[11]在2010年提出了最小化近交(Minimum inbreeding, MI)設計方案，在實驗中設計了5種不同場景實驗，考慮了更復雜的親緣結構(如無親緣、半同胞和親子關系等)。LSTIBREK等[12]在2015年進一步提出將擴展全局(遺傳禁忌)算法(Extended global(genetic-tabu)algorithm, EGA)用于解決經營林業的現實問題，其在MI的基礎上添加了半同胞約束條件。CHALOUPKOV等[13]將提出的最優近鄰算法(Optimum neighborhood algorithm, ONA)應用于平衡和非平衡無性系規模的實驗中。近幾年來，在種子園設計方面的研究更側重于無性系的親緣關系，而其他的影響因素(如花期和花粉量等)考慮較少，僅王強金[14]提及的光澤華橋國有林場配置方案中考慮了物候期，但物候期的劃分不夠明確。LSTIBREK等[12]在2015年提出，可以使用假設情景的方式引入開花同步、育種值、花粉產量等因素，但并沒有提供相關實驗。

種子園中無性系的花期和花粉量會影響種子園子代際的遺傳結構。其中，開花同步數據會影響無性系之間的雜交繁殖，如果無性系之間花期不遇，可能導致大多數無性系之間的授粉期和散粉期錯開，從而嚴重影響種子的產量和品種[15]。花粉量只有達到有效值才能授粉成功，當花粉量傳播不足時，會導致授粉失敗，出種率降低，敗育率升高[16]。可見，花期和花粉量是設計各世代種子園內無性系配置的關鍵因素，只有種子園內的花期基本一致或者相鄰無性系花期一致，并且花粉量足夠，才能授粉成功，進行正常繁殖。目前，尚未見同時引入花期和花粉量作為種子園設計約束條件的相關報道。

齊建東等[17]通過無性系之間的遺傳距離作為親緣關系的衡量標準，利用改進型果蠅算法實現了無性系種子園的遺傳設計，本文在此基礎上，加入花期和花粉量作為新的約束條件，并設計雙種群改進型自適應步長的果蠅優化算法(Two-population improved adaptive step-length fruit fly optimization algorithm, TIASFOA)對無性系進行配置設計。

1 種子園設計問題

1.1 問題描述

依據T株親本之間的遺傳距離、花期和花粉傳播量，合理選擇親本及其分株進行無性系遺傳設計，在栽種過程中重點考慮近鄰位置的近交繁殖現象、同一無性系不同分株的自交情況，以及花粉有效傳播范圍內的無性系花期同步性和花粉接收量。無性系之間的遺傳距離越大，親緣關系越遠，更有利于雜交繁殖，配置時應該選擇遺傳距離較大的無性系作為近鄰；此外，在配置過程中，需要優先考慮花期問題，只有當無性系的花期同步時，散粉和授粉才有意義，應該將同一花期的無性系栽種距離盡量縮小[15]；其次考慮花粉量，根據經驗易知花粉量傳播與樹的高度相關，花粉傳播會隨著距離增加而減少，當花粉量低于有效花粉傳播量時，則無法成功授粉，因此，在有效傳播距離內盡量接受更多的花粉量[18]。

無性系的具體花期時間與花粉量，會受當年的氣候影響(如風、溫度和濕度等)，也會受當地的地形和花粉自身特征等影響[19]，本文研究工作暫不考慮氣候、地形和不同花粉自身特性等因素。

假設種子園是一個規模為M行×N列的規則的種子園，所有的樹高均為h，行間距為s，花粉是在靜風條件下進行傳播，在有效范圍的邊緣處只能接受到10%的花粉量(即有效花粉量不能低于10%)，花粉量在靜風條件下呈線性遞減。種子園栽種示意圖如圖1所示，每一個方格代表一個可以栽種的位置，如1-1表示第1行第1列的位置。

圖1 種子園栽種示意圖

1.2 目標函數

本文設計目標函數為

(1)

其中

R=nh

式中dmin——所有植株的近鄰距離和同一無性系所有分株距離之和與所有植株有效花粉量倒數和相加的最小值

Gij——第i株無性系和第j株無性系之間的遺傳距離，且第i株無性系和第j株無性系為正對近鄰

Gik——第i株無性系和第k株無性系之間的遺傳距離，且第i株無性系和第k株無性系為斜角近鄰

dit——第i株無性系和第t株無性系之間的物理距離，第i株無性系和第t株無性系為同一無性系親本的分株

Piq——第i株無性系與第q株無性系之間傳播的花粉量

q——以第i株樹為中心半徑為R范圍內的所有無性系的數量

w1、w2——限制因子

R——花粉傳播半徑

Pq——第q株樹向第i株樹傳播的花粉量

h——樹高s——行距

n——半徑與樹高的比值，本文設為3

Q——花粉傳播半徑為R的圓范圍內的無性系數量

u——無性系數量

diq——第i株無性系與第q株無性系之間的物理距離

同時考慮花期的影響，添加以下約束條件

(2)

式中Pj——第j株樹向第i株樹傳播的花粉量

Pj-self——第j株樹的花粉量

2 材料與方法

2.1 數據來源

從內蒙古紅花爾基樟子松國家良種基地1代種子園、1.5代種子園、2代種子園采集的當年生針葉3～5針的樟子松無性系材料，由于單核苷酸多態性(Single nucleotide polymorphism，SNP)分子標記技術成本較高、分型技術不太成熟，因此本文通過簡單、成熟、成本低的SSR分子標記法提取樟子松基因組DNA，從15對SSR引物中選擇多態性良好、穩定、清晰的11對SSR引物作為實驗引物，利用Gene Marker V2.2軟件對條帶信息進行比對后，基于等位基因頻率的Nei 1983距離計算得到樟子松無性系材料中不同無性系之間的遺傳距離[12]，作為本文的實驗數據。

2.2 研究方法

王晴等[10]、LSTIBREK等[12]均采用改進的遺傳算法對種子園內的無性系進行配置設計，齊建東等[17]利用改進型果蠅算法實現了種子園設計，結果均表明改進智能優化算法優于傳統設計方法。與較成熟的智能算法(如遺傳算法(Genetic algorithm，GA)、粒子群算法(Particle swarm optimization，PSO)等)相比，果蠅優化算法(Fruit fly optimization algorithm，FOA)作為新型仿生智能算法，具有簡單、易于實現等優點，獲得了國內外眾多學者的廣泛關注和研究，在醫學、生物、工程和科學等領域得到了應用，并且相對于計算量較大的啟發式智能優化算法，如人工蜂群算法(Artificial bee colony algorithm，ABC)、GA算法的多種群策略，FOA算法計算簡單，這使得FOA算法所消耗的資源更少。

2.2.1標準果蠅算法

FOA算法利用果蠅(Drosophilamelanogaster)個體嗅覺優勢獲取食物的味道，并將自身獲取的味道與其他果蠅個體共享；果蠅個體通過視覺比較得出種群中獲得最優食物味道的果蠅個體，然后向具有最優食物味道的果蠅位置聚集，并按照該搜索方式繼續搜索食物，直到找到食物為止[20]。FOA算法具有眾多優點，但無法直接應用于離散問題，且算法穩定性較差[21]。針對種子園遺傳設計問題，齊建東等[17]設計的IFOA算法可用于解決離散問題，但穩定性并沒有得到提高。本文對FOA算法和IFOA算法的不足進行改進，設計了雙種群改進型自適應步長果蠅優化算法(Two-population improved adaptive step-length fruit fly optimization algorithm, TIASFOA)。

2.2.2TIASFOA算法

在IFOA算法[17]的基礎上，本文進行了以下改進。

(1)自適應步長

FOA算法的尋優過程中，步長是一個常量函數，從而導致迭代前期收斂速度慢，迭代后期其尋優精度低；IFOA算法步長變異因子的調節需要人為控制；參考文獻[22-24]，根據種子園特點，設計了一個自適應逐步遞減的步長函數，計算公式為

(3)

其中

式中L——步長L0——初始步長

g——當前迭代次序

maxgen——最大迭代次數

Xorchard——種子園X軸方向可栽種位置數量

Yorchard——種子園Y軸方向可栽種位置數量

(2)多種群策略

在FOA算法和IFOA算法中均采用單種群尋優，種群多樣性降低，算法易陷入局部最優。TIASFOA算法將整個種群劃分為兩個規模相同的子種群，分別對兩個子種群進行獨立尋優操作，利用種群之間的信息交流機制，增加精英個體，保留最優解和次優解，子種群的其他果蠅分別以最優解和次優解為標準值，形成兩個新的子種群，按照該方法迭代尋優，直到滿足最大迭代次數，合并兩個子種群，輸出最優解以及最優解位置上的其他信息。

(3)尋優過程

為了在保留算法較優的收斂速度的同時，加強算法的隨機效果，提高果蠅個體跳出局部最優點的能力，在每次覓食時同時采用最優和隨機兩種覓食行為。增添一個覓食概率，通過多次實驗將其設置為0.8，使得果蠅個體在每次覓食(即尋優)時按照一定的概率采用最優覓食或隨機覓食。即在每次覓食過程中隨機生成一個隨機概率，當隨機概率小于覓食概率時，尋找果蠅個體的最差濃度基因位，進行最優覓食；當隨機概率大于覓食概率時，采用隨機機制。

2.2.3TIASFOA算法流程

TIASFOA算法流程如下：

(1)輪盤賭法初始化種群。設置實驗相關參數：種子園規模M行×N列、種群規模sizepop、最大迭代次數maxgen、覓食概率P。

(2)確定初始濃度和位置。初始化時保留最優解和次優解果蠅的濃度和位置，記為全局最優解和次優解個體濃度和位置。

(3)將種群劃分成2個相同規模的子種群。

(4)動態調整步長，分別對2個子種群進行獨立尋優操作。每次覓食開始前，先隨機生成概率Pi，然后判斷Pi和P的大小，當Pi

(5)對比2個子種群的解，保存所有解中的最優解和次優解個體以及他們對應的位置。

(6)對當前迭代中的最優解、次優解、全局最優解、全局次優解進行比較，判斷是否需要更新全局最優解和次優解個體濃度以及他們對應的位置。

(7)進入迭代尋優，如果迭代次數g小于最大迭代次數，則循環執行步驟(3)～(6)；否則合并2個子種群輸出最優種子園方案、該方案的最佳濃度以及其他相關信息。

TIASFOA算法的基本流程圖見圖2。

圖2 雙種群改進型自適應步長的果蠅優化算法流程圖

3 實驗與結果分析

3.1 實驗設置

種子園規模設定為9行×9列；TIASFOA算法最大迭代次數為1 000，行間距為1 m×1 m，樹高為1 m(本文實驗數據，可以自定義)；將花期設定在5月15—20日，不同單株之間的相互接受花粉量的范圍為50～500單位量之間，在此范圍內，隨機生成每株親本的花期和花粉量；花粉最大的傳播半徑為3倍樹高(即3 m)，在此基礎上與IFOA、GA、PSO算法進行對比，其中IFOA算法的參數設置：迭代前期變異因子為0.4，迭代后期變異因子為0.02；GA算法的參數設置：交配概率為0.8，變異概率為0.2；PSO算法的參數設置：學習因子為2，速度初始化為1。

3.2 結果分析

3.2.1不同種群規模下的適應度對比

對TIASFOA、IFOA、PSO、GA算法在種群20～60范圍內以式(1)為目標函數計算近似最優適應度dmin，并分別執行200次循環后對平均值、最大值、最小值和方差進行比較分析。

觀察表1，在不同種群規模下，PSO、GA算法的最小值、最大值、平均值、方差均大于TIASFOA、IFOA算法的對應值，且PSO、GA算法的最小值始終大于TIASFOA的最大值,說明PSO、GA算法表現較差；隨著種群規模的增加，IFOA算法的最大值減小，在種群規模為50的情況下IFOA算法取得最小值137.097，大于TIASFOA算法的最小值132.733，說明IFOA算法陷入了局部最優；在不同種群規模下，IFOA算法最小值、最大值、平均值、方差均大于TIASFOA算法，說明TIASFOA算法優于IFOA算法。在種群規模為60的情況下，TIASFOA算法取得了方差最小值1.288，小于其他算法。在種群規模為50的情況下，TIASFOA算法取得最小值132.733。TIASFOA算法從適應度和穩定性方面均優于其他3種算法。隨著種群規模的增大，各算法的方差會有所減小，但時間代價更大，TIASFOA算法可以在種群規模20下得到較優值，時間消耗少。

表1 4種算法在不同種群規模下的適應度對比

3.2.2算法收斂比較

為在較短的時間內獲取較優的適應度，本文從種群規模為20的200次實驗中隨機選用6次結果進行收斂速度和效果對比(圖3)。從圖3中可以看到，TIASFOA算法的收斂效果均優于其他算法。IFOA算法在迭代250次左右適應度達到最小值，而TISFOA算法在迭代700次左右達到最小值，說明IFOA算法收斂速度最快，但最優適應度仍大于TIASFOA算法，說明IFOA算法陷入局部最優，TIASFOA算法可以跳出局部最優，得到更優解。

3.2.3不同種群規模下的花期對比

表2給出了在種群規模20～60下，TIASFOA、IFOA、PSO、GA算法分別執行200次，得出的無性系種子園方案在距離中心位置的半徑為1、2、3 m的傳播范圍內，種子園花期相同或相鄰無性系數量的平均值、最大值和最小值。

圖3 各算法收斂情況對比

表2 不同種群規模下4種算法分別執行200次的無性系數量對比

觀察表2，在不同種群規模下距離中心位置半徑為1、2、3 m的傳播范圍內，種子園內花期相同或相鄰的無性系數量，PSO算法和GA算法始終劣于TIASFOA算法和IFOA算法；TIASFOA算法計算出來的無性系數量的平均值均優于IFOA算法，IFOA算法計算出來的無性系數量的最大值和最小值等于或略大于TIASFOA算法，例如表2總計中，TIASFOA算法和IFOA算法在種群規模為20時，種子園內花期相鄰的情況下，無性系數量的最小值都為123；在種群規模為40時，種子園內花期相鄰情況下，無性系數量的最小值IFOA算法結果為134，大于TIASFOA算法的結果130，這是受算法的隨機性影響產生的偶然值。整體來看，TIASFOA算法在花期相同或相鄰時平均值和最大值都更優，可以更好地使整個種子園維持較好的花期一致性。

4 結論

(1)目標函數不僅考慮了無性系間的遺傳距離，而且引入花期和花粉量作為約束條件，并考慮了花粉量在傳播過程中不斷減小的可能性，不局限于實驗變量，靈活性較好。

(2)設計的TIASFOA算法擴大了搜索空間，引入了多種群的信息交流機制，可以進一步跳出局部最優，得到較優的適應度和較好的花期一致性。

(3)在實際應用中，可以根據現實種子園的樹高規定花粉傳播距離、設置真實行間距等，并使用種群規模為20的TIASFOA算法對無性系種子園進行設計。