王士蘭，王 斌

（1．黑龍江農(nóng)墾勘測設計研究院，哈爾濱 150030；2．東北農(nóng)業(yè)大學 水利與建筑學院，哈爾濱 150030）

FAO－56推薦的Penman－Monteith公式 （以下簡稱P－M公式）計算參考作物蒸發(fā)蒸騰量 （ET0）的精度較高，被公認為計算ET0的標準方法［1］。然而，P－M公式需要的氣象數(shù)據(jù)較多，雖然這些數(shù)據(jù)均為常規(guī)氣象資料，但我國氣象觀測資料的實時共享程度并不高，在氣象數(shù)據(jù)不完備時無法直接應用P－M公式；即使在氣象數(shù)據(jù)完備情況下，也只能應用P－M公式計算ET0而不能直接預測ET0，除非公式中涉及的各種氣象數(shù)據(jù)均為預測（報）值，但實際工作中很難同時獲得如此眾多氣象數(shù)據(jù)的預測 （報）值。

國內(nèi)外關于ET0的預測模型較多，比較典型的有：Tracy等［2］和 Marino等［3］各自提出的基于時間序列分析的ET0預測模型；羅玉峰等［4］提出了ET0預測的傅立葉級數(shù)模型；茆智等［5］提出了預測逐日ET0的方法和指數(shù)模型；郭宗樓等［6］提出了ET0灰色預測模型；顧世祥等［7］、霍再林等［8］、崔遠來等［9］、徐俊增等［10］、彭世彰等［11］應用神經(jīng)網(wǎng)絡對ET0進行了預測；段春青等［12］建立了基于混沌遺傳程序設計的ET0預測模型。由于ET0的日變化較劇烈，時間序列分析模型在預測日步長ET0時精度不是很高，而灰色、神經(jīng)網(wǎng)絡方法缺乏機理。本文引入偏最小二乘回歸 （Partial Least－Squares Regression，PLS）算法［13］， 分 析P－M公式計算ET0時涉及的部分計算變量的推求過程，應用與天氣預報對應的幾種常規(guī)氣象指標觀測數(shù)據(jù)，建立逐日ET0預測模型，以期為預測氣象站點可控范圍內(nèi)的ET0提供一種新方法。

1 建模準備

1．1 ET0計算變量分析

FAO－56推薦計算ET0的P－M 公式為［1］：

式中ET0為參考作物蒸發(fā)蒸騰量，mm／d；T為2 m高處日平均氣溫，℃；Δ為溫度－飽和水汽壓關系曲線在T處的斜率，kPa／℃；Rn為作物表面凈輻射，MJ·m－2·d－1；G為土壤熱通量，MJ·m－2·d－1；γ為濕度表常數(shù)，kPa／℃；u2為2m高處風速，m／s；es為飽和水汽壓，kPa；ea為實際水汽壓，kPa。

從P－M公式計算ET0的具體過程可以看出：

1）一些中間變量 （如：日序數(shù)J、日傾角δ、日地相對距離的倒數(shù)dr）可以應用公式直接計算；

2）在計算地點 （即緯度）確定后，另外部分中間變量 （如：大氣邊緣太陽輻射Ra、日照時數(shù)角Ws、最大可能日照時數(shù)N等）也可以應用公式計算；

3）P－M公式的8項基本計算變量 （T、es、Δ、γ、G、ea、Rn、u2）中有6項 （T、es、Δ、γ、G、ea）與溫度有關。因此，對于固定地點，可以考慮選取ET0的部分計算變量構建預測ET0的回歸模型［14］。

1．2 天氣預報提供的氣象信息及其換算

目前，我國天氣預報業(yè)務比較成熟，在時間上可以向前預報1d或＞1d，在空間上可以預報到縣 （市）級甚至鄉(xiāng) （鎮(zhèn)）級。每日的天氣預報項目一般包括天氣狀況、日最高氣溫、日最低氣溫、風力等級等，這些信息可以很方便地通過各種媒體獲得。天氣狀況與天氣指標對應情況見表1，風力等級與風速對應情況見表2。

表1 天氣狀況與天氣指標對應表［15］Table 1 Weather conditions and weather indiators corresponding table［15］

表2 風力等級與風速對應表［16］Table 2 Wind level wind speed corresponding table［16］

為了利用實測的氣象數(shù)據(jù)建模，需使不同天氣狀況下的n／N以及不同風力等級下的風速變化區(qū)間連續(xù)，因此分別調(diào)整了n／N及風速的變化范圍，根據(jù)表1換算出了與天氣狀況對應的αm，見表3；根據(jù)表2換算出了與風力等級對應的風速um，見表4。這樣處理后，當n／N和風速落入某區(qū)間時，建模應用的實際日照時數(shù)和風速分別取該區(qū)間的nm及um值：

式中N為最大可能日照時數(shù)，可應用日序數(shù)和當?shù)鼐暥阮A先計算，h；αm為修正系數(shù)，可根據(jù)n／N參考表3取值。

表3 天氣狀況與αm值對應表Table 3 Weather conditions andαmvalue corresponding table

表4 風力等級與um值對應表Table 4 Wind level with the umvalue corresponding table

1．3 選擇建模變量

選定Ra、日最低氣溫Tmin、日最高氣溫Tmax、T、Tmin時飽和水汽壓es（Tmin）、Tmax時飽和水汽壓es（Tmax）、飽和水汽壓es、Δ、nm和um為建模自變量，P－M公式計算的日ET0為建模因變量。選用這些自變量的原因如下［14］：

1）除nm、um外，其它自變量均為P－M公式計算ET0時的變量；

2）Ra可以直接計算，它不僅綜合了計算ET0過程中的J、δ、dr、當?shù)氐腤s和N 等信息，還兼顧了Rn的部分計算因子；

3）選取 Tmin、Tmax可以計算es（Tmin）、es（Tmax）、es，通過Tmin和Tmax還可以計算出T、Δ；

4）在應用建成的模型時，不能預知的Tmin、Tmax、nm、um可根據(jù)天氣預報情況確定，例如：某日的天氣預報為 “多云，風力3級，溫度10～15℃”，則Tmin＝10℃，Tmax＝15℃，根據(jù)表3有nm＝ （0．55 N）h，根據(jù)表4得um＝4．40m／s。

綜上，所有的建模自變量均可預先計算或根據(jù)天氣預報確定，當應用前期的氣象實測數(shù)據(jù)率定模型參數(shù)以后，就可以根據(jù)天氣預報情況預測逐日的ET0。

2 偏最小二乘回歸方法在逐日ET0預測中的應用

2．1 建模自變量間的相關性分析

應用黑龍江省建三江氣象臺2001～2006年逐日氣象資料，分析了部分自變量之間的相關系數(shù)（表5），可見各自變量間的相關系數(shù)均較高，這一分析結果對于6a是一致的。此時，常規(guī)多元回歸分析方法已不再適用。PLS方法采用信息綜合與篩選技術，它在變量系統(tǒng)中提取若干對系統(tǒng)具有最佳解釋能力的新綜合變量 （成分），利用這些新提取的成分進行回歸建模。因此，PLS具有普通最小二乘回歸方法所不能比擬的優(yōu)點，恰能較好地解決建模自變量間相關性較高的問題。

2．2 基于PLS的逐日ET0預測模型的建立

黑龍江省大田作物主要生長季為5～9月，其它時段的ET0可計算但對大田農(nóng)業(yè)生產(chǎn)意義不大，因此選取涵蓋主要作物生育期的4～10月為研究時段，應用建三江氣象臺2001～2006年逐日氣象資料，采用P－M公式計算了逐年日ET0序列。根據(jù)PLS的建模步驟［13］，采用Matlab語言編制了PLS算法程序，應用2001年數(shù)據(jù)建立ET0預測模型。由于收集不到2001年的天氣預報數(shù)據(jù)，nm和um利用實測的實際日照時數(shù)和風速數(shù)據(jù)，并根據(jù)表3和表4中劃定的范圍轉化為對應的數(shù)值。經(jīng)檢驗提取7個成分已經(jīng)足夠，最后所得原始變量的逐日ET0預測PLS模型為：

表5 自變量之間的相關系數(shù) （2001年4～10月）Table 5 Between the independent variables，the correlation coefficient（April 2001to October）

2．3 模型擬合及預測效果檢驗

將2001～2006年的建模自變量數(shù)據(jù)輸入建立的PLS模型，應用2001年日樣本資料檢驗模型的擬合效果，應用2002～2006年日樣本資料檢驗模型的預測效果。模型擬合及預測值的合格率見表6，模型輸出的ET0值與P－M公式計算的ET0值對比見圖1、圖2（限于篇幅，僅給出2001年的擬合及2002年預測情況）。

表6 PLS模型擬合及預測值相對誤差合格率Table 6 PLS model fit and predictive value relative error qualified rate ／%

由圖1、圖2及表6可見：PLS模型擬合及預測逐日ET0序列的合格率均較高，從預測平均結果看，逐日ET0預測相對誤差的絕對值＜5%、10%、15%、20% 的合格率分別占40．4%、64．8%、78．1%、87．4%。

3 結論與討論

基于P－M方法計算ET0的過程，分析了該方法涉及的各種計算因子的推求規(guī)律，結合天氣預報信息，在盡可能多地利用可以直接計算的變量情況下，引入PLS方法，建立了預測逐日ET0的PLS模型。主要結論如下：

1）PLS模型程序編制簡單，容易實現(xiàn)，建立的模型參數(shù)固定，應用方便，模擬和預測精度較高，基本能夠滿足農(nóng)業(yè)生產(chǎn)實踐的需要。

2）PLS模型需要的Tmin、Tmax、nm和um4項基本氣象資料可以根據(jù)天氣預報數(shù)據(jù)確定，模型的其它自變量可以直接計算，因此，結合天氣預報的日最高氣溫、日最低氣溫、天氣狀況、風力等級等信息，可以應用本文建立的PLS模型預測逐日的ET0。

3）由于收集不到長系列的天氣預報數(shù)據(jù)，PLS模型是在北方寒區(qū)單站點條件下應用氣象觀測數(shù)據(jù)構建的，研究結論有待進一步檢驗，但該模型及其建模思路具有一定的實用價值和借鑒意義。

［1］Richard G Allen，Luis S Pereira，Dirk Raes，et al．Crop evapotranspiration－guidelines for computing crop water equirements ［M］．Rome：FAO Irrigation and drainage paper 56，1998．

［2］Tracy J C，Marino M A，Taghavi S A．Predicting water demand in agricultural regions using time series forecasts of reference crop evapotranspiration ［A］．Karamouz M．Water Resources Planning and Management，Saving a Threatened Resource－In Search of Solutions［C］．New York：American Society of Civil Engineers，1992．

［3］Marino M A，Tracy J C，Taghavi S A．Forecasting of reference crop evapotranspiration ［J］．Agricultural Water Management，1993，24：163－187．

［4］羅玉峰，崔遠來，蔡學良．參考作物騰發(fā)量預報的傅立葉級數(shù)模型 ［J］．武漢大學學報：工學版，2005，38 （6）：45－47，52．

［5］茆 智，李遠華，李會昌．逐日作物需水量預測數(shù)學模型研究 ［J］．武漢水利電力大學學報，1995，28（3）：253－259．

［6］郭宗樓，白憲臺，馬學強．作物需水量灰色預測模型［J］．水電能源科學，1995，13 （3）：186－192．

［7］顧世祥，李遠華，袁宏源．霍泉灌區(qū)作物需水量實時預報 ［J］．武漢水利電力大學學報，1998，31（1）：37－41．

［8］霍再林，史海濱，陳亞新，等．參考作物潛在蒸散量的人工神經(jīng)網(wǎng)絡模型研究 ［J］．農(nóng)業(yè)工程學報，2004，20 （1）：40－43．

［9］崔遠來，馬承新，沈細中，等．基于進化神經(jīng)網(wǎng)絡的參考作物騰發(fā)量預測 ［J］．水科學進展，2005，16（1）：76－80．

［10］徐俊增，彭世彰，張瑞美，等．基于氣象預報的參考作物蒸發(fā)蒸騰量的神經(jīng)網(wǎng)絡預測模型 ［J］．水利學報，2006，37 （3）：376－379．

［11］彭世彰，魏 征，徐俊增，等．參考作物騰發(fā)量主成分神經(jīng)網(wǎng)絡預測模型 ［J］．農(nóng)業(yè)工程學報，2009，24（9）：161－164．

［12］段春青，邱 林，黃 強，等．基于混沌遺傳程序設計的參考作物騰發(fā)量預測模型 ［J］．水利學報，2006，37 （4）：499－503．

［13］王惠文，吳載斌，孟 潔．偏最小二乘回歸的線性與非線性方法 ［M］．北京：國防工業(yè)出版社，2006．

［14］王 斌，王貴作，黃金柏，等．農(nóng)業(yè)旱情評估模型及其應用 ［M］．北京：中國水利水電出版社，2011．

［15］熊運章，宋松柏，胡彥華，等．計算機在農(nóng)業(yè)水土工程中的應用 ［M］．北京：清華大學出版社，1999．

［16］中國氣象局．地面氣象觀測規(guī)范 ［M］．北京：氣象出版社，2003．