融合马尔科夫链_蒙特卡洛算法的改进通用似然不确定性估计方法在流域水文模型中的应用

2009年4月水 利 学 报SHUILI XUEBAO 第40卷 第4期收稿日期:2008-04-23基金项目:国家自然科学基金重点项目(40730632);教育部新世纪优秀人才支持计划(NCET-05-0624);霍英东青年教师基金资助项目(101077)作者简介:卫晓婧(1984-),女,山西阳泉人,硕士生,主要从事水文水资源方面的研究。

E -mail:hellomuki@to 文章编号:0559-9350(2009)04-0464-10融合马尔科夫链-蒙特卡洛算法的改进通用似然不确定性估计方法在流域水文模型中的应用卫晓婧,熊立华,万 民,刘 攀(武汉大学水资源与水电工程科学国家重点实验室,湖北武汉 430072)摘要:本文在Blasone 研究工作的基础上,进一步提出了基于马尔科夫链-蒙特卡洛算法的改进通用似然不确定性估计方法(Markov Chain -Monte Carlo based Modified Generalized Likelihood Uncertainty Esti mation,MMGLUE)。

该方法结合近年来被广泛用于推求参数后验分布的MC MC 方法,对基于Mon te Carlo 随机取样方法的传统GLUE 方法进行改进,并以预测区间性质最优为标准,对可行参数组阈值进行判断与选择,提出了衡量预测区间对称性的标准,并就预测区间性质与可行参数组个数的相关关系进行了探索。

在汉江玉带河流域的实例研究证明,MMGLUE 方法较传统的GLUE 方法能够推求出性质更为优良的预测区间,从而更真实合理地反映水文模型的不确定性。

关键词:MC MC;GLUE;MMGLUE;预测区间;覆盖率;区间宽度;区间对称性中图分类号:P333文献标识码:A1 研究背景近10年来,流域水文模型的不确定性研究逐渐成为当今水文界广泛研究的热点之一,各国的水文学家就此做了大量的工作[1]。

Beven [2-3]于1992年率先提出了流域水文模型/异参同效性0的观点,并针对流域水文模型的不确定性研究问题提出了通用似然不确定性估计(Generalized Likelihood Uncertainty Estimation,GLUE)方法。

该方法结合Monte Carlo 随机取样技术与Bayesian 框架,对由模型结构、参数冗余及相关性、输入输出误差等因素造成的不确定性进行综合分析。

GLUE 方法原理简单,易于操作,但由于其自身理论结构的缺陷,越来越多的研究者就GLUE 方法提出了质疑[4-5],即并非经典的Bayesian 方法、主观判断参数可行域阈值和推求的参数后验概率分布不具有显著的统计特征。

因此,基于不同假设的其他不确定性研究方法,如:基于经典Bayesian 理论的Ba RE(Bayesian Recursive Estimation)方法[6],基于全局卡尔曼滤波理论的EnKF(Ense mble Kalman Filter )方法[7],多目标方法如MOSCE M (Mult-i objective Shuffled Complex Evolution Metropolis)方法[8]等被用于估计模型的不确定性工作中。

然而,上述方法尽管理论结构相对复杂,应用效果与GLUE 方法相比却并没有明显的提高。

同时期另一种基于经典Bayesian 理论的马尔科夫链-蒙特卡洛(Markov Chain Monte Carlo,MC MC)方法也被广泛应用于推求参数后验分布的研究中。

特别是SCE M -UA (The Shuffled Complex E volutionMetropolis Algorithm)方法[9]能够有效地探索参数空间,使Markov Chain 能够朝着高概率密度区进化,从而推导出具有显著统计特征的水文模型参数的后验分布。

因此,Blasone [10]提出将两种方法结合起来,采用SCE M -UA 采样方法替代传统的GLUE 方法中的)464)Monte Carlo随机取样方法,并根据估计的预测区间的覆盖率来控制可行参数组个数的选择,对传统的GL UE方法进行改进。

本文在Blasone所做工作基础之上,进一步提出以预测区间性质最优为指标来控制可行参数组个数的选取。

2方法211贝叶斯统计推断贝叶斯学派是数理统计中的一个重要学派,其重要观点是[11]:任一未知参数H 都可以看作随机变量,因为任一未知量都有不确定性,因此可以用概率分布来描述。

人们根据先验信息对未知参数H的先验分布P(H),通过实验获得样本x1,x2,,,x n,对H的先验分布进行调整,调整的结果是H的后验分布h(H|x1,x2,,,x n)。

在这个过程中,人们的认识由P(H)调整到h(H|x1,x2,,, x n)。

贝叶斯方法中样本x1,x2,,,x n对H的条件密度p(x,x2,,,x n|H)就是经典方法中H已知时样本的联合密度。

一旦样本已知,就只有H在变化,把联合密度看成参数H的似然函数,用l(H|x1,x2,,, x n)来表示。

参数的后验分布表示为[12]h(H|x1,x2,,,x n)=P(H)l(H|x1,x2,,,x n)Q P(H)l(H|x1,x2,,,x n)d H(1)因为参数的后验分布综合了总体信息、样本信息和先验信息,因此对H的统计推断就应建立在后验分布的基础上。

贝叶斯假设:参数的无信息先验分布P(x)所在的取值范围内是/均匀0分布的。

根据最大熵原则,无信息如果意味着不确定性最大,那么,无信息的先验分布应是最大熵的相应分布,因为只有在分布是均匀时,熵才达到最大值,故本文中两种方法所采用的先验分布都是均匀分布。

经典统计学中处理点估计与区间估计方法不同,但在贝叶斯学派却是统一的。

对于贝叶斯统计中的区间估计,只要存在后验分布,就可以用相应分布的分位点给出参数H的置信区间,就模型参数不确定性分析而言,也就是预测区间。

问题就在于评判估计效果的标准。

本文中采用预测区间覆盖率、区间宽度、区间对称性作为最优后验分布判定的标准。

当后验分布已知时,对于给定的置信概率1-A可以求出很多置信区间。

由于参数H的最大后验区域估计集中了分布密度似然函数值取值尽可能最大的点,因此H的最大后验区间一定是在统一置信概率下区间宽度最狭窄的区间。

进而,推求参数的最大后验估计,成为不确定性分析方法研究的最终目的和手段。

212GL UE方法GLUE方法是目前最常用于不确定性估计的经验频率方法,它的原理与步骤如下:首先假设参数服从某一先验分布,通过Monte Carlo取样方法生成一定数目的可行参数组,然后利用流域降雨、蒸发、径流资料,计算各组参数值的对应的似然值。

那些与实际过程越接近的模型参数被认为具有越高的可信度与似然度。

最后主观选定一阈值,对似然度低于该阈值的参数组,令其相应的似然度为0;对高于该阈值的参数组,按照似然函数值由高到低排序,并标准化,再按照其似然值赋予相应的权重。

通过更新样本信息,从而取得参数的后验分布。

213MC MC方法MC MC是为了获得参数后验分布一系列后验量而发展起来的一种行之有效的计算方法,主要适用于多变量,非标准形式,且各变量间相互不独立时的分布模拟。

显然,MC MC方法非常适用于推求流域水文模型各参数的后验量。

Markov链具有如下特性:(1)无后效性:由随机变量序列组成的Markov链{X(0),X(1),X(2),,},在任一时刻t(t\0),序列中下一时刻处的X(t+1)由条件分布产生,它只依赖于时刻t处的当前状态而与时刻t之前的历史状态{X(0),X(1),,X(t-1)}无关;(2)各态遍历性:从不同的X(0)出发,链经过一段时间的迭代后,历经各种状态的Markov链最终收敛于平稳分布[13]。

MC MC方法的基本原理就是基于建立的平稳分布为P(x)的Markov链来获得P(x)的样本。

产生若干条独立并行的Markov链来探索模型参数空间,通过不断更新样本信息而使Markov链收敛于高概率)465)密度区,也就是Bayesian 方法中的最大后验估计。

MC MC 方法中的SCE M -UA 取样方法能够更有效的探索未知参数空间,因此本文采用该种方法推求实验模型参数的后验分布及其预测区间。

214 MC MC -based Modified GLUE (MMGLUE )方法 MMGLUE 方法,采用SCE M -UA 取样方法代替传统的GL UE 方法中的Monte Carlo 随机取样,并采用预测区间性质作为可行参数组数目x 的选取标准。

从而推求出有统计学意义的、性质优良的预测区间。

方法流程如图1所示,具体步骤如下:图1 MMGLUE 方法结构流程(1)选择实验模型以及相应的流域资料。

(2)确定似然函数,本文采用模型效率系数R 2作为似然函数:R 2=1-E M i =1(Q i -Q ^i )2P E M i =1(Q i - Q )2(2)式中:Q i 为实测径流量;Q ^i 为预测径流量; Q 为实测径流序列的均值;M 表示实测系列长度。

(3)根据先验分布随机产生s 个决策变量H t (t =1,2,,,s )。

(4)将s 个样本点划分为q 个区,每个区采用SCEM -UA[14-18]方法独立并行演化L 次以获得L #q 个样本点。

(5)将这些样本点掺混。

(6)将掺混后的样本点按照似然值由高到低排序。

(7)以预测区间对观测值覆盖率最合理、预测区面宽度最窄、区间对称性最优为标准,选取一合理初值x 带入模型进行预测区间性质检验,调整x 的值,使之达到既定标准。

覆盖率CR [19]及区间宽度IW 计算公式为CR =E M t =1J [Q obs,t ]P M ;IW =E M t =1(Q u p,t -Q low,t )P M (3)其中:J [Q obs,t ]=1,Q low,t <Q obs,t <Q up ,t0,其它式中:Q ob s,t 为时段t 实测径流量;Q low,t 为时段t 预测区间下界;Q up,t 为时段t 预测区间上界。

为了更好的反映预测区间的偏移程度,本文提出预测区间对称性IS 计算公式如下:)466)IS =E M t =1I [Q up,t ]P E M t =1I [Q low,t ](4)其中:I [Q low,t ]=1,Q ob s,t <Q low,t 0,其它; I [Q up,t ]=1,Q obs,t >Q up ,t 0,其它. 由式(4)可以看出,当IS =1时,区间对称,当0[IS <1时,区间较实际观测值偏高,IS >1时,区间偏低。