邱 玥,龐 亮,董 勝

(中國(guó)海洋大學(xué)工程學(xué)院,山東 青島 266100)

引 言

海工結(jié)構(gòu)物設(shè)計(jì)標(biāo)準(zhǔn)規(guī)定了根據(jù)不同重現(xiàn)期環(huán)境條件計(jì)算海洋環(huán)境荷載的方法。即使用長(zhǎng)期觀測(cè)或后報(bào)資料,選取合適的概率模型進(jìn)行擬合,計(jì)算各種環(huán)境要素不同重現(xiàn)期的估計(jì)值,然后根據(jù)一定法則組合作用于海工建筑物,推求荷載。

此類(lèi)研究中,以極值分布等模型為工具的概率分析是最主要的方法。常見(jiàn)的極值分布模型有Gumbel分布、Fréchet分布、 Weibull分布,對(duì)數(shù)正態(tài)分布、P-Ⅲ曲線等模型,也廣泛的應(yīng)用于極值分析[1]。以上模型采用的都是年極值取樣方法,但進(jìn)行極端海洋環(huán)境概率分析時(shí)關(guān)心的往往是高分位數(shù)問(wèn)題。同樣的置信水平,樣本容量越小,得出結(jié)果的置信區(qū)間越寬,表明結(jié)果準(zhǔn)確度低,因此樣本容量的大小是影響重現(xiàn)期極值計(jì)算準(zhǔn)確度的最主要因素。而許多工程海域無(wú)法保證長(zhǎng)期、有效的實(shí)測(cè)數(shù)據(jù),使用傳統(tǒng)的年極值取樣方式,會(huì)因信息量不足導(dǎo)致計(jì)算偏差?？紤]到氣候波候波動(dòng)變化帶來(lái)的不確定性影響,如果觀測(cè)資料序列過(guò)短還會(huì)出現(xiàn)使用不同時(shí)間段的序列推求設(shè)計(jì)環(huán)境要素偏差明顯的問(wèn)題,即選取早期發(fā)生的樣本與選取近期發(fā)生的樣本得到的重現(xiàn)期海況計(jì)算結(jié)果會(huì)有很大的差別。解決樣本不足的方法有基于閾值法取樣(POT)的廣義Pareto模型(GPD),劉德輔等提出的采用過(guò)程取樣的復(fù)合極值分布模型[1-3]。

除了上述概率理論模型,借助數(shù)值模擬擴(kuò)充樣本也是解決觀測(cè)數(shù)據(jù)匱乏問(wèn)題的有效方法。Gatey[4]為計(jì)算50年重現(xiàn)期的風(fēng)速,通過(guò)壓力場(chǎng)使用Bratseth方法統(tǒng)計(jì)插值而得到背景風(fēng)場(chǎng)。Dukes等[5]通過(guò)構(gòu)造馬爾科夫鏈進(jìn)行數(shù)據(jù)仿真,生成了與觀測(cè)值相同時(shí)間長(zhǎng)度的序列,解決了長(zhǎng)期觀測(cè)中樣本缺失的問(wèn)題。研究表明當(dāng)所求重現(xiàn)期較短時(shí),理論方法與此方法具有相近的結(jié)果;對(duì)于超長(zhǎng)重現(xiàn)期(1 000年以上),此方法得出的結(jié)果偏低。Sanabria等[6]利用蒙特卡羅法模擬生成了上千年的澳大利亞海域的風(fēng)暴海況數(shù)據(jù),并證明了他們方法的合理性。龐亮[7]、劉德輔等[8]提出的基于隨機(jī)模擬與數(shù)值后報(bào)的雙層嵌套概率分析模式應(yīng)用于中國(guó)沿海風(fēng)暴海況的長(zhǎng)期概率分析得到了良好的效果。

為妥善解決海洋工程中因風(fēng)、浪、流等多種環(huán)境要素的聯(lián)合作用而產(chǎn)生的設(shè)計(jì)標(biāo)準(zhǔn)制定問(wèn)題,多元極值模型成為發(fā)展的趨勢(shì)。Gumbel等[9]最早提出了多元極值理論——對(duì)稱(chēng)Logistic模型,因其復(fù)雜性未能引起足夠的重視。20世紀(jì)90年代,Coles等[10]詳細(xì)研究了多元極值理論,提出了非對(duì)稱(chēng)Logistic模型、負(fù)非對(duì)稱(chēng)Logistic模型、Dirichlet模型,為工程界解決上述問(wèn)題提供了理論依據(jù)。但多元極值模型的表達(dá)式多為隱式形式,不便于工程應(yīng)用。為此,史道濟(jì)等[11]提出具有更廣泛應(yīng)用意義的三元嵌套Logistic模型及其顯式表達(dá)式。劉德輔等[1-3]提出多維復(fù)合極值分布模型。

綜上所述,目前在極端海況長(zhǎng)期概率分析領(lǐng)域,最需要解決的關(guān)鍵問(wèn)題就是樣本數(shù)據(jù)缺乏和不同環(huán)境要素合理組合。本文在復(fù)合極值分布理論的基礎(chǔ)上,提出基于短期觀測(cè)樣本的二項(xiàng)-廣義Pareto分布模型,并應(yīng)用于極端海況要素推算。

1 理論模型構(gòu)建

1.1 模型建立背景

極值理論包括BMM(Block Maximum Method)和POT(Peaks Over Threshold)兩種模型。其中,BMM模型將所有的樣本數(shù)據(jù)分組,選取每組最大值后,對(duì)重構(gòu)的樣本組進(jìn)行擬合估計(jì);POT模型則是將超過(guò)某一較大閾值所有數(shù)據(jù)作為樣本,在建模過(guò)程中將超過(guò)閾值的分布近似為廣義Pareto分布(GPD)[12]。已有學(xué)者對(duì)POT模型和其他數(shù)學(xué)模型在地震等巨災(zāi)事件的擬合進(jìn)行了比較,均得到POT模型在擬合方面效果更佳的結(jié)論[13],且POT模型在觀測(cè)值較少的情況下,估計(jì)的準(zhǔn)確性較高[14]。在實(shí)際問(wèn)題中,原始數(shù)據(jù)資料中超過(guò)閾值的樣本個(gè)數(shù)有限,其取值服從二項(xiàng)分布,預(yù)測(cè)方法可見(jiàn)《概率論與數(shù)理統(tǒng)計(jì)》一書(shū)中的附錄二所示[15]。因此,本文將POT模型篩選超閾值得到的二項(xiàng)分布與GPD模型得到的連續(xù)型極值分布進(jìn)行組合,構(gòu)造復(fù)合極值分布,即二項(xiàng)-廣義Pareto復(fù)合分布。

1.2 模型建立過(guò)程

(1)

Fμ(y)為超閾量的分布函數(shù),表示為:

(2)

對(duì)于充分大的閾值μ,超閾量Yi近似服從廣義Pareto分布,即:

(3)

式中,Gξσ(u)(y)為廣義Pareto分布,用于描述極端海況序列,例如風(fēng)速、波高、波周期等,得其分布函數(shù)為

(4)

式中,位置參數(shù)u∈R,尺度參數(shù)σ>0,形狀參數(shù)ξ∈R。

令x=y+μ,則(4)可變形為:

(5)

將上述離散型分布Y(x)與連續(xù)型極值分布Gξσ(u)(x)構(gòu)造復(fù)合極值分布,即二項(xiàng)-廣義Pareto復(fù)合分布,則其分布函數(shù)為:

由表2可知，各因素對(duì)提取辣椒堿含量影響的主次順序?yàn)椋好附鈺r(shí)間>酶解溫度>酶添加量>料液比>SSL的添加量>粒度；各因素對(duì)提取辣椒二氫堿含量影響的主次順序?yàn)椋好附鈺r(shí)間>酶添加量>酶解溫度>粒度>料液比>SSL的添加量，而最優(yōu)水平組合為A3B1C2D1E3F2，即纖維素酶添加量0.4%(W/W)、酶解時(shí)間4 h、酶解溫度45 ℃、料液比1∶8 (g/mL)、SSL的添加量1%(W/W)、辣椒粒度100目，此時(shí)，辣椒堿的含量可達(dá)7.82 mg/g，辣椒二氫堿的含量為4.94 mg/g。

(6)

代入各項(xiàng),得二項(xiàng)-廣義Pareto復(fù)合分布函數(shù)表達(dá)式為:

(7)

極端海洋工程環(huán)境的概率分析關(guān)心的是高分位數(shù)問(wèn)題,希望通過(guò)計(jì)算預(yù)測(cè)出多年一遇的極端波高,提前做好防范工作,以降低自然災(zāi)害給人類(lèi)帶來(lái)的影響[16]。

(8)

(9)

同時(shí),我們將基于年極值波高數(shù)據(jù),利用Gumbel極值分布、對(duì)數(shù)正態(tài)分布進(jìn)行多年一遇波高值的預(yù)測(cè),并將預(yù)測(cè)結(jié)果與上述二項(xiàng)-廣義Pareto復(fù)合分布預(yù)測(cè)的波高值進(jìn)行比較,以說(shuō)明二項(xiàng)-廣義Pareto復(fù)合分布模型更適用于巨災(zāi)事件的統(tǒng)計(jì)分析。其中,Gumbel極值分布和對(duì)數(shù)正態(tài)分布的設(shè)計(jì)波高計(jì)算公式如下:

Gumbel的分布函數(shù)為:

F(x)=exp{1-exp [-a(x-b)]}

(10)

(11)

對(duì)數(shù)正態(tài)分布的分布函數(shù)為:

(12)

多年一遇的設(shè)計(jì)波高為:

(13)

2 算例

2.1 閾值選取

閾值的合理選取關(guān)系到超閾值樣本的樣本數(shù)量和二項(xiàng)-廣義Pareto復(fù)合分布中的參數(shù)估計(jì),是POT模型成功擬合的關(guān)鍵。若閾值選取過(guò)小,超閾量分布與二項(xiàng)-廣義Pareto復(fù)合分布相差較大,則估計(jì)量會(huì)產(chǎn)生有偏估計(jì);若閾值選取過(guò)大,超閾值樣本數(shù)量個(gè)數(shù)減少,則影響模型的擬合效果[17]。

本文以中國(guó)青島某水文站1985—1989年日最大波高數(shù)據(jù)為樣本,運(yùn)用HILL圖法選取閾值,選取原則為找到圖形尾部相對(duì)穩(wěn)定的線段作為起始點(diǎn),則該點(diǎn)的橫坐標(biāo)對(duì)應(yīng)的數(shù)據(jù)作為閾值μ。

利用MATLAB繪制Hillplot,如圖1。

圖1 Hillplot圖

從圖1可以看出,中尾部首次達(dá)到平穩(wěn)狀態(tài)時(shí)橫坐標(biāo)對(duì)應(yīng)的值約為2.9。為了檢驗(yàn)閾值選在2.9附近是否合理,本文結(jié)合參數(shù)穩(wěn)定性估計(jì)圖來(lái)進(jìn)一步判斷。判斷方法為,在Hillplot圖確定的閾值范圍內(nèi),u=2.9作為初始閾值,在u=2.7到u=3.1之間均勻選取80個(gè)值作為閾值,觀察計(jì)算得到的極大似然估計(jì)值在區(qū)間內(nèi)是否保持相對(duì)穩(wěn)定。在上述閾值的選取范圍內(nèi),尺度參數(shù)和形狀參數(shù)的變化趨勢(shì)圖如圖2所示。

圖2 不同閾值下各參數(shù)的最大似然估計(jì)值

從圖2可以看出,尺度參數(shù)和形狀參數(shù)均在(2.8,2.995)內(nèi)保持相對(duì)穩(wěn)定,為確保POT模型擬合的準(zhǔn)確性,本文選取相對(duì)穩(wěn)定區(qū)間內(nèi)的較大值作為最終閾值,即u=2.995。

2.2 參數(shù)估計(jì)

在得到最佳閾值后,利用極大似然估計(jì)法對(duì)自定義的二項(xiàng)-廣義Pareto分布函數(shù)中的各個(gè)統(tǒng)計(jì)參數(shù)進(jìn)行估計(jì)。該復(fù)合分布的各參數(shù)見(jiàn)值表1。

表1 二項(xiàng)-廣義Pareto復(fù)合分布統(tǒng)計(jì)參數(shù)

(14)

基于上述參數(shù)估計(jì)結(jié)果以及分布函數(shù),為檢驗(yàn)二項(xiàng)-廣義Pareto分布模型是否具有較好的擬合效果,我們繪制了項(xiàng)-廣義Pareto擬合曲線以及樣本數(shù)據(jù)的經(jīng)驗(yàn)分布曲線進(jìn)行診斷。擬合曲線如圖3所示。

圖3 經(jīng)驗(yàn)累積分布曲線及二項(xiàng)-廣義Pareto擬合曲線

觀察圖3可以發(fā)現(xiàn),參數(shù)估計(jì)結(jié)果良好,由此構(gòu)造的二項(xiàng)-廣義Pareto復(fù)合極值分布模型對(duì)日波高最大值序列具有很好的擬合效果。

為預(yù)測(cè)出多年一遇的極端波高,我們將上述各個(gè)統(tǒng)計(jì)參數(shù)代入式(9),得到二項(xiàng)-廣義Pareto復(fù)合分布的分位數(shù)函數(shù),如下:

(15)

2.3 多年一遇波高預(yù)測(cè)以及方法比較

二項(xiàng)-廣義Pareto模型是能夠使用短期觀測(cè)資料推算較長(zhǎng)重現(xiàn)期海況設(shè)計(jì)值的方法。而Gumbel模型、對(duì)數(shù)正態(tài)模型等是基于年極值取樣的概率分析法,在計(jì)算不同重現(xiàn)期的設(shè)計(jì)值時(shí)都要求20～30年不等的長(zhǎng)期海況資料,而這在很多海域無(wú)法做到。

本文利用青島5年(1985—1989)的日最大波高數(shù)據(jù)資料,選用二項(xiàng)-廣義Pareto模型計(jì)算多年一遇波高設(shè)計(jì)值。同時(shí)用青島20年(1970—1989)的年極值波高數(shù)據(jù)資料,選用Gumbel模型、對(duì)數(shù)正態(tài)模型計(jì)算多年一遇波高設(shè)計(jì)值。將這三種方法所得的結(jié)果進(jìn)行比較,如表2。

表2 不同分布下的不同重現(xiàn)期的重現(xiàn)水平

通過(guò)分析表2,對(duì)比50年一遇的波高重現(xiàn)期極值,我們可以發(fā)現(xiàn)自定義的二項(xiàng)-廣義Pareto模型計(jì)算得出的結(jié)果較小,Gumbel模型、對(duì)數(shù)正態(tài)模型得到的結(jié)果近似相等。因此,我們推斷基于較短樣本資料,利用自定義模型計(jì)算出的50年一遇波高極值適用于設(shè)計(jì)標(biāo)準(zhǔn)低、服役年限短的海洋水工建筑物。此外,觀察推求出的百年一遇、150年一遇、200年一遇的波高極值,可以發(fā)現(xiàn)利用自定義的二項(xiàng)-廣義Pareto模型,采用短期數(shù)據(jù)資料預(yù)測(cè)波高值得到的結(jié)果,與采用長(zhǎng)期數(shù)據(jù)資料,利用Gumbel模型、對(duì)數(shù)正態(tài)模型得到的結(jié)果相差不大。因此,自定義的二項(xiàng)-廣義Pareto模型在預(yù)測(cè)波高方面是適用的。

此外可以發(fā)現(xiàn),由于現(xiàn)有的數(shù)據(jù)很難達(dá)到完全精確和充足,導(dǎo)致各種復(fù)雜的統(tǒng)計(jì)模型也無(wú)法達(dá)到完全精確,因此沒(méi)有一種模型可以在所有情況下都適用。

3 結(jié)論

本文基于POT模型,將POT模型篩選超閾值得到的二項(xiàng)分布與GPD模型得到的連續(xù)型極值分布進(jìn)行組合,構(gòu)造新的復(fù)合極值分布,即二項(xiàng)-廣義Pareto復(fù)合分布,并應(yīng)用于極端海況要素推算,得出以下主要結(jié)論:

(1)二項(xiàng)-廣義Pareto復(fù)合分布模型具有良好的擬合效果,可以較好的描述多年一遇波高分布,彌補(bǔ)了傳統(tǒng)方法需要長(zhǎng)期原始海況數(shù)據(jù)的缺陷。

(2)運(yùn)用該模型對(duì)青島海域5年(1985—1989)的日最大波高數(shù)據(jù)進(jìn)行實(shí)例分析,結(jié)果顯示模型擬合較好,并分別預(yù)測(cè)了50、100、150、200年一遇的極值波高。將上述模型擬合結(jié)果與Gumbel模型、對(duì)數(shù)正態(tài)模型利用青島海域20年(1970—1989)數(shù)據(jù)預(yù)測(cè)的50、100、150、200年一遇的極值波高結(jié)果做比較研究,結(jié)果指出:利用自定義的二項(xiàng)-廣義Pareto模型,采用短期數(shù)據(jù)資料預(yù)測(cè)波高值得到的效果,與采用長(zhǎng)期數(shù)據(jù)資料,利用Gumbel、對(duì)數(shù)正態(tài)模型得到的結(jié)果相差不大。得出自定義的二項(xiàng)-廣義Pareto模型在預(yù)測(cè)波高方面是適用的。

(3)由于現(xiàn)有的數(shù)據(jù)很難達(dá)到完全精確,導(dǎo)致各種復(fù)雜的統(tǒng)計(jì)模型也無(wú)法達(dá)到完全精確,因此沒(méi)有一種方法可以在所有情況下適用。