廣義帕雷托分布在波浪極值推算中的應用＊

韓樹宗，劉 昆

（中國海洋大學海洋環境學院，山東 青島266100）

海洋中的波浪是海水運動形式之一，它的產生是外力、重力與海水表面張力共同作用的結果。引起海水波動的外力因素有很多，比如風、大氣壓力的變化、天體的引潮力、海底地震以及人為引起的船體的運動等。

海岸和近海工程建筑物處于嚴酷的海洋環境下，在著手進行海上建筑物的規劃和設計之前，必須得到相應海區的可靠波浪資料，掌握其統計特性。為了掌握工程海域的波浪狀況，最可靠的方法是進行現場觀測。1960年代以來，為了適應建設的需要，國家海洋局在沿海各地陸續建立了一些水文氣象觀測站、臺，進行系統的觀測以累積資料。但是對于像石油平臺之類的遠離岸邊的工程建筑，海浪資料的獲取往往非常困難。因此，運用衛星資料進行海浪的數值模擬一直是海浪研究的重要方面，也是海洋預報和分析的主要手段和工具。

海浪的數值模擬發展到20世紀末已達到比較成熟的階段，由荷蘭科學家開發的SWAN（Simulating Waves Nearshore）模式具有較好的海浪模擬精度。近些年國外較為先進的SWAN海浪模式已被廣泛應用于中國海域的風浪模擬。

帕雷托分布是以意大利經濟學家維弗雷多·帕雷托命名的，起初應用于經濟領域，后來漸漸發現其在氣象和水文極值參數分析中也有很好的適用性。廣義帕雷托分布（GPD）是由Pickands在1975年首次引入的，后來很多學者做了進一步研究，該分布被廣泛應用于極值分析領域。國內研究中，張香云等[1]利用GPD分析了洛杉磯降雨量數據，發現GPD能較好的擬合數據分布的尾部，且隨著樣本容量增加，估計效果越來越好；江志紅等[2]利用GPD擬合了我國東部地區的極端降雨量，并且與廣義極值分布（GEV）相比較，GPD擬合效果更好，且計算更為方便；程炳巖等[3]引進GPD模擬重慶極端降水事件，借助于L－矩估計方法，GPD基本不受原始序列樣本量的影響，具有更好精度的實用性和穩定性。國外研究者們則把GPD的應用范圍擴展到海洋水文極值參數的研究中，Bermudez等[4]對GPD的參數估計做了細致的研究，運用了最大似然估計（ML）、概率權重矩法（PWM）和矩法（MOM）分別進行了估計分析；Kevin等[5]利用GPD對北海北部海域百年一遇有效波高進行研究，分析了不同臨界值下百年一遇極值的變化；Philip等[6]運用GPD估計了GOMOS海浪數據百年一遇有效波高，對GPD參數估計引入了波向，討論了考慮波向的分布模型和不考慮波向的模型之間的穩定性問題，認為考慮方向的模型穩定性更好。

對于波浪重現期極值參數的推算工作，大家所熟知的是年最大值統計法（AM），即每年取1個最大值組成樣本的方法，年最大值統計法是廣義極值分布（GEV）的1個重要前提。部分歷時序列統計法（PDS）則是取合適的臨界值，每年可取多個超過臨界值的值組成樣本。GPD屬于運用部分歷時序列統計法的1種分布，對于1年多次取樣法，目前在波浪重現期極值參數的推算當中的研究較少。本文在前人基礎上繼續深入研究，CCMP風場資料（1988年1月～2010年12月）共23年的風場資料對中國海的海浪場進行模擬，并借助GPD模型對重現期極值參數進行推算。

1 資料簡介

1.1 風場資料

CCMP（Cross－Calibrated Multi－Platform）數據是結合了SSM/I、 TMI、AMSR－E、 QuikSCAT 和ADEOS－II幾種資料通過交叉校驗和同化得到的，其時間分辨率為6h，空間分辨率為0.25（°）×0.25（°），時間范圍從1987年7月～2011年6月，空間范圍為0.125°E～359.875°E，78.375°S～78.375°N。CCMP風場資料有著良好的精度和時空分辨率。

本文選取了CCMP風場資料（1988年1月～2010年12月）共23年的風場資料作為SWAN模式的驅動場，進行海浪場的數值模擬。

1.2 實測海浪資料

實測海浪資料采用日本氣象廳（JMA）22001號海洋觀測浮標資料，22001號浮標坐標為126°20′E，28°10′N。浮標提供了有效波高（1990—2000年）和海表面10m高風速（1980—2000年）資料。觀測間隔為3h，當風速大于18m/s時觀測時間間隔變為1h。

圖1 模擬計算范圍網格圖Fig.1 Grid map of calculation range

2 模擬方法和資料驗證

2.1 模擬方法

SWAN模式屬于第三代淺海海浪數值模式，由荷蘭Delft大學土木工程系開發并維護。從第一個公開發布的版本SWAN30.51開始，不斷進行改進和擴充，性能不斷提高，功能也逐漸增強。本文利用SWAN模式對中國海1988年1月～2010年12月的海浪場進行模擬。

模擬海域計算范圍：4°N～42°N，105°E～136°E。驅動風場數據插值到5（′）×5（′）的網格上，模型計算網格分辨率為4（′）×4（′），計算步長為1h，每24h輸出一次結果，計算時間為1988年1月1日00：00時～2010年12月31日00：00時。

2.2 資料驗證

鄭崇偉[7]利用月 CCMP（Cross－Calibrated Multi－Platform）風場資料，對中國海近22年的海表風場特征進行分析，結果顯示，中國海大部分海域的海表風速呈顯著性逐年線性遞增。李燕[8]利用SWAN模式對黃渤海域波浪場進行計算，經系數訂正后，預報準確率達70%以上，有一定的預報能力。梅嬋娟等[9]利用第三代海浪數值模式 WAVEWATCH和SWAN模式，分別對黃海區域進行了理想模擬計算和實際浪場的模擬計算，發現SWAN模式模擬結果較 WAVEWATCH模式好，只是在高風速的模擬情況下，SWAN模式模擬結果偏大，而WAVEWATCH模式模擬結果偏小。

圖2 模擬計算范圍地形圖Fig.2 Topographic map of calculation range

為了更加直觀的比較SWAN模式模擬的結果與浮標實測結果之間的差異，本文在計算中提取了2000年的驅動風場風速值和有效波高模擬值，與22001號浮標（126°20′E，28°10′N）實測資料相對應，作出散布圖。并且進行了誤差分析，計算了偏差（Bias）、平均絕對誤差（MAE）、均方根誤差（RMSE）和相關系數（CC）。

式中：xi代表浮標觀測數據；yi代表SWAN計算結果；分別代表觀測數據和計算結果的平均值；n為樣本總量。

圖3 實測風速與CCMP混合風場風速散布圖Fig.3 Comparison of wind speed between observed data and CCMP data

圖4 實測有效波高與CCMP混合風場有效波高散布圖Fig.4 Comparison of significant wave height between observed data and CCMP data

由圖3可見，驅動風場和22001號浮標觀測的風速相一致，統計上存在0.54的正偏差，說明驅動風場風速稍大于實測風速，相關系數為0.77，通過了99%的信度檢驗，均方根誤差為2.55m/s，平均絕對誤差為1.95m/s，衛星風場精度較好。由圖4可見，模擬有效波高和22001號浮標觀測的有效波高相一致，統計上存在0.01的負偏差，說明模擬有效波高稍大于實測有效波高，相關系數為0.75，通過了99%的信度檢驗，均方根誤差為0.71m，平均絕對誤差為0.48m，模擬的海浪場可用。

圖5 計算有效波高與實測有效波高極大值驗證圖Fig.5 Comparison of extreme significant wave height between observed data and calculated data

本文所關心的是波浪的極值參數，因此波浪極值參數模擬的準確度就至關重要。鑒于22001號浮標有較為完整的1990—2000年波高資料，根據浮標資料分別補充驗證了年最大值，第二大值和第三大值（見圖5）。圖5中點表示浮標實測有效波高，線表示模擬計算有效波高。其中，1988年計算波高普遍較小，低于4m，最大值、第二大值和第三大值均取為4m。從圖5中可以看出，模擬結果極值有效波高在大部分年份偏低，這與SWAN模式在極端風速下運算不穩定有關，但整體差異不算很大。本文采用此次SWAN模擬的波浪場數據作為樣本，進行波浪極值參數推算新方法的嘗試，還是可行的。

3 極值估計方法

目前國內外常用的極值估計方法總的來說可以分為兩類：一類是經驗型的，如皮爾遜III型曲線法等，一類是以極值分部理論為基礎的，如耿貝爾（Gumbel）曲線法、威布爾（Weibull）曲線法等。皮爾遜III型曲線法雖然有較大的實用性，但是其缺乏嚴格的概率論理論依據，在海洋資料的分析中，最常用的是耿貝爾分布和威布爾分布。耿貝爾分布和威布爾分布都屬于廣義極值分布（GEV）的特殊形式。

GEV的分布形式函數為：

式中：γ成為形狀參數或尾部指數；σ成為尺度參數；u為閾值。當γ＝0時，GEV簡化為Tippett I型分布，即Gumbel分布；當γ＞0時，為Tippett II型分布；放γ＜0時，則為Tippett III型分布，即 Weibull分布。

前面提到廣義極值分布（GEV）的1個重要前提就是年最大值統計法。對于一些年份，可能出現年最大值和年次大值差距很大的情況，也可能出現1年有多個差不多的極大值的情況，如果采取每年抽取1個最大值的方法，其實并不符合實際。因此，這種年最大值統計法，會舍去很多有用的信息，形成的樣本量小，不能充分利用分析資料。而廣義帕雷托分布（GPD）屬于運用部分歷時序列統計法的1種分布，每年可取多個超過臨界值的值組成樣本，一定程度上解決了上述方法的缺陷。

GPD的分布形式函數為：

式中：γ成為形狀參數或尾部指數；σ成為尺度參數；u為閾值。對于給定的重現期T（年），可證明[3]重現期極值xT如下：

式中：λ成為年交叉率[10]，即極值超過給定閾值的個數。對于給定閾值情況下的年交叉率，一般采用多年平均的年交叉率即可。

式中：n為超過給定閾值的總樣本量；A為資料的總年數。在以上的理論基礎上，本文采用L－矩參數估計方法求得GPD對應的參數，進而對重現期極值進行推算。

3.1 L－矩法參數估計

參數估計方法中比較常用的有矩法（MOM）、概率權重矩法（PWM）和最大似然估計（ML）。根據段忠東等[11]對不同極值概率分布參數的研究，雖然最大似然估計具有較高的精確度和穩定性，但在多數情況下，矩法和最大似然估計法都不能獲得參數估計的解析表達，需要數值求解強非線性方程組。L－矩法起源于“概率權重矩（PWM）”，是概率權重矩的線性組合，他的一個顯著優點是可以得到參數的顯示表達式，從而使得參數估計變的簡便，L－矩法的參數估計精度較高，估計值的穩定性也比較強。

首先將樣本由大到小排序（x1≤x2≤ …xn），根據PWM的定義[12]，隨即變量x的第r階概率加權矩為：

線性組合后的樣本L－矩前四階計算式如下：

根據文獻[13－14]，GPD的L－矩估計式為：

根據樣本L－矩和GPD的L－矩估計式，可以推算出GPD的參數估計式如下：

4 GPD方法應用實驗

4.1 不同閾值下的GPD模型參數估計及其效果檢驗

在以上的理論參考下，對不同閾值情況下GPD模型的擬合結果進行檢驗分析。鑒于前面22001號浮標計算結果驗證良好，故選取22001號浮標（126°20′E，28°10′N）點作為實驗點，對該點處SWAN計算得出的23年的有效波高數據進行分析，分析結果見表1。選取科爾莫格洛夫檢驗統計量（K－S）、相關系數和均方根誤差3個指標對擬合見過進行檢驗。

由表1可見，GPD模型的模擬結果較好，K－S統計量很小，擬合的分布函數和經驗分布函數之間的相關系數都在0.9以上，均方根誤差幾乎為零。從表1還可以看出：隨著閾值的增大，樣本量不斷下降，一直到閾值為5.5m時，樣本量為25，仍大于年最大值法的樣本量，利用有限的數據獲取更豐富的信息；隨著閾值的增大，擬合的分布函數和經驗分布函數之間的相關系數不斷減小，這與樣本量的減小是息息相關的，這說明樣本量越大，GPD模型的擬合結果越好。

表1 不同閾值下的GPD模型參數估計及其效果檢驗Table 1 The test of effect for GPD model under different thresholds

4.2 選定閾值（4m）的情況下，對GPD模型和GEV模型的模擬結果進行比較

選定閾值為4m，GPD模型運用部分歷時序列統計法提取樣本。GEV模型采用年最大值統計法提取樣本。分別用GPD分布和GEV分布對擬合累積頻率曲線并進行檢驗，結果如下：

表2 兩種模型擬合結果檢驗Table 2 Comparison of GPD and GEV models

圖6 實驗點累積概率的GPD模擬曲線Fig.6 The simulated curve of GPD for cumulative probability

從以上結果可以看出，GPD模型的模擬結果在各方面都好于GEV模型，檢驗指標上也表現良好，主要原因在于GPD模型更充分的利用了有限的資料，取得了更豐富的樣本，曲線擬合效果更好，在海浪極值參數的推算上有一定應用價值。4.3選定閾值（4m）的情況下，對浮標附近區域尺度參數σ和百年一遇有效波高的求解分析

圖7 實驗點累積概率的GEV模擬曲線Fig.7 The simulated curve of GEV for cumulative probability

尺度參數σ是標準差線性函數，σ的大小標識著極值有效波高的穩定性，σ越大表明極值有效波高之間的差別也越大越不穩定。下面選定浮標附近區域（27°N～29°N，126°E～127°E）為例，對區域內σ的分布進行簡單的分析。由圖可以看出，在選定區域內浮標西部為尺度參數高值區，浮標東部為尺度參數低值區。說明浮標西部區域極大波高值較不穩定，浮標東部區域相反，極大波高值較穩定。

從浮標附近區域百年一遇有效波高空間分布來看，浮標西南部百年一遇有效波高較大，最大值超過15m，浮標東北部百年一遇有效波高較小，多在10m以下。大浪的成長需要足夠的風區，浮標所在的海域近似可以看成由中國大陸、朝鮮半島、日本群島、琉球群島和臺灣島所包圍。從有效波高的分布來看，岸界附近的海域有效波高較小，遠離岸界的海域有效波高較大，同時有利于極值波高的出現。因此圖中浮標左側波高大于右側。

圖8 22001號浮標附近區域σ分布Fig.8 Spatial distribution ofσnear NO.22001buoy

5 結論

圖9 22001號浮標附近區域百年一遇有效波高分布（m）Fig.9 Spatial distribution of significant wave height of 100－year return level near NO.22001buoy

隨著經濟的發展，海洋石油的開采越來越受到國際的關注，海洋石油平臺的選址不僅要考慮油氣資源分布、水深、地質、氣象等條件，工程海域的海洋環境狀況不容忽視。浮標東北部海域百年一遇有效波高較小且極大波高值較穩定，海洋環境狀況明顯優于浮標西部海域。在實際應用中，GPD分布模型有著較強的應用價值。

本文利用CCMP衛星風場資料模擬了中國海域1988—2010年共23年的波浪場，并利用日本氣象廳22001號海洋觀測浮標資料進行了風場和有效波高驗證；采用廣義帕雷托分布（GPD）模型對22001號浮標點計算有效波高極值進行擬合，分析了GPD模型在不同閾值下的估計結果并進行效果檢驗，比較了GPD模型和GEV模型的優劣；計算了浮標附近海域尺度參數σ和百年一遇有效波高分布，得出如下結論：

（1）在不同閾值的情況下，樣本量越大，GPD模型的擬合結果越好。

（2）GEV每年只取一個極大值，得到的樣本量小，資料信息較少，GPD模型采用部分歷時序列統計法采樣，增大了樣本容量，擬合結果在相關程度和效果檢驗方面都表現出色。從前人的研究成果和本文的驗證結果來看，GPD模型不僅在降雨、洪水等方面有較好的應用，在海洋水文極值參數的估計中也有很好的應用價值。

（3）GPD模型的尺度參數σ的大小可以反映極值的穩定性，尺度參數σ和百年一遇有效波高可以在一定程度上表征出海洋環境惡劣狀況，它們的空間分布特征對海洋石油平臺工程的選址有一定參考價值。

[1]張香云，趙旭.廣義Pareto模型統計推斷及其應用 [J].數理統計與管理，2011，30（6）：989－995.

[2]江志紅，丁裕國，朱蓮芳，等.利用廣義帕雷托分布擬合中國東部日極端降水的試驗 [J].高原氣象，2006，28（3）：573－580.

[3]程炳巖，丁裕國，張金鈴，等.廣義帕雷托分布在重慶暴雨強降水研究中的應用 [J].高原氣象，2008，27（5）：1004－1009.

[4]Pierre Ribereau，Philippe Naveau，Armelle Guillou.A note of caution when interpreting parameters of the distribution of excesses[J].Advances in Water Resources，2011，34：1215－1221.

[5]Kevin Ewans，Philip Jonathan.The effect of directionality on Northern North Sea extreme wave design criteria [J].Journal of Offshore Mechanics and Arctic Engineering，2008，130：1－8.

[6]Philip Jonathan，Kevin Ewans.The effect of directionality on extreme wave design criteria [J].Ocean Engineering，2007，34：1977－1994.

[7]鄭崇偉.基于CCMP風場的近22年中國海海表風場特征分析[J].氣象與減災研究，2011，34（3）：41－46.

[8]李燕，薄兆海.SWAN模式對黃渤海海域浪高的模擬能力試驗[J].海洋預報，2005，22（3）：75－82.

[9]梅嬋娟，趙棟梁，史劍.兩種海浪模式對中國黃海海域浪高模擬能力的比較 [J].海洋預報，2008，25（2）：92－98.

[10]郭軍，任國玉，李明財.環渤海地區極端降水事件概率分布特征[J].氣候與環境研究，2010，15（4）：425－432.

[11]段忠東，周道成.極值概率分布參數估計方法的比較研究 [J].哈爾濱工業大學學報，2004，36（12）：1605－1609.

[12]丁裕國，劉吉峰，張耀存.基于概率加權估計的中國極端氣溫時空分布模擬試驗 [J].大氣科學，2004，28（5）：771－782.

[13]Hosking J R M，Wallis J R.Parameter and quantile estimation for the generalized Pareto distribution[J].Technometries，1987，29：339－349.

[14]Zhai Panmao，Sun Anjing，Ren Fumin，et al.Changes of climate extremes in China[J].Climatic Change，1999，42：203－218.

[15]P deZeaBermudeza，SamuelKotzb.Parameter estimationofthegeneralizedParetodistribution[J].Journal of Statistical Planning and Inference，2010，140：1353－1373.