广东西江北江洪水联合概率分布研究*

陈子燊，刘曾美，路剑飞

(1. 中山大学水资源与环境系，广东 广州510275；2. 华南理工大学水利水电工程系，广东 广州510640)

1 研究背景

多变量极端事件之间的相关性普遍存在于水文气象领域。然而，至今对于极端水文气象事件的概率分布研究多数还是基于单变量的统计分析。一些研究尝试利用多变量正态联合分布，但由于要求边缘分布必须是正态或对数正态分布，因此限制了模型使用和选择范围。近些年来借助于Copula函数，多变量联合分布已经在水文气象极值事件分析中得到了较深入的研究[1-9]，为更精确地推算极端事件的重现水平和控制工程风险提供了一个重要的理论基础和技术手段。

流入广东省的西江发源于云南省霑益县马雄山，至广西梧州会桂江后始称西江，东流至三水市思贤滘与北江相通。北江发源于江西省信丰县石碣大茅坑，流入韶关市与武江汇合后始称北江，此后向南流经清远市，经思贤滘与西江干流连通后进入珠江三角洲网河区。位于珠江三角洲西北部顶点思贤滘两侧的马口站和三水站属于国家重点水文观测站，分别监测进入西北江三角洲网河区的流量过程。20世纪80年代以来，珠江三角洲经济和城市化进程高速发展，三角洲出现了大规模无序的河道采沙，河道发生了大规模的下切，而三水测站所处的河道断面下切幅度远远大于同期马口测站所在的河道断面下切幅度。由于思贤滘的连通作用，造成三水断面的水量分配比不断增大[10 ]。西北江洪水受流域特征制约，两江洪水特性差异明显。每年洪水发生时间通常北江早于西江，两江同时发洪的情况也时有发生[11]。通过马口、三水的两江洪水直接威胁珠江西北江三角洲的防洪安全，深入分析两江洪水遭遇的联合概率分布特征对于防御西、北江洪水拱卫广州等珠三角城市的飞来峡水利工程枢纽的洪水调度和洪水资源的规划利用具有重要意义。

2 基于Copula函数的联合概率分布模型

2.1 Copula函数、参数估计和拟合优度评价

根据Sklar定理，令F是具有单变量边缘分布函数F1,...,Fn的n维分布函数，若边缘分布函数F1,...,Fn连续，则存在一个唯一满足F(x1,...,xn)=C(F1(x1),...,Fn(xn))关系的连接函数C。相反，如果C是一个n维Copula，F1,...,Fn为分布函数。作为一新的统计方法，Copula函数为描述随机变量之间的相关结构提供了一个新思路而得到了广泛应用。利用Copula函数建立多变量的联合分布的突出优点是：① Copula函数能够把边缘分布和变量间的相关关系分开处理，因此能够对边缘分布灵活选择，选择能较好地刻画单变量自身的变化规律的边缘分布；② Copula函数能够捕捉随机变量间的非线性相关关系, 且求解比较简单，相对于线性相关提高了适用范围；③ 由Copula函数导出的一系列的相关性度量指标，拓展了变量间的相关性度量范围，在实际工作中有更加广泛的应用。Nelsen[12]系统地总结了Copula函数的性质和这个领域的主要研究成果。郭生练等较全面地介绍了Copula 函数在水文分析计算中的应用[13]。

根据Copula函数性质，构建两变量间的联合概率分布模型可分两步进行：首先确定边缘分布，然后选择一个能够恰当地反映变量间的相关结构的Copula函数。研究表明，阿基米德(Archimedean)Copula是由其生成元唯一确定的单参数函数，是当前应用于水文气象领域的一类非常重要的Copula函数。几种常用二维Archimedean Copula函数为：

① Gumbel-Hougaard(GH) Copula函数：

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

θ∈[1,∞)

(1)

其生成元φ(t)为：φ(t)=(-lnt)θ。GHCopula函数仅适用于变量存在正相关的情形。

②Clayton Copula函数：

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

θ∈(0,∞)

(2)

其生成元φ(t)为：φ(t)=(t-θ-1)/θ。同GH函数，适用于描述正相关的随机变量。

③Ali-Mikhail-Haq (AMH)Copula函数：

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

θ∈[-1,1)

(3)

其生成元φ(t)为：φ(t)=ln(1-θ(1-t)/t)。AMH函数能够描述正或负相关的随机变量，但不适用于正相关或负相关非常高的变量。

④Frank Copula函数：

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函数类似，但对相关程度没有限制。

由于使用不同的Copula函数分析结果可能明显不同，因此对Copula 函数的拟合优度检验以选择合适的Copula函数显得十分重要。择优选用Copula函数的主要检验方法有由图示直观选择Copula 函数的Genest-Rivest方法[14]、均方根误差(RMSE)准则法和AIC 信息准则法等。

2.2 概率分布与重现期

根据Copula函数，两变量联合分布可表示为:

1999年11月，国家经贸委与世界银行召开“现代物流发展国际研讨会”。时任国务院副总理吴邦国提出，“要把现代物流作为国民经济的重要产业和国民经济新的增长点。”

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

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

(5)

变量X和Y的边缘分布重现期：

(6)

变量X、Y的联合重现期To：

(7)

变量X和的Y同现重现期Ta：

(8)

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

当给定X≥x时，Y≥y的条件概率：

(9)

同理，可定义Y≥y时，X≥x的概率。由条件概率分布，可计算相应的条件重现期。

3 西北江洪水联合概率分布

3.1 洪峰流量数据

基本资料为1959-2007年间马口、三水两测站的日平均流量。两江洪水联合概率分布分析所用的样本分别以马口历年最大日均洪峰流量和同步的三水洪峰流量、三水历年最大日均洪峰流量和同步的马口洪峰流量两组抽样数据构成，即：马口-三水洪水(样本1)、三水-马口洪水(样本2)。

3.2 联合分布计算与分析

3.2.1 边缘分布 函数对两样本的边缘分布，都使用了以下3个三参数的概率分布函数择优：

1)广义极值分布(GEV)：

FX(x)=P(X

(10)

2)皮尔逊三型分布(P-III)：

FX(x)=P(X

(11)

3)威布尔分布(WBL)：

FX(x)=P(X

(12)

式中，ξ,β,μ,α,a0为分布函数参数。参数估计使用概率权重矩(PWM)或线性矩(L-矩)方法。经验频率分布Pi使用Gringorten公式计算：Pi=(i-0.44)/(n+0.12)。拟合结果采用均方根误差(RMSE)、经验频率和理论频率拟合误差平方和(Q)和概率点据相关系数(PPCC)检验其拟合优度。根据优度检验结果(表1)综合比较，两样本的马口和三水较优边缘分布分别选用P-III分布和GEV分布(图1)。样本2的边缘分布线型和拟合优度检验结果基本同样本1，为省去篇幅不再列出图与表。

3.2.2 Copula 函数的参数估计及拟合优度评价 由Kendall秩相关系数(τ)计算的马口、三水洪水之间相关性分别为，样本1：τ=0.857；τ=0.847。采用相关性指标法计算了两个联合概率分布样本的copula参数θ，并利用AIC、RMSE检验其拟合优度，结果见表2。根据拟合优度评价指标，选择两样本中AIC和RMSE最小、KT-Ke关系图中点据和理论直线最接近45°对角线的GH Copu1a函数作为联合概率分布的连接函数(图略)。GH Copu1a函数的θ=7.0(明显大于1)，反映两站洪水联合概率分布之间存在很高的上尾相关性，即，当其中一个测站洪水出现小频率的洪峰流量时，另一个测站洪水出现小频率的洪峰流量的概率随之增大，此对于进一步深入分析两测站洪水的空间相关性和联合概率分布特征提供了重要的基础。择优构建马口、三水洪水联合概率分布的GH Copula联合分布函数如下：

表1 样本1的边缘分布参数与优度检验值

图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参数及拟合优度评价指标

3.2.3 联合概率分布与重现期 样本1的联合分布、联合重现期和同现重现期见图2所示(样本2的联合概率分布与重现期图与图2基本相同，为省篇幅略去)不同重现期的洪水设计值见表3。

表3显示，联合重现期小于边缘分布重现期，边缘分布重现期则小于同现重现期。根据单变量同频率方法推算的设计值小于两站洪水联合重现期设计值，差值随重现期减小而有所增大，马口洪水设计值相对差值大约介于0.8%～1.4%，三水洪水设计值相对差值大约介于0.6%～1.5%。受样本抽样方法影响，使得同频率条件下样本1的马口洪水设计值略大于样本2的设计值，样本1的三水洪水设计值则略小于样本2的设计值。在设计频率的低频部分，样本2的同现重现期略大于样本1的同现重现期。以200 a一遇洪水为例，样本1中马口、三水边缘分布重现期设计值分别为55 960和18 790.m3/s，同频率的联合分布重现期设计值分别为564 310.和18 960.m3/s；样本2中马口、三水边缘分布重现期设计值分别为55 150和18 220.m3/s，联合分布重现期设计值分别为55 640和18 360.m3/s；样本1边缘分布200 a一遇的同现重现期为223年，样本1边缘分布200 a一遇的同现重现期为225 a。上述结果表明，与两江洪水遭遇联合概率分布比较，以单变量推算的马口、三水洪水设计值实际上没有达到设计标准。

3.2.4 条件概率分布 考虑到西江和北江流域爆发洪水的时间差异，需要进一步分析马口(三水)洪水出现某特定流量条件下，三水(马口)洪水流量的概率分布。图3表示当马口出现10%、5%、2%、1%、0.5%和0.2%超值概率的洪水时，三水遭遇洪水超过相应标准洪水的条件概率P(Q三水≥q|Q马口≥q|)。图4表示当三水出现超值概率为10%、5%、2%、1%、0.5%和0.2%洪水时，马口遭遇的洪水超过相应标准洪水的条件概率P(Q马口≥q|Q三水≥q|)。表4说明，当马口(三水)洪水大于等于某一特定频率设计值时，三水(马口)出现大于等于该频率设计值的条件概率随着超值概率的减小而减小，以马口出现大于等于概率P10%时的洪水为例，三水洪水出现大于等于概率P10%、P5%、P2%、P1%、P0.5%和P0.2%的条件概率分别为0.901、0.501、0.200、0.100、0.050和0.020，其相应的条件重现期为：1.11、2、4.99、9.98、19.92和50.1。两站相同设计频率洪水相遇的概率都超过88%，充分表明马口与三水洪水之间关系密切，同频率洪水之间遭遇的概率非常高，由马口或三水洪水设计值大致可推算出同频率的三水或马口洪水设计值。

图2 样本1联合分布三维图、联合重现期等值线图、同现重现期等值线图

表3 马口、三水洪水不同重现期的洪水设计值

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

样本样本1:马口-三水洪水样本2:三水-马口洪水边缘分布联合分布同现重现期设计值边缘分布联合分布同现重现期设计值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条件重现期/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条件重现期/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 结 语

利用由思贤滘相通的西江马口和北江三水两测站流量资料，本文分析了广东西江与北江洪水之间的联合概率分布特征，研究表明：两站洪水之间存在高关联性。基于Copula函数方法能够选择更优的边缘分布，构造的联合概率分布和条件概率更全面地反映了马口和三水洪水极值统计特征和两江洪水遭遇概率，为珠江西北江三角洲的防洪规划与风险控制提出了一个较好的解决方案。

参考文献：

[1]闫宝伟，郭生练，肖义，等. 基于两变量联合分布的干旱特征分析[J]. 干旱区研究， 2007， 24(4): 537- 542

[2]方彬，郭生练，肖义，等. 年最大洪水两变量联合分布研究[J].水科学进展，2008， 19( 4):505-511

[3]陆桂华，闫桂霞，吴志勇，等.基于copula函数的区域干旱分析方法[J].水科学进展，2010，21(2)：188-193

[4]刘曾美，陈子燊.区间暴雨和外江洪水位遭遇组合的风险[J]. 水科学进展，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].水利科技与经济，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 函数在多变量水文分析计算中的应用及研究进展[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.