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

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

1 引言

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

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

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

2 計(jì)算方法和資料

2.1 極值Ⅰ型頻率曲線分布

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

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

2.2 離散型Poisson分布

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

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

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

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

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

具體推算見文獻(xiàn)[6]。

2.4 過閾法

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

2.5 資料來源

臺(tái)風(fēng)信息采用中國氣象局的《臺(tái)風(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ù)報(bào)中心，在“908”專項(xiàng)臺(tái)風(fēng)風(fēng)暴潮災(zāi)害歷史數(shù)據(jù)集基礎(chǔ)上進(jìn)行延長和補(bǔ)充，增水值由臺(tái)風(fēng)影響期間測站實(shí)測逐時(shí)水位值去除天文潮得到，實(shí)測水位值根據(jù)國家標(biāo)準(zhǔn)[17]進(jìn)行觀測，數(shù)據(jù)經(jīng)過嚴(yán)格的質(zhì)量控制，排除了由于臺(tái)站遷址、儀器變更和地面沉降等產(chǎn)生的影響。

3 計(jì)算個(gè)例

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

圖1 部分代表驗(yàn)潮站分布圖Fig.1 Location of some typical water gauge

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

圖2 三沙站臺(tái)風(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 三沙站臺(tái)風(fēng)風(fēng)暴潮增水極值累計(jì)頻率曲線Fig.3 Accumulative frequency curve of typhoon storm surge extreme value in Sansha station

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

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

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

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

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

圖5 Gumbel分布概率檢驗(yàn)圖Fig.5 Gumbel distribution probability plot

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

4 比較分析

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

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

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

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

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

4.2 閾值敏感性試驗(yàn)

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

表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 分段計(jì)算結(jié)果比較

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

表5 隨機(jī)分段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 沿海代表站臺(tái)風(fēng)風(fēng)暴潮重現(xiàn)期分布

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

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

6 小結(jié)

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

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