張大偉，張 鵬，郭 珊，吳輝明

(1.廣州珠科院工程勘察設(shè)計(jì)有限公司，廣州 510610；2．珠江水利委員會(huì)珠江水利科學(xué)研究院，廣州 510610)

珠海市位于珠江出海口西岸，瀕臨南海，大部分為沖擊平原區(qū)，地勢較低，受洪潮交匯作用、水文情勢相當(dāng)復(fù)雜，屬于典型的平原感潮河網(wǎng)地區(qū)。珠海市受特有的地理因素影響，天文大潮、風(fēng)暴潮與暴雨等災(zāi)害性天氣頻發(fā)，時(shí)常發(fā)生內(nèi)澇，損失嚴(yán)重。當(dāng)降雨量較大并遭遇外海高潮位時(shí)，澇水不能及時(shí)排出，引發(fā)洪澇災(zāi)害，造成重大經(jīng)濟(jì)損失，嚴(yán)重威脅著人民生命財(cái)產(chǎn)安全。中珠聯(lián)圍內(nèi)澇與暴雨和潮水位密切相關(guān)，因此，研究中珠聯(lián)圍降雨與潮位遭遇聯(lián)合分布，分析內(nèi)澇發(fā)生的頻率，為中珠聯(lián)圍內(nèi)澇防治提供技術(shù)支撐。

目前，關(guān)于珠海市暴雨與潮位遭遇聯(lián)合分布研究甚少，遭遇組合常采用統(tǒng)計(jì)分析和外包法處理[1]。由于暴雨與潮水位的空間相關(guān)性，區(qū)域暴雨與潮水位總有一定的聯(lián)系，因此它們的聯(lián)合分布函數(shù)不能將2者的邊緣分布函數(shù)簡單相乘而得到[2]。傳統(tǒng)的雨潮遭遇統(tǒng)計(jì)分析方法無法考慮雨潮之間遭遇的內(nèi)在聯(lián)系，而Copula函數(shù)理論可以通過構(gòu)造聯(lián)合分布來描述降雨和潮位之間的內(nèi)在聯(lián)系，Copula理論的出現(xiàn)為構(gòu)造多變量事件的相關(guān)關(guān)系和聯(lián)合分布提供了一種新的理論計(jì)算方法[3]。該方法最早于1959年提出，20世紀(jì)末才被運(yùn)用到金融領(lǐng)域，21世紀(jì)初才引入工程領(lǐng)域。近10 a多來，在多變量水文因素分析中，國內(nèi)外[4-6]學(xué)者開始采用Copula函數(shù)來構(gòu)建多變量的聯(lián)合分布函數(shù)，研究多變量水文因素概率分布。林榮等[7]采用正態(tài)變換法分析了吳淞口設(shè)計(jì)潮位與不同頻率區(qū)間降雨及相應(yīng)黃浦江曲江降雨遭遇的概率，然而正態(tài)變換法在對原始數(shù)據(jù)進(jìn)行轉(zhuǎn)換處理之后，往往會(huì)導(dǎo)致一些數(shù)據(jù)失真，并不能保證處理后的數(shù)據(jù)呈正態(tài)分布[8]。劉曾美等[9]運(yùn)用Copula函數(shù)研究了中山市坦洲鎮(zhèn)澇區(qū)的暴雨和外江潮水位遭遇的聯(lián)合分布。任錦亮等[10]采用AMH Copula函數(shù)研究了上海市臺(tái)風(fēng)降雨和潮位遭遇組合概率。

本文采用Copula函數(shù)建立區(qū)域暴雨和潮水位2變量之間的聯(lián)合分布函數(shù)，用概率分布來描述2者遭遇的機(jī)率，并基于聯(lián)合分布研究雨潮遭遇組合重現(xiàn)期的計(jì)算方法，通過實(shí)例探討區(qū)域暴雨和潮水位遭遇的同現(xiàn)概率和內(nèi)澇發(fā)生的概率及重現(xiàn)期。

1 Copula理論及降雨-潮水位聯(lián)合分布模型

1.1 Copula理論簡介

Sklar定理：設(shè)x、y為連續(xù)的隨機(jī)變量，其邊緣分布函數(shù)分別為Fx和Fy，F(xiàn)(x,y)為變量x和y的聯(lián)合分布函數(shù)，如果Fx和Fy連續(xù)，則存在唯一的函數(shù)Cθ(u,v)使得：

F(x,y)=Cθ[Fx(x),Fy(y)],?x,y

(1)

式中：Cθ(u,v)為Copula函數(shù)；θ為待定參數(shù)。

降雨h和相應(yīng)的潮位z均為連續(xù)隨機(jī)變量，F(xiàn)(h，z)為降雨和潮位聯(lián)合分布函數(shù)，因此存在唯一的Copula函數(shù)使得F(h，z)=Cθ[Fh(h),Fz(z)]成立，F(xiàn)h(h)、Fz(z)分別為降雨、潮位的連續(xù)邊緣分布函數(shù)。

Copula函數(shù)是用來描述隨機(jī)變量間相關(guān)性結(jié)構(gòu)的函數(shù)，應(yīng)首先對隨機(jī)變量的相關(guān)性進(jìn)行度量[11]。目前常用于度量水文變量間相關(guān)性的指標(biāo)為皮爾遜線性相關(guān)系數(shù)和Kendall秩相關(guān)系數(shù)[3]。Kendall秩相關(guān)系數(shù)的計(jì)算公式為：

(2)

式中：sign(·)是符號(hào)函數(shù)，取值如下：

(3)

1.2 Copula函數(shù)選用

Copula函數(shù)中，采用阿基米德族Copula函數(shù)描述水文變量的研究越來越多，并已經(jīng)初具理論和應(yīng)用基礎(chǔ)[4,12]。變量聯(lián)合分布函數(shù)的性質(zhì)和與之對應(yīng)的Copula函數(shù)密切相關(guān)，因此，要根據(jù)變量h和z的相關(guān)性程度合理選擇阿基米德族中的Copula函數(shù)。

降雨—潮水位聯(lián)合分布函數(shù)實(shí)質(zhì)上是將一個(gè)二元聯(lián)合分布分為2個(gè)獨(dú)立的部分來分別處理：變量間的相關(guān)性結(jié)構(gòu)和變量的邊緣分布，其中相關(guān)性結(jié)構(gòu)用Copula函數(shù)來描述，邊緣分布采用皮爾遜Ⅲ型曲線來描述。

內(nèi)港站多年平均最高潮位為1.73 m，年最大24 h降雨遭遇最高潮位均值為0.91 m，相差較大，說明年最大24 h降雨基本與外江年最高潮位不遭遇；根據(jù)1960-2018年59 a降雨與潮位遭遇數(shù)據(jù)，僅有1 a相應(yīng)潮位為年最高潮位，年最大24 h降雨量與相應(yīng)最高潮位水位2變量的Kendall秩相關(guān)系數(shù)τ=0.004 1，說明這2變量相關(guān)性較弱，因此不適合用Gumbel Copula函數(shù)和Clayton Copula函數(shù)來描述它們之間的相關(guān)結(jié)構(gòu)，可以選用AMH Copula函數(shù)或Frank Copula函數(shù)[3,13]來描述它們之間的相關(guān)結(jié)構(gòu)。

如用AMH Copula函數(shù)描述降雨h和相應(yīng)潮位z之間相關(guān)結(jié)構(gòu)，則聯(lián)合分布函數(shù)F(h，z)為：

(4)

τ與θ的關(guān)系如下：

(5)

如用Frank Copula函數(shù)描述降雨h和相應(yīng)潮位z之間相關(guān)結(jié)構(gòu)，則聯(lián)合分布函數(shù)F(h，z)為：

(6)

τ與θ的關(guān)系如下：

(7)

2 降雨與潮位遭遇風(fēng)險(xiǎn)分析

降雨與潮位聯(lián)合密度函數(shù)反映區(qū)間暴雨和潮位遭遇的機(jī)率，2者遭遇存在的風(fēng)險(xiǎn)一般用風(fēng)險(xiǎn)率和重現(xiàn)期來表示。風(fēng)險(xiǎn)率指一定條件下，非期望事件發(fā)生的頻率，對于單變量來說，風(fēng)險(xiǎn)率為超值概率，因此對于降雨h和潮位z的風(fēng)險(xiǎn)率P分別表示如下：

P(h)=P(h≥H)=1-F(h

(8)

P(z)=P(z≥Z)=1-F(z

(9)

式中：H為設(shè)計(jì)暴雨；Z為排澇工況設(shè)計(jì)潮位。

超過某一特定值的重現(xiàn)期為風(fēng)險(xiǎn)率的倒數(shù)，因此對于降雨h和潮位z的重現(xiàn)期T分別表示如下：

(10)

(11)

對于中珠聯(lián)圍降雨和潮位遭遇關(guān)系有關(guān)的排澇問題，是指發(fā)生內(nèi)澇風(fēng)險(xiǎn)可能造成對人們生命財(cái)產(chǎn)、社會(huì)公共交通等造成危害或構(gòu)成不利影響。對于珠海市中珠聯(lián)圍發(fā)生內(nèi)澇的風(fēng)險(xiǎn)來說，可采用降雨h與外海潮位z的函數(shù)表示。在中珠聯(lián)圍排澇設(shè)計(jì)暴雨和相應(yīng)外海潮位條件下，以下3種情況不會(huì)發(fā)生內(nèi)澇風(fēng)險(xiǎn):①當(dāng)降雨h≤H(設(shè)計(jì)暴雨)和潮水位z≤Z(排澇工況設(shè)計(jì)潮位)時(shí)，不會(huì)發(fā)生內(nèi)澇；②當(dāng)潮水位z>Z(排澇工況設(shè)計(jì)潮位)且不超過海堤設(shè)計(jì)潮位值Z2時(shí)，因中珠聯(lián)圍內(nèi)河網(wǎng)密集高，調(diào)蓄空間大，只要降雨h≤H1(與河涌調(diào)蓄容積相關(guān))，即使中珠聯(lián)圍不向外排水的情況下，亦不會(huì)發(fā)生內(nèi)澇；③當(dāng)潮水位z≤Z1(多年平均低潮位)時(shí)，因中珠聯(lián)圍內(nèi)河道排水和泄水閘排水能力強(qiáng)，只要降雨h≤H2(河涌過流及水閘泄流能力條件下相應(yīng)降雨)，亦不會(huì)發(fā)生內(nèi)澇。因此發(fā)生內(nèi)澇的風(fēng)險(xiǎn)區(qū)域見圖1。

圖1 內(nèi)澇風(fēng)險(xiǎn)區(qū)示意圖Fig.1 Drawing of waterlogging risk zone

于是，中珠聯(lián)圍內(nèi)澇風(fēng)險(xiǎn)率采用下式表示：

P(h,z)=1-F(H,Z)-F(H1,Z2)-F(H2,Z1)+

F(H,Z1)+F(H1,Z)

(12)

如果采用水利標(biāo)準(zhǔn)重現(xiàn)期來表示內(nèi)澇風(fēng)險(xiǎn)，則內(nèi)澇重現(xiàn)期T(h，z)可采用下式表示：

T(h,z)=

(13)

往往澇區(qū)整治制定內(nèi)澇標(biāo)準(zhǔn)時(shí)比較關(guān)心降雨和潮位均超過某一量級的概率，即P(h≥H,z≥Z)，稱為同現(xiàn)概率，其倒數(shù)為同現(xiàn)重現(xiàn)期，具體表達(dá)式如下：

P(h≥H,z≥Z)=1-Fh(H)-Fz(Z)+F(H,Z)

(14)

(15)

3 實(shí)例分析

珠海市中珠聯(lián)圍位于珠江口磨刀門水道東岸，鄰近澳門。中珠聯(lián)圍主要通過前山河排渠灣仔水道，灣仔水道內(nèi)設(shè)有內(nèi)港潮位站，臨近有神灣雨量站，因此，采用神灣雨量站1960-2018年的歷年最大24 h量以及內(nèi)港潮位站的相應(yīng)潮位來分析確定最大24 h雨和其相應(yīng)潮位的聯(lián)合分布，計(jì)算步驟如下：

(1)確定單變量降雨、相應(yīng)潮位及年最高潮位的P-Ⅲ型曲線。

(2)利用式(2)、式(3)計(jì)算Kendall秩相關(guān)系數(shù)τ。

(3)利用式(16)、式(17)統(tǒng)計(jì)樣本經(jīng)驗(yàn)頻率，判斷其擬合效果。

(4)擇優(yōu)選取聯(lián)合分布函數(shù)，計(jì)算降雨與潮位同現(xiàn)概率及內(nèi)澇風(fēng)險(xiǎn)率。

3.1 單變量和聯(lián)合分布參數(shù)的確定

神灣站年最大24 h降雨、內(nèi)港站相應(yīng)最高潮位和年最高潮位均服從P-Ⅲ型曲線分布，本文采用矩法對降雨、相應(yīng)潮位和年最高潮位的參數(shù)進(jìn)行估算，并結(jié)合適線法優(yōu)化估算參數(shù)，得到年最大24 h降雨、相應(yīng)潮位和年最高潮位的統(tǒng)計(jì)參數(shù)，見表1。

表1 P-Ⅲ型曲線參數(shù)估算成果Tab.1 Estimated results of P-Ⅲ curve parameters

根據(jù)神灣站和內(nèi)港站實(shí)測資料，利用式(2)和式(3)可計(jì)算得到最大24 h降雨和相應(yīng)潮位的Kendall秩相關(guān)系數(shù)τ=0.004 1。說明神灣站年最大24小時(shí)降雨與內(nèi)港相應(yīng)潮位相關(guān)關(guān)系很弱。利用相關(guān)性指標(biāo)法[9]，采用式(5)、式(7)分別計(jì)算得θ為0.018 4和0.036 9。

3.2 擬合效果評價(jià)

為選定最大24 h降雨與相應(yīng)潮位擬合最高的Copula函數(shù)，需對Copula函數(shù)進(jìn)行檢驗(yàn)，本文采用K-S檢驗(yàn)來評價(jià)2變量的理論聯(lián)合概率分布函數(shù)與經(jīng)驗(yàn)聯(lián)合概率分布的擬合程度。

假設(shè)第i對聯(lián)合觀測值(hi，zi)，按降序排列后對應(yīng)為(j，k)，則(hi，zi)的累積經(jīng)驗(yàn)頻率Femp(hi，zi)可由下式求得：

Femp(hi,zi)=P(h≤hi,z≤zi)=P[h≤h(j),z≤z(k)]=

(16)

式中：n為樣本數(shù)量；j為hj的排序號(hào)；k為zj的排序號(hào)；nm,l為符號(hào)函數(shù)，取值如下：

(17)

K-S檢驗(yàn)統(tǒng)計(jì)量D計(jì)算公式如下：

(18)

除K-S檢驗(yàn)外，還需根據(jù)評價(jià)指標(biāo)選取最適合的Copula函數(shù)。本文選取理論聯(lián)合分布概率與經(jīng)驗(yàn)聯(lián)合分布概率適線、離差平方和最小準(zhǔn)則(ΟLS)以及AIC信息準(zhǔn)則來評價(jià)，并選取ΟLS和AIC最小的Copula作為聯(lián)合概率分布函數(shù)[9]。其計(jì)算公式如下：

(19)

AIC=nln (MSE)+2K

(20)

式中：K為模型參數(shù)的個(gè)數(shù)，取59。

(21)

神灣站最大24 h降雨與內(nèi)港相應(yīng)潮位其經(jīng)驗(yàn)聯(lián)合分布值與理論聯(lián)合分布值散點(diǎn)均落在45°對角線附近，見圖2。其誤差來源于經(jīng)驗(yàn)聯(lián)合分布由實(shí)測數(shù)據(jù)統(tǒng)計(jì)，統(tǒng)計(jì)過程中小于一定數(shù)值的計(jì)為1，否則計(jì)為0，而沒有中間數(shù)值，導(dǎo)致經(jīng)驗(yàn)頻率與理論頻率具有一定的誤差。另外Copula函數(shù)僅能近似表達(dá)雨潮相關(guān)關(guān)系用，其誤差主要來源于以上2個(gè)方面。雖然經(jīng)驗(yàn)分布與理論分布存在誤差，但AMH Copula函數(shù)擬合相關(guān)系數(shù)達(dá)0.955，F(xiàn)rank Copula函數(shù)擬合相關(guān)系數(shù)達(dá)0.950，擬合效果良好，雨潮聯(lián)合分布可采用Copula函數(shù)來表示。離差平方和最小準(zhǔn)則(ΟLS)、AIC信息準(zhǔn)則指標(biāo)和K-S檢驗(yàn)統(tǒng)計(jì)結(jié)果見表2。

圖2 降雨與潮位聯(lián)合理論分布與經(jīng)驗(yàn)分布散點(diǎn)圖Fig.2 Scatter plots of theoretical and empirical distributions of rainfall and tide levels

表2 擬合檢驗(yàn)統(tǒng)計(jì)量Tab.2 Fitting test statistical scale

由表2可知，2種Copula函數(shù)均能通過K-S檢驗(yàn)，且與聯(lián)合經(jīng)驗(yàn)頻率擬合很好，說明所采用的Copula理論聯(lián)合分布函數(shù)是合理的，根據(jù)擬合度檢驗(yàn)統(tǒng)計(jì)量結(jié)果，本文選用擬合效果更優(yōu)的AMH Copula函數(shù)來描述神灣站最大24 h降雨與內(nèi)港相應(yīng)潮位的聯(lián)合概率分布，其聯(lián)合分布函數(shù)如下：

(22)

3.3 組合遭遇頻率和內(nèi)澇風(fēng)險(xiǎn)分析

根據(jù)神灣雨量站和內(nèi)港潮位站的最大24 h暴雨、相應(yīng)潮位的邊緣分布，計(jì)算不同重現(xiàn)期下的降雨和潮位組合，對于給定的相應(yīng)潮位，在年最高潮位頻率曲線中對應(yīng)一定的頻率，從而計(jì)算得到相應(yīng)潮位的年重現(xiàn)期。

由式(22)可求得降雨和潮位的聯(lián)合分布概率，由式(12)、式(13)可求得中珠聯(lián)圍內(nèi)澇的風(fēng)險(xiǎn)率及重現(xiàn)期，由式(14)、(15)求得設(shè)計(jì)雨潮遭遇組合同現(xiàn)概率及重現(xiàn)期。表3給出了不同頻率降雨、潮位組合的概率及中珠聯(lián)圍內(nèi)澇的風(fēng)險(xiǎn)率和重現(xiàn)期。

由表3可知，中珠聯(lián)圍發(fā)生20 a一遇暴雨與2 a年一遇、10 a一遇、5 a一遇、多年平均最高潮位遭遇同現(xiàn)概率分別為0.02%、0.04%、0.06%和0.10%；10 a一遇暴雨與20 a一遇、10 a一遇、5 a一遇、多年平均最高潮位遭遇同現(xiàn)概率分別為0.05%、0.08%、0.13%和0.21%，遭遇概率較低，說明中珠聯(lián)圍發(fā)生10 a一遇以上降雨時(shí)，遭遇多年平均最高潮位以上潮位的幾率非常??；根據(jù)《珠海市流域綜合規(guī)劃修編》，中珠聯(lián)圍內(nèi)澇防治標(biāo)準(zhǔn)采用20 a一遇最大24 h設(shè)計(jì)暴雨不成災(zāi)，外江遭遇5 a一遇最高潮位，由雨潮遭遇同現(xiàn)頻率及重現(xiàn)期分析可知，20 a一遇暴雨遭遇5 a一遇高潮位組合的同現(xiàn)頻率為0.06%，重現(xiàn)期為1 548 a，遠(yuǎn)高于內(nèi)澇防治標(biāo)準(zhǔn)的20 a一遇。

然而，中珠聯(lián)圍實(shí)際內(nèi)澇的情況并非暴雨、潮位同時(shí)超過設(shè)計(jì)值時(shí)才發(fā)生內(nèi)澇，因此中珠聯(lián)圍發(fā)生內(nèi)澇的概率要遠(yuǎn)大于雨潮組合同現(xiàn)概率。根據(jù)表3計(jì)算結(jié)果，當(dāng)中珠聯(lián)圍發(fā)生20年一遇暴雨遭遇外海20 a一遇、10 a一遇、5 a一遇、多年平均最高潮位時(shí)，發(fā)生內(nèi)澇的概率分別為5.34%、5.61%、6.09%和6.82%；當(dāng)中珠聯(lián)圍發(fā)生10 a一遇暴雨遭遇外海20 a一遇、10 a一遇、5 a一遇、多年平均最高潮位時(shí)，發(fā)生內(nèi)澇的概率分別為10.17%、10.43%、10.88%和11.57%。發(fā)生內(nèi)澇的概率均高于暴雨頻率，重現(xiàn)期低于暴雨設(shè)計(jì)重現(xiàn)期。

表3 遭遇組合概率和內(nèi)澇風(fēng)險(xiǎn)計(jì)算成果Tab.3 Calculation results of encounter combination probability and waterlogging risk

4 結(jié) 論

(1)采用P-Ⅲ型分布曲線擬合珠海市中珠聯(lián)圍最大24 h降雨、相應(yīng)潮位和年最高潮位，采用矩法估算曲線參數(shù)，結(jié)合適線法優(yōu)化估算參數(shù)，得到不同重現(xiàn)期的設(shè)計(jì)值，并運(yùn)用Copula函數(shù)構(gòu)建最大24 h降雨和相應(yīng)潮位的聯(lián)合分布模型，通過理論分布值與K-S檢驗(yàn)，表明AMH Copula函數(shù)能較好地?cái)M合其聯(lián)合分布。

(2)考慮中珠聯(lián)圍內(nèi)河涌及排水閘泵條件，構(gòu)建內(nèi)澇風(fēng)險(xiǎn)率計(jì)算模型及雨潮同現(xiàn)頻率計(jì)算模型，最后給出不同雨潮組合的內(nèi)澇風(fēng)險(xiǎn)率及同現(xiàn)頻率。

(3)根據(jù)計(jì)算結(jié)果，中珠聯(lián)圍內(nèi)澇風(fēng)險(xiǎn)率大于設(shè)計(jì)暴雨重現(xiàn)期，遠(yuǎn)大于雨潮同現(xiàn)概率；分析中珠聯(lián)圍內(nèi)澇風(fēng)險(xiǎn)率為治澇工程的暴雨和潮水位的遭遇組合確定提供參考，明確發(fā)生內(nèi)澇的風(fēng)險(xiǎn)率。中珠聯(lián)圍20 a一遇最大24 h設(shè)計(jì)暴雨遭遇外海5 a一遇高潮位設(shè)計(jì)工況條件下，內(nèi)澇發(fā)生的概率為6.09%，重現(xiàn)期為16.04 a。

(4)根據(jù)內(nèi)澇發(fā)生概率可知，中珠聯(lián)圍內(nèi)澇發(fā)生的重現(xiàn)期低于暴雨設(shè)計(jì)重現(xiàn)期，中珠聯(lián)圍內(nèi)澇整治工程實(shí)施過程中，應(yīng)適當(dāng)增加工程規(guī)模，使內(nèi)澇發(fā)生重現(xiàn)期與設(shè)計(jì)暴雨重現(xiàn)期相適應(yīng)。

(5)本文構(gòu)建的降雨和相應(yīng)潮位的二維聯(lián)合分布，主要考慮暴雨和外江潮水位頂托可能形成內(nèi)澇的組合風(fēng)險(xiǎn)，若考慮防潮風(fēng)險(xiǎn)可分析潮位與相應(yīng)降雨的聯(lián)合分布函數(shù)來考慮，但對于一個(gè)完整的防洪封閉圈而言，不能簡單采用Copula函數(shù)分析防洪、防潮、防澇綜合風(fēng)險(xiǎn)概率及重現(xiàn)期，需進(jìn)一步研究。