基于粒子群算法的P-III型分布極大似然估計研究

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

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

0 引 言

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

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

1 理論與方法

1.1 傳統極大似然法

傳統極大似然法通過構造所選分布線型的對數似然函數，利用函數取極值的必要條件計算出參數的極大似然估計值。

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

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

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

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

1.2 粒子群極大似然法

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

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

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

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

將第t代種群與第t-1代種群的適應度值標準差之比定義為進化離散度κ(t)，用來描述種群進化過程中適應度值的變化。

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

1.3 矩法與優化適線法

矩法估計（MOM）構造估計量的原理與方法簡單，應用方便，而且，適用于不同的總體分布。可以先利用矩法估計參數的大致區間，便于選取粒子群的優化空間。

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

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

2 蒙特卡羅模擬

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

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

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

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

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

表1 參數估計的優良性評價Tab.1 Optimal evaluation of parameter estimation

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

3 實例分析

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

表2 參數估計結果Tab.2 Parameter estimation results

3.1 參數估計結果檢驗

3.1.1 必要條件檢驗

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

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

圖2 不同方法的參數估計結果Fig.2 Parameter estimation results using different methods

3.1.2 擬合優度檢驗

最小信息量準則（AIC）、均方根誤差（RMSE）和Filliben相關系數能定量評價模型參數估計的合理性。RMSE與AIC越小，Filliben相關系數越大則模型擬合優度越好［10］。模型殘差通常服從正態分布，可對殘差序列進行正態檢驗，利用正態檢驗QQ圖可以判定模型選擇是否合理［12］。

圖3 正態檢驗Q-Q圖Fig.3 Q-Q diagram of normal test

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

3.2 討 論

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

4 結 論

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

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