2.2　暴雨强度公式

暴雨强度公式能够反映一定频率的暴雨在规定时段最不利时程分配的平均强度，它对研究当地降水气候特征（时空分布）以及预防大面积洪涝的发生发展都有一定的指导意义，同时它还是城市雨水排水系统规划与设计的基本依据之一。

2.2.1　暴雨强度公式的编制方法

目前我国各地已积累了完整的自动雨量记录资料，可采用数理统计法计算确定暴雨强度公式。水文统计学的取样方法有年最大值法和非年最大值法两类，国际上的发展趋势是采用年最大值法。我国《室外排水设计规范》（GB 50014）规定，具有20年以上自动雨量记录的地区宜采用年最大值法编制暴雨强度公式，小于20年自动雨量记录的地区可采用年多个样法。统计资料的年限越长，编制出的暴雨强度公式就越能反映当地的暴雨强度规律。下面分别对两种方法进行介绍。

（1）年多个样法

采用年多个样法取样推求暴雨强度公式的工作程序按先后时间顺序可分为3部分。

①选取雨样　a.资料年数。当地的自记雨量记录最少要求连续10年以上。b.降雨历时和重现期。降水历时采用9个历时，分别为5min、10min、15min、20min、30min、45min、60min、90min、120min；重现期采用8个，分别为0.25a、0.33a、0.5a、1a、2a、3a、5a、10a。c.取样方法。采用年多个样法，每年每个历时选择6～8个最大值，然后不论年次，将所有样本按各个历时做降序排列，从中选取资料年数的3～4倍的最大值作为编制公式的基本统计资料。

②频率分析　选取的各历时降雨资料，应用频率曲线加以调整。精度要求不太高时，可采用经验频率曲线；当精度要求较高时，可采用皮尔逊Ⅲ型分布曲线或指数分布曲线等理论频率曲线。根据确定的频率曲线，得出重现期、降雨强度和降雨历时三者的关系，即P-i-t关系表（见表2-2及图2-7）。

③推求公式　根据P、i、t关系值求得b、m、A1、C各个参数，可用解析法、图解与计算结合法或图解法等方法进行。将求得的各参数带入公式（2-17），即得当地的暴雨强度公式。

　　（2-17）

式中　q——设计暴雨强度，L/（s·hm2）；

t——降雨历时，min；

P——设计重现期，a；

A1、b、C、n——参数，根据统计方法计算确定。

（2）年最大值法

本方法适用于具有20年以上自记雨量记录的地区，有条件的地区可用30年以上的雨量系列。

①选取雨样　a.资料年数。当地的自记雨量记录最少要求连续20年以上。b.降雨历时和重现期。降水历时采用11个历时，分别为5min、10min、15min、20min、30min、45min、60min、90min、120min 、150min、180min；重现期采用8个，分别为2a、3a、5a、10a、20a、30a、50a、100a。c.取样方法。采用年最大值取样方法作为编制公式的基本统计资料。

②频率分析　选取的各历时降雨资料，应用频率曲线加以调整。精度要求不太高时，可采用经验频率曲线；当精度要求较高时，可采用皮尔逊III型分布曲线或指数分布曲线等理论频率曲线。根据确定的频率曲线，得出重现期、降雨强度和降雨历时三者的关系，即P-i-t关系值。

③推求公式　根据P-i-t关系值求得b、m、A1、C各个参数，可用解析法、图解与计算结合法或图解法等方法进行。为提高暴雨强度公式的精度，一般采用高斯-牛顿法。同样将求得的各参数带入公式（2-17），即得当地的暴雨强度公式。

2.2.2　暴雨强度公式编制案例

【例2.1】 某市有30年自记雨量记录，每年选择了各历时的最大暴雨强度值1个，然后将历年各历时的暴雨强度不论年次而按大小排列成表2-3。采用年最大值法编制暴雨强度公式。

【解】暴雨强度公式估算。

（1）根据公式计算各强度组的经验频率，式中的m为各强度组的序号数，也就是等于或大于该强度组的暴雨强度出现的次数。N值为观测资料的年数，也是暴雨强度的序号总数，本例的序号总数N为30。再根据公式计算相对应的重现期。

（2）理论频率曲线

经过我国水文工作者的大量拟合和分析工作，认为水文现象中最常用的理论频率曲线是皮尔逊III型分布曲线。本例题选用皮尔逊III型分布曲线进行理论频率曲线的计算，步骤如下。

①由实测资料，统计并计算平均值（）和变差系数（CV）

式中　xi——历年各历时的暴雨强度值，本例为表格中的270个数；

N——样本系列的序号总数，本例为30。

经计算，

②根据经验关系确定偏态系数（CS）

设计暴雨量：CS=3.5CV

设计最大流量：CV<0.5时，CS=（3-4）CV

CV>0.5时，CS=（2-3）CV

年径流及年降水量：CS=2CV

本例题求设计暴雨量，因此CS=3.5CV=1.82

③根据CS查皮尔逊Ⅲ型分布曲线的离均系数ϕP值表（见表2-4），得出不同概率下的ϕP值。

④由KP=ϕPCV+1计算模比系数KP，同时求得xP，列于表2-5。

（3）暴雨强度公式中各参数的推算

本例题采用解析法。地方特性参数A1、C、b、n采用非线性最小二乘估计法求解。非线性最小二乘估计法求暴雨强度公式中未定参数的数学原理为，找到一组未定参数的具体取值，使得基于这组参数的暴雨强度估计值与同一条件下暴雨强度的实际观测值误差的平方和最小。暴雨强度公式的参数辨识问题可转化为数学上典型的优化问题，其目标函数如下：

式中　A1、C、b、n——待优化计算确定的暴雨强度地方特性参数；

——当重现期为Pk，降雨历时为tk时，基于某组地方特性参数（A1、C、b、n）的暴雨强度估计值；

ik——当重现期为Pk，降雨历时为tk时，暴雨强度实际观测值。

①构造优化问题的m函数文件　为了降低编程计算的难度，本例题利用MATLAB非线性优化函数“lsqnonlin”对暴雨强度公式中的待定地方特性参数（A1、C、b、n）进行优化计算。lsqnonlin函数是MATLAB优化工具箱（optimitool）中的一员，能够很好地解决非线性最小二乘或数据拟合的问题，本文使用该算法对暴雨强度公式中的地方特性参数进行提取。

为了利用lsqnonlin函数，首先要根据已知条件（见表2-3），构造优化问题目标函数的m函数文件（Rainstormfitting.m）。该函数的输入变量为待优化的地方特性参数向量，输出为270组（30组重现期，9组降雨历时）不同观测条件下的暴雨强度估计值与观测值ik的残差向量。MATLAB代码如下。

a.function E = Rainstormfitting （A）

b.构造利用最小二乘法提取暴雨公式中地方特性参数（A1、C、b、n）的目标函数；

c.地方特性参数矩阵A为包含4个元素的向量，分别代表A1、C、b和n；

d.根据表格数据构造重现期向量P[重现期P的计算见步骤（1），将各序号各历时对应的270个重现期依次列出]；

P=[31，31，31，31，31，31，31，31，31，15.5，15.5，15.5，15.5，15.5，15.5，15.5，15.5，15.5，10.33，10.33，10.33，10.33，10.33，10.33，10.33，10.33，10.33，7.75，7.75，7.75，7.75，7.75，7.75，7.75，7.75，7.75，6.20，6.20，6.20，6.20，6.20，6.20，6.20，6.20，6.20，5.17，5.17，5.17，5.17，5.17，5.17，5.17，5.17，5.17，4.43，4.43，4.43，4.43，4.43，4.43，4.43，4.43，4.43，3.88，3.88，3.88，3.88，3.88，3.88，3.88，3.88，3.88，3.44，3.44，3.44，3.44，3.44，3.44，3.44，3.44，3.44，3.10，3.10，3.10，3.10，3.10，3.10，3.10，3.10，3.10，2.82，2.82，2.82，2.82，2.82，2.82，2.82，2.82，2.82，2.58，2.58，2.58，2.58，2.58，2.58，2.58，2.58，2.58，2.38，2.38，2.38，2.38，2.38，2.38，2.38，2.38，2.38，2.21，2.21，2.21，2.21，2.21，2.21，2.21，2.21，2.21，2.07，2.07，2.07，2.07，2.07，2.07，2.07，2.07，2.07，1.94，1.94，1.94，1.94，1.94，1.94，1.94，1.94，1.94，1.82，1.82，1.82，1.82，1.82，1.82，1.82，1.82，1.82，1.72，1.72，1.72，1.72，1.72，1.72，1.72，1.72，1.72，1.63，1.63，1.63，1.63，1.63，1.63，1.63，1.63，1.63，1.55，1.55，1.55，1.55，1.55，1.55，1.55，1.55，1.55，1.48，1.48，1.48，1.48，1.48，1.48，1.48，1.48，1.48，1.41，1.41，1.41，1.41，1.41，1.41，1.41，1.41，1.41，1.35，1.35，1.35，1.35，1.35，1.35，1.35，1.35，1.35，1.29，1.29，1.29，1.29，1.29，1.29，1.29，1.29，1.29，1.24，1.24，1.24，1.24，1.24，1.24，1.24，1.24，1.24，1.19，1.19，1.19，1.19，1.19，1.19，1.19，1.19，1.19，1.15，1.15，1.15，1.15，1.15，1.15，1.15，1.15，1.15，1.11，1.11，1.11，1.11，1.11，1.11，1.11，1.11，1.11，1.07，1.07，1.07，1.07，1.07，1.07，1.07，1.07，1.07，1.03，1.03，1.03，1.03，1.03，1.03，1.03，1.03，1.03]；

e.构造降雨历时向量t（将各序号对应的各历时依次列出）；

t=[5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120，5，10，15，20，30，45，60，90，120]；

f.构造暴雨强度观测值向量i（将各序号各历时对应的270个暴雨强度值依次列出）；

i=[3.02，2.89，2.227，2.07，1.93，1.664，1.32，0.999，0.803，2.78，2.40，2.14，1.925，1.773，1.46，1.225，0.884，0.751，2.76，2.18，2.04，1.855，1.523，1.329，1.123，0.831，0.663，2.60，2.14，1.973，1.83，1.513，1.287，1.098，0.801，0.628，2.58，2.12，1.86，1.80，1.507，1.08，1.087，0.794，0.601，2.58，2.10，1.86，1.695，1.363，1.078，0.978，0.764，0.573，2.52，2.05，1.84，1.655，1.337，1.053，0.950，0.743，0.565，2.50，2.03，1.833，1.625，1.317，1.053，0.947，0.673，0.563，2.48，1.96，1.80，1.575，1.303，1.016，0.883，0.641，0.484，2.46，1.93，1.773，1.50，1.30，1.016，0.813，0.553，0.477，2.44，1.85，1.72，1.50，1.26，0.980，0.767，0.552，0.467，2.36，1.84，1.713，1.49，1.237，0.944，0.760，0.546，0.452，2.34，1.84，1.70，1.485，1.16，0.920，0.725，0.539，0.435，2.30，1.83，1.653，1.465，1.143，0.911，0.725，0.532，0.431，2.26，1.80，1.633，1.42，1.14，0.896，0.702，0.514，0.430，2.26，1.79，1.58，1.36，1.067，0.887，0.697，0.500，0.414，2.22，1.76，1.507，1.35，1.067，0.853，0.667，0.473，0.386，2.22，1.75，1.487，1.30，1.053，0.778，0.652，0.469，0.378，2.20，1.72，1.467，1.28，1.04，0.773，0.645，0.468，0.356，2.10，1.70，1.46，1.27，1，0.773，0.633，0.450，0.354，2.08，1.67，1.46，1.255，0.993，0.758，0.585，0.444，0.351，2.08，1.67，1.413，1.235，0.983，0.758，0.580，0.434，0.349，2.06，1.63，1.413，1.215，0.983，0.738，0.572，0.398，0.334，2.06，1.60，1.407，1.20，0.967，0.733，0.570，0.394，0.329，2.04，1.59，1.367，1.17，0.960，0.727，0.567，0.392，0.328，2.04，1.57，1.333，1.155，0.957，0.722，0.560，0.390，0.327，2.02，1.56，1.333，1.155，0.920，0.673，0.555，0.390，0.314，2.02，1.55，1.32，1.15，0.917，0.660，0.515，0.390，0.308，2，1.55，1.273，1.14，0.893，0.656，0.510，0.387，0.305，2，1.52，1.193，1.105，0.883，0.647，0.505，0.387，0.302]；

g.I=zeros（1，270）； 构造暴雨强度估计值向量I的结构；

h.E=zeros（1，270）； 构造残差向量E的结构；

for j=1：270

I（j）=（A（1）*（1+A（2）*log10（T（j））））/（（t（j）+A（3））^A（4））； 计算基于地方参数向量A的暴雨强度估计值向量I；

　　end

E=i-I；% 得到基于地方参数向量A的残差向量

end

②利用MATLAB优化工具——lsqnonlin函数进行地方特性参数的提取

假设地方特性参数（A1、C、b、n）的初始猜测值分别为2 1 1 1；

A=[2 1 1 1]；

lsqnonlin算法参数的设定：使用Levenberg-Marquardt算法并在计算完成后显示迭代过程；

options=optimset（'Algorithm'，'Levenberg-Marquardt'，'Display'，'iter'）；

调用lsqnonlin函数进行计算

输出地方参数（A1、C、b、n）的优化结果（最小二乘估计）

由于算法要求，在起始搜索向量A后面须输入两个空矩阵，表示不对搜索的上下界进行设定；

a=lsqnonlin（@Rainstormfitting，A，[]，[]，options）

MATLAB程序运行结果如下。

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to

its initial value is less than the default value of the function tolerance.

<stopping criteria details>

a =

　 21.7823　　0.5595　 15.1645　　0.8215

可以看到，算法经过35次迭代后达到收敛，残差为1.94894，地方特性参数（A1、C、b、n）的优化计算结果为：21.78、0.56、15.16、0.82，因此最终的暴雨强度拟合公式为：

MATLAB运行图如图2-8所示。