考慮經驗因素的暴雨頻率曲線最優化擬合算法

姬鵬杰 杜坤 馮燕 周明 杜雨

摘要：暴雨頻率曲線擬合是推求暴雨強度公式必不可少的步驟，考慮經驗因素進行暴雨頻率曲線擬合，提出將暴雨強度頻率曲線擬合作為最優化問題，采用加權阻尼高斯牛頓迭代算法求解。與已有方法相比，提出引入權重系數以提高工程常用重現期段擬合精度，避免不同歷時暴雨頻率曲線相交；提出應用有限差分法簡化雅克比矩陣計算，并在海塞矩陣對角添加阻尼系數改進迭代收斂。以云南省保山市隆陽區33 a實測降雨資料為例，證明了算法的可行性及實用性。

關鍵詞：暴雨頻率曲線；經驗因素；最優化；迭代算法

中圖分類號：TU992.02 文獻標志碼：A文章編號：16744764（2018）02007706

收稿日期：20170314

基金項目：國家自然科學基金（51608242）；云南省人才培養計劃項目 （14118943）

作者簡介：姬鵬杰（1992），男，主要從事市政工程研究，Email：443838096@qq.com。

杜坤（通信作者），男，博士， Email：250977426@qq.com。

Received：20170314

Foundation item：National Natural Science Foundation of China（No.51608242）；Personnel Training Program of Yunnan Province （No.14118943）.

Author brief：Ji Pengjie（1992），main research interest： municipal engineering，Email：860655976@qq.com.

Du Kun （correspondence author）， PhD， Email：250977426@qq.com.Study on optimal fitting algorithm of rainstorm frequency

curve considering the empirical factors

Ji Pengjie1， Du Kun1， Feng Yan1， Zhou Ming1， Du Yu2

（1. Faculty of Civil Engineering and Mechanics， Kunming University of Science and Technology， Kunming 650500， P. R. China；

2. The Third Construction Engineer Company LTD， China Construction Second Engineer Bureau， Wuhan 430022， P. R. China）

Abstract：The rainstorm frequency curve fitting is essential for the identification of storm intensity formula， this paper carried out the study of rainstorm frequency curve fitting with considering of experience factors， and put forward to regard the rainstorm intensity frequency curve fitting as an optimization problem， and then to solve it by the weighted damped GaussNewton iterative algorithm. Compared to existing methods， the proposed method introduced weight coefficients to improve the fitting precision of commonly used recurrence period in engineering， and to avoid the intersection problem of different frequency curves. The finite difference method is proposed to simplify the calculation of Jacobian matrix， and the damping coefficient was added in Hesse matrix to improve iterative convergence. Thirtythree years of rainfall data of Longyang District of Baoshan city in Yunnan Province were used as an example to illustrate and demonstrate the feasibility and practicability of the proposed algorithm.

Keywords：rainstorm frequency curve； empirical factors； optimization； iterative algorithm

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

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

筆者將暴雨強度頻率曲線擬合作為最優化問題，提出加權阻尼高斯牛頓迭代算法進行求解。與已有方法相比，該算法引入權重系數提高工程常用重現期段擬合精度，通過調節權重系數避免不同歷時理論頻率曲線相交問題；再者，采用有限差分法簡化雅克比矩陣計算，并在海塞矩陣對角添加阻尼系數改進迭代收斂，避免了反復適線導致的巨大計算工作量。以云南省保山市隆陽區1981—2013年間33 a實測降雨資料為例，證明算法的可行性與實用性。

第2期 姬鵬杰，等：考慮經驗因素的暴雨頻率曲線最優化擬合算法1基于最優化的暴雨頻率曲線參數擬

合框架選擇PIII型曲線作為理論暴雨頻率曲線，其密度函數為y=βαΓ（α）（x-a0）α-1e-β（x-a0）（1）式中：α=4C2s；β=2（CvCs）；a0=（1-2CvCs）；Г（α）為伽瑪函數；x、Cv和Cs分別為均值、偏差系數及偏態系數。將參數擬合作為優化問題，目標函數可寫為f（x，Cs，Cv，tp）=（Xp-p）TW（Xp-p）（2）式中：Xp為實測暴雨強度，為n×1的列向量，n為數據采集年限；W為權重系數矩陣；p為對應的理論暴雨強度。在Matlab環境下，理論暴雨強度p可采用式（3）計算。p=h（x，Cv，Cs，Τp）=

1βgaminv（1-Τp，α，1）-a0（3）式中：Tp為降雨頻率，其中元素tp=m（n+1），m為實測暴雨強度從大到小排列的年序次。目標函數的物理意義是求取參數x、Cv和Cs使理論暴雨強度p與實測值Xp間均方差最小，并引入權重系數矩陣W提高工程常用段降雨頻率下擬合精度。由于理論暴雨強度p與參數x、Cv和Cs非線性相關，故需采用迭代法求解優化問題[9]。筆者采用阻尼高斯牛頓迭代法進行求解[10]，基本思路是采用參數的矩估計值作為初始值，通過計算雅克比矩陣構造搜索向量，沿目標函數減小方向修正參數，并在海塞矩陣對角添加阻尼系數改進迭代收斂性。為便于推導，記κ=[x，Cv，Cs]，第k次迭代的解為f（κk+Δκk，Τp）=[Xp-h（κk+Δκk，Τp）]T·

W[Xp-h（κk+Δκk，Τp）]（4）式（4）的一階線性展開式為f（κk+Δκk，Τp）=[ΔXkp-J（κk，Τp）Δκk]T·

W[ΔXkp-J（κk，Τp）Δκk]（5）式中：ΔXkp=Xp-h（κk，Τp），雅克比矩陣J（κ，Τp）=[pxpCvpCs]，為n×3矩陣。

由于理論暴雨強度計算涉及伽瑪函數和不完全伽瑪函數運算，使得雅克比矩陣解析式推導非常復雜。采用有限差分法計算雅克比矩陣，如式（6）所示。px=h（x+Δx，Cv，Cs，Tp）-h（x，Cv，Cs，Tp）Δx

pCv=h（x，Cv+ΔCv，Cs，Tp）-h（x，Cv，Cs，Tp）ΔCv

pCs=h（x，Cv，Cs+ΔCs，Tp）-h（x，Cv，Cs，Tp）ΔCs（6）理論上，在進行有限差分計算時，式（6）中Δ、ΔCv、ΔCs的取值越小，計算結果越接近解析解，但取值過小會超出計算機計算精度，反而導致不準確的計算結果[11]。經過嘗試，在差分過程中，推薦Δx=0.001，ΔCv=0.01，ΔCs=0.01。根據多元函數極值理論，當目標函數取得極小值時，應有f（κk+Δκk，Τp）Δκk=-2J（κk，Τp）T·

W[ΔXkp-J（κk，Τp）Δκk]=0（7）根據式（7）可得Δκk=[J（κk，Τp）TWJ（κk，Τp）]-1J（κk，Τp）TWΔXkp（8）為改進迭代收斂性，在海塞矩陣的對角添加阻尼系數矩陣Γ[12]，可得Δκk=[J（κk，Τp）TWJ（κk，Τp）+Γ]-1·

J（κk，Τp）TWΔXkp（9）在參數擬合過程中，可設置阻尼系數Γ的取值隨迭代次數的增加而減小，擬合框架如圖1所示。

圖1暴雨頻率曲線參數擬合框架圖

Fig. 12算例分析

算例分析旨在利用實際降雨數據闡明：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如圖2中曲線a所示，當設置權重矩陣W中所有元素為1時，即采用普通的高斯牛頓法求解優化問題[13]，此時殘差Δφ=Xp-p2n=0.088，整體擬合最佳，工程常用重現段擬合殘差Δφ=0062。為進一步提高工程常用段擬合精度，減小其他段權重系數并保持工程常用段權重系數不變，表1給出了調整權重系數時殘差變化情況。表1調整權重系數時殘差變化情況

Table 1權重系數工程常用段其它段整體殘差工程常用

段殘差對應曲線110.0880.062a10.50.0890.060b10.250.0900.059c10.062 50.0920.055d

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

2.2調節權重系數避免理論頻率曲線相交

采用傳統高斯牛頓迭代法進行參數擬合，即設置權重矩陣W中所有元素值為1，將各歷時適線結果繪制在同一海森機率格圖上。如圖3所示，20 min與30 min降雨歷時下理論頻率曲線出現相交趨勢，這有悖于暴雨強度隨歷時增大而減小這一基本前提，明顯不合理。由此可見，如果簡單地將暴雨頻率曲線擬合作為數學問題求解，可能導致理論頻率曲線相交這一不合理結果。針對該問題，可通過改變權重系數對適線結果進行微調，使理論頻率曲線不相交，如圖4所示。

圖3不同歷時下理論頻率曲線相交

Fig. 3圖4調節權重系數使頻率曲線不相交

Fig. 42.3添加阻尼系數改進迭代收斂

高斯牛頓迭代法的收斂性與初值選取相關，一般只有當初值比較靠近真值時才能保證迭代收斂[14]。以5 min降雨歷時下暴雨強度頻率曲線擬合為例，采用矩估計值作為參數x、Cv、Cs的初值，應用高斯牛頓迭代法求解優化問題。如圖5所示，迭代到第6步時，程序提示矩陣奇異、計算不精確，迭代發散、運算終止。

圖5普通高斯牛頓迭代法收斂情況

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

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

圖6添加阻尼系數改進迭代收斂

Fig.62.4基于差分進化算法的優化精度檢驗

優化問題的求解方法分為確定性算法與隨機搜索算法。采用基于梯度信息的確定性搜索算法，其特點是計算效率高，但可能陷入局部最優解。差分進化算法屬于隨機搜索類算法，其特點是計算量大，但通過多次運算能逼近全局最優解[17]。為檢驗所提出算法優化精度，通過多次運行差分進化算法進行對比。限于篇幅以5 min降雨歷時下暴雨強度數據為例，將差分進化算法運行100次，給出目標函數殘差最小的10個解與所提出算法解進行對比。其中，差分進化算法的種群規模取30，3個參數的搜索范圍均為[0 10]，變異因子F=0.6，交叉因子Cr=0.6，最大迭代次數G=100。兩種算法優化結果及目標函數殘差如表2，其中A算法為加權阻尼高斯牛頓迭代算法，B算法為差分進化算法。表2二種算法優化結果及目標函數殘差

Table 2 算法均值mCvCs目標函數殘差AB1.695 8820.284 8012.849 5360.088 2621.695 7580.284 8342.849 4380.088 2621.695 9360.284 8162.848 0230.088 2621.695 8510.284 6972.846 2250.088 2621.695 8550.284 8002.847 4490.088 2621.695 8340.284 8192.850 060.088 2621.695 8680.284 7692.847 7930.088 2621.695 8650.284 7622.848 2600.088 2621.695 868 0.284 758 2.847 301 0.088 262 1.695 948 0.284 763 2.848 609 0.088 262 1.695 806 0.284 841 2.848 003 0.088 262

如表2所示，為對比計算精度，將結果保留6位小數。其中加權阻尼高斯牛頓迭代法優化所得的目標函數殘差為0.008 826 2，差分進化算法逼近的全局最優解的目標函數殘差也為0.008 826 2。雖然參數優化結果有一定差異，但差異均小于10-3，從工程應用角度來看，上述差異可以忽略。因此，可以認為加權阻尼高斯牛頓迭代算法能獲得全局最優解。

3結論

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

1）在利用傳統的高斯牛頓法進行暴雨頻率曲線擬合時，迭代可能不收斂，且存在不同歷時暴雨頻率曲線相交的問題；

2）通過在海塞矩陣對角添加阻尼系數能保證高斯牛頓迭代法的收斂，建議采用阻尼系數二次方衰減進行暴雨頻率曲線擬合。

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

4）由于PIII型曲線涉及伽瑪函數和不完全伽瑪函數運算，雅克比矩陣解析式推導繁復，建議采用有限差分法簡化計算，推薦Δ=0.001，ΔCv=001，ΔCs=0.01。

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

參考文獻：

[1] 車伍， 楊正， 趙楊，等. 中國城市內澇防治與大小排水系統分析[J]. 中國給水排水， 2013， 29（16）： 1319.

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）：1319. （in Chinese）

[2] CHE G W. PearsonIII frequency curve plotting in Excel table [J]. Applied Mechanics & Materials， 2014， 556562： 58295834.

[3] 王正發. MATLAB在PⅢ型分布離均系數p值計算及頻率適線中的應用[J]. 西北水電， 2007（4）： 14.

WANG Z F. MATLAB programming language used to calculate variation coefficient p of the PIII distribution and fit a frequency curve [J]. Northwest Hydropower， 2007（4）： 14. （in Chinese）

[4] BLAIN G C. Standardized precipitation index based on pearson type III distribution [J]. Revista Brasileira De Meteorologia， 2014， 26（2）：167180.

[5] 崔俊蕊， 王政然， 梁爽，等. 城市設計暴雨頻率曲線的擬合及參數優化[J]. 水電能源科學， 2014（11）： 4851.

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）： 4851. （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）：517528.

[7] WU Y C， LIU J J， SU Y F， et al. Establishing acceptance regions for Lmoments based goodnessoffit tests for the Pearson type III distribution [J]. Stochastic Environmental Research and Risk Assessment， 2014， 26（6）： 873885.

[8] 高琳， 周玉文， 唐穎， 等. 城市暴雨強度公式皮爾遜Ⅲ型適線問題研究[J]. 給水排水， 2016（8）： 4751.

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）： 4751. （in Chinese）

[9] RAFIQ A， RAFIULLAH M. Some multistep iterative methods for solving nonlinear equations [J]. Computers & Mathematics with Applications， 2009， 58（8）： 15891597.

[10] FAIRBANK M， ALONSO E. Efficient calculation of the GaussNewton approximation of the Hessian matrix in neural networks [J]. Neural Computation， 2012， 24（3）： 607610.

[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）： 20812085.

[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）： 715749.

[13] ERINA M Y，IZMAILOV A F. The GaussNewton method for finding singular solutions of systems of nonlinear equations [J]. Anz Journal of Surgery， 2015， 58（2）： 395405.

[14] PHAN A H， TICHAVSKY P， CICHOCKI A. Low complexity damped GaussNewton algorithms for CANDECOMP /PARAFAC[J]. Siam Journal on Matrix Analysis& Applications， 2013， 34（1）： 126147.

[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）： 5777.

[16] JIN Q. Further convergence results on the general iteratively regularized GaussNewton methods under the discrepancy principle [J]. Mathematics of Computation， 2013， 82（283）： 16471665.

[17] MACIEL L， GOMIDE F， BALLINI R. A differential evolution algorithm for yield curve estimation [J]. Mathematics & Computers in Simulation， 2016， 129：1030.