，

(1.河北水利電力學院 水利工程系，河北 滄州 061001；2.寧夏水利科學研究院，銀川 750021)

基于隸屬度修正的加權(quán)馬爾可夫鏈的降水預(yù)測

苗正偉1，徐利崗2

(1.河北水利電力學院 水利工程系，河北 滄州 061001；2.寧夏水利科學研究院，銀川 750021)

應(yīng)用Fisher最優(yōu)分割法將榆林地區(qū)1951—2015年的年降水序列劃分為9個狀態(tài)，采用規(guī)范化的各階自相關(guān)系數(shù)為權(quán)重，建立了加權(quán)馬爾可夫鏈模型。以屬于同一狀態(tài)的所有降水量的均值作為聚類中心，應(yīng)用模糊C均值聚類(Fuzzy C-Means，F(xiàn)CM)算法中的隸屬函數(shù)計算隸屬度，以隸屬度向量作為預(yù)測時的初始狀態(tài)向量。該模型逐年預(yù)測了榆林市2006—2015年的降水狀態(tài)，結(jié)果全部與實際情況一致。基于馬爾可夫鏈的預(yù)測結(jié)果，采用模糊集中的級別特征值理論分別預(yù)測了2006—2015年的降水量，所有預(yù)測結(jié)果的相對誤差都在10%以內(nèi)，初步表明基于隸屬度修正的加權(quán)馬爾可夫鏈模型是合理可行的。

降水預(yù)測；隸屬度；加權(quán)馬爾可夫鏈；模糊集；榆林

1 研究背景

如何科學劃分降水狀態(tài)，是用加權(quán)馬爾可夫鏈進行預(yù)測的重要問題，最常用的方法是均值-標準差法，但該方法在參數(shù)選擇上主觀性較強，尤其當標準差較大時，不同參數(shù)值會使狀態(tài)區(qū)間范圍變化幅度很大。Fisher最優(yōu)分割法不但能使狀態(tài)劃分更客觀合理，而且還能為分類數(shù)的選擇提供數(shù)學依據(jù)。

此外，傳統(tǒng)的馬爾可夫鏈是將某時刻所處的實際狀態(tài)作為單位向量，然后以該單位向量為初始向量進行預(yù)測[1]，這種處理方法明顯與實際情況不符：設(shè)某時刻的降水量處于偏豐范圍的上限，此時它已很接近豐水年的水平，而且降水狀態(tài)本就沒有具體的數(shù)值范圍來絕對劃分，這種情況下若毫不考慮豐水狀態(tài)的概率，是不合理的。

鑒于此，本文利用模糊C均值(Fuzzy C-Means，F(xiàn)CM)算法中的隸屬函數(shù)，賦予每個年降水量屬于各狀態(tài)的隸屬度，從而使用隸屬度向量作為初始向量進行預(yù)測。

2 基本原理

2.1 Fisher最優(yōu)分割法

Fisher最優(yōu)分割法是在保持序列順序不變的前提下[2]，以類內(nèi)離差平方和最小為目標的一種聚類算法。基本步驟如下：

(1) 對于數(shù)據(jù)為一維變量的序列X，設(shè)某分法中某一類包含的數(shù)據(jù)為{Xi,Xi+1,…,Xj}，則定義該類的類直徑為

(1)

式中:i,j為數(shù)據(jù)在原序列中的序號。

(2) 用f(n,k)表示將n個有序樣品分為k類的某種分法，該分法的目標函數(shù)定義為

(2)

式中:it表示第t類中第一個數(shù)據(jù)的序號;it+1-1表示第t類中最后一個數(shù)據(jù)的序號。

最優(yōu)分割要求目標函數(shù)最小，即：L[p(n,k)]=min{L[f(n,k)]}。

(3) 可以證明，將n個有序數(shù)據(jù)分成k類的最優(yōu)分割，一定是先將前j-1個數(shù)據(jù)最優(yōu)分割為k-1類，剩下的數(shù)據(jù)自成一類，即[3]

L[p(n,k)]=min{L[p(j-1,k-1)]}+D(j,n),

k≤j≤n。

(3)

可見最優(yōu)分割算法的重點就是上述遞推公式，須編程實現(xiàn)。

(4) 分類數(shù)的確定。在分類數(shù)不能提前確定時，可以通過L[p(n,k)]-k關(guān)系圖，以某k值時曲線是否有明顯轉(zhuǎn)折及大于某k值時曲線是否變得很平緩為標準來選擇合理的分類數(shù)。

2.2 加權(quán)馬爾可夫鏈模型

關(guān)于加權(quán)馬爾可夫鏈的基本原理和建模步驟參見文獻[4]至文獻[6]。

2.3 模糊C均值算法中的隸屬函數(shù)

用模糊集理論來處理聚類問題，即為模糊聚類分析。在模糊聚類算法中，應(yīng)用最廣泛且較成功的是FCM算法，該算法中隸屬函數(shù)的規(guī)定如下[7-8]：

?j=1,2,…,N；

(4)

(5)

式中:uij表示第j個樣本數(shù)據(jù)對第i個聚類中心的隸屬度；dij=‖ci-Xj‖為第i個聚類中心與第j個樣本數(shù)據(jù)間的歐氏距離；ci為第i個模糊聚類的聚類中心；m是一個模糊加權(quán)指數(shù)，通常取為2[9]。

2.4 模糊集理論的級別特征值

由下式定義級別特征值[10-11]h,即

(6)

式中:k為狀態(tài)；S為狀態(tài)空間;φk為狀態(tài)k的權(quán)；pk,pr分別為狀態(tài)k,r的預(yù)測概率；θ為最大概率作用系數(shù)。

降水量預(yù)測值Xk由下式給出，即

(7)

式中:Uk,Dk分別為狀態(tài)k所在區(qū)間的上下限;h為級別特征值。

3 預(yù)測模型的建立步驟

(1) 應(yīng)用Fisher最優(yōu)分割法劃分降水序列的狀態(tài)。

(2) 以各狀態(tài)中所有降水量的算術(shù)平均值作為聚類中心，采用式(5)計算降水量對于每個狀態(tài)的隸屬度。

(3) 權(quán)重定義為

(8)

式中:k表示步長;ωk表示步長為k的轉(zhuǎn)移概率的權(quán)重；b為按預(yù)測需要計算到的最大步長;rk表示第k階自相關(guān)系數(shù)，其表達式為

(9)

(4) 計算轉(zhuǎn)移矩陣：因降水因素極其復(fù)雜且規(guī)律性不明顯，所以很難準確給出降水狀態(tài)的轉(zhuǎn)移概率，當數(shù)據(jù)系列足夠長時，可用轉(zhuǎn)移頻率近似代替轉(zhuǎn)移概率，步長為1的轉(zhuǎn)移頻率為

(10)

式中fij表示年降水量從狀態(tài)i一步轉(zhuǎn)移到j(luò)的頻次。與式(10)同理，由樣本序列統(tǒng)計各狀態(tài)間各步長的轉(zhuǎn)移頻率，各步長的轉(zhuǎn)移頻率構(gòu)成相應(yīng)步長的轉(zhuǎn)移矩陣，它決定了降水狀態(tài)轉(zhuǎn)移的概率法則。

(5) 以距離待預(yù)測年份步長為k的年降水量的隸屬度為初始向量，結(jié)合相應(yīng)的k步轉(zhuǎn)移矩陣和k階權(quán)重，即可得到基于步長k的降水數(shù)據(jù)預(yù)測出的概率向量Tk，k=1，2，…，b。

(6) 將各步長的概率向量求和,則有

(11)

向量T中各元素的最大值所對應(yīng)的狀態(tài)即為預(yù)測狀態(tài)。

(7) 應(yīng)用級別特征值公式計算降水量的預(yù)測值。

4 模型在榆林市年降水量預(yù)測中的應(yīng)用

4.1 研究區(qū)概況

榆林市位于陜西省北部，位置在107°28′E—111°15′E，36°57′N—39°35′N，地勢由西部向東南傾斜，西南部平均海拔1 600～1 800 m，其他各地平均海拔1 000～1 200 m。地貌分為風沙灘區(qū)和黃土丘陵溝壑區(qū)兩大類。氣候為典型的溫帶大陸性季風氣候，區(qū)域多年平均相對濕度57.5%，年平均氣溫為10.7 ℃，極端最高氣溫為38.9 ℃，極端最低氣溫為-24 ℃，氣象災(zāi)害多發(fā)。全市多年平均降水量為397.7 mm，歷年最大年降水量為849.6 mm，最小年降水量為108.6 mm；降水量年內(nèi)分配不均，6—9月份的降水量占全年的75%，降水變差系數(shù)范圍為0.23～0.31；多年平均水面蒸發(fā)量為1 211 mm，多年平均干旱指數(shù)為3.0[12]。

4.2 數(shù)據(jù)來源

從“中國氣象數(shù)據(jù)網(wǎng)”下載了榆林市1951—2015年的逐月降水量，并將之整理為年降水序列，結(jié)果如表1。

表1榆林市1951—2015年降水量分類及聚類中心
Table1ClassificationofannualprecipitationinYulincityfrom1951to2015andcorrespondingclustercenters

狀態(tài)聚類中心/mm降水量區(qū)間/mm1159．1x<203．902271．5203．90≤x<290．853315．9290．85≤x<339．304359．0339．30≤x<387．255414．0387．25≤x<432．806448．4432．80≤x<474．407492．7474．40≤x<512．858552．3512．85≤x<632．059690．5x≥632．05

4.3 隸屬度修正的加權(quán)馬爾可夫鏈模型的建立及應(yīng)用

(1) 首先將數(shù)據(jù)按從小到大的順序排列[13]，然后根據(jù)Fisher最優(yōu)分割的基本原理編制MATLAB程序進行計算。因序列長度為65 a，類直徑及各目標函數(shù)值的計算結(jié)果數(shù)據(jù)較多，限于篇幅，故略去，只給出65個數(shù)據(jù)k類最優(yōu)分割的目標函數(shù)L[p(65,k)]與分類數(shù)k的關(guān)系曲線，如圖1所示。

圖1 榆林市L[p(65,k)]-k曲線Fig.1 The L[p(65,k)]-k curve in Yulin city

普遍把降水狀態(tài)均分為5類，但從圖1可以看出k=5處曲線還有明顯轉(zhuǎn)折；在k=6～8處盡管轉(zhuǎn)折不明顯，但曲線下降的趨勢還是比較大；而在k=9處開始曲線變得很平緩，逐漸逼近橫坐標軸，說明分類數(shù)超過9之后，目標函數(shù)值減小趨勢變得很緩慢，換言之，分類大于9的情況與分成9類相比并無明顯意義。因此，本文將降水序列劃分為9類，并以各狀態(tài)中所有降水量的算術(shù)平均值作為聚類中心，結(jié)果如表1。

(2)按表1確定榆林市1951—2015年各年降水量的狀態(tài)，結(jié)果如表2所示。

(3) 取模糊加權(quán)指數(shù)m=1.8，根據(jù)式(5)計算各年降水量的隸屬度，因數(shù)據(jù)量為9×65=585，限于篇幅，此處不詳述。

(4) 選擇轉(zhuǎn)移矩陣和自相關(guān)系數(shù)的最大步長并計算各階權(quán)重：轉(zhuǎn)移矩陣最大步長普遍根據(jù)經(jīng)驗確定，并無嚴謹理論依據(jù)。本文經(jīng)多次試驗后確定最大步長b=22。計算1—22階自相關(guān)系數(shù)和權(quán)重，結(jié)果如表3所示。

表2 榆林市1951—2015年降水量及其狀態(tài)Table 2 Annual precipitation and states inYulin city from 1951 to 2015

(5) 對降水狀態(tài)序列進行統(tǒng)計計算，得步長為1—22的轉(zhuǎn)移概率矩陣，因篇幅限制，只列出P(1)和P(22)示意，即:

(6) 首先根據(jù)1984—2005年的年降水數(shù)據(jù)預(yù)測2006年的年降水狀態(tài)。

1984—2005年降水量的隸屬度如表4。

求2005年的隸屬度向量u2005與一步轉(zhuǎn)移矩陣P(1)的積，然后再用步長為1的規(guī)范化自相關(guān)系數(shù)ω1進行加權(quán)，即得到基于步長為1的轉(zhuǎn)移矩陣的預(yù)測結(jié)果T1，即：T1=u2005·P(1)·ω1；

求2004年的隸屬度向量u2004與兩步轉(zhuǎn)移矩陣P(2)的積，然后再用步長為2的規(guī)范化自相關(guān)系數(shù)ω2進行加權(quán)，即得到基于步長為2的轉(zhuǎn)移矩陣的預(yù)測結(jié)果T2，即：T2=u2004·P(2)·ω2；

同理：求年份為t的隸屬度向量ut與k步轉(zhuǎn)移矩陣P(k)的積，然后再用步長為k的規(guī)范化自相關(guān)系數(shù)ωk進行加權(quán)，即得到基于步長為k的轉(zhuǎn)移矩陣預(yù)測結(jié)果Tk，即

Tk=ut·P(k)·ωk，

k=2006-t，t=2005,2004,…,1984 。

(12)

將各步長的預(yù)測結(jié)果Tk相加得到

(13)

概率向量T中各元素的最大值對應(yīng)的狀態(tài)即為預(yù)測值,見表5。

表4 1984—2005年降水量隸屬度Table 4 Membership of precipitation from 1984 to 2005

表5 2006年降水狀態(tài)預(yù)測Table 5 Prediction result of precipitation state in 2006

由表5可知，概率向量T中各元素的最大值為0.284 86，對應(yīng)狀態(tài)3，即該模型預(yù)測2006年的降水狀態(tài)編號為3，由表2知，榆林市2006年的降水狀態(tài)編號正是3，預(yù)測正確。

(7) 重復(fù)上述步驟逐年預(yù)測2007—2015年的降水狀態(tài)，結(jié)果如表6。

(8) 在本文狀態(tài)劃分的基礎(chǔ)之上，應(yīng)用傳統(tǒng)加權(quán)馬氏鏈預(yù)測2006—2015年的降水狀態(tài)，結(jié)果如表6。

(9) 狀態(tài)預(yù)測結(jié)果分析：

由表6知，2006—2015年共計10 a的狀態(tài)預(yù)測，本文所建模型的預(yù)測結(jié)果全部與實際一致(準確率達100%)，而傳統(tǒng)加權(quán)馬氏鏈的預(yù)測準確率只有70%。

表6榆林市2006—2015年降水狀態(tài)預(yù)測結(jié)果

Table6PredictionresultofprecipitationstateinYulincityfrom2006to2015

年份實測狀態(tài)編號傳統(tǒng)加權(quán)馬氏鏈預(yù)測隸屬度修正馬氏鏈預(yù)測狀態(tài)預(yù)測值預(yù)測結(jié)果狀態(tài)預(yù)測值預(yù)測結(jié)果200633正確3正確200766正確6正確200866正確6正確200955正確5正確201044正確4正確201163錯誤6正確201288正確8正確201384錯誤8正確201444正確4正確201564錯誤6正確

傳統(tǒng)馬氏鏈對于2013—2015年的預(yù)測結(jié)果全是4，其中有2個是錯誤的，而由表2知狀態(tài)4恰恰是實測序列里出現(xiàn)次數(shù)最多的，達到12次，說明傳統(tǒng)馬氏鏈的預(yù)測結(jié)果會受到狀態(tài)在實測序列里出現(xiàn)頻率大小的影響，從而降低其預(yù)測精度；隸屬度修正的模型克服了這個缺點，使得預(yù)測精度較傳統(tǒng)模型有明顯提高。

同時注意到，2013年的實測降水狀態(tài)編號為8，2014年突然變?yōu)?，跨度達到4個狀態(tài)，體現(xiàn)了降水量波動性大的特點，這也是對降水量精準預(yù)測的難點所在。但在這種情況下，該模型依然預(yù)測準確，說明本文所建模型預(yù)測機制合理，由于隸屬度的修正使其蘊含了更豐富的歷史信息以致于能和降水的大幅度波動性擬合較好。

4.4 降水量預(yù)測及結(jié)果分析

在馬氏鏈預(yù)測結(jié)果基礎(chǔ)之上引入模糊集的級別特征值理論預(yù)測具體降水數(shù)值。因最大概率作用系數(shù)θ的取值沒有成熟的理論依據(jù)，故經(jīng)多次試算，本文確定θ=57，根據(jù)式(7)分別預(yù)測2006—2015年的降水量，結(jié)果如表7。

表7榆林市2006—2015年降水量預(yù)測結(jié)果
Table7PredictionresultofprecipitationinYulincityfrom2006to2015

年份實測值/mm預(yù)測狀態(tài)h值預(yù)測值/mm相對誤差/%2006313．533．0000290．8-7．232007439．465．7491452．42．962008447．566．0000472．15．512009420．854．9997430．32．252010363．944．0000344．2-5．412011445．465．3298419．4-5．842012566．888．0000547．0-3．492013562．587．9141541．2-3．792014378．344．0000344．2-9．012015451．065．3000417．1-7．52

由表7還可知：2006—2015年共計10 a降水量預(yù)測值的相對誤差全在10%以內(nèi)，最大相對誤差也只有-9.01%，相對誤差在5%以內(nèi)的占到了40%。以上結(jié)果表明所建模型合理可行，精度較高，可用于榆林市年降水量的預(yù)測。

4.5 榆林市2016—2021年降水量預(yù)測

用所建模型動態(tài)預(yù)測了榆林市2016—2021年的降水量，結(jié)果如表8。

表8 榆林市2016—2021年降水量預(yù)測結(jié)果Table 8 Prediction result of precipitation inYulin city from 2016 to 2021

由表8可知，2016—2018年連續(xù)3 a時間，榆林市降水量都明顯小于1951—2015年的平均值404.7 mm；2019年降水量較豐； 2020年降水量顯著減少，只有不到300 mm，大致為枯水年水平；2021年降水量大概為平水年水平。

5 結(jié)論與展望

(1) 降水狀態(tài)的劃分沒有絕對的數(shù)值依據(jù)，尤其是區(qū)間兩極的一些降水數(shù)據(jù)，具有明顯的模糊特征，從這個角度講，將降水狀態(tài)模糊化是合理的。本文利用FCM算法中的隸屬函數(shù)賦予每個降水量屬于各狀態(tài)的隸屬度，使降水狀態(tài)模糊化的同時，也使得預(yù)測時的初始向量蘊含了更多的信息，也因此具有了較高的預(yù)測精度。

(2)利用該模型預(yù)測了榆林市2016—2021年的降水量，結(jié)果表明： 6 a之中，除2019年降水較豐、2021年為平水年之外，其余4 a降水量都明顯少于過去65 a的平均水平，尤其是2020年，年降水還不到300 mm,這可能會對榆林市造成一定程度的干旱影響。

(3) 隸屬函數(shù)的構(gòu)造方法有很多，哪一種方法最適于降水序列還需進一步實踐、探討。

(4) 因降水量不能嚴格滿足馬爾可夫模型的一階無后效性，故采用多步長、加權(quán)2種方法進行修正。但對于轉(zhuǎn)移矩陣最大步長的選擇卻一直缺乏嚴謹?shù)睦碚撘罁?jù)，普遍經(jīng)驗性地采用狀態(tài)數(shù)作為最大步長，這樣做在很多情況下都限制了預(yù)測的精度。本文根據(jù)榆林市降水數(shù)據(jù)經(jīng)過多次試驗，確定了最大步長為22，并取得了滿意的預(yù)測結(jié)果。但該步長的選擇具有極強的針對性，榆林之外的地區(qū)未必適用。因此，如何在理論上確定合理的最大步長是馬爾可夫鏈預(yù)測的一個重點問題。

(5) 級別特征值理論中的最大作用系數(shù)θ通常都取為2或4，但這同樣沒有成熟的理論依據(jù)。本文經(jīng)過多次試驗發(fā)現(xiàn)θ>40時才能取得很好的結(jié)果，θ=57時全局結(jié)果相對最佳。因此，最大作用系數(shù)θ的選擇仍需進一步研究。

[1] 黃 華,蔡 仁,穆振俠,等.基于模糊集修正加權(quán)馬爾可夫模型在新疆降水預(yù)測中的應(yīng)用[J].新疆農(nóng)業(yè)科學,2015,52(10):1891-1898.

[2] 李 俊,畢華興,李笑吟,等.有序聚類法在土壤水分垂直分層中的應(yīng)用[J].北京林業(yè)大學學報,2007,29(1):98-101.

[3] 武琳琳.基于Fisher最優(yōu)分割法的聚類分析應(yīng)用[D].鄭州:鄭州大學,2013:13-30.

[4] 溫海彬.馬爾可夫鏈預(yù)測模型及一些應(yīng)用[D].南京:南京郵電大學,2012:7-31.

[5] 楊皓翔,梁 川,崔寧博.基于加權(quán)灰色-馬爾可夫鏈模型的城市需水預(yù)測[J].長江科學院院報,2015,32(7):15-21.

[6] 萬 臣,李建峰,趙 勇,等.基于新維BP神經(jīng)網(wǎng)絡(luò)-馬爾科夫鏈模型的大壩沉降預(yù)測[J].長江科學院院報,2015,32(10):23-27.

[7] 何嘉婧,王晉東,于智勇.基于模糊C-均值的改進人工蜂群聚類算法[J].計算機應(yīng)用研究,2016,33(5):1342-1345.

[8] 李培強,李欣然,陳輝華,等.基于模糊聚類的電力負荷特性的分類與綜合[J].中國電機工程學報,2005,25(24):73-78.

[9] 李 中,苑津莎.不同相似度測量方式的模糊C均值聚類分析[J].計算機工程與應(yīng)用,2011,47(18):17-18.

[10] 孫才志,林學鈺.降水預(yù)測的模糊權(quán)馬爾可夫模型及應(yīng)用[J].系統(tǒng)工程學報,2003,18(4):294-299.

[11] 王 濤,錢 會,李培月.加權(quán)馬爾可夫鏈在銀川地區(qū)降雨量預(yù)測中的應(yīng)用[J].南水北調(diào)與水利科技,2010,8(1):78-81.

[12] 王化齊,張茂省,黨學亞.榆林地區(qū)降水蒸發(fā)時間序列的多尺度特征和突變分析[J].水電能源科學,2010,28(9):1-4.

[13] 姚美娟,陳建平,王 翔,等.基于最優(yōu)分割分級法的月球撞擊坑分級及其演化分析[J].巖石學報,2016,32(1):119-126.

Prediction of Annual Precipitation by Weighted Markov ChainBased on Membership Modification

MIAO Zheng-wei1，XU Li-gang2

(1.Department of Hydraulic Engineering,Hebei University of Water Resources and Electric Engineering,Cangzhou 061001, China; 2. Scientific Research Institute of Water Conservancy of Ningxia, Yinchuan 750021, China)

The annual precipitation series of Yulin city from 1951 to 2015 was divided into 9 states by the Fisher optimal partition method. The weighted Markov chain model was established by taking the standardized autocorrelation coefficients as weights. With the mean value of all precipitation in the same state as the cluster center, the membership function of the Fuzzy C-Means was applied to calculate the membership of annual precipitation, and the membership vector was taken as the initial state vector for the time period. The precipitation state from 2006 to 2015 in Yulin city was predicted year by year. All the results agree with the reality. Based on the prediction results of Markov Chain, the precipitation was predicted respectively from 2006 to 2015 by the level characteristics value of Fuzzy Sets, and the relative error of all the prediction results is less than 10%. The preliminary results show that the model of weighted Markov chain based on membership modification is feasible.

precipitation prediction；membership；weighted Markov chain；fuzzy sets；Yulin

2016-09-01;

2016-09-29

苗正偉(1981-)，男，山東聊城人，講師，碩士，主要從事水文水資源方面的研究。E-mail:myjxbe@126.com

10.11988/ckyyb.20160893

P333

A

1001-5485(2018)01-0040-07

(編輯：王 慰)