宋 帆，楊曉華，武翡翡，劉 童

(北京師范大學(xué)環(huán)境學(xué)院 水環(huán)境模擬國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100875)

降雨量是流域或地區(qū)水資源循環(huán)系統(tǒng)的重要輸入項(xiàng)目，是最主要的水資源補(bǔ)給途徑，在農(nóng)業(yè)生產(chǎn)、工業(yè)建設(shè)、水利防洪等方面起著重要的作用。由于氣象條件具有多變性以及復(fù)雜性，導(dǎo)致降雨過(guò)程存在較大的隨機(jī)性[1]，如何較為準(zhǔn)確地對(duì)降雨量進(jìn)行預(yù)測(cè)也成了一個(gè)熱點(diǎn)研究問(wèn)題。目前較為常用的降雨量預(yù)測(cè)方法有灰色預(yù)測(cè)[2]、馬爾科夫鏈[3]、人工神經(jīng)網(wǎng)絡(luò)[4]等。

馬爾可夫過(guò)程是研究事物狀態(tài)及狀態(tài)轉(zhuǎn)移規(guī)律的理論。它通過(guò)確定系統(tǒng)不同狀態(tài)的初始概率及狀態(tài)之間的轉(zhuǎn)移概率關(guān)系, 來(lái)確定狀態(tài)的變化趨勢(shì), 從而進(jìn)行預(yù)測(cè)[5]。常用的馬爾科夫鏈有基于絕對(duì)分布的馬爾科夫鏈、疊加馬爾科夫鏈、加權(quán)馬爾科夫鏈等，其中加權(quán)馬爾科夫鏈?zhǔn)切Ч詈?、?yīng)用最廣泛的[6]。本次研究在加權(quán)馬爾科夫鏈的基礎(chǔ)上進(jìn)行了改進(jìn)：首先利用歐氏距離聚類(lèi)代替常用的均值分類(lèi)；其次引入隸屬度函數(shù)來(lái)確定初始狀態(tài)向量；最后對(duì)預(yù)測(cè)結(jié)果進(jìn)行分類(lèi)改進(jìn)。從而建立了基于聚類(lèi)分析的模糊馬爾科夫鏈降雨量預(yù)測(cè)模型，并利用中國(guó)氣象站記錄的1951-2013年降雨量數(shù)據(jù)對(duì)模型進(jìn)行驗(yàn)證。

1 模型原理及步驟

1.1 馬爾科夫預(yù)測(cè)原理

馬爾科夫預(yù)測(cè)是根據(jù)馬爾科夫理論與方法來(lái)研究未來(lái)發(fā)展趨勢(shì)的一種動(dòng)態(tài)分析方法，既適用于時(shí)間序列，又適用于空間序列。其基本特征是無(wú)后效性，即系統(tǒng)在未來(lái)某一時(shí)刻的狀態(tài)只與現(xiàn)在的狀態(tài)有關(guān)，而與過(guò)去的狀態(tài)無(wú)關(guān)。因此馬爾科夫鏈在模擬預(yù)測(cè)隨機(jī)波動(dòng)性較大的問(wèn)題上有較好的效果[7]。

1.2 聚類(lèi)-模糊馬爾可夫鏈模型流程

本次研究所采用的聚類(lèi)-模糊馬爾可夫鏈模型如圖1所示，首先通過(guò)聚類(lèi)分析將降雨量序列進(jìn)行分類(lèi)，隨后計(jì)算轉(zhuǎn)移概率矩陣，并進(jìn)行馬氏性檢驗(yàn)。為了確立初始狀態(tài)向量，引入隸屬度的概念，對(duì)參照樣本初始狀態(tài)進(jìn)行測(cè)算。隨后通過(guò)自相關(guān)系數(shù)的計(jì)算來(lái)確定各階權(quán)重。三者結(jié)合即可得到待測(cè)樣本降雨量的初步預(yù)測(cè)值，再經(jīng)過(guò)結(jié)果改進(jìn)即可得到最終預(yù)測(cè)值。

圖1 模型流程Fig.1 The process of model

1.3 模型預(yù)測(cè)步驟

(1)歷史數(shù)據(jù)聚類(lèi)分析。用聚類(lèi)分析代替均值分類(lèi)法對(duì)歷史降雨量數(shù)據(jù)進(jìn)行聚類(lèi)分析，將歷史降雨量數(shù)據(jù)分為5～6組。具體操作在SPSS軟件中進(jìn)行，使用歐氏距離公式測(cè)算各年數(shù)據(jù)間的相似性，運(yùn)用組間聯(lián)結(jié)法進(jìn)行聚類(lèi)[8]。

歐式距離的計(jì)算公式為:

(1)

(2)求轉(zhuǎn)移概率矩陣。根據(jù)聚類(lèi)結(jié)果可將降雨量序列分為N個(gè)狀態(tài)。即狀態(tài)空間E={1,2,…,N}。若本次研究的降雨量序列{Sn}=S1,S2,…,Sn為非負(fù)參數(shù)，且第n+1個(gè)參數(shù)Sn+1的狀態(tài)與其他參數(shù)無(wú)關(guān)，只是Sn的一個(gè)概率函數(shù)P，那么可得到此時(shí)的一步轉(zhuǎn)移概率矩陣：

P(Sn+1=j|Sn=i)=Pij(1) (i,j∈E)

(2)

同理可得步長(zhǎng)為k時(shí)的轉(zhuǎn)移概率矩陣如下[9]：

(3)

其中Pij=fij/fi，fij表示由i狀態(tài)到j(luò)狀態(tài)對(duì)應(yīng)的頻數(shù)，fi表示樣本總數(shù)。在步長(zhǎng)的選擇上，如果選擇的步長(zhǎng)過(guò)大結(jié)果會(huì)明顯失真，為了保證預(yù)測(cè)結(jié)果的準(zhǔn)確性，本次探究步長(zhǎng)k最大取5,即k=1,2,3,4,5[10]。

(3)“馬氏性”檢驗(yàn)。檢驗(yàn)歷史降雨量序列是否具有馬氏性, 是能否應(yīng)用馬爾可夫鏈模型進(jìn)行預(yù)測(cè)的前提，通常用χ2統(tǒng)計(jì)量來(lái)檢驗(yàn)。假設(shè)將一步轉(zhuǎn)移頻數(shù)矩陣的第j列之和與各行各列總和的比值稱(chēng)為“邊際概率”，即：

(4)

統(tǒng)計(jì)量：

(5)

(4)求各參照樣本狀態(tài)向量。模糊馬爾科夫鏈模型與普通馬爾科夫鏈模型相比，不僅僅需要轉(zhuǎn)移概率矩陣，并且需要各參照樣本的初始狀態(tài)向量。在此引入隸屬度的概念，利用隸屬度公式測(cè)算各參照樣本分屬于不同狀態(tài)的隸屬度，從而確定一個(gè)最準(zhǔn)確的初始狀態(tài)[6]。假設(shè)目標(biāo)年份降雨量為A，聚類(lèi)分析的各等級(jí)降雨量下臨界值分別為A0,A1,A2,…,Ax,那么各參照樣本降雨量對(duì)各狀態(tài)的隸屬度計(jì)算如下：

(6)

(5)求待測(cè)樣本各階狀態(tài)向量。將參照樣本年份的狀態(tài)向量與各步轉(zhuǎn)移概率矩陣相乘即可得到待測(cè)樣本年份的各階狀態(tài)向量。例如，本次研究中2011年降雨量數(shù)據(jù)為待測(cè)數(shù)據(jù)，將2010年?duì)顟B(tài)向量與一步轉(zhuǎn)移概率相乘即可得到2011年一階狀態(tài)向量，同理由2009年?duì)顟B(tài)向量以及二步轉(zhuǎn)移概率矩陣可得到2011年二階狀態(tài)向量。

(6)計(jì)算各階自相關(guān)系數(shù)以及權(quán)重。為了正確反映步長(zhǎng)對(duì)馬爾科夫鏈預(yù)測(cè)值的重要性，采用各階自相關(guān)系數(shù)來(lái)計(jì)算權(quán)重，具體公式如下：

(7)

(8)

式中：m為步長(zhǎng)k的最大值，即5。

(8)年降雨量預(yù)測(cè)值。根據(jù)狀態(tài)概率可得到當(dāng)年降雨量預(yù)測(cè)值Xn+1。具體公式如下：

Sn+1=HD/i

(9)

(10)

式中：H為狀態(tài)特征值；Pi為各分類(lèi)狀態(tài)的預(yù)測(cè)概率；D為各預(yù)測(cè)狀態(tài)對(duì)應(yīng)的下臨界值。

(9)結(jié)果改進(jìn)。上述步驟在應(yīng)用時(shí)結(jié)果并不十分理想，因此本次研究在實(shí)際計(jì)算的基礎(chǔ)上對(duì)結(jié)果進(jìn)行了分類(lèi)改進(jìn)。具體公式如下：

(11)

式中：T是馬爾科夫預(yù)測(cè)的待測(cè)樣本所處的狀態(tài)等級(jí)；H為狀態(tài)特征值。

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

本次研究所用數(shù)據(jù)全部來(lái)自中國(guó)氣象數(shù)據(jù)網(wǎng)中的《中國(guó)地面氣候年值數(shù)據(jù)集》。選用了16個(gè)氣象數(shù)據(jù)站的1951-2013年降雨量數(shù)據(jù)，將1951-2010年數(shù)據(jù)作為參照樣本，2011-2013年為待測(cè)樣本。采用基于聚類(lèi)分析的模糊馬爾科夫鏈來(lái)進(jìn)行預(yù)測(cè)。

3 實(shí)例研究

3.1 研究步驟

以江蘇省為例進(jìn)行模型驗(yàn)證如下：

(1)用SPSS軟件將降雨量數(shù)據(jù)進(jìn)行聚類(lèi)分析，共分為枯水年、偏枯年、正常年、偏豐年和豐水年5類(lèi)，結(jié)果如表1所示(1951-2010年數(shù)據(jù)按順序標(biāo)為1～60)。

(2)根據(jù)公式(3)對(duì)分類(lèi)結(jié)果進(jìn)行統(tǒng)計(jì)計(jì)算，得到步長(zhǎng)1～5時(shí)的轉(zhuǎn)移概率矩陣如下，它決定了降雨量指標(biāo)狀態(tài)轉(zhuǎn)移過(guò)程的概率法則。

表1 降雨量序列分類(lèi)Tab.1 Classification of rainfall sequence

(4)根據(jù)表1中分類(lèi)以及公式(6)可得2006-2010年降雨量狀態(tài)向量，再分別與相應(yīng)的概率轉(zhuǎn)移矩陣相乘可得2011年歸一化狀態(tài)向量如下：

(5)根據(jù)公式(7)、(8)可求得各階自相關(guān)系數(shù)以及相應(yīng)的權(quán)重，如表2所示。

表2 各階自相關(guān)系數(shù)以及權(quán)重Tab.2 Autocorrelation coefficients and weights

(7)根據(jù)改進(jìn)后的公式(11)對(duì)2011年降雨量數(shù)值進(jìn)行預(yù)測(cè)，同理重復(fù)步驟(2)～(6)對(duì)2012和2013年降雨量進(jìn)行預(yù)測(cè)，結(jié)果見(jiàn)表3。

表3 南京市降雨量預(yù)測(cè)結(jié)果Tab.3 Rainfall prediction results of Nanjing

3.2 結(jié)果討論

3.2.1 精度分析

使用基于結(jié)果改進(jìn)的聚類(lèi)—模糊馬爾科夫鏈對(duì) 2011-2013年3年降水量進(jìn)行預(yù)測(cè)，結(jié)果顯示3年降雨量的相對(duì)誤差全部都在10%之內(nèi)，平均誤差為8.37%。說(shuō)明建立的馬爾科夫鏈模型可用于降雨量的預(yù)測(cè)，并且結(jié)果改進(jìn)后具有較高精度。

3.2.2 遍歷性與平穩(wěn)分布

由于降雨量的5個(gè)狀態(tài)枯水年、偏枯年、正常年、偏豐年和豐水年是互通的，并且都沒(méi)有周期性，因此全部狀態(tài)構(gòu)成的空間是一個(gè)閉集，那么此馬爾科夫鏈?zhǔn)遣豢杉s的。因此本次研究的馬爾科夫鏈?zhǔn)且粋€(gè)正常返鏈，具有遍歷性，當(dāng)其達(dá)到平穩(wěn)狀態(tài)時(shí)的分布為平穩(wěn)分布，也就是其極限分布。此時(shí)極限分布{πj,j∈E}以及各狀態(tài)出現(xiàn)的周期Ti公式如下：

(12)

(13)

(14)

本次研究中步長(zhǎng)為3時(shí)的模型擬合性最好，因此采用三步轉(zhuǎn)移概率矩陣求平穩(wěn)分布，結(jié)果如表4所示。

表4 平穩(wěn)分布與狀態(tài)周期Tab.4 Stationary distribution and period of state

由表4可知，根據(jù)南京市1951-2013年歷史降雨量數(shù)據(jù)推算，偏枯年降雨?duì)顟B(tài)出現(xiàn)的幾率最大，平均2.47年出現(xiàn)一次，而豐水年出現(xiàn)的幾率最小，平均每隔20.62年才出現(xiàn)一次。這個(gè)數(shù)據(jù)與各狀態(tài)實(shí)際占比幾乎相同，說(shuō)明該周期預(yù)測(cè)合理，可為城市用水規(guī)劃提供依據(jù)。

3.3 拓展預(yù)測(cè)

為了驗(yàn)證模型的正確性以及適用范圍，在數(shù)據(jù)可獲得且完整的基礎(chǔ)上，隨機(jī)取了江蘇省另外兩個(gè)城市降雨測(cè)站以及其他省13個(gè)降雨量測(cè)站數(shù)據(jù)進(jìn)行進(jìn)一步驗(yàn)證。各城市2011-2013年降雨量預(yù)測(cè)結(jié)果以及相對(duì)誤差見(jiàn)表5。

表5 拓展預(yù)測(cè)結(jié)果Tab.5 The results of expanded prediction

由表5以及南京市預(yù)測(cè)數(shù)據(jù)可知，16個(gè)城市共48個(gè)預(yù)測(cè)數(shù)據(jù)中，絕對(duì)誤差超過(guò)15%的年份有15個(gè)，占31.25%；絕對(duì)誤差低于10%的年份有27個(gè)，占56.25%。48個(gè)預(yù)測(cè)數(shù)據(jù)總體相對(duì)誤差為12.4%，總的來(lái)說(shuō)預(yù)測(cè)結(jié)果相對(duì)較好。但由于氣象條件的復(fù)雜性和偶然性，有一些年份會(huì)出現(xiàn)較大誤差，這是不可避免的。從各城市平均誤差來(lái)看，超過(guò)15%的城市有5個(gè)，占31.25%，與高誤差年份占比一致。但誤差低于10%的城市只有7個(gè)，占43.75%，可見(jiàn)個(gè)別極端年份對(duì)整體誤差的影響較大，而這可作為今后降雨量預(yù)測(cè)研究的重要方向。

4 結(jié) 語(yǔ)

本文以降雨量預(yù)測(cè)為研究目標(biāo)，建立了基于聚類(lèi)分析的模糊馬爾科夫鏈降雨量預(yù)測(cè)模型，并以南京市降雨量數(shù)據(jù)為實(shí)例進(jìn)行分析，結(jié)果顯示2011-2013年降雨量預(yù)測(cè)平均誤差為8.37%，預(yù)測(cè)精度較高。隨后將結(jié)果改進(jìn)后的模型推廣到全國(guó)范圍內(nèi)的其他15個(gè)城市，結(jié)果顯示誤差低于10%的城市占43.75%，而誤差低于10%的年份占56.25%。表明模型在預(yù)測(cè)降雨量方面效果較好，可為水資源的規(guī)劃提供一定依據(jù)。但模型在異常值的預(yù)測(cè)上能力較差，整體預(yù)測(cè)精度也因此下滑，如何進(jìn)一步改善極端值年份的影響將是今后的研究重點(diǎn)。