趙嘉靜，馮 曦，馮衛(wèi)兵，李慧超

(河海大學(xué) 港口海岸與近海工程學(xué)院，南京 210098)

海岸工程的規(guī)劃和建設(shè)對(duì)一個(gè)國(guó)家的經(jīng)濟(jì)發(fā)展起到帶頭作用。為了保證海岸基礎(chǔ)設(shè)施的設(shè)計(jì)與建設(shè)的安全穩(wěn)定，如碼頭與港口的設(shè)計(jì)與施工[1-2]，需要對(duì)波浪參數(shù)有確切的掌握。波高極值資料對(duì)于工程的設(shè)計(jì)與施工都不可或缺。

波高的長(zhǎng)期變化趨勢(shì)可以通過累積分布函數(shù)圖和重現(xiàn)期分布來(lái)表征，此方法廣泛應(yīng)用于沿海極端波高的研究[3-5]。極值分析(EVA)在海洋、氣象、金融、交通等領(lǐng)域的應(yīng)用十分廣泛，極值分析理論應(yīng)用于極端波高的分析可追溯于Goda[6]、Mathiesen等人[7]和Menendez等人[8]，由于工程的需要，在中國(guó)近海海域已有基本的研究[9-12]。廣義極值分布函數(shù)(GEV)經(jīng)常被應(yīng)用于極值的統(tǒng)計(jì)分布，最常用的方法是對(duì)于長(zhǎng)期數(shù)據(jù)，采用每年取一個(gè)極大值的方法進(jìn)行擬合，缺點(diǎn)是往往會(huì)忽略每年發(fā)生的其他極端情況。對(duì)此有2種解決辦法：分塊取極值法和超閾值(peak-over-threshold，POT)方法[13-14]。分塊取極值法最常用的是采用取月極值波高的方法，即每個(gè)月取一個(gè)極值，盡可能統(tǒng)計(jì)到發(fā)生在全年的極端波高情況；POT方法即選取一個(gè)特定的閾值，對(duì)超過該閾值的值進(jìn)行統(tǒng)計(jì)并分析。極端波高的長(zhǎng)期變化趨勢(shì)一直以來(lái)備受關(guān)注，Komar等人[15]使用NOAA浮標(biāo)測(cè)量數(shù)據(jù)對(duì)美國(guó)西海岸等地也做了類似的研究。

南黃海地區(qū)位于中國(guó)江蘇、浙江省沿岸，由于輻射沙洲的存在，近岸地形復(fù)雜(圖1-a)，同時(shí)受到季風(fēng)影響，多臺(tái)風(fēng)浪和風(fēng)暴潮，出現(xiàn)極端水位情況較多[16]?，F(xiàn)有的測(cè)站多位于沿岸且數(shù)據(jù)有缺失，不能完全表示出整個(gè)南黃海地區(qū)的長(zhǎng)期波浪分布情況，因此不僅需要合適的極值計(jì)算模型，還需進(jìn)一步分析數(shù)據(jù)長(zhǎng)度對(duì)于重現(xiàn)波高計(jì)算的影響[13]。本文研究目的在于利用3種方法計(jì)算南黃海地區(qū)的百年重現(xiàn)波高空間分布，并進(jìn)行比選，得到不同情況下的最優(yōu)方法。同時(shí)，計(jì)算3種方法受到不同時(shí)間跨度的影響，為工程建設(shè)提供指導(dǎo)。

1 模型及極值計(jì)算方法介紹

1-a 區(qū)域地形 1-b 網(wǎng)格區(qū)域

1.1 模型設(shè)置

本研究采用SWAN模型進(jìn)行模擬，此波浪模型已廣泛應(yīng)用于中國(guó)近海海域，模擬結(jié)果良好[17-19]。網(wǎng)格經(jīng)度范圍119°E～124°E，緯度范圍31°N～35°N(圖1-b)，共有10 429個(gè)節(jié)點(diǎn)和20 097個(gè)單元，精度范圍從0.02°到0.11°(1 322～12 454 m)。利用歐洲天氣預(yù)報(bào)中心(ECMWF)1979～2018年共40 a的風(fēng)場(chǎng)數(shù)據(jù)作為初始條件，時(shí)間分辨率為6 h，空間分辨率為12 km。物理過程考慮了包括白帽、破波、底摩擦、波-波相互作用在內(nèi)的波能耗散機(jī)制。SWAN模型的控制方程如下[20]

(1)

S=Sin+Swc+Sbrk+Sbot+Snl4+Snl3

(2)

式中：N為波作用量密度，方程左邊第1項(xiàng)為動(dòng)密度譜對(duì)時(shí)間的偏導(dǎo)，第2～5項(xiàng)分別為能量在x、y、σ、θ空間的傳播，Ci(i=x、y、σ、θ)為各個(gè)空間的能量傳播速度。方程右邊S為源項(xiàng)，包括能量產(chǎn)生項(xiàng)(風(fēng)能輸入Sin)，耗散項(xiàng)(白帽耗散Swc，深度誘導(dǎo)引起的波浪破碎Sbrk，底摩阻耗散Sbot)及轉(zhuǎn)化過程(四波相互作用Snl4, 三波相互作用Snl3)，σ為相對(duì)頻率。

1.2 模型驗(yàn)證

表1 有效波高的驗(yàn)證參數(shù)

為了證明模擬結(jié)果的可信性，選取大豐(120.81°E，33.285°N)和蠣蚜山(121.568°E，32.147°N)2個(gè)浮標(biāo)測(cè)點(diǎn)2013年的實(shí)測(cè)波高(Hs)數(shù)據(jù)進(jìn)行對(duì)比。相關(guān)系數(shù)R，均方根誤差RMSE及偏差Bias如表1所示，大豐測(cè)站吻合較好，蠣蚜山測(cè)站由于浮標(biāo)放置在狹窄的潮道中，潮流的影響較明顯，故在單純考慮風(fēng)場(chǎng)的情況下模擬準(zhǔn)確性不高。圖2為大豐測(cè)站2013年部分模擬值與實(shí)測(cè)值，可以看出除了實(shí)測(cè)數(shù)據(jù)有缺失的部分，模擬波高和實(shí)測(cè)數(shù)據(jù)擬合較好。圖3中給出了大豐測(cè)站模擬值與實(shí)測(cè)值兩者的線性比例的概率密度分布圖，實(shí)測(cè)數(shù)據(jù)與模擬值貼近45°線，說(shuō)明SWAN模型能夠較準(zhǔn)確地模擬本地區(qū)波高的時(shí)間變化。

圖2 大豐測(cè)站模擬波高和實(shí)測(cè)數(shù)據(jù)對(duì)比

Fig.2 Comparison of simulatedHsand measuredHsat Dafeng station

圖3 大豐測(cè)站模擬波高和實(shí)測(cè)數(shù)據(jù)線性比例的概率密度分布圖

Fig.3 Linearly scaled probability density function(pdf)of observedHsversus modeledHsat Dafeng station

1.3 極值分布函數(shù)

1.3.1 GEV模型

在極值理論中，已經(jīng)證明，對(duì)于足夠長(zhǎng)的獨(dú)立和相同分布的隨機(jī)變量序列，大小為n的樣本的最大值可以擬合到廣義極值(GEV)分布中[21]，該分布具有以下累積分布函數(shù)

(3)

μ、σ和ξ是位置、尺度和形狀參數(shù)，與重現(xiàn)期T相對(duì)應(yīng)的重現(xiàn)波高的求解可以使用以下方程

(4)

GEV經(jīng)常被應(yīng)用于極值的統(tǒng)計(jì)分布，最常用的方法是對(duì)于長(zhǎng)期數(shù)據(jù)，采用每年取一個(gè)極大值的方法進(jìn)行擬合，缺點(diǎn)是往往會(huì)忽略每年發(fā)生的其他極端情況。

1.3.2 基于POT方法的廣義帕累托分布模型(GP模型)

對(duì)于超過一定閾值的數(shù)據(jù)，滿足如下廣義帕累托分布函數(shù)

(5)

與重現(xiàn)期T相對(duì)應(yīng)的重現(xiàn)波高的求解如下

(6)

POT方法比每月取一個(gè)極值的方法更能反映出極端天氣的發(fā)生情況，但是需要慎重給定取樣時(shí)間間隔和閾值。時(shí)間間隔的選取應(yīng)結(jié)合不同地區(qū)的極端天氣發(fā)生情況，要保證盡可能多地考慮到極端天氣發(fā)生次數(shù)，且不至于在同一次極端天氣中選取多個(gè)極值。閾值的選取會(huì)直接決定取值的數(shù)量，過低的閾值會(huì)造成數(shù)據(jù)取樣數(shù)量大進(jìn)而導(dǎo)致模型預(yù)測(cè)結(jié)果的偏低，過高的閾值會(huì)造成分布參數(shù)計(jì)算結(jié)果的劇烈變動(dòng)[22]。

對(duì)于時(shí)間間隔的選取，已經(jīng)有學(xué)者研究發(fā)現(xiàn)，極端天氣事件的發(fā)生間隔是隨著不同地區(qū)、不同年限而相互獨(dú)立的[22-23]，因此只能確定一個(gè)常用的合理值，常用3 d[13-14]。

1.3.3 極值分布模型的適用性

本文以南黃海地區(qū)劃分的20個(gè)區(qū)域的40 a風(fēng)浪推算數(shù)據(jù)為基礎(chǔ)，采用基于取年極值與月極值的GEV分布模型和基于POT的GP分布模型，分析兩種模型的適用性，并推算對(duì)應(yīng)100 a重現(xiàn)期的重現(xiàn)波高。限于篇幅，本文僅選取了20個(gè)劃分區(qū)域中具有代表性5個(gè)區(qū)塊的進(jìn)行詳細(xì)說(shuō)明：輻射沙洲北部區(qū)塊N2、輻射沙洲中部區(qū)塊R2、輻射沙洲邊緣區(qū)塊P4共3個(gè)區(qū)塊(圖1-b)。

通常采用累積分布函數(shù)(CDF)圖和Q-Q圖來(lái)判斷數(shù)據(jù)與分布函數(shù)擬合的優(yōu)劣。GEV分布的CDF圖如圖4所示，當(dāng)各區(qū)域取年極值波高(圖4-a～圖4-c)和月極值波高(圖4-d～圖4-f)時(shí)，模擬值(虛線)與理論值(實(shí)線)均符合分布趨勢(shì)，說(shuō)明GEV分布適用于整個(gè)南黃海區(qū)域的波高擬合。在Q-Q圖(圖5)上可以看出年極值和月極值取樣基本符合分布趨勢(shì)，基本分布在45°線(虛線)附近，說(shuō)明兩者吻合良好。但在曲線尾部會(huì)有出入，預(yù)測(cè)曲線比實(shí)際偏低，即預(yù)測(cè)波高會(huì)小于實(shí)際值，這種現(xiàn)象在取月極值(圖5-d～圖5-f)時(shí)更明顯。

4-a N2區(qū)域取年極值 4-b R2區(qū)域取年極值 4-c P4區(qū)域取年極值

4-d N2區(qū)域取月極值 4-e R2區(qū)域取月極值 4-f P4區(qū)域取月極值

圖4 代表區(qū)塊GEV分布的CDF圖

Fig.4 The CDF plot of theGEVdistribution in the representative block

5-a N2區(qū)域取年極值 5-b R2區(qū)域取年極值 5-c P4區(qū)域取年極值

5-d N2區(qū)域取月極值 5-e R2區(qū)域取月極值 5-f P4區(qū)域取月極值

圖5 代表區(qū)塊GEV分布的Q-Q圖

Fig.5 The Quantile-Quantile(Q-Q)plot derived fromGEVdistribution in the representative block

基于POT方法的GP模型中時(shí)間間隔選為3 d，對(duì)于模型中閾值的選取，則根據(jù)超額均值函數(shù)確定。圖6所示為3個(gè)代表性區(qū)塊的超額均值函數(shù)圖。圖7所示為3個(gè)代表性區(qū)塊的平均剩余值圖，圖中實(shí)線代表形狀參數(shù)值Xi，上下兩條虛線為形狀參數(shù)估計(jì)值的95%置信區(qū)間的上下限。在超額均值函數(shù)和形狀參數(shù)滿足穩(wěn)定線性狀態(tài)后來(lái)看，3個(gè)區(qū)域閾值結(jié)果分別為1.9 m、1.85 m和2.45 m。閾值確定之后，繪制滿足超過閾值GP分布的擬合圖，如圖8所示，真實(shí)值由點(diǎn)組成，而黑色線表示分布趨勢(shì)，縱軸表示GP分布函數(shù)?？梢?，選擇最佳閾值后，在3個(gè)代表區(qū)塊中，超出相應(yīng)閾值的樣本均符合GP分布，除了中間部分稍微有所出入，頭尾兩部分預(yù)測(cè)狀況較好。

圖6 代表區(qū)塊超額均值函數(shù)圖(時(shí)間間隔為3 d)

Fig.6 The Mean Excess Function(MEF)in the representative block(Time span is △t=3 d)

圖7 代表區(qū)塊95%置信區(qū)間的平均剩余值圖(時(shí)間間隔為3 d)

Fig.7 The mean residual life plot with 95% confidence intervals in the representative block(Time span is △t=3 d)

圖8 代表區(qū)塊滿足超過閾值的GP分布(時(shí)間間隔為3 d)

Fig.8 TheGPdistribution function plot meeting the peak over threshold in the representative block(Time span is △t=3 d)

圖9 20個(gè)區(qū)塊取年極值、月極值方法與POT方法計(jì)算的百年重現(xiàn)波高對(duì)比

綜合之前采用取年、月極值的GEV分布和采用POT方法的GP分布來(lái)看，2種分布擬合均符合南黃海地區(qū)的波高分布趨勢(shì)。在所研究的分布均適用的情況下，分別用取年極值、月極值、POT方法獲得的數(shù)據(jù)，可以進(jìn)一步計(jì)算得到對(duì)應(yīng)分布函數(shù)下的重現(xiàn)波高。

2 結(jié)果

2.1 基于不同方法的百年重現(xiàn)波高空間分布

使用本文方法利用40 a的數(shù)據(jù)計(jì)算得到20個(gè)區(qū)塊的100 a重現(xiàn)波高統(tǒng)計(jì)如圖9所示。可見，采用GEV分布取年極值的方法比取月極值方法計(jì)算的百年重現(xiàn)波高明顯高，這是因?yàn)槟陿O值僅考慮了每年最大波高的情況，忽略了其他極端情況。POT方法采用的GP分布擬合出的百年重現(xiàn)波高和采用取年極值的GEV方法相差不大，互有高低，顯然也比取月極值的GEV方法計(jì)算值高，這是由于POT方法平均每年取的極值數(shù)約為3個(gè)，并沒有像取月極值那樣每年取12個(gè)之多，忽略相對(duì)小的極端波高。

為了直觀地表述南黃海地區(qū)的100 a重現(xiàn)波高的分布特征，圖10給出了3種不同取極值方法所計(jì)算的百年重現(xiàn)波高分布圖。其共同特征是，重現(xiàn)波高以輻射沙洲為中心向呈輻射狀外圍遞增。不同點(diǎn)是，重現(xiàn)波高的較大值集中分布區(qū)域有所不同，取年極值方法的重現(xiàn)波高最大值出現(xiàn)于東北部深水區(qū)，POT方法位于東南部區(qū)域但范圍較廣，月極值方法大致分布區(qū)域則較為均勻。究其原因，可能是不同區(qū)域出現(xiàn)極端波況的頻率和極端波高值各有不同。東北部地區(qū)更偏向于每年發(fā)生最極端情況，其他深水區(qū)更多發(fā)生次極端情況且發(fā)生次數(shù)更為頻繁。

10-a 年極值 10-b 月極值 10-c POT

圖10 三種方法分別計(jì)算的百年重現(xiàn)波高分布

Fig.10 The 100-year return wave height distribution calculated by three methods

圖11 每個(gè)區(qū)域重現(xiàn)波高計(jì)算方法的最保守選擇

在圖11中給出了3種方法中，每個(gè)區(qū)域得到百年重現(xiàn)波高最大值時(shí)，即工程上最保守時(shí)采用的方法，可見在南北部地區(qū)取年極值方法計(jì)算的百年重現(xiàn)波高大，其他中部地區(qū)和輻射沙洲區(qū)域則基本以POT方法為大，說(shuō)明這些地區(qū)在個(gè)別年份出現(xiàn)了多次高于正常年份的極端波況。

2.2 不同時(shí)間跨度對(duì)百年重現(xiàn)波高計(jì)算的影響

遵循工程上的取值思路，從最近的10 a數(shù)據(jù)開始，逐次增加數(shù)據(jù)長(zhǎng)度，使用對(duì)應(yīng)數(shù)據(jù)長(zhǎng)度為10 a、20 a、30 a和40 a的2009～2018年、1999～2018年、1989～2018年和1979～2018年期間的數(shù)據(jù)。采用基于年、月極值的GEV分布和基于POT方法的GP分布計(jì)算了輻射沙洲內(nèi)部的區(qū)域的百年重現(xiàn)波高，同時(shí)計(jì)算20～40 a相對(duì)10 a計(jì)算結(jié)果的增長(zhǎng)率(表2～表4)。

表2 基于年極值的不同時(shí)間跨度的百年重現(xiàn)波高及增長(zhǎng)率

表3 基于月極值的不同時(shí)間跨度的百年重現(xiàn)波高及增長(zhǎng)率

表4 基于POT方法的不同時(shí)間跨度的百年重現(xiàn)波高及增長(zhǎng)率

從不同取值長(zhǎng)度計(jì)算的百年重現(xiàn)波高增長(zhǎng)率看，取月極值的GEV分布計(jì)算結(jié)果對(duì)時(shí)間跨度的改變不敏感，而取年極值和POT方法均受制于時(shí)間跨度的選取?；谀陿O值的方法取近10 a的數(shù)據(jù)計(jì)算的百年重現(xiàn)波高最大，且不同數(shù)據(jù)長(zhǎng)度計(jì)算出來(lái)的值差別很大，這是因?yàn)槊磕耆∫粋€(gè)極值的情況偶然性大，所記錄的每一個(gè)極端波況對(duì)百年重現(xiàn)波高結(jié)果都會(huì)產(chǎn)生明顯的影響。取年極值的方法計(jì)算的重現(xiàn)波高值隨著取值年段長(zhǎng)度增加而減小，即較長(zhǎng)的數(shù)據(jù)集計(jì)算產(chǎn)生較小的重現(xiàn)波高，這與Mazas等[25]的發(fā)現(xiàn)相同，一方面可能是短期數(shù)據(jù)集的偶然性導(dǎo)致，另一方面可能是近10 a的極端天氣情況發(fā)生較之前更為頻繁。

取月極值的方法計(jì)算結(jié)果最為穩(wěn)定，但是由于取樣數(shù)量較多，影響相對(duì)小的波高也被考慮到，導(dǎo)致計(jì)算結(jié)果偏于不安全。POT方法的計(jì)算結(jié)果的穩(wěn)定程度介于前兩者之間。

取月極值方法計(jì)算的百年重現(xiàn)波高隨著取值年段長(zhǎng)度變化雖然穩(wěn)定，但是明顯偏小，取年極值計(jì)算結(jié)果最大，而POT方法計(jì)算的大小介于前二者之間。如果考慮到工程上偏于安全的情況，即在POT方法和取年極值方法中在不同地區(qū)選擇各自最安全的重現(xiàn)波高值。

3 結(jié)論

本文利用第三代近岸波浪模型SWAN在中國(guó)南黃海地區(qū)模擬了長(zhǎng)達(dá)40 a的波高數(shù)據(jù)，模擬波高和實(shí)測(cè)數(shù)據(jù)擬合較好。

基于年極值和月極值的廣義極值分布(GEV)和超閾值取值方法(POT)的廣義帕累托分布模型(GP)在南黃海地區(qū)適用。利用40 a模擬波高計(jì)算得到的百年重現(xiàn)波高在輻射沙洲北部地區(qū)計(jì)算的結(jié)果差別最大，采用月極值所得極值偏小，輻射沙洲南北外圍地區(qū)采用年極值計(jì)算的重現(xiàn)波高大，其余地區(qū)則用POT方法計(jì)算出的結(jié)果較大。如果考慮到工程上偏于安全的情況，即在POT方法和取年極值方法中在不同地區(qū)選擇各自最安全的重現(xiàn)波高值。

不同時(shí)間跨度的計(jì)算結(jié)果表明，取月極值的GEV分布計(jì)算結(jié)果對(duì)時(shí)間跨度的改變不敏感，而取年極值受之影響最大，POT方法介于兩者之間。取月極值方法計(jì)算的百年重現(xiàn)波高隨著取值時(shí)間跨度變化雖然穩(wěn)定，但是在不同時(shí)間跨度下百年重現(xiàn)波高值明顯偏小，取年極值的百年重現(xiàn)波高值最大，而POT方法計(jì)算的大小介于前二者之間。