復(fù)權(quán)馬爾可夫鏈及其在怒江水沙預(yù)測中的應(yīng)用

劉新有+彭海英+吳捷+謝飛帆

摘要：[HJ18mm]針對(duì)傳統(tǒng)馬爾可夫鏈及其改進(jìn)的預(yù)測方法只能進(jìn)行狀態(tài)預(yù)測的局限，根據(jù)相依隨機(jī)變量的特點(diǎn)，在以傳統(tǒng)馬爾可夫鏈預(yù)測方法求得各狀態(tài)預(yù)測概率的基礎(chǔ)上，進(jìn)一步以狀態(tài)預(yù)測概率為權(quán)重與狀態(tài)平均值加權(quán)求和，實(shí)現(xiàn)了馬爾可夫鏈預(yù)測方法從狀態(tài)預(yù)測到數(shù)值預(yù)測的關(guān)鍵性改進(jìn)。利用我國西南國際大河怒江干流道街壩水文站1957-2010年徑流和1964-2010年懸移質(zhì)輸沙序列為分析期，2011-2015年徑流和懸移質(zhì)輸沙為驗(yàn)證期，對(duì)所建立的復(fù)權(quán)馬爾可夫鏈預(yù)測方法步驟進(jìn)行驗(yàn)證表明，復(fù)權(quán)馬爾可夫鏈預(yù)測方法具有較高的數(shù)值預(yù)測精度，能夠滿足隨機(jī)時(shí)間序列短期數(shù)值預(yù)測的需要。

關(guān)鍵詞：復(fù)權(quán)馬爾可夫鏈；數(shù)值預(yù)測；徑流；懸移質(zhì)輸沙；怒江

中圖分類號(hào)：P333文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：

16721683（2017）06002607

Abstract：In view of the limitations of traditional Markov chain and its improved prediction methods which can only predict the state，in this paper we realized a critical improvement of the Markov chain forecasting method to being able to conduct numerical predictionWe did so by using weighted summation of the average value of each state multiplied by the corresponding predicted probability，on the basis of obtaining the predicted probability of each state with the traditional Markov chain forecasting method according to the characteristics of dependent stochastic variablesThe data of this study were collected from Daojieba hydrological station on the Nujiang river，which is a famous international river in southwest ChinaWe used the runoff series from 1957 to 2010 and the suspended sediment series from 1964 to 2010 for analysis，and used the runoff and suspended sediment series from 2011 to 2015 for validationResults showed that the reweighted Markov chain forecasting had a high accuracy in numerical prediction and could meet the demand of shortterm numerical prediction in stochastic time series

Key words：reweighted Markov chain；numerical prediction；runoff；suspended sediment；Nujiang river

馬爾可夫鏈?zhǔn)嵌砹_斯數(shù)學(xué)家馬爾可夫1906-1912年間提出的一種隨機(jī)事件預(yù)測的重要方法，在教育、經(jīng)濟(jì)、生物、農(nóng)業(yè)、災(zāi)害、水文氣象、環(huán)境預(yù)測等眾多領(lǐng)域得到了廣泛應(yīng)用。尤其在水文氣象預(yù)測中，馬爾可夫鏈預(yù)測方法應(yīng)用非常廣泛，并在應(yīng)用過程中不斷得以改進(jìn)，加權(quán)馬爾可夫鏈[112]、灰色馬爾可夫鏈[1317]、疊加馬爾可夫鏈[18]、時(shí)間序列馬爾可夫模型[19]、基于多重轉(zhuǎn)移概率的馬爾可夫模型[20]均取得了較好的預(yù)測精度。夏樂天等[2122]系統(tǒng)研究了各種馬爾可夫鏈預(yù)測方法在水文預(yù)測中的應(yīng)用，并對(duì)比了三種常用馬爾可夫鏈預(yù)測方法的優(yōu)劣，認(rèn)為加權(quán)馬爾可夫鏈預(yù)測方法精度最高。這些研究為馬爾可夫鏈預(yù)測方法的應(yīng)用和發(fā)展起到了積極作用，但這些改進(jìn)方法仍然沒有超出對(duì)隨機(jī)事件狀態(tài)預(yù)測的范疇。因此，如何根據(jù)馬爾可夫鏈預(yù)測狀態(tài)概率分布得到預(yù)測值仍然有待解決[12]。本文在加權(quán)馬爾可夫鏈預(yù)測、基于絕對(duì)分布的馬爾可夫鏈預(yù)測和疊加馬爾可夫鏈預(yù)測方法的基礎(chǔ)上，進(jìn)一步以狀態(tài)預(yù)測概率為權(quán)重，結(jié)合狀態(tài)平均值進(jìn)行加權(quán)求和，實(shí)現(xiàn)了馬爾可夫鏈預(yù)測方法從狀態(tài)預(yù)測到數(shù)值預(yù)測的關(guān)鍵性改進(jìn)，并通過怒江水沙預(yù)測實(shí)例對(duì)復(fù)權(quán)馬爾可夫鏈預(yù)測方法的數(shù)值預(yù)測精度進(jìn)行驗(yàn)證。

1復(fù)權(quán)馬爾可夫鏈預(yù)測方法

馬爾可夫鏈通過統(tǒng)計(jì)隨機(jī)事件過去一定時(shí)期內(nèi)的狀態(tài)轉(zhuǎn)移概率來預(yù)測將來狀態(tài)變化的概率，其中時(shí)間參數(shù)集T={0，1，2，…}及狀態(tài)參數(shù)集E={0，1，2，…}稱為馬爾可夫鏈。在實(shí)際應(yīng)用中，一般采用齊次馬爾可夫鏈，即對(duì)任意參數(shù)u，k∈T，有

Pij（u；k）∈E（1）

式中：Pij（u；k）表示隨機(jī)事件u時(shí)段所處的狀態(tài)i，經(jīng)過k步狀態(tài)轉(zhuǎn)移后變?yōu)闋顟B(tài)j的概率。

傳統(tǒng)齊次馬爾可夫鏈的狀態(tài)轉(zhuǎn)移步長一般取1，即利用初始分布推算未來狀態(tài)的絕對(duì)分布，沒有考慮各種步長馬爾可夫鏈的絕對(duì)分布在預(yù)測中所起的作用。為彌補(bǔ)這一缺陷，一些學(xué)者將各種步長馬爾可夫鏈求得的狀態(tài)絕對(duì)分布疊加起來進(jìn)行狀態(tài)預(yù)測，但在疊加過程中沒有考慮各種步長在權(quán)重上的差異。因此，利用各種步長自相關(guān)性的強(qiáng)弱確定不同步長權(quán)重的加權(quán)馬爾可夫鏈進(jìn)行狀態(tài)預(yù)測更符合實(shí)際[9]。但由于加權(quán)馬爾可夫鏈得到的預(yù)測結(jié)果仍然是狀態(tài)，在實(shí)際應(yīng)用中受到一定的限制。復(fù)權(quán)馬爾可夫鏈在之前的研究基礎(chǔ)上，進(jìn)一步以各狀態(tài)的預(yù)測概率為權(quán)重，結(jié)合其對(duì)應(yīng)狀態(tài)均值進(jìn)行加權(quán)求和，從而實(shí)現(xiàn)從狀態(tài)預(yù)測到數(shù)值預(yù)測的跨越。endprint

2復(fù)權(quán)馬爾可夫鏈預(yù)測方法步驟

復(fù)權(quán)馬爾可夫鏈以馬爾科夫鏈求得的各狀態(tài)的預(yù)測概率為基礎(chǔ)，因此步驟（1）至（9）與加權(quán)馬爾可夫鏈基本一致。但為方便對(duì)復(fù)權(quán)馬爾可夫鏈預(yù)測的理解，本研究以加權(quán)馬爾科夫鏈為基礎(chǔ)，完整介紹復(fù)權(quán)馬爾可夫鏈預(yù)測方法步驟。

（1）初步判斷對(duì)象序列是否是隨機(jī)變量。若受大型水利工程等人為控制則不適用于馬爾可夫鏈，反之則可能適用于馬爾可夫鏈，最終確定是否適用于馬爾可夫鏈有待馬氏性檢驗(yàn)結(jié)果。

（2）建立序列狀態(tài)分級(jí)標(biāo)準(zhǔn)，確定資料序列的對(duì)應(yīng)狀態(tài)。常用的狀態(tài)分級(jí)方法有聚類分析法、樣本均值標(biāo)準(zhǔn)差分級(jí)法、頻率曲線法等。水文分析中常用PIII型頻率曲線法來確定各年份的豐枯狀態(tài)，且為使樣本序列具有代表性，一般要求樣本序列不應(yīng)少于30年。根據(jù)狀態(tài)分級(jí)標(biāo)準(zhǔn)，即可確定資料序列所對(duì)應(yīng)的狀態(tài)。

（3）用fij表示指標(biāo)值序列x1，x2，…xn中從狀態(tài)i經(jīng)過一步或多步轉(zhuǎn)移到達(dá)狀態(tài)j的頻數(shù)，i，j∈E。對(duì)資料序列所對(duì)應(yīng)的狀態(tài)進(jìn)行統(tǒng)計(jì)計(jì)算，得到各狀態(tài)的轉(zhuǎn)移規(guī)律，進(jìn)而建立各階（步長）的狀態(tài)轉(zhuǎn)移頻數(shù)矩陣。

（4）將狀態(tài)轉(zhuǎn)移頻數(shù)矩陣（fij）i，j∈E的第i行第j列元素fij除以各行的總和所得的值稱為“轉(zhuǎn)移概率”，記為Pij，i，j∈E，即

Pij=[SX（]fij[]∑[DD（]m[]j=1[DD）]fij（2）

式中：m為指標(biāo)值序列包含的可能的狀態(tài)。

（5）對(duì)隨機(jī)變量進(jìn)行馬氏檢驗(yàn)。

將轉(zhuǎn)移概率矩陣（pij）的第j列之和除以各行各列的總和所得的值稱為“邊際概率”，記為Pj，即

Pj=[SX（]∑[DD（]m[]i=1[DD）]fij[]∑[DD（]m[]i=1[DD）]∑[DD（]m[]j=1[DD）]fij（3）

則當(dāng)序列長度n充分大時(shí)，統(tǒng)計(jì)量

X2=2∑[DD（]m[]i=1[DD）]∑[DD（]m[]j=1[DD）]fij[JB（|]lg[SX（]Pij[]Pj[JB）|]（4）

給定顯著性水平α，查表可得分位點(diǎn)X2α（（m-1）2）的值，計(jì)算后得統(tǒng)計(jì)量X2的值，若X2>X2α·（（m-1）2），則可以認(rèn)為{Xi}符合馬氏性，否則可以認(rèn)為該序列不可作為馬氏鏈來處理。

[JP3]（6）計(jì)算各階（步長）自相關(guān)系數(shù)。計(jì)算公式如下：

rk=∑[DD（]n-k[]l=1[DD）]（xl-[AKx-]）（xl+k-[AKx-]）∑[DD（]n[]l=1[DD）]（xl-[AKx-]）2（5）

式中：rk為第k步長自相關(guān)系數(shù)；xl為序列的第l個(gè)值；[AKx-]為序列均值；n為序列長度。

（7）各步長自相關(guān)系數(shù)規(guī)范化。計(jì)算公式如下：

wk=|rk|∑[DD（]c[]k=1[DD）]|rk|（6）

式中：wk為規(guī)范化后的各步長自相關(guān)系數(shù)，即各步長的馬爾可夫鏈權(quán)重；c為按預(yù)測需要的最大步長。

（8）以各種步長為初始狀態(tài)，結(jié)合其對(duì)應(yīng)的轉(zhuǎn)移概率矩陣，預(yù)測其狀態(tài)概率Pki。

（9）將同一狀態(tài)的各預(yù)測概率加權(quán)求和，得到該狀態(tài)的預(yù)測概率，即

Pi=∑[DD（]m[]k=1[DD）]wkPki（7）

（10）以各狀態(tài)的預(yù)測概率Pi為權(quán)重，與其對(duì)應(yīng)狀態(tài)的均值[AKx-]i加權(quán)求和，得到預(yù)測值d，即

d=∑Pi[AKx-]i（8）

將預(yù)測值加入原序列，再重復(fù)以上步驟，即可進(jìn)行下一步的數(shù)值預(yù)測。

基于絕對(duì)分布的復(fù)權(quán)馬爾可夫鏈預(yù)測和疊加復(fù)權(quán)馬爾可夫鏈預(yù)測方法與基于加權(quán)馬爾可夫鏈的復(fù)權(quán)馬爾可夫鏈預(yù)測方法相似，即在各自的狀態(tài)預(yù)測概率基礎(chǔ)上[21]加上步驟（10）求得預(yù)測值。

3怒江水沙預(yù)測應(yīng)用實(shí)例

怒江薩爾溫江是全球最典型的南北向發(fā)育國際大河，其上游中國境內(nèi)稱為怒江。怒江流域?qū)賺{谷地形，南北跨度大，獨(dú)特地理環(huán)境和氣候條件使其成為全球生物多樣性最突出的地區(qū)之一，怒江也蘊(yùn)藏了極為豐富的水能資源。但由于多種原因，怒江干流水電開發(fā)一直未能實(shí)施，其水文過程至今沒有受到水利工程等人類活動(dòng)的控制。本文以怒江干流道街壩水文站1957-2015年徑流和1964-2015年懸移質(zhì)輸沙序列為數(shù)據(jù)基礎(chǔ)，并將1957-2010年徑流和1964-2010年懸移質(zhì)輸沙序列作為預(yù)測方法分析期，將2011-2015年徑流和懸移質(zhì)輸沙作為預(yù)測方法的驗(yàn)證期，以說明復(fù)權(quán)馬爾可夫鏈預(yù)測方法的具體應(yīng)用并檢驗(yàn)預(yù)測精度。道街壩水文站控制流域面積1102[KG-7]萬[KG-9]km2，占中國境內(nèi)怒江干流流域面積的883%，該站徑流和懸移質(zhì)輸沙變化基本能代表怒江干流徑流和懸移質(zhì)輸沙變化特征。

以怒江道街壩站1957-2010年54年徑流量和1964-2010年47年懸疑質(zhì)輸沙為例，預(yù)測2011年徑流量和懸疑質(zhì)輸沙量，以基于加權(quán)馬爾可夫鏈為基礎(chǔ)的復(fù)權(quán)馬爾可夫預(yù)測為例，詳細(xì)介紹其計(jì)算過程。

（1）初步判斷道街壩站年徑流和年懸移質(zhì)輸沙序列是否是隨機(jī)變量。怒江干流水電梯級(jí)開發(fā)尚未實(shí)施，徑流和懸移質(zhì)輸沙沒有受到人為控制。同時(shí)，怒江干流流域云南段涉及5個(gè)縣區(qū)，但2014年末總?cè)丝趦H15965萬人，社會(huì)經(jīng)濟(jì)發(fā)展落后，加之山高水低，耕地少且分散，區(qū)域內(nèi)農(nóng)業(yè)以自然耕種為主，產(chǎn)流產(chǎn)沙條件基本保持天然狀態(tài)。因此，可初步判斷怒江徑流和懸移質(zhì)輸沙序列屬隨機(jī)變量。

（2）建立道街壩站年徑流和年懸移質(zhì)輸沙序列分級(jí)標(biāo)準(zhǔn)。徑流和年懸移質(zhì)輸沙序列長度超過30年，樣本具有代表性，宜采用PIII型分布頻率曲線法來確定其所處狀態(tài)。分別以保證率0～125%、>125%～375%、>375%～625%、>625%～875%、>875%～100%將年徑流和年懸移質(zhì)輸沙分為豐、偏豐、平、偏少、少5級(jí)，對(duì)應(yīng)狀態(tài)E={1，2，3，4，5}。年徑流和年懸移質(zhì)輸沙PIII型分布各保證率對(duì)應(yīng)的數(shù)值見表1。endprint

（3）按照分級(jí)標(biāo)準(zhǔn)，確定年徑流和年懸移質(zhì)輸沙序列對(duì)應(yīng)的狀態(tài)（表2）。

（4）據(jù)表2進(jìn)行統(tǒng)計(jì)分析，得到1至5階（步長）狀態(tài)轉(zhuǎn)移頻數(shù)矩陣（表3、表4）。

注：矩陣a，b，c，d，e分別為步長1，2，3，4，5的馬爾科夫轉(zhuǎn)移頻數(shù)矩陣，下同。

（5）對(duì)1至5階（步長）狀態(tài)轉(zhuǎn)移頻數(shù)矩陣進(jìn)行統(tǒng)計(jì)分析，得到各個(gè)步長的馬爾科夫鏈轉(zhuǎn)移概率矩陣（表5、表6）。

（6）結(jié)合步長為1的轉(zhuǎn)移概率矩陣和式（3）、式（4），求得怒江道街壩站64年徑流量和47年懸疑質(zhì)輸沙量序列對(duì)應(yīng)的邊際概率和統(tǒng)計(jì)量x2，計(jì)算得到x2值分別為3029和3471，大于α=005顯著性

水平下分位點(diǎn)X2α（16）的值26296，因此該徑流量和懸疑質(zhì)輸沙序列滿足馬氏性。

（7）按照式（5）、式（6）分別計(jì)算各步長自相關(guān)系數(shù)和馬爾可夫鏈權(quán)重，結(jié)果如表7所示。

（8）以各種滯時(shí)為初始狀態(tài)，結(jié)合相應(yīng)的轉(zhuǎn)移概率矩陣預(yù)測其狀態(tài)概率。依據(jù)2010、2009、2008、2007、2006年的年徑流和年懸移質(zhì)輸沙量及其相應(yīng)的狀態(tài)轉(zhuǎn)移概率矩陣，結(jié)合式（7）將同一狀態(tài)的各預(yù)

測概率加權(quán)求和即可對(duì)2011年的年徑流和年懸移質(zhì)輸沙狀態(tài)概率進(jìn)行預(yù)測（表8、表9）。

（9）將各狀態(tài)的預(yù)測概率作為權(quán)重，與其對(duì)應(yīng)狀態(tài)的均值依據(jù)式（8）進(jìn)行加權(quán)求和，即可得到2011年徑流和年懸移質(zhì)輸沙的預(yù)測值（表10）。

（10）由表10可知，2011年徑流預(yù)測值1 766 m3s與實(shí)測值對(duì)比，相對(duì)預(yù)測誤差為208%；2011年懸移質(zhì)輸沙預(yù)測值1 234 kgs與實(shí)測值對(duì)比，相對(duì)預(yù)測誤差為508%。

將預(yù)測值加入原序列，重復(fù)以上步驟，即可得到2012-2015年徑流和懸疑質(zhì)輸沙量的預(yù)測值?；诮^對(duì)分布的馬爾可夫鏈預(yù)測和疊加馬爾可夫鏈預(yù)測方法的復(fù)權(quán)馬爾可夫預(yù)測方法與基于加權(quán)馬爾科夫鏈的復(fù)權(quán)馬爾科夫預(yù)測方法相似，即先求得各狀態(tài)的預(yù)測概率，再以預(yù)測概率為權(quán)重，結(jié)合數(shù)據(jù)序列中各對(duì)應(yīng)狀態(tài)的均值加權(quán)求和，即求得數(shù)值預(yù)測結(jié)果。

表11是在加權(quán)馬爾可夫鏈預(yù)測、基于絕對(duì)分布的馬爾可夫鏈預(yù)測和疊加馬爾可夫鏈預(yù)測方法的基礎(chǔ)上，以各狀態(tài)預(yù)測概率為權(quán)重，結(jié)合狀態(tài)平均值進(jìn)行加權(quán)求和的復(fù)權(quán)馬爾科夫預(yù)測方法的數(shù)值預(yù)測結(jié)果。由表11可知，2011-2015年預(yù)測值的徑流和懸疑質(zhì)輸沙量序列的馬爾可夫檢驗(yàn)統(tǒng)計(jì)量X2均大于26926，說明預(yù)測所用的時(shí)間序列在005顯著性水平下均滿足馬氏性。預(yù)測值與實(shí)測值對(duì)比表明，徑流數(shù)值預(yù)測精度總體高于懸移質(zhì)輸沙數(shù)值預(yù)測精度，這可能是由于相對(duì)于徑流量而言，懸移質(zhì)輸沙受人類活動(dòng)影響更大，導(dǎo)致怒江干流懸移質(zhì)輸沙狀態(tài)之間的數(shù)值跨度較大，1964-2015年期間怒江輸沙極值比達(dá)719。2011-2014年怒江干流徑流數(shù)值預(yù)測精度相對(duì)較高，而2015年徑流數(shù)值預(yù)測精度較低；2011-2013年怒江干流輸沙數(shù)值預(yù)測精度相對(duì)較高，而2014-2015年輸沙數(shù)值預(yù)測精度較低，這可能與復(fù)權(quán)馬爾科夫鏈更適合短期預(yù)測有關(guān)，隨著預(yù)測時(shí)間的延長，預(yù)測誤差可能被逐步放大。本研究短期數(shù)值預(yù)測結(jié)果精度與馬占青等[23]基于馬爾可夫鏈預(yù)測模型的杭州市降水量數(shù)值預(yù)測結(jié)果精度相當(dāng)，高于馬建琴等[24]的改進(jìn)型灰色馬爾可夫鏈模型對(duì)三門峽入庫年徑流的預(yù)測精度，為較少受人類活動(dòng)控制的河流的徑流和輸沙的數(shù)值預(yù)測提供了一條值得探索的途徑。

4結(jié)論

已有的馬爾可夫鏈預(yù)測方法多限于進(jìn)行狀態(tài)預(yù)測，而本文建立的復(fù)權(quán)馬爾可夫鏈預(yù)測方法能夠進(jìn)行數(shù)值預(yù)測，實(shí)現(xiàn)了對(duì)馬爾可夫鏈預(yù)測方法的關(guān)鍵性改進(jìn)，不僅提高了預(yù)測精度，也擴(kuò)展了該方法的應(yīng)用范圍。不受人為控制的隨機(jī)性序列和足夠的序列長度，是適用于馬爾馬爾可夫鏈的前提條件。復(fù)權(quán)馬爾可夫鏈在馬爾可夫鏈前期研究的基礎(chǔ)上，進(jìn)一步以各狀態(tài)的預(yù)測概率為權(quán)重，結(jié)合狀態(tài)平均值進(jìn)行二次加權(quán)求和，從而實(shí)現(xiàn)數(shù)值預(yù)測。與其他馬爾可夫鏈改進(jìn)方法相比，復(fù)權(quán)馬爾可夫鏈能更充分地挖掘隨機(jī)序列的信息。怒江干流水沙預(yù)測實(shí)例表明，所建立的復(fù)權(quán)馬爾可夫鏈預(yù)測方法思路清晰、物理概念明確、計(jì)算簡便，為提高隨機(jī)變量的數(shù)值預(yù)測精度提供了一種可行的途徑。

參考文獻(xiàn)（References）：

[1]賀娟，王曉松，王彩云加權(quán)馬爾可夫鏈模型在密云水庫入庫流量中的應(yīng)用[J]南水北調(diào)與水利科技，2015，13（4）：618621（HE J，WANG X S，WANG C YApplication of the weighted Markov chain model in the inflow prediction of the Miyun Reservoir[J]SouthtoNorth Water Transfers and Water Science & Technology，2015，13（4）：618621（in Chinese）） DOI：1013476jcnkinsbdqk201504003

[2]王濤，錢會(huì)，李培月加權(quán)馬爾可夫鏈在銀川地區(qū)降雨量預(yù)測中的應(yīng)用[J]南水北調(diào)與水利科技，2010，8（1）：7881（WANG T，QIAN H，LI P YPrediction of precipitation based on the weighted Markov chain in Yinchuan area[J]SouthtoNorth Water Transfers and Water Science & Technology，2010，8（1）：7881（in Chinese）） DOI：103969jissn16721683201001021

[3]王亞雄，黃淑嫻，劉祖發(fā)，等變化環(huán)境下北江下游年徑流量的加權(quán)馬爾可夫鏈預(yù)測[J]生態(tài)環(huán)境學(xué)報(bào)，2011，20（4）：754760（WANG Y X，HUANG S X，LIU Z F，et alForecast of yearly river runoff in lower reaches of Beijiang River by Weighted MarkovChain Method in changing environments[J]Ecology and Environmental Sciences，2011，20（4）：754760（in Chinese）） DOI：103969jissn16745906201104030endprint