韓樹宗，劉 昆

（中國海洋大學(xué)海洋環(huán)境學(xué)院，山東 青島266100）

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

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

海浪的數(shù)值模擬發(fā)展到20世紀(jì)末已達(dá)到比較成熟的階段，由荷蘭科學(xué)家開發(fā)的SWAN（Simulating Waves Nearshore）模式具有較好的海浪模擬精度。近些年國外較為先進(jìn)的SWAN海浪模式已被廣泛應(yīng)用于中國海域的風(fēng)浪模擬。

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

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

1 資料簡介

1.1 風(fēng)場資料

CCMP（Cross－Calibrated Multi－Platform）數(shù)據(jù)是結(jié)合了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風(fēng)場資料有著良好的精度和時空分辨率。

本文選取了CCMP風(fēng)場資料（1988年1月～2010年12月）共23年的風(fēng)場資料作為SWAN模式的驅(qū)動場，進(jìn)行海浪場的數(shù)值模擬。

1.2 實測海浪資料

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

圖1 模擬計算范圍網(wǎng)格圖Fig.1 Grid map of calculation range

2 模擬方法和資料驗證

2.1 模擬方法

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

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

2.2 資料驗證

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

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

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

式中：xi代表浮標(biāo)觀測數(shù)據(jù)；yi代表SWAN計算結(jié)果；分別代表觀測數(shù)據(jù)和計算結(jié)果的平均值；n為樣本總量。

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

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

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

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

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

3 極值估計方法

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

GEV的分布形式函數(shù)為：

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

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

GPD的分布形式函數(shù)為：

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

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

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

3.1 L－矩法參數(shù)估計

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

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

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

根據(jù)文獻(xiàn)[13－14]，GPD的L－矩估計式為：

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

4 GPD方法應(yīng)用實驗

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

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

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

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

4.2 選定閾值（4m）的情況下，對GPD模型和GEV模型的模擬結(jié)果進(jìn)行比較

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

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

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

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

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

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

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

圖8 22001號浮標(biāo)附近區(qū)域σ分布Fig.8 Spatial distribution ofσnear NO.22001buoy

5 結(jié)論

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

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

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

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

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

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

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

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

[3]程炳巖，丁裕國，張金鈴，等.廣義帕雷托分布在重慶暴雨強(qiáng)降水研究中的應(yīng)用 [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風(fēng)場的近22年中國海海表風(fēng)場特征分析[J].氣象與減災(zāi)研究，2011，34（3）：41－46.

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

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

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

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

[12]丁裕國，劉吉峰，張耀存.基于概率加權(quán)估計的中國極端氣溫時空分布模擬試驗 [J].大氣科學(xué)，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.