逐日參考作物蒸發(fā)蒸騰量預(yù)測的偏最小二乘回歸模型

王士蘭，王 斌

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

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

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

1 建模準(zhǔn)備

1．1 ET0計(jì)算變量分析

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

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

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

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

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

3）P－M公式的8項(xiàng)基本計(jì)算變量 （T、es、Δ、γ、G、ea、Rn、u2）中有6項(xiàng) （T、es、Δ、γ、G、ea）與溫度有關(guān)。因此，對于固定地點(diǎn)，可以考慮選取ET0的部分計(jì)算變量構(gòu)建預(yù)測ET0的回歸模型［14］。

1．2 天氣預(yù)報(bào)提供的氣象信息及其換算

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

表1 天氣狀況與天氣指標(biāo)對應(yīng)表［15］Table 1 Weather conditions and weather indiators corresponding table［15］

表2 風(fēng)力等級與風(fēng)速對應(yīng)表［16］Table 2 Wind level wind speed corresponding table［16］

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

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

表3 天氣狀況與αm值對應(yīng)表Table 3 Weather conditions andαmvalue corresponding table

表4 風(fēng)力等級與um值對應(yīng)表Table 4 Wind level with the umvalue corresponding table

1．3 選擇建模變量

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

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

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

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

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

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

2 偏最小二乘回歸方法在逐日ET0預(yù)測中的應(yīng)用

2．1 建模自變量間的相關(guān)性分析

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

2．2 基于PLS的逐日ET0預(yù)測模型的建立

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

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

2．3 模型擬合及預(yù)測效果檢驗(yàn)

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

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

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

3 結(jié)論與討論

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

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

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

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

［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］羅玉峰，崔遠(yuǎn)來，蔡學(xué)良．參考作物騰發(fā)量預(yù)報(bào)的傅立葉級數(shù)模型 ［J］．武漢大學(xué)學(xué)報(bào)：工學(xué)版，2005，38 （6）：45－47，52．

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

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

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

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

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

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

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

［12］段春青，邱 林，黃 強(qiáng)，等．基于混沌遺傳程序設(shè)計(jì)的參考作物騰發(fā)量預(yù)測模型 ［J］．水利學(xué)報(bào)，2006，37 （4）：499－503．

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

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

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

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