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

齊建東 買晶晶 劉春霞 李 偉

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

0 引言

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

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

種子園中無性系的花期和花粉量會影響種子園子代際的遺傳結(jié)構(gòu)。其中，開花同步數(shù)據(jù)會影響無性系之間的雜交繁殖，如果無性系之間花期不遇，可能導(dǎo)致大多數(shù)無性系之間的授粉期和散粉期錯開，從而嚴(yán)重影響種子的產(chǎn)量和品種[15]?；ǚ哿恐挥羞_到有效值才能授粉成功，當(dāng)花粉量傳播不足時，會導(dǎo)致授粉失敗，出種率降低，敗育率升高[16]。可見，花期和花粉量是設(shè)計各世代種子園內(nèi)無性系配置的關(guān)鍵因素，只有種子園內(nèi)的花期基本一致或者相鄰無性系花期一致，并且花粉量足夠，才能授粉成功，進行正常繁殖。目前，尚未見同時引入花期和花粉量作為種子園設(shè)計約束條件的相關(guān)報道。

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

1 種子園設(shè)計問題

1.1 問題描述

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

無性系的具體花期時間與花粉量，會受當(dāng)年的氣候影響(如風(fēng)、溫度和濕度等)，也會受當(dāng)?shù)氐牡匦魏突ǚ圩陨硖卣鞯扔绊慬19]，本文研究工作暫不考慮氣候、地形和不同花粉自身特性等因素。

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

圖1 種子園栽種示意圖

1.2 目標(biāo)函數(shù)

本文設(shè)計目標(biāo)函數(shù)為

(1)

其中

R=nh

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

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

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

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

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

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

w1、w2——限制因子

R——花粉傳播半徑

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

h——樹高s——行距

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

Q——花粉傳播半徑為R的圓范圍內(nèi)的無性系數(shù)量

u——無性系數(shù)量

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

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

(2)

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

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

2 材料與方法

2.1 數(shù)據(jù)來源

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

2.2 研究方法

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

2.2.1標(biāo)準(zhǔn)果蠅算法

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

2.2.2TIASFOA算法

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

(1)自適應(yīng)步長

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

(3)

其中

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

g——當(dāng)前迭代次序

maxgen——最大迭代次數(shù)

Xorchard——種子園X軸方向可栽種位置數(shù)量

Yorchard——種子園Y軸方向可栽種位置數(shù)量

(2)多種群策略

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

(3)尋優(yōu)過程

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

2.2.3TIASFOA算法流程

TIASFOA算法流程如下：

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

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

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

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

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

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

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

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

圖2 雙種群改進型自適應(yīng)步長的果蠅優(yōu)化算法流程圖

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

3.1 實驗設(shè)置

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

3.2 結(jié)果分析

3.2.1不同種群規(guī)模下的適應(yīng)度對比

對TIASFOA、IFOA、PSO、GA算法在種群20～60范圍內(nèi)以式(1)為目標(biāo)函數(shù)計算近似最優(yōu)適應(yīng)度dmin，并分別執(zhí)行200次循環(huán)后對平均值、最大值、最小值和方差進行比較分析。

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

表1 4種算法在不同種群規(guī)模下的適應(yīng)度對比

3.2.2算法收斂比較

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

3.2.3不同種群規(guī)模下的花期對比

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

圖3 各算法收斂情況對比

表2 不同種群規(guī)模下4種算法分別執(zhí)行200次的無性系數(shù)量對比

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

4 結(jié)論

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

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

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