基于聚類分析的模糊馬爾科夫鏈在降雨量預測中的應用

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

(北京師范大學環境學院 水環境模擬國家重點實驗室，北京 100875)

降雨量是流域或地區水資源循環系統的重要輸入項目，是最主要的水資源補給途徑，在農業生產、工業建設、水利防洪等方面起著重要的作用。由于氣象條件具有多變性以及復雜性，導致降雨過程存在較大的隨機性[1]，如何較為準確地對降雨量進行預測也成了一個熱點研究問題。目前較為常用的降雨量預測方法有灰色預測[2]、馬爾科夫鏈[3]、人工神經網絡[4]等。

馬爾可夫過程是研究事物狀態及狀態轉移規律的理論。它通過確定系統不同狀態的初始概率及狀態之間的轉移概率關系, 來確定狀態的變化趨勢, 從而進行預測[5]。常用的馬爾科夫鏈有基于絕對分布的馬爾科夫鏈、疊加馬爾科夫鏈、加權馬爾科夫鏈等，其中加權馬爾科夫鏈是效果最好、應用最廣泛的[6]。本次研究在加權馬爾科夫鏈的基礎上進行了改進：首先利用歐氏距離聚類代替常用的均值分類；其次引入隸屬度函數來確定初始狀態向量；最后對預測結果進行分類改進。從而建立了基于聚類分析的模糊馬爾科夫鏈降雨量預測模型，并利用中國氣象站記錄的1951-2013年降雨量數據對模型進行驗證。

1 模型原理及步驟

1.1 馬爾科夫預測原理

馬爾科夫預測是根據馬爾科夫理論與方法來研究未來發展趨勢的一種動態分析方法，既適用于時間序列，又適用于空間序列。其基本特征是無后效性，即系統在未來某一時刻的狀態只與現在的狀態有關，而與過去的狀態無關。因此馬爾科夫鏈在模擬預測隨機波動性較大的問題上有較好的效果[7]。

1.2 聚類-模糊馬爾可夫鏈模型流程

本次研究所采用的聚類-模糊馬爾可夫鏈模型如圖1所示，首先通過聚類分析將降雨量序列進行分類，隨后計算轉移概率矩陣，并進行馬氏性檢驗。為了確立初始狀態向量，引入隸屬度的概念，對參照樣本初始狀態進行測算。隨后通過自相關系數的計算來確定各階權重。三者結合即可得到待測樣本降雨量的初步預測值，再經過結果改進即可得到最終預測值。

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

1.3 模型預測步驟

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

歐式距離的計算公式為:

(1)

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

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

(2)

同理可得步長為k時的轉移概率矩陣如下[9]：

(3)

其中Pij=fij/fi，fij表示由i狀態到j狀態對應的頻數，fi表示樣本總數。在步長的選擇上，如果選擇的步長過大結果會明顯失真，為了保證預測結果的準確性，本次探究步長k最大取5,即k=1,2,3,4,5[10]。

(3)“馬氏性”檢驗。檢驗歷史降雨量序列是否具有馬氏性, 是能否應用馬爾可夫鏈模型進行預測的前提，通常用χ2統計量來檢驗。假設將一步轉移頻數矩陣的第j列之和與各行各列總和的比值稱為“邊際概率”，即：

(4)

統計量：

(5)

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

(6)

(5)求待測樣本各階狀態向量。將參照樣本年份的狀態向量與各步轉移概率矩陣相乘即可得到待測樣本年份的各階狀態向量。例如，本次研究中2011年降雨量數據為待測數據，將2010年狀態向量與一步轉移概率相乘即可得到2011年一階狀態向量，同理由2009年狀態向量以及二步轉移概率矩陣可得到2011年二階狀態向量。

(6)計算各階自相關系數以及權重。為了正確反映步長對馬爾科夫鏈預測值的重要性，采用各階自相關系數來計算權重，具體公式如下：

(7)

(8)

式中：m為步長k的最大值，即5。

(8)年降雨量預測值。根據狀態概率可得到當年降雨量預測值Xn+1。具體公式如下：

Sn+1=HD/i

(9)

(10)

式中：H為狀態特征值；Pi為各分類狀態的預測概率；D為各預測狀態對應的下臨界值。

(9)結果改進。上述步驟在應用時結果并不十分理想，因此本次研究在實際計算的基礎上對結果進行了分類改進。具體公式如下：

(11)

式中：T是馬爾科夫預測的待測樣本所處的狀態等級；H為狀態特征值。

2 數據來源

本次研究所用數據全部來自中國氣象數據網中的《中國地面氣候年值數據集》。選用了16個氣象數據站的1951-2013年降雨量數據，將1951-2010年數據作為參照樣本，2011-2013年為待測樣本。采用基于聚類分析的模糊馬爾科夫鏈來進行預測。

3 實例研究

3.1 研究步驟

以江蘇省為例進行模型驗證如下：

(1)用SPSS軟件將降雨量數據進行聚類分析，共分為枯水年、偏枯年、正常年、偏豐年和豐水年5類，結果如表1所示(1951-2010年數據按順序標為1～60)。

(2)根據公式(3)對分類結果進行統計計算，得到步長1～5時的轉移概率矩陣如下，它決定了降雨量指標狀態轉移過程的概率法則。

表1 降雨量序列分類Tab.1 Classification of rainfall sequence

(4)根據表1中分類以及公式(6)可得2006-2010年降雨量狀態向量，再分別與相應的概率轉移矩陣相乘可得2011年歸一化狀態向量如下：

(5)根據公式(7)、(8)可求得各階自相關系數以及相應的權重，如表2所示。

表2 各階自相關系數以及權重Tab.2 Autocorrelation coefficients and weights

(7)根據改進后的公式(11)對2011年降雨量數值進行預測，同理重復步驟(2)～(6)對2012和2013年降雨量進行預測，結果見表3。

表3 南京市降雨量預測結果Tab.3 Rainfall prediction results of Nanjing

3.2 結果討論

3.2.1 精度分析

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

3.2.2 遍歷性與平穩分布

由于降雨量的5個狀態枯水年、偏枯年、正常年、偏豐年和豐水年是互通的，并且都沒有周期性，因此全部狀態構成的空間是一個閉集，那么此馬爾科夫鏈是不可約的。因此本次研究的馬爾科夫鏈是一個正常返鏈，具有遍歷性，當其達到平穩狀態時的分布為平穩分布，也就是其極限分布。此時極限分布{πj,j∈E}以及各狀態出現的周期Ti公式如下：

(12)

(13)

(14)

本次研究中步長為3時的模型擬合性最好，因此采用三步轉移概率矩陣求平穩分布，結果如表4所示。

表4 平穩分布與狀態周期Tab.4 Stationary distribution and period of state

由表4可知，根據南京市1951-2013年歷史降雨量數據推算，偏枯年降雨狀態出現的幾率最大，平均2.47年出現一次，而豐水年出現的幾率最小，平均每隔20.62年才出現一次。這個數據與各狀態實際占比幾乎相同，說明該周期預測合理，可為城市用水規劃提供依據。

3.3 拓展預測

為了驗證模型的正確性以及適用范圍，在數據可獲得且完整的基礎上，隨機取了江蘇省另外兩個城市降雨測站以及其他省13個降雨量測站數據進行進一步驗證。各城市2011-2013年降雨量預測結果以及相對誤差見表5。

表5 拓展預測結果Tab.5 The results of expanded prediction

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

4 結 語

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