供水不足地區園藝植物類溫室自動灌溉系統的設計與試驗

國俊保， 高愛玲

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

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

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

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

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

分交換過程的各個環節，并且能夠仿真出溫室內蒸騰速率與微氣候的變化過程。

1 材料與方法

1.1 試驗環境與方法

試驗于2017年6月在山東省濰坊地區10 m×10 m的Venlo型玻璃溫室進行。溫室由4 mm厚的園景玻璃構成，在屋脊兩側裝備了連續的屋頂通風口，通風口在室外溫度超過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，花盆內填充均勻分布的泥炭，密度為0.12 g/cm3。

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

1.2 溫室微氣候的數據采集

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

1.3 基質栽培的參數

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

1.4 植物活動的數據

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

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

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

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

2.1 模型的基礎方程

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

(1)

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

(2)

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

2.2 輻射子模型

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

(3)

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

(4)

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

(5)

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

2.3 植物子模型

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

(6)

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

(7)

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

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

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

Rgabs-Trd-Qs=0。

(8)

可根據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是輻射的消光系數；供水充足的情況下該值為0.95，供水不足的情況下該值為0.64，這2個值是使用式(9)計算冠層上面與下面的短波輻射而來。這個差值的原因是在植物缺水的情況下，植物的葉片會自動地旋轉，這個生理機制能夠降低植物吸收的太陽輻射量，因此降低了蒸騰作用。

冠層的可感知熱量定義為

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

(11)

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

一個給定區域的潛熱密度定義為

(12)

式中：γ是溫濕常數；VPDa是空氣的飽和水汽壓差，Pa；Δ是飽和水蒸氣壓曲線的斜率；根據文獻[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 水平衡模型 根據文獻[13]模型，可通過下式計算基質水勢的減少量ψw：

(15)

式中：θres與θsat分別是土壤剩余含水量與土壤飽和含水量；α(kPa)與N是2個待校準的參數，這2個參數來自于拉力計與含水量傳感器的測量數據。根據多線性擬合的結果，可獲得以下參數值：θ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數量級。

3 模型驗證

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

3.1 試驗結果與分析

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

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

3.2 空氣溫度與濕度

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

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

3.3 泥炭土的基質勢

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

表1 溫室內微氣候各個參數的測量結果

3.4 觀賞植物的相關參數

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

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

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

4 結論

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