Weibull分布參數(shù)的粒子群算法估計＊

董 勝，韓 意，陶山山，樊敦秋

（1.中國海洋大學工程學院，山東青島266100；2.勝利石油管理局鉆井工藝研究院，山東東營，257017）

在海洋技術、水利工程和城市建筑設計中，Weibull分布以其較強的適應性和靈活性被廣泛地運用于多年一遇極值風速、波高、水位、降水量等的長期預測中［1］。目前，Weibull分布的參數(shù)估計方法有矩法、極大似然法、概率權重矩法、相關系數(shù)法及Bayes估計法等［1－3］。嚴曉東等人［4］對三參數(shù)Weibull分布的多種估計方法進行比較研究，并指出各種方法的適用范圍。極大似然法需要迭代求解聯(lián)立的3個超越方程，計算相當復雜［3］，而且對初值的要求比較高［2］；概率權重矩法（PWM）不適合于子樣較少的情況，當子樣少于10個時，計算精度很差［3］；相關系數(shù)法的計算比較復雜［5］；Bayes方法將先驗分布和后驗分布結合起來［5］，但后驗分布比較復雜［3］，并且一般需要和其他方法聯(lián)合使用［5］。2001年鄭紅星和李麗娟將最速下降法應用于水質模型參數(shù)的優(yōu)化中［6］，并將結果與模擬退火算法和遺傳算法進行比較，得出最速下降法的計算結果易受初始搜索條件的影響，尋優(yōu)過程容易陷入局部最優(yōu)，而且精度低于非數(shù)值隨機優(yōu)化算法，但是最速下降法的計算速度明顯高于非數(shù)值隨機優(yōu)化算法。

近幾年，粒子群算法的研究引起了許多學者的關注。該算法簡潔、收斂速度快等優(yōu)勢被廣泛地應用于參數(shù)的非線性擬合中。江燕等［7］將其應用于新安江模型參數(shù)的優(yōu)選中。郭建青等［8］在求解河流水質參數(shù)的函數(shù)優(yōu)化問題時引入了該算法，取得了良好的效果。路志強等［9］運用該法對暴雨強度公式參數(shù)進行了優(yōu)化。郭建青等［10］研究了含水層參數(shù)的函數(shù)優(yōu)化問題，認為粒子群算法給出的結果優(yōu)良。羅航等［11］結合雙線性回歸估計和極大似然估計方法求解三參數(shù)Weibull分布的參數(shù)，即首先基于最小二乘法的雙線性回歸法求出初始解，再運用基于粒子群算法的極大似然估計法迭代求解。

上述成果說明粒子群算法正在被廣泛應用于實際工程中。本文將粒子群算法引入對Weibull分布的參數(shù)估計中。運用Weibull分布參數(shù)估計常用的3種方法，即最小二乘法、矩估計法和最速下降法，與新引入的粒子群算法，以廣西潿洲島海洋監(jiān)測站3個方向的年極值波高觀測數(shù)據(jù)為例，對Weibull分布的參數(shù)進行了估計，并對各種方法所得結果進行了比較分析。

1 Weibull分布參數(shù)的估計方法

三參數(shù)Weibull分布的分布函數(shù)可表示如下：

式中：a，b和c分別表示位置參數(shù)、尺度參數(shù)和形狀參數(shù)，且b＞0，c＞0。其對應的密度函數(shù)為

Weibull分布的參數(shù)估計方法有很多，比如最小二乘法、矩估計法、最速下降法等，在引入Weibull分布的粒子群算法前，本文將簡單介紹一下這些方法。以下假設x1，x2，…，xn為來自滿足式（1）的Weibull分布的獨立同分布樣本。

1.1 最小二乘法

在式（1）中令F（x）＝p，則式（1）變形為：

然后對式（3）變形，得

對式（4），令

則式（4）變形為AZ＋B＝Y，即為線性方程。由于在Z中a值未知，故首先對a循環(huán)賦值。假設x（1）≤x（2）≤…≤x（n）為樣本x1，x2，…，xn的順序樣本，假設已取定某一a值，令zi＝ln（x（i）－a），（i＝1，…，n），可以得到一列z1，z2，…，zn。另外假設對應于任一x（i）（i＝1，…，n）的累積頻率為i／（n＋1），根據(jù)式y(tǒng)i＝ln［－ln（1－i／（n＋1））］可以容易得到對應于順序樣本x（1）≤x（2）≤…≤x（n）的一列值y1，y2，…，yn。對兩列數(shù)據(jù)z1，z2，…，zn及y1，y2，…，yn運用最小二乘法，根據(jù)AZ＋B＝Y進行一次擬合，求出系數(shù)A，B，然后按照式（5），反求參數(shù)b，c。在遍歷所有滿足條件的a值的條件下，找出滿足離差平方和最小，即滿足的參數(shù)組a，b和c，則此時的參數(shù)a，b和c即為所求Weibull分布3個參數(shù)的最小二乘估計值。

1.2 矩估計法

Weibull分布的一階原點矩為：

Weibull分布的二階中心矩和三階中心距為：

其偏態(tài)系數(shù)為

1.3 最速下降法

在運籌學中已知負梯度方向為最速下降方向，所以在求解目標函數(shù)的非線性約束極小化問題時，可以使用負梯度方向作為尋優(yōu)方向［12］。也就是通過在負梯度方向的一維搜索，來確定使目標函數(shù)最小的步長，這種方法即為最速下降法。步驟如下：

（1）給定初始點X0＝（a0，b0，c0）以及精度η＞0，令k＝0，若梯度，則X0即為近似的極小點。

在Weibull分布參數(shù)的最速下降法求解中目標函數(shù)為

式中：yi可由i／（n＋1）估計。Weibull分布的最速下降法的梯度為

由于負梯度方向僅在一點附近才具有最速下降的性質，對于全局極小化求解過程，運用梯度法極易得到局部最優(yōu)解，所以初值的選取非常關鍵。

1.4 粒子群算法

粒子群算法是1種基于群體智能的算法，它在解決非線性問題上具有非常優(yōu)良的特性。在其迭代搜索過程中，通過局部最優(yōu)解pbest和全局最優(yōu)解gbest來動態(tài)調整自己的位置及速度，最終完成尋優(yōu)。

粒子群算法對初值的要求不是很高。使用以下迭代公式對粒子的速度和位置進行迭代［13］：

用粒子群算法求解Weibull分布參數(shù)的估計，可以歸結為以下優(yōu)化問題：

式中：N表示粒子個數(shù)。初始粒子的速度為：

更新粒子速度，計算V2，按照公式

更新粒子位置，計算X2。檢驗是否滿足或者k≥Nmax（Nmax為最大迭代次數(shù)）。若滿足，則此時的gbest1即為擬合Weibull分布的a，b，c 3個參數(shù)；若不滿足則令k＝k＋1，同時按照式（13）進行迭代，按照公式

圖1 粒子群算法流程圖Fig.1 Flowchart of particle swarm optimization

2 實例分析

廣西沿海潿洲島海洋監(jiān)測站3個方向（NNE、E、SSW）1960－1993年的年最大波高見圖2～4［14］。

用三參數(shù)Weibull分布對圖2～4的數(shù)據(jù)進行擬合，采用最小二乘法、矩估計法、最速下降法及粒子群算法分別對Weibull分布的未知參數(shù)進行了估計。

圖2 潿洲島站NNE方向年最大波高曲線Fig.2 Curve of annual maximum wave heights in NNE at Weizhoudao Station

對于矩法，先根據(jù)各樣本數(shù)據(jù)求出樣本的均值，方差和偏度系數(shù)。分別利用式（9），（7）和（6），將其中的總體統(tǒng)計量用樣本統(tǒng)計量代換，得到Weibull分布3個參數(shù)的解。

在應用粒子群算法時，c1取0.9，c2取0.1，w取0.6，令最大迭代次數(shù)Nmax＝100，并且要求前后2次全局最優(yōu)代入Weibull分布后概率結果的離差平方和小于η＝10－7，取群體中粒子個數(shù)popsize＝1 000。由于分布中含有3個未知參數(shù)，所以群體維數(shù)dimsize＝3。

通過以上方法，擬合Weibull分布中的a，b，c三參數(shù)，最終3個方向（NNE，E，SSW）的結果如表1～3。其中表中最速下降法1和最速下降法2是指對于初值不同的最速下降方法的2次計算結果，粒子群算法1和粒子群算法2是指對于沒有改變任何參數(shù)情況下重復試驗后的不同的2次結果。

圖5～7是采用最小二乘法、矩估計法、最速下降法及粒子群算法對潿洲島3個方向年極值波高進行Weibull分布擬合的結果。由圖可見，實測數(shù)據(jù)分布在理論曲線兩側，擬合結果較好。

比較表1～3及圖5～7，可以看出：

（1）采用的4種擬合方法中，粒子群算法的概率離差平方和最小。針對本算例，粒子群算法的精度高于最小二乘法，矩估計法，以及最速下降法。

（2）矩估計法和最小二乘法在擬合分布的參數(shù)時無需迭代初值，粒子群算法對初值的要求不高，而最速下降法對初值有較高的要求。表1～3中的數(shù)據(jù)說明：即使初值的選擇只有微小的差別，最速下降法得到的結果仍存在很大差異。采用最速下降法實際計算時，若初值設定得不好，計算容易陷入局部最優(yōu)。

（3）由于粒子群算法在生成初始矩陣時具有隨機性，因此，擬合結果不穩(wěn)定，然而，每次結果相差不大，且與其他3種算法相比，發(fā)生概率的離差平方和最小，對觀測數(shù)據(jù)的擬合效果更優(yōu)。

（4）最速下降法與粒子群算法都需要給定初始迭代值，二者相比，最速下降法耗時更大。最小二乘法和矩估計法無需提供初始迭代值。粒子群算法與它們比較，能夠看出：粒子群算法明顯快于最小二乘法，但是稍慢于矩估計法。

圖7 潿洲島站SSW向極值波高Weibull分布Fig.7 Weibull distribution of extreme H in SSW at Weizhoudao station

表1 潿洲島站NNE向年極值波高擬合結果Table 1 Fitting results of maximum wave height in NNE at Weizhoudao station

表2 潿洲島站E向年極值波高擬合結果Table 2 Fitting results of maximum wave height in E at Weizhoudao station

表3 潿洲島站SSW向年極值波高擬合結果Table 3 Fitting results of maximum wave height in SSW at Weizhoudao station

3 結語

本文將粒子群算法應用于Weibull分布的參數(shù)擬合中，提高了分布參數(shù)擬合的速度?；趶V西沿海潿洲島海洋監(jiān)測站3個方向（NNE、E、SSW）1960—1993年的年最大波高數(shù)據(jù)，采用粒子群算法，最小二乘法，矩估計法，以及最速下降法4種方法，擬合了Weibull分布中的參數(shù)。計算結果表明：粒子群算法具有對初值要求不高，計算速度快，計算耗時不高等優(yōu)點。綜上可知，粒子群算法是海洋工程環(huán)境要素統(tǒng)計中估計Weibull分布參數(shù)的1種有效方法。

［1］ Goda Y.Random seas and design of maritime structures［M］.Singapore：World Scientific，2010.

［2］ 徐以鋒，翁永剛，王明星，等.威布爾分布三參數(shù)的擬合特性［J］.鄭州大學學報：理學版，2003，35（3）：30－34.

［3］ 湯銀才，侯道燕.三參數(shù)Weibull分布參數(shù)的Bayes估計［J］.系統(tǒng)科學與數(shù)學，2009，29（1）：109－115.

［4］ 嚴曉東，鄭榮越，吳亮，等.三參數(shù)威布爾分布參數(shù)估計方法比較［J］.寧波大學學報：理工版，2005，18（3）：301－305.

［5］ 胡恩平，羅興柏，劉國慶.三參數(shù)Weibull分布幾種常用的參數(shù)估計方法［J］.沈陽工業(yè)學院學報，2000，19（3）：88－94.

［6］ 鄭紅星，李麗娟.水質模型參數(shù)的非數(shù)值隨機優(yōu)化［J］.地理研究，2001，20（1）：97－102.

［7］ 江燕，胡鐵松，武夏寧，等.粒子群算法在新安江模型參數(shù)優(yōu)選中的應用［J］.武漢大學學報：工學版，2006，39（4）：14－17.

［8］ 郭建青，王洪勝，李彥，等.粒子群優(yōu)化算法在確定河流水質參數(shù)中的應用［J］.水利水電科技進展，2007，27（6）：1－5.

［9］ 路志強，楊文海，康迎賓，等.粒子群算法及其在暴雨強度公式參數(shù)優(yōu)化中的應用［J］.南水北調與水利科技，2008，6（3）：72－73.

［10］ 郭建青，王洪勝，李彥，等.粒子群優(yōu)化算法在確定河流水質參數(shù)中的應用［J］.水利水電科技進展，2007，27（6）：1－5.

［11］ 羅航，王厚軍，黃建國，等.基于PSO的三參數(shù)威布爾分布參數(shù)的聯(lián)合估計方法［J］.儀器儀表學報，2009，30（1）：160－161.

［12］ 甘應愛，田豐，李維錚，等.運籌學［M］.北京：清華大學出版社，2005：151－155.

［13］ 紀震，廖惠連，吳青華.粒子群算法及應用［M］.北京：科學出版社，2009：16－21.

［14］ 夏華永，李樹華.廣西沿海年極值波高分析［J］.熱帶海洋學報，2001，20（2）：1－6.