通过迭代重复上述步骤。这样对足够大的k值，序列日41,02，…可以认为是后验分布(2)的抽 样。和独立样本一样，可以用这个序列来估计一些数字特征，比如后验均值。这种方法可以完全不 考虑所选的g(当然必须服从某些正则条件)。从极限意义上讲，式(7)给出的拒绝步骤确保了模拟序 列的平稳分布具有所需的边缘分布。 4实例分析 丹江口水库是汉江中下游防洪和水资源开发的关键水利枢纽工程，也是南水北调中线工程的水 源工程。枢纽控制流域面积9.52万km2,约占汉江集水面积的60%，多年平均入库径流量约363亿m。 根据该水库建库后1969一2008年的入库径流资料，分别选取设计时段(1、3、5、7d)的年最大洪量系 列和年最大洪峰流量系列作为洪水频率计算的样本数据。 4.1GEV分布参数与设计洪水的估计GEV分布的分布函数为： ≠0 F(x)= (8) x-) exp-exp =0 其中：μ、σ、分别为GEV分布模型的位置参数、尺度参数和形状参数，并满足4∈R,σ>0，∈R, 1+(x-u)/o>0。 为计算方便，洪量系列的单位统一为万m/s,以下类同。分别对丹江口水库设计时段(1、3、5、7d) 的年最大洪量值系列和年最大洪峰流量系列建立GEV分布模型。令中=Iog0,以保证σ>0。选择先验 密度函数为π(u,中，)=π(u)π（中）π()。其中π（日、（日和π，（均值为零，方差分别为4、6和E 的正态密度函数，参数、中和之间是独立的。选择足够大的方差使密度函数相当平坦，取u=。=10,三 100,这样就确定了先验密度。 根据Metropolis--Hastings方法原理，对向量(w,中，)的每个分量分别应用式(7)，且每个建议分 布g选为各自坐标轴上的随机游动，即u=4+e4,中=中+e。,=(+e,其中，ee、e。、e是均值 为0，方差分别为w4、"。、”的正态变量。"。、"b、0,的数值根据试验确定。 限于篇幅，仅给出年最大1日洪量GEV分布模型每个参数20O00次迭代产生的MCMC序列，见 图2。其中，令σ，=，得到尺度参数σ后验分布的模拟样本。在除去前面不稳定的模拟值后，其余 比较稳定的模拟值就认为是要求的后验分布的观测。 从图2可以看出，左边的序列大约在3000次迭代后趋于稳定，右边的序列大约在1000次迭代后 趋于稳定。但左边图中显示各参数的收敛值与右边图中显示各参数的收敛值非常接近。这说明参数 0a、06、w,的选择并不影响计算的模型参数，只是影响算法的效率。 将向量(4,0，)的每个模拟值代入式(9) 4--(g-p以门5≠0 X1-p= (9) u-log-log(1-p),=0 就得到相应的1p年重现期的设计洪量值的后验分布样本。图3给出了选取参数w.=0.01、心。= 025、心=0.25的情况下，年最大1日洪量GEV分布模型的参数u、o、和百年一遇设计洪量xo9的后 验密度估计。其中，百年一遇设计洪量的单位为万m/s。 表1列出了贝叶斯估计方法和极大似然估计方法估计的不同极值系列GEV模型的分布参数值及 百年一遇设计洪量值。 -945-通过迭代重复上述步骤。这样对足够大的 k 值，序列θk+1，θk+2，…可以认为是后验分布（2）的抽 样。和独立样本一样，可以用这个序列来估计一些数字特征，比如后验均值。这种方法可以完全不 考虑所选的 q（当然必须服从某些正则条件）。从极限意义上讲，式（7）给出的拒绝步骤确保了模拟序 列的平稳分布具有所需的边缘分布。 4 实例分析 丹江口水库是汉江中下游防洪和水资源开发的关键水利枢纽工程，也是南水北调中线工程的水 源工程。枢纽控制流域面积 9.52 万 km2 ，约占汉江集水面积的 60%，多年平均入库径流量约 363 亿 m3 。 根据该水库建库后 1969—2008 年的入库径流资料，分别选取设计时段（1、3、5、7d）的年最大洪量系 列和年最大洪峰流量系列作为洪水频率计算的样本数据。 4.1 GEV 分布参数与设计洪水的估计 GEV 分布的分布函数为： F ( x ) = ì í î ï ï ï ï ï ï ï ï exp ì í î ï ï ü ý þ ï ï - é ë êê ù û úú 1 + ξ ( x - μ ) σ -1 ξ ξ ≠ 0 exp ì í î ü ý þ -exp é ë êê ù û - úú ( x - μ ) σ ξ = 0 （8） 其中：μ、σ、ξ分别为 GEV分布模型的位置参数、尺度参数和形状参数，并满足 μ ∈ R，σ > 0，ξ ∈ R ， 1 + ξ ( x - μ ) σ > 0。 为计算方便，洪量系列的单位统一为万 m3 /s，以下类同。分别对丹江口水库设计时段（1、3、5、7d） 的年最大洪量值系列和年最大洪峰流量系列建立 GEV 分布模型。令ϕ = logσ ，以保证σ > 0。选择先验 密度函数为π ( μ，ϕ，ξ ) = πμ ( μ )πϕ (ϕ )πξ (ξ )。其中πμ (⋅)、πϕ (⋅)和πξ (⋅)均值为零，方差分别为 vμ、vϕ和 vξ 的正态密度函数，参数μ、ϕ和ξ之间是独立的。选择足够大的方差使密度函数相当平坦，取vμ=vϕ=104 ，vξ = 100，这样就确定了先验密度。 根据 Metropolis-Hastings 方法原理，对向量（μ，ϕ，ξ）的每个分量分别应用式（7），且每个建议分 布 q 选为各自坐标轴上的随机游动，即 μ* = μ + εμ，ϕ* = ϕ + εϕ，ξ * =ξ + εξ，其中，εμ、εϕ、εξ是均值 为 0，方差分别为wμ、wϕ、wξ的正态变量。wμ、wϕ、wξ的数值根据试验确定。 限于篇幅，仅给出年最大 1 日洪量 GEV 分布模型每个参数 20 000 次迭代产生的 MCMC 序列，见 图 2。其中，令σi = e ϕi ，得到尺度参数σ后验分布的模拟样本。在除去前面不稳定的模拟值后，其余 比较稳定的模拟值就认为是要求的后验分布的观测。 从图 2 可以看出，左边的序列大约在 3 000 次迭代后趋于稳定，右边的序列大约在 1000 次迭代后 趋于稳定。但左边图中显示各参数的收敛值与右边图中显示各参数的收敛值非常接近。这说明参数 wμ、wϕ、wξ的选择并不影响计算的模型参数，只是影响算法的效率。 将向量（μi ，σi ，ξi ）的每个模拟值代入式（9） x1 - p = ì í î ï ï μ - σ ξ é ë ù û 1 - {-log (1 - p )} -ξ ，ξ ≠ 0 μ - σlog {-log (1 - p )}， ξ = 0 （9） 就得到相应的 1/p 年重现期的设计洪量值的后验分布样本。图 3 给出了选取参数wμ=0.01、wϕ = 0.25、wξ=0.25 的情况下，年最大 1 日洪量 GEV 分布模型的参数μ、σ、ξ和百年一遇设计洪量 x0.99的后 验密度估计。其中，百年一遇设计洪量的单位为万 m3 /s。 表 1 列出了贝叶斯估计方法和极大似然估计方法估计的不同极值系列 GEV 模型的分布参数值及 百年一遇设计洪量值。 — 945 —