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

姬鵬杰 杜坤 馮燕 周明 杜雨

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

關(guān)鍵詞：暴雨頻率曲線；經(jīng)驗(yàn)因素；最優(yōu)化；迭代算法

中圖分類號(hào)：TU992.02 文獻(xiàn)標(biāo)志碼：A文章編號(hào)：16744764（2018）02007706

收稿日期：20170314

基金項(xiàng)目：國家自然科學(xué)基金（51608242）；云南省人才培養(yǎng)計(jì)劃項(xiàng)目 （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

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

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

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

第2期 姬鵬杰，等：考慮經(jīng)驗(yàn)因素的暴雨頻率曲線最優(yōu)化擬合算法1基于最優(yōu)化的暴雨頻率曲線參數(shù)擬

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

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

由于理論暴雨強(qiáng)度計(jì)算涉及伽瑪函數(shù)和不完全伽瑪函數(shù)運(yùn)算，使得雅克比矩陣解析式推導(dǎo)非常復(fù)雜。采用有限差分法計(jì)算雅克比矩陣，如式（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）理論上，在進(jìn)行有限差分計(jì)算時(shí)，式（6）中Δ、ΔCv、ΔCs的取值越小，計(jì)算結(jié)果越接近解析解，但取值過小會(huì)超出計(jì)算機(jī)計(jì)算精度，反而導(dǎo)致不準(zhǔn)確的計(jì)算結(jié)果[11]。經(jīng)過嘗試，在差分過程中，推薦Δx=0.001，ΔCv=0.01，ΔCs=0.01。根據(jù)多元函數(shù)極值理論，當(dāng)目標(biāo)函數(shù)取得極小值時(shí)，應(yīng)有f（κk+Δκk，Τp）Δκk=-2J（κk，Τp）T·

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

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

圖1暴雨頻率曲線參數(shù)擬合框架圖

Fig. 12算例分析

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

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

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

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

Fig. 2如圖2中曲線a所示，當(dāng)設(shè)置權(quán)重矩陣W中所有元素為1時(shí)，即采用普通的高斯牛頓法求解優(yōu)化問題[13]，此時(shí)殘差Δφ=Xp-p2n=0.088，整體擬合最佳，工程常用重現(xiàn)段擬合殘差Δφ=0062。為進(jìn)一步提高工程常用段擬合精度，減小其他段權(quán)重系數(shù)并保持工程常用段權(quán)重系數(shù)不變，表1給出了調(diào)整權(quán)重系數(shù)時(shí)殘差變化情況。表1調(diào)整權(quán)重系數(shù)時(shí)殘差變化情況

Table 1權(quán)重系數(shù)工程常用段其它段整體殘差工程常用

段殘差對(duì)應(yīng)曲線110.0880.062a10.50.0890.060b10.250.0900.059c10.062 50.0920.055d

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

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

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

圖3不同歷時(shí)下理論頻率曲線相交

Fig. 3圖4調(diào)節(jié)權(quán)重系數(shù)使頻率曲線不相交

Fig. 42.3添加阻尼系數(shù)改進(jìn)迭代收斂

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

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

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

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

圖6添加阻尼系數(shù)改進(jìn)迭代收斂

Fig.62.4基于差分進(jìn)化算法的優(yōu)化精度檢驗(yàn)

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

Table 2 算法均值mCvCs目標(biāo)函數(shù)殘差A(yù)B1.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所示，為對(duì)比計(jì)算精度，將結(jié)果保留6位小數(shù)。其中加權(quán)阻尼高斯牛頓迭代法優(yōu)化所得的目標(biāo)函數(shù)殘差為0.008 826 2，差分進(jìn)化算法逼近的全局最優(yōu)解的目標(biāo)函數(shù)殘差也為0.008 826 2。雖然參數(shù)優(yōu)化結(jié)果有一定差異，但差異均小于10-3，從工程應(yīng)用角度來看，上述差異可以忽略。因此，可以認(rèn)為加權(quán)阻尼高斯牛頓迭代算法能獲得全局最優(yōu)解。

3結(jié)論

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

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

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

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

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

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

參考文獻(xiàn)：

[1] 車伍， 楊正， 趙楊，等. 中國城市內(nèi)澇防治與大小排水系統(tǒng)分析[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] 王正發(fā). MATLAB在PⅢ型分布離均系數(shù)p值計(jì)算及頻率適線中的應(yīng)用[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] 崔俊蕊， 王政然， 梁爽，等. 城市設(shè)計(jì)暴雨頻率曲線的擬合及參數(shù)優(yōu)化[J]. 水電能源科學(xué)， 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] 高琳， 周玉文， 唐穎， 等. 城市暴雨強(qiáng)度公式皮爾遜Ⅲ型適線問題研究[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.