苗正偉 徐利崗
(1.河北水利電力學(xué)院,河北 滄州 061001; 2.寧夏水利科學(xué)研究院,寧夏 銀川 750021)
降水的科學(xué)預(yù)測,不僅可為制定水資源開發(fā)利用策略提供依據(jù),又可為防災(zāi)、抗災(zāi)、減災(zāi)提供參考。由物理成因的定性分析及大量降水序列資料的統(tǒng)計分析得知,降水量可視為一列相依的隨機變量,其相依程度的強弱,可采用自相關(guān)系數(shù)作為其定量的測度[1],基于此,可建立馬爾可夫鏈模型進行預(yù)測分析。但實際上,降水是一種非常復(fù)雜的自然現(xiàn)象,具有大量的不確定性因素,很難嚴格滿足“齊次性”[2];同時,有研究表明:考慮隨機變量的中間過程狀態(tài)從而充分利用信息,可提高預(yù)測精度[3],因此本文以降水量序列的規(guī)范化后的各階自相關(guān)系數(shù)為權(quán),建立加權(quán)馬爾可夫模型進行預(yù)測。
其原理參見文獻[4],基本步驟如下:
1)建立年降水量的分級標準[5],從而確定降水序列的狀態(tài)空間??梢詷颖镜木怠獦藴什顬闃藴?,將年降水序列劃分為五級,此時年降水序列的狀態(tài)空間E={1,2,3,4,5}。
2)按照上述分級標準確定資料序列中各年降水量的狀態(tài),構(gòu)成年降水狀態(tài)序列。
3)計算一步轉(zhuǎn)移概率:對于降水狀態(tài)這一隨機變量,很難確定其準確的轉(zhuǎn)移概率,但樣本序列足夠長時,可用轉(zhuǎn)移頻率近似代替轉(zhuǎn)移概率,其計算公式如下:
(1)
其中,m為樣本序列的狀態(tài)個數(shù);fij為樣本序列中年降水量從狀態(tài)i一步轉(zhuǎn)移到j(luò)的頻次;pij為由樣本序列統(tǒng)計得出的狀態(tài)i到j(luò)的一步轉(zhuǎn)移概率(頻率)。
4)對樣本序列做馬氏性檢驗[6]:當(dāng)樣本序列長度較大時,構(gòu)造服從自由度為(m-1)2的χ2分布的統(tǒng)計量:
(2)
其中,p*j為邊際頻率,計算公式如下:
(3)
給定顯著性水平α,查表可得分位點χ2((m-1)2)的值,若χ2>χ2((m-1)2),則認為樣本序列具有馬氏性,可作為馬爾可夫鏈來處理。
5)計算各階自相關(guān)系數(shù):
(4)
6)對各階自相關(guān)系數(shù)規(guī)范化:
(5)
其中,ωk為步長為k的轉(zhuǎn)移概率的權(quán)重;b為按預(yù)測需要計算到的最大步長。
7)計算轉(zhuǎn)移矩陣:與式(4)同理,由樣本序列統(tǒng)計各狀態(tài)間各步長的轉(zhuǎn)移概率,各步長的轉(zhuǎn)移概率構(gòu)成相應(yīng)步長的轉(zhuǎn)移矩陣,它決定了降水狀態(tài)轉(zhuǎn)移的概率法則。
關(guān)于擴大灌溉面積。第一,還是要堅持外延與內(nèi)涵相結(jié)合。內(nèi)涵就是對老的配套不完善的這些灌區(qū)進行挖潛改造,特別是著力解決好“最后一公里”的問題,這樣可以使有效灌溉面積得到恢復(fù)和改善。第二,在水土資源條件組合比較好、有開發(fā)利用潛力的地方,結(jié)合千億斤糧食生產(chǎn)規(guī)劃,大力發(fā)展灌溉面積。第三,大力推行節(jié)水措施,包括北方地區(qū)的節(jié)水、東北四省區(qū)的節(jié)水、西北地區(qū)的節(jié)水。第四,充分地調(diào)動農(nóng)民的積極性,調(diào)動社會的資金。比如像剛才郝益東同志提出,內(nèi)蒙古搞的大型噴灌,基本上很多都是用社會資金來搞的。還有水權(quán)置換,工業(yè)企業(yè)要上新項目,但是沒有用水指標,那么就支持河套灌區(qū)的農(nóng)業(yè)節(jié)水,通過節(jié)出來的水量發(fā)展工業(yè)項目。
9)將同一狀態(tài)的各預(yù)測概率加權(quán)后求和:
(6)
則,max{pi,i∈E}所對應(yīng)的i即為該年年降水的預(yù)測狀態(tài)。若該年年降水狀態(tài)可確定,將其加入原序列,重復(fù)以上步驟進行下一年度的預(yù)測。
10)進一步對該馬爾可夫鏈的特征(非周期、遍歷性、平穩(wěn)分布等)進行分析。
于“中國氣象數(shù)據(jù)網(wǎng)”下載了榆林地區(qū)1951年—2010年的逐月降水?dāng)?shù)據(jù),并將之整理為年降水序列。首先以1951年—2008年的資料建立模型預(yù)測2009年的年降水狀態(tài)。
表1 榆林地區(qū)1951年—2008年降水量分級表
表2 榆林地區(qū)1951年—2008年降水量及其狀態(tài)
3)一步轉(zhuǎn)移概率及馬氏性檢驗。經(jīng)統(tǒng)計計算:
由表1知榆林地區(qū)年降水序列的狀態(tài)數(shù)m=5,給定顯著性水平α=0.05,查χ2分布分位數(shù)表可得χ2((5-1)2)=χ2(16)=26.296??梢姡?>χ2(16),因此,該時間序列具有馬氏性,即榆林地區(qū)1951年—2008年年降水狀態(tài)序列可視為馬爾可夫鏈。
4)計算各階自相關(guān)系數(shù)和權(quán)重,結(jié)果如表3所示。
表3 各階自相關(guān)系數(shù)及權(quán)重
5)對降水狀態(tài)序列進行統(tǒng)計計算,得步長從1~5的轉(zhuǎn)移概率矩陣:
6)根據(jù)2004年—2008年的年降水狀態(tài)預(yù)測2009年的年降水狀態(tài),結(jié)果見表4。
表4 2009年年降水狀態(tài)預(yù)測
由表4可知max{pi,i∈E}=0.449,對應(yīng)的狀態(tài)為3,即榆林地區(qū)2009年年降水量狀態(tài)為3,平水年。榆林地區(qū)2009年的實際降水量為420.8 mm,由表1知,正屬于狀態(tài)3。預(yù)測結(jié)果與實際結(jié)果吻合。
7)以1951年—2009年的實測資料重復(fù)上述步驟預(yù)測2010年的降水狀態(tài),結(jié)果見表5。
表5 2010年年降水狀態(tài)預(yù)測
由表5知,榆林地區(qū)2010年年降水量預(yù)測狀態(tài)為3,其實際降水量為363.9 mm,屬于平水年,預(yù)測與實際相吻合。
1)按照均值—標準差方法將榆林地區(qū)實測年降水序列劃分為豐水、偏豐、平水、偏枯、枯水年五種狀態(tài),利用加權(quán)馬爾可夫鏈分別對2009年,2010年降水狀態(tài)進行預(yù)測,結(jié)果均與實測情況吻合,初步說明用加權(quán)馬爾可夫鏈預(yù)測榆林地區(qū)的年降水狀態(tài)是可行的。
2)狀態(tài)標準的劃分有一定的隨意性,對結(jié)果會產(chǎn)生一定影響,如何更科學(xué)的劃分降水狀態(tài)有待于進一步探討。
3)馬爾可夫鏈的預(yù)測結(jié)果是一個區(qū)間,如何結(jié)合其他理論或方法進一步預(yù)測出較為精確的具體數(shù)值,從而拓展預(yù)測結(jié)果的適用范圍,值得進一步研究。
4)用馬爾可夫鏈進行預(yù)測的基本假設(shè)是:年降水狀態(tài)序列是齊次馬爾可夫鏈,但實際上這很難嚴格滿足:齊次馬氏鏈的統(tǒng)計特性完全取決于初始分布和一步轉(zhuǎn)移概率,但對于年降水狀態(tài)序列卻很難獲得準確的轉(zhuǎn)移概率,只能用由足夠長的實際資料進行統(tǒng)計得出的轉(zhuǎn)移頻率代替轉(zhuǎn)移概率;齊次馬爾可夫鏈的n步轉(zhuǎn)移概率矩陣等于一步轉(zhuǎn)移概率矩陣的n次方,這對于由統(tǒng)計計算得出的各步轉(zhuǎn)移矩陣顯然是難以成立的。這些都會對預(yù)測精度產(chǎn)生影響,盡管引入了加權(quán)馬氏鏈,但這不能從根本上消弭不足。因此該方法對于水文時間序列預(yù)測的適用程度有待進一步檢驗。
參考文獻:
[1] 韓璞璞,張 生,李暢游,等.基于權(quán)馬爾可夫鏈模型的廬江縣降水量預(yù)測[J].水文,2012,32(3):38-42.
[2] 韓合忠,高淑會,國偉華.基于加權(quán)馬爾可夫模型的濟南市降水豐枯狀況預(yù)測研究[J].水文,2012,32(5):72-75.
[3] 王鑫東.加權(quán)馬爾科夫鏈模型在農(nóng)業(yè)灌溉用水預(yù)測中的應(yīng)用研究[J].中國農(nóng)村水利水電,2016(5):58-64.
[4] 劉次華.隨機過程[M].第5版.武漢:華中科技大學(xué)出版社,2014.
[5] 王宇博,梁秀娟,喬 雨,等.基于疊加馬爾科夫鏈和BP神經(jīng)網(wǎng)絡(luò)模型的降水量預(yù)測研究[J].中國農(nóng)村水利水電,2014(9):80-86.
[6] 王 濤,錢 會,李培月.加權(quán)馬爾可夫鏈在銀川地區(qū)降雨量預(yù)測中的應(yīng)用[J].南水北調(diào)與水利科技,2010,8(1):78-81.