非線性回歸模型參數估計的輪回選擇算法

陳庭木 方兆偉 王寶祥 劉艷 邢運高 徐波 徐大勇

摘要：生物學領域規律更多呈現不可線性化的純非線性方程關系，對此類方程的參數準確估計能促進生物學領域的應用數學研究。但純非線性方程最優擬合一直是應用數學沒有完全解決的課題;受水稻輪回選擇育種啟發，以育種進化方法構建新的進化算法，求解復雜非線性方程的擬合問題。結果表明，通過Richard方程、房室模型中的二室模型試算，前者與NIST結果相同，后者明顯優于原文結果;通過奶牛泌乳曲線中的Dijkstra、Wood方程各14組數據最優擬合計算，表明本算法準確有效，且計算結果優于統計軟件SAS。本算法為生物學研究中更復雜的非線性方程擬合及其應用提供了可能，為生物數學的深入發展提供了有力計算工具。

關鍵詞：最優擬合;輪回選擇;非線性回歸

中圖分類號： S11+5文獻標志碼： A

文章編號：1002-1302（2019）18-0253-07

收稿日期：2018-05-04

基金項目：現代農業技術體系建設專項（編號：CARS-01-61）;江蘇省科技計劃重點項目（編號：SBE2017310472）。

作者簡介：陳庭木（1977—），男，安徽蕪湖人，副研究員，主要從事水稻品質遺傳育種與生物統計研究。E-mail：chentingmu@139.com。

通信作者：徐大勇，博士，研究員，主要從事水稻遺傳育種研究。E-mail：xudayong3030@sina.com。

一般線性回歸關系有成熟算法，能應用最小二乘法準確無偏估計回歸參數及其標準差，在農業科學研究領域，變量間更多呈現非線性關系，非線性關系只有經線性代換化為線性方程后才能應用成熟的線性算法，且參數估計只是在線性化條件下最優，以原非線性方程考量不是最優[1-2]。對于不能化為線性方程的非線性方程，稱為純非線性方程。純非線性方程在生物學研究中更為普遍，如生物生長衰老方程、病害發展動態模型[3-7]、混合正態分布方程、復雜房室模型[8-11]、植物光合光響應方程[12-14]、奶牛泌乳曲線機理方程[15]、土壤水分模型[16-18]、農藥殘留降解模型[19-21]等。對于純非線性方程最優擬合難度極大，傳統方法有最速下降法、牛頓法、共軛梯度法，近代的有麥夸法、模擬退火法、遺傳算法、縮張算法、人工神經網絡、蟻群算法等。傳統最優擬合方法存在對初值依賴性強、易陷于局優陷阱[1-2]的缺點，麥夸法對簡單非線性模型初值要求較低[22]，但對于復雜非線性方程（如Richard方程）極難擬合，對于房室模型中一室模型相對容易，對二室模型則會得出不恰當的解[23]。其他算法多是模擬生物自然演化過程的進化類算法。進化算法屬于計算智能，現有進化算法以遺傳算法為代表，但遺傳算法存在早熟收斂的問題，至今未能很好克服[24]。縮張算法是一種優秀的最優擬合算法，不依賴導數支持，在已報道實例中均能達到全局最優，但對于參數過多，如參數20個以上，則按最少試算點計算，每次試算量達320個試算點，計算難以實現。故應當進一步研究開發新的算法，以克服遺傳算法與縮張算法的固有缺陷。

1?材料與方法

1.1?數據來源

1.1.1?Richard方程擬合

美國國家標準網站（NIST）提供了Richard方程的最優擬合計算案例，殘差平方和（http：//www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml）由15對數據組成，最小離回歸平方和Qe=8 786.404 908 0。

yi=A（1+Be-kti）1/N。（1）

式中：yi為第i個觀察值;A為生長極限值;B為初始生長量參數;k為生長速度參數;ti為第i個時刻;N為曲線形狀參數。

Qe=∑n-1i=0（yi-yei）2。（2）

式中：Qe為殘差平方和;yi為第i個觀察值;yei為第i個預測值。

1.1.2?二室模型

二室模型共5個參數，且滿足Km>K2>K4。摘自袁志發《多元統計分析》中二室模型擬合實例演算[22]，文獻[25-26]以麥夸法求解，初值以退層法求得。

二室模型：A2（e-k2t-e-kmt）+A4（e-k4t-e-kmt）。（3）

1.1.3?奶牛泌乳方程Dijkstra與Wood方程[25-26]

分別用下面2個公式表示。

y=Aeb（1-e-ct）/c-dt;（4）

y=Atbe-ct。（5）

式中：各字母含義見文獻[25-26]。

以2種數學模型分別對羅清堯、能本海等《中國荷斯坦奶牛第二泌乳期泌乳曲線模型的研究》《中國荷斯坦奶牛第三泌乳期泌乳曲線模型的研究》文中第2、第3胎的各7組數據重算。

1.2?輪回選擇算法

輪回選擇是經典的育種方法，工作重點是選種、鑒定與自由互交。輪回選擇在有限的遺傳背景中，能最大程度地淘汰不良基因，聚合優良基因，以產生更高產量、抗性與品質的新種質。在輪回選擇時，每輪選擇大量選種材料，鑒定出10%的優良選種材料，作為自由互交親本，親本間自由雜交產生下輪選種材料。如此周而復始，不斷提高親本遺傳表現水平，多代選擇后必會獲得一群高水平的個體。

筆者受輪回選擇育種程序啟發，將殘差平方和作為目標函數值，看作遺傳表現值，約束條件不滿足看作致死基因（只選擇存活個體），各自變量看作染色體，染色體以實數編碼，先按各參數給定區間隨機生成染色體組，構成若干親本，親本間進行自由雜交、自由突變產生一定量存活選種群體，對存活選種群體依目標函數值排序，選擇若干優良選種材料作下輪親本，親本間再自由雜交、自由突變產生下代選種群體，如此周而復始，直到親本間沒有遺傳變異與或目標函數值收斂到迭代精度要求為止。經歷多個連續世代進化，收斂達到精度要求，稱為1輪輪回選擇進化，此進化算法稱為輪回選擇算法。

當對參數取值區間完全不知情，可以設定1個很寬的區間，首輪（第1輪）進化結束，以最優個體自變量（染色體）取值為中心，以其絕對值為鄰域半徑，進行第2輪輪回選擇進化計算，之后每輪輪回選擇進化計算，均以最優個體自變量（染色體值）取值為中心，將區間縮小到0.618倍，重新生成隨機親本，并加入上輪最優個體，重新作輪回選擇進化計算，一般進行10輪計算方能達到全局最優解。如上組織多輪輪回選擇進化計算，放寬初始區間要求及提高計算精度的算法稱為多重輪回選擇算法。

1.2.1?基本輪回選擇算法

輪回選擇算法基于親本自由雜交及自由突變產生選種群體，再從選種群體中優選存活個體組成新一輪親本群體，每一個由初始親本群體形成下輪親本群體的進化過程稱為1個進化世代。不同進化世代間，只要上代親本選種群體存在差異且有優劣區別，進化就會發生，反之進化則會停止。每輪進化由下列算子構成：初始親本生成算子、雜交算子、突變算子、選擇算子、迭代終止算子。

不斷地應用雜交算子、突變算子、選擇算子、迭代終止算子進行世代更替，當世代間進化達到收斂要求，則完成了一輪輪回選擇計算。計算實踐表明，1輪輪回選擇計算結果一般不差于遺傳算法，一般表現更優。

針對本算法，以C++編程語言，以面向對象的編程方法設計RSBase類，將常用參數設計成類的屬性，各算子設計成類的私有方法，參數輸入與輸出方法設計為公有方法，以利于使用者改寫與調用。類中將Fx（目標函數定義）、Gx（約束函數定義）設定為純虛函數，每次調用必須根據使用者特有的統計模型改寫，同樣，參數輸出方法Para也設計成了虛函數，方便改寫。另外設計了函數求極值及積分的公有方法，利于對使用者定義的復雜數學模型求函數極值點、極值及區間積分。軟件著作權《連農最優生長擬合類庫軟件V1.0》（登記號：2147408）包含了本算法全部C++代碼。算法流程詳見圖1。

1.2.2?多重輪回選擇算法

一輪輪回選擇計算多由于區間不恰當或區間過大，一般不能達到全局最優解，如要求更高計算精度，可在上輪計算結果基礎上精算。每輪初始親本生成邊界值由空間收斂算子給出，算子算法描述如下：

tmp=0.618repeat

repeatrepeat+1

if（tmp<0.001）tmp0.001

xa=optX-tmp×|optX|

xb=optX+tmp×|optX|。

對每個參數獨立計算取值區間，由xa、xb計算產生初始親本群體及表征突變程度上限。收斂倍率0.618是筆者經過數百次試算由經驗給出，一般取值區間（0.5，1）。repeat初始值為0，optX為參數上輪最優估計值，tmp最小值暫定為0.001，最佳值有待深入研究確定。

每輪最優親本由上輪最優選種個體遺傳，能保障輪次間收斂與持續進化，參數取值區間的縮小，能更大程度地減少早熟收斂的可能性。多重輪回選擇算法流程見圖2。

1.2.3?輪回選擇算法屬性構造與各算子設計

1.2.3.1?屬性構造

本算法以C++語言設計，程序執行前，要求定義一些基本的參數，以利程序安排必要的內存空間用于計算。重要參數如下：

Nx為參數變量個數;Repeat為多重輪回選擇計算輪數計數器;nP、nS為親本與選種群體規模;iG、iS、maxG為進化世代計數、選種個體計數及最大進化代數;isint為整型變量標志指針，規定變量取實數還是整數;xa、xb參數取值區間指針，初始群體由之計算產生，突變程度也由之表征;optX、optValue最優個體染色體取值指針及其目標函數值，作為算法輸出;Parent、valueP親本群體染色體取值及對應目標函數值指針;Select、valueS選種群體染色體取值及對應目標函數值指針。

1.2.3.2?初始群體構造算子

第i個變量依據（xa，i，xb，i）區間隨機生成親本參數，如整型標志變量isinti=0將相應變量取整，每個變量的區間、整型約束獨立處理。初始群體由存活個體構成。初始群體的構造，不像其他算法那樣過度依賴初始值，給研究者一個寬松的條件，之后在第1輪最優結果獲得后，再據其修正取值區間，區間長為最優個體相應自變量取值絕對值的2倍，以防初始區間不恰當造成的隨機生成樣本分布不均勻，之后每輪區間長度縮小到上輪0.618倍，當區間長達到原區間長的0.001倍時不再縮小。每輪區間收斂能保證最優解區間的收斂，減少進化難度，區間愈小，進化愈快。

1.2.3.3?雜交算子

每條染色體獨立雜交，每次雜交隨機選擇2個不同的親本，以[0，1]間隨機數α為因子一次雜交形成2個選種個體。有整型約束的變量變換為整型。雜交如同輪回選擇中的自由互交，個體都由自由交配形成，保障了遺傳多樣性。雜交采用實數形式表示，由式（6）、式（7）計算：

x1=αPi+（1-α）Pj;（6）

x2=（1-α）Pi+αPj。（7）

1.2.3.4?誘變算子?因為誘變產生變異及雜交產生變異是本算法的2種變異方式，誘變概率過低會導致雜交變異產生的變異個體顯著多于誘變產生的變異個體，為相對均衡起見，不能將誘變概率定太低。αi∈[0，1]間隨機數，用于計算第i個變量突變程度，βi∈[0，1]間隨機數，表示第i個自變量是否突變，表征誘變發生的概率。當βi<0.1時進行突變處理，否則不作突變處理，整體突變概率很高。計算實踐表明，突變在進化過程中起重要作用。如果αi<0.5，采用式（8）。

x=P-αi[（xi）max-xai]。（8）

式中：（xi）max為當前世代親本群體第i個自變量的最大值。如果αi≥0.5，采用式（9）。

x=P+αi[xbi-（xi）max]。（9）

式中：a、b表示屬性值的左、右邊界值。

1.2.3.5?選擇算子

選擇分為2個部分。一部分選擇是致死選擇，即約束函數Gx=0，為致死個體，不能參與進化計算;另一部分選擇對存活個體選優汰劣，由選種群體的排序功能完成，優選nP個個體加入下輪親本群體。選擇算子優選出的解總是可行解，這與運籌學中的外點法不同，同時又避免了內點法初始可行點難以選擇的矛盾。迭代過程中不會因為函數性態差而產生大幅度偏差。

1.2.3.6?空間收斂算子

當解空間過大時，會造成進化計算獲取的采樣點不一定能覆蓋最優解空間，容易造成早熟收斂，如對不包含解空間的空間進行割舍，能集中計算能力尋找最優解。具體方法見“1.1.2”節。

1.2.3.7?迭代終止算子

迭代終止算子是判斷計算能否終止的方法，當進化代數達到上限時，終止計算;或親本沒有遺傳變異時，終止計算;或親本群體目標函數值變異標準差小于指定精度值，提前結束本輪計算。當群體變異存在，但目標函數值已經基本沒有變化時，應提前終止;如有整數約束時，常存在多個最優解。

1.3?經典算法介紹

高斯-牛頓法應用各一階導數構造一個正交矩陣，建立一個最優搜索方向，當目標函數全局性態單一、有類似正定二次型的曲面形態且初始搜索點處對應正交矩陣可逆，則計算過程很快，速度優于所有現代算法。

麥夸法是在高斯-牛頓法基礎上發展起來的一種近代算法。當初始搜索點處對應正交矩陣不可逆，通過將正交矩陣主對角線元素加上一足夠小正數，可以改變矩陣性態，可以適當放寬高斯-牛頓法的初始值條件，計算量相對其成倍提高。

牛頓法以初始點處黑塞矩陣計算優化點，當黑塞矩陣可逆時，迭代計算非常快，如目標函數為正定二次型可一次迭代成功，當黑塞矩陣不可逆時，類似于麥夸法，將黑塞矩陣主對角線元素加上一足夠小正數，此法稱為修正牛頓法。牛頓法與修正牛頓法要計算各種二階導數，計算復雜，不利于編程計算。二者與前述兩法相比，存在同樣的問題，當目標函數性態很差時，存在相當多的局部最優，則計算很難達到全局最優，有的甚至沒有可行解。

縮張算法是由揚州大學顧世梁等提出的最優化算法，通過擴張步與收縮步組成一輪循環，初始值由隨機數生成，經過多輪的縮張循環，可以較快地計算得到全局最優解，在低維問題處理上本算法目前是最優方法，但面對高維問題時，其試算點數過多，計算時間相當長，計算難以實現。

1.4?各算法間比較指標

1.4.1?誤差方差

各算法比較最關鍵指標為殘差平方和，好的算法獲得估計參數更加適合所研究模型，應變量變異平方和中誤差平方和占比最小化。均方誤差（mean squared error，MSE）是衡量“平均誤差”的一種較方便的方法，表征實際值與預測值的差異程度，可由殘差平方和計算得到。均方根誤差（root mean squared error，RMSE）是均方誤差的算術平方根，也稱作標準誤差（standard error）。

1.4.2?AIC指標

AIC是赤池弘次提出的衡量統計模型擬合優良性的一種標準，可以權衡所估計模型的復雜度和此模型擬合數據的優良性。

AIC=2k+nln（Qe/n）。（10）

式中：k是參數的數量;n為觀察數;Qe為殘差平方和。

1.4.3?計算時間

計算時間是衡量算法效率的重要指標，但算法有效性往往比算法效率更重要，如牛頓法是非線性最優化問題的最快算法，但計算往往不能得到正確結果，現代遺傳算法計算時間遠多于牛頓法，但獲得了廣泛的應用。縮張算法計算時間也遠多于牛頓法，但它是處理低維非線性問題的最佳算法。可見，計算時間在算法間比較意義相對較小，更重要的是比較計算可行性、準確性與魯棒性。

1.4.4?對初值的依賴性?好的算法不僅要有精確的計算結果，還不能依賴初值的精確性，如牛頓法對初值要求較高，只有好的初值才會收斂到精確結果，麥夸法相對于牛頓法則對初值依賴性有所下降，但不當的初值還會造成不收斂[22-23]。縮張算法對初值基本沒有要求，本研究算法對初值也沒有要求，對初值沒有依賴的算法稱為沒有初值依賴性。

2?結果與分析

2.1?Richard方程最優擬合結果

應用輪回選擇算法經歷11輪計算，第0輪經歷10代進化，第1輪經歷17代進化，…，第10輪經歷7代進化（表1），殘差平方和Qe達到本算法最優值8 786.404 913 187 29，與NIST網站公布的全球最優結果8 786.404 908 0前8位有效數字完全相同，AIC分別為49.515 799 32、49.515 799 31。經多次計算均達到相近計算結果（前9位有效數字一直相同），本次計算第6輪已達到計算精度要求，之后RMSE變化很小，且參數估計值變化也相當小，有的輪次間不能產生進化。本次嘗試隨機生成初值再作麥夸法求解，當隨機點數達5 000個時，其最優結果仍不及本算法。本算法屬性值構造：親本群體規模100，每代雜交規模10 000，誘變規模10 000，迭代收斂精度要求1.0×10-10。A初始區間[0，1 000]，B初始區間[0，1 000]，K初始區間[0，10]，N初始區間[0，10]。建議10輪為佳，如提前進入全局最優值，之后輪次的進化世代數很少，進化很快。平均每個數據集單個模型計算用時不超過120 s。

擬合殘差總體為可接受水平，擬合決定系數0.991 8，相關系數0.995 9。殘差分析見表2。

2.2?房室模型中二室模型擬合

二室模型共5個參數，且滿足要求Km>K2>K4。對于有約束的擬合，傳統非線性擬合難度極大。給定初值不當，不能產生符合約束的解，而初值選擇難度又很大。輪回選擇算法可將約束作為致死因子，在存活個體中優選親本雜交、突變，產生新的選種群體，遴選最優個體，求解最優擬合參數。以式（3）模型與數據試算（表3）。本算法屬性值構造：親本群體規模100，每代雜交規模10 000，誘變規模10 000，迭代收斂精度要求1.0×10-10。初始區間A2∈[0，5 000]、K2∈[0，10]、Km∈[0，50]、A4∈[0，50]、K4∈[0，50]，經11輪計算，得 Qe=6 731.829 722 431 19（AIC=47.548 071 21），明顯低于文獻[22]中的最優結果6 761.911 38（AIC=47.575 180 15），且完全符合約束條件。

以上擬合數據說明，K4=0，不適合二室模型擬合，更適合一室模型擬合，但原文中K4明顯不為0，不能說明一室模型不合適;A4與原文差異較大，其他3個參數差異較小。本算法計算更準確。同Richard模型計算相似，第6輪基本不再有大的進化。最優擬合殘差見表4，決定系數0.996 4，相關系數0.998 2。

2.3?最優數學模型選擇

奶牛泌乳方程，國內外已報道有100多種，以“1.1.3”節模型與數據計算。

第2胎泌乳數據擬合殘差平方和（表5）作模型、組別兩因素無重復方差分析，表明模型間差異顯著（F值9.083 321 499，無差別顯著概率0.023 583 05），Dijkstra顯著優于Wood模型。第3胎泌乳數據擬合殘差平方和（表5）作模型、組別兩因素無重復方差分析，也表明模型間差異顯著（F值11.093 4，無差別顯著概率0.015 793），Dijkstra也顯著優于Wood模型。這與文獻[25-26]結論完全相反，對比表7數據發現，Wood模型參數文獻[25-26]SAS計算結果與本研究輪回選擇計算結果相似，但本算法所有殘差平方和均小于SAS結果，方差分析表明，2種方法間差異F值4.282 694 075，無差別概率0.058 986 506，接近顯著水平，2種算法殘差平方和差異越大，則參數估計值差異越大，表明輪回選擇算法在處理簡單非線性模型方面優于SAS。表6顯示，Dijkstra模型數據差異較大，少數組別相近，多數組別差異大，且擬合殘差14組數據均小于SAS結果，表明輪回選擇算法在處理復雜非線性模型方面也優于SAS軟件系統的非線性分析模塊（一般采用高斯-牛頓法）。文獻[25-26]14組試驗數據分別用2種模型計算，形成28組擬合數據，擬合殘差平方和Qe均表明本算法較SAS非線性擬合優秀，但AIC結果2個模型間差別不明顯，可能AIC指標結合了參數個數信息，實質上泌乳數據需要更精確的預測結果，以Qe比較更為恰當。

3?討論與結論

3.1?輪回選擇算法的優點

3.1.1?魯棒性

輪回選擇算法采用親本自由雜交、自由突變產生選種群體，從選種群體中優選個體組成新的親本群體，不斷地選優汰劣，實現全局尋優計算。本算法不采用導數計算，不要求目標函數及約束函數連續可導，另可以針對整數約束求解，不受一般最優化問題的實數連續約束的制約。當模型添加了區間限制（如非負），不等式約束限制（如房室模型中的二室、三室模型）等時，有特別大的靈活性。

3.1.2?適合解混合非線性模型與高難度生物學模型

如非線性模型中部分自變量存在整數約束，則稱為混合非線性模型，解此類問題，傳統優化算法無能為力。運籌學中采用罰函數法中的外點法求解，初始值易找，但所得常為不可行解，內點法初始點不易尋找，所得解均為可行解，但可行域不一定唯一，不能保證解為所有可行域中的最優解，實踐應用不便，如采用輪回選擇算法，易于實現，如將Richard方程中N要求取整數，求解中僅加整數標志即可。

由于非線性模型擬合一般采用導數支持的傳統算法，當模型復雜時，會在最優解附近存在很多局部優化陷阱，擬合難度極大，如Richard方程，在美國NIST網站中將其歸為高難度數組。植物生長過程不是一個單調生長過程，還包括衰老過程，僅用生長方程來擬合，不能解析衰老過程，衰老過程還干擾生長過程的參數擬合。如設計生長衰老方程，如式（11），將生長與衰老過程相疊加，則可以解析生長過程與衰老過程，還能找到二者平衡點，揚州大學孫成明在博士論文中提及此想法[27]，但一直未有應用報道，主要原因是參數擬合難度遠高于Richard方程，筆者研究發現，輪回選擇算法可以高效準確擬合此類方程，將另文報道。

W=A11+B1e-kt-A21+e-k（t-t0）;（11）

A2（e-k2t-e-kmt）;（12）

A2（e-k2t-e-kmt）N。（13）

式（12）方程用于描述生物對藥物的吸收清除過程，稱為一室模型，在應用中一室模型擬合精度一般不高，但給方程括號項加一指數成分后變為四參數方程，形如式（13），擬合精度顯著提高，也極大地提高了擬合難度，用輪回選擇算法也能準確擬合，將另文報道。此類復雜的生物學方程還有很多，意義很大，難度很高，限制了在生物學科研中的應用，輪回選擇算法的成功將拓展對復雜生物學模型的開發與應用，促進生物學理論研究的深入。

3.1.3?每輪重新初始化，防止了解樣本分布不均勻、早熟收斂問題

遺傳算法與模擬退火算法在本算法提出前早有應用，但其早熟收斂缺陷不易克服，主要是進化計算過程中對最優參數空間的不自覺舍棄[22]。輪回選擇算法因采用了類似輪回選擇育種中自由互交的變異方式，單輪計算出現早熟收斂的概率較低，采用多輪優化技術，每輪優化均對最優解空間作分割，區間長縮小到0.618倍，當早期輪次未計算出全局最優解時，可在以后的輪次出現全局最優解，在足夠的計算輪后總能找到全局最優解。如提高計算精度（高精度計算支持）還可以突破現有最優記錄。從表1和表3結果可知，隨著計算輪次的增加，進化增益很小，一般沒必要超過10輪，完全可以處理一般科研應用問題的求解。

3.1.4?隱并行性

如同遺傳算法，本算法避開了最優方向及最優步長的搜索，也有計算的隱并行性，可實現多線程的并行計算，大幅提高計算時間效率。

3.2?算法優化

輪回選擇算法計算，對不同的模型，最優算法參數不同，采用相同的參數也完全可以計算出各種類型模型的全局最優解，但計算效率有區別，針對不同的模型可以探討最優參數配置問題，如以全局最優解解出計算時間為目標函數值進行復雜的非線性規劃，可以計算最優參數配置問題，此問題有待深入研究。本算法有待改進，促進算法的完善與推廣。

參考文獻：

[1]顧世梁，惠大豐，莫惠棟. 非線性方程最優擬合的縮張算法[J]. 作物學報，1998，24（5）：513-519.

[2]顧世梁，徐辰武，蒯建敏. 用縮張算法最優擬合非線性方程的檢驗[J]. 揚州大學學報（自然科學版），2001，1（3）：16-19.

[3]劉維紅，林茂松，李紅梅，等. 人工接種測定水稻干尖線蟲在水稻上的病害發展動態[J]. 中國農業科學，2007，40（12）：2734-2740.

[4]劉維紅. 水稻干尖線蟲在“小穗頭”上的危害及人工接種測定病害發展動態[D]. 南京：南京農業大學，2007.

[5]彭懷俊，顧?鋼，紀成燦，等. 烤煙根系土壤中青枯病菌動態與田間病害發生發展的關系[J]. 湖南農業大學學報（自然科學版），2005，31（4）：384-387.

[6]曹?靜. 設施條件下馬鈴薯晚疫病流行主導因素及病害防治的初步研究[D]. 保定：河北農業大學，2002.

[7]柏自琴. 中國柑橘黃龍病發生動態及其病原菌亞洲種群分化研究[D]. 重慶：西南大學，2012.

[8]曾文藝，孫曉穎. 藥物動力學房室模型的改進[J]. 北京師范大學學報（自然科學版），2012，48（1）：6-10.

[9]王曉麗，趙?晶. 傳染病房室模型的建立及求解方法[J]. 山東輕工業學院學報（自然科學版），2006，20（3）：88-90.

[10]何紹雄，曹鑒萍，黃團華. 非線性房室模型參數計算法[J]. 數值計算與計算機應用，1986（3）：147-152.

[11]楊博文，劉?萍. 藥物動力學二房室模型[J]. 哈爾濱師范大學（自然科學學報），2016，32（1）：13-15，66.

[12]韓?剛，趙?忠. 不同土壤水分下4種沙生灌木的光合光響應特性[J]. 生態學報，2010，30（15）：4019-4026.

[13]閆小紅，尹建華，段世華，等. 四種水稻品種的光合光響應曲線及其模型擬合[J]. 生態學雜志，2013，32（3）：604-610.

[14]王榮榮，夏江寶，楊吉華，等. 貝殼砂生境干旱脅迫下杠柳葉片光合光響應模型比較[J]. 植物生態學報，2013，37（2）：111-121.

[15]賈先波，陳仕毅，王?杰，等. 中國荷斯坦牛泌乳曲線數學模型擬合及應用[J]. 南京農業大學學報，2016，39（1）：133-138.

[16]李洪文，高煥文. 保護性耕作土壤水分模型[J]. 中國農業大學學報，1996，1（2）：25-30.

[17]于?浕. 土壤水分特征曲線Gardner模型參數的預報模型研究[D]. 太原：太原理工大學，2017.

[18]張?巖，朱?巖，張建軍， 等. 林地土壤水分模型SWUF在晉西黃土高原的適用性[J]. 林業科學，2012，48（5）：8-14.

[19]吳?鵬，秦智偉，周秀艷，等. 蔬菜農藥殘留研究進展[J]. 東北農業大學學報，2011，42（1）：138-144.

[20]陳振德，陳雪輝，馮明祥，等. 毒死蜱在菠菜中的殘留動態研究[J]. 農業環境科學學報，2005，24（4）：728-731.

[21]張韓杰，閆艷春. 農藥殘留及微生物在農藥降解中的應用與展望[J]. 湖北植保，2004（1）：31-35.

[22]袁志發，周靜芋. 多元統計分析[M]. 北京：科學出版社，2002：247-248.

[23]周靜芋，宋世德，袁志發，等. 用麥夸法進行室分析曲線擬合[J]. 西北農業大學學報，1996，24（1）：75-78.

[24]范青武. 遺傳算法動態演化過程的研究[D]. 北京：北京工業大學，2007.

[25]羅清堯，熊本海，馬?毅，等. 中國荷斯坦奶牛第二泌乳期泌乳曲線模型的研究[J]. 中國農業科學，2010，43（23）：4910-4916.

[26]熊本海，馬?毅，羅清堯，等. 中國荷斯坦奶牛第三泌乳期泌乳曲線模型的研究[J]. 中國農業科學，2011，44（2）：402-408.

[27]孫成明. FACE水稻生長發育模擬模型研究[D]. 揚州：揚州大學，2006.