(浙江工業(yè)大學 機械工程學院，浙江 杭州 310014)

設施農(nóng)業(yè)是現(xiàn)代農(nóng)業(yè)的基礎，而溫室栽培是其重要組成部分，其中番茄以其廣闊的市場前景和較大的需求，是溫室中主要栽培蔬菜之一。以往溫室種植者擁有的反季蔬菜提供能力使得其種植的溫室蔬菜在市場上擁有較強的競爭力。然而近年來溫室的數(shù)量逐漸增多，國外蔬菜的進口量越來越大，市場的競爭越來越激烈，用科學的方法對溫室栽培的過程進行管理，增加溫室栽培的經(jīng)濟效益成為在競爭中站穩(wěn)腳跟的唯一方法。作物產(chǎn)量是溫室生產(chǎn)中衡量經(jīng)濟效益的關(guān)鍵指標。因此，建立精確的番茄產(chǎn)量預測模型具有重要意義。

作物產(chǎn)量可以通過干物質(zhì)分配模型計算，干物質(zhì)分配模型主要有經(jīng)驗模型與機理模型。國內(nèi)學者對經(jīng)驗模型的研究較多，倪紀恒等[1]利用累積輻熱積擬合分配指數(shù)，并結(jié)合采收指數(shù)來預測溫室番茄的產(chǎn)量。齊維強等[2]采用積溫法研究溫室番茄的生長發(fā)育模型。然而這類經(jīng)驗模型在建立溫室番茄生長發(fā)育模型時主要依靠經(jīng)驗，缺乏機理性。國外相關(guān)的研究開展較早，目前能夠較為準確地描述番茄生長的產(chǎn)量模型有TOMSIM(Tomato simulator)番茄生長發(fā)育動態(tài)模型[3]和Vanthoor模型[4]。這類機理模型能夠模擬各器官的生長過程，描述作物的生長狀態(tài)，如當前果實數(shù)量、葉面積指數(shù)及干物質(zhì)分配比例等，同時還能夠反映室內(nèi)氣候?qū)Ψ焉L的影響，如溫度、CO2含量和光照強度的影響。然而機理模型一般都有一個通用的缺點，雖然能夠精確地模擬作物的生長過程，但是由于含有大量的作物生長狀態(tài)，很多參數(shù)難以計算，在溫室優(yōu)化控制中無法直接使用。因此出現(xiàn)了簡化TOMGRO(Tomato growth)番茄生長模型[5-6]。簡化TOMGRO模型機理性強，在降低模型復雜度的同時還不失精確性，能夠準確預測未來產(chǎn)量變化，但是在實際應用時需要根據(jù)實際種植環(huán)境對描述番茄生長狀態(tài)的部分參數(shù)進行辨識[7]。參數(shù)辨識依據(jù)實驗數(shù)據(jù)和所建立的模型，結(jié)合優(yōu)化算法來確定模型中待辨識的參數(shù)，使模擬值最大程度上擬合實際值。魯建廈等利用改進粒子群算法對混流混合車間生產(chǎn)調(diào)度問題[8]和裝配線平衡[9]進行優(yōu)化。王衛(wèi)紅等[10]提出PCA-FOA-SVR方法對股票價格進行預測。張娟等[11]運用RBF神經(jīng)網(wǎng)絡和遺傳算法建立了番茄莖干的生長模型。文歡等[12]利用灰色神經(jīng)網(wǎng)絡模型對加工番茄產(chǎn)量進行預測。Vazquez-Cruz等[13]用遺傳算法校準TOMGRO模型。然而神經(jīng)網(wǎng)絡模型需要采集大量數(shù)據(jù)進行訓練，遺傳算法中遺傳操作參數(shù)的選擇大部分依靠經(jīng)驗，算法搜索速度慢。因此需要一種快速精確的優(yōu)化算法對簡化TOMGRO模型進行參數(shù)辨識。結(jié)合PSO算法收斂快速[14]和SA算法不易陷入局部最優(yōu)[15]的特性，筆者提出了IPSO-SA算法。通過分析算法求解待辨識參數(shù)的過程以及模型產(chǎn)量預測的精確度，比較PSO，SA和IPAO-SA算法的優(yōu)劣，并得到適合本地種植條件的番茄產(chǎn)量預測模型，為溫室番茄的生產(chǎn)提供參考。

1 材料與方法

1.1 試驗材料

試驗溫室位于上海市崇明示范基地，主要用于栽培無限生長型番茄。溫室控制系統(tǒng)自帶的Priva氣象站、CO2檢測儀和光照傳感器等儀器每5 min進行一次數(shù)據(jù)采集，對溫室進行全方位監(jiān)控[16]。無限生長型番茄單個果實質(zhì)量不大，由于生產(chǎn)周期長，具有高產(chǎn)量的特點。根據(jù)崇明溫室的種植經(jīng)驗，串番茄單果鮮重質(zhì)量為100～200 g，根據(jù)一次收獲的果實總產(chǎn)量和總數(shù)，可以估算出單果果實鮮重質(zhì)量為160 g左右。番茄的生長周期總天數(shù)為300 d，集中育苗時間約為30 d，溫室種植的時間約270 d，結(jié)合崇明溫室實際的種植情況及已有數(shù)據(jù)，選擇2015年9月30日—2016年5月12日(共32 周，224 d)為番茄在溫室內(nèi)的種植周期。Vanthoor模型是一種機理比較充分的番茄生長模型，能夠比較精確地模擬番茄的生長發(fā)育過程，已經(jīng)在歐洲和北美等地區(qū)得到了廣泛驗證。本研究使用的數(shù)據(jù)包括上海崇明試驗溫室2015—2017 年兩個種植周期的室內(nèi)溫度數(shù)據(jù)、CO2濃度數(shù)據(jù)、光照強度數(shù)據(jù)和番茄生長數(shù)據(jù)，其中番茄幼苗移植到溫室中時的初始數(shù)據(jù)如表1所示，番茄實際產(chǎn)量數(shù)據(jù)由Vanthoor模型仿真得到。

表1 溫室番茄初始狀態(tài)Table 1 Greenhouse tomato initial state

1.2 番茄產(chǎn)量預測模型構(gòu)建方法

筆者采用機理模型與參數(shù)辨識方法[17]，建立番茄產(chǎn)量預測模型，構(gòu)建方法如圖1所示。番茄的生長模型采用簡化TOMGRO模型，結(jié)合前人對TOMGRO模型的研究，提取模型中的待辨識參數(shù)，得到待辨識的番茄產(chǎn)量預測模型。將溫室內(nèi)的溫度、CO2濃度、光照強度以及番茄實際產(chǎn)量數(shù)據(jù)輸入待辨識的模型，通過比較模擬值與實際產(chǎn)量值，采用優(yōu)化算法計算求解不確定參數(shù)。利用不同日期的室內(nèi)外環(huán)境數(shù)據(jù)和產(chǎn)量值對優(yōu)化后的番茄產(chǎn)量模型進行有效性校驗，若模型驗證合理則可用來預測番茄產(chǎn)量，否則重新進行上述步驟，直到模型驗證合理。

圖1 番茄產(chǎn)量預測模型流程圖Fig.1 Tomato yield forecasting model flow chart

2 番茄生長模型

在簡化TOMGRO模型中，主要用果實干物質(zhì)總量WF和成熟果實干物質(zhì)總量WM兩個變量來描述果實生長狀態(tài)。地表干物質(zhì)總量W、植株莖節(jié)數(shù)N和葉面積指數(shù)LAI被用來描述莖、葉等營養(yǎng)性器官的生長狀態(tài)。簡化TOMGRO模型通過植株莖節(jié)發(fā)育率、光合作用速率和呼吸速率計算植株地表干物質(zhì)日生長率和不同器官的干物質(zhì)分配比例，并引入新的函數(shù)用于模擬計算葉面積增長、生長性器官、果實和成熟的干物質(zhì)分配比例。

溫室無限生長型番茄生長發(fā)育過程在簡化TOMGRO模型中分為生長期、結(jié)果期和成熟期。當植株上的莖節(jié)數(shù)達到NFF時，果實進入結(jié)果期。之后當莖節(jié)數(shù)增長KF時，認為果實進入成熟期。植株莖節(jié)發(fā)育率dN/dt忽略了CO2濃度對植株莖節(jié)發(fā)育的影響，具體模型描述為

(1)

式中：Nmax為最大莖節(jié)生長率，在最適宜溫度下一般為每天0.5 節(jié)；FN(T)為非最適宜溫度下莖節(jié)發(fā)育率的下降趨勢函數(shù)；T為室內(nèi)溫度，℃；Th為當天第h個小時的室內(nèi)平均溫度，℃。

番茄生長時的葉面積指數(shù)隨時間的變化率可以表示為

(2)

式中：ρ為植株密度，株/m2；δ為相對于莖節(jié)的葉面積最大增長量，m2/莖節(jié)；β為莖節(jié)數(shù)在式中的相關(guān)系數(shù)；Nb為相關(guān)莖節(jié)數(shù)。

果實在一個特定的時間或者主莖節(jié)數(shù)達到NFF時開始發(fā)育，NFF指第一個果實直徑超過10 mm時植株的主莖節(jié)數(shù)，WF在模型中可以用來表示果實開始發(fā)育后植株地上部分向果實輸送的生物量。在求WF之前，首先要計算植株地上部分生物量增長率GRnet，即

GRnet=E·(A-Rm)·[1-fR(N)]

(3)

式中：A是日光合作用總量，以CH2O計，g/(m2·d)；Rm是日呼吸作用總量，以CH2O計，g/(m2·d)；函數(shù)fR(N)為每日地表生物量分配至根的部分；E為光合同化速率至干物質(zhì)積累速率的轉(zhuǎn)換系數(shù)。

從莖節(jié)數(shù)達到NFF后每天的果實干物質(zhì)增長量dWF/dt可以表示為

(4)

式中：函數(shù)fF(Td)表示正常的日平均溫度Td對地上干物質(zhì)總增長量分配至果實比例的影響函數(shù)；函數(shù)g(Tdaytime)表示白天平均溫度超過26 ℃的高溫環(huán)境對果實生長的影響，與Tdaytime有關(guān)；αF表示地表總生物量分配至果實的最大比例；θ表示營養(yǎng)性器官與果實生長率的轉(zhuǎn)換系數(shù)。

葉面積指數(shù)達到最大值LAImax后，植物葉片會衰老或者被剪除，所以在計算地表總生物量時應該去除被剪除或者衰老的那一部分，地表植株干物質(zhì)增長率的公式可表示為

(5)

式中P1表示葉面積指數(shù)LAI達到LAImax后因為日常的剪除葉面或植株衰老而產(chǎn)生的葉片干重損失率，g/莖節(jié)。先假設P1在LAI未達到最大時為0，在LAI達到最大后P1恒等于2 g/莖節(jié)。

高溫條件下果實的生長率會下降，但植株內(nèi)有較多碳水化合物未轉(zhuǎn)化。假設營養(yǎng)性器官的最大生長率Vmax和最大地表總干物質(zhì)增長率(dW/dt)max，用來限制果實的生長，公式為

(6)

移栽后，番茄在某天的植株地表干物質(zhì)生長率dW/dt取上述兩式計算結(jié)果的較小值。

簡化模型中采用從第一顆果實出現(xiàn)到其成熟這段時間內(nèi)的植株莖節(jié)差KF表示發(fā)育期時間延遲，為了描述第一個果實成熟后的所有果實發(fā)育情況，引入函數(shù)DF(Td)來表示所有果實成熟前的平均發(fā)育率，即

(7)

式中DF(Td)是日平均氣溫在第一顆果實成熟后對后續(xù)果實的發(fā)育影響函數(shù)。

結(jié)合以上公式，可得到溫室番茄產(chǎn)量WM與時間t之間的關(guān)系為

WM(t)=f(t,P(t))

(8)

式中P(t)是簡化TOMGRO模型中的環(huán)境變量，t是時間。

在種植地區(qū)、溫室小氣候條件、番茄品種發(fā)生改變時，簡化TOMGRO模型中需要重新估計的參數(shù)包括用于描述成熟果實WM的兩個參數(shù)：第一顆果實出現(xiàn)到成熟植株上增長的莖節(jié)數(shù)KF和成熟果實最大發(fā)育率DF max，以及用于描述番茄果實總干重W的4 個參數(shù)：地表總生物量分配至果實的最大比例αF、營養(yǎng)性器官與果實生長率的轉(zhuǎn)換系數(shù)θ、出現(xiàn)第一顆果實時植株上的莖節(jié)數(shù)NFF和果實發(fā)生敗育現(xiàn)象的白天氣溫臨界值TCRIT。

3 改進型PSO-SA算法

粒子群算法[18]是一種通過研究鳥類和魚類群體性覓食行為并結(jié)合演化計算理論提出的具有群體智能的全局優(yōu)化算法。PSO算子在上一代粒子的基礎上，更新速度和位置信息，得到新種群。其下一代的粒子的位置和速度信息向量分別為

Vi(t+1)=wi(t)Vi(t)+c1(t)r1(gi best-
xi(t))+c2(t)r2(Pbest-xi(t))

(9)

x(t+1)=x(t)+V(t+1)

(10)

式中：Vi(t)表示第i個粒子在第t次進化中的速度；r1和r2是(0, 1)之間隨機數(shù)；gi best表示第i個粒子經(jīng)過第t次進化之后產(chǎn)生的局部最優(yōu)解；Pbest表示全局最優(yōu)解；x(t)表示第i個粒子在經(jīng)歷第t次進化時所在的位置；w(t)表示第t次進化中的慣性權(quán)重；c1(t)和c2(t)分別表示第t次進化中的學習因子以及社會因子。慣性權(quán)重、學習因子和社會因子在種群進化過程中會不斷變化以適應種群的變化。

SA算法思想最早是由Steinbrunn等[19]于1953年提出，稱為Metropolis準則。SA算法是一種隨機尋優(yōu)算法，基于蒙特卡洛迭代求解策略，其基本根據(jù)為對物理中固體物質(zhì)退火過程的模擬，以某個特定的概率接受優(yōu)化問題的求解結(jié)果。SA算法結(jié)合概率突跳特性，在限制區(qū)域空間內(nèi)采用隨機搜尋的方法得到目標函數(shù)的全局最優(yōu)解。由于一般需要解決的優(yōu)化問題與固體退火過程有很大的相似性，使得SA算法成為一種通用優(yōu)化算法，并得到廣泛應用。

PSO算法公式簡單和收斂快速，但是容易過早成熟而陷入局部收斂。而SA算法雖然全局尋優(yōu)能力不強，但是能夠在求解陷入局部收斂時概率性地跳出局部收斂，從而保證可以求得全局最優(yōu)解。綜合考慮上述兩種算法的優(yōu)缺點，筆者提出了改進型退火粒子群優(yōu)化算法(IPSO-SA)。

為了加快IPSO-SA算法中PSO部分的搜索速度，加入了自適應算子操作，其目的是提高收斂速度和收斂精度，減少陷入局部較優(yōu)的可能，包括自適應權(quán)重系數(shù)、社會因子和學習因子組成。權(quán)重系數(shù)ω、社會因子c1和學習因子c2分別為

(11)

(12)

(13)

式中：ωmax，ωmin分別為最大、最小慣性權(quán)重；(為加速因子；gi best表示第i個粒子在第t次進化后的局部最優(yōu)解；Pbest表示全局最優(yōu)解，c1 max，c1 min分別為最大、最小學習因子；c2 min是最小社會因子。

IPSO-SA算法的流程圖如圖2所示，具體過程為

1) 參數(shù)首次初始化：辨識粒子總數(shù)pop_size，模型中待辨識參數(shù)的個數(shù)以及其取值范圍，PSO中的學習因子c1和c2的初始值，粒子最大速度向量以及最小速度向量，IPSO-SA迭代的總次數(shù)max_p，內(nèi)嵌的種群代數(shù)max_s，降溫系數(shù)以及初始溫度T。

2) 對粒子種群進行初始化，單獨一個種群中有pop_size 個粒子，每個種群都有一組辨識參數(shù)。確定適應性函數(shù)。比較適應度向量，可以挑選出全局最優(yōu)粒子、全局最優(yōu)適應度值、局部最優(yōu)粒子種群和局部最優(yōu)適應度向量。

3) 進行PSO計算。通過更新權(quán)重因子、速度因子和社會因子，接著更新粒子的速度和位置向量，得到新的種群。

4) 將更新后的粒子重新代入模型中，得到新一輪的值，之后依據(jù)適應性函數(shù)計算得到種群新的適應度向量。

5) 進入SA操作。將新種群通過SA規(guī)則替換舊種群。

6) 重選全局最優(yōu)粒子，全局最優(yōu)適應性值，局部最優(yōu)粒子種群、局部最優(yōu)適應度值。

7) 將當代粒子群作為下一代的初始種群。更新溫度T、慣性權(quán)重、速度因子和社會因子后，算法跳至步驟3) 繼續(xù)迭代，直到gen_max達到預設值。

圖2 IPSO-SA算法的程序流程圖Fig.2 IPSO-SA algorithm program flow chart

4 結(jié)果與分析

4.1 參數(shù)辨識

分別在Matlab中編寫PSO，SA，IPSO-SA以及簡化TOMGRO模型的M文件代碼。采用2015年9月30日—2016年5月12日之間的室內(nèi)溫度、CO2濃度、光照強度和番茄生長數(shù)據(jù)進行參數(shù)辨識?？紤]到不同算法的收斂速度，將算法總代數(shù)初始化為300，保證所有算法都能完成求解。優(yōu)化算法在求解過程中調(diào)用簡化TOMGRO模型計算適應度函數(shù)值，辨識出模型中的不確定參數(shù)。參數(shù)辨識的結(jié)果如表2所示。

表2 參數(shù)辨識結(jié)果Table 2 Parameter identification result

根據(jù)IPSO-SA，PSO，SA每一代的均方根誤差，可分別得到3 種算法的收斂情況，如圖3所示。該圖顯示SA在第32 代時陷入局部收斂，直到第240 代時又跳出了局部收斂，表明SA具有良好的全局求解能力。PSO運行至第20 代時就已完成收斂。而IPSO-SA在第3 代時達到收斂，此后分別在第65 代時跳出局部收斂。表明IPSO-SA不僅繼承了PSO在前期快速收斂的特性，而且融合了SA全局搜索的優(yōu)點，具備了較強的跳出局部收斂的能力。

根據(jù)優(yōu)化算法辨識獲得的參數(shù)(表2)，結(jié)合溫室番茄產(chǎn)量預測模型獲得3 種優(yōu)化算法模擬番茄的產(chǎn)量，如圖4所示。在移植后第77 d左右和種植周期結(jié)束時，IPSO-SA算法跟隨實際番茄產(chǎn)量的能力要高于其他兩種算法，可知IPSO-SA算法的魯棒性更強，對于果實成熟的響應更快。PSO和SA最終收斂時的均方根誤差分別比IPSO-SA要大90.8%和55.5%，可見IPSO-SA模擬和實際產(chǎn)量的偏差要小于其他兩種算法。

圖3 3 種算法收斂情況Fig.3 Three kinds of algorithm convergence

圖4 辨識過程中3 種優(yōu)化算法產(chǎn)量預測值和實測值Fig.4 Experimental and forecasted yield with three algorithms in identification process

4.2 產(chǎn)量模型驗證

采用2016年9月30日—2017年5月12日期間的溫室環(huán)境數(shù)據(jù)和番茄生長數(shù)據(jù)對溫室番茄產(chǎn)量預測模型進行驗證。為排除偶然誤差，先利用2015年9月30日—2017年5月12日期間的相應數(shù)據(jù)進行10 次參數(shù)辨識，IPSO-SA，PSO，SA等3 種算法分別運行10 次得到10 組參數(shù)，即10 個模型。然后利用得到的模型和2016年9月30日—2017年5月12日期間的相應數(shù)據(jù)對該生長周期的番茄產(chǎn)量進行預測，對預測值取平均，結(jié)果如圖5所示。

圖5 驗證過程中3 種優(yōu)化算法產(chǎn)量預測值和實測值Fig.5 Experimental and forecasted yield with three algorithms in forecasting process

雖然3 種算法校準后的模型均能較好地預測番茄產(chǎn)量，但是根據(jù)196～224 d內(nèi)的產(chǎn)量曲線可知：IPSO-SA預測的產(chǎn)量更接近實際值。IPSO-SA預測的總產(chǎn)量與實際總產(chǎn)量的相對誤差為2.2%，比PSO和SA分別要小2.1%和0.7%，PSO和SA均方根誤差的平均值分別比IPSO-SA高66.1%和17.7%，表明IPSO-SA的參數(shù)辨識能力要優(yōu)于PSO和SA，預測的總產(chǎn)量精度高達97.8%。

5 結(jié) 論

筆者將PSO算法與SA算法結(jié)合并引入PSO算法的自適應算子操作加速算法的收斂，辨識得到產(chǎn)量模型能很好地預測溫室番茄產(chǎn)量。通過實際模擬分析，將IPSO-SA算法用于簡化TOMGRO模型參數(shù)辨識時，搜索最優(yōu)解的收斂速度和收斂精度均有較大提升。校準后的簡化TOMGRO模型可預測溫室番茄不同溫度條件下一個生產(chǎn)周期的總產(chǎn)量，對溫室番茄栽培具有一定指導意義。由于溫室環(huán)境以及植物生長過程的復雜性，為更精確地預測溫室番茄產(chǎn)量，需考慮晝夜溫度分開設定對番茄產(chǎn)量的影響，并結(jié)合室外氣候變化進行下一步的研究。