确定含有重现期修正系数的暴雨强度公式的一种方法

摘要

本发明确定含有重现期修正系数的暴雨强度公式的一种方法，对于选配4个参数的暴雨强度公式i＝(A

说明书

一、技术领域

本发明属于市政工程技术领域。暴雨强度公式是城市排水工程设计中确定雨水设计流量 的基本依据之一，该公式的合理性与精度直接影响城市排水工程的规模。

二、背景技术

(一)常用的暴雨强度公式的形式

目前，我国城市常用的暴雨强度公式为^[1]：

单一重现期的分公式$i=\frac{A}{{(t+b)}^n}---\left(1\right)$

各重现期的总公式$i=\frac{A_1(1+\mathrm{ClgP})}{{(t+b)}^n}=\frac{A_1+C_1\mathrm{lgP}}{{(t+b)}^n}---\left(2\right)$

式中i为暴雨强度，简称雨强，mm/min；P为重现期，a；t为降雨历时，min。而A，A₁， C，b，n是与地方暴雨特性有关的参数，其中A为雨力，即t+b＝1时的雨强，mm/min；A₁为重现期为1a的设计降雨的雨力，mm/min；C为雨力变动参数，b为降雨历时修正参数， min；n为暴雨衰减指数；C₁＝A₁×C。

单一重现期的分公式对降雨历时t，雨强i的关系值(t，i)的拟合精度较高，但只能应用 于一个重现期的情况，故目前我国常用式(2)形式的公式。

由于公式(2)为纯非线性回归方程，可采用高斯-牛顿法或麦夸尔特法求解式(2)的拟合参 数。

(二)高斯-牛顿法与麦夸尔特法

1.高斯牛顿法

设一般的非线性方程为

y＝f(X，θ)+ε

式中f为一般函数；X可以是单个自变量，也可以是r个自变量X=(x₁，x₂，…，x_r)；θ为P维参 数向量，即θ＝(θ₁，θ₂，…，θp)′；ε为随机误差项，且∈N(0，σ²)。设y和X具有m组观察 值(x_1j，x_2j，…，x_rj；y_j)，j=1～m。求非线性模型参数的“最小二乘”估计，就是求θ的待估计值 使残差平方和

$Q\left(\theta \right)=\underset{j=1}{\overset{m}{\Sigma }}{ϵ}_i^2=\underset{j=1}{\overset{m}{\Sigma }}{[y_j-f(X_j,\theta )]}^2---\left(3\right)$

为最小。

为方便起见，用记号f_j(θ)代替f(X_j，θ)，高斯-牛顿法求解满足式(3)Q(θ)为最小的参数递 推公式，并写成矩阵形式为^[2，3]

θ_(k+1)＝θ_(k)+[J′(θ_(k))J(θ_(k))]^-1J′(θ_(k))[(y-f(θ_(k))]  (4)

$J\left({\theta }_{\left(k\right)}\right)={\left(\begin{array}{cccc}\frac{{\partial f}_1\left(\theta \right)}{{\partial \theta }_1}& \frac{{\partial f}_1\left(\theta \right)}{{\partial \theta }_2}& ...& \frac{{\partial f}_1\left(\theta \right)}{{\partial \theta }_p}\\ \frac{{\partial f}_2\left(\theta \right)}{{\partial \theta }_1}& \frac{{\partial f}_2\left(\theta \right)}{{\partial \theta }_2}& ...& \frac{{\partial f}_2\left(\theta \right)}{{\partial \theta }_p}\\ .& .& & .\\ .& .& ...& .\\ \frac{{\partial f}_m\left(\theta \right)}{{\partial \theta }_1}& \frac{{\partial f}_m\left(\theta \right)}{{\partial \theta }_2}& ...& \frac{{\partial f}_m\left(\theta \right)}{{\partial \theta }_p}\end{array}\right)}_{\theta ={\theta }_{\left(k\right)}}---\left(5\right)$

式中k为递推次数；y＝(y₁，y₂，…，y_j…，y_m)′，y_j为因变量的第j个观察值；f(θ_(k))＝[f₁(θ_(k))，f₂(θ _(k))，…，f_j(θ_(k))，…，f_m(θ_(k))]′，f_j(θ_(k))为由非线性方程及第k次迭代参数计算的第j个因变 量y的估计值，j＝1～m。

利用式(4)从参数初值θ₍₀₎开始，一步步递推下去，直到θ_(k)收敛稳定，即|θ_(k+1)-θ_(k)|小于 或等于预先指定的小正数δ(例如δ＝0.0001)，从而得到θ的估计值。

2.麦夸尔特法

麦夸尔特法，放宽了迭代计算时对参数初值的要求，其求解式(3)中暴雨公式参数θ的递 推公式，写成矩阵形式为^[4，5]

θ_(k+1)＝θ_(k)+[J′(θ_(k))J(θ_(k))+λI]^-1J′(θ_(k))[(y-f(θ_(k))]  (6)

式中：λ为阻尼因子；I为P阶单位阵；其他符号含义与式(4)中相应符号含义相同；偏导数 矩阵J(θ_(k))的计算式见式(5)。

关于阻尼因子λ的确定方法详见文献[5]。利用式(6)从参数初值θ₍₀₎开始，一步步递推下 去，直到θ_(k)收敛稳定，即|θ_(k+1)-θ_(k)|小于或等于预先指定的小正数δ(例如δ＝0.0001)，从而 得到θ的估计值。

(三)暴雨强度公式的推求路径及其拟合误差

推求暴雨强度公式中的参数，有两种不同路径。

路径一：利用理论频率计算所得的雨强选配公式，称其为传统方法。该法涉及两个环节 ^[1，6]。(1)理论频率曲线拟合：首先选配理论频率分布模型，例如，皮尔逊III型分布曲线、耿 贝尔分布曲线和指数分布曲线^[6]，并对各历时相应的实测雨强样本系列进行频率曲线拟合； 然后，根据各历时t的理论频率曲线，确定各计算重现期P的设计雨强i_p，进而得计算重现 期P、历时t、设计雨强i_p关系表。此环节存在理论频率曲线与实测雨强样本的拟合误差。(2) 暴雨公式参数拟合：根据计算重现期、历时、设计雨强关系表，选配暴雨强度公式。此环节 存在由暴雨强度公式拟合由理论频率曲线计算的设计雨强的误差，对于此误差， GB50014-2006《室外排水设计规范》^[6](2014版)附录A指出，对于采用年最大值法取样时， “计算重现期在2年～20年时，在一般强度的地方，平均绝对均方差不宜大于0.05mm/min(注： 原文为“平均绝对方差不宜大于0.05mm/min”。显然是错误的，均方差才与暴雨强度具有相 同量纲)。在较大强度的地方，平均相对均方差不宜大于5％”。

路径二：直接拟合法^[4]。文献[4]采用非线性回归方法，根据重现期，历时，实测雨强的 关系值(P_j，t_j；i_j)，j＝1～m，直接拟合暴雨强度公式参数，该法可实现拟合实测雨强样本的误 差平方和为最小，求得的暴雨强度公式参数为整体上的最佳参数，该法不仅减少了理论频率 分布模型的选择与拟合环节，从而避免了选择何种理论频率分布模型的纠结及避免推求暴雨 强度总公式参数的任意性，而且拟合实测雨强的均方差小于传统方法的结果，提高了拟合精 度。

需要指出，实测雨强样本是暴雨强度公式选配与拟合的基础，应与其拟合最佳。直接拟 合法的拟合误差则是拟合实测雨强样本的误差，其含义不同于规范^[6]中规定的误差。因此， 直接拟合法的误差不应该按规范^[6]中规定的误差上限来衡量。

三、发明内容

一些地区的计算成果表明，雨力A(即t+b＝1时的雨强)与重现期P的对数lgP不是直线关 系^[1]，即A≠A₁+C₁lgP，此种情况，即使采用非线性最小二乘法或其他最优化方法求解式(2) 的拟合参数，所得暴雨强度公式拟合某些重现期的雨强的绝对均方差或平均绝对均方差(即标 准差)仍超过0.05mm/min^[1，7，8]。有鉴于此，为了更密切的吻合当地多年统计的降雨强度子样 点的分布规律，提高拟合精度，不应再拘泥于4参数公式，而应采用式(7)5参数的公式

$i=\frac{A_1[1+\mathrm{Clg}(P+d)]}{{(t+b)}^n}=\frac{A_1+C_1\lg (P+d)}{{(t+b)}^n}---\left(7\right)$

式中：d为重现期修正系数，简称修正系数；其他符号含义与式(2)中相应符号含义相同。

文献[1]给出了确定式(7)中参数的手工算法(第932-939页)。到目前为止，尚无采用解析 法确定式(7)中参数的方法。

为了克服手工算法确定式(7)中参数效率低、所求参数具有一定的任意性、拟合精度不如 解析法高的缺陷，本发明提出考虑修正系数d的判别方法及含有该系数的暴雨强度公式(7)中 参数的确定方法。

本发明采用的技术方案是：对于4参数暴雨强度公式(2)拟合重现期P、历时t以及由理 论频率曲线确定的计算重现期相应的设计雨强i_p(或实测雨强i)关系值效果不理想的情况，首 先，根据求得的式(2)的参数n，b，利用最小二乘法计算各重现期的雨力A；并点绘A与lgP 关系点，若不为直线趋势，即A≠A₁+C₁lgP，则要考虑修正系数d，选配暴雨强度公式(7)； 然后，采用非线性回归方法确定包含修正系数d的暴雨强度公式中的5个参数。

本项发明的有益效果是，当A与lgP关系不为直线趋势时，采用含有修正系数d的5参 数暴雨强度公式更能客观地反映当地的暴雨特性；采用高斯-牛顿法或麦夸尔特法使所求暴雨 强度公式的参数为最佳拟合参数，且计算效率高；拟合精度显著提高。本项发明已利用我国 有关城市的暴雨特征值资料进行了验证，效果很好。

四、附图说明

下面结合附图和实施例对本发明进一步说明。

附图1为实例1的A与lgT关系图。

附图2为实例1修正系数初值d₍₀₎＝-0.6时A与lg(T-0.6)关系图。

附图3为实例2的A与lgT关系图。

附图4为实例2修正系数初值d₍₀₎＝-0.98时A与lg(T-0.98)关系图。

五、具体实施方式

(一)路径一——利用理论频率计算所得的雨强选配公式

实例1：文献[8]中昆明市松华坝雨量站，具有1964～2000年年最大10、20、30、45、 60、90、120、180min暴雨系列，将各时段样本系列采用皮尔逊III型分布频率计算后，得到 重现期P、降雨历时t、设计暴雨强度i_p的72组样本点(P_j，t_j；i_p，j)，j＝1～72，见表1。

表1实例1重现期历时设计雨强P-t-i_p关系值

利用本发明确定暴雨强度公式，其具体操作如下。

1.根据表1数据，采用文献[2，4]方法，推求4参数暴雨强度公式：

$i=\frac{9.8857+9.5129\mathrm{lgP}}{{(t+11.4082)}^{0.7342}},1\le P\le 100---\left(8\right)$

并计算由式(8)拟合表1数据的误差平方和Q＝0.2078。

进一步计算由式(8)拟合表1中各重现期相应的不同历时雨强的均方差σ：

$\sigma =\sqrt{\frac{\underset{j=1}{\overset{m_1}{\Sigma }}{(i_{p,i}-i_{g,j})}^2}{m_1}}---\left(9\right)$

式中i_p，j为根据计算重现期由理论频率曲线得到的雨强，mm/min，本例即为表1中某一重现期 相应的不同历时的雨强；i_g，j为根据重现期、降雨历时由暴雨强度公式计算的雨强，mm/min； m₁为降雨历时个数，本例m₁＝8。

计算结果见表2。

表2实例1四参数公式拟合雨强的均方差

可见，4参数暴雨强度公式对于重现期较大或较小的情况拟合效果不理想，拟合均方差 超过0.05mm/min，需判断暴雨公式中是否需要考虑修正系数d。

2.A与lgT关系趋势的判别方法

对于4参数暴雨强度公式中的参数n，b一定时，本发明计算某一重现期P相应的雨力A 的方法是：

根据某一重现期P的m₁组实测关系点(t_j，i_p，j)，j＝1～m₁，基于式(1)，求满足

$Q\left(A\right)=\underset{j=1}{\overset{m_1}{\Sigma }}{[i_{p,j}-A{(t_j+b)}^{-n}]}^2---\left(10\right)$

为最小的A值。即对A求偏导数：

$\frac{\partial Q}{\partial A}=-2\underset{j=1}{\overset{m_1}{\Sigma }}(i_{p,j}-A{(t_j+b)}^{-n}){(t_j+b)}^{-n}=0---\left(11\right)$

整理得：

$A=\frac{\underset{j=1}{\overset{m_1}{\Sigma }}{i}_{p,j}·{(t_j+b)}^{-n}}{\underset{j=1}{\overset{m_1}{\Sigma }}{(t_j+b)}^{-2n}}---\left(12\right)$

因此，利用式(12)即可计算各重现期P相应的雨力A，进而点绘(A，lgP)关系点，若不为 直线趋势，则要考虑修正系数d，使A＝A₁+C₁lg(P+d)，进而选配5参数暴雨强度公式(7)，否 则，可通过其他方法，提高拟合精度。

对实例1，根据表1数据、式(8)中参数n＝0.7342，b＝11.4082，利用式(12)计算各重现期 P相应的雨力A，并点绘(A，lgP)关系点，如图1，可见关系不为直线趋势，则要考虑修正系 数d。

3.含有重现期修正系数d的暴雨强度公式参数的确定方法

(1)对式(7)中参数A₁，C₁，d，b，n分别求偏导数

$\partial i/{\partial A}_1=1/{(t+b)}^n---\left(13\right)$

$\partial i/{\partial C}_1=\lg (P+d)/{(t+b)}^n---\left(14\right)$

$\partial i/\partial d=C_1/\left[(P+d)\ln 10{(t+b)}^n\right]---\left(15\right)$

$\partial i/\partial b=-n[A_1+C_1\lg (P+d)]/{(t+b)}^{n+1}---\left(16\right)$

$\partial i/\partial n=-\ln (t+b)[A_1+C_1\lg (P+d)]/{(t+b)}^n---\left(17\right)$

(2)确定暴雨强度公式(7)中参数的初值θ₍₀₎＝(A₁₍₀₎，C₁₍₀₎，d₍₀₎，b₍₀₎，n₍₀₎)′

将求得的4参数暴雨强度公式中的参数值作为初值A₁₍₀₎，C₁₍₀₎，b₍₀₎，n₍₀₎；根据A与lgT 关系点趋势，在d＞-P_min(P_min为表1样本点中的最小重现期)范围初拟修正系数d₍₀₎，使A与lg(P +d₍₀₎)近似为直线趋势。本发明实例1中初拟d₍₀₎＝-0.6，A＝A₁+C₁(lgP-0.6)近似为直线关系，如 图2所示。故本发明实例1参数的初值θ₍₀₎＝(9.8857，9.5129，-0.6，11.4082，0.7342)′

(3)由参数的初值θ₍₀₎、式(13)～式(17)以及各个重现期的M组观察值(P_j，t_j；i_p，j)，j＝1～M 及式(5)计算偏导数矩阵J(θ₍₀₎)及矩阵f(θ₍₀₎)。

(4)根据式(4)或式(6)计算θ₍₁₎。

(5)再以θ₍₁₎作为参数的初始值θ₍₀₎，重复步骤(3)，步骤(4)，根据给定的δ值(例如， δ＝0.0001)，经若干次递推迭代，可求得公式(7)中参数的估计值θ＝(A₁，C₁，d，b，n)′。

本发明实例1的计算结果：θ＝(A₁，C₁，d，b，n)′＝(12.3965，7.9181，-0.7194，11.4300， 0.7347)′，则

$i=\frac{12.3965+7.9181\lg (P-0.7194)}{{(t+11.4300)}^{0.7347}},1\le P\le 100---\left(18\right)$

由式(18)拟合表1雨强的误差平方和Q＝0.0309。进一步由式(18)计算拟合表1中各重现期不 同历时雨强的均方差σ，结果见表3。

表3实例1五参数公式拟合雨强的均方差

对实例1，比较4参数、含有修正系数d的5参数暴雨强度公式的拟合表1全部雨强数 据的误差平方和Q及拟合各重现期雨强的均方差(表2、表3)，可见，考虑修正系数d的5参 数暴雨强度公式显著提高了拟合精度。

(二)路径二——直接拟合法

根据暴雨资料选样得到的重现期P、实测暴雨强度i和降雨历时t的m组关系值(P_j，t_j； i_j)，j＝1～m，即为直接拟合法确定暴雨强度公式(7)中拟合参数的样本点。下面结合实例2说 明采用直接拟合法确定含有修正系数d的暴雨强度公式的方法。

实例2：A城市雨量站为国家基本气象站，利用该站1980～2011年的短历时自记雨量计 观测的雨量资料，采用年最大值选样法，分别得到降雨历时5min、10min、15min、20min、 30min、45min、60min、90min、120min、180min的10个历时的实测雨强样本，对其进行经 验频率及重现期计算后，得重现期历时实测雨强P-t-i关系表，如表4所示，共320组数据(限 于篇幅，仅列出部分数据)。

表4实例2重现期历时实测雨强P-t-i关系值

1.根据表4数据，采用文献[2，4]方法，推求4参数暴雨强度公式：

$i=\frac{18.4451+18.1057\mathrm{lgP}}{{(t+14.5130)}^{0.7921}},1.03\le P\le 33---\left(19\right)$

并计算由式(19)拟合表4中全部实测雨强的误差平方和为Q＝5.6316。进一步计算由式(19)拟 合表4中各重现期不同历时相应雨强的均方差σ，结果见表5。

表5实例2直接法四参数暴雨强度公式的拟合均方差(单位：mm/min)

由表5可见，4参数暴雨强度公式对于重现期较大或较小的情况拟合效果不理想，需判 断暴雨公式中是否需要考虑修正系数d。

2.A与lgT关系趋势的判别方法

根据式(19)中b＝14.5130，n＝0.7921及表4中各重现期33.00a，16.50a，…，1.03a的雨 强，利用式(12)计算各重现期的雨力A，并绘制A与lgP关系点如图3所示，关系点不是直线 趋势，而应考虑修正系数d，采用5参数暴雨强度公式(7)。

3.确定式(7)中参数的初值θ₍₀₎＝(A₁₍₀₎，C₁₍₀₎，d₍₀₎，b₍₀₎，n₍₀₎)′

将式(19)中的各参数值分别作为初值A₁₍₀₎，C₁₍₀₎，b₍₀₎，n₍₀₎；在d＞-P_min(P_min为表4样本点 中的最小重现期)范围初拟修正系数d₍₀₎＝-0.98，绘制A与lg(P-0.98)关系点如图4所示，关系 点近似为直线趋势。故实例2暴雨强度公式参数的初值θ₍₀₎＝(18.4451，18.1057，-0.98，14.5130， 0.7921)′。

4.利用式(4)或式(6)确定含有修正系数d的5参数暴雨强度公式，具体方法与路径一相应 环节的计算方法类似，不再赘述。对实例2，得

$i=\frac{26.2714+11.0536\lg (P-0.9728)}{{(t+14.7778)}^{0.7980}},1.03\le P\le 33---\left(20\right)$

由式(20)拟合表4中全部实测雨强的误差平方和为Q＝1.6593。进一步计算由式(20)拟合各重 现期不同历时实测雨强的均方差及全部重现期的均方差，结果见表6。

强调指出，表6为直接拟合各重现期实测雨强的均方差，不同于传统方法拟合由理论频 率曲线计算雨强的误差，故不能用规范^[6]中的误差限制值0.05mm/min来衡量。

表6实例2五参数暴雨强度公式的拟合均方差(单位：mm/min)

5.实例2不同公式拟合效果的比较。4参数公式(19)与5参数公式(20)拟合实测雨强样本 的误差平方和分别为5.6316、1.6593；比较表5、表6拟合各重现期实测雨强的均方差，在 32个重现期中，5参数公式有28个重现期的拟合均方差均显著小于4参数公式的结果，仅有 4个重现期的拟合均方差与4参数公式的结果相等；5参数公式拟合全部重现期实测雨强的均 方差也显著小于4参数的拟合结果。故5参数公式显著提高拟合精度，特别是能显著降低重 现期较小与较大时的拟合均方差。

有关文献

[1]北京市市政工程设计研究总院.给水排水设计手册，第5册，城镇排水，第二版[M].北京： 中国建筑工业出版社，2002.932-939.

[2]张子贤.用高斯-牛顿法确定暴雨公式参数[J].河海大学学报，1995，23(5)：106-111.

[3]方开泰.实用回归分析[M].北京：科学出版社，1988.168-172.

[4]张子贤，孙光东，孙建印，等.城市暴雨强度公式拟合方法的研究[J].水利学 报，2013，44(11)：1263-1271.

[5]袁志发，周静芋.多元统计分析[M].北京：科学出版社，2002.99-102.

[6]GB50014-2006.室外排水设计规范(2014版)[s].

[7]樊建军，王峰，陈鹏飞.利用MATLAB推导城市暴雨强度公式[J].中国给水排 水，2010，26(11)：114-115.

[8]柏绍光，黄英，方绍东，等.城市暴雨总公式与分公式精度分析[J].人民长 江，2006，37(8)：34-35.