Weibull分布参数的粒子群算法估计＊

2012-01-10陶山山樊敦秋

中国海洋大学学报（自然科学版） 2012年6期

董胜，韩意，陶山山，樊敦秋

（1.中国海洋大学工程学院，山东青岛266100；2.胜利石油管理局钻井工艺研究院，山东东营，257017）

在海洋技术、水利工程和城市建筑设计中，Weibull分布以其较强的适应性和灵活性被广泛地运用于多年一遇极值风速、波高、水位、降水量等的长期预测中［1］。目前，Weibull分布的参数估计方法有矩法、极大似然法、概率权重矩法、相关系数法及Bayes估计法等［1－3］。严晓东等人［4］对三参数Weibull分布的多种估计方法进行比较研究，并指出各种方法的适用范围。极大似然法需要迭代求解联立的3个超越方程，计算相当复杂［3］，而且对初值的要求比较高［2］；概率权重矩法（PWM）不适合于子样较少的情况，当子样少于10个时，计算精度很差［3］；相关系数法的计算比较复杂［5］；Bayes方法将先验分布和后验分布结合起来［5］，但后验分布比较复杂［3］，并且一般需要和其他方法联合使用［5］。2001年郑红星和李丽娟将最速下降法应用于水质模型参数的优化中［6］，并将结果与模拟退火算法和遗传算法进行比较，得出最速下降法的计算结果易受初始搜索条件的影响，寻优过程容易陷入局部最优，而且精度低于非数值随机优化算法，但是最速下降法的计算速度明显高于非数值随机优化算法。

近几年，粒子群算法的研究引起了许多学者的关注。该算法简洁、收敛速度快等优势被广泛地应用于参数的非线性拟合中。江燕等［7］将其应用于新安江模型参数的优选中。郭建青等［8］在求解河流水质参数的函数优化问题时引入了该算法，取得了良好的效果。路志强等［9］运用该法对暴雨强度公式参数进行了优化。郭建青等［10］研究了含水层参数的函数优化问题，认为粒子群算法给出的结果优良。罗航等［11］结合双线性回归估计和极大似然估计方法求解三参数Weibull分布的参数，即首先基于最小二乘法的双线性回归法求出初始解，再运用基于粒子群算法的极大似然估计法迭代求解。

上述成果说明粒子群算法正在被广泛应用于实际工程中。本文将粒子群算法引入对Weibull分布的参数估计中。运用Weibull分布参数估计常用的3种方法，即最小二乘法、矩估计法和最速下降法，与新引入的粒子群算法，以广西涠洲岛海洋监测站3个方向的年极值波高观测数据为例，对Weibull分布的参数进行了估计，并对各种方法所得结果进行了比较分析。

1 Weibull分布参数的估计方法

三参数Weibull分布的分布函数可表示如下：

式中：a，b和c分别表示位置参数、尺度参数和形状参数，且b＞0，c＞0。其对应的密度函数为

Weibull分布的参数估计方法有很多，比如最小二乘法、矩估计法、最速下降法等，在引入Weibull分布的粒子群算法前，本文将简单介绍一下这些方法。以下假设x1，x2，…，xn为来自满足式（1）的Weibull分布的独立同分布样本。

1.1 最小二乘法

在式（1）中令F（x）＝p，则式（1）变形为：

然后对式（3）变形，得

对式（4），令

则式（4）变形为AZ＋B＝Y，即为线性方程。由于在Z中a值未知，故首先对a循环赋值。假设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），根据式yi＝ln［－ln（1－i／（n＋1））］可以容易得到对应于顺序样本x（1）≤x（2）≤…≤x（n）的一列值y1，y2，…，yn。对两列数据z1，z2，…，zn及y1，y2，…，yn运用最小二乘法，根据AZ＋B＝Y进行一次拟合，求出系数A，B，然后按照式（5），反求参数b，c。在遍历所有满足条件的a值的条件下，找出满足离差平方和最小，即满足的参数组a，b和c，则此时的参数a，b和c即为所求Weibull分布3个参数的最小二乘估计值。

1.2 矩估计法

Weibull分布的一阶原点矩为：

Weibull分布的二阶中心矩和三阶中心距为：

其偏态系数为

1.3 最速下降法

在运筹学中已知负梯度方向为最速下降方向，所以在求解目标函数的非线性约束极小化问题时，可以使用负梯度方向作为寻优方向［12］。也就是通过在负梯度方向的一维搜索，来确定使目标函数最小的步长，这种方法即为最速下降法。步骤如下：

（1）给定初始点X0＝（a0，b0，c0）以及精度η＞0，令k＝0，若梯度，则X0即为近似的极小点。

在Weibull分布参数的最速下降法求解中目标函数为

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

由于负梯度方向仅在一点附近才具有最速下降的性质，对于全局极小化求解过程，运用梯度法极易得到局部最优解，所以初值的选取非常关键。

1.4 粒子群算法

粒子群算法是1种基于群体智能的算法，它在解决非线性问题上具有非常优良的特性。在其迭代搜索过程中，通过局部最优解pbest和全局最优解gbest来动态调整自己的位置及速度，最终完成寻优。

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

用粒子群算法求解Weibull分布参数的估计，可以归结为以下优化问题：

式中：N表示粒子个数。初始粒子的速度为：

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

更新粒子位置，计算X2。检验是否满足或者k≥Nmax（Nmax为最大迭代次数）。若满足，则此时的gbest1即为拟合Weibull分布的a，b，c 3个参数；若不满足则令k＝k＋1，同时按照式（13）进行迭代，按照公式

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

2 实例分析

广西沿海涠洲岛海洋监测站3个方向（NNE、E、SSW）1960－1993年的年最大波高见图2～4［14］。

用三参数Weibull分布对图2～4的数据进行拟合，采用最小二乘法、矩估计法、最速下降法及粒子群算法分别对Weibull分布的未知参数进行了估计。

图2 涠洲岛站NNE方向年最大波高曲线Fig.2 Curve of annual maximum wave heights in NNE at Weizhoudao Station

对于矩法，先根据各样本数据求出样本的均值，方差和偏度系数。分别利用式（9），（7）和（6），将其中的总体统计量用样本统计量代换，得到Weibull分布3个参数的解。

在应用粒子群算法时，c1取0.9，c2取0.1，w取0.6，令最大迭代次数Nmax＝100，并且要求前后2次全局最优代入Weibull分布后概率结果的离差平方和小于η＝10－7，取群体中粒子个数popsize＝1 000。由于分布中含有3个未知参数，所以群体维数dimsize＝3。

通过以上方法，拟合Weibull分布中的a，b，c三参数，最终3个方向（NNE，E，SSW）的结果如表1～3。其中表中最速下降法1和最速下降法2是指对于初值不同的最速下降方法的2次计算结果，粒子群算法1和粒子群算法2是指对于没有改变任何参数情况下重复试验后的不同的2次结果。

图5～7是采用最小二乘法、矩估计法、最速下降法及粒子群算法对涠洲岛3个方向年极值波高进行Weibull分布拟合的结果。由图可见，实测数据分布在理论曲线两侧，拟合结果较好。

比较表1～3及图5～7，可以看出：

（1）采用的4种拟合方法中，粒子群算法的概率离差平方和最小。针对本算例，粒子群算法的精度高于最小二乘法，矩估计法，以及最速下降法。

（2）矩估计法和最小二乘法在拟合分布的参数时无需迭代初值，粒子群算法对初值的要求不高，而最速下降法对初值有较高的要求。表1～3中的数据说明：即使初值的选择只有微小的差别，最速下降法得到的结果仍存在很大差异。采用最速下降法实际计算时，若初值设定得不好，计算容易陷入局部最优。

（3）由于粒子群算法在生成初始矩阵时具有随机性，因此，拟合结果不稳定，然而，每次结果相差不大，且与其他3种算法相比，发生概率的离差平方和最小，对观测数据的拟合效果更优。

（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分布的参数拟合中，提高了分布参数拟合的速度。基于广西沿海涠洲岛海洋监测站3个方向（NNE、E、SSW）1960—1993年的年最大波高数据，采用粒子群算法，最小二乘法，矩估计法，以及最速下降法4种方法，拟合了Weibull分布中的参数。计算结果表明：粒子群算法具有对初值要求不高，计算速度快，计算耗时不高等优点。综上可知，粒子群算法是海洋工程环境要素统计中估计Weibull分布参数的1种有效方法。

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

［2］徐以锋，翁永刚，王明星，等.威布尔分布三参数的拟合特性［J］.郑州大学学报：理学版，2003，35（3）：30－34.

［3］汤银才，侯道燕.三参数Weibull分布参数的Bayes估计［J］.系统科学与数学，2009，29（1）：109－115.

［4］严晓东，郑荣越，吴亮，等.三参数威布尔分布参数估计方法比较［J］.宁波大学学报：理工版，2005，18（3）：301－305.

［5］胡恩平，罗兴柏，刘国庆.三参数Weibull分布几种常用的参数估计方法［J］.沈阳工业学院学报，2000，19（3）：88－94.

［6］郑红星，李丽娟.水质模型参数的非数值随机优化［J］.地理研究，2001，20（1）：97－102.

［7］江燕，胡铁松，武夏宁，等.粒子群算法在新安江模型参数优选中的应用［J］.武汉大学学报：工学版，2006，39（4）：14－17.

［8］郭建青，王洪胜，李彦，等.粒子群优化算法在确定河流水质参数中的应用［J］.水利水电科技进展，2007，27（6）：1－5.

［9］路志强，杨文海，康迎宾，等.粒子群算法及其在暴雨强度公式参数优化中的应用［J］.南水北调与水利科技，2008，6（3）：72－73.

［10］郭建青，王洪胜，李彦，等.粒子群优化算法在确定河流水质参数中的应用［J］.水利水电科技进展，2007，27（6）：1－5.

［11］罗航，王厚军，黄建国，等.基于PSO的三参数威布尔分布参数的联合估计方法［J］.仪器仪表学报，2009，30（1）：160－161.

［12］甘应爱，田丰，李维铮，等.运筹学［M］.北京：清华大学出版社，2005：151－155.