高俊國, 翟金金, 董 勝

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

青島地處北溫帶季風(fēng)區(qū)域，受來自洋面上的東南季風(fēng)及海流、水團(tuán)的影響，空氣濕潤，年降水量充沛。降水對(duì)青島港口營運(yùn)影響很大，尤其是夏季的降雨。由于降雨停工、天晴突擊裝卸，加重了港口生產(chǎn)和船舶運(yùn)行的不均衡性，降低了港航效益。為了加強(qiáng)青島港口的營運(yùn)水平，提高港口的通過能力，有必要對(duì)青島地區(qū)歷年日最大降水量及其出現(xiàn)日期和歷年各月日最大降水量及其出現(xiàn)日期進(jìn)行深入研究。

迄今，專家學(xué)者對(duì)降雨量的極值估計(jì)進(jìn)行了研究：陳創(chuàng)買等[1]采用Gumbel分布和指數(shù)分布對(duì)廣東省日最大降雨量重現(xiàn)期進(jìn)行了計(jì)算。林巧美等[2]利用指數(shù)分布對(duì)揭陽市日最大降雨量和月最大降雨量進(jìn)行擬合，分析降雨量的年際演變特征及月際變化，進(jìn)而計(jì)算在不同重現(xiàn)期的極值。Cho等[3]利用Gamma分布和對(duì)數(shù)正態(tài)分布對(duì)衛(wèi)星得到的降雨率數(shù)據(jù)擬合，比較得知對(duì)數(shù)正態(tài)分布的結(jié)果更優(yōu)。Olofintoye等[4]選用Gumbel、對(duì)數(shù)正態(tài)、正態(tài)、對(duì)數(shù)-Gumbel、皮爾遜和對(duì)數(shù)皮爾遜分布對(duì)Nigeria地區(qū)的日最大降雨量特征進(jìn)行描述，檢驗(yàn)結(jié)果表明對(duì)數(shù)皮爾遜分布的擬合最優(yōu)。石先武等[5]以溫州站為例，利用5種極值模型模擬了臺(tái)風(fēng)降水概率，結(jié)果表明Weibull分布的整體擬合效果最佳。

降雨作為一種常見的水文事件，可以用降雨量、降雨天數(shù)和降雨強(qiáng)度等相互關(guān)聯(lián)的特征量來描述[6]。而傳統(tǒng)的單變量分布只能描述單一降雨指標(biāo)的概率變化特征，多變量聯(lián)合分布對(duì)研究降雨的綜合統(tǒng)計(jì)特征有重要意義。Yue[7]采用兩變量正態(tài)模型描述降雨特征量的聯(lián)合分布，實(shí)現(xiàn)非正態(tài)變量的正態(tài)化；Yue[8]采用二維Gumbel邏輯模型分析降雨峰值和降雨總量的聯(lián)合分布。以上成果要求對(duì)隨機(jī)變量進(jìn)行一定的假設(shè)，Copula函數(shù)的出現(xiàn)克服了這些局限，Salvadori[9]運(yùn)用二維Frank Copula函數(shù)模擬降雨歷時(shí)和降雨強(qiáng)度的聯(lián)合概率分布，推求聯(lián)合重現(xiàn)期等值線，引入第二重現(xiàn)期的概念，為風(fēng)險(xiǎn)評(píng)價(jià)提供科學(xué)依據(jù)。Lennartsson等[10]采用二維Gaussian Copula函數(shù)表征連續(xù)天數(shù)降雨之間的相關(guān)關(guān)系，結(jié)合馬爾科夫鏈，構(gòu)建了瑞典地區(qū)的降雨模型。Kao和Govindaraju[11]基于Plackett Copula函數(shù)構(gòu)造降雨量、歷時(shí)和峰值強(qiáng)度的二維和三維聯(lián)合分布模型，并與Frank Copula函數(shù)建立的降雨模型進(jìn)行對(duì)比，得出Plackett Copula函數(shù)可以克服Frank Copula不能描述相關(guān)性偏低的隨機(jī)變量和難以構(gòu)造高維聯(lián)合分布的缺點(diǎn)，能夠更全面地分析降雨特性。Zhang和Singh采用了三維Gumbel-Hougaard Copula函數(shù)建立了降雨歷時(shí)、降雨強(qiáng)度和雨深的三維聯(lián)合分布，計(jì)算了降雨事件的重現(xiàn)期[12]。Grimaldi和Serinaldi[13]采用7種三維Copula函數(shù)，建立了最大降雨強(qiáng)度、總降雨深和臨界雨深的聯(lián)合分布，確定了在已知重現(xiàn)期和臨界雨深時(shí)的總降雨深和最大降雨強(qiáng)度。Bardossy和Pegram[14]運(yùn)用非對(duì)稱Copula函數(shù)描述降雨量及其發(fā)生的內(nèi)在相關(guān)性，建立了多站降雨量的隨機(jī)模型。Copula函數(shù)在降雨事件中得到了廣泛的應(yīng)用。

對(duì)于青島地區(qū)降水的研究，多數(shù)集中在降水變化[15]以及日/月最大降水量的概率統(tǒng)計(jì)分析[16]上，缺乏對(duì)日最大降水量及其最可能出現(xiàn)日期(以下簡稱出現(xiàn)日期)的研究，而對(duì)于給定聯(lián)合重現(xiàn)期，且日最大降水量給定某一重現(xiàn)值時(shí)，能確定該值在一年中(或某個(gè)月份中)的可能發(fā)生時(shí)間，能夠有效地避開災(zāi)害風(fēng)險(xiǎn)，為提高港口的營運(yùn)水平，擴(kuò)大港口的通過能力有重要意義。因此，本文針對(duì)青島地區(qū)近30 a的歷年日最大降水量(RY)及其出現(xiàn)日期(LY)，以及各月的日最大降水量(RM)及其出現(xiàn)日期(LM)，基于Copula函數(shù)，分別建立二維聯(lián)合分布模型，對(duì)(RY,LY)和各月(RM,LM)進(jìn)行了分析討論，得到對(duì)工程有意義的結(jié)論。

1 二維統(tǒng)計(jì)模型的構(gòu)建

Copula函數(shù)是多個(gè)定義在[0,1]區(qū)間上均勻分布變量的聯(lián)合概率分布函數(shù)，它可以描述多維邊緣間的統(tǒng)計(jì)關(guān)系[9]?；赟klar定理，若F(x,y)是邊緣分布為FX(x)和FY(y)的二維聯(lián)合分布函數(shù)，則必然存在一個(gè)二維Copula函數(shù)C(u,v)滿足對(duì)任意的(x,y) ([-∞, +∞]2，有[17-19]

F(x,y)=C(u,v)=C(FX(x),FY(y))。

(1)

若FX(x)和FY(y)連續(xù)，則C(u,v)唯一；否則，C(u,v)由RanFX(RanFY唯一確定(Ran表示值域)。反之，若C(u,v)為二維Copula函數(shù)，u=FX(x)和v=FY(y)為邊緣分布，則由式(1)確定的函數(shù)F(x,y)即為X和Y的二維聯(lián)合分布函數(shù)[18]。

Copula函數(shù)主要分為橢圓型Copula和Archimedes Copulas等。其中，橢圓型Copula中的Gaussian Copula函數(shù)與常用的Nataf變換一致，其分布函數(shù)C(u,v)為[20]

(2)

式中：Φ(x)為標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)，Φ-1(x)為其逆；-1≤θ≤ 1為Φ-1(U)和Φ-1(V)的線性相關(guān)系數(shù)。θ= 0時(shí)，U和V獨(dú)立；|θ| = 1時(shí)，U和V完全相關(guān)。

Archimedean Copulas由其生成元確定的單參數(shù)函數(shù)，是目前水文統(tǒng)計(jì)中運(yùn)用較為廣泛的一類Copula函數(shù)。其中，最常用的二維Archimedean Copula主要有以下幾種[20-21]：

式中：u=FX(x)和v=FY(y)為邊緣分布；θ為Copula函數(shù)的待定參數(shù)。

聯(lián)合分布的參數(shù)估計(jì)分為兩步[19,21]：第1步，根據(jù)日最大降水量(RY,RM)及其出現(xiàn)日期(LY,LM)的觀測值序列，采用單變量分布的參數(shù)估計(jì)方法估計(jì)邊緣分布u和v的參數(shù)；第2步，采用極大似然法估計(jì)Copula函數(shù)的參數(shù)[22-23]。

2 歷年日最大降水量與出現(xiàn)日期的聯(lián)合統(tǒng)計(jì)分析

2.1 歷年日最大降水量及其出現(xiàn)日期的邊緣分布

青島地區(qū)歷年日最大降水量(Maximum daily rainfall:RY)與其在該年中的出現(xiàn)日期(Occurrence date:LY)的散點(diǎn)圖見圖1，其中LY指日最大降水量RY在該年中發(fā)生時(shí)，從1月1日開始累加到該天的總天數(shù)。采用Pearson-Ⅲ型分布[24-25]和最大熵分布[26-27]分別擬合日最大降水量及其出現(xiàn)日期的邊緣分布，

P-Ⅲ型分布的概率密度函數(shù)為

(7)

式中:a0為位置參數(shù)，且0

最大熵分布的概率密度函數(shù)為

(8)

式中:a0> 0為位置參數(shù)；β、γ與ε為待定參數(shù)。

首先采用式(7)對(duì)日最大降水量及其出現(xiàn)日期序列進(jìn)行Kolmogorov-Smirnov檢驗(yàn)(以下簡稱K-S檢驗(yàn))

(9)

圖1 青島地區(qū)RY和LY的散點(diǎn)圖Fig.1 The scatter plots ofRYandLY

經(jīng)過計(jì)算，歷年日最大降水量及其出現(xiàn)日期的K-S檢驗(yàn)統(tǒng)計(jì)量DRY=0.130 4和DLY=0.097 7均小于D30(0.05)=0.241 7，說明所選取的分布均可以擬合日最大降水量及其出現(xiàn)日期序列(見圖2)，并分別推算給定重現(xiàn)期時(shí)日最大降水量及其出現(xiàn)日期的重現(xiàn)值(見表1)。

2.2 歷年日最大降水量及其出現(xiàn)日期的聯(lián)合分布

日最大降水量及其出現(xiàn)日期分別采用Pearson-Ⅲ型分布和最大熵分布作為邊緣分布，基于二維Clayton Copula、Frank Copula、G-H Copula和Gaussian Copula函數(shù)，建立歷年日最大降水量(RY)及其出現(xiàn)日期(LY)的二維概率模型，采用Pearson’s χ2檢驗(yàn)對(duì)4種Copula函數(shù)進(jìn)行假設(shè)檢驗(yàn)，其檢驗(yàn)統(tǒng)計(jì)量M的估計(jì)值為[28-29]

(10)

圖2 日最大降水量及其出現(xiàn)日期的擬合分布Fig.2 The fitting distribution forRYandLY

重現(xiàn)期 Return period/a510202550100日最大降水量Maximum daily rainfall RY/mm138168196205232258出現(xiàn)日期Occurrence date LY/d233248260263272281

表2 Copulas函數(shù)的假設(shè)檢驗(yàn)和擬合優(yōu)度評(píng)價(jià)Table 2 Hypothesis testing and goodness-of-fit evaluation of Copulas

將4種Copula函數(shù)計(jì)算得到的理論聯(lián)合概率分布和經(jīng)驗(yàn)聯(lián)合聯(lián)合概率分別點(diǎn)繪在圖中(見圖3)，從圖中可以看出，數(shù)據(jù)點(diǎn)均分布在45°線附近，可以直觀地看出4種Copula函數(shù)擬合效果均較好。為了獲得最優(yōu)的擬合分布函數(shù)，結(jié)合均方根誤差法(RMSE)選擇一種Copula函數(shù)：

圖3 二維Copulas擬合優(yōu)度評(píng)價(jià)Fig.3 Goodness-of-fit test of bivariate Copulas

(11)

式中:n為樣本容量；F為理論分布計(jì)算值；P為經(jīng)驗(yàn)概率。RMSE越小，說明該Copula函數(shù)擬合的越好。

根據(jù)表2的擬合優(yōu)度評(píng)價(jià)，選取RMSE最小的Frank Copula函數(shù)作為聯(lián)合概率分布的連結(jié)函數(shù)，其計(jì)算得到的聯(lián)合重現(xiàn)期(RY和LY同時(shí)發(fā)生大于某閾值時(shí)聯(lián)合概率的倒數(shù))等值線見圖4。

(點(diǎn)為歷年日最大降水量及其出現(xiàn)日期的組合。Points denote the combinations ofRYandLY.)

圖4 日最大降水量及其出現(xiàn)日期的聯(lián)合重現(xiàn)期
Fig.4 Joint return periods ofRYandLY

由等值線(見圖4)可知，當(dāng)聯(lián)合重現(xiàn)期為某個(gè)值時(shí)，有無數(shù)組RY和LY與其對(duì)應(yīng)。此時(shí)給定一個(gè)日最大降水量，就可以得到該降水量相應(yīng)所處該年的確切位置。文中給出了當(dāng)聯(lián)合重現(xiàn)期分別為50和100年時(shí)，日最大降水量取5、10、20、25年一遇的重現(xiàn)值時(shí)的出現(xiàn)日期，具體結(jié)果見表3。

表3 (RY,LY)聯(lián)合重現(xiàn)期為50 a或100 a條件下，RY取不同重現(xiàn)值時(shí)的出現(xiàn)日期LYTable 3 LYcorresponding to a given univariate return valueRYunder joint return period (50 a or 100 a) of (RY,LY)

3 歷年各月的日最大降水量及出現(xiàn)日期的統(tǒng)計(jì)分析

歷年各月日最大降水量及其在該月中的出現(xiàn)日期分別為RMi和LMi(i=1～12)，選取擬合較優(yōu)的Pearson-Ⅲ型分布和最大熵分布分別擬合12個(gè)樣本的日最大降水量及其出現(xiàn)日期的邊緣分布，然后基于第2節(jié)中的4種二維Copula函數(shù)，建立聯(lián)合概率分布模型，通過Pearson′s χ2假設(shè)檢驗(yàn)和擬合優(yōu)度評(píng)價(jià)(RMSE)，均選取Frank Copula函數(shù)構(gòu)建歷年各月日最大降水量及其在該月中的出現(xiàn)日期的二維聯(lián)合概率分布模型。

表4 (RMi,LMi)聯(lián)合重現(xiàn)期為50、100 a，RMi取單因素重現(xiàn)值時(shí)對(duì)應(yīng)該月的出現(xiàn)日期Table 4 LMifor different joint return period: 50、100 a when theRMiis the given univariate return value

圖5 日最大降水量及其出現(xiàn)日期的邊緣擬合分布圖和聯(lián)合重現(xiàn)期等值線圖Fig.5 Marginal fitting distributions and joint return period contours of maximum daily rainfall and the associate occurrence date

由于篇幅的限制，文中僅給出春夏秋冬四個(gè)季節(jié)各一個(gè)月結(jié)果，RMi和LMi的邊緣擬合分布和聯(lián)合重現(xiàn)期等值線如圖5。由等值線圖可見，當(dāng)聯(lián)合重現(xiàn)期為某個(gè)值時(shí)，有無數(shù)組RMi和LMi與其對(duì)應(yīng)。此時(shí)給定一個(gè)日最大降水量，就可以得到該降水量在該月的出現(xiàn)日期(見表4)，其結(jié)果對(duì)避開降水大值區(qū)，合理安排港口作業(yè)時(shí)間，避免怠工損失具有十分重要的參考價(jià)值。

圖6給出了當(dāng)聯(lián)合重現(xiàn)期分別為50和100 a時(shí)，歷年各月日最大降水量取重現(xiàn)期為5、10、20、25、50 a的重現(xiàn)值時(shí)所對(duì)應(yīng)的其在該月的出現(xiàn)日期。

圖6 給定日最大降水量時(shí)，不同聯(lián)合重現(xiàn)期下的出現(xiàn)日期(12個(gè)月)Fig.6 LMifor different joint return periodswhenRMis the given univariate return value (12 months)

4 結(jié)論

收集青島地區(qū)近30年的日最大降水量與其在該年中的出現(xiàn)日期，以及各月的日最大降水量及其在該月中的出現(xiàn)日期，共計(jì)13個(gè)樣本序列為例，通過構(gòu)建二維聯(lián)合分布模型進(jìn)行統(tǒng)計(jì)分析，得到如下結(jié)論：

(1)針對(duì)青島地區(qū)歷年日最大降水量(RY)與其在該年中的出現(xiàn)日期(LY)，以及各月的日最大降水量(RM)及其在該月中的出現(xiàn)日期(LM)，分別建立了基于4種二維Copula函數(shù)(Clayton、Frank、Gumbel-Hougaard和Gaussian)的二維聯(lián)合分布模型，根據(jù)擬合優(yōu)度檢驗(yàn)，二維Frank Copula聯(lián)合分布模型最佳。

(2)通過二維聯(lián)合分布模型，對(duì)日最大降水量(RY)與其在該年中的出現(xiàn)日期(LY)的相關(guān)關(guān)系以及各月的日最大降水量(RM)與其在該月中的出現(xiàn)日期(LM)相關(guān)關(guān)系進(jìn)行研究，當(dāng)聯(lián)合重現(xiàn)期取定，且日最大降水量給定某一重現(xiàn)值時(shí)，可以估計(jì)該值在一年中(或某個(gè)月份中)的最可能出現(xiàn)的日期，因此能夠有效地避開降水大值區(qū)，選擇其他時(shí)段完成任務(wù)，從而為港口施工、建設(shè)和營運(yùn)提供指導(dǎo)性建議，為防范極端氣象情況做出合理、充分的準(zhǔn)備。