基于SARIMA模型的贛江流域近50年參考作物需水量預(yù)測(cè)

曾招財(cái)，梁藉，吳豪，曾志強(qiáng)

(華中科技大學(xué)水電與數(shù)字化工程學(xué)院，湖北武漢430074)

中國(guó)是一個(gè)農(nóng)業(yè)大國(guó)，保證糧食安全的一大前提就是要使農(nóng)業(yè)用水得到充分保障。隨著中國(guó)農(nóng)業(yè)現(xiàn)代化的不斷發(fā)展，“農(nóng)業(yè)經(jīng)濟(jì)用水量[1]”這個(gè)概念被相關(guān)研究人員廣泛地提及。在中國(guó)的農(nóng)業(yè)用水量中，很重要的組成部分就是農(nóng)作物需水量，在整個(gè)國(guó)民經(jīng)濟(jì)的總耗水量中占據(jù)著非常大的比重。參考作物需水量(以下簡(jiǎn)稱ET0)用以確定具體地區(qū)的農(nóng)作物需水量，從而得出區(qū)域的農(nóng)業(yè)灌溉用水量，為相關(guān)部門制定作物灌溉制度提供可靠的現(xiàn)實(shí)基礎(chǔ)，同時(shí)也是流域布局規(guī)劃、地區(qū)水利工程及灌溉排水工程的規(guī)劃、設(shè)計(jì)和運(yùn)營(yíng)管理的基本依據(jù)[2]。作物需水量的計(jì)算可以分為直接測(cè)定法和間接推求法。間接推求作物需水量的方法最早由Penman在1948年提出。通過(guò)引入表面阻力的相關(guān)概念，Monteith于1963年在Penman等人對(duì)于作物騰發(fā)量研究的相關(guān)理論的基礎(chǔ)上得出了后來(lái)被FAO(國(guó)際糧農(nóng)組織)確定為參考作物需水量的標(biāo)準(zhǔn)方法的彭曼-蒙特斯方法[3]。Doorenboos和Pruit于1977年在Makkink公式中引入調(diào)整系數(shù)來(lái)解決空氣動(dòng)力學(xué)因素的影響問(wèn)題，提出了適用于干旱地區(qū)的FAO-Radiation法[4]。時(shí)間序列分析法最開始由英國(guó)統(tǒng)計(jì)學(xué)家G.U.Yule于1972年提出，并成功預(yù)測(cè)了市場(chǎng)的變化趨勢(shì)和規(guī)律[5]。經(jīng)過(guò)不斷的應(yīng)用和發(fā)展，時(shí)間序列分析的相關(guān)研究方法已經(jīng)廣泛應(yīng)用于于經(jīng)濟(jì)學(xué)、工程學(xué)、氣象學(xué)、水文學(xué)、醫(yī)學(xué)、信息學(xué)、環(huán)境科學(xué)、社會(huì)科學(xué)、數(shù)據(jù)挖掘等各個(gè)領(lǐng)域[6]。

利用時(shí)間序列分析方法對(duì)江西省境內(nèi)贛江流域需水量變化特征及其預(yù)測(cè)的相關(guān)研究，可以充分利用歷史氣象資料，研究流域內(nèi)參考作物的需水特性及規(guī)律，有利于當(dāng)?shù)毓?jié)水灌溉的發(fā)展。在建設(shè)節(jié)約型社會(huì)的今天，水資源高效利用一直是一個(gè)熱點(diǎn)話題，掌握區(qū)域作物需水規(guī)律正是為了更好地利用水資源，提高水資源的利用效率。論文從分析贛江流域參考作物需水量年際、年內(nèi)變化特征出發(fā)，利用SPSS統(tǒng)計(jì)分析軟件對(duì)ET0時(shí)間序列建立季節(jié)性ARIMA模型，對(duì)ET0的變化進(jìn)行了擬合及預(yù)測(cè)。

1 研究區(qū)域概況

作為長(zhǎng)江的主要支流之一，贛江也是江西省最大的河流。贛江位于長(zhǎng)江中下游南岸，其源頭是江西省和福建省交接處的武夷山脈西邊，自南向北縱貫全省，從九江湖口匯入長(zhǎng)江。贛江共有13條主要支流匯入，總長(zhǎng)766 km，流域面積83 500 km2。從河源至贛州市為上游，稱貢水，在贛州市城西和章水匯流后開始被稱為贛江。貢水為贛江東源全長(zhǎng)255 km，發(fā)源于江西省贛州市石城縣橫江鎮(zhèn)。章水是贛江西源，又名章江，發(fā)源于贛州市崇義縣的聶都山。贛州至新干為中游，長(zhǎng)303 km，穿行于丘陵之間。新干至吳城為下游，長(zhǎng)208 km。

2 數(shù)據(jù)來(lái)源與研究方法

2.1 數(shù)據(jù)來(lái)源

論文原始數(shù)據(jù)為選自中國(guó)氣象數(shù)據(jù)網(wǎng)的12個(gè)氣象站1961—2010年50 a的逐日實(shí)測(cè)氣象數(shù)據(jù)(氣溫、水汽壓、濕度、日照時(shí)數(shù)、風(fēng)速、太陽(yáng)輻射)。采用FAO-Penman-Monteith公式計(jì)算得到ET0的逐日數(shù)據(jù)。根據(jù)流域內(nèi)站點(diǎn)的分布特點(diǎn)，按照站點(diǎn)所在的地理位置，將其分為上游站點(diǎn)(龍南、贛縣)、中游站點(diǎn)(井岡山、遂川、吉安、寧都、永豐、蓮花)和下游站點(diǎn)(宜春、宜豐、樟樹、南昌)3個(gè)子研究區(qū)域。選取的12個(gè)氣象站點(diǎn)(其中11個(gè)選取自贛江流域內(nèi)，1個(gè)位于流域外)的空間分布見(jiàn)圖1。

圖1 贛江流域氣象站空間分布

2.2 研究方法

參考作物是一種假想的草類參照作物，其高度為0.12 m,具有固定的表面阻力70 s/m和反射率0.23。根據(jù)《參考作物需水量計(jì)算指南》書中的描述，此假想?yún)⒄兆魑锩骖愃埔粋€(gè)面積很大、高度均勻、生長(zhǎng)旺盛、完全遮蔽地面且供水充分的綠色草地。國(guó)際糧農(nóng)組織推薦在濕潤(rùn)地區(qū)采用Penman-Monteith[7]公式來(lái)計(jì)算參考作物需水量(也稱參考作物騰發(fā)量)，該方法已經(jīng)得到了廣泛的應(yīng)用，其表達(dá)式如下：

ET0=

(1)

式中ET0——參照騰發(fā)速率，mm/d；Rn——凈太陽(yáng)輻射，MJ/(m2·d)；G——土壤熱通量，MJ/(m2·d)；T——2 m高度處的日平均氣溫,℃；u2——2 m高度處的風(fēng)速,m/s；es——飽和水汽壓，kPa；eq——實(shí)際水汽壓，kPa；γ——濕度計(jì)常數(shù)，kPa/℃；Δ——飽和水汽壓曲線斜率。

對(duì)于季節(jié)性ARIMA(p,d,q)(P,D,Q)S模型[8]，可以用下式來(lái)表示:

(2)

3 結(jié)果與討論

3.1 ET0年際變化分析

將計(jì)算得到的各個(gè)站點(diǎn)ET0日數(shù)據(jù)進(jìn)行整理得到了ET0年平均值數(shù)據(jù)并在Matlab軟件中繪制序列，得到圖2。從圖中可以看出全流域內(nèi)ET0的年際變化較大，全流域ET0最大值為1963年的1 228.7 mm/a，最小值1997年的937.5 mm/a，極差為291.2 mm/a。由于氣候變化，1986年到2005年這連續(xù)20 a內(nèi)ET0經(jīng)歷了從1986—始波動(dòng)下降而后波動(dòng)上升的變化過(guò)程。從上游、中游、下游3個(gè)子研究區(qū)域的ET0年際變化來(lái)看，上游的年ET0值最大，三者的變化特征和趨勢(shì)和全流域的ET0變化規(guī)律相似。從50 aET0年平均值來(lái)看，全流域平均值為1 061.0 mm/a，上游為1 116.5 mm/a，中游為1 036.1 mm/a，下游為1 038.5 mm/a，整體變化并不大。

圖2 贛江流域參考作物需水量年際變化

3.2 ET0年內(nèi)變化分析

對(duì)于全流域ET0日均值的年內(nèi)變化情況，圖3給出了其日時(shí)序和月時(shí)序年內(nèi)分布。從圖中可以發(fā)

a) 日時(shí)序圖3 贛江流域參考作物需水量年內(nèi)變化

b)月時(shí)序續(xù)圖3 贛江流域參考作物需水量年內(nèi)變化

現(xiàn)，ET0在一年中的中期即6、7月的時(shí)候達(dá)到最大值然后開始下降。在6、7月份處于夏季溫度較高，降雨也多，導(dǎo)致作物騰發(fā)量較大，所以作物參考需水量也大。過(guò)了夏季，氣溫降低使得作物騰發(fā)量開始顯著減小，所以其分布曲線在6、7月份之后開始呈現(xiàn)下降趨勢(shì)。

3.3 時(shí)間序列特性分析

時(shí)間序列建模要求數(shù)據(jù)必須為平穩(wěn)序列[9]。利用SPSS軟件作出研究序列的時(shí)間序列圖和自相關(guān)系數(shù)圖。從圖4可以看出，全流域ET0月數(shù)據(jù)序列為非平穩(wěn)序列，其數(shù)值圍繞平均值88.42 mm/月上下波動(dòng)。同時(shí)，圖5自相關(guān)圖也表現(xiàn)為非截尾性和非拖尾性，可以認(rèn)為序列非平穩(wěn)。因此，需要對(duì)數(shù)據(jù)進(jìn)行差分運(yùn)算。經(jīng)研究發(fā)現(xiàn)，對(duì)其進(jìn)行一階季節(jié)性差分可以使序列趨于平穩(wěn)，序列差分后序列的時(shí)序見(jiàn)圖6，其自相關(guān)系數(shù)見(jiàn)圖7。

圖4 贛江流域月參考作物需水量時(shí)序

圖5 月ET0序列自相關(guān)系數(shù)

圖6 月ET0序列一階季節(jié)性差分后時(shí)序

圖7 ET0序列一階季節(jié)性差分后自相關(guān)

3.4 模型的建立

3.4.1模型階數(shù)識(shí)別

時(shí)間序列模型階數(shù)的確定是建立模型的重要一步。隨著統(tǒng)計(jì)軟件和計(jì)算機(jī)技術(shù)的不斷發(fā)展，可以運(yùn)用統(tǒng)計(jì)軟件根據(jù)擬合效果的信息準(zhǔn)則法來(lái)確定時(shí)間序列模型的階數(shù)。常見(jiàn)的判別準(zhǔn)則有FPE準(zhǔn)則、AIC準(zhǔn)則和BIC準(zhǔn)則。論文結(jié)合SPSS軟件，采用BIC準(zhǔn)則來(lái)確定模型的階數(shù)。BIC準(zhǔn)則表達(dá)式如下：

(3)

如果在某一階數(shù)p0下滿BIC(p0)=minBIC(p)(1≤p≤L)，式中L是預(yù)先假定的所研究的模型階數(shù)的上限值，則認(rèn)為模型的最佳階數(shù)為p0。利用贛江全流域的參考作物需水量數(shù)據(jù)，經(jīng)過(guò)不斷試驗(yàn)和調(diào)整并結(jié)合BIC準(zhǔn)則法，選定合適的時(shí)間序列模型[10]。論文采用1961—2010年50 a月時(shí)序數(shù)據(jù)進(jìn)行模型的建立。首先以前49 a作為模型的訓(xùn)練期，2010年作為檢驗(yàn)期，利用SPSS軟件得到贛江流域的季節(jié)性ARIMA模型形式為SARIMA(2,0,1)(0,1,1)。

3.4.2模型參數(shù)估計(jì)

ARIMA模型參數(shù)的估計(jì)有比較多的方法,例如矩估計(jì)、極大似然估計(jì)最小二乘估計(jì)等。利用SPSS軟件的分析工具可以得到模型的各項(xiàng)參數(shù)。對(duì)于贛江流域參考作物月需水量時(shí)間序列模型。其模型統(tǒng)計(jì)量和模型的相關(guān)參數(shù)見(jiàn)表1，可以得到贛江流域ET0預(yù)測(cè)模型的表達(dá)式為：

(1+0.8B-0.195B2)(1-B12)xt=

(1+0.85B12)(1-0.999B)wt

表1 模型SARIMA(2,0,1)(0,1,1)參數(shù)

3.4.3模型的診斷檢驗(yàn)

模型的顯著性檢驗(yàn)，即模型的有效性檢驗(yàn)，驗(yàn)證模型是否充分地提取了數(shù)據(jù)的有關(guān)信息。從圖8中，可以看出殘差序列的自相關(guān)系數(shù)非常小，接近于0。這說(shuō)明，模型殘差是白噪聲序列，建立的模型已經(jīng)充分提取了數(shù)據(jù)的有效信息。表2給出了模型的統(tǒng)計(jì)量和LB(Ljung-Box)值。可以知道，當(dāng)LB值較大時(shí)，拒絕原假設(shè)，說(shuō)明模型擬合不顯著，當(dāng)LB值較小時(shí)，說(shuō)明模型擬合顯著顯著。經(jīng)計(jì)算得到LB統(tǒng)計(jì)量的值15.412，其 Sig.值大于顯著性水平0.05，所以不能拒絕原假設(shè)，模型的殘差序列是白噪聲序列，模型擬合顯著有效，可以用于預(yù)測(cè)。

圖8 殘差自相關(guān)

模型擬合統(tǒng)計(jì)量平穩(wěn)R2統(tǒng)計(jì)量Ljung-Box Q(18)R2BICSig.0.420.914.9115.41140.35

3.4.4預(yù)測(cè)結(jié)果分析

一般認(rèn)為時(shí)間序列模型的短期預(yù)測(cè)較為準(zhǔn)確。利用1961—2009年月數(shù)據(jù)來(lái)預(yù)測(cè)2010年各月的參考作物需水量，結(jié)果見(jiàn)表3。結(jié)果顯示最小相對(duì)誤差為2.25%，最大相對(duì)誤差為23.04%，平均相對(duì)誤差為10.02%，在正常允許范圍之內(nèi)。為了避免模型預(yù)測(cè)結(jié)果的偶然性，利用2001—2010年作為檢驗(yàn)期，進(jìn)行滾動(dòng)預(yù)測(cè)，發(fā)現(xiàn)月尺度下ARIMA模型預(yù)測(cè)的平均相對(duì)誤差最小為2006年的4%，最大為2002年的15%，10 a平均相對(duì)誤差為9.1%，年尺度下最小相對(duì)誤差0.41%，最大相對(duì)誤差8.61%，10 a平均相對(duì)誤差4.35%，均在合理范圍之內(nèi)。年尺度預(yù)測(cè)結(jié)果見(jiàn)表4。2001—2010年月尺度實(shí)際值和預(yù)測(cè)值擬合見(jiàn)圖9。從圖9中可以看出，實(shí)際值和預(yù)測(cè)值擬合較好，模型的預(yù)測(cè)精度滿足要求，說(shuō)明考慮季節(jié)性因素的ARIMA模型運(yùn)用于贛江流域ET0預(yù)測(cè)具有一定合理性。

表3 2010年月參考作物需水量預(yù)測(cè)結(jié)果

表4 2001—2010年參考作物需水量預(yù)測(cè)結(jié)果

圖9 實(shí)際值和預(yù)測(cè)值對(duì)比

4 結(jié)論

a) 贛江流域ET0在1961—2010年50 a間經(jīng)歷了先波動(dòng)遞減后波動(dòng)上升的過(guò)程。且3個(gè)子研究區(qū)域的ET0年際變化和全流域的變化規(guī)律大體一致，說(shuō)明贛江流域ET0在整個(gè)流域內(nèi)分布較為均勻。ET0年內(nèi)變化則呈現(xiàn)明顯的季節(jié)性特征。贛江流域參考作物需水量時(shí)間序列為非平穩(wěn)序列，在考慮季節(jié)性因素的情況下進(jìn)行一階差分可使序列平穩(wěn)。

b) 預(yù)測(cè)結(jié)果表明，月尺度下ARIMA模型預(yù)測(cè)的平均相對(duì)誤差最小為2006年的4%，最大為2002年的15%，10 a平均相對(duì)誤差為9.1%；年尺度下最小相對(duì)誤差0.41%，最大相對(duì)誤差8.61%，10 a平均相對(duì)誤差4.35%，均在合理范圍之內(nèi)。說(shuō)明ARIMA模型在考慮季節(jié)性因素的情況下運(yùn)用于贛江流域參考作物需水量預(yù)測(cè)，具有實(shí)際意義，可以為農(nóng)業(yè)生產(chǎn)和相關(guān)研究提供參考。