引言

全球变暖增加了城市极端暴雨的发生概率^[1]，随着我国城镇化速度的加快，也使得城市降雨频率与强度增大^[2-4]。近年来暴雨对于城市的威胁越来越严峻，不仅造成了惨重的人员伤亡和财产损失，也严重影响了城市经济的正常发展^[5]。2012年7月21日北京市大部分地区出现了大暴雨、局地特大暴雨天气过程，强降水造成北京市78人死亡，经济损失达116.4亿元^[6]。

2007年7月18日下午至19日凌晨，山东省济南市出现的强降雨造成市区街道严重积水，洪灾造成34人死亡，并造成重大财产损失^[7]。2007年7月29日山西南部强降水仅运城市受灾人口就达20万余人，死亡13人、失踪3人、受伤14人，直接经济损失达14.96亿元^[8]。2010年住房和城乡建设部组织开展的全国范围内351个城市调研发现，在2008—2010年的3 a间，全国有62%的城市均曾发生过内涝事件，内涝发生3次以上的城市有137个^[9]。因此，系统有效地解决城市排水防涝难题已经成为当前国家面临的严峻问题。确定城市暴雨强度公式是城镇给排水系统的规划与室外排水设计的重要基础性工作。为此许多学者对暴雨强度公式确定方法进行了研究，现行的国家规范的暴雨强度公式为一个非线性模型^[10]。过去常用图解法、最小二乘法、或两者相结合等传统方法确定^[11]，随着计算技术的发展，许多数理计算方法被应用于暴雨强度公式参数的优化中，呈现出了各种形式的推求方法。许拯民等^[12]采用麦夸尔特法、倍比搜索法对暴雨强度公式进行推求；陈正洪等^[13]先将公式线性化，确定出未知参数取值范围，优化了算法；还有SCE-UA算法^[14]、遗传算法^[15]、高斯牛顿法^[16]、切比雪夫多项式逼近计算伽玛函数法^[17]、GEV分布模型法^[18]、Lingo-BFGS算法^[19]、粒子群算法^[20]、微分进化算法^[21]、曲面搜索法^[22]等方法也被应用在暴雨强度公式参数优化中。但对区域暴雨强度公式参数的空间分布特征研究比较少见，本文对适合山西省降水特点的暴雨强度公式推算方法及空间分布特征进行研究，首先采用无约束的非线性优化求解方法，同时拟合4个参数，得到误差最小的暴雨强度公式参数，然后对暴雨强度公式参数进行分组拟合试验，得到另外一组暴雨强度公式参数，以最小误差的暴雨强度公式结果为标准，最终给出合理的暴雨强度公式参数空间分布。这对目前我国各地都在进行的新一轮暴雨强度公式修编工作具有重要的现实指导意义。

1 资料与方法 1.1 资料

资料为山西省107个站1981—2012年5 min、10 min、15 min、20 min、30 min、45 min、60 min、90 min、120 min、180 min十个历时的年最大降水量，来源于山西省气象信息中心。

1.2 暴雨强度公式形式及误差计算

根据我国现行室外排水设计规范，暴雨强度公式形式为^[10]：

$ i = \frac{{A(1 + C\lg p)}}{{{{(t + b)}^n}}} $

(1)

暴雨强度公式均方差计算公式为：

$ \sigma = \sqrt {\frac{{\sum {{{({i_g} - {i_j})}^2}} }}{m}} $

(2)

式(2)中，i_g为i-t-p关系表中雨强(单位：mm·min^-1)；i_j为所求的暴雨强度公式计算的雨强(单位：mm·min^-1)；m为降雨强度的项数(m=历时总数×重现期总数)。国家规范规定，在一般降雨强度的地方，均方差不宜大于0.05 mm·min^-1。

1.3 频率调整

在推算暴雨强度公式参数时，需要确定降雨强度i、降雨历时t和重现期p三者的关系，即i-t-p关系表。由于实测数据的时间序列长度有限，直接用经验频率点绘出的曲线往往不能满足推求稀遇频率的要求，目估或外延会产生较大的误差。气象要素极值作为随机变量在一般数学意义上是不稳定的，但在概率意义上却可以是稳定的^[23]，可以采用理论频率曲线对经验频率加以调整，作为定线和外延的依据。理论频率曲线具有一定的数学形式，故能探讨其统计性质，减少主观的目估或外延弊病，并能对有限样本得到的经验频率存在的偏差有一定程度的矫正。邓培德等^[24]研究表明，年最大降水量原始分布为偏态铃形。王振华等^[25]也对山西省年最大降水量原始分布进行研究，指出其分布也为偏态铃形。本文采用的选样法为年最大值法，其频率分布函数满足皮尔逊Ⅲ型和耿贝尔型分布，为了简便起见，利用耿贝尔函数分析雨强频率。

耿贝尔分布函数为^[26]：

$ F({x_P}) = P(x < {x_P}) = {e^{ - {e^{ - \alpha ({x_P} - \beta )}}}}\;\;\left( {\alpha > 0, - \infty < \beta < \infty } \right) $

(3)

公式(3)中，α为尺度参数，与均方差成反比；β是位置参数，反映频率分布集中在数轴上的位置。

需要推求设计频率f = F(x ≥ x_f)所对应极值x_f，通过推算，可得到：

$ {x_f} = \beta - \frac{1}{n}\ln [ - \ln (1 - f)] $

(4)

耿贝尔法是一种直接与经验概率相结合的参数估计方法。假定数据有序列x₁ ≤ x₂… ≤ x_i… ≤ x_n，设Y_i =-ln[-ln(1-f_i)]，参数Y_i =-ln[-ln(1-p)]，α、β可根据下式得到：

$ \sigma \left( x \right) = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{({x_i} - \bar x)}^2}} } $

(5)

$ \sigma \left( y \right) = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{({y_i} - \bar y)}^2}} } $

(6)

$ \alpha = \frac{{\sigma \left( y \right)}}{{\sigma \left( x \right)}} $

(7)

$ \beta = \bar x - \frac{{\bar y}}{\alpha } $

(8)

1.4 暴雨强度公式的推算

暴雨强度公式参数的求解过程是整个暴雨强度公式编制的核心。由暴雨强度公式的基本形式可知，参数求解就是暴雨强度公式拟合出i-t-p关系表中的暴雨强度误差最小，这其实就是一个非线性模型参数的优化过程。即：

有一组数据(x1_i，x2_i，x3_i，...，y_i)(i = 1，2，...，n)，要用已知函数y = F(params，x)拟合，x为一组变量，params =(pa₁，pa₂，...，pa_l)(l ≤ n)为待定参数，为实现参数优化，需构造目标函数$R = \sum\limits_{i = 1}^n {{{[{y_i} - F(params, {x_i})]}^2}} $，于是问题就成为求解minR条件下的params。为了达到此目的，可以利用多种数理算法。一些软件也开发了相应的功能，如Matlab中就有为求解优化问题专门开发的优化工具箱，excel中规划求解也能够得到最优参数。

本文借助于Matlab中无约束的非线性优化求解方法，其工作原理是采用Nelder-Mead下山单纯形(Downhill Simplex)法收缩局域最小值点^[27]。无约束的非线性优化求解方法对函数的形式和参数的个数没有限制，可广泛运用于各类模型的参数估计，实现多个参数的同时寻优，避免了图解类方法试凑和反复调整等繁琐工作，使得推算误差最小的参数成为可能并且其结果不会由于推算人员变化而有大的变化。在推算山西省各地暴雨强度公式过程中，其结果都是收敛的，没有参数溢出情况。

对于求解暴雨强度公式中的参数优化过程，定义目标函数为：

$ R = \sum\limits_{l = 1}^m {{{\left[ {{i_l} - \frac{{A(1 + c\lg {P_l})}}{{{{({t_l} + b)}^n}}}} \right]}^2}} $

(9)

式(9)中, i_l、t_l、P_l为i-t-p关系表中的第l点数据的相应雨强、历时、重现期。m为表中的数据个数，求解minR条件下的A、b、C、n即为最优解。本文在推求暴雨强度公式参数时，降雨历时采用5 min、10 min、15 min、20 min、30 min、45 min、60 min、90 min、120 min、180 min共十个历时，降雨重现期采用2 a、3 a、5 a、10 a、20 a、30 a、50 a、100 a共8个重现期。下面以一个个例来说明暴雨强度公式参数优化求解的过程和结果。

表 1为大同县站根据耿贝尔分布函数调整频率得到的i-t-p关系表，根据关系表，以minR最小为目标，得到的暴雨强度公式均方差为0.028 mm·min^-1。暴雨强度公式如下：

$ i = \frac{{10.511\;1(1 + 1.206\;6\lg p)}}{{{{(t + 9.505\;6)}^{0.864\;2}}}} $

(10)

2 山西不同重现期的雨强和暴雨强度公式参数空间分布 2.1 山西省自然概况

山西省省境轮廓大体呈平行四边形，地形复杂，山地、高原、丘陵占全省面积的72%，按地形起伏特点，可以将山西省分为东部山区、西部高原区和中部盆地区3大部分，总的地势是东西两侧为山地和丘陵的隆起，中部为一列串珠式断陷盆地。山西省降水的总体特征是从东南向西北递减，山地降水量普遍多于川谷。以24 h降水量大于等于50 mm为暴雨标准，山西省暴雨日数的分布特点也是东南多、西北少^[28]。

2.2 山西省重现期的空间分布

根据上述频率调整结果，采用同时拟合4个参数的求解方法得到山西省各站绝对误差最小的暴雨强度公式，进而根据暴雨强度公式可推算出不同历时各重现期的暴雨强度。图 1给出山西省20 min、60 min历时5 a、20 a、50 a重现期的暴雨雨强空间分布，从中可见，它们的分布趋势一致，雨强总体特征是从东南向西北递减，山区大于盆地，高值区在晋东南一带。这与山西省年降水量及暴雨日数分布基本一致。对于20 min历时，重现期5 a的雨强介于0.8~1.5 mm·min^-1之间，重现期20 a的雨强介于1.3~2.1 mm·min^-1之间，重现期50 a的雨强介于1.5~2.5 mm·min^-1之间；对于60 min历时，重现期5 a的雨强介于0.4~0.9 mm·min^-1之间，重现期20 a的雨强介于0.6~1.2 mm·min^-1之间，重现期50 a的雨强介于0.7~1.4 mm·min^-1之间。

2.3 暴雨强度公式参数分布

综上可知，由暴雨强度公式计算的不同重现期的雨强空间分布有一定规律，但将上述参数点绘制出来(图略)，却发现其空间分布没有规律性。究其原因，是因为同时拟合4个参数暴雨强度公式时，只追求绝对误差最小，有时误差减小0.001 mm·min^-1，A和b可能会增大10倍以上。为克服这一问题，使参数也能够和地形、雨强有一定的匹配性，根据以上方法推算出全省107站暴雨强度公式，计算4个参数的全省平均值，并将其作为初始参数，把4个参数分成两组，如第一组参数是1个(或2个)，第二组参数就是3个(或2个)，两组参数分多步进行拟合，其步骤如下：

步骤1：以全省各站最优参数的平均值为各站初始参数；

步骤2：在第一组参数不变的条件下，求解第二组参数的最优解；

步骤3：在上一步确定的参数基础上，在第二组参数不变条件下，求解第一组参数的最优解；

步骤4:在上一步参数基础上，重复步骤2，如此反复。

经过大量的实验最终发现，把4个参数分成两组，一组参数是b，另一组参数是A、C、n，反复求各组参数的最优解，反复10次以后误差变化很小，且已逼近同时拟合4个参数的暴雨强度公式误差，参数的空间分布也有一定的规律。图 2为山西省暴雨强度公式参数空间分布，从中可见，参数A和C基本上为经向分布，b和n为纬向分布；参数A的分布是东部大于西部，其值范围在4.2~15.4之间；参数C的分布基本上是东部小、西部大，最大值在太原盆地以西，其值变化范围比较小，在0.7~1.9之间；参数b的分布基本是南部大于北部，其值变化范围在7.5~8.7之间；参数n的分布基本上是南部小、北部大，其值范围在0.50~0.98之间。这些特征与山西省地形及气候分布有一定的相似性，这些相似性在推广应用暴雨强度公式时具有重要意义。

2.4 参数空间分布的合理性分析

分析分组拟合参数和同时拟合4个参数得到的暴雨强度公式的误差可知，全省平均误差仅增大了0.013 mm·min^-1。表 2为分组拟合参数的暴雨强度公式误差变化情况，从中可见，全省一半站点以上的误差增大不超过0.01 mm·min^-1。根据分组拟合参数得到暴雨强度公式进而求得的重现期暴雨雨强空间分布(图略)与图 1基本相同，说明分组拟合参数得到的暴雨强度公式与同时拟合4个参数得到的暴雨强度公式是一致的。

3 结论与讨论

根据山西省107个气象站32 a的10个历时年最大降水量资料，利用Gumbel分布函数调整频率，得出雨强(i)-历时(t)-重现期(p)关系表。在此基础上，首先采用无约束的非线性优化求解方法拟合得到误差最小的暴雨强度公式参数，然后对暴雨强度公式参数进行分组拟合试验，得到另外一组暴雨强度公式，以最小误差的暴雨强度公式结果为标准，最终给出合理的暴雨强度公式参数空间分布。

(1) 山西省20 min、60 min历时的重现期5 a、20 a、50 a的暴雨雨强空间分布一致，总体特征是从东南向西北递减，山区大于盆地，雨强高值区在晋东南一带，这与山西省降水分布基本一致。

(2) 山西省暴雨强度公式中的地区参数空间分布特征是：参数A和C为经向分布，b和n为纬向分布；参数A的分布是东部大于西部，其值范围在4.2~15.4之间；参数C的分布基本上是东部小、西部大，最大值在太原盆地以西，其值变化范围较小，在0.7~1.9之间；参数b的分布基本是南部大于北部，其值变化范围在7.5~8.7之间；参数n的分布基本上是南部小、北部大，其值范围在0.50~0.98之间。

(3) 无约束的非线性优化求解方法对函数的形式和参数的个数没有限制，可广泛运用于各类模型的参数估计，其结果可用于缺乏气象观测站的区域根据下垫面和气候分区特点构建暴雨强度公式。

山西省暴雨强度公式中参数的空间分布特征与山西省地形特征及气候分布有一定的相似性，这些相似性的内在联系有待进一步分析探讨。