齊建東 劉春霞 崔曉暉 李 偉
(1.北京林業(yè)大學(xué)信息學(xué)院, 北京 100083; 2.北京林業(yè)大學(xué)生物科學(xué)與技術(shù)學(xué)院, 北京 100083)
當(dāng)前我國林業(yè)發(fā)展中,營林造林主要依靠林木育種。其中種子園是林業(yè)生產(chǎn)中使用的主要手段,種子園是由優(yōu)樹無性系或家系按營建設(shè)計(jì)要求,實(shí)現(xiàn)集約經(jīng)營,以生產(chǎn)遺傳品質(zhì)和播種品質(zhì)優(yōu)良的林木良種為目的的特種人工林[1]。在種子園設(shè)計(jì)中,為了避免近交衰退的不利影響,需要盡量減少無性系之間的近交繁殖,以達(dá)到種子園遺傳效益最大化的目標(biāo)。
現(xiàn)今,我國大部分種子園已經(jīng)結(jié)束了初級改良工作,正在向高世代育種發(fā)展,主要考慮的問題有種子園親本材料選擇、種子園設(shè)計(jì)、種子園植株管理、育種群體組建及長期育種計(jì)劃制定等[2]。其中,種子園設(shè)計(jì)的常用方法有:系統(tǒng)設(shè)計(jì)(System design,SD)(又稱為順序錯(cuò)位排列)、完全隨機(jī)設(shè)計(jì)(Complete random design,CR)、隨機(jī)完全區(qū)組設(shè)計(jì)(Randomized complete block,RCB)等[3]。國內(nèi)學(xué)者在種子園設(shè)計(jì)領(lǐng)域的研究較少,代表性的工作有:許魯平[4]、申文輝等[5]和袁虎威等[6]將無性系材料分為不同區(qū)組,分別采用約束變換區(qū)組設(shè)計(jì)、約束分組結(jié)合隨機(jī)完全分組設(shè)計(jì)和不平衡、不完全固定區(qū)組設(shè)計(jì)的配置方式;程祥等[7]采用無性系順序錯(cuò)位排列的固定配置法進(jìn)行設(shè)計(jì),以上研究工作均未充分考慮親本之間的不同親緣關(guān)系以及最大化近交距離。國外學(xué)者對種子園設(shè)計(jì)研究更為深入,BELL等[8]采用重置近鄰設(shè)計(jì),LSTIBUREK等[9]提出了最小化近交(Minimum-inbreeding, MI)設(shè)計(jì)方案,EI-KASSABY等[10]提出了隨機(jī)、重復(fù)、交錯(cuò)的無性系行(Randomized, replicated, staggered clonal-row, R2SCR)種子園設(shè)計(jì)方案,LSTIBUREK等[11]通過遺傳禁忌算法以及合并獨(dú)立的MI算法改進(jìn)了原有的MI方案,CHALOUPKOVA等[12]提出最優(yōu)近鄰算法(ONA),以上學(xué)者在種子園設(shè)計(jì)方法上更注重考慮不同的親緣關(guān)系和親緣比例,研究工作的共同點(diǎn)是:無性系的分株數(shù)量以及非平衡設(shè)計(jì)比例均為已知條件;在固定的無性系種類和無性系分株比例中完成種子園設(shè)計(jì)。
在現(xiàn)有研究工作中,種子園無性系的親緣關(guān)系、數(shù)量和比例都是研究者實(shí)驗(yàn)假設(shè)的前提條件,如:LSTIBUREK等[9]在實(shí)驗(yàn)中假設(shè)種子園的無性系只存在半同胞親緣關(guān)系或無親緣關(guān)系等不同親緣關(guān)系的場景,但是在現(xiàn)實(shí)種子園中親緣關(guān)系和遺傳多樣性往往會(huì)比假設(shè)條件更復(fù)雜、更混亂。因此,明確種子園中的親緣關(guān)系,可以有效控制近交,提高子代雜合率。
本文僅基于無性系種類和無性系之間的遺傳距離,利用改進(jìn)型果蠅算法來使種子園中的無性系在滿足最小近親繁殖的同時(shí)得到較高的遺傳效益,以及種子園的較優(yōu)配置方案和合理的親本分株比例。最后將改進(jìn)的果蠅算法與傳統(tǒng)的種子園配置方法(CR、RCB)、遺傳算法(GA)以及其他兩個(gè)改進(jìn)果蠅算法進(jìn)行比較,以評估該設(shè)計(jì)方案的優(yōu)越性,并為高世代種子園的設(shè)計(jì)提供參考。
假設(shè)存在一個(gè)規(guī)模為M×N(M行、N列)的種子園,需要從已知遺傳距離的T株親本中挑選出合適的親本比例進(jìn)行栽種。在栽種過程中,除種子園的四周位置以外,每個(gè)中心位置由8個(gè)位置包圍,這8個(gè)位置作為中心位置的近鄰位置,分為兩種情況:正對近鄰和斜角近鄰,易知種子園四角位置只存在3個(gè)近鄰位置,剩下的周邊位置都存在5個(gè)近鄰位置。其中近鄰位置作為中心位置計(jì)算近鄰距離的重點(diǎn)考慮對象,因?yàn)橛H本之間的遺傳距離越小,則說明親本之間的親緣關(guān)系越近,所以應(yīng)該選擇遺傳距離較大的植株作為近鄰,更利于提高種子園的遺傳效益。除此之外還需要考慮同一親本的不同分株之間的距離影響,即同一親本的不同分株應(yīng)該栽種在適當(dāng)遠(yuǎn)的距離,以避免自交現(xiàn)象以及優(yōu)先交配親本組合不斷循環(huán)的現(xiàn)象。本文實(shí)驗(yàn)中,不考慮種子園地形、氣候等客觀因素。
該問題的目標(biāo)函數(shù)為
(i≠j,k,t且j≠k)
(1)
式中fmin——所有植株的近鄰距離與同一無性系所有分株距離之和的最小值
Gij——第i株無性系和第j株無性系之間的遺傳距離, 第i株無性系和第j株無性系為正對近鄰
Gik——第i株無性系和第k株無性系之間的遺傳距離,第i株無性系和第k株無性系為斜角近鄰
dit——第i株無性系和第t株無性系之間的物理距離,第i株無性系和第t株無性系為同一無性系親本的分株
將本文設(shè)計(jì)的目標(biāo)函數(shù)分別與現(xiàn)有MI算法[9]中的目標(biāo)函數(shù)
(2)
以及改進(jìn)MI算法[11]中存在半同胞(half-siblings, h-s)關(guān)系時(shí)的目標(biāo)函數(shù)
(3)
式中dil,jl——第l株無性系的第i株和第j株分株之間的距離
Nc——無性系的數(shù)量
NR——同一無性系的分株數(shù)量
dmin——所有種類的無性系中同一無性系的任何兩個(gè)分株的所有可能距離的平方倒數(shù)和
dmin(h-s)——所有具有半同胞關(guān)系的成對分株距離之和
進(jìn)行對比。
本文提出的算法以及兩個(gè)對比算法的最終目標(biāo)都是為了最小化近親繁殖,其中,式(2)只將同一無性系的不同分株之間的距離作為近親繁殖的衡量標(biāo)準(zhǔn),而同一無性系的不同分株間授粉稱為自交,所以式(2)其實(shí)只考慮了自交繁殖現(xiàn)象,并沒有考慮到不同無性系之間的遺傳關(guān)系以及它們之間的交配現(xiàn)象。式(3)加入了半同胞情況,但是在現(xiàn)實(shí)的種子園中特別是高世代種子園中,無性系之間具有的親緣關(guān)系往往更復(fù)雜,算法設(shè)計(jì)不能單純只考慮同一無性系之間的自交和半同胞現(xiàn)象。
本文提出的目標(biāo)函數(shù)不僅考慮了同一無性系的不同分株之間的影響,還引入了遺傳距離作為不同無性系之間的親緣關(guān)系的衡量標(biāo)準(zhǔn),利用遺傳距離計(jì)算種子園中所有不同無性系之間的近交繁殖距離。
實(shí)驗(yàn)材料是從內(nèi)蒙古紅花爾基樟子松(Pinussylvestrisvar.mongolica)國家良種基地中的初級種子園、1.5代種子園、2代種子園采集的當(dāng)年生針葉3~5針的樟子松無性系材料,由于SNP分子標(biāo)記技術(shù)成本較高、分型技術(shù)不太成熟,因此本實(shí)驗(yàn)材料是通過簡單、成熟、成本低的SSR分子標(biāo)記法提取的樟子松基因組DNA,從15對SSR引物中選擇其中多態(tài)性良好、穩(wěn)定、清晰的11對SSR引物作為實(shí)驗(yàn)所用引物,利用GeneMarker V2.2軟件對條帶信息進(jìn)行比對后,再基于等位基因頻率的Nei 1983距離計(jì)算得到樟子松無性系材料中不同無性系之間的遺傳距離[13],作為本文的實(shí)驗(yàn)數(shù)據(jù)。
種子園設(shè)計(jì)本質(zhì)上是一個(gè)二次分配問題(Quadratic assignment problem,QAP),作為一個(gè)NP難問題,QAP的解決方法有窮舉法和智能算法等。受計(jì)算能力和時(shí)間制約,窮舉法只適用于解決小規(guī)模的QAP案例,智能算法是解決大規(guī)模QAP案例的可行、有效手段。
較為成熟的智能算法有遺傳算法(Genetic algorithms,GA)、禁忌搜索(Tabu search,TS)、粒子群算法(Particle swarm optimization,PSO)等,其中,LSTIBUREK等[9]的MI算法利用了禁忌搜索算法,且LSTIBUREK等[11]在改進(jìn)的MI算法中采用了遺傳禁忌搜索算法,進(jìn)一步優(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é)、計(jì)算量小、尋優(yōu)精度較高等優(yōu)點(diǎn),易于實(shí)現(xiàn),更適合解決結(jié)構(gòu)復(fù)雜的現(xiàn)實(shí)問題[14-15]。
2.2.1果蠅優(yōu)化算法
果蠅優(yōu)化算法主要分為嗅覺覓食環(huán)節(jié)和視覺覓食環(huán)節(jié),然后對這兩個(gè)覓食環(huán)節(jié)進(jìn)行不斷的迭代最終實(shí)現(xiàn)該果蠅種群的進(jìn)化,獲得該算法對問題的最優(yōu)解。但是,到目前為止,果蠅算法多用于解決連續(xù)型目標(biāo)優(yōu)化問題,對于離散型目標(biāo)函數(shù)研究還較少,用該算法解決的離散型目標(biāo)函數(shù)最常見的有旅行商問題(Traveling salesman problem, TSP)[16-18]和調(diào)度問題[19-20],在解決TSP問題中王克甫等[16]引入局部最優(yōu)概念并采用自適應(yīng)步長策略抑制了早熟現(xiàn)象并提高了搜索效率,段艷明等[18]將果蠅算法和遺傳算法以及C2Opt算子相結(jié)合加快了解決問題的局部搜索能力和收斂速度?;诓煌瑢W(xué)者對果蠅算法的改進(jìn)方法,本文針對無性系種子園配置問題對傳統(tǒng)果蠅優(yōu)化算法的不足進(jìn)行了改進(jìn),提出改進(jìn)型果蠅算法(Improved fruit fly optimization algorithm,IFOA)。
2.2.2改進(jìn)果蠅優(yōu)化算法
無性系種子園設(shè)計(jì)問題不同于連續(xù)型優(yōu)化問題,對T株無性系按隨機(jī)比例栽種到規(guī)模為M×N的種子園的配置方案為一只果蠅,配置方案中的任意一個(gè)位置為一個(gè)基因位,其中,改變果蠅的任意一個(gè)基因位,稱為一次覓食;在覓食過程中對基因位的選擇行為稱為覓食的方向選擇,每只果蠅的濃度判定函數(shù)為式(1)。
為了使FOA能夠更好地解決本文的無性系種子園配置問題,對FOA做了以下改進(jìn):
(1)輪盤賭法初始化種群。果蠅算法的初始值會(huì)在一定程度上影響果蠅算法的收斂速度,本實(shí)驗(yàn)采用輪盤賭法初始化果蠅種群,將所選的T株親本的遺傳距離的倒數(shù)作為輪盤賭法中的概率p(i),即
(4)
式中 sum(G)——候選親本中所有遺傳距離的和
G(i,j)——第i株無性系到第j株無性系的遺傳距離
(2)過濾無效樹種選擇。在每次覓食過程中,獲取種子園當(dāng)前基因位的正對近鄰的樹種,為節(jié)約覓食搜索時(shí)間,加快搜索速度,過濾與正對近鄰相同無性系的選擇,避免無效的覓食。
(3)步長選擇。FOA的步長會(huì)影響算法的效率和精度,當(dāng)步長較大時(shí)有利于全局搜索,加快搜索速度,但精度較低;步長較小時(shí)有利于局部搜索,提高算法精度,但搜索速度會(huì)下降。傳統(tǒng)的FOA是采用固定步長進(jìn)行覓食,因此很難平衡算法的效率和精度。為兼顧兩者平衡,本文引入遺傳算法中的基因變異對步長進(jìn)行調(diào)整,在迭代前期可設(shè)置較大的變異因子,有利于加快算法的搜索和收斂速度;在迭代后期減小變異因子,在保證算法精度的同時(shí)也有助于算法跳出局部最優(yōu)。
(4)覓食方向。在迭代前期為了加快搜索速度,將覓食方向定為種子園中濃度最差的基因位,在迭代后期為了避免算法陷入局部最優(yōu),將覓食方向改為隨機(jī)方向。
在求解無性系種子園配置問題時(shí)需要把FOA的連續(xù)解空間對應(yīng)到種子園配置方案,其主要步驟見圖1。
圖1 IFOA算法流程圖Fig.1 Flow chart of IFOA algorithm
3.1.1數(shù)據(jù)材料
從內(nèi)蒙古紅花爾基樟子松國家良種基地中的1代種子園、1.5代種子園、2代種子園遺傳距離的Excel文件中,隨機(jī)選取其中13個(gè)無性系和它們的遺傳距離作為實(shí)驗(yàn)數(shù)據(jù)。
3.1.2參數(shù)設(shè)置
本實(shí)驗(yàn)中,種子園的規(guī)模假設(shè)為9×9(即9行、9列),初始種群大小為25,迭代次數(shù)為800,迭代前期變異因子為0.4,迭代后期變異因子為0.02,在此基礎(chǔ)上將本實(shí)驗(yàn)與傳統(tǒng)的種子園配置方法(CR、RCB)、遺傳算法(GA)以及其他兩個(gè)改進(jìn)果蠅算法:MFFA[16]、ASFOA[18]進(jìn)行了比較,以上所有算法均
采用式(1)作為目標(biāo)函數(shù)。
3.1.3實(shí)驗(yàn)說明
該算法并不只針對數(shù)據(jù)材料中選取的13個(gè)無性系或參數(shù)設(shè)置中的9×9固定大小的種子園,它適用于任意無性系的種類和數(shù)量以及任意規(guī)模的種子園。
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次試驗(yàn)中隨機(jī)選擇6次運(yùn)行結(jié)果,作為各算法的收斂速度,見圖2,IFOA的收斂速度以及效果明顯優(yōu)于其他算法。
3.2.3不同方法設(shè)計(jì)對比
將本文的IFOA算法得出的設(shè)計(jì)方案,分別與效果較好的智能算法GA設(shè)計(jì)方案以及傳統(tǒng)算法CR設(shè)計(jì)方案進(jìn)行對比,如圖3所示。
圖3a存在多處同一無性系的分株作為直接鄰居,如第1行第3個(gè)位置和第4個(gè)位置都是無性系9的分株;圖3b、3c沒有同一無性系的分株作為直接鄰居出現(xiàn);圖3c無性系3和無性系7作為直接鄰居出現(xiàn)的概率明顯高于圖3b,原因是無性系3和無性系7之間的遺傳距離是實(shí)驗(yàn)數(shù)據(jù)中遺傳距離最大的組合,因此會(huì)優(yōu)先交配,有利于雜交育種;圖3c沒有無性系2,因?yàn)闊o性系2與其他無性系之間的遺傳距離都較小,即無性系2與其他無性系都具有較近的遺傳距離,因此不建議栽種。
圖2 各算法收斂情況對比Fig.2 Comparison of convergence of different algorithms
圖3 13個(gè)無性系在9×9種子園的設(shè)計(jì)方案Fig.3 Design of 13 clones in 9×9 orchard
將無性系的分株數(shù)量作為未知條件,基于無性系之間的遺傳距離,不僅考慮了同一親本的分株影響,也考慮了近鄰之間的影響,使得種子園最終能夠達(dá)到較優(yōu)的遺傳效益,通過改進(jìn)型果蠅算法實(shí)現(xiàn)該目標(biāo),并得到較優(yōu)的種子園設(shè)計(jì)方案和親本分株比例。計(jì)算結(jié)果表明IFOA的收斂速度和效果優(yōu)于其他算法,并且IFOA的設(shè)計(jì)方案中同一無性系分株不會(huì)作為近鄰出現(xiàn),當(dāng)存在親本組合的遺傳距離明顯大于其他組合時(shí),即親緣關(guān)系最遠(yuǎn)組合,該組合通常會(huì)作為優(yōu)先交配組合。