城市降雨频率分析数学方法与计算软件——改进

城市降雨频率分析数学方法与计算软件4.1 概述目前，对城市降雨资料进行频率分析的方法仍然使用P-Ⅲ分布模型和Weibull分布模型。

在推求城市暴雨强度公式时，可以根据降雨基础资料画出直方图，直方图的横坐标为暴雨强度，纵坐标为某一强度区间范围内出现的暴雨场数与总场数的比值，由直方图可以直观地看出基础资料分布密度函数的分布形态是呈乙形还是铃形，从而确定选用哪一种分布模型。

如果直方图呈铃形分布，就可以采用P-Ⅲ分布模型，如果直方图呈乙形分布，就可以采用Weibull分布模型。

我们在编制雨样频率分析计算软件的过程中，对P-Ⅲ分布模型和Weibull分布模型的各种计算方法进行了深入的研究，其中P-Ⅲ分布模型统计方法中采用了矩法、适线法、极大似然法、遗传算法和加速遗传算法。

4.2 城市降雨频率P-Ⅲ分布模型与分析计算软件4.2.1 P-Ⅲ分布模型1895年，英国生物学家K.Pearson在统计分析了大量随机现象后，发现概率密度曲线均为类似于铃形的曲线（见下图），这种曲线有两个特点：y（1）只有一个众数在众数处曲线的斜率等于零。

若把纵坐标移到均值处，即当dx -=0=dxdy（2）曲线的两端或一端以横轴为渐进线，即当0=y0=dxdy根据这两点，K.Pearson 建立了概率密度曲线微分方程式2210)(xb x b b y d x dxdy +++=根据上式，所得出的Pearson 曲线簇共有13条曲线，P-Ⅲ型曲线是其中的一种，其02=b ，微分方程的形式为xb b y d x dxdy 10)(++=经移轴、参数代换、分离变量积分、整理得dx d aeax y y -+=)1(0式中 0y —— 众值处纵坐标；a —— 系列起点到众值的距离； d —— 均值到众值的距离；经移轴、参数代换，利用概率分布特性最后得出的皮尔逊-Ⅲ型曲线方程式的另一种形式)(1)()()(αβγγαγβ----Γ=x ex x f （4-1）式中 γ —— 代换参数 d a /1=-γ；β —— 代换参数 d /1=β；α —— 系列起点到坐标原点的距离；e —— 自然对数的底；)(γΓ —— 伽玛函数，)(γΓ =⎰∞--01dx exxγ24sC=γsv C C x 2=β)21(sv C C x -=α其中x 、v C 、s C 分别为均值、离差系数、偏差系数，它们可以由矩法、适线法、极大似然法和遗传算法等方法求得。

若已知P-Ⅲ型频率密度曲线，则理论累积频率曲线可以求得：=≥)(p x x P dx ex dx x f ppx x x ⎰⎰∞---∞-Γ=)(1)()()(αβγγαγβ（4-2）令)(αβ-=x t ，则 )(αβ-=p p x t ，β/dx dt =。

因为)(γΓ与x 无关，所以dt etP px t⎰∞--Γ=)(1γγ （4-3）即 dt etP pt t⎰Γ-=--01)(1γγ （4-4）dt etP pt t⎰Γ='--01)(γγ （4-5）公式(4-5)是著名的不完全伽玛函数，可由数学方法算出，所以理论频率P P '-=1可以求得。

4.2.2 P-Ⅲ模型参数估计与计算软件─矩法1．矩法基本原理矩法是用样本矩估计总体矩，并通过矩和参数之间的关系，来估计频率曲线参数的一种方法。

一阶原点矩的计算公式就是均值x ，均方差σ的计算式为二阶中心矩开方，偏态系数s C 计算式中的分子则为三阶中心矩。

设x 为由总体随机变数X 抽出的一个样本，当∞→n 时，样本的分布)(x Fn 将趋近于总体的分布)(x F ，因此，样本的各阶矩也必然相应地趋向于总体的各阶矩。

参数估计的矩法就是根据这个性质，采用样本分布函数的相应数字特征去估计总体的同一种数字特征。

具体来说，对于常用的三个数字特征x 、v C 、s C ，采用如下的估计量。

对于Ex ，采用样本分布的数学期望n xx i∑=作为估计量。

当有了一个具体样本x 时，就假定∑=ix nx 1（4-6） 对于x C v /σ=，采用样本的离差系数作为估计量，即令xx x nx S C i v 2)(1-∑==（4-7）因此常用1)1(1)(122--∑=--∑=n k n x x xC i i v （4-8）作为v C 的估计量。

其中S 为样本的方差，x x k i i =称为模比系数。

对于s C ，采用样本的偏差系数作为估计量。

由于2/32333μμσμ==s C （4-9）所以习惯上常令33)2)(1()1(vi s C n n k n C ---∑=（4-10）式中2μ、3μ分别为二阶中心矩和三阶中心矩。

矩法估计的优点是使用简便，且计算公式与总体分布无关。

但是该方法具有如下两个缺点：(1) 除对x 的估计外，一般EU 要小于待估参数0μ，不满足不偏性。

(2) 有效率较低，即抽样误差较大，特别是对三阶矩s C 更是如此，而对于包含两参数的分布，矩法确实是一种很好的方法。

2．参数估计与计算步骤使用矩法对降雨资料进行频率分析的具体步骤如下：第一步：对降雨的原始数据进行处理，得到符合要求的工作数据，并对其依次读取。

第二步：计算降雨强度的均值x ，离差系数v C 和偏差系数s C∑=ix nx 11)1(2--∑=n k C i v33)2)(1()1(vi s C n n k n C ---∑=其中i x 为各历时的降雨强度， x x k i i =为模比系数。

第三步：计算公式参数α、β和γ)21(sv C C x -=αsv C C x 2=β24sC=γ第四步：将参数α、β和γ代入公式(4-2)，得到理论分布曲线，进而完成频率分析。

3. 程序流程图矩法对降雨资料进行频率分析的程序流程图如下所示：开 始 ↓ 输入工作数据↓ 计算x ,v C ,s C↓ 计算α,β,γ↓ 数据存盘 ↓结 束4.2.3 P-Ⅲ模型参数估计与计算软件─适线法1．适线法基本原理设所研究的随机变数为X ，通过观测有了一个样本，将这个样本按大小重新排列，得到：n m x x x x ≥≥≥≥≥ 21，然后采用下式计算其经验频率：1+=n m P n （4-11）对每一个样本i x ，求出一个经验频率i n P ,，并设n P x P x i n i i n i ~1,),,(,,=，这n 对数值点在几率格纸上，称为经验点据。

当根据某种理由把分布函数的形式确定以后，假设一组参数值1321),,(u u u ，就可通过积分{}dx u u u x f x X P x G P x),,;()(321⎰∞=≥== （4-12）求出)(x G P =的关系。

这个积分一般都有制成的表格可查，不需进行积分。

把)(x G P =的关系也画在上述几率格纸上，这样画出的曲线通常称为理论曲线。

检查这条曲线与经验点据的配合情况，如配合不好，就再假设另外一组参数1321),,(u u u ，再重画一条理论曲线。

如果仍然符合不好，就再重新假设，直到曲线与点据拟合良好为止。

适线法本来是一种图解方法，为了采用计算机计算，必须用解析形式表达出来。

采用某种经验频率公式，对每一个i x ，都规定一个i n p ,，与之对应，将各对数值ni P x i n i ,,2,1),(, =，逐个点绘在几率格纸上，共得n 个经验点据。

先由（4-6）式算出x ，由于数学期望的这个估计具有不偏性而且方差最小，因此不需要通过适线法来确定。

再假定一组v C ，s C ，就可以求出相应于这一组v C ，s C 的频率曲线：),,(s v C C x x G P ；=（4-13）求出这条曲线与经验点据的纵向离差i x ∆及横向离差i n P ,∆，此处)(,i n i i P x x x -=∆（4-14）其中)(,i n p x 为满足下式的x 值。

i n p x s v P dx C C x x f i n ,)(,),,(=⎰∞； （4-15）)(i i i x P P P -=∆n i ,,2,1 = （4-16）其中 ),(),,()(s v i xs v i C C x x G dx C C x x f x P ；；==⎰∞（4-17）i P ∆相当于在均匀分格的格纸上点据与曲线的横向离差，也可采用在海森几率格纸上点据与曲线的横向离差，记为i n P ,∆。

（1）适线准则作为拟合曲线与经验点据间离差的量度，可以有多种不同的指标，如21i x ∆∑=Φ，i x ∆∑=Φ2，213ix ∆∑=Φ，24i p ∆∑=Φ，i p ∆∑=Φ5，26)(i p ∆∑=Φ，m p ∆∑=Φ7，ip x ∆∆∑=Φ8等。

取不同的指标达到最小来求解v C ，s C 就构成了不同的适线准则。

如果根据min 21=∆∑=Φi x 来求解v C ，s C ，就称为1Φ准则，以此类推。

一般对于同一样本，不同的准则，所得的结果也不相同。

根据经验，采用8Φ准则不如1Φ准则，对以上准则所做的计算比较发现1Φ准则和4Φ准则比较好。

（2）经验频率公式适线法中使用的经验频率公式有很多种类，对于简单的样本来说，经验频率公式的一般形式为：cn b m P n +-=（4-18）为了使上式所确定的n 个n P 对于%50=n P 能成为对称的，可以证明b 、c 必须满足12=+c b ，于是(4-18)式可以写成bn b m P n 21-+-=（4-19）又由b 的三种取值情况，使（4-19）式有三种变形，当b =0时成为威布尔公式。

一种适线准则结合一种绘点位置公式，就可以构成一种适线方法。

于是上述的七种准则、三种公式就可以组合成二十一种不同的适线方法。

2．参数估计与计算步骤适线法对降雨资料进行频率分析的步骤如下：第一步：选取合适的降雨资料，对每一降雨强度计算其经验频率n P1+=n m P n其中m 为降雨资料作为工作数据由大到小排序后的序号，n 为工作数据的总数。

第二步：计算降雨工作数据的均值x 和离差系数v C∑=ix nx 11)1(2--∑=n k C i v第三步：在倍数因子范围内取对应的s C 值。

由公式（4-1）中1)(--γαx 可知min 0x ≤≤α，又因为)21(s v C C x -=α，所以)(22min x x x C C s v -≤=Φ≤，在此范围内按步长为1依次取整数Φ值，并计算v s C C Φ=的值。

第四步：根据所得x ，v C 和s C 值，计算对应的α，β，γ值，并计算出理论频率P 。

第五步：由上可得72≤Φ≤，用黄金分割法搜索Φ值，以21i x ∆∑=Φ和24i p ∆∑=Φ为适线准则。

具体为，先找出2~7之间的两个黄金分割点1m ，2m ，且382.0)27(1⨯-=m ,618.0)27(2⨯-=m 。