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

（1.中國(guó)海洋大學(xué)工程學(xué)院，山東青島266100；2.勝利石油管理局鉆井工藝研究院，山東東營(yíng)，257017）

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

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

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

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

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

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

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

1.1 最小二乘法

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

然后對(duì)式（3）變形，得

對(duì)式（4），令

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

1.2 矩估計(jì)法

Weibull分布的一階原點(diǎn)矩為：

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

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

1.3 最速下降法

在運(yùn)籌學(xué)中已知負(fù)梯度方向?yàn)樽钏傧陆捣较?，所以在求解目?biāo)函數(shù)的非線性約束極小化問(wèn)題時(shí)，可以使用負(fù)梯度方向作為尋優(yōu)方向［12］。也就是通過(guò)在負(fù)梯度方向的一維搜索，來(lái)確定使目標(biāo)函數(shù)最小的步長(zhǎng)，這種方法即為最速下降法。步驟如下：

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

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

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

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

1.4 粒子群算法

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

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

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

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

更新粒子速度，計(jì)算V2，按照公式

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

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

2 實(shí)例分析

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

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

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

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

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

通過(guò)以上方法，擬合Weibull分布中的a，b，c三參數(shù)，最終3個(gè)方向（NNE，E，SSW）的結(jié)果如表1～3。其中表中最速下降法1和最速下降法2是指對(duì)于初值不同的最速下降方法的2次計(jì)算結(jié)果，粒子群算法1和粒子群算法2是指對(duì)于沒(méi)有改變?nèi)魏螀?shù)情況下重復(fù)試驗(yàn)后的不同的2次結(jié)果。

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

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

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

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

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

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

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

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

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

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

3 結(jié)語(yǔ)

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

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

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

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

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

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

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

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

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

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

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

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

［12］ 甘應(yīng)愛(ài)，田豐，李維錚，等.運(yùn)籌學(xué)［M］.北京：清華大學(xué)出版社，2005：151－155.

［13］ 紀(jì)震，廖惠連，吳青華.粒子群算法及應(yīng)用［M］.北京：科學(xué)出版社，2009：16－21.

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