王志旭，陳子燊

（中山大學(xué) 水資源與環(huán)境系，廣東 廣州 510275）

在海岸工程建設(shè)設(shè)計(jì)中，最重要的又難以決定的問題之一是設(shè)計(jì)波浪標(biāo)準(zhǔn)的選擇，選用的設(shè)計(jì)波高有少量變化就會(huì)明顯影響到工程的費(fèi)用、涉及工程安全和維修方案。在得到波高長期樣本數(shù)據(jù)后，可通過概率分布函數(shù)來擬合實(shí)測數(shù)據(jù)，并推算設(shè)計(jì)波高。Goda等 （1990）、Ferreira等 （2000）和Todd（2000）對各種分布函數(shù)做過分析比較，發(fā)現(xiàn)各種函數(shù)有其適用性。目前常用的分布函數(shù)有對數(shù)正態(tài)分布、Gumbel分布、皮爾遜Ⅲ型分布、Weibull分布。我國海港水文規(guī)范規(guī)定，對于年極值波高及其對應(yīng)的周期的理論頻率曲線，一般采用皮爾遜Ⅲ型曲線，也可以實(shí)測經(jīng)驗(yàn)累積頻率點(diǎn)擬合最佳為原則，選用其他理論頻率曲線。

1927年，F(xiàn)réchet首先給出威布爾分布的定義，即極值Ⅲ型分布。Rosin等（1933）在研究碎末的分布時(shí)，第一次應(yīng)用了威布爾分布。1951年，瑞典工程師、數(shù)學(xué)家Weibull通過研究軸承壽命、結(jié)構(gòu)強(qiáng)度和疲勞等問題從而詳細(xì)解釋了這一分布，于是，該分布便以他的名字命名為Weibull distribution。隨后，威布爾分布廣泛應(yīng)用于質(zhì)量管理上的元件和材料壽命可靠性研究，是可靠性分析和壽命檢驗(yàn)的理論基礎(chǔ)。1954年Gumbel將其應(yīng)用到年最枯流量的水文統(tǒng)計(jì)分析中，Neelamani等（2007）和Seguro等（2000）分別將其應(yīng)用于風(fēng)速和波高的極值推算。

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

威布爾曲線的概率密度函數(shù)為：

分布函數(shù)為：

式中，a≥0，為位置參數(shù)；b＞0，為尺度參數(shù)；c＞0，為形狀參數(shù)；x≥a。

三參數(shù)威布爾分布型式具有較大靈活性，形狀參數(shù)的存在使其在數(shù)據(jù)擬合上極富彈性，比二參數(shù)威布爾分布有更高的擬合精度，然而其參數(shù)的優(yōu)良估計(jì)比較復(fù)雜。對于威布爾分布的參數(shù)估計(jì)，國內(nèi)外有眾多學(xué)者進(jìn)行相關(guān)研究，傳統(tǒng)的參數(shù)估計(jì)方法有極大似然估計(jì)法、最小二乘法、矩估計(jì)法等。我國學(xué)者相繼提出了雙線形回歸法（莊渭峰，1999）、概率權(quán)重矩法（張秀芝，1994）、相關(guān)系數(shù)優(yōu)化法（傅惠民等，1990；史景釗等，1995）和灰色估計(jì)法（鄭榮躍等，1989）。不過，各方法存在著相應(yīng)的優(yōu)缺點(diǎn)，各自的使用范圍和精度也都不同。本文利用蒙特卡羅法對常用的4種參數(shù)方法（極大似然法、概率權(quán)重矩法、相關(guān)系數(shù)優(yōu)化法和灰色估計(jì)法）進(jìn)行模擬，討論各方法的適用性及其精度，并用上述各種參數(shù)估計(jì)方法對潿洲站34年實(shí)測年極值波高進(jìn)行設(shè)計(jì)波高的推算，分析其優(yōu)良性。

1.1 極大似然估計(jì)法

極大似然估計(jì)（Maximum Likelihood Estimation，MLE）的基本思想是選擇待定參數(shù)使樣本出現(xiàn)在觀測值的領(lǐng)域內(nèi)的概率最大，并以這個(gè)值作為未知參數(shù)的點(diǎn)估計(jì)值（朱勇華等，1999）。

對于此超越方程組有多種算法，本文采用王華勝等（2004）提出的求解方法。

1.2 相關(guān)系數(shù)優(yōu)化法

相關(guān)系數(shù)優(yōu)化法（Correlation Coefficient Method，CCM）是我國學(xué)者史景釗等（1995）提出的一種估計(jì)威布爾分布位置參數(shù)的新方法，對式（2）進(jìn)行變換得：

上式可以簡化為線性方程

將獲取的極值波高序列從小到大排列得到序列xi，通過極值波高累計(jì)率估計(jì)量F（xi）計(jì)算得到序列Y，由式（4）、（5）可知，如果位置參數(shù)a估計(jì)正確，X、Y將呈線性關(guān)系，即X、Y之間有著最大的相關(guān)系數(shù)。

根據(jù)這個(gè)思路，可以給定a一個(gè)初始值（略小于波高序列最小值）和一個(gè)步長（步長可根據(jù)實(shí)際數(shù)據(jù)進(jìn)行設(shè)定），a按步長遞增可得一系列的X與Y的相關(guān)系數(shù)R2。a與R2的關(guān)系為：隨著a的增加，R2先增后減，當(dāng)R2達(dá)到最大值時(shí)，此時(shí)的a為最佳理論估計(jì)值。

另外可以通過推導(dǎo)證明，在下式成立時(shí)，相關(guān)系數(shù)R2最大：

用相關(guān)系數(shù)法求出位置參數(shù)a后，c和B的參數(shù)估計(jì)可用最小二乘法或其他方法進(jìn)行估計(jì)。

1.3 灰色估計(jì)法

灰色估計(jì)法是鄭榮躍等（1989）基于鄧聚龍?zhí)岢龅幕疑到y(tǒng)理論發(fā)展的一種較新的參數(shù)估計(jì)方法。在已知數(shù)列（如波高、風(fēng)速序列）前提下，以單數(shù)列的微分方程GM（1，1）為基礎(chǔ)，得到各類灰色預(yù)測方法，通過預(yù)測方法進(jìn)行參數(shù)估計(jì)。

式 （1-2）可變換為：

式 （8）對應(yīng)于灰色估計(jì)GM（1，1）的時(shí)間響應(yīng)模型：

從而可用GM（1，1）模型對參數(shù)a，b，e，d進(jìn)行估計(jì)，再由此得到威布爾分布三參數(shù)的估計(jì)值（其中，c=-1/d）?；疑烙?jì)的參數(shù)為：

式中：

1.4 概率權(quán)重矩法

Greenwood等（1979）提出了概率權(quán)重矩法（Probability Weighted Moments，PWM），并應(yīng)用于Gumbel、GEV和Wakeby分布的參數(shù)估計(jì)。張秀芝（1994）將其運(yùn)用于三參數(shù)Weibull分布的參數(shù)估計(jì)，并推算了極值風(fēng)速。樣本的概率權(quán)重矩估計(jì)值為：

其中，Pi為對應(yīng)于xi的累積頻率。用概率權(quán)重矩表示的威布爾分布參數(shù)估計(jì)值為：

2 參數(shù)估計(jì)方法分析比較

利用蒙特卡羅模擬方法對以上4種參數(shù)估計(jì)方法進(jìn)行參數(shù)估計(jì)并分析其優(yōu)劣性、適用性。為使計(jì)算估計(jì)結(jié)果更加精確，本文中蒙特卡羅模擬運(yùn)行次數(shù)J=5 000，樣本數(shù)分別取20、50、100、200、500，選取參數(shù)估計(jì)量的偏差、標(biāo)準(zhǔn)差和均方誤差來評價(jià)估計(jì)結(jié)果。估計(jì)參數(shù)的偏差、標(biāo)準(zhǔn)差和均方誤差越小，估計(jì)方法越好。威布爾分布三參數(shù)的估計(jì)值為，令待估參數(shù) a=20，b=30，c=2.5，模擬結(jié)果見表1：

（1）通過偏差、標(biāo)準(zhǔn)差和均方誤差的對比，對于不同樣本大小，極大似然估計(jì)法參數(shù)估計(jì)結(jié)果較好，相關(guān)系數(shù)法次之，灰色估計(jì)法的效果最差。當(dāng)樣本數(shù)越來越大時(shí)，各種參數(shù)估計(jì)方法的結(jié)果差異逐漸縮小。

（2）灰色估計(jì)法在小樣本時(shí)即可建模；概率權(quán)矩法和灰色估計(jì)法計(jì)算簡單，避免了試定和迭代等復(fù)雜計(jì)算，然而在小樣本時(shí)難以得到良好參數(shù)估計(jì)；相關(guān)系數(shù)法假定位置參數(shù)是獨(dú)立的，雖能得到最佳的位置參數(shù)估計(jì)，但忽略了位置參數(shù)與形狀參數(shù)、尺度參數(shù)的相關(guān)性，盡管有此缺陷，其估計(jì)所得的參數(shù)仍有較高精度；極大似然估計(jì)法算法最為復(fù)雜，需要進(jìn)行試定和迭代計(jì)算。

（3）概率權(quán)矩法、灰色估計(jì)法、相關(guān)系數(shù)法均受樣本元素經(jīng)驗(yàn)概率Pi影響。Pi取值不同，估算精度也有所差別。對于概率權(quán)重矩法，Hosking（1985）認(rèn)為可采用估算值：Pi=（i-0.35）/n。表1各種參數(shù)估計(jì)方法均采用此公式。

為分析不同參數(shù)估計(jì)方法適用的經(jīng)驗(yàn)頻率，采用以下多種經(jīng)驗(yàn)頻率公式（黃振平，等，2002）進(jìn)行蒙特卡羅模擬：海森公式：Pi=（i-0.5）/n；數(shù)學(xué)期望公式：Pi=i/（n+1）；中值公式：Pi=（i-0.3）/（n+0.4）；Blom公式：Pi=（i-0.375）/（n+0.25）；Tukey公式：Pi=（3i-1）/（3n+1）；Gringorten公式：Pi=（i-0.44）/（n+0.12）。結(jié)果表明：灰色估計(jì)法采用海森公式時(shí)，計(jì)算結(jié)果精度較好，也相對穩(wěn)定；概率權(quán)重矩法仍采用上述公式；采用海森公式和Gringorten公式進(jìn)行相關(guān)系數(shù)優(yōu)化法參數(shù)估計(jì)所得的結(jié)果相差很小，但均優(yōu)于其他經(jīng)驗(yàn)頻率公式所得結(jié)果。

表1 蒙特卡羅模擬不同樣本長度與不同參數(shù)估計(jì)的誤差

3 設(shè)計(jì)波高推算

以潿洲站東（E）向浪34a的年最大波高為樣本（夏華永等，2001）其中最大波高為5.3m，最小波高為0.9m，平均波高為2.85m，利用上述各種三參數(shù)威布爾分布的參數(shù)估計(jì)方法進(jìn)行設(shè)計(jì)波高的推算，其中相關(guān)系數(shù)優(yōu)化法和灰色估計(jì)法的經(jīng)驗(yàn)頻率公式使用海森公式，概率權(quán)重矩法采用Pi=（i-0.35）/n。為使實(shí)測波高與Weibull分布擬合有著更直觀的線性對比，對式（2）進(jìn)行變換，得Weibull簡化變量W：

各種參數(shù)估計(jì)方法擬合的結(jié)果見圖1。利用相關(guān)系數(shù)R、均方根誤差和Q統(tǒng)計(jì)值來分析比較各種參數(shù)估計(jì)結(jié)果的優(yōu)良性。1）均方根誤差：

2）Q統(tǒng)計(jì)值：

式中：xi為實(shí)測樣本序列排序后觀測值和經(jīng)驗(yàn)累積頻率分布Pei計(jì)算的與威布爾分布相應(yīng)的分位數(shù)；Pi：理論累積頻率；n：樣本容量。在樣本容量相同的情況下均方根誤差和Q值越小，擬合越好。

由表2可知，極大似然估計(jì)法的相關(guān)系數(shù)最大且Q統(tǒng)計(jì)值最小，相關(guān)系數(shù)法的均方根誤差最小。各估計(jì)方法均估算得良好的參數(shù)，使得波高數(shù)據(jù)的擬合具有最好的線性度（圖1），綜合而言，極大似然估計(jì)法所得結(jié)果更加精確。相關(guān)系數(shù)優(yōu)化法與灰色估計(jì)法所得的三參數(shù)大小均相差較大，推算而得的H1%相差達(dá)0.22m；而概率權(quán)矩法和極大似然估計(jì)法參數(shù)差異較小。采用極大似然估計(jì)法計(jì)算結(jié)果，潿洲島50年一遇極值波高為4.98m，100年一遇極值波高為5.27m。1992年潿洲站年極值波高為5.3m，達(dá)到100年一遇。

圖1 三參數(shù)W eibu ll分布線性擬合結(jié)果

表2 不同參數(shù)估計(jì)方法的設(shè)計(jì)波高推算結(jié)果

4 結(jié)論

（1）對于各種容量的樣本，極大似然估計(jì)法參數(shù)估計(jì)結(jié)果精度較好，然而在小樣本時(shí)（n=20），各種估計(jì)方法所得結(jié)果均不理想，偏離實(shí)際參數(shù)較大。隨著樣本容量的增加，各種估計(jì)方法之間的差異越來越小，當(dāng)樣本容量大于100時(shí)，各方法所得結(jié)果已很接近。

（2）需選用合適的樣本元素經(jīng)驗(yàn)概率Pi，以保證取得較好的參數(shù)估計(jì)，其中相關(guān)系數(shù)法和灰色估計(jì)法宜采用海森公式，概率權(quán)矩法采用Pi=1-（i-0.35）/n。

（3）一般年極值波高序列樣本長度有限（少于50年），進(jìn)行設(shè)計(jì)波高推算時(shí)可采用精度較高的極大似然估計(jì)法，也可根據(jù)實(shí)際情況進(jìn)行估計(jì)方法分析對比，選擇合適的估計(jì)方法。

（4）利用極大似然法估計(jì)威布爾分布三參數(shù)，推算而得潿洲島50年一遇極值波高為4.98m，100年一遇極值波高為5.27m。在34年的年極值波高序列中，1992年的年極值波高達(dá)到了100年一遇。

Ferreira J A，Guedes S C，2000.Modelling distributions of significant waveheight.CoastalEngineering，4（40）：361-374.

Goda Y，KobuneK，1990.Distribution function fitting forstorm wave data.Proceedingsof the22nd InternationalConferenceon Coastal Engineering，ASCE，18-31.

Greenwood J A，Landwehr JM，Matalas N C，et al，1979.Probability Weighted Moments：Definition and Relation to Parametersof Several Distributions Expressible in Inverse Form.Water Resources Research，15（5）：1049-1054.

Hosking JR M，Wallis JR，Wood E F，1985.Estimation of the generalized extreme value distribution by themethod of probability weight-ed moments.Technometrics，27（3）：251-261.

Neelamani S，Al-Salem K，Rakha K，2007.Extreme waves for Kuwaiti territorial waters.Ocean Engineering，34（10）：1496-1504.

Rosin P，Ramm ler E，1933.The laws governing the fineness of powdered coal.Journalof the Instituteof Fuel，7：29-36.

Seguro JV，Lambert TW，2000.Modern estimation of the parameters of the Weibullwind speed distribution for wind energy analysis.Journal of Wind Engineeringand Industrial Aerodynamics，1（85）：75-84.

Todd L，2000.Distributions for storm surge extremes.Ocean Engineering，27（12）：1279-1293.

傅惠民，高鎮(zhèn)同，1990.確定威布爾分布三參數(shù)的相關(guān)系數(shù)優(yōu)化法.航空學(xué)報(bào)，11（7）：323-327.

黃振平，薩迪伊，王春霞，等，2002。關(guān)于適線法中經(jīng)驗(yàn)頻率計(jì)算公式的對比研究.水利水電科技進(jìn)展，22（5）：5-8

史景釗，蔣國良，1995.用相關(guān)系數(shù)法估計(jì)威布爾分布的位置參數(shù).河南農(nóng)業(yè)大學(xué)學(xué)報(bào)，29（2）：167-171.

王華勝，李忠厚，林榮文，等，2004.耗損故障的三參數(shù)Weibull分布極大似然估計(jì)方法.中國鐵道科學(xué)，25（5）：39-42.

夏華永，李樹華，2001.廣西沿海年極值波高分析.熱帶海洋學(xué)報(bào)，20（2）：1-7.

張秀芝，1994.概率權(quán)重矩法及其在Weibull分布參數(shù)估計(jì)中的應(yīng)用.海洋預(yù)報(bào)，11（3）：55-61.

鄭榮躍，秦子增，1989.Weibull分布參數(shù)估計(jì)的灰色方法.強(qiáng)度與環(huán)境，（2）：34-40.

朱勇華，邰淑彩，1999．應(yīng)用數(shù)理統(tǒng)計(jì)．武漢：武漢水利電力大學(xué)出版社.

莊渭峰，1999.用微機(jī)實(shí)現(xiàn)威布爾參數(shù)的雙線性回歸最小二乘估計(jì).電子產(chǎn)品可靠性與環(huán)境試驗(yàn)，5：2-7.