基于规划求解的暴雨强度公式推求方法研究
2019-08-06王子奇邹浩东洪毅怡晖符荣松
李 莉， 冷 艺， 黄 俊， 张 赛， 王子奇， 邹浩东， 洪毅怡晖， 符荣松
(1.重庆大学 三峡库区生态环境教育部重点实验室, 重庆 400045； 2.中国交通建设股份有限公司， 北京 100088； 3.广东海洋大学 海洋与气象学院， 广东 湛江 524088； 4.重庆工业职业技术学院， 重庆 401120)
1 研究背景
近5年来，我国300多个城市均发生了不同程度的内涝灾害，而暴雨是形成洪涝灾害的主要原因[1]。 研究反映降雨规律的暴雨强度公式和参数对于城镇防洪、海绵城市建设[2]有着重要的意义。目前我国暴雨强度公式的推求基本按照《室外排水设计规范》(GB50014-2006)推荐的基本步骤和方法进行，另外还有不少文献也对暴雨强度公式编制中的一些技术问题进行了研究[2]，使得暴雨强度公式的编制逐步规范化。但目前对编制过程中基础数据异常值的剔除、Cv、Cs值对拟合曲线的影响、有约束的规划求解法在暴雨强度公式推求的应用等问题的研究较少。
现有的皮尔逊-III分布适线研究[6]和公式推求方法[8-11]很少涉及异常值的剔除，随着数据的更新、资料年限的增加，降雨资料中出现异常值的概率也随之增加，这直接影响暴雨强度公式的准确度。Cv、Cs值的大小不仅反映实测数据的分布规律，还影响暴雨强度-降雨历时-重现期(q-t-P)的关系，只有在明确把握Cv、Cs的意义和影响的基础上，才能根据样本分布情况正确估计Cv、Cs的大小关系，得到准确可靠的(q-t-P)表。此外，皮尔逊-III型曲线适线过程中推荐采用的目估适线法具有主观性，工作量大，精度不够高，通过约束条件采用规划求解不仅能解决研究者发现参数出现负值的问题[12]，还能提高拟合速度。因此本文以重庆市酉阳县1993-2013年共21 a的年最大降雨资料为例，研究了异常值的剔除方法、Cv和Cs值对拟合曲线的影响规律以及目标规划求解的应用，实现对公式的快速精确的推求。
2 降雨数据资料
本文收集了重庆市酉阳县1993-2013年共计21 a的逐分钟降雨资料，按年最大值法取样后将不同年份5～120 min共9个降雨历时下的最大雨量统计如下(表1)。
表1 降雨数据资料
3 暴雨强度公式拟合与分析
3.1 异常值的剔除
我国制定的城市暴雨强度公式编制和设计暴雨雨型确定技术导则(以下简称《导则》)[13]中提到了数据筛选，但是并未提出数据筛选的方法。异常值可能是由于雪融水、飓风、人为等非正常因素造成的，将它们纳入计算范围会导致拟合曲线与普遍降雨规律偏离较远，因此需要根据气象历史资料剔除高异常值和低异常值。美国关于洪水控制的技术导则[14](以下简称Bulletin 17B)中提出，选取样本时会出现遇到特大值的情况，需要进行专门分析确认。一般特大值是指比相应历史资料序列的平均值(计入特大值后的均值)大2倍以上的稀遇暴雨值，合理取舍特大值有利于提升理论频率曲线拟合的准确度。通过分析研究，当数据服从对数皮尔逊分布且偏态值在-3～3之间时，采用如下的方法剔除异常值较为合理：
(1)
(2)
小于qL的所有样本数据应被视为异常值而被剔除，而大于qH的数据只能在有气象资料支撑的基础上剔除最大值。当样本偏态系数大于0.4时，Bulletin 17B建议先剔除高异常值，偏态系数小于0.4时，先剔除低异常值，当偏态系数介于-0.4～0.4时，对高、低异常值的剔除顺序不作要求。
鉴于本文中酉阳的样本量仅有21 a，异常值出现的概率较低，同时修订公式时样本量不应小于20 a，对于每一个降雨历时在本次拟合中仅需检验是否存在一个异常值。同时，酉阳城市防涝仅需达到20年一遇的等级，故采用0.1显著水平下的KN值。在拟合酉阳暴雨强度公式时，按照上述方法首先分析数据是否存在异常值，然后根据计算分别剔除了原始资料中45、60、90、120 min的最大暴雨强度值。
3.2 Cv、Cs值对拟合曲线的影响规律
皮尔逊-III型分布的概率密度函数为：
(3)
由于累积分布函数表示的是小于等于某一暴雨强度q的概率，而理论频率是指大于等于某一暴雨强度q的概率,故按下列积分式计算理论频率：
Ptheoretical=P(q≥qtheoetical)=
(4)
对于公式(4)的求解方法可以参考罗雅文[17]的方法，令t=β(q-b),公式(4)变为：
Ptheoretical=P(t>ttheoretical)
(5)
通过以上分析，对Cv、Cs赋不同的数值计算皮尔逊-III型概率密度曲线及其积分后的累积概率曲线，曲线作图详见图1所示。
由图1(a)可知，当Cs不变时，Cv越大，概率密度曲线矮且宽，累积分布曲线越来越缓，由此推导Cv越大，皮尔逊-III型曲线可拟合的范围也越大，此时曲线能拟合出离散程度高的数据。当实际降雨数据计算的概率集中在均值附近时，在进行适线时可适当减小Cv，而当数据离散程度较大时，可以适当增大Cv。由图1(b)可知，当Cv一定时，Cs越大，概率密度曲线越高且窄，左偏程度越高,累积分布曲线中部越向左偏，上段越陡，下段越平缓,由此可推导，当数据正偏态程度越大，即小于平均值的数据占总数据的比例越大时，应该用更大的Cs/Cv曲线去拟合数据。
综上所述，在皮尔逊-III型曲线适线过程中，可通过计算样本均值，分析数据的离散与偏态程度，结合图1中数据离散程度与Cv、Cs变化的相关关系，来调整Cv/Cs比值，提高适线拟合速度。
3.3 基于规划求解的分布曲线拟合及适线
为得到指定重现期及指定降雨历时下的暴雨强度，需按照理论频率分布拟合曲线，将分散的原始数据拟合得出内在趋势，为暴雨强度公式的外延计算创造条件。理论频率分布曲线包括：皮尔逊-III型分布曲线(P-Ⅲ型曲线)、耿贝尔分布曲线(Gumbel曲线)和指数分布曲线。由于指数分布曲线适用于非年最大值法取样的数据[18]，耿贝尔曲线是皮尔逊-III型曲线Cs=1.140的一个特例，其弹性不足，对样本适应性较差[19]，因此本研究采用皮尔逊-III型曲线来拟合暴雨强度-降雨历时-重现期(q-t-P)的关系。
将通过规划求解优化拟合的皮尔逊-III型分布曲线与经验频率点共同体现在海森频率格纸[20]上，得到如图2所示的暴雨强度-降雨历时-重现期(q-t-P)关系图。
根据以上方法按照皮尔逊-III型分布拟合曲线，拟合的绝对均方误差为0.0557 mm/min，相对均方误差为6.27%。同时对未剔除异常值时误差进行了对比分析，发现未剔除异常值时理论频率曲线的绝对均方误差为0.0704 mm/min，相对均方误差为7.03%；剔除异常值后绝对均方误差值减小0.0147 mm/min，相对均方误差值减小0.76%。
3.4 暴雨强度公式推求
通过计算qtheoretical或者由图2查询得到一组包含5个重现期(酉阳县管网按20年一遇防涝等级设计，故只需20、10、5、3、2年的重现期)、9个降雨历时(5、10、15、20、30、45、60、90、120 min)的数据，分别采用最小二乘法和高斯牛顿法推导暴雨强度公式。
依据《室外排水设计规范》(GB50014-2006)[21]，暴雨强度公式定义为：
(6)
式中：q为暴雨强度，L/(s·hm2)；P为重现期，a；t为降雨历时，min；A1、b、C、n为与地方暴雨特性有关且需求解的参数，其中A1为雨力参数，即重现期为1 a时的1 min设计雨量，mm ；C为雨力变动参数；b为降雨历时修正参数，即对暴雨强度公式两侧求对数后能使曲线转为直线所加的时间参数，min；n为暴雨衰减指数，与重现期有关。
3.4.1 最小二乘法推求结果 为求得暴雨强度公式，对于公式(6)两侧取对数得：
lnq=ln 167A1+ln(1+ClgP)-nln(t+b)
(7)
令y=lnq，b0=ln 167A1，x1=ln(1+ClgP)，b2=-n,x2=ln (t+b)，公式(7)变为：
y=β0+x1+β2x2
(8)
(9)
(10)
应用最小二乘法求解暴雨强度公式曲线与皮尔逊-III频率分布设计值的比较如图2所示。最小二乘法求解暴雨强度公式的绝对均方误差为0.0654 mm/min，相对均方误差为3.96%。最小二乘法是以残差平方和最小为目标，对估量值求偏导从而求解参数，这种方法对异常值特别敏感，在计算中收敛速度不是很高，此外对于非线性函数需要进行线性转化，故最小二乘法在暴雨强度公式的拟合上具有一定局限性。
图1 不同Cv、Cs对概率密度函数和累积分布函数的影响
图2 暴雨强度-降雨历时-重现期(q-t-P)关系的经验频率点和理论频率曲线
(11)
高斯牛顿法采用Matlab程序求解，得暴雨强度公式为：
(12)
应用高斯牛顿法求解暴雨强度公式曲线与皮尔逊-III频率分布设计值的比较如图4所示。高斯牛顿法求解暴雨强度公式的绝对均方误差为0.0502 mm/min，相对均方误差为5.38%。
图3最小二乘法求解暴雨强度公式推求值 图4高斯牛顿法求解暴雨强度公式推求值
与P-III拟合值的比较 与P-III拟合值的比较
4 结 论
针对暴雨强度公式推求中存在的问题，以实际暴雨强度数据结合理论知识对推求过程的几个问题进行了研究，得到以下几个结论：
(1)通过研究分析提出了异常值的剔除方法，结合分析Cv、Cs的意义总结其对适线的影响规律。当数据服从对数皮尔逊分布且偏态值在-3～3之间时，采用对数形式的高低限值方法剔除异常值，合理取舍异常值有利于提升理论频率曲线拟合的准确度；当数据集中分布在均值左右时适当减小Cv，当数据正偏程度较大时，适当增大Cs/Cv，为方便准确地拟合出暴雨强度-降雨历时-重现期(q-t-P)曲线提供了支撑。
(2)通过约束条件，设置目标函数采用Excel规划求解Cv、Cs，与目估适线法相比，明显提高了拟合速度和精度。
(3)基于规划求解的最小二乘法推导的暴雨强度公式绝对均方误差为0.0654 mm/min，相对均方误差为3.96%；高斯牛顿法推导的暴雨强度公式绝对均方误差为0.0502 mm/min，相对均方误差为5.38%，通过高斯牛顿法推求的酉阳暴雨强度公式较合理。高斯牛顿法可采用最小二乘法的计算值作为初始值进行迭代计算，相比将而言高斯牛顿法具有更快的收敛速度和更高的拟合精度。
