蒙特卡洛與馬爾可夫方法在降水預(yù)測中的應(yīng)用

黃鶴 張維江 李娟

摘 要：為了更深入地分析原州區(qū)的降水特征，為該區(qū)域水資源預(yù)測提供依據(jù)，基于原州區(qū)1957—2016年60 a降水資料，采用蒙特卡洛方法推求降水分布，采用K-S檢驗(yàn)對模型進(jìn)行顯著性檢驗(yàn)，采用基于歐氏距離的層次聚類方法進(jìn)行狀態(tài)劃分，確定了原州區(qū)的降水分布，建立了適用于原州區(qū)的滑動(dòng)平均加權(quán)馬爾可夫預(yù)測模型。根據(jù)已有數(shù)據(jù)驗(yàn)證了預(yù)測結(jié)果的有效性，再結(jié)合已確定的降水分布，通過K-S檢驗(yàn)，檢驗(yàn)了未來5 a降水預(yù)測的準(zhǔn)確性。結(jié)果表明：原州區(qū)降水分布符合P-Ⅲ型分布;馬爾可夫模型適用于原州區(qū)降水預(yù)測，且未來5 a的降水預(yù)測結(jié)果是準(zhǔn)確的，分別為508.5、520.8、554.9、451.0、466.6 mm。

關(guān)鍵詞：蒙特卡洛方法;馬爾可夫模型;隨機(jī)模擬;降水預(yù)測;原州區(qū);K-S檢驗(yàn)

中圖分類號：TV11 文獻(xiàn)標(biāo)志碼：A

doi：10.3969/j.issn.1000-1379.2020.05.004

Abstract： In order to analyze the characteristics of precipitation in YuanzhouDistrict and provide a basis for water resources prediction in the region， based on the precipitation data from 1957 to 2016 of the district， the Monte Carlo method was used to estimate the precipitation distribution and the K-S test was used to model the significant test. It used the hierarchical clustering method based on Euclidean distance to divide the state， determined the precipitation distribution of the region and established a sliding average weighted Markov prediction model which was suitable for Yuanzhou District. Based on the existing data， the true validity of the prediction results was verified. Combined with the determined precipitation distribution， the accuracy of precipitation prediction for the next 5 years was tested by K-S test. Comparing with Pearson three-type distribution， the results show that the precipitation distribution in Yuanzhou District is more consistent with the log-normal distribution. The Markov model can be applied to the precipitation forecast in the region， and the prediction results of the sliding average precipitation in the next 5 years are real and effective， respectively 508.5， 520.8， 554.9， 451.0 and 466.6 mm.

Key words： Monte Carlo method; Markov model; stochastic simulation; precipitation prediction; Yuanzhou District; K-S test

1 引 言

由于客觀世界中的一些現(xiàn)象可能與另一種現(xiàn)象存在著某種相似性，因此我們經(jīng)常從一種現(xiàn)象出發(fā)來研究另一種現(xiàn)象。當(dāng)某個(gè)概率模型可以描述隨機(jī)系統(tǒng)并可以基于此概率模型進(jìn)行實(shí)驗(yàn)時(shí)，這種實(shí)驗(yàn)方法即為隨機(jī)模擬方法[1]。近年來，隨機(jī)模擬在水文系統(tǒng)預(yù)測中發(fā)展迅速，取得了顯著進(jìn)展[2]。劉新立[3]研究了隨機(jī)過程情況下隨機(jī)模擬和馬爾可夫鏈在水災(zāi)風(fēng)險(xiǎn)管理中所起的作用，研究表明兩者結(jié)合可以評估未來若干個(gè)時(shí)間單位水災(zāi)所造成的風(fēng)險(xiǎn);尹正杰等[4]提出了一個(gè)同時(shí)含有趨勢、季節(jié)、隨機(jī)3個(gè)分量的時(shí)間序列模型，并通過此時(shí)間序列模型對灌區(qū)灌溉需水量進(jìn)行隨機(jī)模擬;溫季等[5]探討了集中作物需水量的隨機(jī)模擬及預(yù)測模型，結(jié)果表明隨機(jī)模擬在作物灌溉管理中具有很強(qiáng)的應(yīng)用性及普適性。

蒙特卡洛方法與馬爾可夫方法在水文系統(tǒng)方面的應(yīng)用發(fā)展十分迅速。ZHANG等[6]研究開發(fā)了基于馬爾可夫鏈-蒙特卡洛的多級因子分析方法，從而更好地對水文模型參數(shù)不確定性進(jìn)行評估;KNIGHTON等[7]通過蒙特卡洛方法將建立的概率分布應(yīng)用于已知模型，用來估計(jì)水文模型的不確定性;劉悅憶等[8]在淮河流域水動(dòng)力學(xué)-水質(zhì)模型的基礎(chǔ)上，使用蒙特卡洛方法隨機(jī)模擬了大量入流數(shù)據(jù)并將其用作模型的輸入條件進(jìn)行計(jì)算，建立了基于蒙特卡洛模擬的水質(zhì)概率預(yù)報(bào)模型;李娟等[9]應(yīng)用滑動(dòng)平均-馬爾可夫模型對固原市隆德縣的降水進(jìn)行預(yù)測，研究證明改進(jìn)后的馬爾可夫模型預(yù)測精度較高;王艷等[10]在傳統(tǒng)的馬爾可夫方法上用最優(yōu)分割法優(yōu)化了分級標(biāo)準(zhǔn);岳遙等[11]提出了一種基于投影距離的處理級別特征值的方法來代替?zhèn)鹘y(tǒng)的處理方法，并將其引入馬爾可夫模型，從而應(yīng)用于對水質(zhì)的定量預(yù)測;李亞斌等[12]用樣本均值-均方差方法對銅川地區(qū)降水量進(jìn)行分級，建立相應(yīng)的馬爾可夫模型進(jìn)行降水預(yù)測。

降水受氣候和人類活動(dòng)兩者的共同影響，且由于氣候因素本身具有很強(qiáng)的變異性、復(fù)雜性以及多樣性，因此降水系統(tǒng)呈現(xiàn)出十分復(fù)雜的行為特征。過于復(fù)雜的系統(tǒng)難以建立精確的數(shù)學(xué)模型進(jìn)行準(zhǔn)確預(yù)測。筆者運(yùn)用隨機(jī)模擬原理，在確定地區(qū)降水的分布模型后，結(jié)合蒙特卡洛方法以及馬爾可夫方法對固原市原州區(qū)的降水進(jìn)行隨機(jī)模擬，以期預(yù)測未來原州區(qū)降水的發(fā)展變化趨勢，為水資源合理利用及調(diào)控提供依據(jù)。

2 研究方法

2.1 蒙特卡洛方法

蒙特卡洛方法又稱為概率統(tǒng)計(jì)法，是一種基于概率論思想，對隨機(jī)變量進(jìn)行數(shù)理統(tǒng)計(jì)實(shí)驗(yàn)及分布概率模擬，從而近似求解得到預(yù)測值的方法[13]。只要能構(gòu)建出適當(dāng)?shù)哪Ｐ?，此方法都能夠進(jìn)行模擬應(yīng)用，其基本框架：①假設(shè)變量X服從某一概率分布;②用隨機(jī)抽樣的方法對概率分布進(jìn)行抽樣從而得到樣本值，一般選擇常用的均勻分布模型，產(chǎn)生（0，1）區(qū)間內(nèi)的隨機(jī)數(shù)，再結(jié)合原始數(shù)據(jù)產(chǎn)生服從特定分布的隨機(jī)序列，即樣本值;③確定和選取統(tǒng)計(jì)值;④由統(tǒng)計(jì)量的算術(shù)平均值得到統(tǒng)計(jì)量的估計(jì)值，從而近似求解出預(yù)測值。

2.2 馬爾可夫方法

2.2.1 馬爾可夫基本原理

馬爾可夫模型是基于馬爾可夫鏈建立起來的，馬爾可夫鏈?zhǔn)且环N特殊的隨機(jī)過程[14]。設(shè)有一隨機(jī)運(yùn)動(dòng)的系統(tǒng)，它可能處的狀態(tài)記為E0、E1、…、En。這個(gè)系統(tǒng)只可能在時(shí)刻t（t=1，2，…，n）上改變它的狀態(tài)。隨著隨機(jī)運(yùn)動(dòng)的進(jìn)行，定義一列隨機(jī)變量Xn（n=0，1，2，…），其中Xn=k，表示在t=n時(shí)，系統(tǒng)狀態(tài)為Ek。

2.2.2 模糊集理論中的級別特征值

使用傳統(tǒng)的馬爾可夫只能夠預(yù)測到某個(gè)區(qū)間，無法預(yù)測較為準(zhǔn)確的降水值，在實(shí)際應(yīng)用中作用有限，筆者選用模糊集理論中的級別特征值方法有效地解決了這個(gè)問題。

假定最大概率的狀態(tài)為i，當(dāng)H>i時(shí)，年降水量的預(yù)測值X預(yù)報(bào)=TiHi+0.5;當(dāng)H

2.3 K-S分布檢驗(yàn)

K-S檢驗(yàn)是用于檢測一組數(shù)據(jù)的分布與已知分布之間相似程度的一種檢驗(yàn)方法，通過顯著性水平來判別是否相似。

設(shè)累積頻率曲線的理論分布形式為Fe（x），樣本的累積頻率F1（x）=k/n。其中：n為總觀測數(shù)，k為不大于x的次數(shù)。

3 原州區(qū)降水隨機(jī)模擬分析

3.1 原州區(qū)降水量理論分布函數(shù)

根據(jù)我國水文計(jì)算的相關(guān)規(guī)范，我國主要流域的降水量可假定服從P-Ⅲ型分布或?qū)?shù)正態(tài)分布[15]。通過計(jì)算機(jī)軟件用原始數(shù)據(jù)畫頻率曲線發(fā)現(xiàn)，當(dāng)變差系數(shù)CV=0.26、偏態(tài)系數(shù)CS=0.52時(shí)，曲線模擬效果很好。又由P-Ⅲ型分布特性可知，當(dāng)CS=2CV時(shí)，參數(shù)a0=0，此時(shí)的P-Ⅲ型分布就是伽瑪分布，故筆者采用a0=0的P-Ⅲ型分布，即伽瑪分布進(jìn)行隨機(jī)模擬。因?yàn)楣P者主要研究降水的隨機(jī)模擬，需要建立最為適合的模型，所以對兩種分布進(jìn)行對比，選出最貼近原始數(shù)據(jù)的分布，從而建立最優(yōu)模型。

3.1.1 伽瑪分布

3.1.2 對數(shù)正態(tài)分布

3.1.3 參數(shù)確定

通過原始降水?dāng)?shù)據(jù)可求得式（7）和式（8）的未知參數(shù)。由原始降水?dāng)?shù)據(jù)可知，原州區(qū)近60 a降水量的期望值為442.0 mm，方差為11 659.2 mm，從而求得對數(shù)正態(tài)分布參數(shù)μ為6.064、σ為0.234，伽瑪分布中的形狀參數(shù)α為14.80、尺度參數(shù)β為0.03。

3.1.4 對比分析

將原州區(qū)近60 a降水量的經(jīng)驗(yàn)分布分別與對數(shù)正態(tài)函數(shù)理論分布及伽瑪函數(shù)理論分布作對比，如圖1所示?？梢钥闯?，樣本數(shù)據(jù)與對數(shù)正態(tài)分布、伽瑪分布的擬合效果均較好，但無法進(jìn)一步判斷最適合的理論分布，故筆者應(yīng)用K-S檢驗(yàn)對兩者進(jìn)行概率分布判斷。

為了使結(jié)果更加精確，取置信度α=0.01，n=60，其K-S臨界值Dα（n）為0.206 7。當(dāng)假設(shè)數(shù)據(jù)服從對數(shù)正態(tài)分布時(shí)，最大偏差D（n）為0.094 1;假設(shè)數(shù)據(jù)服從伽瑪分布時(shí)，最大偏差D（n）為0.107 6。通過對比分布的最大偏差與K-S檢驗(yàn)臨界值，發(fā)現(xiàn)兩者均不能拒絕原假設(shè)，而且在99%的置信區(qū)間上符合對數(shù)正態(tài)分布及伽瑪分布。由此可知，原州區(qū)降水?dāng)?shù)據(jù)既服從對數(shù)正態(tài)分布，又服從P-Ⅲ型分布。此處選取一般水文規(guī)范所用到的P-Ⅲ型分布進(jìn)行隨機(jī)模擬及預(yù)測研究。

3.2 降水量蒙特卡洛隨機(jī)模擬

借助計(jì)算機(jī)編程，引用均勻分布模型，通過蒙特卡洛方法產(chǎn)生（0，1）區(qū)間的隨機(jī)數(shù)，將產(chǎn)生的結(jié)果視為隨機(jī)變量概率。為了使模擬更具可靠性，使用伽瑪分布的反函數(shù)，根據(jù)原始降水序列的均值和方差，產(chǎn)生1 000 000組隨機(jī)降水序列值。用平方誤差作為評判標(biāo)準(zhǔn)，選取與原始降水序列數(shù)據(jù)平方誤差最小的一組作為隨機(jī)模擬降水量結(jié)果，見圖2。

將原始降水序列依次與通過伽瑪分布函數(shù)隨機(jī)模擬出的降水序列進(jìn)行對比分析?？梢钥闯?，在伽瑪分布隨機(jī)模擬出的新降水序列中，整體波動(dòng)較大，部分極值點(diǎn)在時(shí)間趨勢上保持一致。部分重要指標(biāo)的對比見表1，伽瑪分布所模擬出來的均值、標(biāo)準(zhǔn)差等指標(biāo)與原始數(shù)據(jù)相差不大，模擬效果良好，進(jìn)一步證明了原始數(shù)據(jù)服從P-Ⅲ型分布的結(jié)論。

3.3 馬爾可夫降水預(yù)測驗(yàn)證

當(dāng)前對馬爾可夫預(yù)測準(zhǔn)確度的評定大都是通過原始數(shù)據(jù)進(jìn)行對比驗(yàn)證，沒有考慮未來數(shù)據(jù)驗(yàn)證的準(zhǔn)確性。而從隨機(jī)模擬結(jié)果得到的P-Ⅲ型分布就可以很好地解決這個(gè)問題。筆者先通過對原始數(shù)據(jù)的預(yù)測及驗(yàn)證判別馬爾可夫方法的可行性，之后再進(jìn)行5 a預(yù)測，并進(jìn)行K-S檢驗(yàn)，判別其是否依舊符合P-Ⅲ型分布，從而確定預(yù)測結(jié)果的準(zhǔn)確性。

3.3.1 分級標(biāo)準(zhǔn)及狀態(tài)確定

考慮到降水序列的連續(xù)性及降水特性，采取降水序列的3 a滑動(dòng)平均值進(jìn)行降水預(yù)測?？紤]水文現(xiàn)象的本身特性及序列數(shù)據(jù)的結(jié)構(gòu)合理性，將降水量序列分為5類，即將降水量劃分為5個(gè)區(qū)間，分別為枯水年、偏枯年、平水年、偏豐年、豐水年，采用基于歐氏距離的層次聚類法進(jìn)行分級，見表2，其中x為年降水量。

3.3.3 相關(guān)系數(shù)與權(quán)重的確定

3.3.4 模型檢驗(yàn)及分析

基于隨機(jī)模擬出的2007—2011年的年降水?dāng)?shù)據(jù)，采用加權(quán)馬爾可夫模型對2012年降水量進(jìn)行預(yù)測，結(jié)果見表5，2012年的預(yù)測值為狀態(tài)3，根據(jù)模糊集理論，級別特征值為2.824，2012年的降水量預(yù)測值為500.8 mm，相對誤差為13.9%，在中長期水文預(yù)報(bào)允許誤差20%之內(nèi)。2013年、2014年預(yù)測值分別為521.4、499.5 mm，相對誤差分別為6.5%、5.8%，都在允許誤差以內(nèi)。

通過以上對原始數(shù)據(jù)的驗(yàn)證，能在一定程度上證明加權(quán)馬爾可夫模型在降水預(yù)測應(yīng)用上的可行性。但是，其只能對未來降水進(jìn)行預(yù)測，并不能對未來降水預(yù)測的準(zhǔn)確性進(jìn)行檢驗(yàn)，要確保對未來預(yù)測的準(zhǔn)確性，還需要進(jìn)一步的驗(yàn)證。由于前文已經(jīng)得到寧南地區(qū)降水符合P-Ⅲ型分布，因此筆者通過對寧南山區(qū)降水進(jìn)行接下來3 a滑動(dòng)平均的預(yù)測，再通過K-S檢驗(yàn)，看是否滿足P-Ⅲ型分布，預(yù)測結(jié)果見表6。

此時(shí)，n=63，取置信度α=0.01，其K-S臨界值Dα（n）為0.201 8，假設(shè)數(shù)據(jù)服從P-Ⅲ型分布，最大偏差D（n）為0.080 6。通過對比分布的最大偏差與K-S檢驗(yàn)臨界值，發(fā)現(xiàn)不能拒絕原假設(shè)，而且在99%的置信區(qū)間上符合P-Ⅲ型分布?？梢姡瑢ξ磥? a的降水預(yù)測結(jié)果是可靠的。

4 結(jié) 論

（1）依據(jù)原州區(qū)1957—2016年實(shí)測降水資料，運(yùn)用蒙特卡洛方法進(jìn)行隨機(jī)模擬，通過K-S檢驗(yàn)，得出原州區(qū)降水符合P-Ⅲ型分布。

（2）通過對降水?dāng)?shù)據(jù)進(jìn)行3 a滑動(dòng)平均處理，采用基于歐氏距離的層次聚類法進(jìn)行狀態(tài)劃分，確定了枯水年、偏枯年、正常年、偏豐年和豐水年5個(gè)狀態(tài)，建立了滑動(dòng)平均-加權(quán)馬爾可夫鏈預(yù)測模型。

（3）通過對2012年、2013年、2014年滑動(dòng)平均降水量的預(yù)測，得出馬爾可夫降水預(yù)測模型可以應(yīng)用于原州區(qū)降水預(yù)測。

（4）對未來5 a降水進(jìn)行預(yù)測，預(yù)測值分別為508.5、520.8，554.9、451.0、466.6 mm，通過K-S檢驗(yàn)進(jìn)行驗(yàn)證，證明了預(yù)測的準(zhǔn)確性。

參考文獻(xiàn)：

[1] 魏艷華，王丙參.概率論與數(shù)理統(tǒng)計(jì)[M].成都：西南交通大學(xué)出版社，2013：264-269.

[2] 丁晶，鄧育仁.隨機(jī)水文學(xué)[M].成都：成都科技大學(xué)出版社，1988：9-10.

[3] 劉新立.隨機(jī)過程與隨機(jī)模擬在水災(zāi)風(fēng)險(xiǎn)管理中的應(yīng)用研究[J].經(jīng)濟(jì)科學(xué)，2003（1）：114-119.

[4] 尹正杰，袁宏源，崔遠(yuǎn)來，等.灌區(qū)灌溉需水量的隨機(jī)模擬[J].中國農(nóng)村水利水電，2001（11）：19-22.

[5] 溫季，郭樹龍，盧聞航.作物灌溉隨機(jī)模擬技術(shù)研究[J].人民黃河，2004，26（5）：39-41.

[6] ZHANG J L， LI Y P， HUANG G H， et al. Assessment of Parameter Uncertainty in Hydrological Model Using a Markov-Chain-Monte-Carlo-Based Multilevel-Factorial-Analysis Method[J].Journal of Hydrology， 2016， 538（7）：471-486.

[7] KNIGHTON J， WHITE E， LENNON E， et al. Development of Probability Distributions for Urban Hydrologic Model Parameters and a Monte Carlo Analysis of Model Sensitivity[J].Hydrological Processes， 2015， 28（19）：5131-5139.

[8] 劉悅憶，趙建世，黃躍飛，等.基于蒙特卡洛模擬的水質(zhì)概率預(yù)報(bào)模型[J].水利學(xué)報(bào)，2015，46（1）：51-57.

[9] 李娟，張維江，馬軼.滑動(dòng)平均-馬爾可夫模型在降水預(yù)測中的應(yīng)用[J].水土保持研究，2005，12（6）：200-202，209.

[10] 王艷，毛明志，范晶，等.最優(yōu)分割法確定的加權(quán)馬爾可夫鏈在降雨量預(yù)測中的應(yīng)用[J].統(tǒng)計(jì)與決策，2009，25（11）：17-18.

[11] 岳遙，李天宏.基于模糊集理論的馬爾可夫模型在水質(zhì)定量預(yù)測中的應(yīng)用[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，2011，19（2）：231-242.

[12] 李亞斌，徐盼盼，錢會，等.加權(quán)馬爾可夫鏈在銅川地區(qū)降水量預(yù)測中的應(yīng)用[J].灌溉排水學(xué)報(bào)，2017，36（5）：96-102.

[13] GARCA-ALONSO C R， ARENAS-ARROYO E， PREZ-ALCAL G M. A Macro-Economic Model to Forecast Remittances Based on Monte-Carlo Simulation and Artificial Intelligence[J].Expert Systems with Applications， 2012， 39（9）：7929-7937.

[14] 夏樂天.馬爾可夫鏈預(yù)測方法及其在水文序列中的應(yīng)用研究[D].南京：河海大學(xué)，2005：11-12.

[15] 夏樂天，朱元甡.馬爾可夫鏈預(yù)測方法的統(tǒng)計(jì)試驗(yàn)研究[J].水利學(xué)報(bào)，2007，38（增刊1）：372-378.

【責(zé)任編輯 張 帥】