畢詩詠 馬 靜

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

1 研究背景

水文預(yù)報(bào)精度在很大程度上取決于流域水文模型參數(shù)的優(yōu)選和檢驗(yàn)，現(xiàn)有的水文預(yù)報(bào)的模型都是概念性水文模型，水文模型的參數(shù)大多具有明確的物理意義，從理論上應(yīng)該是可以通過實(shí)測(cè)來確定的。但是由于缺乏實(shí)測(cè)值，現(xiàn)有的模型主要是對(duì)實(shí)際水文物理過程的概化，很多待定的模型參數(shù)并不是反映單一過程或者單一影響因素的結(jié)果，因此不大可能通過實(shí)測(cè)來進(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)的問題，難以確保得到的參數(shù)全局最優(yōu)，而且求解步驟比較繁瑣。因此，流域水文模型的參數(shù)優(yōu)選工作一直都是一個(gè)大家關(guān)注的重點(diǎn)以及難點(diǎn)問題[2]。本文將結(jié)合Mcm C方法對(duì)新安江模型參數(shù)不確定性進(jìn)行分析。

2 McmC方法

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

MCMC 方法的核心思想是構(gòu)造一個(gè)概率轉(zhuǎn)移矩陣，建立一個(gè)以分布π（x）為平穩(wěn)分布的Markov 鏈來得到π（x）的樣本，產(chǎn)生若干條獨(dú)立并行的Markov 鏈來探索模型參數(shù)空間，通過不斷更新樣本信息而使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 理論，通過建立以π（x）為平衡分布的馬爾可夫鏈，并對(duì)其分布進(jìn)行采樣，通過不斷的更新樣本信息使得馬爾可夫鏈能夠充分的搜索模型的參數(shù)空間，最終收斂到高概率密度區(qū)域。因此，MCMC 方法是對(duì)理想的Bayesian 推斷過程的一種近似[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的不確定性分析

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

4 結(jié)論

針對(duì)新安江模型參數(shù)不確定性分析過程中出現(xiàn)的參數(shù)多、收斂慢、計(jì)算負(fù)擔(dā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.