付翔，梁森棟，李濤，吳少華

（國家海洋環(huán)境預(yù)報中心，北京 100081）

1 引言

風(fēng)暴潮災(zāi)害是我國最主要的海洋災(zāi)害，特別是臺風(fēng)風(fēng)暴潮。自1989年有統(tǒng)計數(shù)據(jù)以來，風(fēng)暴潮（含近岸浪）災(zāi)害造成的直接經(jīng)濟損失占全部海洋災(zāi)害直接經(jīng)濟損失的90%以上，嚴(yán)重的臺風(fēng)風(fēng)暴潮災(zāi)害幾乎每年都會在不同地區(qū)發(fā)生[1]。因此，在涉及海洋工程建設(shè)、海岸環(huán)境損害評價、保險的災(zāi)害風(fēng)險計算以及氣候變化影響評估時，需考慮風(fēng)暴潮的致災(zāi)危險性。風(fēng)暴潮致災(zāi)因子的危險性體現(xiàn)在風(fēng)暴潮的強度、頻率以及疊加浪、涌和天文潮后引起的高潮位上，極值發(fā)生頻率（P）是常用的定量化指標(biāo)之一[2-3]。

極值發(fā)生頻率和重現(xiàn)期（T= 1/P）傳遞的是極端事件最大可能的信息。它有兩層含義：一是對于某頻率p的極值，下次出現(xiàn)該值的平均等待時間為Ta；二是在Ta 里出現(xiàn)該極值的次數(shù)平均是1 次[4]。對于臺風(fēng)風(fēng)暴潮，由于臺風(fēng)影響具有隨機性，某一區(qū)域并不是每年都遭受顯著的臺風(fēng)風(fēng)暴潮，但在某些年卻可能遭受一次以上的臺風(fēng)風(fēng)暴潮災(zāi)害，例如2014 年連續(xù)遭受臺風(fēng)“威馬遜”和臺風(fēng)“海鷗”襲擊的雷州半島，以及2015 年臺風(fēng)“蘇迪羅”和臺風(fēng)“杜鵑”連續(xù)登陸的福建中部。因此，如果使用年極值法來計算風(fēng)暴潮重現(xiàn)期，每年只取一個極值，有可能會取到無臺風(fēng)或弱臺風(fēng)影響年的小值而舍去了頻發(fā)年的次大值，使得計算數(shù)據(jù)失真，估算的重現(xiàn)期極值出現(xiàn)較大偏差?？紤]隨機事件發(fā)生次數(shù)的復(fù)合極值分布可以有效地解決這一問題。

引入一個描述無規(guī)律事件發(fā)生次數(shù)的離散型概率分布——Poisson 分布，使其與一個連續(xù)型的極值分布復(fù)合，組成的復(fù)合極值分布不僅能擬合每年實測變量隨機的數(shù)據(jù)序列，對于資料年限較短的情況也能可靠應(yīng)用[5]。復(fù)合極值分布提出后，不少學(xué)者在臺風(fēng)波高、極值風(fēng)速和風(fēng)暴潮的計算預(yù)測中進行了研究和應(yīng)用[6-10]。目前Poisson-Gumbel 復(fù)合極值分布已被推薦用于臺風(fēng)波高的設(shè)計值計算[11-12]。本文將該分布法用于臺風(fēng)風(fēng)暴潮的重現(xiàn)期計算，分析其優(yōu)缺點及應(yīng)用效果。

2 計算方法和資料

2.1 極值Ⅰ型頻率曲線分布

港口與海岸工程設(shè)計中，常采用極值Ⅰ型分布律推求指定頻率或重現(xiàn)期的風(fēng)暴潮和極端高低潮位值[11-12]。它是廣義極值分布（Fisher-Tippett）的一種，最早由Gumbel[13]用于水文統(tǒng)計。其概率分布函數(shù)為雙指數(shù)形式，表示為：

式中，μ為位置參數(shù)，即分布的眾數(shù)；α為尺度參數(shù)，代表分布的離散性。這兩個參數(shù)可用多種方法進行估計[14-15]。

2.2 離散型Poisson分布

連續(xù)的極值Ⅰ型分布律僅能用于規(guī)律的極值樣本，如年極值或月極值等。而對于某一地區(qū)，每年臺風(fēng)風(fēng)暴潮的發(fā)生次數(shù)和強度都是隨機的，它們構(gòu)成的是離散型分布。描述稀有事件發(fā)生次數(shù)的Poisson 分布非常適合用以描述臺風(fēng)風(fēng)暴潮的發(fā)生概率。假定某地一定強度以上的臺風(fēng)風(fēng)暴潮的發(fā)生頻次n符合Poisson分布，其概率質(zhì)量函數(shù)表示為：

式中，λ為年平均臺風(fēng)風(fēng)暴潮發(fā)生次數(shù)。

2.3 Poisson-Gumbel復(fù)合極值分布

根據(jù)復(fù)合極值分布式[5-6]，給定設(shè)計頻率p′的極值分布為：

式中，p= 1 -p′。定義稱為重現(xiàn)期。將式（1）、（2）代入（3），則某T年一遇的臺風(fēng)風(fēng)暴潮值xp可表示為：

具體推算見文獻[6]。

2.4 過閾法

過閾法（Peak Over Threshold，POT）是極值理論下確定計算樣本序列的一種方法，即先確定某一觀測值下限（閾值），凡超過此值的均列入統(tǒng)計序列，而后再進行極值分布擬合。閾值的確定十分關(guān)鍵，閾值過高會使可利用的樣本數(shù)過少，從而導(dǎo)致分布函數(shù)估計的參數(shù)方差偏高；閾值過低又會使分布的漸進性得不到滿足，產(chǎn)生有偏估計。閾值的選取有多種方法，如平均剩余壽命圖法（Mean Residual Life Plot）、平均超越期望圖法（Mean Excess Plot）、Hill 圖（Hill plot）等。計算法包括峰度法和漸進均方誤差法等。

2.5 資料來源

臺風(fēng)信息采用中國氣象局的《臺風(fēng)年鑒》（1980—1988 年）和《熱帶氣旋年鑒》（1989—2016 年），以及熱帶氣旋資料中心最佳路徑數(shù)據(jù)集[16]（網(wǎng)址：http：//tcdata.typhoon.org.cn）中的資料（2017—2019年）。風(fēng)暴潮數(shù)據(jù)來自國家海洋環(huán)境預(yù)報中心，在“908”專項臺風(fēng)風(fēng)暴潮災(zāi)害歷史數(shù)據(jù)集基礎(chǔ)上進行延長和補充，增水值由臺風(fēng)影響期間測站實測逐時水位值去除天文潮得到，實測水位值根據(jù)國家標(biāo)準(zhǔn)[17]進行觀測，數(shù)據(jù)經(jīng)過嚴(yán)格的質(zhì)量控制，排除了由于臺站遷址、儀器變更和地面沉降等產(chǎn)生的影響。

3 計算個例

三沙驗潮站位于福建省東北部（見圖1），幾乎每年都會受到臺風(fēng)影響。以其為例，利用Poisson-Gumbel復(fù)合極值分布計算其多年一遇臺風(fēng)風(fēng)暴潮。

圖1 部分代表驗潮站分布圖Fig.1 Location of some typical water gauge

第一步：確定閾值。采用平均剩余壽命圖法[17]和百分位法確定閾值范圍，挑選1980—2019 年間進入距三沙站約500 km 范圍以內(nèi)的臺風(fēng)，共計171個，統(tǒng)計每次過程的最大風(fēng)暴潮值，剔除其中的減水過程后建立經(jīng)驗頻率序列。根據(jù)平均剩余壽命圖（見圖2），考慮到近似線性和置信區(qū)間范圍[18]，閾值可在20～80 cm 之間選取。取其中值，同時也是增水極值累積頻率的中位數(shù)50 cm 為統(tǒng)計閾值（見圖3），使其發(fā)生次數(shù)符合Poisson分布。

圖2 三沙站臺風(fēng)風(fēng)暴潮增水極值序列平均剩余壽命圖（虛線為95%置信區(qū)間）Fig.2 Mean residual life of typhoon storm surge extreme value in Sansha station（The dotted line is the 95%confidence interval）

圖3 三沙站臺風(fēng)風(fēng)暴潮增水極值累計頻率曲線Fig.3 Accumulative frequency curve of typhoon storm surge extreme value in Sansha station

第二步：Poisson 分布檢驗。統(tǒng)計三沙站超過閾值的臺風(fēng)風(fēng)暴潮過程，由式（2）計算其發(fā)生次數(shù)的理論頻率分布（見圖4），采用卡方分布χ2=進行檢驗。取顯著水平0.05，查自由度為5（k- 1）的卡方分布臨界表= 11.07，由χ2＜χ02.05可知三沙站最大增水50 cm以上的臺風(fēng)風(fēng)暴潮發(fā)生頻數(shù)符合Poisson分布。

圖4 三沙站臺風(fēng)風(fēng)暴潮發(fā)生頻次及概率密度分布Fig.4 Distribution of probability density and typhoon storm surge occurrence frequencies in Sansha station

如果發(fā)生頻數(shù)擬合不能通過檢驗，則應(yīng)該適當(dāng)提高閾值大小減少樣本量，再進行表1計算，直至通過檢驗。

表1 臺風(fēng)風(fēng)暴潮發(fā)生次數(shù)擬合檢驗Tab.1 Fitting test of the number of typhoon storm surge events

第三步：Gumbel 分布檢驗。P-P 圖是變量的累計比例與指定分布的累積比例之間的關(guān)系圖，可直觀地檢驗樣本數(shù)據(jù)是否符合指定概率分布。當(dāng)數(shù)據(jù)符合指定分布時，P-P圖中各點近似呈一條直線，即代表樣本數(shù)據(jù)的點基本在代表理論分布的對角線上。根據(jù)P-P 圖，三沙站過閾臺風(fēng)風(fēng)暴潮極值序列符合Gumbel 分布（見圖5），也可采用Kolmogorov-Smirnov 檢驗進行數(shù)據(jù)序列的Gumbel分布檢驗。

圖5 Gumbel分布概率檢驗圖Fig.5 Gumbel distribution probability plot

第四步：復(fù)合極值分布擬合。根據(jù)式（4）擬合三沙站過閾風(fēng)暴潮增水極值序列，利用有限樣本分布函數(shù)代替總體分布函數(shù)的Gumbel 法[13]估算參數(shù)μ和α，然后計算不同重現(xiàn)期的臺風(fēng)風(fēng)暴潮極端增水值。

4 比較分析

4.1 與年極值序列計算結(jié)果比較

建立三沙站1980—2019 年的臺風(fēng)風(fēng)暴潮年極值增水序列，利用Gumbel 分布擬合計算三沙站臺風(fēng)風(fēng)暴潮多年一遇重現(xiàn)期增水值，與復(fù)合極值分布擬合的結(jié)果進行比較（見表2）。結(jié)果顯示，對于10年一遇以下增水，Gumbel計算的增水值比復(fù)合分布擬合的結(jié)果小，而對于10 年一遇以上增水，Gumbel計算值比復(fù)合分布擬合的結(jié)果大，且隨著重現(xiàn)期增大，差比也增大。在中高頻次發(fā)生范圍內(nèi)，兩種分布的計算值比較接近，但在高頻次和低頻次發(fā)生范圍內(nèi)的計算差異較大。究其原因，采用年極值序列擬合計算時，由于無臺風(fēng)影響年低值的存在會導(dǎo)致Gumbel 曲線的斜率增大，使得高頻發(fā)生值更小，低頻發(fā)生值更大；而過閾法序列不僅包含了在相同年份發(fā)生的高值，還剔除了弱過程，避免了增水小值帶來的較大觀測誤差和潮汐預(yù)報誤差，但由于相對低值樣本增多，導(dǎo)致低頻次發(fā)生的風(fēng)暴潮值計算偏小。

表2 不同分布擬合的重現(xiàn)期增水值比較（單位：cm）Tab.2 Comparison of surge return values between different distribution（unit:cm）

根據(jù)三沙站臺風(fēng)風(fēng)暴潮原始統(tǒng)計序列，40 a 間增水大于79 cm 的過程為33 次，約為1.2 年一遇；增水大于90 cm 以上的過程為23 次，約為1.7 年一遇；增水大于95 cm 的過程為20 次，正為2 年一遇。可見Poisson-Gumbel復(fù)合極值分布的高頻次發(fā)生值的結(jié)果更為合理。另外，由圖6也可看出，對于過閾原始序列的擬合，復(fù)合極值分布對極大值的擬合結(jié)果明顯好于Gumbel分布。

圖6 分布擬合比較Fig.6 Comparison between curves of different distribution

4.2 閾值敏感性試驗

由于接近閾值的小值數(shù)量的多少對擬合估算結(jié)果影響較大，因此需要調(diào)整閾值大小，比較不同閾值下樣本序列的復(fù)合分布擬合結(jié)果。取最大增水閾值60 cm、65 cm、70 cm 和最高75 cm，分別截取臺風(fēng)風(fēng)暴潮過程樣本序列進行復(fù)合分布計算。所有序列的發(fā)生頻數(shù)分布均通過0.05 顯著性檢驗的Poisson 分布。結(jié)果顯示（見表3），閾值選取增加，所有重現(xiàn)期的計算值都會升高，在閾值為70 cm時達到最高值，后又開始下降。文獻資料顯示[19]，三沙站1949 年以來臺風(fēng)風(fēng)暴潮極值為190～200 cm（6614 號臺風(fēng)風(fēng)暴潮），與閾值為70 cm 和75 cm 的計算結(jié)果相近。由此可見，雖然高頻次發(fā)生值的估計偏高，但對于低頻范圍的高年遇值，高閾值序列的擬合效果更好。

表3 擬合不同閾值樣本的Possion-Gumbel復(fù)合分布重現(xiàn)期增水值比較（單位：cm）Tab.3 Comparison of surge return values of Poisson-Gumbel distributions fitting samples with different threshold（unit:cm）

4.3 分段計算結(jié)果比較

根據(jù)相關(guān)標(biāo)準(zhǔn)[11-12]，在進行極端水位值確定時，應(yīng)用不少于連續(xù)20 a的潮位實測資料進行計算。將三沙站40 a 的資料序列分為前20 a 和后20 a，以閾值70 cm 為例（其余結(jié)果略），分別對兩種方法的擬合計算進行比較，結(jié)果見表4。我們發(fā)現(xiàn)，兩種方法兩個時間段的差值均自高頻向低頻增大，復(fù)合分布結(jié)果僅在高頻處差值較小，10 年一遇以上低頻結(jié)果兩個時間段的差值的復(fù)合分布遠高于Gumbel 年極值分布。其原因在于三沙站這兩個時間段的臺風(fēng)風(fēng)暴潮極值序列差別較大，前20 a 極值偏小且分布均勻，以閾值70 cm 為例，前20 a 的極值樣本僅為18個，均值為89.8 cm，標(biāo)準(zhǔn)差為16.8 cm；而后20 a的樣本為26 個，極值均值為99.7 cm，標(biāo)準(zhǔn)差為27.3 cm。樣本數(shù)據(jù)的較大差異導(dǎo)致了復(fù)合極值分布對低頻發(fā)生值的估算差別顯著，前20 a 數(shù)據(jù)的計算值嚴(yán)重偏低。

表4 連續(xù)分段20 a樣本數(shù)據(jù)重現(xiàn)期增水值比較（單位：cm）Tab.4 Comparison of surge return values between different distribution of consecutive 20 a subsection data(unit:cm)

在包含有40 a間最大增水的原始過閾樣本和年極值樣本中，任意挑選兩組20 a 的數(shù)據(jù)序列進行兩種分布擬合計算，由于觀測年代較短，復(fù)合分布應(yīng)取低閾值以增加樣本數(shù)量，這里取50 cm。由兩種分布的兩組結(jié)果比較可見（見表5），在觀測年代較短且樣本較少的情況下，復(fù)合分布的擬合結(jié)果比單一的年極值Gumbel 分布結(jié)果更為穩(wěn)定，兩組數(shù)據(jù)結(jié)果的差值更小，低頻值也更為合理。

表5 隨機分段20 a樣本數(shù)據(jù)重現(xiàn)期增水值比較（單位：cm）Tab.5 Comparison of surge return values between different distribution of random 20 a subsection data(unit:cm)

5 沿海代表站臺風(fēng)風(fēng)暴潮重現(xiàn)期分布

根據(jù)第2 部分的計算步驟，選取部分代表驗潮站（見圖1），利用Poisson-Gumbel 復(fù)合分布擬合1980—2019 年40 a 的臺風(fēng)風(fēng)暴潮增水極值序列，計算出各站50 年一遇和100 年一遇的臺風(fēng)風(fēng)暴潮增水值，閾值選取均在累積分布中位數(shù)附近并保證序列通過分布檢驗。由表6 可見，塘沽、吳淞、閘坡和秀英站可能遭受的臺風(fēng)風(fēng)暴潮強度較強，不僅增水值大，100年一遇和50年一遇增水的差值也大，發(fā)生頻次低的風(fēng)暴潮強度更強，說明臺風(fēng)在這些地方可能引起的極端風(fēng)暴潮更高，危險性也更大。

表6 各站臺風(fēng)風(fēng)暴潮重現(xiàn)期值（單位：cm）Tab.6 Return values of typhoon storm surge in several typical stations(unit:cm)

6 小結(jié)

考慮隨機事件發(fā)生頻率的Poisson-Gumbel復(fù)合極值分布是解決觀測資料年限較短或是出現(xiàn)無年極值（無臺風(fēng)影響）情況下，風(fēng)暴潮極值重現(xiàn)期計算的有效方法。復(fù)合極值分布的計算結(jié)果對樣本數(shù)據(jù)較為敏感，特別是標(biāo)準(zhǔn)偏差較小的樣本數(shù)據(jù)，低頻次發(fā)生值的計算結(jié)果顯著偏低。閾值的選取對結(jié)果影響較大，對于較長時間序列數(shù)據(jù)，復(fù)合分布的閾值取值偏低則高頻次發(fā)生值更準(zhǔn)，閾值取值偏高則低頻次發(fā)生值更準(zhǔn)；而對于短時間序列數(shù)據(jù)，復(fù)合分布的計算結(jié)果更穩(wěn)定，且對低頻發(fā)生值的估算也更準(zhǔn)，此時由于樣本數(shù)量較少，應(yīng)取偏小閾值以增加樣本數(shù)量。對幾個典型潮位站的重現(xiàn)期計算表明，塘沽、吳淞、閘坡和秀英站臺風(fēng)引起的可能極端風(fēng)暴潮較強，臺風(fēng)風(fēng)暴潮致災(zāi)危險性較高。

復(fù)合分布在風(fēng)暴潮極值計算中的特征可延伸至潮位極值計算。但由于其對數(shù)據(jù)樣本的敏感性，在工程應(yīng)用中，仍建議以差比法延長資料年限后再行計算。若考慮無年極值或一年多極值的情況，則需根據(jù)對高頻值或低頻值的需求，確定具體閾值。