供水不足地區(qū)園藝植物類溫室自動灌溉系統(tǒng)的設計與試驗

國俊保， 高愛玲

(濰坊科技學院，山東壽光 62700)

隨著全球氣候變暖的加劇，水資源日益減少，農(nóng)業(yè)生產(chǎn)過程中須要消耗大量的水資源，因此，提高農(nóng)業(yè)生產(chǎn)過程中水資源利用率成為了當前農(nóng)業(yè)領域的研究重點[1]。在我國西北地區(qū)，大量區(qū)域目前仍處于水資源匱乏狀態(tài)，所以適用于水資源緊缺情況的溫室灌溉系統(tǒng)是一個有實際意義的研究內(nèi)容[2]。并且我國西南部區(qū)域時而發(fā)生旱災，如果能夠設計出水資源緊缺情況的溫室灌溉系統(tǒng)，則能夠幫助旱災地區(qū)減少經(jīng)濟損失[3]。

許多研究人員采用人工智能技術與機器學習技術對溫室的歷史數(shù)據(jù)進行學習，使用訓練獲得的模型預測溫室內(nèi)植物未來的變化趨勢[4-5]。在溫室供水充足的情況下，溫室內(nèi)微氣候的變化趨勢較為穩(wěn)定，所以人工智能與機器學習訓練的模型能夠準確地預測溫室內(nèi)植物參數(shù)的變化趨勢。但是如果發(fā)生突發(fā)性旱災，采用人工智能或者機器學習技術則難以獲得理想的預測與控制效果。

計算流體力學能夠有效地模擬各種流場的過程[6]，溫室中土壤-植物-空氣的水分交換過程是一種典型的流場過程[7]，并且計算流體力學不須要通過對歷史數(shù)據(jù)的學習即可建立模型，因此計算流體力學建立的預測模型實時性較高，能夠隨著問題的實時變化而自適應地改變模型的參數(shù)，從而能夠應對一些突發(fā)性狀況的發(fā)生。

為了在不影響植物蒸騰速率的前提下降低灌溉的水量，采用計算流體力學對溫室微氣候、 植物活動參數(shù)與土壤基質參數(shù)均進行建模與預測， 該模型考慮了土壤-植物-空氣水

分交換過程的各個環(huán)節(jié)，并且能夠仿真出溫室內(nèi)蒸騰速率與微氣候的變化過程。

1 材料與方法

1.1 試驗環(huán)境與方法

試驗于2017年6月在山東省濰坊地區(qū)10 m×10 m的Venlo型玻璃溫室進行。溫室由4 mm厚的園景玻璃構成，在屋脊兩側裝備了連續(xù)的屋頂通風口，通風口在室外溫度超過20 ℃的時候打開，一般通風口在白天打開，在夜晚關閉。為了防止溫室溫度過高，采用1個遮陽幕，遮陽幕的成分為鋁(16%)與聚乙烯(84%)，導熱率為0.15 W/mK，透射率為50%。溫室地面為混凝土，地面放置4個架子(3 m×1.5 m)，架子高度為0.8 m，水仙花花盆均勻地放置于4個架子上，花盆容積為0.74 L，高度為8.7 cm，花盆內(nèi)填充均勻分布的泥炭，密度為0.12 g/cm3。

架子上植物的冠層面積約為18 m2，葉面積指數(shù)(leaf area index，LAI)為2.36，植物密度為15盆/m2，植物平均高度為24 cm。每天使用滴灌設備為每盆植物澆灌2次完全營養(yǎng)液，其中3個架子的植物進行供水充足的試驗，這些花盆中泥炭的水勢保持在-2 kPa以上；第4個架子的植物進行供水不足的試驗，周期性地停止灌溉，使其處于供水不足的狀態(tài)。

1.2 溫室微氣候的數(shù)據(jù)采集

使用傳感器采集溫室內(nèi)的微氣候數(shù)據(jù)。試驗溫室的總體結構如圖1所示，使用溫度傳感器(T)與濕度傳感器(R)測量空氣的溫度與相對濕度，傳感器分別放于3個位置：植物上方15 cm處(T1，RH1)，植物內(nèi)部(T2，RH2)，架子下方(T3，RH3)。此外測量葉片的溫度(TI)，在地面放置PT100鉑探針，測量架子下方的地面溫度。使用CNR4型凈全輻射傳感器采集輻射的長波、短波數(shù)據(jù)，在1.29 m的高度放置1個熱球式電風速計，測量冠層的氣流速度。每隔3 s采集1次上述所有參數(shù)數(shù)據(jù)，每隔10 min將數(shù)據(jù)保存至數(shù)據(jù)庫中。

1.3 基質栽培的參數(shù)

在6個花盆中埋了負壓計與含水量傳感器，分別測量花盆的泥炭基質勢ψ(kPa)與含水量θ(體積分數(shù))，其中3個花盆為供水充足的花盆，另外3個為供水不足的花盆。每隔3 s采集1次上述所有參數(shù)數(shù)據(jù)，每隔10 min將數(shù)據(jù)保存至數(shù)據(jù)庫中。

1.4 植物活動的數(shù)據(jù)

通過有孔石膏在6張植物葉片下面固定Cu-Cs熱電偶，每隔3 s測量1次6張葉片的溫度，將10 min總數(shù)據(jù)的平均值作為結果保存至數(shù)據(jù)庫中。使用氣孔計測量葉片的氣孔阻力。在冠層的中間放置1個高分辨率Melter-Toledo電子秤，測量冠層蒸騰所流失的水量，電子秤的托盤(56 cm×56 cm)保持水平，每隔1 h記錄花盆流水量。

2 建立計算流體力學(CFD)模型

為了建立計算流體力學模型，在感興趣區(qū)域周圍定義1個計算域，然后將計算域離散化為表征體元，最終求解計算域中質量(空氣、水分)與能量的守恒方程量。根據(jù)文獻分析[8]，溫室的湍流一般位于溫室的高處，為了記錄湍流的運動變化，將變量分為1個平均分量與1個波動分量，并且將等式平均化處理，獲得RANS等式(雷諾平均等式)。

使用ANSYS Fluent 15.0流體計算軟件仿真CFD模型，該軟件使用有限體積法求解二維非穩(wěn)態(tài)對流擴散方程。

2.1 模型的基礎方程

因為溫室中連續(xù)通風口的主風向一般垂直于通風口，則氣流的二維平面垂直于屋脊，因此模型僅須要考慮二維的情況。二維湍流動能輸運方程的一般形式為

(1)

式中：Φ表示非量綱輸運量的集合，即動量、質量與能量；U與V是速度向量的元素；Γ是擴散系數(shù)；SΦ是輸運方程的源項(source term)。將k-ε湍流模型設為多組分湍流封閉模型，因為空氣密度(ρ)不僅依賴溫度還依賴空氣的含水量，所以無法應用Boussinesq模型。將空氣密度定義為一個關于溫度與質量的函數(shù)：

(2)

式中：T是局部溫度，K；R是理想氣體常數(shù)，R=8.31 J/(mol·K)；Pop為工作壓力，Pa；ωi是組分i的質量值；Mi為組分i的分子量，kg/mol。

2.2 輻射子模型

為了區(qū)分輻射中短波(0.1～3 μm)與長波(3～100 μm)的貢獻度，設計了一個輻射子模型。對于有限數(shù)量的離散立體角，可用式(3)求解單色亮度方程，每個離散立體角對應的向量方向為

(3)

聚集整個空間的每個單色亮度，即可獲得全局的輻射通量qs，表示為

(4)

熱平衡中每個單位體積的輻射貢獻度是能量守恒方程的一個源項(source term)，表示了輻射qs(s)的散度。

(5)

在當前的研究中，短波長與長波長之間存在差異[9]。與短波輻射相比，冠層吸收長波輻射的量可忽略不計，在白天冠層吸收的長波輻射小于短波輻射的4%。

2.3 植物子模型

2.3.1 多孔介質 將植物抽象為均勻分布的多孔介質，該多孔介質對氣流具有阻力(對應植物的氣孔阻力)，表示為能量守恒方程中的一個源項SΦ。根據(jù)Darcy-Forchheimer方程，植物每個單位體積的拉力可表示為黏性丟失項與內(nèi)部丟失項2個項的總和，如下式所示。

(6)

式中：K是多孔介質的滲透性，%；CF是非線性動量的衰減系數(shù)；ρa是空氣密度，kg/m3；μ是空氣的動態(tài)黏度，(N·s)/m2；V是空氣的速度值，m/s。根據(jù)文獻[10]分析，實際的空氣動態(tài)黏度非常低，是10-5(N·s)/m2數(shù)量級，所以可以忽略式(6)的第1項，可得：

(7)

式中：LAD是葉密度，定義為LAI/冠層高度(H)；CD是拉力系數(shù)。

從式(7)可以看出，僅須要獲得CD值即可求解植物每個單位體積的拉力，鳳仙花的CD值約為0.32。

2.3.2 能量平衡模型 冠層的太陽輻射為Rg，植物通過蒸騰作用以及與周圍空氣的熱交換進行能量傳遞，假設蒸騰作用的熱通量密度為Trd，周圍空氣的熱通量密度為Qs。最終，葉片的能量平衡公式寫為

Rgabs-Trd-Qs=0。

(8)

可根據(jù)Beer法則[式(9)]直接推導出冠層每個細胞吸收的輻射量Rgabs，W/m2：

Rg(y)=Rg0exp[-KcLAD(H-y)]；

(9)

Rgabs=Rg(yi)-Rg(yi+1)。

(10)

式中：yi與yi+1分別是細胞i的上一坐標與下一坐標。Rg0是冠層上的全局輻射，W/m2；Kc是輻射的消光系數(shù)；供水充足的情況下該值為0.95，供水不足的情況下該值為0.64，這2個值是使用式(9)計算冠層上面與下面的短波輻射而來。這個差值的原因是在植物缺水的情況下，植物的葉片會自動地旋轉，這個生理機制能夠降低植物吸收的太陽輻射量，因此降低了蒸騰作用。

冠層的可感知熱量定義為

Qs=2LADρaCp[(TI-Ta)/Ra]。

(11)

式中：ρa是空氣密度，kg/m3；Cp是指定壓強對應的空氣熱量，J/(kg·K)；TI與Ta分別是葉片與空氣的溫度。

一個給定區(qū)域的潛熱密度定義為

(12)

式中：γ是溫濕常數(shù)；VPDa是空氣的飽和水汽壓差，Pa；Δ是飽和水蒸氣壓曲線的斜率；根據(jù)文獻[11]，將氣動阻力Ra設為一個常量，Ra=271 s/m。Rs是氣孔阻力，s/m，采用FFD(全因子試驗設計)[12]方法可推導出Rs，表示為

Rs=(-115×rg-139×rh-39×τ+139×rg×rh+43×rg×τ+11×rh×τ+661×rg2-368×rg3)×[1+(ψ-11.41)1.05]。

(13)

式中：rg=(Rg-75)/75，rh=(RH-65)/10，τ=(T-20.5)/5.5，rg是輻射量，rh是濕度值，τ是溫度因子，ψw是土壤基質勢(kPa)的減少量，RH為相對濕度。將式(8)進行轉換，可計算出細胞的葉片溫度(TI)：

(14)

2.3.3 水平衡模型 根據(jù)文獻[13]模型，可通過下式計算基質水勢的減少量ψw：

(15)

式中：θres與θsat分別是土壤剩余含水量與土壤飽和含水量；α(kPa)與N是2個待校準的參數(shù)，這2個參數(shù)來自于拉力計與含水量傳感器的測量數(shù)據(jù)。根據(jù)多線性擬合的結果，可獲得以下參數(shù)值：θsat=0.887 cm3/cm3，θres=0.1 cm3/cm3，α=0.134 kPa-1，N=1.469。通過式(16)計算每個時間點i的含水量，式(16)定義了蒸騰作用與初始化含水量(θ0)的關系。

(16)

式中：θi是當前時間點的含水量i，cm3/cm3；θ0是初始含水量，cm3/cm3；ts是時段長度，本研究ts等于3 600 s；ss是架子的總面積；Travg，j是在時間點j時的冠層平均蒸騰速率；λ是水的潛熱量，J/kg；m是植物總量(m=264 kg)；ρw是水密度，kg/m3；Vp是花盆的容積，m3。

2.3.4 模擬程序 對于動量與熱量的輸運方程采用二階迎風差分格式[14]處理，提高輸運方程的計算準確率。對于壓力連接等式采用半隱式格式[15]來求解壓力-動量方程。所有變量的收斂條件是10-6數(shù)量級。

3 模型驗證

采用2個準確率性能指標評估本模型的預測準確率，分別為RMSE(均方根)與R2(確定系數(shù))，RMSE也稱為預測系統(tǒng)的擬合標準差，值越小表示預測越精確；R2是通過數(shù)據(jù)的變化來表征一個擬合的好壞，R2取值范圍為[0，1]，越接近1，擬合效果越好。

3.1 試驗結果與分析

首先對供水充足與供水不足2種情況分別進行分析與研究。將二維非穩(wěn)態(tài)仿真的時間步長設為2 h，周期為從6月14日01：00—23：00，從11：00—20：00打開通風口，6月14日對供水不足的花盆停止灌溉。試驗中測量的數(shù)據(jù)作為預測模型的輸入數(shù)據(jù)，測量的結果如圖2所示。從圖2可以看出，太陽輻射的變化曲線與氣溫、地面溫度均具有相關性，在13：00達到最大值。

在6月14日01：00使用6月14日00：00的數(shù)據(jù)進行一次初步仿真，對于供水不足的情況，初步仿真結果中的初始水量為0.652 cm3/cm3，使用基于水平衡的子模型計算每個時間點的土壤基質水勢。對于供水充足的情況，將土壤基質水勢設為恒定值，該值固定為-1 kPa。

3.2 空氣溫度與濕度

2種供水條件下3個測量位置的氣溫隨著時間的變化曲線如圖3與圖4所示。氣溫與葉片溫度的測量值與預測值表現(xiàn)出相同的趨勢，2種供水條件下氣溫的預測值與測量值能夠較好地擬合。由表1可以看出，對于供水充足情況，所有的結果中R2結果均大于或等于0.97，RMSE均小于或等于1.14；對于供水不足情況，所有的結果中R2結果均大于或等于0.94，RMSE均小于或等于 1.62。因為2種供水條件的溫度差異較小，所以供水是否充足對觀察周圍溫度的影響較小。2種供水條件的差異主要由初始的環(huán)境參數(shù)引起，供水不足情況的葉片溫度比供水充足情況約高1K，其原因是供水不足情況的蒸騰作用減弱，而蒸騰作用能夠降低葉片的溫度。2種供水條件的葉片底部溫度幾乎相同，因為葉片底部的蒸騰作用均較弱。

此外，研究了2種供水條件下3個測量位置的相對濕度，分別表示為RH1、RH2與RH3，相對濕度隨著時間的變化曲線如圖5所示。對于供水充足的情況，結果顯示冠層內(nèi)部相對濕度的預測值與測量值高于冠層的上面與下面。對于供水不足的情況，3個位置相對濕度的測量值與預測值較為接近?？傮w而言，2種供水條件的預測值與測量值實現(xiàn)了較好的擬合。從表1可以看出，對于供水充足的情況，R2大于0.88，RMSE小于8.66；對于供水不足的情況，R2大于0.78，RMSE小于6.69。

3.3 泥炭土的基質勢

在停止灌溉5 d后，統(tǒng)計了1 d內(nèi)測量與預測的泥炭土的基質勢值(圖6)。因為白天植物的視覺特征較為明顯，因此試驗主要關注白天。因為停止了灌溉，土壤內(nèi)的水量隨著植物的蒸騰作用而減少。從圖6可以看出，泥炭土基質勢的預測值較好地擬合了測量值，并且獲得了較高的預測準確率，預測性能指標分別為r2=0.99，RMSE=-3.38。

表1 溫室內(nèi)微氣候各個參數(shù)的測量結果

3.4 觀賞植物的相關參數(shù)

3.4.1 植物的氣孔阻力 為了分析光照對植物氣孔阻力的影響，對陰影區(qū)域葉片與光照區(qū)域葉片均進行了試驗分析，供水充足與供水不足2種情況的氣孔阻力變化曲線如圖7所示。從圖7-b可以看出，從大約10：00開始，氣孔阻力開始上升，因為此時土壤基質勢開始降低，植物通過增加氣孔阻力來限制植物的蒸騰作用。從圖7-a可以看出，供水充足情況的氣孔阻力也從10:00開始上升，此時陽光輻射與空氣氣壓差均開始上升，所以植物通過增加氣孔阻力來限制植物的蒸騰作用，但其上升幅度遠小于供水不足的情況。

CFD模型能夠有效地預測出陰影區(qū)域葉片與光照區(qū)域葉片的氣孔阻力變化曲線，供水充足情況的預測性能指標為r2>0.4，RMSE<168；供水不足情況的預測性能指標為r2>0.51，RMSE<527。

3.4.2 植物的蒸騰速率 供水充足與供水不足2種情況的植物蒸騰速率變化曲線如圖8所示。供水不足情況的蒸騰速率小于供水充足情況。供水不足情況的蒸騰速率最小測量值為1.43 g/(m2·h)，低于供水充足情況的最小蒸騰速率?？傮w而言，本模型能夠有效地預測供水充足與供水不足2種情況的植物蒸騰速率，供水充足情況的預測性能指標為r2=0.98，RMSE=5.87；供水不足情況的預測性能指標為r2=0.97，RMSE=5.64。此外，計算了24 h蒸騰速率的總預測誤差，供水充足情況的總預測誤差為3.32%，供水不足情況的總預測誤差為4.82%。

4 結論

目前溫室自動灌溉系統(tǒng)大多針對供水充足的穩(wěn)定環(huán)境，為了提高供水不足情況的灌溉效率，設計了一種采用物聯(lián)網(wǎng)的園藝植物類溫室自動灌溉系統(tǒng)。為了在不影響植物蒸騰速率的前提下降低灌溉的水量，采用計算流體力學對溫室微氣候、植物活動參數(shù)與土壤環(huán)境參數(shù)均進行建模與預測，該模型考慮了土壤-植物-空氣水分交換過程的各個環(huán)節(jié)，并且能夠仿真出溫室內(nèi)蒸騰速率與微氣候的變化過程。在山東濰坊地區(qū)采用水仙花盆景進行了試驗，對供水充足與供水不足2種情況均進行了試驗，該模型對溫室微氣候、植物活動參數(shù)與泥炭土基質勢均獲得了準確的預測效果。