http://www.paper.edu.cn
基于贝叶斯理论的TOPMODEL参数不确定性分析1
梁忠民， 李彬权，余钟波，华家鹏，刘金涛
河海大学水文水资源与水利工程科学国家重点实验室，南京（210098）
E-mail: libinquan@gmail.com
摘 要：应用贝叶斯理论探讨了TOPMODEL模型参数不确定性问题。首先通过Monte Carlo途径确定了模型中的敏感参数，结果表明，模型中的参数SZM、RV较为敏感，它们的微小改变都将影响模拟结果；其次用Markov chain Monte Carlo(MCMC)算法得到其后验分布的抽样，进而获得流域出口断面流量过程的抽样分布，并据此对模型参数的不确定性及其对预报结果的影响进行评价。提供了参数不确定性分析的流程，并通过模型预报结果的抽样分布，构造了90%的置信预报区间。
关键字：贝叶斯理论；MCMC；不确定性分析；参数敏感性分析；TOPMODEL；置信区间
水文模型参数反映的是流域水文特征，是一些原则上可以实测的物理量[1]，但实际中难以做到，所以往往都需要通过实测资料率定一套“最佳”模型参数[2] [3]。由于测量误差、资料的代表性不足及率定过程中目标函数、计算方法的局限性，一般只能得到局部最优解或多组不同的参数值（即所谓的“异参同效”），所以，实际应用中模型参数不可避免地存在着不确定性。模型参数不确定性的存在，必将导致对水文过程模拟和预测结果的不确定性。目前，关于这种不确定性的定量描述及其对水文模拟和预报不确定性的影响评价，在国内外已成为研究热点。
本文在贝叶斯理论框架下研究水文模型参数的不确定性。通过Monte Carlo途径，分析模型参数变差系数与模型模拟的确定性系数或Nash-Sutcliffe（NS）系数之间的关系，设定确定性系数变化阈值并挑选模型的敏感参数。对敏感参数，根据贝叶斯公式构建其后验分布，并采用MCMC抽样技术获得“集束”样本，运行水文模型可以得到对应的模型预报的抽样分布，根据模型参数的后验分布及预报值的抽样分布即可对参数的不确定性及其对预报的影响进行定量分析。作为实例，本文选择TOPMODEL参数的不确定性问题进行研究，以浙江密赛流域为典型流域，提供了具体的分析流程和结论。
1. 模型参数敏感性分析
大量研究表明，在水文模型的识别中，灵敏度较低的参数的取值对目标函数值影响很小，而灵敏度较高的参数的影响较大。灵敏度较低的参数对模型预报不确定性的影响也较小，所以，为减少计算量，本文只对灵敏度较高的模型参数进行不确定性分析，因此，首先需进行参数的敏感性分析。具体步骤如下：
ˆ(i=1,2,L,n)，n代① 通过任何传统的模型识别方法可以得到模型参数θi的率定值θi
表模型参数个数；
② 对第i个参数θi，将其作为随机变量并假设其概率分布，本文假设其服从对数正态分布：
(lnθi−a)22σ2
f(θi)=
1
12πθiσe
− （1）
基金项目：教育部博士点基金资助项目（20070294018） 、国家自然科学基金资助项目（50779013）、长江学者和创新团队发展计划资助（IRT0717）
- 1 -
http://www.paper.edu.cn
ˆ；σ=Cv∗a，Cv为θ的变差系数。从该分布中随机抽取θ，代入到水文 式中，a=lnθiii
ˆ(j≠i)），得到预报流量过程，计算确定性系数NS(i)； 模型中（其它参数采用θj
③ 改变Cv值，重复步骤②，则可以建立第i个参数的NS(i)~Cv关系； ④ 重复步骤②和③,得到所有参数的NS(i)~Cv关系；
⑤ 分析各参数的确定性系数NS(i)随Cv值的变化情况，给出衡量参数是否敏感的阀值λ。若NS(i)的变化区间大于λ，则认为参数θi为敏感性参数。
2. 模型参数不确定性分析的贝叶斯方法
贝叶斯理论的基本原则是利用数据获得变量的条件分布，即变量的后验概率密度函数，公式为：
h(θ/x)=
L(x/θ)g(θ)
∫L(x/θ)g(θ)dθ (2)
式中θ和x分别代表随机变量和与之相关的统计数据，h(θ/x)为后验概率密度函数，
L(x/θ)为似然函数，g(θ)为先验分布。
对本文研究的水文模型，θ为模型的敏感性参数，x为模型的输出变量，即流量Q；并将进行水文模型模拟或预报所需要的其它要素，如实测降雨、流量等次洪资料等记为w。则公式（2）改写为：
h(θ/Q,w)=
L(Q/θ,w)g(θ/w)
∫L(Q/θ,w)g(θ)dθ2−1
j
∝L(Q/θ,w)g(θ) （3）
似然度函数L(Q/θ,w)的度量有多种方法[4]，本文取如下形式： L(Q/θ,w)=(
∏σj=1
N
)
1m2
， σ=∑(Qcal,k−Qobs,k) (4)
mk=1
2j
2
式中，Qcal,j、Qobs,j分别为计算与实测的流量值，σj为第k场洪水误差的方差，m为第i场洪水的总时段数，N为次洪个数。
20世纪90年代，研究人员将马尔科夫链蒙特卡罗法（Markov Chain Monte Carlo, MCMC）[5]引入到参数不确定性研究中，用于待估参数的贝叶斯分布采样，以估计参数的后验分布。由于这种方法的应用，使得随机模拟在很多领域的计算中，显示出巨大的优越性。本文将应用MCMC法进行模型参数的不确定性研究。
3. AM-MCMC方法
本文采用自适应的Metropolis算法(Adaptive Metropolis, AM)估计参数的后验分布。AM是一种改进的MCMC采样器[6]，相比传统的M-H与Gibbs采样，AM不再需要事先确定变量的转移分布，而是决定于初始抽样的协方差。将转移分布定义为参数空间的正态分布形式，以此为依据对参数空间运行随机抽样，其初始协方差可根据先验分布确定。在抽样过程中根据马尔可夫链的历史抽样信息自适应地调整转移密度（即协方差矩阵），且可并行运算，从而大大提高算法的收敛速度。这样，第i步参数的推荐分布为均值θi ，协方差Ci的多元正态分布。协方差迭代计算公式如(5)所示，在初始i0次迭代中,协方差矩阵Ci取固定值
- 2 -
http://www.paper.edu.cn
C0，之后自适应更新。
−−TSd−−Ti−1
Ci+(iθi−1θi−1−(i+1)θiθi+θiθiT+εId) (5) Ci+1=ii
式中，ε为一个较小的数(在这里取ε=10-5)，以确保Ci不成为奇异矩阵；Sd为一个比例因子，
以确保接受率在一个合适范围内(在这里取Sd=(2.4)/d)；θi−1依赖于参数的空间维度d，
和θi为前i−1和i次迭代参数的均值。Id 为d维单位矩阵。
MCMC 研究的一个重要任务是判断并行采样序列是否收敛到后验分布。理论上一个各向同性的采样器在t→∞时一定收敛，而在实际应用中并非如此。在实际应用中总是希望找到最小的t值。为了实现该目标，本文采用Geman 等[7]在1992年提出了一种定量的比例缩小得分(scale reduction score)R用以诊断收敛性。计算方法如下:
−
2
−
R=
g−1q+1B
+ (6) gq∗gW
式中，g为每一参数采样序列的迭代次数；q为用于评价的序列数；B/g为q个序列的平均值的方差；W为q个序列的方差的平均值。计算每一参数的比例缩小得分R，若接近于1表示参数收敛到了后验分布上。Gelman和Rubin建议将R<1.2作为多序列抽样算法收敛判断条件。针对单序列是否稳定，则可采用平均值法和方差法，考察迭代过程中的平均值和方差是否稳定。
4. 算例
（1）流域及资料情况
选择TOPMODEL模型在浙江省开化县密赛流域的应用作为分析实例。该流域面积为700.6km2，年降雨量在1500~2000mm之间，属湿润地区。共选用了从1982年到1988年的9场次洪资料。 （2）参数敏感性分析
在本研究中，计算时段Δt为1小时，并将马氏京根法参数Ke取为Δt，所以需要进行敏感性分析的参数共7个，表1给出了这些参数的取值范围。
表1 TOPMODEL模型参数取值范围
参数 最小值 最大值
SZM Ln(T0) Td SRmax SR0 RV CHV 0 1 0 0.001 0 4000 100 0.05 10 30 0.1 0.05 8000 10000
在对模型参数敏感性分析过程中，由于根带土壤饱和缺水量的初值SR0 在不同场次洪水中取值不同，而其它参数在同一单元流域上各场次洪的取值可以认为是相同的，因此SR0的参数敏感性需要单独分析。
以确定性系数NS为目标函数，敏感性阀值λ取0.2，分析结果如图1所示。从图1(a)可以看出参数SZM和RV取值的微小改变对模拟结果影响较大，而无论参数T0、Td和CHV的取值如何变化，目标函数NS始终处于0.9～0.95之间，同样地，参数SRmax变化幅度
- 3 -
http://www.paper.edu.cn
也小于敏感性阈值λ；图1(b)使用了9场次洪对参数SR0进行敏感性分析，从图中可以看出
SR0并不敏感。分析结果表明，土壤下渗率呈指数衰减的速率参数SZM和地表坡面汇流的
有效速度RV为最敏感的2个参数。
（a）
SZM等6个参数NS
~Cv关系图
（b）参数SR0的NS
~Cv关系图
图1 模型参数的NS~Cv关系图
（3） MCMC算法配置
对SZM和RV，采用AM算法对其后验分布进行估计，其中公式（3）的先验分布都取为均匀分布。AM算法的配置为：参数初始协方差矩阵C0为对角矩阵，方差取参数取值范围(如表1所示)的1/20；初始迭代次数i0=1000。算法平行运行3次，每次采样5000个样本。
以参数SZM的分析结果为例，图2、图3和图4分别为单一采样序列的参数值、参数平均值和方差的抽样采样变化过程图。从图可看出，当i>500后，参数SZM的平均值和方差基本达到稳定，因此单一序列是收敛的。
，R利用公式(6)计算的R演化过程如图5所示。由图可知，在迭代初期（i<200）变化剧烈；当i>200后，R趋于稳定，并接近于1.0，说明参数SZM和RV的MCMC采样序列均能稳定收敛到其参数的后验分布，且算法全局收敛。综合考虑上述单序列和多序列参数的收敛情况，将每一序列的前500次舍去，这样3次平行试验共采集了13500个样本用于参数后验分布的统计分析。
- 4 -
http://www.paper.edu.cn
图2 参数SZM值抽样变化过程 图3 参数SZM均值变化过程
图4 参数SZM方差变化过程 图5 收敛指标值的变化过程
（4）不确定性分析结果
根据上述MCMC抽样得到参数SZM和RV的后验分布，见图6～图7。由图可知，RV比SZM具有更大的不确定性。分别取后验分布的众数作为SZM和RV的点估值，并与传统的Ronsenbrock 优化结果[8]一同列于表2。由表2可知，SZM的贝叶斯点估值与优化值较为接近，而RV却相差较大。由此说明在参数的识别中，具有较小不确定性的参数更可能得到合理的率定结果；相反，若参数不确定较大，则更可能得到具有相同效果的参数率定值。
图6 参数SZM后验分布直方图 图7 参数RV后验分布直方图
表2 MCMC法和Ronsenbrock优化结果对比
参数名称
后验分布的众数
Ronsenbrock
0.02 5000
相对误差(%)
SZM 0.019 RV 6300
5 26
将SZM和RV的MCMC抽样值输入TOPMODEL模型中（其它输入变量及模型参数取自文献[8]），得到预报变量的抽样分布，根据此分布可以构造不同置信水平的预报置信区间；同时，采用SZM和RV的贝叶斯点估值也可以得到一条预报过程线。图8提供了贝叶斯点估值、优化参数的预报结果及90％置信度的不确定性范围，采用参数后验分布众数的贝叶斯点估值可以得到比较好的预报结果（图8中贝叶斯预报值），而且，可以提供置信区
- 5 -
http://www.paper.edu.cn
间预报。从图中可以看出，不确定性范围随流量而变，在高流量区较大，在低流量区较小。小部分时段模拟预报与实测流量有偏差，这可能是由多种因素共同作用所致，具体可能是由于19820619次洪的多峰特性造成的。除此外，大多数实测数据都包含于90%的置信区间内，表明在TOPMODEL模型中，考虑了敏感参数的不确定性后，就可以较有效地定量描述预报结果的不确定性。但是，预报置信区间没能包含所有时段的实测数据的事实也表明，模型并不能完全模拟流域出口断面的流量过程，说明仅通过对敏感性参数的不确定性分析并不能对整个预报过程的不确定性做出全面的评价。
图8 19820619次洪预报与实测流量过程对比图
5. 主要结论
本文利用贝叶斯理论研究了水文模型参数的不确定性及其对预报结果的影响，主要结论如下：
（1）本文采用确定性系数（NS系数）并结合随机抽样，建立NS与参数Cv的变化关系，以此确定模型的敏感参数，结果表明，模型中的参数SZM、RV较为敏感，它们的微小改变都将影响模拟结果。
（2）MCMC算法可以用于模型参数后验分布的抽样，根据抽样样本可以直观地构建各参数的后验分布，进而对参数的不确定性进行评价。根据后验分布抽样，可以得到水文模型预报的抽样分布，由此可以提供诸如参数众数对应的预报结果，也可以获得某一置信水平的区间预报，这为模型预报不确定性的定量评价及水文过程的概率预报提供了有效途径。
（3）尽管MCMC采样方法可以克服自动搜索、随机搜索、试错搜索等寻优方法在水文模型高维参数空间的缺陷，但是对于多参数的模型结构，参数的组合方式太多，常常需要几万或几十万的参数取样，因此MCMC法需要消耗大量的计算资源，一般的模拟洪水过程程序运行时间都较长，因此一定程度了MCMC法在取样上的优点。
（4）水文模型的应用包含着输入、模型参数和结构等诸多的不确定性因素，本文仅对参数的不确定性问题进行了初步探讨，更完整的结论需要进一步的深入研究。
- 6 -
http://www.paper.edu.cn
参考文献
[1] 赵人俊. 流域水文模拟[M]. 北京：水利电力出版社. 1984.
[2] 叶守泽，夏军. 水文系统识别（原理与方法）[M]. 北京：水利电力出版社. 19.
[3] Duan Q., Sorooshian, Gupa V.K. Effective and efficient global optimization for conceptual rainfall-runoff models[J]. Water Resources Research,1992, 28(4), 265～284.
[4] Beven K Binley A. The future of distributed models:model calibration and uncertainty prediction [J]. Hydrological processes,1992,6:279～298
[5] 龚光鲁，钱敏平.应用随机过程教程及其在算法与智能计算中的应用[M].北京：清华大学出版社，2003. [6] Haario H,Saksman E,Tamminem J. An adaptive Metropolis algorithm[J]. Bernoulli, 2001,7(2):223～242
[7] Gelman A,Rubin DB. Inference from iterative simulation using multiple sequences[J]. Statistics Science,1992,7(4):457～511
[8] 姚成. 基于栅格的分布式新安江模型构建与分析[D]，河海大学硕士学位论文，2007，pp47.
Bayesian Theory for Parameter Uncertainty Analysis in
TOPMODEL
Liang Zhongmin, Li Binquan, Yu Zhongbo, Hua Jiapeng, Liu Jintao
State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai
University, Nanjing (210098) Abstract
Bayesian theory is applied to the parameter uncertainty analysis of TOPMODEL in this paper. The sensitive parameters of TOPMODEL are determined firstly by the Monte Carlo Method. The results show that the parameters SZM and RV are sensitive, any a slight change of them will influence the simulated results; Moreover, the sampling distribution of flow process at the basin outlet can be obtained by the later distribution sampling, which can be calculated with the algorithm of Markov chain Monte Carlo (MCMC). And the paper analyzes the parameter uncertainty and its influence on forecast based on MCMC method. The paper also provides the analysis process of the parameter uncertainty, and constructs the 90% confidence interval with the sampling distribution of the model forecast results.
Keywords: Bayesian theory; MCMC; uncertainty analysis; parameter sensitive analysis; TOPMODEL; confidence interval
作者简介:梁忠民(1962-),男,辽宁省凤城人,教授、博导,水文水资源专业。
- 7 -
