考慮經(jīng)驗因素的暴雨頻率曲線最優(yōu)化擬合算法

姬鵬杰，杜坤，馮燕，周明，杜雨

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

近年來，中國城鎮(zhèn)內(nèi)澇災害頻發(fā)，極大危害了人們的生命財產(chǎn)安全[1]。在開展雨水管網(wǎng)設(shè)計、調(diào)蓄工程規(guī)劃時，許多城鎮(zhèn)目前仍采用1987版《室外排水設(shè)計規(guī)范》規(guī)定的暴雨強度公式。隨著城市化進程加速及全球氣候持續(xù)變暖，各地區(qū)降雨特性發(fā)生了較大變化，1987版規(guī)范的公式存在推源數(shù)據(jù)過舊、無法合理反映降雨特征的問題，因此，有必要對暴雨強度公式進行修編。暴雨頻率曲線擬合是推求暴雨強度公式必不可少的步驟，其通過假定降雨強度與頻率服從某一理論函數(shù)分布，采用實測數(shù)據(jù)對函數(shù)中的參數(shù)進行擬合，進而實現(xiàn)降雨強度與降雨頻率關(guān)系外延計算，并消除測量誤差影響、防止參數(shù)過擬合。

常用的暴雨頻率曲線包括皮爾遜III型分布曲線(P-III型曲線)、耿貝爾和指數(shù)分布曲線。耿貝爾和指數(shù)分布曲線是P-III型曲線在Cs=1.14和Cs=2時的特例，計算相對簡單，但擬合精度不高；P-III型頻率曲線擬合精度高，但計算復雜。針對P-III型曲線離均系數(shù)計算速度慢、精度低的問題，Che[2]應用Excel軟件簡化計算，王正發(fā)[3]指出Excel計算時存在數(shù)字發(fā)散區(qū)，提出應用Matlab軟件計算離均系數(shù)。針對暴雨頻率曲線擬合，Balin[4]通過標準化降水指數(shù)模型提高計算效率，崔俊蕊等[5]利用水文頻率分析軟件簡化試線過程，Mandal等[6]引入馬爾可夫鏈模型進行擬合，Wu等[7]引入置信區(qū)間法提高擬合精度。上述研究在提高計算效率及擬合精度兩個方面取得了一定成果，但高琳等[8]的最新研究指出，目前,大多數(shù)研究都忽視了實踐工程經(jīng)驗因素，將暴雨頻率曲線擬合作為單純的參數(shù)求解問題。高琳等[8]還提出，對于暴雨強度頻率曲線擬合，并非誤差越小越好，而應更多地照顧工程實際要求重現(xiàn)期段下的樣本，這樣得到的暴雨強度公式更加符合實際工程需求。值得注意的是，高琳等雖然給出了參數(shù)擬合準則，但未提出相應算法，其通過反復適線擬合參數(shù)，導致了巨大計算工作量。在考慮經(jīng)驗因素的情況下，如何高效實現(xiàn)暴雨頻率曲線參數(shù)擬合是值得研究的問題。

1 基于最優(yōu)化的暴雨頻率曲線參數(shù)擬合框架

選擇P-III型曲線作為理論暴雨頻率曲線，其密度函數(shù)為

(1)

(2)

(3)

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

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

(4)

式(4)的一階線性展開式為

(5)

由于理論暴雨強度計算涉及伽瑪函數(shù)和不完全伽瑪函數(shù)運算，使得雅克比矩陣解析式推導非常復雜。采用有限差分法計算雅克比矩陣，如式(6)所示。

(6)

(7)

根據(jù)式(7)可得

(8)

為改進迭代收斂性，在海塞矩陣的對角添加阻尼系數(shù)矩陣Γ[12]，可得

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

(9)

在參數(shù)擬合過程中，可設(shè)置阻尼系數(shù)Γ的取值隨迭代次數(shù)的增加而減小，擬合框架如圖1所示。

圖1 暴雨頻率曲線參數(shù)擬合框架圖Fig. 1 Frame diagram of parameters fitting for frequency curve of rainstor

2 算例分析

算例分析旨在利用實際降雨數(shù)據(jù)闡明：1)調(diào)節(jié)權(quán)重系數(shù)提高工程常用段擬合精度；2)調(diào)節(jié)權(quán)重系數(shù)避免暴雨頻率曲線相交；3)添加阻尼系數(shù)改進迭代收斂性。值得說明的是，該工程是短歷時排水系統(tǒng)，但所提出方法同樣適用于長歷時排澇系統(tǒng)的暴雨頻率曲線擬合。

2.1 調(diào)節(jié)權(quán)重系數(shù)提高工程常用段擬合精度

通過收集云南省保山市隆陽區(qū)1981—2013年33 a間實測降雨數(shù)據(jù)，整理出5、10、15、20、30、45、60、90、120、150、180 min共11個降雨歷時下的暴雨強度。對擬合結(jié)果進行分析發(fā)現(xiàn)，不同暴雨強度的擬合精度不同，降雨歷時越小，擬合精度越低。例如，5、10 min的擬合精度遠低于150、180 min擬合精度，其原因是降雨歷時越小，暴雨強度離均系數(shù)越大，尤其對最大值與最小值，往往偏離擬合曲線較遠，如圖2所示。以圖2中5 min降雨歷時下暴雨強度為例，闡明通過調(diào)節(jié)權(quán)重系數(shù)，提高工程常用重現(xiàn)期段擬合精度。

圖2 調(diào)節(jié)權(quán)重系數(shù)提高工程常用段擬合精度示意圖Fig. 2 Adjusting weight coefficient to improve the fitting precision of common section of Engineerin

表1 調(diào)整權(quán)重系數(shù)時殘差變化情況Table 1 Change of residual difference when adjusting weight coefficient

如表1所示，隨著其他段權(quán)重系數(shù)減小，整體殘差增大，工程常用段殘差減??；對應于圖2，曲線逐步向下偏移，使得適線結(jié)果與工程常用段樣本點更貼近。由此可見，通過改變權(quán)重系數(shù)能對適線結(jié)果進行微調(diào)，有效提高工程常用重現(xiàn)期段擬合精度。但值得說明的是，提高工程常用段擬合精度時，整體精度會不可避免的下降，且不同案例的精度變化不同。如果要定量給出精度取值或取值范圍，需要收集多個城市降雨數(shù)據(jù)進行綜合分析，工作量巨大。鑒于筆者的目的在于證明所提出的算法能高效調(diào)整二者精度，故不對上述問題進行深入分析。

2.2 調(diào)節(jié)權(quán)重系數(shù)避免理論頻率曲線相交

圖3 不同歷時下理論頻率曲線相交Fig. 3 Intersecting theoretical frequency curves at different diachronic times

圖4 調(diào)節(jié)權(quán)重系數(shù)使頻率曲線不相交Fig. 4 Adjusting the weight coefficient to disjoint the frequency curv

2.3 添加阻尼系數(shù)改進迭代收斂

圖5 普通高斯牛頓迭代法收斂情況Fig. 5 Convergence of ordinary Gauss Newton

在海塞矩陣對角添加阻尼系數(shù)，如式(9)所示，阻尼系數(shù)初值取1，分別采用一次方衰減(1/k)及二次方衰減(1/k2)進行算法測試，其中,k為迭代次數(shù)[15]。

如圖6所示，當阻尼系數(shù)一次方衰減時，迭代49次達到收斂精度要求(ε<10-4)；當阻尼系數(shù)二次方衰減時，迭代33次達到收斂精度要求。雖然二次方衰減法能更快的達到收斂精度，但當初值偏離真值較遠時，不能保證迭代收斂[16]。建議先采用二次方衰減法進行初算，若迭代不收斂，則采用一次方或更低衰減方式，可加大阻尼系數(shù)初值進一步改進迭代收斂性。

圖6 添加阻尼系數(shù)改進迭代收斂Fig.6 Adding damping coefficient to improve iterative convergenc

2.4 基于差分進化算法的優(yōu)化精度檢驗

表2二種算法優(yōu)化結(jié)果及目標函數(shù)殘差
Table2Theoptimizationresultsoftwoalgorithmsandthedifferenceoftheresidualnumberoftheobjectivefunction

算法均值mCvCs目標函數(shù)殘差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 結(jié)論

研究了考慮經(jīng)驗因素時的暴雨頻率曲線擬合算法，以云南省保山市隆陽區(qū)33 a實測降雨資料為例論證了算法的可行性，得到如下結(jié)論：

3)通過調(diào)節(jié)權(quán)重系數(shù)能方便對適線結(jié)果進行微調(diào)，避免不同歷時暴雨頻率曲線相交的問題，并提高工程常用重現(xiàn)期段擬合精度，但如何平衡工程常用段與整體精度需要進一步研究。

5)采用差分進化算法搜索全局最優(yōu)解，驗證了所提出的加權(quán)阻尼高斯牛頓算法同樣能獲得全局最優(yōu)解。文中所有涉及運算都進行了編程，程序運算時間小于10 s，使繁復的適線工作能在5～10 min完成。

[1] 車伍, 楊正, 趙楊,等. 中國城市內(nèi)澇防治與大小排水系統(tǒng)分析[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] 王正發(fā). MATLAB在P-Ⅲ型分布離均系數(shù)φ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] 崔俊蕊, 王政然, 梁爽,等. 城市設(shè)計暴雨頻率曲線的擬合及參數(shù)優(yōu)化[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.