基于Penman-Monteith模型的低丘紅壤區(qū)稻田蒸散模擬?

文建川，景元書**，韓麗娟

（1．南京信息工程大學(xué)氣象災(zāi)害預(yù)報(bào)預(yù)警與評(píng)估協(xié)同創(chuàng)新中心/應(yīng)用氣象學(xué)院，南京 210044；2．國家氣象中心，北京 100081）

生態(tài)系統(tǒng)中植被發(fā)育以及系統(tǒng)生產(chǎn)力都受到蒸散的直接影響，同時(shí)它不僅是土壤-植被-大氣連續(xù)體中水分傳輸中的重要環(huán)節(jié)，也是生態(tài)系統(tǒng)中能量平衡的主要組成部分[1]。據(jù)統(tǒng)計(jì)，以蒸散的形式消耗農(nóng)業(yè)用水占到90%[2-3]，蒸散的準(zhǔn)確估算有利于作物水分生產(chǎn)力提升以及改善農(nóng)田或流域灌溉管理[4-5]，因此，蒸散觀測(cè)和模擬一直受到生態(tài)氣象領(lǐng)域?qū)W者的廣泛關(guān)注。

近幾十年來，蒸散估算模型眾多，在該研究領(lǐng)域較為有名的是以能量平衡和微氣候?qū)W方法為原理的Penman-Monteith模型（簡稱PM模型），該模型模擬精度主要取決于冠層阻力這一關(guān)鍵性參數(shù)[7]。國內(nèi)外研究者基于冠層阻力公式，對(duì)其與PM模型結(jié)合后模擬蒸散量開展了不少研究。Irmak 采用非線性回歸方法設(shè)計(jì)了7個(gè)形式相同，氣象因子數(shù)量不同導(dǎo)致復(fù)雜程度不一樣的模型，之后利用該模型對(duì)大豆蒸散量進(jìn)行模擬，取得了較好的結(jié)果[8-9]。吳林等[10]將Irmak 公式應(yīng)用到PM模型中，可以很好地估算黑河綠洲區(qū)玉米半小時(shí)尺度上的蒸散量。趙華等[11]采用Jarvis 公式計(jì)算了南京地區(qū)2012年和2013年水稻拔節(jié)-成熟期的冠層阻力，并將其耦合到PM模型中模擬水稻該時(shí)期的蒸散，兩年的擬合度都超過0.95。Li 等[12]采用PM_Jarvis模型時(shí)，基于玉米葉面積指數(shù)大小將整個(gè)生育期劃分為兩個(gè)階段，分段擬合模型中的經(jīng)驗(yàn)系數(shù)，對(duì)比后發(fā)現(xiàn)以葉面積指數(shù)為0.5 分段擬合的系數(shù)進(jìn)行蒸散模擬的效果最佳。Srivastava等[13]比較了Monteith、Katerji-Perrier、Todorovic和Jarvis 四種冠層阻力公式應(yīng)用在PM模型中的模擬效果，結(jié)果表明Jarvis和Todorovic 都較好地估算了半濕潤地區(qū)玉米的蒸散。Jarvis和Irmak 公式中包含了較多經(jīng)驗(yàn)系數(shù)，有學(xué)者分別將其它研究的經(jīng)驗(yàn)系數(shù)值和通過非線性擬合得到的系數(shù)值代入到模型中估算蒸散量，比較后發(fā)現(xiàn)借用其它研究的經(jīng)驗(yàn)系數(shù)得到的模擬值與實(shí)測(cè)值相關(guān)性較差[14]。

因?yàn)椴煌貐^(qū)氣候差異大，這兩種模型同屬為經(jīng)驗(yàn)?zāi)Ｐ停琂arvis和Irmak 公式應(yīng)用在PM模型中估算蒸散量，能否取得較高的模擬精度，是否適用于低丘紅壤區(qū)稻田蒸散研究，有待進(jìn)一步探索。因此，需要針對(duì)試驗(yàn)區(qū)及作物，挑選出最為合適的冠層阻力公式，得到模擬精度高的蒸散模型。另外，Jarvis模型中各脅迫函數(shù)存在多種表達(dá)式，導(dǎo)致有多種組合形式，較多學(xué)者在使用該模型時(shí)，基本上都是選取其中一種組合形式進(jìn)行模擬，但是不同的表達(dá)式所涉及的參數(shù)數(shù)量不同，可能導(dǎo)致最后的模擬精度有差別。趙華[15]采用了Jarvis模型中的幾種組合形式，對(duì)冬小麥和水稻蒸散量進(jìn)行研究，但只選取了各生育期中典型晴天數(shù)據(jù)，未對(duì)作物全生育期進(jìn)行驗(yàn)證。因此，本研究以波文比儀監(jiān)測(cè)該地區(qū)稻田的蒸散量作為實(shí)測(cè)值，對(duì)比分析PM_Jarvis1-8（脅迫函數(shù)的8種組合形式）與PM_Irmak 共9個(gè)模型的蒸散模擬效果，旨在選出更適用于低丘紅壤區(qū)稻田的實(shí)際蒸散模型，準(zhǔn)確估算蒸散量，以利于制定更合理的灌溉計(jì)劃和提高稻田的水分利用效率。

1 材料與方法

1.1 試驗(yàn)站點(diǎn)概況

試驗(yàn)站點(diǎn)位于江西省鷹潭市中國科學(xué)院紅壤生態(tài)試驗(yàn)站(116°55'E，28°15'N)，該站地處武夷山至潘陽湖平原的過渡地帶鷹潭盆地，屬亞熱帶濕潤季風(fēng)氣候，年平均氣溫17.6℃，降水量1788.8mm，日照時(shí)數(shù)1739.4h。試驗(yàn)觀測(cè)時(shí)間為2014年和2015年晚稻生長季（7-10月），水稻移栽后，依照當(dāng)?shù)卮筇锷a(chǎn)進(jìn)行田間管理。觀察期內(nèi)水稻生育期見表1。

表1 試驗(yàn)區(qū)2014年和2015年水稻生育期Table1 Growth stage of rice in experimental area in 2014 and 2015

1.2 數(shù)據(jù)來源

波文比儀架設(shè)在稻田中部，該系統(tǒng)由NR-Lite凈輻射傳感器(Kipp&Zonen，Netherland)、兩層空氣溫濕度傳感器( HMP155A)、CR1000 數(shù)據(jù)采集器、風(fēng)速儀( 010C-1）組成，自動(dòng)采集1.5m 高度差之間的溫濕度差，計(jì)算得到波文比值。土壤熱通量由埋深為 5cm的土壤熱通量板(Hukseflux，HFP01，Netherland)測(cè)得。氣象數(shù)據(jù)由布設(shè)在農(nóng)田中部的自動(dòng)氣象站（HOBO U30，ONSET，USA）監(jiān)測(cè)，觀測(cè)內(nèi)容包括相對(duì)濕度、氣溫、2m 風(fēng)速、氣壓等小氣候數(shù)據(jù)，觀測(cè)頻率為20min。

1.3 數(shù)據(jù)處理方法

波文比測(cè)定水熱通量是以地表能量平衡方程為依據(jù)，即

式中，Rn 為凈輻射通量（W·m-2），λET 為潛熱通量（W·m-2），H 為感熱通量（W·m-2），G 為土壤熱通量（W·m-2）。其中，Rn和G可以實(shí)測(cè)得到，λET和H 通過計(jì)算得到。

波文比（β）為感熱通量和潛熱通量的比值，根據(jù)莫寧-奧布霍夫理論，假定水汽和熱量的湍流擴(kuò)散系數(shù)Kh與Kw相等，則有

式中，γ 為干濕表常數(shù)，取值為0.067kPa·℃；ΔT、Δe 分別為波文比系統(tǒng)的兩層觀測(cè)高度的溫度差和水汽壓差。

波文比儀測(cè)定出兩層高度的溫度差和水汽壓差、凈輻射通量以及土壤熱通量板測(cè)定的土壤熱通量數(shù)據(jù)，結(jié)合式（1）和式（2），可得潛熱通量λET 為

式中，λ 為水的汽化潛熱，取值為2.45MJ·kg-1，ET 為蒸散量(mm)。

因?yàn)槿粘?、日落以及降水天氣，可能造成誤差，所以根據(jù)Perez 等[16]提供的方法確定無效β 值，將其剔除并進(jìn)行插補(bǔ)。計(jì)算式為

式中，δ1和δ2分別為溫、濕度傳感器的測(cè)量精度。

1.4 模型介紹

1.4.1 冠層阻力計(jì)算模型

冠層阻力rc，采用兩種模型進(jìn)行計(jì)算。

（1）PM_Jarvis模型

Jarvis 認(rèn)為冠層阻力與環(huán)境因子具有協(xié)同關(guān)系，而且不同環(huán)境因子對(duì)冠層阻力的影響?yīng)毩ⅲ虼?，提出了冠層阻力與環(huán)境因子之間的關(guān)系式，即

式中，rmin為最小冠層阻力，根據(jù)前人[17]研究結(jié)果，水稻在拔節(jié)期、孕穗期、抽穗期、乳熟期的rmin分別取為40、45、50、60s·m-1，LAIe為有效葉面積指數(shù)，F(xiàn)(Rn)、F(D)、F(T)分別為太陽輻射脅迫函數(shù)、飽和水汽壓差脅迫函數(shù)、溫度脅迫函數(shù)，因以上3種脅迫函數(shù)均有多種表達(dá)式，導(dǎo)致存在不同組合形式，不同的表達(dá)形式有可能影響模擬精度，因此，每個(gè)脅迫函數(shù)均選取其中兩種表達(dá)式。

有效葉面積指數(shù)LAIe為[18]

太陽輻射脅迫函數(shù)分別采用線性關(guān)系式[19]和指數(shù)關(guān)系式[12]，即

飽和水汽壓差脅迫函數(shù)分別采用雙曲線關(guān)系式[19]和線性關(guān)系式[12]，即

空氣溫度脅迫函數(shù)分別采用以下兩種關(guān)系式[12，19]，即

式中，Rn為凈輻射通量（W·m-2），D 為飽和水汽壓差（kPa），T 為空氣溫度(℃)，a1、a2、a3、a4、b1、b2、b3、c1、c2、c3、c4為模型經(jīng)驗(yàn)系數(shù)。

基于上述各脅迫函數(shù)的響應(yīng)表達(dá)式，得到8種不同的Jarvis模型組合形式（表2）。

（2）PM_Irmak模型

Irmak 等采用廣義非線性回歸方法，構(gòu)建了關(guān)于氣象因子、葉面積指數(shù)與空氣動(dòng)力學(xué)阻力的冠層阻力模型，即

式中，rc為冠層阻力（s·m-1），Rn 為凈輻射（W·m-2），T 為空氣溫度（℃），RH 為相對(duì)濕度，U 為風(fēng)速（m·s-1），LAI 為葉面積指數(shù)，ra為空氣動(dòng)力學(xué)阻力（s·m-1），按Rana 等提出的公式計(jì)算[20]，即

表2 PM_Jarvis 冠層阻力模型中3種脅迫函數(shù)表達(dá)式的不同組合Table2 Different combinations of three stress function expressions of the PM_Jarvis canopy resistance model

式中，z 為參考高度（m），取值為1.5m；h 為冠層高度（m），以實(shí)測(cè)株高代替；k 為Von Karman常數(shù)，取值為0.41；u 為風(fēng)速（m·s-1），z0=0.13h，d=0.63h。

1.4.2 蒸散量計(jì)算模型

Monteith 在Penman 基礎(chǔ)上，引入了空氣動(dòng)力學(xué)阻力和冠層阻力等參數(shù)，提出以能量平衡和水汽擴(kuò)散理論為基礎(chǔ)的適用于作物蒸騰計(jì)算的阻力模式，即Penman-Monteith模型。

式中，λET 為潛熱通量（W·m-2）；Rn 為凈輻射通量（W·m-2）；G 為土壤熱通量（W·m-2）；Δ 為飽和水汽壓曲線斜率（kPa·K-1）；ρ 為空氣密度，取值為 1.29kg·m-3；Cp 為空氣定壓比熱，取值為1004J·kg-1·K-1；D 為飽和水汽壓差（kPa）；γ 為干濕表常數(shù)，取值為0.067kPa·K-1；rc為冠層阻力（s·m-1）；ra為空氣動(dòng)力學(xué)阻力（s·m-1）。

1.5 模型精度分析

以波文比儀測(cè)定的蒸散量ET 作為參照，分別把不同冠層阻力模型應(yīng)用到Penman-Monteith 公式中，再將模擬得到的ET 與實(shí)測(cè)值比較。模型的模擬精度評(píng)價(jià)指標(biāo)包括均方根誤差（RMSE)、平均絕對(duì)誤差（MAE）、納什效率系數(shù)（NSE）和一致性指數(shù)（IA）、確定系數(shù)（R2）。

式中，Pi為模擬蒸散量(mm·h-1)，Oi為實(shí)測(cè)蒸散量(mm·h-1)，為實(shí)測(cè)蒸散量的平均值(mm·h-1)，為模擬蒸散量的平均值(mm·h-1)，n 為個(gè)數(shù)。

RMSE和MAE 反映模擬值與實(shí)測(cè)值之間的差異，數(shù)值越接近0，表明誤差越?。籒SE 越接近1，表明模型可信度越高；IA 反映模擬值與實(shí)測(cè)值的符合程度，數(shù)值越接近1，說明模型的預(yù)測(cè)效果越好。R2越接近1，表明擬合程度越好。

2 結(jié)果與分析

2.1 PM_Jarvis和PM_Irmak 冠層阻力模型模擬效果

基于2014年波文比儀監(jiān)測(cè)的水稻全生育期實(shí)測(cè)蒸散量與同階段的氣象數(shù)據(jù)，將表2中組合好的8種Jarvis模型及Irmak模型分別與Penman-Monteith模型結(jié)合，簡稱PM_Jarvis1-8模型（脅迫函數(shù)8種組合形式）、PM_Irmak模型，利用1stOpt 軟件進(jìn)行非線性擬合(http: //www.7d-soft.com/en/)，通過全局優(yōu)化算法，得到各模型的經(jīng)驗(yàn)系數(shù)值（表3）。統(tǒng)計(jì)這9種蒸散模型的誤差及相關(guān)性統(tǒng)計(jì)參量，結(jié)果見表4。從表4中可知，不同組合形式的PM_Jarvis模型精度指標(biāo)都存在一定的差異，其中PM_Jarvis6模型精度最低，其RMSE、MAE 最高，NSE、IA 及R2最低。PM_Irmak模型的RMSE、MAE 在9種模型中最低，NSE、IA 及R2較其余模型更接近1，表明該模型模擬精度最高。就PM_Jarvis模型而言，8種組合形式模型的模擬精度存在差別，將模型進(jìn)行兩兩對(duì)比，結(jié)果表明，對(duì)于太陽輻射脅迫函數(shù)，F(xiàn)1(Rn)的表達(dá)式在應(yīng)用于該模型中較F2(Rn)效果更佳；對(duì)于飽和水汽壓差脅迫函數(shù)和空氣溫度脅迫函數(shù)，F(xiàn)1(D)、F1(T)的表達(dá)式分別優(yōu)于F2(D)和F2(T)。

表3 PM_Jarvis（式7-12）和PM_Irmak（式13）模型中各參數(shù)的擬合值Table3 Simulation value of parameters in PM_Jarvis（Equation7-12） and PM_Irmak models（Equation13）

表4 PM_Jarvis和PM_Irmak 模擬2014年逐時(shí)蒸散的統(tǒng)計(jì)參量Table4 Statistical parameters of hourly evapotranspiration simulated by PM _Jarvis and PM_Irmak models

2.2 PM_Jarvis和PM_Irmak 模擬結(jié)果驗(yàn)證

利用2015年7月23日-10月31日的氣象資料和水稻生理數(shù)據(jù)，代入9種模型計(jì)算蒸散量，并與波文比儀監(jiān)測(cè)的該時(shí)段蒸散量比較，進(jìn)一步驗(yàn)證模型。從圖1看出，PM_Irmak模型模擬的蒸散量與實(shí)測(cè)蒸散量基本沿1:1 線分布，且散點(diǎn)集中，表明模擬值與實(shí)測(cè)值較接近；其余8種PM_Jarvis模型散點(diǎn)大都位于1:1 線下方，總體上表現(xiàn)出偏低趨勢(shì)。為了更全面地體現(xiàn)模型模擬效果，在全生育期及各生育階段進(jìn)行精度分析，結(jié)果見表5和圖2。由表5可以看出，PM_Irmak模型在估算水稻全生育期蒸散量中，各評(píng)價(jià)指標(biāo)整體均優(yōu)于其余模型，同樣，就PM_Jarvis模型而言，脅迫函數(shù)表達(dá)式不同導(dǎo)致8種PM_Jarvis模型精度不一，與2014年結(jié)果一致。

圖1 PM-Jarvis和PM-Irmak模型模擬逐時(shí)蒸散量與實(shí)測(cè)蒸散量對(duì)比Fig.1 Comparison of hourly evapotranspiration( ET) between measured and simulated by PM-Jarvis and PM-Irmak models

表5 PM-Jarvis和PM-Irmak模型模擬的稻田全生育期蒸散量精度驗(yàn)證（2015年）Table5 Accuracy of evapotranspiration during the whole growth stage of rice paddy simulated by PM-Jarvis and PM-Irmak models (2015)

將水稻的全生育期分為4個(gè)階段，由圖2可見，8個(gè)PM_Jarvis模型的4個(gè)評(píng)價(jià)指標(biāo)在返青-分蘗期與后 3個(gè)生育期有明顯差距，在該生育階段，PM_Irmak模型的均方根誤差（RMSE）和平均絕對(duì)誤差（MAE）顯著低于8個(gè)PM_Jarvis模型，納什效率系數(shù)（NSE）和一致性指數(shù)（IA）明顯較高；在后3個(gè)生育時(shí)期，PM_Irmak模型的RMSE和MAE雖然不是最低，NSE和IA 也不是最高，但與其余模型差距較小。所以認(rèn)為，PM_Irmak模型總體上對(duì)水稻整個(gè)生長階段蒸散量的模擬效果較好。

圖2 PM-Jarvis和PM-Irmak模型模擬的稻田不同生育期蒸散量精度驗(yàn)證Fig.2 Accuracy of evapotranspiration in the different growth stages of rice paddy simulated by PM-Jarvis and PM-Irmak models

2.3 PM-Jarvis和PM-Irmak模型模擬蒸散量的日動(dòng)態(tài)變化

在水稻4個(gè)生育時(shí)期中，每一階段均選擇連續(xù)9d 用來觀察9個(gè)模型估算的蒸散與實(shí)測(cè)蒸散的日變化，結(jié)果見圖3（下頁）。由圖可以看出，在返青-分蘗期（7月23-31日），PM_Jarvis1-8模型呈現(xiàn)出低估現(xiàn)象，并且低估程度嚴(yán)重，蒸散估算值不及實(shí)際蒸散量的50%；而PM_Irmak模型在該生育時(shí)期，除了在個(gè)別時(shí)間點(diǎn)（7月27日）蒸散模擬值稍偏高，其余時(shí)段均與實(shí)測(cè)值較接近，模擬效果明顯優(yōu)于PM_Jarvis1-8模型。在后3個(gè)生育時(shí)期，各模型模擬的日變化動(dòng)態(tài)趨勢(shì)同步，數(shù)值上存在差異，但差異不大，均與實(shí)測(cè)蒸散的日動(dòng)態(tài)曲線基本吻合?？偟膩碚f，PM_Irmak模型在水稻整個(gè)生長階段都能較好地估算蒸散量，所以更適用于紅壤地區(qū)稻田蒸散量模擬。

3 結(jié)論與討論

3.1 結(jié)論

本研究利用不同冠層阻力公式的PM模型估算蒸散量，與波文比儀監(jiān)測(cè)的實(shí)測(cè)蒸散量對(duì)比，PM_Jarvis 1-8和PM_Irmak 共9種模型中，PM_Irmak模型的均方根誤差、平均絕對(duì)誤差比8種PM_Jarvis模型的均偏小，確定系數(shù)、納什效率系數(shù)和一致性指數(shù)相比更接近 1，整體模擬效果最好。8種PM_Jarvis模型在水稻返青分蘗期的初始階段存在明顯低估現(xiàn)象，而PM_Irmak模型在水稻所有生育階段的模擬值均與蒸散日變化動(dòng)態(tài)值接近。

3.2 討論

圖3 PM_Jarvis和PM_Irmak模型模擬值與實(shí)測(cè)值的蒸散日變化動(dòng)態(tài)Fig.3 Daily variation of observed and simulated evapotranspiration by PM_Jarvis and PM_Irmak models

因Jarvis 公式脅迫函數(shù)表達(dá)式不止一種，本研究每個(gè)函數(shù)各選取兩種，組成了8種組合形式，發(fā)現(xiàn)PM_Jarvis1-8 模擬效果存在不同程度的差異，其中PM_Jarvis6模型（3個(gè)脅迫函數(shù)分別采用F2(Rn)、F2(D)、 F2(T)的表達(dá)式）精度最低。雖然這8種模型為同一形式，但脅迫函數(shù)的表達(dá)式不同，對(duì)最后蒸散量的模擬產(chǎn)生了較大影響。趙華[15]在進(jìn)行水稻和冬小麥蒸散量研究時(shí)，也將不同組合形式的Jarvis 公式應(yīng)用到PM模型中，發(fā)現(xiàn)不同模型的均方根誤差和擬合度相差較大。8種PM_Jarvis模型，在返青期均呈現(xiàn)出明顯低于實(shí)測(cè)蒸散的現(xiàn)象。Li 等[21]在利用PM_Jarvis模型對(duì)葡萄園蒸散進(jìn)行模擬研究時(shí)，同樣指出在葡萄生長前期階段模擬蒸散也明顯低于實(shí)測(cè)值，隨著葡萄逐漸生長發(fā)育，二者的差異在減小，擬合效果才得以提高。因?yàn)镴arvis 把每個(gè)環(huán)境變量對(duì)作物氣孔產(chǎn)生的影響假定為一種相互獨(dú)立的行為，并認(rèn)為作物所有葉片的氣孔開閉都同時(shí)具有一致性，基于這樣一種假設(shè)提出了該模型[22]，但實(shí)際上環(huán)境變量之間存在著互相影響的關(guān)系，并且較為復(fù)雜，所以這有可能導(dǎo)致了該模型在作物部分生育階段蒸散模擬效果較差。

PM_Irmak模型在返青-分蘗期各評(píng)價(jià)指標(biāo)明顯優(yōu)于8種PM_Jarvis模型，后3個(gè)生育時(shí)期該評(píng)價(jià)指標(biāo)雖不是最優(yōu)，但9種模型之間差異較小，所以最終 PM_Irmak模型仍在水稻全生育期中較 8種PM_Jarvis模型有更好的模擬效果。Li 等在研究干旱地區(qū)玉米蒸散中發(fā)現(xiàn)，PM_Irmak模型在玉米葉面積指數(shù)小于2的時(shí)期模擬蒸散量的誤差要小于葉面積指數(shù)大于2的時(shí)期，雖然在LAI 大于2時(shí)期中模型的精度略低于該研究中的其它冠層阻力模型，但是就整個(gè)生育期而言，PM_Irmak模型的均方根誤差和平均相對(duì)誤差都明顯要小[23]，這與本研究結(jié)果一致。PM_Irmak模型在水稻4個(gè)生育階段估算的蒸散量與實(shí)測(cè)值的日變化動(dòng)態(tài)均較為吻合且數(shù)值接近，精度總體來說較高，說明該模型適用于本試驗(yàn)區(qū)，能夠更好地模擬蒸散量。