基于McmC的新安江模型參數不確定性分析及應用

畢詩詠 馬 靜

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

1 研究背景

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

2 McmC方法

MCMC（Markov Chain Monte Carlo）方法源自于物理學研究，MCMC 方法是根據Bayesian推斷為中心的后驗分布來模擬隨機樣本的一種動態的蒙特卡洛方法[3]。

MCMC 方法的核心思想是構造一個概率轉移矩陣，建立一個以分布π（x）為平穩分布的Markov 鏈來得到π（x）的樣本，產生若干條獨立并行的Markov 鏈來探索模型參數空間，通過不斷更新樣本信息而使Markov 鏈收斂于高概率密度區，也就是Bayesian 方法中的最大后驗估計[4]。Mcm C方法中的Metropolis-Hastings取樣方法能夠更有效的探索未知參數空間，因此本文采用Mcm C方法推求新安江模型參數的后驗分布及其預測區間。

2.1 常用的MCMC采樣方法

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

2.2 McmC方法的實施步驟

實施步驟概括為如下三步：

Step1：在狀態空間D上建立一個以 π（x）為平穩分布，轉移核為(p·，·)的Markov鏈；

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

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

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

表1 參數的取值范圍

表2 參數后驗分布

3 實例分析

將Mcm C方法應用到新安江模型，采用東洋河流域1963年～1988年10場歷史洪水進行研究，應用Mcm C方法對由模型參數引起的模型輸出不確定性進行估計。

本例中Mcm C方法采用的參數為：新安江模型參數個數n=15；種群規模pops=400；線程數p=8；初始進化代數n=10；最大循環次數L=50000。

3.1 確定參數的初始范圍和先驗分布函數

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

3.2 基于MCMC的不確定性分析

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

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

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

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

圖2 新安江模型各參數邊緣后驗分布直方圖

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

圖3 東洋河流域19820801場次洪水的不確定性范圍

圖4 東洋河流域19830907場次洪水的不確定性范圍

圖5 東洋河流域19880814場次洪水的不確定性范圍

圖3給出了198208號洪水實測流量過程、水文預報90%不確定性區間。圖4給出了198309號洪水實測流量過程、水文預報90%不確定性區間。圖5給出了198808號洪水實測流量過程、水文預報90%不確定性區間。由圖3、圖4和圖5可以看出，不確定性范圍在不同時間段是不同的，隨流量而變，在高流量區較大，在低流量區較小。模擬得到的流量過程線的上、下邊界并不能完全包含實測流量過程線，總有一些實測流量落在90%的置信區間之外，并不能完全模擬流域的流量過程。但不確定性范圍基本包含了實測流量過程的多數，說明新安江模型在東洋河流域的應用是可行的，模型結構本身誤差所產生的不確定性影響在可接受范圍內。

4 結論

針對新安江模型參數不確定性分析過程中出現的參數多、收斂慢、計算負擔重等問題，提出了Mcm C方法并對東洋河流域新安江模型參數不確定性進行分析。

（1）Mcm C方法可以用于新安江模型參數后驗分布的抽樣，對參數的不確定性進行評價。根據Mcm C方法得到的后驗分布抽樣，得出新安江模型預報值的經驗分布，由此可以得到新安江模型參數分布的平均值、眾數、標準差等估計值，為水文模型不確定性的定量評價以及水文預報提供有效途徑。

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

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

[2]N.K.Goel.隨機水文學[M].黃河水利出版社,2001.

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

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

[5]黎光明,張敏強.先驗信息對Mcm C方法估計概化理論方差分量變異量的影響[J].統計與決策，2012，7：27-29.