考虑经验因素的暴雨频率曲线最优化拟合算法

姬鹏杰，杜坤，冯燕，周明，杜雨

(1.昆明理工大学 建筑工程学院，昆明 650500；2.中建二局第三建筑工程有限公司，武汉430022)

近年来，中国城镇内涝灾害频发，极大危害了人们的生命财产安全[1]。在开展雨水管网设计、调蓄工程规划时，许多城镇目前仍采用1987版《室外排水设计规范》规定的暴雨强度公式。随着城市化进程加速及全球气候持续变暖，各地区降雨特性发生了较大变化，1987版规范的公式存在推源数据过旧、无法合理反映降雨特征的问题，因此，有必要对暴雨强度公式进行修编。暴雨频率曲线拟合是推求暴雨强度公式必不可少的步骤，其通过假定降雨强度与频率服从某一理论函数分布，采用实测数据对函数中的参数进行拟合，进而实现降雨强度与降雨频率关系外延计算，并消除测量误差影响、防止参数过拟合。

常用的暴雨频率曲线包括皮尔逊III型分布曲线(P-III型曲线)、耿贝尔和指数分布曲线。耿贝尔和指数分布曲线是P-III型曲线在Cs=1.14和Cs=2时的特例，计算相对简单，但拟合精度不高；P-III型频率曲线拟合精度高，但计算复杂。针对P-III型曲线离均系数计算速度慢、精度低的问题，Che[2]应用Excel软件简化计算，王正发[3]指出Excel计算时存在数字发散区，提出应用Matlab软件计算离均系数。针对暴雨频率曲线拟合，Balin[4]通过标准化降水指数模型提高计算效率，崔俊蕊等[5]利用水文频率分析软件简化试线过程，Mandal等[6]引入马尔可夫链模型进行拟合，Wu等[7]引入置信区间法提高拟合精度。上述研究在提高计算效率及拟合精度两个方面取得了一定成果，但高琳等[8]的最新研究指出，目前,大多数研究都忽视了实践工程经验因素，将暴雨频率曲线拟合作为单纯的参数求解问题。高琳等[8]还提出，对于暴雨强度频率曲线拟合，并非误差越小越好，而应更多地照顾工程实际要求重现期段下的样本，这样得到的暴雨强度公式更加符合实际工程需求。值得注意的是，高琳等虽然给出了参数拟合准则，但未提出相应算法，其通过反复适线拟合参数，导致了巨大计算工作量。在考虑经验因素的情况下，如何高效实现暴雨频率曲线参数拟合是值得研究的问题。

1 基于最优化的暴雨频率曲线参数拟合框架

选择P-III型曲线作为理论暴雨频率曲线，其密度函数为

(1)

(2)

(3)

f(κk+Δκk,Τp)=[Xp-h(κk+Δκk,Τp)]T·

W[Xp-h(κk+Δκk,Τp)]

(4)

式(4)的一阶线性展开式为

(5)

由于理论暴雨强度计算涉及伽玛函数和不完全伽玛函数运算，使得雅克比矩阵解析式推导非常复杂。采用有限差分法计算雅克比矩阵，如式(6)所示。

(6)

(7)

根据式(7)可得

(8)

为改进迭代收敛性，在海塞矩阵的对角添加阻尼系数矩阵Γ[12]，可得

Δκk=[J(κk,Τp)TWJ(κk,Τp)+Γ]-1·

(9)

在参数拟合过程中，可设置阻尼系数Γ的取值随迭代次数的增加而减小，拟合框架如图1所示。

图1 暴雨频率曲线参数拟合框架图Fig. 1 Frame diagram of parameters fitting for frequency curve of rainstor

2 算例分析

算例分析旨在利用实际降雨数据阐明：1)调节权重系数提高工程常用段拟合精度；2)调节权重系数避免暴雨频率曲线相交；3)添加阻尼系数改进迭代收敛性。值得说明的是，该工程是短历时排水系统，但所提出方法同样适用于长历时排涝系统的暴雨频率曲线拟合。

2.1 调节权重系数提高工程常用段拟合精度

通过收集云南省保山市隆阳区1981—2013年33 a间实测降雨数据，整理出5、10、15、20、30、45、60、90、120、150、180 min共11个降雨历时下的暴雨强度。对拟合结果进行分析发现，不同暴雨强度的拟合精度不同，降雨历时越小，拟合精度越低。例如，5、10 min的拟合精度远低于150、180 min拟合精度，其原因是降雨历时越小，暴雨强度离均系数越大，尤其对最大值与最小值，往往偏离拟合曲线较远，如图2所示。以图2中5 min降雨历时下暴雨强度为例，阐明通过调节权重系数，提高工程常用重现期段拟合精度。

图2 调节权重系数提高工程常用段拟合精度示意图Fig. 2 Adjusting weight coefficient to improve the fitting precision of common section of Engineerin

表1 调整权重系数时残差变化情况Table 1 Change of residual difference when adjusting weight coefficient

如表1所示，随着其他段权重系数减小，整体残差增大，工程常用段残差减小；对应于图2，曲线逐步向下偏移，使得适线结果与工程常用段样本点更贴近。由此可见，通过改变权重系数能对适线结果进行微调，有效提高工程常用重现期段拟合精度。但值得说明的是，提高工程常用段拟合精度时，整体精度会不可避免的下降，且不同案例的精度变化不同。如果要定量给出精度取值或取值范围，需要收集多个城市降雨数据进行综合分析，工作量巨大。鉴于笔者的目的在于证明所提出的算法能高效调整二者精度，故不对上述问题进行深入分析。

2.2 调节权重系数避免理论频率曲线相交

图3 不同历时下理论频率曲线相交Fig. 3 Intersecting theoretical frequency curves at different diachronic times

图4 调节权重系数使频率曲线不相交Fig. 4 Adjusting the weight coefficient to disjoint the frequency curv

2.3 添加阻尼系数改进迭代收敛

图5 普通高斯牛顿迭代法收敛情况Fig. 5 Convergence of ordinary Gauss Newton

在海塞矩阵对角添加阻尼系数，如式(9)所示，阻尼系数初值取1，分别采用一次方衰减(1/k)及二次方衰减(1/k2)进行算法测试，其中,k为迭代次数[15]。

如图6所示，当阻尼系数一次方衰减时，迭代49次达到收敛精度要求(ε<10-4)；当阻尼系数二次方衰减时，迭代33次达到收敛精度要求。虽然二次方衰减法能更快的达到收敛精度，但当初值偏离真值较远时，不能保证迭代收敛[16]。建议先采用二次方衰减法进行初算，若迭代不收敛，则采用一次方或更低衰减方式，可加大阻尼系数初值进一步改进迭代收敛性。

图6 添加阻尼系数改进迭代收敛Fig.6 Adding damping coefficient to improve iterative convergenc

2.4 基于差分进化算法的优化精度检验

表2二种算法优化结果及目标函数残差
Table2Theoptimizationresultsoftwoalgorithmsandthedifferenceoftheresidualnumberoftheobjectivefunction

算法均值mCvCs目标函数残差AB1.6958820.2848012.8495360.0882621.6957580.2848342.8494380.0882621.6959360.2848162.8480230.0882621.6958510.2846972.8462250.0882621.6958550.2848002.8474490.0882621.6958340.2848192.850060.0882621.6958680.2847692.8477930.0882621.6958650.2847622.8482600.0882621.6958680.2847582.8473010.0882621.6959480.2847632.8486090.0882621.6958060.2848412.8480030.088262

3 结论

研究了考虑经验因素时的暴雨频率曲线拟合算法，以云南省保山市隆阳区33 a实测降雨资料为例论证了算法的可行性，得到如下结论：

3)通过调节权重系数能方便对适线结果进行微调，避免不同历时暴雨频率曲线相交的问题，并提高工程常用重现期段拟合精度，但如何平衡工程常用段与整体精度需要进一步研究。

5)采用差分进化算法搜索全局最优解，验证了所提出的加权阻尼高斯牛顿算法同样能获得全局最优解。文中所有涉及运算都进行了编程，程序运算时间小于10 s，使繁复的适线工作能在5～10 min完成。

[1] 车伍, 杨正, 赵杨,等. 中国城市内涝防治与大小排水系统分析[J]. 中国给水排水, 2013, 29(16): 13-19.

CHE W, YANG Z, ZHAO Y, et al. Analysis of urban flooding control and major and minor drainage systems in China [J]. China Water & Wastewater, 2013, 29(16):13-19. (in Chinese)

[2] CHE G W. Pearson-III frequency curve plotting in Excel table [J]. Applied Mechanics & Materials, 2014, 556-562: 5829-5834.

[3] 王正发. MATLAB在P-Ⅲ型分布离均系数φp值计算及频率适线中的应用[J]. 西北水电, 2007(4): 1-4.

WANG Z F. MATLAB programming language used to calculate variation coefficientφpof the P-III distribution and fit a frequency curve [J]. Northwest Hydropower, 2007(4): 1-4. (in Chinese)

[4] BLAIN G C. Standardized precipitation index based on pearson type III distribution [J]. Revista Brasileira De Meteorologia, 2014, 26(2):167-180.

[5] 崔俊蕊, 王政然, 梁爽,等. 城市设计暴雨频率曲线的拟合及参数优化[J]. 水电能源科学, 2014(11): 48-51.

CUI J R, WANG Z R, LIANG S, et al. Fitting and parameter optimization of urban design storm frequency curve [J]. Water Resources and Power, 2014(11): 48-51. (in Chinese)

[6] MANDAL K G, PADHI J, KUMAR A, et al. Analyses of rainfall using probability distribution and Markov chain models for crop planning in Daspalla region in Odisha, India [J]. Theoretical and Applied Climatology, 2015, 121(3):517-528.

[7] WU Y C, LIU J J, SU Y F, et al. Establishing acceptance regions for L-moments based goodness-of-fit tests for the Pearson type III distribution [J]. Stochastic Environmental Research and Risk Assessment, 2014, 26(6): 873-885.

[8] 高琳, 周玉文, 唐颖, 等. 城市暴雨强度公式皮尔逊Ⅲ型适线问题研究[J]. 给水排水, 2016(8): 47-51.

GAO L, ZHOU W Y, TANG Y, et al. Research on the fitting of pearson type III in urban storm water intensity equation [J]. Water & Wastewater Engineering, 2016(8): 47-51. (in Chinese)

[9] RAFIQ A, RAFIULLAH M. Some multi-step iterative methods for solving nonlinear equations [J]. Computers & Mathematics with Applications, 2009, 58(8): 1589-1597.

[10] FAIRBANK M, ALONSO E. Efficient calculation of the Gauss-Newton approximation of the Hessian matrix in neural networks [J]. Neural Computation, 2012, 24(3): 607-610.

[11] ZHANG J, GENG X, DAI R. Analysis on two approaches for high order accuracy finite difference computation [J]. Applied Mathematics Letters, 2012, 25(12): 2081-2085.

[12] HAN Q, ZHANG Q S. An upper bound for hessian matrices of positive solutions of heat equations [J]. The Journal of Geometric Analysis, 2016, 26(2): 715-749.

[13] ERINA M Y,IZMAILOV A F. The Gauss-Newton method for finding singular solutions of systems of nonlinear equations [J]. Anz Journal of Surgery, 2015, 58(2): 395-405.

[14] PHAN A H, TICHAVSKY P, CICHOCKI A. Low complexity damped Gauss-Newton algorithms for CANDECOMP /PARAFAC[J]. Siam Journal on Matrix Analysis& Applications, 2013, 34(1): 126-147.

[15] SHEHU Y. Iterative method for fixed point problem, variational inequality and generalized mixed equilibrium problems with applications [J]. Journal of Global Optimization, 2013, 52(1): 57-77.

[16] JIN Q. Further convergence results on the general iteratively regularized Gauss-Newton methods under the discrepancy principle [J]. Mathematics of Computation, 2013, 82(283): 1647-1665.

[17] MACIEL L, GOMIDE F, BALLINI R. A differential evolution algorithm for yield curve estimation [J]. Mathematics & Computers in Simulation, 2016, 129:10-30.