高俊國, 翟金金, 董 勝

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

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

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

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

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

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

Copula函數(shù)是多個定義在[0,1]區(qū)間上均勻分布變量的聯(lián)合概率分布函數(shù)，它可以描述多維邊緣間的統(tǒng)計關(guān)系[9]。基于Sklar定理，若F(x,y)是邊緣分布為FX(x)和FY(y)的二維聯(lián)合分布函數(shù)，則必然存在一個二維Copula函數(shù)C(u,v)滿足對任意的(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)為標準正態(tài)分布的分布函數(shù)，Φ-1(x)為其逆；-1≤θ≤ 1為Φ-1(U)和Φ-1(V)的線性相關(guān)系數(shù)。θ= 0時，U和V獨立；|θ| = 1時，U和V完全相關(guān)。

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

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

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

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

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

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

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

(7)

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

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

(8)

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

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

(9)

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

經(jīng)過計算，歷年日最大降水量及其出現(xiàn)日期的K-S檢驗統(tǒng)計量DRY=0.130 4和DLY=0.097 7均小于D30(0.05)=0.241 7，說明所選取的分布均可以擬合日最大降水量及其出現(xiàn)日期序列(見圖2)，并分別推算給定重現(xiàn)期時日最大降水量及其出現(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檢驗對4種Copula函數(shù)進行假設(shè)檢驗，其檢驗統(tǒng)計量M的估計值為[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ōu)度評價Table 2 Hypothesis testing and goodness-of-fit evaluation of Copulas

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

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

(11)

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

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

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

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

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

表3 (RY,LY)聯(lián)合重現(xiàn)期為50 a或100 a條件下，RY取不同重現(xiàn)值時的出現(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)計分析

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

表4 (RMi,LMi)聯(lián)合重現(xiàn)期為50、100 a，RMi取單因素重現(xiàn)值時對應(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

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

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

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

4 結(jié)論

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

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

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