廣東西江北江洪水聯(lián)合概率分布研究*

陳子燊，劉曾美，路劍飛

(1. 中山大學(xué)水資源與環(huán)境系，廣東 廣州510275；2. 華南理工大學(xué)水利水電工程系，廣東 廣州510640)

1 研究背景

多變量極端事件之間的相關(guān)性普遍存在于水文氣象領(lǐng)域。然而，至今對于極端水文氣象事件的概率分布研究多數(shù)還是基于單變量的統(tǒng)計(jì)分析。一些研究嘗試?yán)枚嘧兞空龖B(tài)聯(lián)合分布，但由于要求邊緣分布必須是正態(tài)或?qū)?shù)正態(tài)分布，因此限制了模型使用和選擇范圍。近些年來借助于Copula函數(shù)，多變量聯(lián)合分布已經(jīng)在水文氣象極值事件分析中得到了較深入的研究[1-9]，為更精確地推算極端事件的重現(xiàn)水平和控制工程風(fēng)險(xiǎn)提供了一個(gè)重要的理論基礎(chǔ)和技術(shù)手段。

流入廣東省的西江發(fā)源于云南省霑益縣馬雄山，至廣西梧州會桂江后始稱西江，東流至三水市思賢滘與北江相通。北江發(fā)源于江西省信豐縣石碣大茅坑，流入韶關(guān)市與武江匯合后始稱北江，此后向南流經(jīng)清遠(yuǎn)市，經(jīng)思賢滘與西江干流連通后進(jìn)入珠江三角洲網(wǎng)河區(qū)。位于珠江三角洲西北部頂點(diǎn)思賢滘兩側(cè)的馬口站和三水站屬于國家重點(diǎn)水文觀測站，分別監(jiān)測進(jìn)入西北江三角洲網(wǎng)河區(qū)的流量過程。20世紀(jì)80年代以來，珠江三角洲經(jīng)濟(jì)和城市化進(jìn)程高速發(fā)展，三角洲出現(xiàn)了大規(guī)模無序的河道采沙，河道發(fā)生了大規(guī)模的下切，而三水測站所處的河道斷面下切幅度遠(yuǎn)遠(yuǎn)大于同期馬口測站所在的河道斷面下切幅度。由于思賢滘的連通作用，造成三水?dāng)嗝娴乃糠峙浔炔粩嘣龃骩10 ]。西北江洪水受流域特征制約，兩江洪水特性差異明顯。每年洪水發(fā)生時(shí)間通常北江早于西江，兩江同時(shí)發(fā)洪的情況也時(shí)有發(fā)生[11]。通過馬口、三水的兩江洪水直接威脅珠江西北江三角洲的防洪安全，深入分析兩江洪水遭遇的聯(lián)合概率分布特征對于防御西、北江洪水拱衛(wèi)廣州等珠三角城市的飛來峽水利工程樞紐的洪水調(diào)度和洪水資源的規(guī)劃利用具有重要意義。

2 基于Copula函數(shù)的聯(lián)合概率分布模型

2.1 Copula函數(shù)、參數(shù)估計(jì)和擬合優(yōu)度評價(jià)

根據(jù)Sklar定理，令F是具有單變量邊緣分布函數(shù)F1,...,Fn的n維分布函數(shù)，若邊緣分布函數(shù)F1,...,Fn連續(xù)，則存在一個(gè)唯一滿足F(x1,...,xn)=C(F1(x1),...,Fn(xn))關(guān)系的連接函數(shù)C。相反，如果C是一個(gè)n維Copula，F(xiàn)1,...,Fn為分布函數(shù)。作為一新的統(tǒng)計(jì)方法，Copula函數(shù)為描述隨機(jī)變量之間的相關(guān)結(jié)構(gòu)提供了一個(gè)新思路而得到了廣泛應(yīng)用。利用Copula函數(shù)建立多變量的聯(lián)合分布的突出優(yōu)點(diǎn)是：① Copula函數(shù)能夠把邊緣分布和變量間的相關(guān)關(guān)系分開處理，因此能夠?qū)吘壏植检`活選擇，選擇能較好地刻畫單變量自身的變化規(guī)律的邊緣分布；② Copula函數(shù)能夠捕捉隨機(jī)變量間的非線性相關(guān)關(guān)系, 且求解比較簡單，相對于線性相關(guān)提高了適用范圍；③ 由Copula函數(shù)導(dǎo)出的一系列的相關(guān)性度量指標(biāo)，拓展了變量間的相關(guān)性度量范圍，在實(shí)際工作中有更加廣泛的應(yīng)用。Nelsen[12]系統(tǒng)地總結(jié)了Copula函數(shù)的性質(zhì)和這個(gè)領(lǐng)域的主要研究成果。郭生練等較全面地介紹了Copula 函數(shù)在水文分析計(jì)算中的應(yīng)用[13]。

根據(jù)Copula函數(shù)性質(zhì)，構(gòu)建兩變量間的聯(lián)合概率分布模型可分兩步進(jìn)行：首先確定邊緣分布，然后選擇一個(gè)能夠恰當(dāng)?shù)胤从匙兞块g的相關(guān)結(jié)構(gòu)的Copula函數(shù)。研究表明，阿基米德(Archimedean)Copula是由其生成元唯一確定的單參數(shù)函數(shù)，是當(dāng)前應(yīng)用于水文氣象領(lǐng)域的一類非常重要的Copula函數(shù)。幾種常用二維Archimedean Copula函數(shù)為：

① Gumbel-Hougaard(GH) Copula函數(shù)：

C(u,v)=exp{-[(-lnu)θ+(-lnv)θ]1/θ},

θ∈[1,∞)

(1)

其生成元φ(t)為：φ(t)=(-lnt)θ。GHCopula函數(shù)僅適用于變量存在正相關(guān)的情形。

②Clayton Copula函數(shù)：

C(u,v)=(u-θ+v-θ-1)-1/θ,

θ∈(0,∞)

(2)

其生成元φ(t)為：φ(t)=(t-θ-1)/θ。同GH函數(shù)，適用于描述正相關(guān)的隨機(jī)變量。

③Ali-Mikhail-Haq (AMH)Copula函數(shù)：

C(u,v)=uv/[1-θ(1-u)(1-v)],

θ∈[-1,1)

(3)

其生成元φ(t)為：φ(t)=ln(1-θ(1-t)/t)。AMH函數(shù)能夠描述正或負(fù)相關(guān)的隨機(jī)變量，但不適用于正相關(guān)或負(fù)相關(guān)非常高的變量。

④Frank Copula函數(shù)：

C(u,v)=-ln{1+[(e-θu-1)(e-θv-1)/

(e-θ-1)]}/θ,θ∈R

(4)

其生成元φ(t)為：φ(t)=-ln[(e-θt-1)/(e-θ-1)]。與AMH Copula函數(shù)類似，但對相關(guān)程度沒有限制。

由于使用不同的Copula函數(shù)分析結(jié)果可能明顯不同，因此對Copula 函數(shù)的擬合優(yōu)度檢驗(yàn)以選擇合適的Copula函數(shù)顯得十分重要。擇優(yōu)選用Copula函數(shù)的主要檢驗(yàn)方法有由圖示直觀選擇Copula 函數(shù)的Genest-Rivest方法[14]、均方根誤差(RMSE)準(zhǔn)則法和AIC 信息準(zhǔn)則法等。

2.2 概率分布與重現(xiàn)期

根據(jù)Copula函數(shù)，兩變量聯(lián)合分布可表示為:

1999年11月，國家經(jīng)貿(mào)委與世界銀行召開“現(xiàn)代物流發(fā)展國際研討會”。時(shí)任國務(wù)院副總理吳邦國提出，“要把現(xiàn)代物流作為國民經(jīng)濟(jì)的重要產(chǎn)業(yè)和國民經(jīng)濟(jì)新的增長點(diǎn)?！?/p>

F(x,y)=P(X≤x,Y≤y)=

C(FX(x),FY(y))

(5)

變量X和Y的邊緣分布重現(xiàn)期：

(6)

變量X、Y的聯(lián)合重現(xiàn)期To：

(7)

變量X和的Y同現(xiàn)重現(xiàn)期Ta：

(8)

由以上可得：To(x,y)≤Min(T(x),T(y))≤Max(T(x),T(y))≤Ta(x,y)

當(dāng)給定X≥x時(shí)，Y≥y的條件概率：

(9)

同理，可定義Y≥y時(shí)，X≥x的概率。由條件概率分布，可計(jì)算相應(yīng)的條件重現(xiàn)期。

3 西北江洪水聯(lián)合概率分布

3.1 洪峰流量數(shù)據(jù)

基本資料為1959-2007年間馬口、三水兩測站的日平均流量。兩江洪水聯(lián)合概率分布分析所用的樣本分別以馬口歷年最大日均洪峰流量和同步的三水洪峰流量、三水歷年最大日均洪峰流量和同步的馬口洪峰流量兩組抽樣數(shù)據(jù)構(gòu)成，即：馬口-三水洪水(樣本1)、三水-馬口洪水(樣本2)。

3.2 聯(lián)合分布計(jì)算與分析

3.2.1 邊緣分布 函數(shù)對兩樣本的邊緣分布，都使用了以下3個(gè)三參數(shù)的概率分布函數(shù)擇優(yōu)：

1)廣義極值分布(GEV)：

FX(x)=P(X

(10)

2)皮爾遜三型分布(P-III)：

FX(x)=P(X

(11)

3)威布爾分布(WBL)：

FX(x)=P(X

(12)

式中，ξ,β,μ,α,a0為分布函數(shù)參數(shù)。參數(shù)估計(jì)使用概率權(quán)重矩(PWM)或線性矩(L-矩)方法。經(jīng)驗(yàn)頻率分布Pi使用Gringorten公式計(jì)算：Pi=(i-0.44)/(n+0.12)。擬合結(jié)果采用均方根誤差(RMSE)、經(jīng)驗(yàn)頻率和理論頻率擬合誤差平方和(Q)和概率點(diǎn)據(jù)相關(guān)系數(shù)(PPCC)檢驗(yàn)其擬合優(yōu)度。根據(jù)優(yōu)度檢驗(yàn)結(jié)果(表1)綜合比較，兩樣本的馬口和三水較優(yōu)邊緣分布分別選用P-III分布和GEV分布(圖1)。樣本2的邊緣分布線型和擬合優(yōu)度檢驗(yàn)結(jié)果基本同樣本1，為省去篇幅不再列出圖與表。

3.2.2 Copula 函數(shù)的參數(shù)估計(jì)及擬合優(yōu)度評價(jià) 由Kendall秩相關(guān)系數(shù)(τ)計(jì)算的馬口、三水洪水之間相關(guān)性分別為，樣本1：τ=0.857；τ=0.847。采用相關(guān)性指標(biāo)法計(jì)算了兩個(gè)聯(lián)合概率分布樣本的copula參數(shù)θ，并利用AIC、RMSE檢驗(yàn)其擬合優(yōu)度，結(jié)果見表2。根據(jù)擬合優(yōu)度評價(jià)指標(biāo)，選擇兩樣本中AIC和RMSE最小、KT-Ke關(guān)系圖中點(diǎn)據(jù)和理論直線最接近45°對角線的GH Copu1a函數(shù)作為聯(lián)合概率分布的連接函數(shù)(圖略)。GH Copu1a函數(shù)的θ=7.0(明顯大于1)，反映兩站洪水聯(lián)合概率分布之間存在很高的上尾相關(guān)性，即，當(dāng)其中一個(gè)測站洪水出現(xiàn)小頻率的洪峰流量時(shí)，另一個(gè)測站洪水出現(xiàn)小頻率的洪峰流量的概率隨之增大，此對于進(jìn)一步深入分析兩測站洪水的空間相關(guān)性和聯(lián)合概率分布特征提供了重要的基礎(chǔ)。擇優(yōu)構(gòu)建馬口、三水洪水聯(lián)合概率分布的GH Copula聯(lián)合分布函數(shù)如下：

表1 樣本1的邊緣分布參數(shù)與優(yōu)度檢驗(yàn)值

圖1 樣本1馬口洪水和三水洪水兩變量邊緣分布圖

樣本1：C1(FX(x),FY(y))=exp{-

[(-lnFX(x))7.00+(-lnFY(y))7.00]1/7.00}

(13)

樣本2：C2(FX(x),FY(y))=exp{-

[(-lnFX(x))6.53+(-lnFY(y))6.53]1/6.53}

(14)

表2 兩樣本的Copula參數(shù)及擬合優(yōu)度評價(jià)指標(biāo)

3.2.3 聯(lián)合概率分布與重現(xiàn)期 樣本1的聯(lián)合分布、聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期見圖2所示(樣本2的聯(lián)合概率分布與重現(xiàn)期圖與圖2基本相同，為省篇幅略去)不同重現(xiàn)期的洪水設(shè)計(jì)值見表3。

表3顯示，聯(lián)合重現(xiàn)期小于邊緣分布重現(xiàn)期，邊緣分布重現(xiàn)期則小于同現(xiàn)重現(xiàn)期。根據(jù)單變量同頻率方法推算的設(shè)計(jì)值小于兩站洪水聯(lián)合重現(xiàn)期設(shè)計(jì)值，差值隨重現(xiàn)期減小而有所增大，馬口洪水設(shè)計(jì)值相對差值大約介于0.8%～1.4%，三水洪水設(shè)計(jì)值相對差值大約介于0.6%～1.5%。受樣本抽樣方法影響，使得同頻率條件下樣本1的馬口洪水設(shè)計(jì)值略大于樣本2的設(shè)計(jì)值，樣本1的三水洪水設(shè)計(jì)值則略小于樣本2的設(shè)計(jì)值。在設(shè)計(jì)頻率的低頻部分，樣本2的同現(xiàn)重現(xiàn)期略大于樣本1的同現(xiàn)重現(xiàn)期。以200 a一遇洪水為例，樣本1中馬口、三水邊緣分布重現(xiàn)期設(shè)計(jì)值分別為55 960和18 790.m3/s，同頻率的聯(lián)合分布重現(xiàn)期設(shè)計(jì)值分別為564 310.和18 960.m3/s；樣本2中馬口、三水邊緣分布重現(xiàn)期設(shè)計(jì)值分別為55 150和18 220.m3/s，聯(lián)合分布重現(xiàn)期設(shè)計(jì)值分別為55 640和18 360.m3/s；樣本1邊緣分布200 a一遇的同現(xiàn)重現(xiàn)期為223年，樣本1邊緣分布200 a一遇的同現(xiàn)重現(xiàn)期為225 a。上述結(jié)果表明，與兩江洪水遭遇聯(lián)合概率分布比較，以單變量推算的馬口、三水洪水設(shè)計(jì)值實(shí)際上沒有達(dá)到設(shè)計(jì)標(biāo)準(zhǔn)。

3.2.4 條件概率分布 考慮到西江和北江流域爆發(fā)洪水的時(shí)間差異，需要進(jìn)一步分析馬口(三水)洪水出現(xiàn)某特定流量條件下，三水(馬口)洪水流量的概率分布。圖3表示當(dāng)馬口出現(xiàn)10%、5%、2%、1%、0.5%和0.2%超值概率的洪水時(shí)，三水遭遇洪水超過相應(yīng)標(biāo)準(zhǔn)洪水的條件概率P(Q三水≥q|Q馬口≥q|)。圖4表示當(dāng)三水出現(xiàn)超值概率為10%、5%、2%、1%、0.5%和0.2%洪水時(shí)，馬口遭遇的洪水超過相應(yīng)標(biāo)準(zhǔn)洪水的條件概率P(Q馬口≥q|Q三水≥q|)。表4說明，當(dāng)馬口(三水)洪水大于等于某一特定頻率設(shè)計(jì)值時(shí)，三水(馬口)出現(xiàn)大于等于該頻率設(shè)計(jì)值的條件概率隨著超值概率的減小而減小，以馬口出現(xiàn)大于等于概率P10%時(shí)的洪水為例，三水洪水出現(xiàn)大于等于概率P10%、P5%、P2%、P1%、P0.5%和P0.2%的條件概率分別為0.901、0.501、0.200、0.100、0.050和0.020，其相應(yīng)的條件重現(xiàn)期為：1.11、2、4.99、9.98、19.92和50.1。兩站相同設(shè)計(jì)頻率洪水相遇的概率都超過88%，充分表明馬口與三水洪水之間關(guān)系密切，同頻率洪水之間遭遇的概率非常高，由馬口或三水洪水設(shè)計(jì)值大致可推算出同頻率的三水或馬口洪水設(shè)計(jì)值。

圖2 樣本1聯(lián)合分布三維圖、聯(lián)合重現(xiàn)期等值線圖、同現(xiàn)重現(xiàn)期等值線圖

表3 馬口、三水洪水不同重現(xiàn)期的洪水設(shè)計(jì)值

Table 3 Flood design values for different return periods of Makou and Sanshui Stations

樣本樣本1:馬口-三水洪水樣本2:三水-馬口洪水邊緣分布聯(lián)合分布同現(xiàn)重現(xiàn)期設(shè)計(jì)值邊緣分布聯(lián)合分布同現(xiàn)重現(xiàn)期設(shè)計(jì)值t/aQ馬口 Q三水(104m3·s-1)T0(x,y) Ta(x,y)aQ馬口 Q三水(104m3·s-1)Q馬口 Q三水(104m3·s-1)T0(x,y)Ta(x,y)aQ馬口 Q三水(104m3·s-1)5006.0311.9184535586.0771.9291.9415.9274505631.9545.9742005.5961.8031812235.6431.8161.8225.5151802251.8365.5641005.2551.706911125.3041.7211.7215.191901131.7375.241504.9001.59845564.9511.6141.6104.85245561.6284.905204.4041.43618224.4581.4551.4444.37618221.4644.432103.9991.2969114.0561.3161.3013.9829111.3234.042

圖3 P(Q三水≥q|Q馬口≥q|)的條件概率分布圖 圖4 P(Q馬口≥q|Q三水≥q|)的條件概率分布圖

表4 兩站洪水遭遇條件概率

Table 4 Conditional probabilities of two station floods

樣本1超值概率P/%條件概率:P(Q三水≥q|Q馬口≥q|)105210.50.2條件重現(xiàn)期/a105210.50.2100.9010.9981.0001.0001.0001.0001.111.001.001.001.001.0050.5010.8990.9991.0001.0001.0002.001.111.001.001.001.0020.2000.4000.8970.9981.0001.0004.992.501.111.001.001.0010.1000.2000.4990.8970.9981.0009.985.002.001.111.001.000.50.0500.1000.2510.5010.8970.99919.929.973.992.001.121.000.20.0200.0400.1000.2000.3990.89450.1025.0710.035.012.511.12樣本2超值概率P/%條件概率:P(Q馬口≥q|Q三水≥q|)105210.50.2條件重現(xiàn)期/a105210.50.2100.8940.9971.0001.0001.0001.0001.121.001.001.001.001.0050.4990.8910.9991.0001.0001.0002.011.121.001.001.001.0020.2000.4000.8890.9971.0001.0005.002.501.121.001.001.0010.1000.2000.4980.8890.9971.00010.005.002.011.131.001.000.50.0500.1000.2500.4980.8880.99920.0010.004.002.011.131.000.20.0200.0400.1000.2000.3990.89850.0025.0310.015.002.511.11

4 結(jié) 語

利用由思賢滘相通的西江馬口和北江三水兩測站流量資料，本文分析了廣東西江與北江洪水之間的聯(lián)合概率分布特征，研究表明：兩站洪水之間存在高關(guān)聯(lián)性。基于Copula函數(shù)方法能夠選擇更優(yōu)的邊緣分布，構(gòu)造的聯(lián)合概率分布和條件概率更全面地反映了馬口和三水洪水極值統(tǒng)計(jì)特征和兩江洪水遭遇概率，為珠江西北江三角洲的防洪規(guī)劃與風(fēng)險(xiǎn)控制提出了一個(gè)較好的解決方案。

參考文獻(xiàn)：

[1]閆寶偉，郭生練，肖義，等. 基于兩變量聯(lián)合分布的干旱特征分析[J]. 干旱區(qū)研究， 2007， 24(4): 537- 542

[2]方彬，郭生練，肖義，等. 年最大洪水兩變量聯(lián)合分布研究[J].水科學(xué)進(jìn)展，2008， 19( 4):505-511

[3]陸桂華，閆桂霞，吳志勇，等.基于copula函數(shù)的區(qū)域干旱分析方法[J].水科學(xué)進(jìn)展，2010，21(2)：188-193

[4]劉曾美，陳子燊.區(qū)間暴雨和外江洪水位遭遇組合的風(fēng)險(xiǎn)[J]. 水科學(xué)進(jìn)展，2009, 20(5): 619- 625.

[5]ZHANG L, SINGH V P. Bivariate flood frequency analysis using the copula method [J]. Journal of Hydrologic Engineering, 2006, 11(2):150- 164.

[6]ZHANG L, SINGH V P. Bivariate rainfall frequency distributions using Archimedean copulas [J]. Journal of Hydrology, 2007, 332:93- 109.

[7]SHIAU J T, SONG F, NADARAJAH S. Assessment of hydrological droughts for the Yellow River, China, using copulas [J].Hydrological Processes, 2007, 21: 2157-2163.

[8]FAVRE A C, ADLOUNI S E, PERRAULT L, et al. Multivariate hydrological frequency analysis using Copulas [J].Water Resources Research, 2004, 40(11):1-12.

[9]JOE H. Multivariate Models and Dependence Concepts [M]. London:Chapman & Hall, 1997.

[10]童娟.珠江流域概況及水文特性分析[J].水利科技與經(jīng)濟(jì)，2007，13(1)：31-33

[11]姚章民，王永勇，李愛鳴. 珠江三角洲主要河道水量分配比變化初步分析[J]. 人民珠江，2009(2)：43-45，

[12]NELSEN R B. An Introduction to Copulas. Springer Series in Statistics [M]. NewYork: Springer,1999:216

[13]郭生練，閆寶偉，肖義，等.Copula 函數(shù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展[J].水文，2008，28(3)：1-7.

[14]GENEST C, RIVEST L. Statistical inference procedures for bivariate Archimedean copulas [J]. Journal of American Statistical Association, 1993, 88: 1034- 1043.