基于McmC的新安江模型參數(shù)不確定性分析及應(yīng)用

畢詩(shī)詠 馬 靜

（吉林省水文水資源局延邊分局 吉林 延吉 133001）

1 研究背景

水文預(yù)報(bào)精度在很大程度上取決于流域水文模型參數(shù)的優(yōu)選和檢驗(yàn)，現(xiàn)有的水文預(yù)報(bào)的模型都是概念性水文模型，水文模型的參數(shù)大多具有明確的物理意義，從理論上應(yīng)該是可以通過(guò)實(shí)測(cè)來(lái)確定的。但是由于缺乏實(shí)測(cè)值，現(xiàn)有的模型主要是對(duì)實(shí)際水文物理過(guò)程的概化，很多待定的模型參數(shù)并不是反映單一過(guò)程或者單一影響因素的結(jié)果，因此不大可能通過(guò)實(shí)測(cè)來(lái)進(jìn)行確定，只能依靠系統(tǒng)分析方法，采用最優(yōu)化技術(shù)進(jìn)行求解[1]。人們常采用的試算法以及客觀優(yōu)選法主要依賴于人的經(jīng)驗(yàn)，因此，所選的參數(shù)常因人而異，且費(fèi)時(shí)費(fèi)力。在求解結(jié)構(gòu)復(fù)雜且參數(shù)比較多的模型時(shí)，采用傳統(tǒng)的最優(yōu)化方法會(huì)出現(xiàn)局部最優(yōu)的問(wèn)題，難以確保得到的參數(shù)全局最優(yōu)，而且求解步驟比較繁瑣。因此，流域水文模型的參數(shù)優(yōu)選工作一直都是一個(gè)大家關(guān)注的重點(diǎn)以及難點(diǎn)問(wèn)題[2]。本文將結(jié)合Mcm C方法對(duì)新安江模型參數(shù)不確定性進(jìn)行分析。

2 McmC方法

MCMC（Markov Chain Monte Carlo）方法源自于物理學(xué)研究，MCMC 方法是根據(jù)Bayesian推斷為中心的后驗(yàn)分布來(lái)模擬隨機(jī)樣本的一種動(dòng)態(tài)的蒙特卡洛方法[3]。

MCMC 方法的核心思想是構(gòu)造一個(gè)概率轉(zhuǎn)移矩陣，建立一個(gè)以分布π（x）為平穩(wěn)分布的Markov 鏈來(lái)得到π（x）的樣本，產(chǎn)生若干條獨(dú)立并行的Markov 鏈來(lái)探索模型參數(shù)空間，通過(guò)不斷更新樣本信息而使Markov 鏈?zhǔn)諗坑诟吒怕拭芏葏^(qū)，也就是Bayesian 方法中的最大后驗(yàn)估計(jì)[4]。Mcm C方法中的Metropolis-Hastings取樣方法能夠更有效的探索未知參數(shù)空間，因此本文采用Mcm C方法推求新安江模型參數(shù)的后驗(yàn)分布及其預(yù)測(cè)區(qū)間。

2.1 常用的MCMC采樣方法

MCMC 方法基于Bayesian 理論，通過(guò)建立以π（x）為平衡分布的馬爾可夫鏈，并對(duì)其分布進(jìn)行采樣，通過(guò)不斷的更新樣本信息使得馬爾可夫鏈能夠充分的搜索模型的參數(shù)空間，最終收斂到高概率密度區(qū)域。因此，MCMC 方法是對(duì)理想的Bayesian 推斷過(guò)程的一種近似[5]。常用的MCMC 采樣方法有：Gibbs 采樣方法和Metropolis-Hastings 方法。

2.2 McmC方法的實(shí)施步驟

實(shí)施步驟概括為如下三步：

Step1：在狀態(tài)空間D上建立一個(gè)以 π（x）為平穩(wěn)分布，轉(zhuǎn)移核為(p·，·)的Markov鏈；

Step2：由D中某一點(diǎn)X（0）出發(fā)，用Step1中的Markov鏈產(chǎn)生點(diǎn)序列X（1），L，X（n）；

Step3：對(duì)某個(gè)m和足夠大的n，用式（1）估計(jì)任一函數(shù)f(x)：

圖1 McmC的水文模型參數(shù)不確定性分析流程

表1 參數(shù)的取值范圍

表2 參數(shù)后驗(yàn)分布

3 實(shí)例分析

將Mcm C方法應(yīng)用到新安江模型，采用東洋河流域1963年～1988年10場(chǎng)歷史洪水進(jìn)行研究，應(yīng)用Mcm C方法對(duì)由模型參數(shù)引起的模型輸出不確定性進(jìn)行估計(jì)。

本例中Mcm C方法采用的參數(shù)為：新安江模型參數(shù)個(gè)數(shù)n=15；種群規(guī)模pops=400；線程數(shù)p=8；初始進(jìn)化代數(shù)n=10；最大循環(huán)次數(shù)L=50000。

3.1 確定參數(shù)的初始范圍和先驗(yàn)分布函數(shù)

一般情況下，參數(shù)的先驗(yàn)分布形式不容易確定，因此常常采用均勻采樣或者對(duì)數(shù)采樣等方式進(jìn)行代替，本文采用均勻分布。

3.2 基于MCMC的不確定性分析

通過(guò)對(duì)敏感性分析選定的新安江模型的SM、KG、KSS、KKSS、GS、WDM這六個(gè)參數(shù)進(jìn)行分析，由Mcm C方法獲取30000個(gè)參數(shù)組，設(shè)置目標(biāo)函數(shù)值的閾值為0.7，選擇高于閥值的所有參數(shù)組，并設(shè)置為有效參數(shù)組，然后對(duì)有效參數(shù)組分別進(jìn)行計(jì)算，得到計(jì)算流量值，按照似然值的大小進(jìn)行排序，并設(shè)定置信水平為90%（置信度上限95%，下限為5%）的模型計(jì)算不確定性區(qū)間。

表2中的最小值以及最大值與表1中的參數(shù)搜索范圍基本一致，說(shuō)明采用Mcm C方法進(jìn)行不確定性研究，參數(shù)空間得到了有效的搜索。

圖2給出由MCMC方法得到的東洋河流域新安江模型的參數(shù)SM、KG、KSS、KKSS、CS、WDM的邊緣分布，每個(gè)刻度對(duì)應(yīng)的縱坐標(biāo)值表示該參數(shù)值落在該刻度值和前一刻度值之間的概率。

圖2給出了由MCMC方法抽樣得到的東洋河流域新安江模型的六個(gè)參數(shù)的邊緣分布，每刻度對(duì)應(yīng)的縱坐標(biāo)值表示該參數(shù)值落在該刻度與前一刻度值之間的概率。由圖2可以看出，SM、KG、KSS、KKSS、CS、WDM這六個(gè)參數(shù)都有明顯的分布規(guī)律。第二章中的優(yōu)選值SM=17.517、KG=0.2817、KSS=0.4134、KKSS=0.6157、CS=0.2001、WDM=78.5762大致落在圖2中各對(duì)應(yīng)參數(shù)概率最大的區(qū)域，由此可以看出Mcm C方法得到的結(jié)果和D E算法是一致的。MCMC方法不僅可以給出模型各參數(shù)的后驗(yàn)分布，而且可以識(shí)別參數(shù)空間內(nèi)概率最大的區(qū)域。

圖2 新安江模型各參數(shù)邊緣后驗(yàn)分布直方圖

采用Mcm C方法達(dá)到收斂以后抽取的新安江模型的SM、KG、KSS、KKSS、CS、WDM六個(gè)參數(shù)的2000個(gè)參數(shù)組樣本對(duì)東洋河流域的19820801號(hào)、19830907號(hào)、19880814號(hào)歷史洪水分別進(jìn)行模擬，可生成三場(chǎng)洪水的2000個(gè)模擬流量數(shù)據(jù)，再根據(jù)求得的模擬數(shù)據(jù)求出流量分布函數(shù)，并且求出該分布的5%和95%分位數(shù)作為其水文預(yù)報(bào)的90%不確定性置信區(qū)間。圖3到圖5給出了東洋河流域新安江模型模擬三場(chǎng)洪水的90%置信度的不確定性區(qū)間。

圖3 東洋河流域19820801場(chǎng)次洪水的不確定性范圍

圖4 東洋河流域19830907場(chǎng)次洪水的不確定性范圍

圖5 東洋河流域19880814場(chǎng)次洪水的不確定性范圍

圖3給出了198208號(hào)洪水實(shí)測(cè)流量過(guò)程、水文預(yù)報(bào)90%不確定性區(qū)間。圖4給出了198309號(hào)洪水實(shí)測(cè)流量過(guò)程、水文預(yù)報(bào)90%不確定性區(qū)間。圖5給出了198808號(hào)洪水實(shí)測(cè)流量過(guò)程、水文預(yù)報(bào)90%不確定性區(qū)間。由圖3、圖4和圖5可以看出，不確定性范圍在不同時(shí)間段是不同的，隨流量而變，在高流量區(qū)較大，在低流量區(qū)較小。模擬得到的流量過(guò)程線的上、下邊界并不能完全包含實(shí)測(cè)流量過(guò)程線，總有一些實(shí)測(cè)流量落在90%的置信區(qū)間之外，并不能完全模擬流域的流量過(guò)程。但不確定性范圍基本包含了實(shí)測(cè)流量過(guò)程的多數(shù)，說(shuō)明新安江模型在東洋河流域的應(yīng)用是可行的，模型結(jié)構(gòu)本身誤差所產(chǎn)生的不確定性影響在可接受范圍內(nèi)。

4 結(jié)論

針對(duì)新安江模型參數(shù)不確定性分析過(guò)程中出現(xiàn)的參數(shù)多、收斂慢、計(jì)算負(fù)擔(dān)重等問(wèn)題，提出了Mcm C方法并對(duì)東洋河流域新安江模型參數(shù)不確定性進(jìn)行分析。

（1）Mcm C方法可以用于新安江模型參數(shù)后驗(yàn)分布的抽樣，對(duì)參數(shù)的不確定性進(jìn)行評(píng)價(jià)。根據(jù)Mcm C方法得到的后驗(yàn)分布抽樣，得出新安江模型預(yù)報(bào)值的經(jīng)驗(yàn)分布，由此可以得到新安江模型參數(shù)分布的平均值、眾數(shù)、標(biāo)準(zhǔn)差等估計(jì)值，為水文模型不確定性的定量評(píng)價(jià)以及水文預(yù)報(bào)提供有效途徑。

（2）Mcm C方法使樣本多樣性得以豐富和保持，減少了收斂于局部最優(yōu)區(qū)域的可能性，從而提高了求解質(zhì)量和計(jì)算速度。抽樣結(jié)果對(duì)東洋河流域進(jìn)行不確定性預(yù)報(bào)，并給出預(yù)報(bào)90%不確定性區(qū)間，結(jié)果表明該區(qū)間較好的反映了與模型參數(shù)不確定性相關(guān)的洪水預(yù)報(bào)不確定性范圍。

[1]丁晶,鄧育仁.隨機(jī)水文學(xué)[M].成都科技大學(xué)出版社,1988.

[2]N.K.Goel.隨機(jī)水文學(xué)[M].黃河水利出版社,2001.

[3]邢貞相,芮孝芳,崔海燕,等.基于A M-Mcm C算法的貝葉斯概率洪水預(yù)報(bào)模型[J].水利學(xué)報(bào),2007，38(12)：1500-1506.

[4]王建平,程聲通,賈海峰.基于Mcm C法的水質(zhì)模型參數(shù)不確定性研究[J].環(huán)境科學(xué),2006,27(1)：24-30.

[5]黎光明,張敏強(qiáng).先驗(yàn)信息對(duì)Mcm C方法估計(jì)概化理論方差分量變異量的影響[J].統(tǒng)計(jì)與決策，2012，7：27-29.