李 群 汪 星 汪有科，3 納文娟 唐 燕 吳久江

(1.西北農(nóng)林科技大學水利與建筑工程學院， 陜西楊凌 712100； 2.寧夏大學農(nóng)學院， 銀川 752201；3.西北農(nóng)林科技大學水土保持研究所， 陜西楊凌 712100)

0 引言

大量研究表明Penman-Monteith(PM)模型可以應用在溫室環(huán)境蒸散估算[1-2]，并且目前是溫室蒸散模擬應用最廣泛的機理模型[3-4]。當前國內(nèi)外學者對溫室PM蒸散模擬的研究多以智能自控溫室為載體，簡易大棚應用PM模型估算蒸散量少有報道。我國90%以上的溫室屬于簡易大棚[5-6]，如小跨度塑料大棚的蒸散模擬鮮有報道[7-9]。由于溫室環(huán)境與大田環(huán)境相比風速明顯更低，PM模型的關鍵參數(shù)空氣動力學阻力ra使用常用的Perrier公式往往會被過高估計，進而導致PM模型的模擬效果較差[10-13]。為解決這一問題，BOULARD等[14]提出一種基于溫室外部環(huán)境估算溫室蒸散的方法，李軍等[7]在我國長江下游單棟塑料大棚中使用此方法模擬塑料大棚小白菜蒸散結果較好。但這種方法要求溫室空氣與室外環(huán)境耦合緊密，即需長期保持溫室自然通風，不適合在我國北方地區(qū)冬季使用。陳新明等[15]與王健等[12]在日光溫室環(huán)境中將風速取0并使用THOM等[16]的公式計算ra，這種方法有效解決了溫室低風速條件下ra被過高估計的問題，但對于大棚有通風情況下的蒸散估算仍未涉及。孫懷衛(wèi)等[8]建立了華東地區(qū)塑料大棚內(nèi)部風速與外部風速的線性關系，并基于能量平衡和紊流擴散理論建立了韭菜蒸騰模型。這種方法較為繁瑣，需要使用計算流體動力學(Computational fluid dynamics，CFD)分析棚內(nèi)氣流場特征。此外，上述研究使用的PM模型其約束函數(shù)形式均是唯一的，這會導致模型模擬精度缺乏比較，也沒有考慮模型的不確定性。

PM模型結構較為復雜，子模型涉及的待率定參數(shù)較多，模型的不確定性可能較強，基于上述考慮選擇使用Bayesian方法作為模型參數(shù)率定方法。Bayesian方法將參數(shù)的概率分布、實測數(shù)據(jù)、參數(shù)波動范圍和不確定性的相關假設相結合，既考慮了模型參數(shù)、輸入變量與輸出結果的不確定性，也考慮到了已知信息以及參數(shù)之間的互相影響[17]，是一種適合高維模型參數(shù)率定的優(yōu)秀方法[18-19]。近20年來，Bayesian方法在大田環(huán)境的蒸騰、蒸散模擬及模型不確定性分析中都得到了較好的應用[18-21]，但在溫室蒸散模擬研究中很少使用Bayesian方法，對溫室蒸散模型不確定性的分析更是鮮有報道。

本研究在關中地區(qū)簡易大棚中使用溫室中應用最為廣泛的PM蒸散機理模型，進行簡易大棚日尺度下章姬草莓蒸散模擬研究，以期為簡易大棚蒸散模擬提供參考。

1 材料與方法

1.1 研究區(qū)概況

試驗區(qū)位于陜西省渭南市白水縣雷牙鎮(zhèn)劉家卓村(35°12′N，109°36′E，海拔809 m)，年平均太陽總輻射5 363.52 MJ/m2，年平均氣溫11.4℃，年平均降水量577.8 mm。試驗地0～0.3 m土層土壤質地均一，土壤滲透性良好，砂礫、粉粒和黏粒3種粒級質量分數(shù)分別為70.1%、21.8%、8.1%，屬于沙質土。0～0.3 m土層土壤容重為1.39 g/cm3，環(huán)刀法測得田間持水率為35%，計劃濕潤層深度取章姬草莓根系平均深度0.3 m。定植時間為2018年8月31日傍晚，定植期共12 d。使用寬窄行起壟方式栽植，9壟，下寬0.45 m，上寬0.3 m，壟高0.35 m，間距0.25 m，每壟660株，栽植2行，行距0.3 m，株距0.2 m，每公頃約栽90 000株。栽植時同一行植株的花序朝同一方向，使草莓苗弓背朝花序預定生長方向，苗心露出畦面，根系平展埋入疏松土層。草莓苗栽植后及時灌溉定植水，灌至0.3 m深度土層使土壤含水率達到田間持水率，定植水每2 d灌溉一次，持續(xù)13 d。使用國產(chǎn)迷宮式滴灌帶，滴灌帶內(nèi)徑16 mm，額定流量2 L/h，滴孔間距200 mm，濕潤比取0.35，灌溉水利用效率取0.9。使用黑色PVE地膜以保溫保墑，定植后全地膜覆蓋栽培。

1.2 數(shù)據(jù)采集與處理

試驗區(qū)簡易大棚中央距離地面1 m處設有小型氣象站(RR-9100型)自動監(jiān)測記錄氣象數(shù)據(jù)(步長1 h)，測量指標包括空氣溫度、空氣相對濕度、風速、太陽輻射，試驗數(shù)據(jù)采集時間段為2018年8月31日至2020年3月15日，包含章姬草莓兩個完整生育階段。

章姬草莓冠層高度和葉面積指數(shù)每5 d測量一次，待測樣本植株每次隨機選取50株取平均值。測量所選取植株所有有效長度和有效寬度，根據(jù)實測折算系數(shù)計算每株草莓的葉面積并取平均值，最后根據(jù)種植密度得到葉面積指數(shù)。實測蒸散量根據(jù)水量平衡原理，采用稱量法測量草莓植株蒸散量。使用電子天平對桶栽草莓進行連續(xù)稱量，試驗桶內(nèi)徑0.5 m，高0.5 m，表面覆蓋和棚內(nèi)草莓相同的地膜。桶栽草莓共6桶，均勻分布在一個簡易大棚中。實測蒸散量為6桶草莓蒸散量的平均值。

模型精度評價使用平均相對誤差(MAE)、決定系數(shù)(R2)和威爾莫特一致性指數(shù)(D)，模型精度可靠的標準為R2和D均大于0.8且MAE小于20%。使用Bayesian方法估算PM模型中參數(shù)的后驗分布，并分析參數(shù)的不確定性。

1.3 PM模型

PM模型是一種基于能量平衡原理的蒸散機理模型，是目前溫室蒸散模擬應用最廣泛的機理模型[3-4]，PM模型一般形式為

(1)

式中λET——潛熱通量，W/m2

Δ——飽和水汽壓-溫度關系曲線在平均氣溫T處的切線斜率，kPa/K

Rn——凈輻射通量，W/m2

G——土壤熱通量，W/m2

ρa——空氣密度，取1.29 kg/m3

Cp——空氣定壓比熱容，取1 004 J/(kg·K)

es——飽和水汽壓，kPa

ea——實際水汽壓，kPa

rc——冠層邊界層阻力，s/m

γ——濕度計常數(shù)，取0.067 kPa/K

ra——空氣動力學阻力，s/m

其中，ra和rc是兩個關鍵參數(shù)。

1.3.1空氣動力學阻力模型

1.3.1.1Perrier公式法

ra模型包括PERRIER[22]提出的Perrier公式，此公式已在通風條件較好的大田環(huán)境下得到普遍認可[18,23]。多數(shù)研究認為該公式在溫室低風速環(huán)境中并不適用[10-13]。最終Perrier公式可以表達為僅與冠層高度hc和參考高度處風速uz相關，計算式為

(2)

1.3.1.2Thom風速零值法

(3)

1.3.1.3熱傳輸系數(shù)法

溫室內(nèi)水汽以渦流方式擴散，ra可以通過熱傳輸系數(shù)計算[26]。溫室的對流類型可分為自由對流、強迫對流和混合對流，計算熱傳輸系數(shù)hs時先判定出對流類型，之后根據(jù)不同對流類型的公式計算對應的hs[27-28]。表1列出了熱傳輸系數(shù)法中3種對流類型出現(xiàn)的比例，以及各類對流類型條件下與hs相關的因子。由表1可知，溫室通風時絕大部分時間溫室內(nèi)的對流都屬于強迫對流，強迫對流條件下冠層高度風速u與hs正相關，葉片特征長度dc與hs負相關；不通風時絕大部分時間屬于混合對流，冠層溫度與空氣溫度的差值Δt、u分別與hs正相關。

表1 3種對流出現(xiàn)比例以及與hs相關的因子Tab.1 Proportion of three kinds of convection and hs related factors

將3類對流條件下計算出的熱傳輸系數(shù)hs代入可計算得到ra[26]，即

(4)

式中LAI——葉面積指數(shù)

1.3.2冠層邊界層阻力模型

1.3.2.1基于葉片氣孔阻力的冠層邊界層阻力模型

BAILEY等[26]通過分析葉片氣孔阻力rs與氣象因子的響應，構建了基于有效葉面積指數(shù)LAIe的冠層邊界層阻力rc模型，計算式為

(5)

(6)

rs是確定rc的關鍵參數(shù)，由氣象因子、水分條件所決定[29-30]。本研究簡易大棚章姬草莓實施充分灌溉，對rs與簡易大棚內(nèi)氣象因子的相關關系分析結果顯示，太陽輻射Rs為影響rs的主要氣象因子，這一結果與文獻[31-32,10]對溫室番茄和甜椒的研究結果相似。圖1為最終通過非線性擬合得出的Rs與rs的關系，其關系式為

圖1 簡易大棚章姬草莓rs與Rs的關系Fig.1 Relationship between rs and Rs of Akihime strawberry in basic greenhouse

(7)

綜合式(5)～(7)，最終構建出包含氣象因子Rs和作物LAIe的雙因子rc模型。

1.3.2.2基于冠層導度的冠層邊界層阻力模型

冠層邊界層阻力rc是冠層導度gc的倒數(shù)，通過構建氣象因子、土壤含水率與冠層導度gc的關系

(8)

式中gmax——冠層頂部葉片最大氣孔導度，m/s

KQ——短波輻射消光系數(shù)

Qh——冠層頂部可見光通量密度，W/m2

Q50——冠層導度達到最大值一半時的可見光通量，W/m2

VPD——空氣飽和水汽壓差，kPa

D50——冠層導度達到最大值一半時的飽和水汽壓差，kPa

F2——水分脅迫系數(shù)

現(xiàn)場觀摩結束后，新洋豐農(nóng)藝師在室內(nèi)會議室，為參會人員詳細地講解茄子高產(chǎn)高效栽培技術要點及洋豐百倍邦套餐肥的優(yōu)勢與特點，農(nóng)戶們聽后激動不已。海南荊島公司張斌經(jīng)理發(fā)布鉅惠訂貨政策后，農(nóng)戶們一算賬，用洋豐百倍幫套餐肥不僅茄子長勢比用進口肥更好，而且按活動政策現(xiàn)場訂肥，一畝地投入還要省幾十甚至上百元。

求倒數(shù)后即可確立冠層邊界層阻力rc2模型[19,33]。本研究中不存在水分脅迫，因此F2取值恒為1。此模型共4個待定參數(shù)，分別為gmax、KQ、Q50、D50。

1.3.2.3基于葉片氣孔導度尺度提升的冠層邊界層阻力模型

冠層邊界層阻力rc與冠層導度gc是倒數(shù)關系，因此構建基于氣孔導度gs的rc模型的一個基本思路是將gs尺度轉換到gc，兩種常見的gs模型是Jarvis模型[34]和Ball模型[35]。本研究參照張寶忠等[36]的方法，基于gs實測值建立gs與氣象因子(光合有效輻射PAR、VPD和空氣溫度Ta[37-43])的Jarvis模型，再以光合有效輻射PAR為尺度轉換因子，構建出基于氣象因子PAR、VPD和Ta的冠層導度gc模型。

gc模型可表示為gs在葉面積指數(shù)上的積分，即

(9)

通過PARa作為尺度轉換因子得到gc模型為

(10)

式中a1、a2、a3、a4——待率定經(jīng)驗參數(shù)

PARa——葉片截獲的光合有效輻射，μmol/(m2·s)

K——消光系數(shù)

PARh——冠層頂部高度的光合有效輻射，μmol/(m2·s)

最后求倒即可得到rc3模型，此模型共5個待定參數(shù)，分別為K、a1、a2、a3、a4。

1.3.3PM模型匯總

將上述PM模型中2個關鍵參數(shù)ra和rc的模型兩兩組合，共得到9種模型(表2)。其中涉及rc2和rc3的6種模型(*標記)使用Bayesian方法進行參數(shù)估計。

表2 PM模型匯總Tab.2 Penman-Monteith models summary

1.4 Bayesian參數(shù)估計方法

PM模型結構較為復雜，子模型rc2和rc3涉及的待率定參數(shù)較多(分別為4個和5個)，模型的不確定性較強，基于上述考慮選擇使用Bayesian方法作為模型參數(shù)率定方法。

子模型rc2中的參數(shù)gmax、KQ、Q50、D50和子模型rc3中的參數(shù)K、a1、a2、a3、a4作為參數(shù)向量β的組成部分，先驗信息見表3。

表3 PM模型參數(shù)先驗分布Tab.3 Prior distribution of parameters in Penman-Monteith models

假設模型模擬值S與實測值O之間的差異即模型誤差獨立且服從方差σ2為未知常數(shù)的正態(tài)分布，模型誤差均值為0[47-48]，則包含n個觀測值的實測數(shù)列似然函數(shù)為[19]

(11)

式中Oi——第i個蒸散實測值

xi——模型輸入向量

S(xi,β)——參數(shù)向量β下模型輸出值

σ——模型誤差的標準差

后驗分布抽樣使用基于Gibbs抽樣的馬爾科夫-蒙特卡羅(MCMC)方法，使用WinBUGS軟件實現(xiàn)，運行3條平行MCMC鏈，迭代次數(shù)30萬次并酌情增加。

2 結果與討論

2.1 簡易大棚日尺度空氣動力學阻力ra模型篩選

由表4可知，ra模型中ra1表現(xiàn)最差，ra2與ra3表現(xiàn)相當；分析ra1-rc3、ra2-rc3、ra3-rc3，可得出相同的結論。ra1表現(xiàn)差的原因是ra1模型中風速在分母，溫室環(huán)境中日均風速較低，因此會高估ra從而影響最終模型的估計效果，因此溫室使用ra1結果較好的報道很少[49]，一些研究甚至建議不在溫室中使用這種方法[10-13]。但本研究中ra1-rc3模型完全符合3個精度指標，ra1-rc2模型僅率定年D指標稍微超出要求。因此，簡易大棚日尺度條件下雖然ra1較ra2和ra3表現(xiàn)差，但并非不可用，這一結論有別于現(xiàn)有溫室環(huán)境的報道。Perrier公式法能在研究區(qū)使用的原因是棚栽草莓通風時間較長，另外簡易大棚的密閉性較高端溫室差，存在漏風情況，這些原因導致研究區(qū)簡易大棚內(nèi)部日均風速較一般溫室高，進而使Perrier公式計算結果不會被過高估計。本研究統(tǒng)計的棚內(nèi)平均風速為0.36 m/s(2018率定年與2019檢驗年2年平均值)，接近已有文獻中日光溫室風速日均值上限(0～0.36 m/s)[50]，使用Perrier方法計算出的ra1為253.41，與已有的ra報道差別不大[29]。但由于ra1模型與ra2模型的復雜程度相當且遠小于ra3，ra2模型與ra3模型精度相當，因此從模型精度與復雜程度兩個角度考慮，研究區(qū)簡易大棚日尺度條件下推薦使用ra2模型。

表4 PM模型蒸散模擬精度評價Tab.4 Accuracy evaluation of evapotranspiration simulation in Penman-Monteith models

2.2 簡易大棚日尺度冠層邊界層阻力rc模型篩選

簡易大棚日尺度蒸散量模擬值與實測值相關關系如圖2所示，可直觀看出除ra1-rc1和ra2-rc1擬合效果較差以外，其余各PM蒸散模型的日蒸散量模擬值與實測值能較好地沿1∶1直線分布。ra3-rc1雖然能較好沿直線1∶1分布，但散點明顯更分散，決定系數(shù)R2也較低(0.75)。分析PM蒸散模型精度評價(表4)，發(fā)現(xiàn)包含rc1子模型的PM蒸散模型表現(xiàn)最差?？赡艿脑蚴莚c1子模型考慮的作物-環(huán)境因子不全面，即模型缺少關鍵的作物-環(huán)境因子，進而導致一些關鍵的生理-物理過程被忽略，最終影響模型精度。rc1涉及的因子僅有太陽有效輻射和葉面積指數(shù)，rc2較rc1增加了描述冠層結構的參數(shù)項并增加了棚內(nèi)水汽壓虧缺這一環(huán)境因子，rc3較rc1增加了溫度、飽和水汽壓虧缺等環(huán)境因子。rc1較rc2、rc3缺少的這些因子都被證明是估算冠層導度的關鍵因子[51]。

圖2 簡易大棚日尺度蒸散量模擬值與實測值相關關系Fig.2 Correlations between simulated values of daily evapotranspiration and measured values in basic greenhouse

分析表4中rc子模型rc2和rc3，發(fā)現(xiàn)rc3較rc2性能更好(ra2-rc3優(yōu)于ra2-rc2，ra3-rc3優(yōu)于ra3-rc2)。對比rc3與rc2兩個子模型，發(fā)現(xiàn)rc3較rc2多考慮了Ta這一環(huán)境因子，但子模型rc3優(yōu)于rc2的原因并不一定是因為Javrvis模型中rc3多考慮了Ta因子。CHEN等[52]通過對包含不同數(shù)量關鍵環(huán)境因子(1～4個)、不同約束函數(shù)形式共計191個Javrvis模型進行了蒸騰模擬研究，發(fā)現(xiàn)Javrvis模型中嵌套增加關鍵環(huán)境因子并不總是能提高模型性能，模型性能還取決于Javrvis模型的約束函數(shù)形式。因此，rc3精度更高的一個可能原因是本文在組建Javrvis模型時所選擇的關鍵環(huán)境因子以及其約束函數(shù)形式綜合起來更加適合研究區(qū)簡易大棚的環(huán)境。此外，rc3子模型較rc2子模型包含有更多的待率定經(jīng)驗參數(shù)，與之對應的是包含有更多的先驗分布信息，這意味著Bayesian分析方法在參數(shù)率定的過程中可以利用更多的先驗分布信息，這對減小模型不確定性、提高模型性能具有積極意義。因此，rc3精度更高的另一個可能原因是Bayesian參數(shù)率定過程中先驗信息更加充分。

rc3子模型中以光合有效輻射PARa為尺度轉換因子，以光合有效輻射PARa和飽和水汽壓差VPD構建的冠層導度估算模型被證明能較好地實現(xiàn)葉片氣孔導度向冠層導度的尺度轉化提升[38]，但此模型中包含有難以測量的因子——消光系數(shù)K，K在測量過程中容易產(chǎn)生測量誤差，誤差累積到達一定數(shù)量必然會影響模型精度。因此，當模型中包含難以實測的因子時，在先驗信息充分的條件下建議嘗試將此因子當作未知參數(shù)，結合Bayesian方法計算出這些難以實測因子的估計值。

2.3 簡易大棚日尺度不同生育階段的參數(shù)不確定性分析

同時滿足3個精度指標的PM模型共有5個(ra1-rc3、ra2-rc2、ra2-rc3、ra3-rc2和ra3-rc3)，圖3所示5個模型中涉及的大部分參數(shù)的95%置信度后驗分布區(qū)間在不同生育階段與先驗分布區(qū)間相比明顯縮小，這說明Bayesian方法在參數(shù)率定過程中減小了模型參數(shù)的不確定性，參數(shù)率定效果較好。某些參數(shù)僅在特定的生育階段率定效果較好，在膨果期中參數(shù)D50、KQ、Q50和gmax率定效果均較好。在模型ra2-rc2和ra3-rc2中，參數(shù)D50的后驗分布區(qū)間較先驗分布區(qū)間分別縮小了96.44%和93.98%，參數(shù)KQ的后驗分布區(qū)間較先驗分布區(qū)間分別縮小了56.08%和24.15%，參數(shù)Q50的后驗分布區(qū)間較先驗分布區(qū)間分別縮小了97.78%和99.50%，參數(shù)gmax的后驗分布區(qū)間較先驗分布區(qū)間分別縮小了99.75%和99.68%。僅特定生育階段參數(shù)率定效果好的原因是先驗信息不夠充分[53]。本研究參照已有文獻中參數(shù)的先驗信息均沒有區(qū)分生育階段(表3)，這是造成先驗信息不充分的一個原因。例如葉面積LAI在草莓生長初期隨著植株的快速生長葉面積增大較快，花期時葉面積才趨于穩(wěn)定，如果使用全生育期LAI均值反推出的參數(shù)作為先驗信息，必然會影響到Bayesian參數(shù)率定結果。本研究中參數(shù)gmax隨生育期變化明顯，在草莓生長初期參數(shù)值較大，在膨果期驟降至極小值，這一變化規(guī)律與青藏高原高山草甸[19]、西北地區(qū)胡楊林[54]和美國威斯康星州草地[55]相同。因此，D50、KQ、Q50和gmax這類參數(shù)與植株生理特性和環(huán)境輻射關系密切，對作物生育階段的變化比較敏感，使用區(qū)分生育階段的先驗信息進行參數(shù)率定效果會更好。

圖3 2018年PM蒸散模型日尺度不同生育階段參數(shù)最優(yōu)估計值和95%置信區(qū)間Fig.3 Optimal estimation and 95% confidence interval of parameters in different growth stages of PM models in 2018

參數(shù)a2率定效果最好，在草莓全生育階段中參數(shù)a2后驗分布區(qū)間都較先驗分布區(qū)間明顯縮小。ra1-rc3模型中參數(shù)a2的后驗分布區(qū)間較先驗分布區(qū)間在4個生育階段分別縮小了72.32%、88.59%、56.10%和86.32%，ra2-rc3模型分別縮小了97.65%、92.38%、93.31%和98.24%，ra3-rc3模型分別縮小了97.65%、92.38%、93.31%和98.24%。參照對D50、KQ、Q50和gmax這類參數(shù)有關生育階段敏感性的分析，同樣使用不夠充分的先驗信息卻得到了在全生育階段都較好的參數(shù)率定結果，說明參數(shù)a2對作物生育階段的變化不敏感。此外，還可以從rc3模型分析得出相同的結論。本研究使用的rc3模型中參數(shù)a2僅與VPD相關，VPD反映了空氣的干燥程度，很明顯這一環(huán)境因子在簡易大棚環(huán)境中不隨生育階段發(fā)生顯著變化。

綜上，結合模型模擬精度評價與模型參數(shù)不確定性分析，篩選出ra2-rc3和ra3-rc3兩種適宜研究區(qū)簡易大棚日尺度下的草莓PM蒸散模型，由于ra3-rc3模型涉及的熱傳輸系數(shù)法計算繁瑣，故推薦優(yōu)先使用ra2-rc3模型。

3 結論

(1)Bayesian參數(shù)估計方法使PM模型在簡易大棚日尺度蒸散估算精度較高，應用Bayesian參數(shù)估計方法的6種PM模型中有5種在模型率定年和模型檢驗年同時滿足3個精度指標。

(2)Bayesian參數(shù)估計方法能有效減小簡易大棚日尺度PM模型中部分參數(shù)的不確定性，在模型ra2-rc3中參數(shù)a2的后驗分布區(qū)間較先驗分布區(qū)間在4個生育階段分別縮小了97.65%、92.38%、93.31%和98.24%，在模型ra2-rc2中參數(shù)D50、KQ、Q50和gmax在膨果期后驗分布區(qū)間較先驗分布區(qū)間分別縮小了96.44%、56.08%、97.78%和99.75%。

(3)篩選出ra2-rc3和ra3-rc32種適宜簡易大棚日尺度蒸散模擬的PM模型，其中最優(yōu)模型為ra2-rc3。