裴昊川,倪 倩,甘富萬,3

(1.廣西大學(xué) 土木建筑工程學(xué)院,南寧 530004;2.廣西珠江委員會南寧勘測設(shè)計(jì)院,南寧 530007;3.廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室,南寧 530004)

近些年來,Copula函數(shù)在解決洪水、降雨、干旱等多變量水文聯(lián)合分布問題上應(yīng)用廣泛。潘國勇等[1]采用Copula函數(shù)構(gòu)建了武澄錫虞區(qū)的雨、潮聯(lián)合風(fēng)險(xiǎn)模型,并分析該地區(qū)雨、潮組合風(fēng)險(xiǎn)率,為武澄錫虞區(qū)的防洪治澇提供了新的科學(xué)依據(jù);許翰卿等[2]應(yīng)用二元Copula函數(shù)構(gòu)建沿海地區(qū)極端降水和風(fēng)暴增水兩變量聯(lián)合風(fēng)險(xiǎn)模型,為防汛工程抵御復(fù)合洪澇災(zāi)害的設(shè)計(jì)標(biāo)準(zhǔn)提供了新的設(shè)計(jì)方法;藺文慧等[3]采用混合Copula函數(shù)構(gòu)建馬達(dá)瓦斯克流域洪水歷時(shí)、洪量及洪峰流量混合模型,對比單一Copula函數(shù)擬合效果更優(yōu),能更加全面表征變量間相關(guān)關(guān)系。本文以北海市鐵山港區(qū)入海明渠為研究對象,采用北海站雨量資料和石頭埠潮位站潮位資料進(jìn)行兩變量分析,首先分析雨、潮聯(lián)合概率分布特征,其次引用了一種考慮雨、潮與調(diào)洪耦合的防洪重現(xiàn)期,并將其與常見的OR、AND重現(xiàn)期進(jìn)行對比分析,最后基于蒙特卡羅法分析抽樣不確定性對于雨、潮共同作用下明渠設(shè)計(jì)水深的影響,此研究為雨、潮作用下的入海明渠的設(shè)計(jì)水深確定提供了一條新的設(shè)計(jì)思路。

1 研究方法

1.1 Copula函數(shù)

Copula函數(shù)可以將多個具有不同形式邊緣分布的隨機(jī)變量聯(lián)結(jié)起來構(gòu)造聯(lián)合分布。令雨量X的概率分布為FX(x),潮位Y的概率分布為FY(y)。依據(jù)Sklar定理,X和Y的聯(lián)合分布函數(shù)可以用Copula函數(shù)C表示成f(x,y)=C(u,v)[4],其中u=FX(x),v=FY(y)。

f(x,y)相應(yīng)的聯(lián)合概率密度函數(shù)可表示為

f(x,y)=C(u,v)fX(x)fY(y)

(1)

式中：C(u,v)=?C(u,c)/?u?v為Copula函數(shù)的密度函數(shù);fX、fY分別為X和Y的概率密度函數(shù)。

在現(xiàn)階段水文氣象領(lǐng)域中,通常采用Archimedean Copula函數(shù)族函數(shù)構(gòu)建多變量水文聯(lián)合分布模型[5-6]。4種常用的Archimedean Copula函數(shù)見表1。

表1 4種常用的Archimedean Copula函數(shù)Tab.1 Four commonly used Archimedean Copula functions

1.2 OR、AND重現(xiàn)期和防洪重現(xiàn)期

AND重現(xiàn)期指所研究區(qū)域降雨量和潮位同時(shí)超過某一特定閾值。

OR重現(xiàn)期指所研究區(qū)域降雨量或潮位其中一個超過某一特定閾值。

定義如下

TOR=1/P(X>x∪Y>y)=1/[1-C(u,v)]

(2)

TAND=1/P(X>x∩Y>y)=1/{1-u-v+C(u,v)}

(3)

根據(jù)OR和AND重現(xiàn)期定義,在同一重現(xiàn)期水平下,有無數(shù)種設(shè)計(jì)雨、潮組合;這表明設(shè)計(jì)雨、潮重現(xiàn)期和設(shè)計(jì)雨、潮組合之間并沒有一一對應(yīng)的關(guān)系,這將給水利工程的運(yùn)行帶來了巨大風(fēng)險(xiǎn)。因此,針對類似雨、潮組合等二變量水文事件的設(shè)計(jì)重現(xiàn)期問題,甘富萬等[7]在研究干、支流洪水二變量共同作用下的防洪設(shè)計(jì)時(shí),提出了防洪重現(xiàn)期的概念,即由洪水(雨、潮)等作用所確定的防洪參數(shù)F超過工程防洪(防潮)設(shè)計(jì)的給定值f的平均間隔時(shí)間,在本文中可以表示為

Tf=1/[F>f設(shè)=f(x,y)]

(4)

式中：F為雨、潮設(shè)計(jì)的防洪參數(shù);x和y分別為雨量和潮位。

黃強(qiáng)等[8]認(rèn)為最可能組合模式可以有效節(jié)省工程費(fèi)用,在同時(shí)滿足工程防洪標(biāo)準(zhǔn)的條件下,其模式代表雨量、潮位聯(lián)合設(shè)計(jì)值取最大值的情況,其表達(dá)式如下

(um,vm)=argmaxf(u,v)

(5)

f(u,v)=C(u,v)f(x)f(y)

(6)

式中：C(u,v)為雨量、潮位聯(lián)合分布的概率密度函數(shù);f(x)和f(y)為邊緣分布的概率密度函數(shù)。

1.3 不確定性度量指標(biāo)

本文利用尹家波等[9]提出的平均橫向偏移程度DX、平均縱向偏移程度DY、置信區(qū)域面積S和平均歐氏距離d作為兩變量估計(jì)不確定性的度量指標(biāo)。DX、DY分別用于度量雨量、潮位與實(shí)測樣本系列推求的設(shè)計(jì)值在一維空間的估計(jì)偏差;S和d用于度量模擬設(shè)計(jì)值點(diǎn)據(jù)與實(shí)測樣本系列設(shè)計(jì)值的空間距離。計(jì)算式如下

(7)

(8)

(9)

1.4 基于防洪重現(xiàn)期的設(shè)計(jì)水深計(jì)算

(1)分別建立雨、潮組合的邊緣分布函數(shù)u=F(x)和v=F(y),并依此構(gòu)建并選取雨、潮組合的最優(yōu)Copula函數(shù)模型C(u,v);

(2)產(chǎn)生樣本容量為n的范圍在(0,1)的隨機(jī)數(shù)n1和n2,設(shè)n1=u,由條件概率公式求解出每一個與u相關(guān)的v,即令C(v|u)=?C(u,v)/?v=n2;

(3)根據(jù)雨、潮最優(yōu)邊緣分布函數(shù),反推每個隨機(jī)組合(u,v)所對應(yīng)的雨、潮組合的設(shè)計(jì)值(x,y),帶入調(diào)洪模型計(jì)算出每個設(shè)計(jì)斷面的設(shè)計(jì)水深;

(4)對每個斷面的n個水深值進(jìn)行排序,根據(jù)水文計(jì)算方法計(jì)算出不同設(shè)計(jì)重現(xiàn)期水平年的設(shè)計(jì)水深h設(shè)。

1.5 基于Monte Carlo Method的抽樣不確定性分析

(1)與上述(1)和(2)步驟相同,獲得樣本容量為n的二元雨、潮組合樣本;

(2)重復(fù)上述步驟N次,便可得到N組與樣本容量n相同的二元雨、潮組合樣本。對N組二元雨、潮組合樣本使用實(shí)測樣本序列的估參方法即可得到N組參數(shù)值;

(3)在每一組參數(shù)下,給定OR重現(xiàn)期TOR,利用式(5)和式(6)得到特定重現(xiàn)期下N組參數(shù)對應(yīng)的N組最可能組合聯(lián)合設(shè)計(jì)值。對于N組最可能聯(lián)合設(shè)計(jì)值,采用核密度估計(jì)方法得到給定某一顯著性水平α的(1-α)置信區(qū)域,置信區(qū)域越大,表明聯(lián)合設(shè)計(jì)值不確定性越大;

(4)在每一組參數(shù)下,給定防洪標(biāo)準(zhǔn)T,計(jì)算OR重現(xiàn)期標(biāo)準(zhǔn)下的設(shè)計(jì)水深。在OR重現(xiàn)期標(biāo)準(zhǔn)下,將步驟(3)中得到的N組最可能設(shè)計(jì)組合值代入調(diào)洪演算模型可得到N個水深設(shè)計(jì)值zOR。從而可求得重現(xiàn)期標(biāo)準(zhǔn)的單變量z點(diǎn)估計(jì)和置信區(qū)間。置信區(qū)間寬度越小,設(shè)計(jì)水深不確定性越小。

2 案例研究

鐵山港區(qū)是廣西地區(qū)面向東南亞的重要出海港口,該區(qū)在夏季易受臺風(fēng)的侵襲,臺風(fēng)帶來強(qiáng)降雨以及洪暴潮,且該區(qū)域降雨主要集中5月至9月,故通過規(guī)劃入海明渠實(shí)現(xiàn)區(qū)域內(nèi)洪澇積水自流入海,以應(yīng)對洪潮災(zāi)害。北海鐵山港工業(yè)區(qū)排水明渠起點(diǎn)在北海至鐵山港一級公路和規(guī)劃的鐵山港鐵路支線的交叉口附近,終點(diǎn)至鐵山港港口,全程自流排放,最后注入大海。線路全長約7 996.602 m,干渠設(shè)計(jì)雨水排泄流量為166.09 m3/s。本文以北海市鐵山港區(qū)入海明渠為研究對象,降雨量資料采用北海站實(shí)測連續(xù)42 a年最大降雨量序列,潮位資料采用降雨量對應(yīng)石頭埠潮位站實(shí)測年最大潮位序列。

2.1 邊緣分布及聯(lián)合分布的確定

降雨量X代表北海站實(shí)測連續(xù)42 a年最大降雨量序列,潮位Y代表降雨量對應(yīng)石頭埠潮位站實(shí)測年最大潮位序列。采用AIC和RMSE準(zhǔn)則[10]對皮爾遜III型(PE3)、Weibull分布、廣義Logistic分布(GLO)和廣義極值分布(GEV)4種雨量、潮位邊緣分布函數(shù)進(jìn)行評估,AIC和RMSE準(zhǔn)則評估數(shù)值越小,則雨、潮邊緣分布模型效果越優(yōu),由表2可知皮爾遜III型(PE3)函數(shù)構(gòu)建雨量和潮位序列邊緣分布模型效果最優(yōu)。

表2 雨、潮邊緣分布計(jì)算結(jié)果Tab.2 Calculated results of rainfall and tidal level margin distribution

根據(jù)皮爾遜III型(PE3)函數(shù)所構(gòu)建雨量和潮位序列邊緣分布模型可以得到不同設(shè)計(jì)頻率下的雨量和潮位計(jì)算結(jié)果,表3為設(shè)計(jì)頻率為 1%、2%、5%、10%及20%所對應(yīng)的雨量和潮位設(shè)計(jì)值。

表3 雨量、潮位資料PE3分布擬合計(jì)算結(jié)果Tab.3 Rainfall and tide level data PE3 distribution fitting calculation

采用Archimedean Copula函數(shù)族中的G-H、Clayton、Frank和AMH函數(shù)構(gòu)建降雨量X和潮位Y聯(lián)合分布模型,Copula函數(shù)參數(shù)θ由Kendall秩相關(guān)性系數(shù)法求得,采用BIC準(zhǔn)則、RMSE準(zhǔn)則和AIC準(zhǔn)則可以有效篩選出擬合結(jié)果最優(yōu)的Copula函數(shù)。BIC準(zhǔn)則、RMSE準(zhǔn)則和AIC準(zhǔn)則數(shù)值越小,則雨、潮聯(lián)合分布模型效果越優(yōu),從表4可以看出,擬合結(jié)果最優(yōu)Copula函數(shù)為G-H Copula函數(shù)。

表4 4種候選Copula函數(shù)參數(shù)估計(jì)與擬合檢驗(yàn)結(jié)果Tab.4 Parameter estimation and fitting test results of four candidate Copula functions

因此,G-H Copula函數(shù)構(gòu)建北海站降雨量對應(yīng)石頭埠潮位聯(lián)合分布模型如下

CGH=exp{-[(-lnu)1.277 5+(-lnv)1.277 5]1/1.277 5}

(10)

綜上所述,雨、潮組合的聯(lián)合分布模型如圖1。

圖1 雨、潮聯(lián)合概率分布模型Fig.1 Joint probability distribution model of rainfall and tidal level

2.2 不同重現(xiàn)期特征分析

如圖2中T分別為100 a、50 a、20 a、10 a和5 a的OR和AND重現(xiàn)期等值線所示,重現(xiàn)期等值線與45°線呈對稱分布;這表明雨量和潮位可能對工程的設(shè)計(jì)水平計(jì)算具有相同的影響。然而,在實(shí)際工程中,雨量和潮位對工程的設(shè)計(jì)水平計(jì)算影響卻不盡相同。因此,雖然圖2在一定程度上顯示了雨量與潮位的相關(guān)性,但未能考慮雨、潮組合與工程調(diào)洪的耦合關(guān)系。

2-a AND重現(xiàn)期等值線圖 2-b OR重現(xiàn)期等值線圖圖2 重現(xiàn)期等值線圖Fig.2 Contour calculation results of design periods

如圖3的防洪重現(xiàn)期等值線所示,防洪重現(xiàn)期等值線上尾部分的斜率絕對值小于1,說明當(dāng)雨量較大,潮位較低時(shí),此時(shí)對明渠而言雨量起主要作用;隨著雨量逐漸減小,潮位不斷升高,防洪重現(xiàn)期等值線斜率絕對值逐漸大于1,說明此時(shí)潮位對調(diào)洪起主要作用,雨量的影響相對較小。

圖3 防洪重現(xiàn)期等值線圖Fig.3 Contour calculation results of flood control return period

圖4描述了OR、AND和防洪重現(xiàn)期為20 a一遇的設(shè)計(jì)情況。圖4中防洪重現(xiàn)期選取的是斷面K7+340,按照防洪重現(xiàn)期計(jì)算其設(shè)計(jì)水深為4.02 m。通過對比發(fā)現(xiàn),防洪重現(xiàn)期等值線處于OR、AND重現(xiàn)期等值線之間。常用的明渠水面線計(jì)算方法有均勻流法[11]、非均勻流法[12]和數(shù)值模擬方法[13-14]等,本文采用明渠非均勻流法,圖5將防洪重現(xiàn)期20 a一遇水面線與OR、AND重現(xiàn)期20 a一遇水面線進(jìn)行對比發(fā)現(xiàn),防洪重現(xiàn)期20 a一遇設(shè)計(jì)水面線也處于OR、AND重現(xiàn)期之間,其中OR重現(xiàn)期對應(yīng)的水面線比防洪重現(xiàn)期對應(yīng)的水面線高約0.20 m,而AND重現(xiàn)期對應(yīng)的水面線則要低約0.22 m,這表明防洪重現(xiàn)期可以避免出現(xiàn)OR、AND重現(xiàn)期設(shè)計(jì)結(jié)果偏高或偏低的問題,能有效提高工程防洪、防潮設(shè)計(jì)水平。

圖4 三種設(shè)計(jì)重現(xiàn)期20 a一遇等值線圖Fig.4 Contour calculation results of the recurrence period of 20 years of three designs

圖5 三種設(shè)計(jì)重現(xiàn)期20 a一遇設(shè)計(jì)水面線計(jì)算成果圖Fig.5 Calculation results of the design water surface line with the recurrence period of 20 years of three designs

2.3 雨、潮不確定性分析

SERINALDIF[15]采用小于樣本容量長度的OR重現(xiàn)期分析兩變量設(shè)計(jì)值不確定性問題。因此本文以小于實(shí)測樣本序列(42 a)的10 a、20 a OR重現(xiàn)期作為分析對象,基于實(shí)測數(shù)據(jù)建立的聯(lián)合分布模型為總體,分別設(shè)置樣本容量n=42、n=100、n=200分析不同重現(xiàn)期下序列長度對聯(lián)合設(shè)計(jì)值的影響。

聯(lián)合分布的不確定性是由于參數(shù)的不確定性引起的。表5給出基于蒙特卡羅法不同樣本容量條件下重現(xiàn)期為20 a,聯(lián)合分布參數(shù)在95%水平下的置信區(qū)間?？梢钥闯?參數(shù)的變幅隨著樣本容量的增加而不斷減小。

表5 不同樣本容量聯(lián)合分布參數(shù)95%置信區(qū)間Tab.5 95% Confidence interval of joint distribution parameters of different sample size

如圖6所示,同一OR重現(xiàn)期水平條件下,樣本容量的增加導(dǎo)致最可能組合聯(lián)合設(shè)計(jì)值95%置信區(qū)間逐漸減小。在同一樣本容量n下,聯(lián)合設(shè)計(jì)值二元置信區(qū)間隨著OR重現(xiàn)期水平的增加而增大。此外,在T=20和n=42時(shí),聯(lián)合設(shè)計(jì)值最可能組合主要分布在T=5到T=100曲線之間。相反,在T=10和n=200時(shí),聯(lián)合設(shè)計(jì)值最可能組合主要集中在T=5和T=20曲線之間。這些結(jié)果表明,當(dāng)樣本容量n較小且設(shè)計(jì)重現(xiàn)期水平較大時(shí),雨、潮聯(lián)合設(shè)計(jì)值的不確定性也更大。

6-a n=42,TOR=10 6-b n=42,TOR=20 6-c n=100,TOR=10

采用本文1.3節(jié)所介紹的4個不確定性評價(jià)指標(biāo)(DX、DY、S、d)對雨、潮組合不確定性進(jìn)行了評估,由表6可知,DX、DY隨著樣本數(shù)量n增加而減小;當(dāng)T=10 a時(shí),d和95%置信區(qū)間的面積S分別減小38.9%、29.1%、38.3%和27.1%,在T=20 a時(shí),分別減少59.9%、49.0%、54.5%和52.8%。

表6 OR重現(xiàn)期設(shè)計(jì)值不確定性評價(jià)指標(biāo)計(jì)算結(jié)果Tab.6 Calculating results of uncertainty evaluation index of OR design period

2.4 設(shè)計(jì)水深不確定性分析

對于明渠的防洪安全而言,最主要的因素是設(shè)計(jì)水深,因此聯(lián)合設(shè)計(jì)雨、潮組合的不確定性問題就轉(zhuǎn)化為單變量的不確定性問題,即設(shè)計(jì)水深的不確定性,根據(jù)調(diào)洪規(guī)則可分別得到每一組參數(shù)下OR重現(xiàn)期最可能組合對應(yīng)的N個設(shè)計(jì)水深以及防洪重現(xiàn)期標(biāo)準(zhǔn)下的設(shè)計(jì)水深。表7給出了不同重現(xiàn)期標(biāo)準(zhǔn)下明渠K7+340處設(shè)計(jì)水深的計(jì)算結(jié)果,可以看出在不同樣本容量下,以防洪重現(xiàn)期為設(shè)計(jì)標(biāo)準(zhǔn)推求得到的20 a一遇設(shè)計(jì)水深的95%置信區(qū)間寬度和標(biāo)準(zhǔn)差均小于OR重現(xiàn)期,例如在雨、潮組合序列長度為實(shí)測序列長度(42 a)時(shí),防洪重現(xiàn)期比OR重現(xiàn)期標(biāo)準(zhǔn)差偏小14.3%,區(qū)間寬度偏小21.4%,這說明了以防洪重現(xiàn)期為設(shè)計(jì)標(biāo)準(zhǔn)可減少設(shè)計(jì)水深的不確定性。

表7 不同重現(xiàn)期標(biāo)準(zhǔn)下明渠設(shè)計(jì)水深估計(jì)結(jié)果Tab.7 Estimation results of open channel design water depth under different return period standards

從表7還可以看出,當(dāng)樣本容量小于100時(shí),兩種重現(xiàn)期標(biāo)準(zhǔn)下的設(shè)計(jì)水深95%置信區(qū)間寬度均超過了10 mm,標(biāo)準(zhǔn)差均超過了0.01,也就是說當(dāng)雨、潮樣本序列長度較短時(shí),設(shè)計(jì)水深估計(jì)結(jié)果在數(shù)值上波動性較大,會給明渠的設(shè)計(jì)水深帶來較大的不確定性。因此,在實(shí)際工程中,有必要延長雨、潮序列使得設(shè)計(jì)結(jié)果更加可靠。

3 結(jié)論

以雨、潮聯(lián)合作用下的鐵山港入海明渠為研究對象,采用Copula函數(shù)構(gòu)建兩種重現(xiàn)期風(fēng)險(xiǎn)模型分析雨、潮遭遇及組合規(guī)律,并針對實(shí)際工程所采用的設(shè)計(jì)方法對雨、潮變量間的相關(guān)性考慮不足情況,提出了一種考慮雨、潮與工程調(diào)洪耦合關(guān)系的防洪重現(xiàn)期,并將該重現(xiàn)期與較為常見的兩種二變量重現(xiàn)期進(jìn)行對比分析,最后基于蒙特卡羅法分析抽樣不確定性對于雨、潮共同作用下明渠設(shè)計(jì)水深的影響,得出以下結(jié)論：(1)與OR、AND重現(xiàn)期相比,防洪重現(xiàn)期可以有效避免最可能設(shè)計(jì)組合值設(shè)計(jì)標(biāo)準(zhǔn)偏高或偏低的情況;(2)在雨、潮聯(lián)合作用下,最可能組合聯(lián)合設(shè)計(jì)值95%置信區(qū)間隨著重現(xiàn)期水平的遞增而逐漸增加,隨著樣本容量的增加而逐漸減小;(3)明渠水深的不確定性和重現(xiàn)期的選取有關(guān),防洪重現(xiàn)期與傳統(tǒng)的OR重現(xiàn)期相比更能有效減少明渠設(shè)計(jì)水深的不確定性。