年降水量概率分布的贝叶斯推断及其数值求解

孔 强

(石家庄市滹沱河生态工程运维服务中心,河北 石家庄 050000)

在水文频率计算方法上,一直面临的一个问题是:基于不同样本容量水文系列的频率计算,会得出不同的结果,反映在实际中,就是很多已建水利工程,当初的设计水文指标值与现实情势不符,须进行重新核算。分析原因,较为普遍的观点主要有:ⓐ自然因素和人类活动使得水文事件产生条件发生变化,水文资料存在非平稳、非一致性;ⓑ水文实测样本容量偏小,不能满足所需的充分样本长度;ⓒ选用线型参数本身存在误差,影响了计算结果精度。 为此,众多学者已进行了多项研究,如由子系列组成混合分布分模型[1]以解决非一致条件的水文频率问题,建立变参数PDS/GP 模型[2]以提高参数稳定性;鲁帆等[3]提出了贝叶斯MCMC 洪水频率分析方法,用以提高分布参数的估计精度;戈立婷等[4]和黄华平等[5]分析了GEV 分布所需的充分样本长度和应用历史调查洪水资料保证充分样本容量的方法。 这些研究对促进水文频率计算发展进行了有益的探索。

本次研究通过计算年降水量样本参数的估计区间,推断水文总体正态分布的可能分布域,采用贝叶斯推断,以数值法估计给定水文样本的后验概率组合,计算一定概率下的年降水量值,旨在将正态分布作为水文事件的可能总体分布。

1 水文样本序列的正态性检验

1.1 频率推断与抽样

水文频率计算的基本假定是将水文事件作为随机事件,认为样本序列是独立随机抽取自某一总体分布,其变化规律服从概率分布律[1]。 通常采用适线法,即选用对水文序列经验点据拟合较好的线型,如极值Ⅰ型和Ⅱ型、广义极值分布(GEV)、对数正态分布(LN)、皮尔逊Ⅲ型(P-Ⅲ)等,计算频率设计值。 对于频率线型存在的误差,有一种观点认为:水文频率推断的非平稳、非一致性假设已不满足,一些“极端”水文事件可能产生系列突变[6]。 从大尺度历史和未来趋势来说,这些已观测的“极端”事件,更是一种自然规律[7]。因此,以几十个实测数据得出的水文系列长周期突变结论,还有待跨学科深入探讨。

从统计抽样视角,由样本估计的确定参数值,只是所有可能的参数取值的其中之一,反映在推断结果上就是,即便对样本拟合再好的线型,也是所有可能线型之一,只是概率不同,而取用基于样本的某个确定的线型来代表整体分布,存在巨大理论缺陷,即便样本拟合度再完美,也是不一定能正确反映总体分布。

1.2 水文事件总体分布正态性检验

水文事件的产生过程受多重原因影响,包括自然环境、人类生产、星系运动等,对于水文事件总体分布是正态分布还是其他某种分布尚无定论,若将水文事件看成许多随机因素叠加的结果,定性分析或可认为水文事件应是正态分布或近似正态分布,这也是对水文事件总体分布的合理猜测。 对于样本的正态性检验方法有很多,下面采用W 检验(Shapiro-Wilk 检)与D检验(D’Agostino 检验)法[8],对地区年降水量等水文样本进行正态性检验。 水文样本容量(3≤n≤50)时,采用W 检验:

设X1,X2,…,Xn是从总体X中抽取的容量为n的样本,若X服从正态分布,应满足

统计量W为

对不同的n,系数αi(i=1,…,n)可查W 检验统计量系数表得到;Wα是在给定的显著性水平α下使P(W≤Wα)=α的临界值。

假设H0:X服从正态分布,水文样本容量(n＞50)时, 采用D检验,检验统计量D为

在H0为真时,E(D) ≈0.28209479,令

D’Agostino 用随机模拟法得到了Y的分位数表,在给定显著性水平α后,用统计量Y进行检验的拒绝域为:,采用D 检验,由于样本Y未落入拒绝域中,故在α= 0.05 时,可认为服从正态分布。

给定显著性水平α=0.05 时,由式(1)计算了资料中43 个城市样本容量为22 的W值(n= 22,Wα=0.911),有41 个通过了正态检验,5 个水文站或入库径流量中,有3 个通过了正态检验;16 个城市样本容量为36 的W值(n= 36,Wα= 0.935),有15 个通过了正态检验,5 个水文站或入库径流量中,有2 个通过了正态检验;由式(2)计算了7 个城市样本容量为60 的D值(n=60,=1.13),全都通过了正态检验。 可见降水量无论是W检验还是D检验,都较好地满足正态性。 而水文站或入库径流量正态性较差,可能是受人为影响因素较大,计算结果见表1。

表1 中,南昌、武汉、长沙、海口、重庆、成都、贵阳、昆明、拉萨、西安、兰州、西宁、乌鲁木齐、天津、石家庄、太原、呼和浩特、沈阳、哈尔滨、上海、杭州、合肥、福州、南京、银川、北京、济南、郑州、广州、南宁30个城市的1998—2020 年降水量数据均取自1998—2020 年的《中国统计年鉴》;榆林[9]、庐江[10]、永吉[11]、张家口[12]、石嘴山[13]、宝鸡[14]、南京[15]、银川[16]、咸阳[17]、密云水库入库径流[18]、滦河滦县站年径流[19]、北京[20]、济南[21]、郑州[21]、广州[23]、南宁[24]、韶关[25]、陇南[26]、温州[27]、洪家渡水文站构皮滩水文站[28]20 个地区的年降水量或径流量均来自公开发表论文或技术资料。

2 推断方法

2.1 均值和方差的估计

对正态分布的均值和方差的估计,采用极大似然估计量[8],即总体均值μ=X,总体方差样本的极大似然值不代表所有可能取值,需要对、S2可能取值作区间估计。

2.1.1 总体均值μ的区间估计

由于总体方差σ2未知,可用σ2的无偏估计S2代替,则得随机变量为含有未知参数μ的枢轴量[8],服从T～tn-1分布,由于t分布的概率密度函数是单峰对称的,查t分布表可得自由度为n-1 的t分布的上侧分位数,使得

故μ的置信度为1-α的置信区间为

2.1.2 总体方差σ的区间估计2

样本方差S2是σ2的无偏估计,得随机变量

由于χ2分布的概率密度函数图形是单峰不对称的,由χ2分布的上侧分位数定义,有

故σ2的置信度为1-α的置信区间为

则标准差σ的置信度为1-α的置信区间为

取显著性水平α=0.05,则置信区间为1-α=0.95。 对于正态分布的样本均值置信区间取t分布,方差取χ2分布,分别计算样本均值和方差置信区间。

2.2 可能的正态总体

正态分布样本的均值和方差S2相互独立[8],即值与S2值无关。 若水文样本来自均值、方差估计区间的边界值范围内的正态分布,则水文样本分布必然会收敛于均值、方差估计区间内的正态分布域。 因此样本均值置信区间内的所有可能值与方差置信区间内所有可能值,组合成的是有界正态分布簇,其边缘为置信区间上下限的4 个组合:(均值下限;标准差下限)、(均值下限;标准差上限)、(均值上限;标准差下限)、(均值上限;标准差上限),1 个最大似然正态分布是(样本均值;样本标准差)。

2.3 后验概率分布

2.3.1 贝叶斯公式

贝叶斯公式可以将先验知识和样本信息结合起来进行推断,基本公式为

式中:θ代表参数;x代表样本;f(θ|x) 为后验概率密度函数;π(θ) 为先验分布;f(x|θ) 为似然函数。

2.3.2 先验概率分布

前述已证实,区域年降水量等非人为控制水文样本总体趋向于正态分布,因此,可以水文事件正态分布作为贝叶斯推断的先验分布:

式中:N(x;μ,σ) 为正态分布概率密度函数;x为水文序列值;μ为均值;σ为标准差。

2.3.3 似然函数

由已知样本拟合得到的分布曲线,如P-Ⅲ分布、GEV 分布等都是基于已知样本的可能分布,其分布是总体分布的似然分布,其分布函数为似然函数,若采用常用的P-Ⅲ型曲线作为拟合线型[29],则式(1)中似然函数f(x|θ) 可写为

式中:Г(·) 为gamma 函数;α0,β,α分别为分布的位置、尺度、形状参数,且有α≥α0、α ＞0、β ＞0,此3 个参数与总体参数Ex、Cv、Cs(期望值、变差系数、偏态系数)关系如下:

2.3.4 后验概率分布

若先验分布π(θ) 和似然函数f(x|θ) 已知,则可以求出二者的乘积函数g(θ|x):

对乘积函数g(θ|x) 在Θ 上求解积分值k,则可由式(8)得出后验密度函数。 贝叶斯公式中,分母为两个函数的积函数在积分域Θ 上,对参数θ的积分,理论上,先验分布π(θ) 和似然函数f(x|θ) 的积分域是[- ∞,+ ∞],一般采用马尔可夫链蒙特卡罗(MCMC)方法求解贝叶斯公式。 而对于特定水文事件,定义积分域[- ∞,+ ∞] 是理论极限,在可遇见的未来,或存在实际的降水或径流的上限和下限,故可取特定的计算区间[a,b], 下限为可能最小值,上限为可能最大值,选取计算步长,以两个函数值乘积函数为分子目标函数值,则式(5)可写为

3 后验概率分布的数值求解

以南宁市1961—2020 年降水量样本均值(μ=1302.14)、标准差(σ=227.241)生成正态分布随机数为算例,对南宁市年降水量后验概率进行数值求解。

3.1 确定正态分布组合

由式(3)和式(4)计算南宁市年降水量均值、方差的95%置信水平的取值范围边界值,可组合出5 个正态总体分布,见表2。

表2 南宁市降水量均值、方差置信区间边界值正态分布组合

作为可能的正态先验密度函数簇边界概率,由式(6)可计算5 个正态总体分布先验概率密度曲线,见图1。

图1 组合1 ～5 正态可能先验概率密度曲线

3.2 计算似然函数

对样本数据以P-Ⅲ曲线进行拟合,由式(7)求得P-Ⅲ拟合曲线密度函数(三参数gamma 分布)[29],作为贝叶斯公式的似然函数,P-Ⅲ拟合频率与样本计算频率拟合情况见图2。

图2 样本频率与P-Ⅲ拟合曲线对比

样本频率与P-Ⅲ拟合曲线拟合度及误差计算见表3。

表3 南宁市年降水量P-Ⅲ拟合曲线拟误差计算表

表4 傅里叶多项式参数

3.3 计算乘积函数

正态先验概率密度函数式(6)与似然函数式(7)相乘的解析解求解较为复杂,可采用数值法,积分域Θ取[400,2500],计算步长设置为1mm,计算似然函数值与正态密度函数值的乘积值,以matlab 曲线拟合工具中8 参数傅里叶多项式对乘积值曲线进行拟合,得出似然函数与组合1 ～5 可能正态先验密度函数乘积的函数拟合公式:

拟合公式的R-square 和Adjusted R-square 值均为1,SSE 和RSME 值见表5。 以组合1 正态先验值与P-Ⅲ似然函数数值计算,先验概率密度曲线、似然函数曲线、后验概率密度曲线见图3,组合2 ～5 计算过程相同。

图3 组合1 概率密度数值计算对比

表5 组合1 ～5 拟合优度及k 值

3.4 计算k 值

在给定积分域Θ[400,2500]与乘积函数g(θ|x)曲线顶点值max(g(θ|x)) 围成的矩形范围内,采用MC 法计算5 个乘积函数g(θ|x) 在Θ 上的积分值,随机取10 万个数值(1mm 单位),计算得g(θ|x) 所围图形面积,由式(8)可求得k值,见表5。

将k值代入式(9)可得5 个后验概率密度函数。组合1 ～5 后验概率密度曲线见图4。

图4 组合1 ～5 后验概率密度曲线对比

3.5 计算给定频率的降水量值

给定频率值,带入组合1 ～5 的后验概率密度函数,即可求出0.95 置信水平的降水量取值区间。 对图4 中5 个后验概率密度取累计值,可得相应的概率质量函数曲线,见图5。 对于给定的某个频率值,其降水量取值为组合1 ～5 对应的取值范围:频率0.10 =[1394,1492],频率0.05=[1451,1562],频率0.01=[1567,1693]。 由于正态分布总体均值和方差置信区间的对称估计,本方法数值计算结果相较于常用的水文频率设计值的估计量抽样分布的标准差计算方法[9]和含义都不同。

图5 组合1 ～5 后验概率质量曲线对比

4 结 论

4.1 推断结果是1 个取值域

本方法将对水文事件的先验知识融入了频率推断过程,通过正态先验分布参数边缘值组合,利用贝叶斯方法,灵活运用数值模拟计算和MC 法,使水文频率计算符合对事物的认知,在水文频率推断方向上是正确的。 对给定样本的水文序列频率计算结果是1 个取值范围,这个结果含义与教材中设计值的估计量抽样分布的标准差计算方法和含义是完全不同的,在一定程度上解决了常用的拟合曲线可能很适合样本数据,但预测仍欠佳的问题,推断结果所确定的值域有较好的概率理论依据,丰富了水文频率计算方法。

4.2 对样本的拟合是误差主要来源

对样本数据的拟合有多种方法,每种方法会得出不同的拟合函数,故而拟合函数对计算结果产生的误差是不可避免的,而且是误差的主要来源。 相较而言,计算过程中对组合1 ～5 的先验和后验的数值拟合函数,以及MC 方法在设定积分域上对k值的模拟计算,所产生的误差会随着增加计算量而大幅减小。

4.3 采用数值法对贝叶斯公式直接求解是可行的

理论上,水文事件概率分布的两端为[- ∞,+ ∞],根据可遇见的实际情况,将不能发生情况的概率认定为0,对贝叶斯公式中分母积分域设置为固定值,采用MC 方法计算定积分,进而求解k值。 采用数值法计算贝叶斯公式,直接模拟了贝叶斯公式分子乘积函数,为贝叶斯公式在水文计算中的应用提供了新方法。