雷慶文，閆 磊，魯東陽，卜嘉駿，羅 云

（1.河北工程大學(xué)水利水電學(xué)院，河北 邯鄲 056038；2.河北工程大學(xué)河北省智慧水利重點實驗室，河北 邯鄲 056038）

0 引 言

水文頻率分析是水利工程中的一個重要環(huán)節(jié)，其計算過程主要是對分布線型的選擇以及統(tǒng)計參數(shù)的估計。1980年以來我國制定的不同版本的水利水電工程水文計算規(guī)范中均規(guī)定皮爾遜III型曲線（P-III）為備選分布線型［1］。確定分布線型后，統(tǒng)計參數(shù)估計的好壞將直接影響到設(shè)計值的合理性。

水文統(tǒng)計中認(rèn)為極大似然法（Maximum Likelihood Estima‐tion，MLE）有良好的參數(shù)估計能力。然而，采用MLE 估計P-III型分布參數(shù)時，解析求解困難，參數(shù)估計結(jié)果不靈敏，而且似然方程還存在無解等問題都很大程度限制了其在實際水文頻率分析中的應(yīng)用［2］。前人對MLE 估計水文分布參數(shù)進(jìn)行了有益的探索，周玉文等［3］利用極大似然法研究了杭州市暴雨強(qiáng)度統(tǒng)計資料，證明采用極大似然法估計的P-Ⅲ型分布較矩估計法與適線法估計的分布擬合程度更好；余泱悅等［4］借助MATLAB軟件分析黃家港觀測站年平均徑流數(shù)據(jù)認(rèn)為極大似然法具有很好的擬合效果；Yan 等［5］為估計時變兩組分混合分布模型中眾多的分布參數(shù)，提出了一種結(jié)合模擬退火算法和極大似然的元啟發(fā)式參數(shù)估計方法；劉桐愷［6］在對福建地區(qū)降雨頻率分析的研究中，將似然方程中無法解析求解的Digamma 函數(shù)通過斯特林展開計算其近似值。極大似然法在P-III 線型參數(shù)估計中的應(yīng)用目前還沒得到有效解決，相關(guān)理論研究并不成熟，因此，在實際的水文頻率計算中使用較少［7］。Python語言中的科學(xué)計算庫Scipy 提供了P-III 型分布極大似然估計的程序，但在實際的應(yīng)用中問題較多。為將極大似然法有效地應(yīng)用于水文頻率分析中，采用一種結(jié)合粒子群優(yōu)化算法（Particle Swarm Optimi‐zation，PSO）與極大似然法的PSO-MLE 算法，通過數(shù)值計算求解P-III型分布參數(shù)的極大似然估計值，以實例分析和蒙特卡羅統(tǒng)計模擬探究此算法的參數(shù)估計效果。

1 理論與方法

1.1 傳統(tǒng)極大似然法

傳統(tǒng)極大似然法通過構(gòu)造所選分布線型的對數(shù)似然函數(shù)，利用函數(shù)取極值的必要條件計算出參數(shù)的極大似然估計值。

公式（1）中L(θ)為將觀測序列xi(i=1，2，…，n)代入所選分布線型f(x；θ)得到的似然函數(shù)，通常為計算方便，構(gòu)造對數(shù)似然函數(shù)Gf，此函數(shù)即為極大似然估計的優(yōu)化目標(biāo)函數(shù)。

公式（2）中的α、β、a0分別對應(yīng)P-III 線型概率密度函數(shù)的形狀參數(shù)、尺度參數(shù)和位置參數(shù)。實際應(yīng)用中更多的是以均值EX、偏態(tài)系數(shù)Cs和變差系數(shù)Cv表示這三個參數(shù)。

P-III型分布的極大似然估計必要條件方程組：

公式（4）為似然方程，是含Digamma 函數(shù)的非線性方程組，解析求解非常困難，在工程應(yīng)用中需要試算，計算繁瑣，參數(shù)估計精度不高。而且當(dāng)≤1（Cs≥2）時，其中3 式第一項≥0，第二項顯然也大于0，因此，似然方程還存在無解問題。傳統(tǒng)極大似然法的缺陷使得MLE在P-III型分布參數(shù)估計中很少使用。

1.2 粒子群極大似然法

概率分布模型參數(shù)的極大似然估計本質(zhì)上是參數(shù)優(yōu)化問題，鑒于PSO 在參數(shù)優(yōu)化中的良好應(yīng)用效果，將其作為研究PIII型分布極大似然估計的優(yōu)化算法。

由于P-III 型分布的對數(shù)似然函數(shù)對參數(shù)變化的反應(yīng)不夠靈敏，粒子尋優(yōu)過程中容易局部收斂。為此，通過動態(tài)調(diào)整慣性權(quán)重ω來提高算法的局部尋優(yōu)能力，周期性更新ω促使粒子可以跳出局部最優(yōu)解，提高算法的收斂精度。PSO-MLE法的適應(yīng)度函數(shù)為公式（1）表示的對數(shù)似然函數(shù)，計算流程見圖1。

圖1 PSO-MLE 計算流程圖Fig.1 Flowchart of PSO-MLE method

P-III 型分布的PSO-MLE 法粒子搜索空間為3 維，慣性權(quán)重采用一種聯(lián)合進(jìn)化離散度κ(t)和Sigmoid 函數(shù)的非線性動態(tài)自適應(yīng)慣性權(quán)重［8］。

將第t代種群與第t-1代種群的適應(yīng)度值標(biāo)準(zhǔn)差之比定義為進(jìn)化離散度κ(t)，用來描述種群進(jìn)化過程中適應(yīng)度值的變化。

公式（6）中，ωmax、ωmin、ω(t)分別為最大、最小、t代種群慣性權(quán)重；Tmax為最大迭代次數(shù)；b為阻尼因子。大量實驗表明ωmax=0.9、ωmin=0.4、Tmax=100、b=0.5時計算效果較好。

1.3 矩法與優(yōu)化適線法

矩法估計（MOM）構(gòu)造估計量的原理與方法簡單，應(yīng)用方便，而且，適用于不同的總體分布?？梢韵壤镁胤ü烙媴?shù)的大致區(qū)間，便于選取粒子群的優(yōu)化空間。

優(yōu)化適線法（FIT）是基于離差平方和、離差絕對值和以及相對離差平方和最小準(zhǔn)則，構(gòu)建水文頻率參數(shù)優(yōu)化模型。最小二乘得出的曲線通過點群中心，所得參數(shù)與目估適線法的結(jié)果比較接近。因此，文中的優(yōu)化適線法以離差平方和為優(yōu)度指標(biāo)求解P-III線型“累計頻率曲線”的最小二乘擬合［9］。

公式（7）中pm為實測水文序列的經(jīng)驗頻率，為累計概率密度函數(shù)的反函數(shù)，以實測序列求解的最小二乘擬合。

2 蒙特卡羅模擬

為證明PSO-MLE 法估計參數(shù)的統(tǒng)計特性，采用蒙特卡羅模擬（Monte Carlo Simulation）檢驗不同方法所得參數(shù)估計值的無偏性和有效性［公式（8）］。

實際水文頻率計算時，序列通常不長，因此，采樣長度選擇80。本研究僅對形狀參數(shù)進(jìn)行改變，分別考慮Cs＜＜Cs＜2、2＜Cs＜∞三種形狀線型，Cs/Cv=3 和均值EX=100 保持不變，隨機(jī)模擬次數(shù)N=1 000次。

從蒙特卡羅的模擬結(jié)果來看，PSO-MLE法的參數(shù)估計效果要明顯優(yōu)于其他方法，而且可有效解決Python提供的MLE 算法對Cs≥2線型的參數(shù)估計表現(xiàn)不佳的問題。

似然方程無解（Cs≥2）時，對數(shù)似然函數(shù)并不存在嚴(yán)格意義上的極大似然值。由于α≤1 時，?ln L ?a0＞0，對數(shù)似然函數(shù)隨a0遞增，參數(shù)a0的估計值是計算機(jī)浮點數(shù)精度允許下最大逼近實測序列中最小值的值。因此，a0對其他參數(shù)的估計結(jié)果影響很大，對數(shù)似然函數(shù)此時

顯然無法衡量模型的擬合優(yōu)劣（并非對數(shù)似然函數(shù)越大模型擬合越好），表1中Cs=3的模擬結(jié)果表明Python中的MLE程序并沒有很好地解決這個問題。雖然，不少學(xué)者認(rèn)為這種分布線型不宜在水文中應(yīng)用，但為確保PSO-MLE 算法的普遍適用性，提出一種改進(jìn)的極大似然估計算法：表1中矩法對均值的估計效果較好，而Cs≥2 時，Python 中的MLE 算法對均值估計很差。改進(jìn)后的算法在粒子優(yōu)化參數(shù)的過程中，直接使用矩法估計均值EX，只對Cs和Cv進(jìn)行優(yōu)化，此時，對數(shù)似然函數(shù)值雖然不是最大，但效果卻更好。

表1 參數(shù)估計的優(yōu)良性評價Tab.1 Optimal evaluation of parameter estimation

當(dāng)Cs＞2 時，即α＜1 的概率密度曲線呈特殊的乙字形，公式（9）表明如果參數(shù)a0的估計值為觀測樣本中的最小值，那么無論其他參數(shù)取何值，對數(shù)似然函數(shù)都將趨于無窮，這顯然是不合理的。因此，在粒子參數(shù)優(yōu)化過程中，引入δ約束參數(shù)：δ=10-λ（λ=2，…，5）。一般水文頻率計算中λ對結(jié)果影響不大，但當(dāng)出現(xiàn)Cs≥2時，可通過對比不同λ模型的均方誤差，選擇適當(dāng)?shù)摩酥担谑莂0的取值為：a0≤(1-δ) × min{xi}。PSO-MLE算法只對Cs和Cv優(yōu)化，粒子的優(yōu)化參數(shù)為：Cs和K=，由公式（3）得粒子搜索空間范圍：的倍比系數(shù)K為[Kmin，10Kmin]，偏態(tài)系數(shù)Cs根據(jù)矩法估計控制為[0.5 ×Cs，5 ×Cs]。PSO-MLE 法具有良好的全局搜索能力，應(yīng)用中若發(fā)現(xiàn)粒子收斂在約束空間邊界，則可繼續(xù)擴(kuò)大搜索范圍。

3 實例分析

為更好地檢驗PSO-MLE 算法在實際水文頻率分析中的應(yīng)用效果，本研究以渭河流域1951-2011年61 a間的降雨資料，通過面積加權(quán)平均計算得到的年總降雨量數(shù)據(jù)為例，分別利用PSO-MLE法、FIT和MOM計算參數(shù)估計值。

表2 參數(shù)估計結(jié)果Tab.2 Parameter estimation results

3.1 參數(shù)估計結(jié)果檢驗

3.1.1 必要條件檢驗

檢驗PSO-MLE 法參數(shù)估計結(jié)果是否滿足極大似然的必要條件，可進(jìn)一步驗證算法的有效性。將求解的參數(shù)代入公式（4），可以證明數(shù)值計算結(jié)果滿足極大似然的必要條件。

當(dāng)參數(shù)估計結(jié)果滿足極大似然必要條件時，由必要條件方程組公式（4）中的2式可以證明極大似然法獲得的P-III型分布均值估計量是無偏估計量。此時，矩法估計的均值和極大似然法相同。

圖2 不同方法的參數(shù)估計結(jié)果Fig.2 Parameter estimation results using different methods

3.1.2 擬合優(yōu)度檢驗

最小信息量準(zhǔn)則（AIC）、均方根誤差（RMSE）和Filliben相關(guān)系數(shù)能定量評價模型參數(shù)估計的合理性。RMSE與AIC越小，F(xiàn)illiben相關(guān)系數(shù)越大則模型擬合優(yōu)度越好［10］。模型殘差通常服從正態(tài)分布，可對殘差序列進(jìn)行正態(tài)檢驗，利用正態(tài)檢驗QQ圖可以判定模型選擇是否合理［12］。

圖3 正態(tài)檢驗Q-Q圖Fig.3 Q-Q diagram of normal test

表3 P-III分布參數(shù)估計誤差比較Tab.3 Comparison of parameter estimation error of P-III distribution

3.2 討 論

在實際的水文頻率分析中，F(xiàn)IT 應(yīng)用最廣，因為其對洪水頻率計算中的上部點據(jù)有較好的擬合效果，但應(yīng)明確的是這種擬合效果是建立在經(jīng)驗頻率假設(shè)的基礎(chǔ)上。關(guān)于計算經(jīng)驗頻率的相關(guān)理論很多，我國常用的經(jīng)驗公式假定也未必合理。而且，F(xiàn)IT 的Cs通常偏大，這勢必會造成設(shè)計值過大，影響工程建設(shè)的經(jīng)濟(jì)效益。PSO-MLE法充分考慮總體分布信息，參數(shù)估計結(jié)果統(tǒng)計特性更好，所得設(shè)計值更為合理。當(dāng)然，極大似然法在目前的工程設(shè)計中使用較少，更多的問題還有待實踐的檢驗。

4 結(jié) 論

（1）基于PSO 的MLE 算法可將極大似然法有效地應(yīng)用于水文頻率分析中，而且參數(shù)估計效果具有較好的無偏性和有效性。PSO-MLE法在實際的水文頻率分析中，具有很好的借鑒和使用意義，此方法也適用于水文頻率分析中的其他分布線型，只需將適應(yīng)度函數(shù)換作所需分布的對數(shù)似然函數(shù)。

（2）針對Cs≥2的分布線型，提出了一種直接利用矩法估計的均值和引入δ參數(shù)約束位置參數(shù)a0的改進(jìn)極大似然算法，取得了較好的應(yīng)用效果。相比于Python語言中提供的MLE 程序，PSO-MLE 算法具有更高的理論價值和實用價值。文中僅對比了FIT和MOM，還需進(jìn)一步與其他參數(shù)估計方法的比較。