基于溫室內氣象數據的適用于冬季日光溫室番茄蒸騰量各子模型對比研究

王 淼，羅新蘭，尹佳琪，蘇 慧,常姝婷

(1.錦州市氣象局,遼寧 錦州 121000; 2.沈陽農業大學農學院，遼寧 沈陽 110866;3. 遼陽市氣象局，遼寧 遼陽，111000； 4. 葫蘆島市氣象局， 遼寧 葫蘆島，125000)

隨著水資源的匱乏及科學有效灌溉的需求，節水灌溉及作物需水量的準確模擬計算成為熱門課題。例如在溫室或大棚等果蔬生產環境下，對節約優質適量灌溉的需求更加迫切。

北方寒區日光溫室冬季僅在中午前后通風，其生產環境相對封閉，蒸騰作用消耗潛熱降低植株溫度，釋放水汽增加溫室內濕度，室內濕度的增加又會一定程度地抑制植株蒸騰作用，增加植株病蟲害發生頻率。Penman公式是世界糧農組織推薦的一種計算作物需水量的方法[1]，Monteith在1963年對此改進，提出以能量平衡和水汽擴散為基礎的Penman-Monteith方程(以下簡稱P-M方程)，適用于多種環境下作物蒸發蒸騰量的直接計算[2-3]。Moritlle等[4]試驗驗證了P-M方程可以用于溫室環境下作物蒸騰量的計算。然而，由于影響蒸騰的因素較多且很多參數很難精準測定，尤其是不同環境下邊界層空氣動力學阻力(ra)、氣孔平均阻力(rc)，目前還缺乏通用、精準、可用于不同類型溫室的子模型計算式。大多數學者認為ra與植株葉片和溫室內空氣之間對流熱交換形式密切相關，不同風速下應采用不同的對流熱交換模型[1]。另外，溫室內風速是影響ra的主要因素，大部分研究應用P-M方程計算作物蒸騰量時，測量微風速后直接帶入半經驗公式即求出ra值[5-19]。冬季室內空氣流通不均勻，微風速測量難度很大時，也可采用不含風速項的模擬式來計算ra值[20-22]。溫室內環境(室內總輻射Rg、飽和水汽壓差VPD、室內空氣溫度Ta)是影響子模型rc大小的重要因素，常采用3種影響因素的半經驗式或植株整個生長周期的特征值來模擬rc值[5,16-18，20]，鮮有利用P-M反演式來計算。國內外研究中常將ra、rc子模型直接應用于蒸騰量試驗中，而少有對各子模型對比分析擇優后帶入P-M方程對植株蒸騰速率模擬計算。

本文分別選出4組適用于北方寒區冬季日光溫室內的ra、rc子模型，試驗證明不同子模型對P-M方程計算植株蒸騰速率有明顯影響，通過初步篩選、檢驗各項指標擇優的方式選出基于此環境下模擬植株蒸騰量準確率最高的ra、rc子模型。北方寒區日光溫室內具有弱光、高濕、低溫及低風速(甚至為0 m·s-1)的環境特點，同時，本文還對多株植株同時進行試驗，這為應用P-M反演式計算rc值提供了條件，子模型ra、rc的選取對P-M方程的準確模擬及冬季日光溫室園藝作物灌溉需水量的確定具有重要參考意義。

1 材料與方法

1.1 試驗概況

試驗于2016年11月—2019年2月在沈陽農業大學北山科研基地日光溫室(123°5′E，41°8′N)內進行。溫室方位南偏西7°。東西長度60 m，跨度7.5 m，脊高3.5 m，北墻高2.3 m，墻體厚45 cm，北墻外側設有冷棚，試驗期間揭簾時間8∶30，蓋簾時間14∶30。

供試番茄品種為適于北方日光溫室長季節栽培的圣羅蘭。植株采用南北成行基質袋培，每行鋪有一條滴灌管，滴頭與植株相對應，滴頭出水口插于栽培袋中間位置，晴天的9∶00、12∶00、14∶00 分別灌溉100 ml·株-1·次-1,基質水分控制在田間持水量80%以上，地面覆蓋塑料薄膜(無裸露土壤)，行距1 m，株距20 cm，共定植920株，定植密度為3.38株·m-2。選擇長勢良好的3株番茄作為測試樣本。溫室頂部通風口可人工開啟，晴天中午室內溫度超過30℃時通風口打開10 cm進行通風，平均通風時間為15 min。其余時間處于密閉狀態，試驗期間植株處于成熟采摘期，莖粗1.41～1.48 cm、株高188.3～193.6 cm。

1.2 觀測項目

室外及室內中部安裝HOBO自動氣象站以獲取試驗所需基本氣象數據(輻射、溫度、相對濕度、土壤含水量等)；番茄群體中間安裝北京雨根公司的環境及莖流監測儀，測定項目包括：植株上方的太陽凈輻射、植株的莖流值、空氣溫濕度、葉溫(用與葉片相同顏色的綠色別針將熱電偶與葉片固定)；采用熱球風速儀測量植株上方微風速。試驗期間24 h連續測試，取樣間隔10 min，蒸騰速率估算中用到的相關氣象要素觀測項目及儀器位置如表1所示。基質含水量探針置于栽培袋中間位置。

表1 冬季日光溫室內氣象要素觀測項目及儀器位置Table 1 Meteorological elements tested items and instrument position in winter solar greenhouse

單株植株實測蒸騰速率[4,9]計算公式：

(1)

式中，Tf為單株植株蒸騰速率(mm·h-1)；La為植株監測探針上方的葉面積(m2)；LAI為單株葉面積指數(植株葉面積/植株占地面積)；Fs為質量形式莖稈液流速率(g·h-1)，根據莖流計流速公式計算。

1.3 分析方法

1.3.1 葉面積指數(LAI) 葉面積采用主干上主葉長寬系數法計算，測量主葉的葉長(包括葉柄長度)和最大葉寬,計算公式[23]：

(2)

式中，LA為葉面積(cm2)；a為葉片長度(cm)；b為葉片寬度(cm)。LAI=LA/單株植株占地面積。

1.3.2 消光系數(K) 根據比爾-朗伯特(Beer-Lambert)定律得到：

(3)

式中，LAI為單株葉面積指數；I，I0分別為單株尺度的地面和植株上方的總輻射(W·m-2)。計算每日9∶00—14∶00每10 min的消光系數，根據馮新靈等[24]研究建議剔除不合理的消光系數數據(即K>1的數據)，再以日為單位將合理數據均值后作為每日的消光系數。

1.3.3 采用Penman-Monteith方程模擬日光溫室基質袋培單株植株蒸騰量 在考慮土壤熱通量、并結合適用于本環境的各邊界層空氣動力學阻力及氣孔平均阻力計算式，Penman-Monteith方程計算溫室內單株植株蒸騰量的基本式可轉換[14,25-26]為:

(4)

聯合國糧食及農業組織(簡稱FAO)推薦土壤熱通量(G)白天采用0.5Rn，夜間采用0.1Rn[4]，本文根據北方寒區日光溫室揭、蓋簾時間對應的溫室內日照時數對G進行修正，白天采用0.35Rn，夜間采用0.13Rn，已在前文對取值做了敏感性檢驗[11],本文不再贅述。

1.3.4 邊界層空氣動力學阻力模擬 于冬季日光溫室內進行，僅在晴天中午室溫超過30℃時通風(通風口開度10 cm，平均通風15 min)，適用于此環境下ra的4種計算方法如下:

(1)根據Campbell[22]研究結果，由葉片溫度與周圍空氣溫度差值引起的自由對流對ra起決定作用時，計算公式：

(5)

式中，ra1為單一葉片空氣動力學阻力(s·m-1)；d為葉片特征長度(m)(本試驗取值0.15 m)；Tl為葉片溫度(℃)；Ta為周圍空氣溫度(℃)；ra為單株植株邊界層空氣動力學阻力(s·m-1)；LAI為單株葉面積指數。

(2)根據康紹忠等[14]和劉賢趙等[19]研究結果，邊界層空氣動力學阻力是不同作物風速廓線的函數，由湍流傳輸阻力與葉片層流邊界層阻力組成，可表示為：

ra=rg+rb

(6)

式中，ra為邊界層空氣動力學阻力(s·m-1)，rg為湍流傳輸阻力(s·m-1)，rb為葉片層流邊界層阻力(s·m-1)。

根據Thom等[22]提出的適用于溫室內計算湍流傳輸阻力的公式，當Uz=0時，公式簡化為：

式中，Z為風速測量高度(m)，Z0為地面粗糙度(m)，d為零平面位移長度(m)。

根據康紹忠等[14]研究，對于一般作物：

Z0=0.13hc，d=0.64hc

式中，hc為番茄冠層高度(m)，其余符號意義同前。

根據此公式可避免由于溫室內風速過低，出現湍流傳輸阻力無窮大的情況。

葉片層流邊界層阻力rb與葉片特征寬度、風速有關, 當溫室內風速不為零時，葉片層流邊界層阻力不能忽略，此時rb可表示為：

式中，c為經驗常數(0.01 m·s-10.5)，u為冠層上方風速(m·s-1)，w為葉片特征寬度(m)。

(3)根據Bailey等[1]關于溫室內作物葉片的研究表明，邊界層空氣動力學阻力(ra)與葉面和空氣之間的對流熱交換系數(hp,a)的關系式為：

(7)

不同對流情況下，hp,a的計算式為：

其中

式中，Kc為空氣的熱導率(0.0267 W·m-1·K-1)，Gr為格拉曉夫數，Pr為普朗特數，β為空氣的熱膨脹系數(0.00375 K-1)，g為重力加速度(9.8 m·s-2)，Tl為葉面溫度(℃)，Ta為室內空氣溫度(℃)，Re為雷諾數，v為空氣運動粘性系數(1.6×10-5m2·s-1)，a為空氣熱擴散率(2.29×10-5m2·s-1)，u為室內空氣流速(m·s-1)，l為作物葉片特征長度，番茄取0.15 m。

(4)Stanghellini[27]應用葉片與空氣之間溫度差值及溫室內空氣流速的經驗公式計算ra值，公式可總結為：

ra=1174L0.5(L|Tl-Ta|+207u2)-0.25

(8)

式中，L為葉片特征長度(m)，Tl為葉片溫度(℃)，Ta為室內空氣溫度(℃)，u為作物冠層上方微風速(m·s-1)。

1.3.5 氣孔平均阻力模擬 作物通過葉片氣孔與周圍環境進行顯熱交換和水分傳輸，交換和傳輸過程中受氣孔阻力控制，氣孔阻力又受諸多氣象因子影響，主要包括溫室內植株上方總輻射(Rg)、溫室內飽和水汽壓差(VPD)、溫室內空氣溫度(Ta)。本文對處于成熟采摘期的番茄植株假定其在大葉模式下進行氣孔平均阻力的模擬，根據Boulard等[28]的研究表明，作物氣孔平均阻力可由以下各統計方程來計算：

(1)根據Pollet等[15]的研究結果此方程可總結為：

×[1+0.016(Ta-16.4)2]

(9)

式中，Rg為溫室內植株上方總輻射(W·m-2)，Ta為溫室內空氣溫度(℃)，D為溫室內飽和水汽壓差(kPa)。

(2)北方寒區日光溫室冬季VPD變化范圍很小，根據Pollet等[20]及Boulard等[28]研究，表達式可總結為Rg、Ta的方程，即：

(10)

(3)Demrati等[29]認為冬季日光溫室處于密閉狀態下，可直接將室內總輻射帶入下式計算出作物冠層對水汽的阻力(即氣孔平均阻力)：

(11)

(4)由于影響溫室內氣孔平均阻力的主要因素包括溫室內總輻射、溫度、飽和水汽壓差等氣象因子，因此可以通過反演的方式計算出植株樣本a的rc值作為植株樣本b、c的rc值并帶入P-M方程計算另兩株植株的蒸騰速率模擬值，rc計算表達式如下：

(12)

式中，λE為植株a的實測蒸騰潛熱值(W·m-2)，ra為植株a的邊界層空氣動力學阻力(s·m-1)，其余物理量為植株a單株尺度對應的數值且意義同前。

2 結果與分析

2.1 不同子模型模擬日光溫室內各邊界層空氣動力學阻力(ra)與氣孔平均阻力(rc)對比分析

基于模型模擬蒸騰速率的準確程度與ra、rc的取值有直接關系[2]，目前對于不同類型溫室中ra、rc計算公式的選取仍沒有一致的結論。為了驗證本文提出的幾種計算ra、rc的方法對P-M模型模擬精確程度的影響，將12月20日及12月21日的試驗數據計算的各ra、rc值帶入P-M方程計算逐時蒸騰速率，以對比各子模型在不同時刻的模擬效果。以1.3.4、1.3.5節中的公式為基礎，計算單株植株蒸騰速率的相對變化率，白天各時刻蒸騰速率計算結果如圖1所示。

由圖1(a)可見不同的ra計算式對白天逐時蒸騰速率影響最大的時間是11∶00—13∶00，蒸騰速率平均變化量為8.01%，影響最小的時間為9∶00及15∶00左右，蒸騰速率平均變化量為0.76%，其余時間蒸騰速率平均變化量為5.11%；由圖1(b)可見不同的rc計算式對白天逐時蒸騰速率影響最大的時間是11∶00—13∶00，蒸騰速率平均變化量為42.21%，影響最小的時間段為9∶00—10∶00，蒸騰速率平均變化量為27.48%。

由于北方寒區冬季嚴寒，日光溫室于夜間覆蓋保溫被以保證室內溫度滿足植株正常生長需要，而此時蒸騰量很小。為了進一步明確子模型選取對植株蒸騰作用的顯著影響，本試驗對夜間番茄蒸騰速率作簡要計算說明，對植株夜間的蒸騰速率值比較分析，給出不同ra、rc組合形式對夜間植株蒸騰量的影響情況，結果如表2所示，不同的ra計算式對夜間蒸騰速率影響較大，且影響程度相近，平均變化量為46.13%；各rc計算式對夜間蒸騰速率影響程度較一致，平均變化量為30.25%。

表2 基于不同子模型的冬季日光溫室內夜間番茄各邊界層空氣動力學阻力與氣孔平均阻力模擬的蒸騰速率/(mm·h-1)Table 2 Transpiration rate of different aerodynamic resistance (ra)、average stomatal resistance (rc)on different submodel in solar greenhouse during night winter

總體上看，不同的ra、rc對P-M方程模擬值的影響程度均是夜間大于白天，白天rc對P-M方程模擬值的影響程度明顯大于ra，夜間則相反，因此在應用P-M方程計算植株蒸騰速率時子模型ra、rc的選用是非常重要的。

2.2 日光溫室內單株尺度子模型運行結果

選取試驗期間有典型天氣的2017-12-13—2017-12-25(其中13日與17日為霧霾天，14日與24日為陰天，其余時間為晴天)連續13 d的測試值。根據1.3.4與1.3.5節各公式計算邊界層空氣動力學阻力(ra)與氣孔平均阻力(rc)，各公式模擬值的對比見圖2。

由圖2(a)可以看出，基于1.3.4節此環境下公式Campbell(1977)與公式Stanghellini(1995)的變化趨勢相同，波動范圍不同，公式Campbell(1977)計算的ra平均日波動范圍在晴天為185～489 s·m-1，陰天為249～357 s·m-1；公式Stanghellini(1995)計算的ra平均日波動范圍在晴天為450～1 403 s·m-1，陰天為649～1 011 s·m-1。而公式劉賢趙(2001)、公式Bailey(1993)計算的ra值波動振幅過于平穩。

由圖2(b)可以看出，冬季日光溫室除午間通風時段，其余時間處于密閉狀態，室內濕度一直維持較高值，VPD變化范圍很小，考慮VPD時，rc的波動范圍反而不那么明顯，所以如公式Pollet(2000)所示僅考慮Rg與Ta時rc的計算值在全天均呈現明顯的波動性。鑒于此，1.3.5節中公式Pollet(2000)及反演式計算的rc平均日波動范圍在晴天分別為261～1 359、253～1 356 s·m-1；陰天分別為250～1 291、235～1 260 s·m-1。此環境下考慮影響rc的3個因素(植株上方總輻射,Rg;飽和水汽壓差,VPD;室內平均空氣溫度,Ta)及一個因素(Rg)時在白天時段rc模擬值基本相同，這與實際情況不符。

由于ra、rc變化具有波動性，由圖2分析結果初步確定1.3.4節中由空氣溫度、葉溫及植株生理指標(葉片特征長度、LAI)代入公式Campbell(1977)、公式Stanghellini(1995)計算的ra子模型、1.3.5節中根據室內植株上方總輻射(Rg)、空氣溫度(Ta)代入公式Pollet(2000)及P-M方程反演式計算出的rc模擬值符合其隨時間波動的特征。

2.3 不同子模型模擬日光溫室內邊界層空氣動力學阻力(ra)與氣孔平均阻力(rc)的驗證與篩選

為了進一步說明子模型對P-M方程模擬植株蒸騰速率的影響程度以及模型的通用性，分別將ra于12月18日—1月1日及rc于12月11日—12月25日的子模型計算值代入P-M方程并求得蒸騰速率模擬值與實測值的Pearson相關系數、平均相對誤差(MRE)、平均絕對誤差(MAE)、均方根誤差(RMSE)及模型效率指數(EF)[30-31]，結果見表3。由表3分析得出，北方寒區冬季日光溫室生產環境中各ra子模型的模擬值與實測值相關性最強的是Campbell與Bailey模擬式，其與實測值對比的相關系數均為0.826(劉賢趙0.813；Stanghellini 0.822)；平均絕對誤差(MAE)Campbell模型為0.02，其余3個子模型均為0.03；Campbell模型的均方根誤差(RMSE)最小；模型效率指數(EF)Campbell>Bailey=Stanghellini>劉賢趙。因為本試驗測試時間為冬季，溫室中午前后開窗通風時室內風速垂直流通且并不均勻，其余時間處于密閉狀態，微風速的測量偏差較大。通過各檢驗系數結合子模型整體波動性，綜合分析得出選用包含空氣溫度、植株葉片溫度及植株生理指標(葉片特征長度、單株植株葉面積指數LAI)等參數模擬單株植株邊界層空氣動力學阻力值；各rc子模型計算得出的模擬值中反演式與實測值對比的相關系數為0.955(Pollet 0.861；Boulard 0.863，Demrati 0.860)，反演式的平均相對誤差(MRE)、平均絕對誤差(MAE)、均方根誤差(RMSE)均低于其余三者，反演式的模型效率指數(EF)為0.95，其準確率高于其余三者，所以本文選用反演式模擬單株植株氣孔平均阻力值。

表3 各子模型模擬的植株蒸騰速率模擬值與實測值的分析檢驗Table 3 Inspection analysis between simulated on different submodel and measured of transpiration rate

3 討論與結論

3.1 討論

基于模型模擬蒸騰速率的準確程度與邊界層空氣動力學阻力、氣孔平均阻力的取值有直接關系[2]，而二者數值的精確獲得很難。冠層內外風速及溫度的變化會產生不同的氣流模式類型[1],氣孔平均阻力的取值不僅與太陽輻射、VPD、空氣溫度、作物水分等有關[32],且還受冠層結構、冠層上方氣流及土壤條件等影響[12]。已有部分研究于空氣流通較均勻的塑料大棚[5]、玻璃溫室[6]、防蟲遮陽網溫室[7]、日光溫室[8-13,24]中直接采用運用微風速計算的ra值、半經驗公式總結的rc值或總結出的特征值模擬整個試驗周期的子模型運行結果，本試驗驗證的4組ra、rc計算式對P-M方程模擬植株蒸騰速率有明顯影響，驗證各組子模型準確率后再運用其計算植株蒸騰量值是非常必要的。

Bailey等[1]在英國(5月)及西班牙(7月)空氣流通較均勻的溫室內利用葉片邊界層空氣流通模式模擬ra值。結果表明：開窗通風時，強迫對流模擬式更接近實測值，兩地區的模擬值與實測值平均差值分別為2.8%、3.2%；康紹忠等[14]和劉賢趙等[19]在Thom等[22]提出的溫室內湍流傳輸阻力計算式基礎上提出葉片層流邊界層阻力值，提高了冬季大棚作物需水量計算的準確率；Stanghellini[27]利用葉片與空氣溫差及室內空氣微風速的半經驗公式也得到了較好的模擬效果；Campbell等[25]總結出基于葉溫、空氣溫度、植株生理指標(LAI、葉片特征長度)為基礎的模擬式，避免了由于冬季溫室內空氣流通不均勻時誤差較大的問題，經檢驗在冬季日光溫室空氣流通不均勻時采用此模擬式可提高試驗準確率。經對比分析，本文選用Campbell等[25]總結出的ra模擬式，并將LAI換算為單株尺度，當室內微風速測量困難時可采用不含有風速項的ra模擬式。

國內外相關研究中多是利用溫室內影響氣孔平均阻力的氣象因子(Rg、VPD、Ta)總結的半經驗公式來模擬rc值[20,28-29]，本文利用反演式求得rc值，經驗證利用反演式模擬植株蒸騰速率的誤差較小，且Pearson相關系數達0.955，準確率明顯提高。

本試驗在北方寒區很少通風的日光溫室內進行，此條件下空氣流速不均勻，微風速測量難度較大。本文僅對冬季日光溫室內ra、rc篩選擇優，并未在夏季日光溫室內驗證。在今后試驗中將增加夏季空氣流通較均勻時各子模型的篩選驗證。試驗番茄為基質袋培，地面覆蓋塑料薄膜，無裸露土壤。番茄根部水熱狀況與盆栽或土壤栽培有區別，各組子模型的適用性對蒸騰結果的影響有待進一步研究及優化。

總體來看，本文主體采用應用較廣泛的P-M方程模擬溫室單株植株蒸騰量，對邊界層空氣動力學阻力及氣孔平均阻力不同子模型模擬式作對比分析；本試驗選取的ra、rc子模型模擬植株蒸騰速率的準確率較高，特別是運用反演式計算的rc值達到了較優的模擬效果。試驗對比分析結果可為北方寒區冬季日光溫室空氣流通并不均勻(微風速測量難度較大)、選取多株植株同時進行試驗時子模型的選取提供參考。

3.2 結 論

(1)總體來看，ra、rc隨時間波動變化，利用不同ra、rc模擬式計算北方寒區日光溫室蒸騰速率的相對變化率有差別。不同子模型模擬的ra在白天、夜間平均變化量分別為5.73%、46.13%；不同子模型模擬的rc在白天、夜間平均變化量分別為37.10%、30.25%，rc模擬式對白天蒸騰速率的影響程度明顯大于ra的選擇，雖然夜間不同的ra計算式對蒸騰速率模擬值影響程度大于rc約15.88%，但由于北方寒區日光溫室冬季番茄蒸騰作用主要集中于白天，因此模擬植株蒸騰速率時子模型rc的選擇尤為關鍵。

(2)本試驗條件下，選取的子模型計算得出的晴天、陰天番茄植株邊界層空氣動力學阻力(ra)的變化范圍分別為185～489、249～357 s·m-1，氣孔平均阻力rc變化范圍分別為253～1 356、235～1 260 s·m-1，應用P-M方程模擬日光溫室作物蒸騰量時應根據溫室內特定的氣象環境選擇合適的ra、rc模擬式。

(3)本研究驗證了ra、rc(各4組)子模型組合式分別代入P-M方程模擬植株蒸騰速率并與實測值對比分析，結果表明：Campbell模型的Pearson相關系數、MRE、MAE、RMSE、EF指標均優于其余3種模型，又結合子模型的波動特征，綜合考慮本試驗的ra值選取Campbell模型模擬；反演式的各指標均優于其他模擬式且其Pearson相關系數為0.955、MRE為7.16%、MAE為0.01、RMSE為8.56、EF為0.95，因此，本文選取反演式計算rc值。建議在溫室內微風速測量難度較大時，可采用測試葉溫、空氣溫度、LAI及葉片特征長度的方式來計算ra值；由于rc值不受植株自身條件影響，僅與室內環境要素有關，因此當試驗同時選取多株植株時，可以通過其中某一植株的反演式計算出的rc值同樣作為其余植株的rc值代入P-M方程計算相應單株植株的蒸騰量可達到較優的模擬效果。