李 揚 秦 聰 牛永華
(太原理工大學水利科學與工程學院,山西 太原 030024)
降水是水文循環(huán)的最基本環(huán)節(jié),可利用降水量的多少直接影響農(nóng)業(yè)生產(chǎn)與人民生活。山西省屬顯著的大陸性季風氣候,大氣降水是水資源的主要補給來源,全省多年平均年降水量為400~600mm,主要集中在7月、8月、9月,且呈東南向西北遞減態(tài)勢。山西省降水時空分布不均,“十年九旱,旱澇交錯”是其主要特點,水資源供需矛盾突出。
汾河灌區(qū)是山西省最大的自流灌區(qū),設計灌溉面積為149.55萬畝,占山西省總灌溉面積的10%。對灌區(qū)而言,灌溉用水量的影響因素既有氣候條件、灌溉方式、種植結(jié)構(gòu)等自然因素,也有水價、管理水平等人為因素。而其中最為關(guān)鍵的自然因素是氣候條件,即降水量及其年內(nèi)分配過程。相關(guān)研究表明,越干旱的地區(qū),其降水量變化對農(nóng)業(yè)灌溉用水量的影響越大[1]。因此,準確預測灌區(qū)降水量,對農(nóng)業(yè)灌溉用水量的估計以及農(nóng)業(yè)水價的制定具有重要意義。
近年來,已有學者采用線性回歸分析法、M-K檢驗法、小波分析法等多種方法對山西省不同地區(qū)的降水量進行分析和預測[2]。如:鈔錦龍等[3]采用Mann-Kendall突變檢驗法、經(jīng)驗正交函數(shù)(EOF)以及馬爾可夫模型等對山西北部地區(qū)的降水時空分布格局進行了研究和預測。解建強等[4]采用距平分析法、線性傾向估計法以及Mann-Kendall檢驗法,對山西省右玉地區(qū)的降水量進行了系統(tǒng)分析,結(jié)果表明,右玉縣降水量呈增加趨勢。針對灌區(qū)的研究尚不多見,因此,本文采用小波分析方法,以汾河灌區(qū)介休灌溉區(qū)域為研究對象,分析研究區(qū)1954—2022年降水量在多時間尺度上的變化特征,并構(gòu)建ARMA模型對介休灌溉區(qū)域2023年的降水量進行預測,以期為汾河灌區(qū)農(nóng)業(yè)水資源利用以及農(nóng)業(yè)水價的制定提供參考依據(jù)。
汾河灌區(qū)位于山西省中部的太原盆地,其中介休灌溉區(qū)域位于灌區(qū)中南部,地跨東經(jīng)111°44′10″~112°10′14″、北緯36°50′1″~37°11′4″,總面積為744km2,屬暖溫帶季風氣候,四季分明,春季多風干燥,夏季多雨、炎熱,秋季多晴,冬季少雪、寒冷。全市降水一般集中在7—9月,年無霜期平均為175天。
介休灌溉區(qū)域1954—2022年月降水量資料來源于國家氣象科學數(shù)據(jù)中心(http://data.cma.cn/),采用MATLAB軟件對降水量變化趨勢進行小波分析,并采用ARMA模型對降水量進行預測。
1.3.1 小波變換
小波變換是分析降水時間序列的典型方法,由于降水存在趨勢性、周期性、隨機性、突變性等特征,通過信號的伸縮平移,可獲得不同時間尺度的小波函數(shù),提取出不同時間尺度的降水分布特征,CWT(連續(xù)小波變換)也被廣泛應用于其中。母小波的選擇是連續(xù)小波變換的基礎,Morlet小波不僅包含多種振動信息,而且可以平衡頻率與時間的局部化問題[5-7]??紤]到降水時間序列的平滑連續(xù)性以及所需小波函數(shù)的非正交性,選取Morlet小波作為母小波。
ψ0(t)=π-1/4eiw0te-t2/2
(1)
式中:w0為無量綱頻率,t為時間。
1.3.2 ARMA模型
ARMA模型(自回歸滑動平均模型)是研究平穩(wěn)隨機序列的典型方法之一[8],主要步驟為序列平穩(wěn)性檢驗、參數(shù)確定、模型診斷與預測,其表達式為
Yt=β1Yt-1+β2Yt-2+…+βpYt-p+εt+α1εt-1+
α2εt-2+…+αqεt-q
(2)
式中:Yt為預測序列;Yt-i(i=1,2,…,p)為滯后序列;εt-i(i=1,2,…,q)為隨機序列;εt為白噪聲序列;βi(i=1,2,…,p)為自回歸系數(shù);αj(j=1,2,…,q)為移動平均系數(shù)。
從介休灌溉區(qū)域降水量年際變化圖(見圖1)可看出,年降水量曲線大致在420mm上下波動,年降水量的最大值和最小值分別出現(xiàn)在2021年(714.4mm)和1957年(248.6mm),極值比為2.87。根據(jù)1954—2022年降水量的年內(nèi)分配(見圖2)可知,介休灌溉區(qū)域6—9月降水較豐,占到全年的70%以上,其中7—8月為降水集中期,約占全年降水的46%,而1月、2月和12月為枯季,降水量占比僅為全年的3%左右。
圖1 介休灌溉區(qū)域降水量年際變化過程
圖2 介休灌溉區(qū)域1954—2022年降水量年內(nèi)分配情況
為消除數(shù)據(jù)不同特征的尺度差異對分析結(jié)果的影響,在進行小波分析之前,先采用SPSS軟件對降水量數(shù)據(jù)進行Z值標準化處理,使其成為無量綱時序數(shù)據(jù)。之后,根據(jù)降水量標準化后的序列,使用MATLAB軟件對介休灌溉區(qū)域1954—2022年降水量數(shù)據(jù)進行小波分析,得到研究區(qū)降水序列的Morlet小波系數(shù)實部的時頻變化圖(見圖3)。
圖3 介休灌溉區(qū)域降水序列的Morlet小波系數(shù)實部時頻變化
圖3中,橫、縱坐標分別表示年份和時間尺度,實部值的大小表示降水量的多少,實線表示實部值為正,降水增多;虛線表示實部值為負,降水減少。由圖3可知,較長時間尺度的振蕩是低頻振蕩,位于圖形的上半部分,等值線較為稀疏,較短時間尺度的振蕩是高頻振蕩,位于圖形的下半部分,等值線相對密集。小波系數(shù)實部虛實線交替出現(xiàn),表明介休灌溉區(qū)域降水量變化過程存在周期性,主要出現(xiàn)在以5年、10年、15年、25年為中心的時間尺度上。
在圖4小波系數(shù)實部時頻分布圖上,取主要時間尺度(5年、10年、15年、20年)的小波實部值,分析Morlet小波系數(shù)實部在主要時間尺度下隨年份的變化過程。在5年時間尺度下,小波系數(shù)實部振幅在1990年之前呈增大趨勢,在1990—2000年之間振幅減小,2000年以后又開始增大;在10年時間尺度下,小波系數(shù)實部振幅整體呈現(xiàn)減小趨勢,2010年以后,其振幅有增大趨勢;在15年時間尺度下,小波系數(shù)實部振幅在1990年前呈減小趨勢,在1990年后呈增大趨勢;在25年時間尺度下,小波系數(shù)實部振幅近69年來一直呈現(xiàn)緩慢增大趨勢。根據(jù)介休灌溉區(qū)域降水量小波實部在不同時間尺度下的變化過程圖可以推測,2023年介休灌溉區(qū)域降水量將呈現(xiàn)增加趨勢。
圖4 介休灌溉區(qū)域降水量小波實部在不同時間尺度下的變化過程
從以上分析中可以看出,不同時間尺度下,介休灌溉區(qū)域降水量呈現(xiàn)不同的變化趨勢,時間尺度越小,其降水量的變化越劇烈,當時間尺度增大時,則呈現(xiàn)相對穩(wěn)定的變化趨勢,因此,降水量的變化趨勢與時間尺度存在密切的關(guān)系。
借助EViews軟件,用ARMA模型對介休灌溉區(qū)域1954—2022年降水量進行分析,得出降水量的變化規(guī)律,進而對介休灌溉區(qū)域2023年的年降水量作出預測,以期為介休灌溉區(qū)域氣候變化的預測和農(nóng)業(yè)水價的制定提供參考依據(jù)。
2.3.1 時間序列平穩(wěn)化
通過對介休灌溉區(qū)域1954—2022年的降水量數(shù)據(jù)進行時間序列平穩(wěn)性檢驗和白噪聲檢驗(見圖5)得知,ADF檢驗的P值為0,且自相關(guān)圖(見圖6)中Q統(tǒng)計量對應的P>0.05,表明原降水量序列為平穩(wěn)白噪聲序列,無分析價值。將原序列進行一階差分,ADF檢驗(見圖7)的P值為0,自相關(guān)圖(見圖9)中Q統(tǒng)計量對應的P<0.05,即一階差分序列為平穩(wěn)非白噪聲序列,且降水量一階差分后的序列(見圖8)在0附近上下波動。
圖5 介休灌溉區(qū)域1954—2022年的降水量數(shù)據(jù)ADF檢驗
圖6 介休灌溉區(qū)域1954—2022年的降水量的自相關(guān)圖和偏自相關(guān)圖
圖7 一階差分序列的ADF檢驗
圖8 介休灌溉區(qū)域1954—2022年降水量一階差分序列
圖9 一階差分序列的自相關(guān)圖和偏自相關(guān)圖
2.3.2 模型參數(shù)確定
介休灌溉區(qū)域降水量時間序列在一階差分后為平穩(wěn)非白噪聲序列,故采用ARMA(p,q)進行降水量預測,其自相關(guān)圖和偏自相關(guān)圖如圖9所示。
圖9中,左圖為自相關(guān)圖(AC圖),右圖為偏自相關(guān)圖(PAC圖),兩虛線之間為隨機區(qū)間,模型存在4種組合:ARMA(1,0),ARMA(1,1),ARMA(2,0),ARMA(2,1)。使用EViews軟件對以上4種模型進行回歸操作,逐步剔除不顯著的滯后項或移動平均項后,發(fā)現(xiàn)只有ARMA(1,0)和ARMA(2,0)符合顯著性條件,需從ARMA(1,0)和ARMA(2,0)兩模型中選出較優(yōu)的一個作為預測模型。
分別計算ARMA(1,0)模型和ARMA(2,0)模型的AIC、SC、HQ值。由圖10可知,ARMA(1,0)模型的AIC值為12.66147,SC值為12.69411,HQ值為12.67440。由圖11可知,ARMA(2,0)模型的AIC值為12.56654,SC值為12.63182,HQ值為12.59241,均小于ARMA(1,0)模型的相應值。因此,根據(jù)AIC、SC、HQ最小為最優(yōu)的信息準則,最終選擇ARMA(2,0)模型作為研究區(qū)降水量模擬預測模型。
圖10 ARMA(1,0)模型擬合結(jié)果
圖11 ARMA(2,0)模型擬合結(jié)果
2.3.3 白噪聲檢驗
構(gòu)建模型時,假設模型的誤差是序列無關(guān)的,因此,模型構(gòu)建完成后有必要進行誤差的序列相關(guān)性檢驗,以證明模型構(gòu)建良好。如果序列之間不存在序列相關(guān)性,則稱之為白噪聲,即認為該序列已經(jīng)沒有任何可供提取的信息。采用Q統(tǒng)計量的方法進行白噪聲檢驗,結(jié)果見圖12。
圖12 模型殘差圖
由模型殘差圖(見圖12)可看出,AC值和PAC值所對應的P值基本大于顯著性水平0.01,通過了白噪聲檢驗,表明所構(gòu)建的模型有效。
2.3.4 預測結(jié)果
選用ARMA(2,0)模型Static Forecast方法對介休灌溉區(qū)域2023年降水量進行預測,預測結(jié)果見圖13:實線代表降水量預測值,兩條虛線代表2倍標準差置信區(qū)間,Theil不相等系數(shù)為0.1357,這表明模型的預測能力較為理想。因此,利用ARMA(2,0)模型Static Forecast預測方法,可以較好地展現(xiàn)未來降水量的變化趨勢,從而預測出短期降水量。將實測降水量與模擬降水量進行擬合,由圖14可知,ARMA(2,0)模型的模擬結(jié)果基本能夠反映介休灌溉區(qū)域1954—2022年間降水量的波動特征和演變規(guī)律,得到介休灌溉區(qū)域2023年的降水量預測值為641.21mm。
圖14 ARMA(2,0)模型模擬降水量擬合圖
本文以汾河灌區(qū)介休灌溉區(qū)域為研究對象,采用小波分析方法分析了研究區(qū)1954—2022年降水量在多時間尺度上的變化特征,并構(gòu)建ARMA模型對介休灌溉區(qū)域2023年的降水量進行了預測,取得如下結(jié)論:
a.介休灌溉區(qū)域1954—2022年降水序列具有5年、10年、15年、25年的主周期,且降水量的變化趨勢與時間尺度存在密切的關(guān)系。時間尺度越小,其降水量的變化越劇烈,當時間尺度增大時,變化趨勢相對穩(wěn)定。根據(jù)介休灌溉區(qū)域降水量小波實部在不同時間尺度下的變化過程,預測2023年介休灌溉區(qū)域降水量將呈現(xiàn)增加趨勢。
b.ARMA(2,0)模型的模擬結(jié)果基本能夠反映介休灌溉區(qū)域1954—2022年間降水量的波動特征和演變規(guī)律,采用Static Forecast方法求得的介休灌溉區(qū)域2023年降水量預測值為641.21mm。