基于copula函数的联合分布的治涝标准分析方法

基于copula函数的联合分布的治涝标准分析方法

0涝区暴雨与承灾体水位组合的致灾机理防止洪水的目的是防止雨水造成的地面游泳池，并立即消除农业、工业、商店、林业、渔业、交通和其他基础设施损失。然而,暴雨是随机的,不可能及时排除每一场暴雨形成的积水,从而彻底消除涝灾,只能以某一标准(重现期)的暴雨为基准,使发生不超过该标准的暴雨时不产生涝灾。因此,相关规范[1~3]都指出:治涝设计标准一般应以涝区发生一定重现期的暴雨不受灾为准。所谓不受灾,是指暴雨造成的积水状况(包括积水深度和积水时间)不能超过承灾体的承灾能力,显然,最终是否受灾以及成灾度的大小与承灾体的承灾能力和积水状况紧密相关。对确定的涝区、确定的承灾体,是否受灾只与积水状况有关,而积水状况不仅与降雨有关,而且与孕灾环境息息相关。孕灾环境包括孕灾内环境和孕灾外环境,孕灾内环境包括涝区的地形地貌,调蓄能力、排水系统等因素;孕灾外环境主要指承泄区水位。孕灾内环境在一段时间内一般是相对确定的,而孕灾外环境(承泄区水位)则是随机变化的。如果把涝区内环境看作为承灾体致灾的内因,暴雨和承泄区水位则是致灾的外因。因此,导致涝灾的因子不只是暴雨,还有承泄区水位。故对一个确定的涝区,治涝设施的组成和规模还与承泄区水位紧密相关,实际的治涝能力达不到设计暴雨标准。可见,对于具有两种致灾因子的治涝事件,单变量(暴雨)分析不能描述两种致灾因子相互之间的相关关系,无法全面反映涝灾的真实特征。因此,必须用基于涝区暴雨和承泄区水位两变量频率分析来研究涝区暴雨与承泄区水位的遭遇规律,从而研究选定的的涝区暴雨与承泄区水位组合下的真正治涝标准。多变量联合分析成为水文计算领域的一个研究热点[6~7]。闫宝伟、郭生练等对多变量水文分析计算的正态变换法、非参数方法和Copula函数法的适用性和合理性做了比较,Copula函数作为一种用于构造灵活的多变量联合分布的工具,无偏性和有效性等统计性能最好。Copula函数已经在水科学领域得到了越来越多的应用[8~12]。本文采用Copula函数建立涝区暴雨和承泄区水位的联合分布,并基于联合概率密度函数和联合分布来研究涝区暴雨和承泄区水位遭遇组合规律和治涝标准分析方法,并以位于广东省连江中上游的阳山涝区和位于珠江河口区的坦洲涝区为实例来对比分析。1流域内同时确定的洪水和承载区水位的联合分布1.1设计暴雨时段的确定治涝是及时排除由于暴雨所产生的地面积水,因此,应该采用年最大设计历时暴雨和相应承泄区水位(简称相应水位)来构建联合分布。从理论上,设计暴雨历时应等于涝区出口控制断面的最大汇流时间[1~3](即流域汇流时间),实践中选取与涝灾关系较好的降雨时段作为暴雨历时。设计暴雨历时一般采用1~3天。对于很大面积的涝区可以采取分片防治,也可避免出现排水路径长,治涝效果不理想。如对广东省而言,设计暴雨历时一般为24h。为了便于分析,此处设计暴雨历时长取为24h,结合实际资料情况,采用最大一天暴雨和相应水位来进行遭遇分析。若承泄区为非感潮河段或湖泊,相应水位可以在年最大暴雨当日和次日中选取最高日均水位;若承泄区为感潮河段或海域,相应水位可以在年最大暴雨当日和次日中选取最高潮水位。特别值得一提的是,由于实际一天雨量的计时情况与最大24小时存在的差异性,为保险起见,相应水位可在年最大暴雨当日、次日甚至前一日中选取。1.2阿基米德族中copula函数的选取及对其参数估计的意义年最大一天暴雨H、相应水位Z均为连续的随机变量,F(h,z)为变量H和Z的联合分布函数,其边缘分布函数为FH(h)和FZ(z),根据著名的Sklar定理,一定存在唯一的Copula函数C(u,v),使得:实质上是将二元联合分布分为两个独立的部分来分别处理:变量间的相关性结构和变量的边缘分布,其中相关性结构用Copula函数来描述。变量联合分布函数的性质和与之对应的Copula函数密切相关。Copula函数中,采用阿基米德族Copula函数描述水文变量的研究越来越多,并已经初具理论和应用基础[8~12]。每种Copula函数的参数θ都与用于度量水文变量相关性的指标Kendall秩相关系数τ存在着一定的关系。因此,可根据变量H和Z的相关程度合理选择阿基米德族中的Copula函数。Kendall秩相关系数τ的计算式如下:式中(xi,yi)是观测点据;sign(·)是符号函数,当(xi-xj)(yi-yj)>0时,sign(·)=1;当(xi-xj)(yi-yj)=0时,sign(·)=0;当(xi-xj)(yi-yj)<0时,sign(·)=-1。一般而言,承泄区为感潮河段和海域时,涝区暴雨h和承泄区水位Z的相关性较弱,可以选用AMHCopula函数或FrankCopula函数来进行探讨;承泄区为非感潮河段和湖泊时,涝区暴雨和承泄区水位一般存在一定的正相关关系,关系密切程度主要与承泄区流域大小、流域形状以及河道形态有关。可以根据区间暴雨量H和外江洪水位Z的秩相关系数τ的值来确定选用Copula函数,一般可以选用ClaytonCopula函数、GHCopula函数甚至FrankCopula函数来进行探讨。联合分布的参数估计采用两步法,第1步,根据区间暴雨H和外江洪水位Z的独立观测值系列,采用单变量分布的参数估计方法估计边缘分布FH(h)和FZ(z)的参数;第2步,根据联合观测值系列估计Copula函数的参数θ。联合分布的参数确定后,尚需根据评价指标选取最适合的Copula函数,可采用离差平方和最小准则(OLS)来评价Copula方法的有效性。2我国最高潮水位供给侧交通特征当涝区发生某一标准的暴雨时,可能遭遇到承泄区各种水位,但是遭遇的机率是不同的。设年最大一天暴雨H和相应水位Z的联合密度函数为f(h,z),边缘密度函数分别为fH(h)和fZ(z),相应的Copula函数C(u,v)的密度函数为c(u,v),则有:式中为便于不同地区的遭遇规律对比,通常分析涝区一定重现期的年最大暴雨与承泄区某一重现期的年最高水位遭遇的机率。对任一年最大暴雨h的重现期TH,可以按下式(4)、(5)分别求出重现期为TH的h的FH(h)和fH(h)。对任一年最高水位zM的重现期TM,先由随机变量年最高水位ZM的分布FM(zM)反求出年最高水位zM,然后令z=zM,便可求出与年最高潮水位zM相等的相应水位的FZ(z)和fZ(z)。zM的计算式为:对重现期为TM的最高水位zM,其相应于相应水位系列的重现期(简称为相应重现期)为:因此,利用年最大一天暴雨和相应潮位的边缘分布、概率密度函数,以及年最高潮位的分布函数便可以得到各种重现期的年最大一天暴雨与各种重现期的年最高潮位的遭遇机率。3设计背景和设计组合的确定如对于单一致灾因子的堤防防洪工程而言,防洪标准一般是指能防御的最大洪水(以洪峰流量来表征)的重现期,该洪水为防洪工程设计的依据,称为设计洪水。只要出现的洪水不超过设计洪水,则防洪工程是安全的,出现的大于该设计洪水的概率为防洪风险率,风险率即为重现期的倒数。若防洪标准为T年一遇,T年一遇的洪峰流量为设计洪峰流量q设,则防洪风险率Prisk=P(Q>q设)=1/T。对于有两种致灾因子(暴雨和承泄区水位)的治涝工程而言,由于设计暴雨可能与各重现期的承泄区水位相遭遇,只是遭遇的几率不同而已,因此若以涝区发生一定重现期的暴雨不受灾为准,则相组合的设计承泄区水位(简称治涝设计外水位)从理论上应该是无穷高,而实际上是不可能的,通常选取与设计暴雨遭遇可能性小的高水位。无论遭遇可能性多么小,只要不为零,则实际的治涝标准就小于设计暴雨标准。若以T年一遇的设计暴雨hT和与之相组合的治涝设计外水位z配为设计组合,则实际治涝风险率Prisk为:从式(8)可见,若要保证发生设计暴雨时不受灾,必须使P[(Z>z配)(H≤hT)]=0,即涝区暴雨不超过设计暴雨hT的情况下,承泄区水位超过治涝设计外水位z配的概率为零。P[(Z>z配)∩(H≤hT)]=P(H≤h,Z>z配)=FH(hT)-F(hT,z配)。则治涝风险率的计算式为:式中F(hT,z配)为联合分布函数F(h,z)在点(hT,z配)的函数值,F(hT,z配)=C(FH(hT),FZ(z配)),C(FH(hT),FZ(z配))为Copula函数C(u,v)在点(FH(hT),FZ(z配))的函数值。可以证明,治涝风险率等于联合风险率P[(H>hT)∪(Z>z配)]。故治涝标准T治为:显然,T治<TH(hT)=T,且T治<TZ(z配),T为设计暴雨hT的重现期;TZ(z配)为设计治涝外水位z配的相应重现期。因此,对任一假设的最大暴雨重现期TH(h)和最高水位重现期TM(z)设计组合,便可由式(4)计算该最大暴雨的不足概率FH(h),由式(6)计算出最高水位后再由相应水位的边缘分布求出FZ(z),然后求出其联合分布F(h,z),进而由式(10)计算该设计组合下的治涝标准。反之,对任一要求的治涝标准T治,也可反推出各种年最大暴雨重现期TH(h)和最高水位重现期TM(z)组合,即可得到T治的等值线。此外,对一给定的治涝设计标准,可在该等值线上确定TH(h)和TM(z)组合方案。4涝区和坦洲涝区的年最高潮位联合概率密度分析分别以广东省阳山县的阳山涝区和中山市的坦洲涝区为例来对比分析。阳山涝区位于连江流域的中上游,承泄区为连江;而中山市坦洲涝区位于珠江三角洲下游,承泄区为磨刀门水道。阳山涝区采用阳山气象站1967-2005年的历年最大一天雨量和阳山水位站相应日均水位来分析确定年最大一天雨量和相应水位的联合分布。其最大一天雨量H1和相应日均洪水位Z1的边缘分布FH1(h)和FZ1(z)均采用P-Ⅲ型分布能拟合得较好,其统计参数分别为:μH1=97mm,Cv=0.38,Cs=1.33以及μZ1=60.57m,Cv=0.024,Cs=0.63。由式(2)计算得最大一天雨量H1与相应日均洪水位Z1的Kendall秩相关系数τ=0.848,说明这两者存在较大的正相关性,可采用ClaytonCopula函数或GumbelCopula函数来描述其相关性结构。经分析检验,ClaytonCopula能更好地描述相关性结构,其联合分布F1(h,z)为:若阳山涝区最大一天雨量H1的边缘概率密度为fH1(h),其承泄区相应日均洪水位Z1的边缘概率密度为fZ1(z),则它们的联合概率密度函数f1(h,z)为:对阳山涝区承泄区年最高日均水位序列进行频率分析,用P-Ⅲ型分布能拟合得较好,得到其统计参数为:=61.10m,Cv=0.025,Cs=0.675。坦洲涝区采用神湾雨量站1960~2007年的历年最大一天雨量和灯笼山站的相应潮位来分析确定最大一天暴雨和其相应潮水位的联合分布。其最大一天雨量H2和相应潮位Z2的边缘分布FH2(h)和FZ2(z)均采用P-Ⅲ型分布能拟合得较好,其统计参数分别为:μH2=180mm,Cv=0.46,Cs=2.065以及μZ2=1.08m,Cv=0.31,Cs=0.165。由式(2)计算得涝区暴雨H2与承泄区潮水位Z2的Kendall秩相关系数τ=0.0352,说明两者相关性很弱,可采用AMHCopula函数或FrankCopula函数来描述其相关性结构。参照文献求得其联合分布F2(h,z)和联合概率密度函数f2(h,z)分别为:对坦洲涝区承泄区年最高潮位序列进行频率分析,用P-Ⅲ型分布也能拟合得较好,得到其统计参数为:=1.65m,Cv=0.18,Cs=1.77。为便于不同地区的遭遇规律对比分析,对涝区暴雨和承泄区水位均采用重现期来表征。对任一假定的最大暴雨重现期TH(h)和最高水位重现期TM(z)设计组合,可由式(4)和式(5)计算该最大暴雨的FH(h)和fH(h);由式(6)求出年最高水位z,然后求出与年最高潮水位z相等的相应水位的FZ(z)和fZ(z)。由式(11)、式(12)或式(13)、式(14)可分别求得阳山涝区和坦洲涝区年最大暴雨量和承泄区水位的联合分布概率和联合概率密度函数值,再由式(10)可求得各设计组合下的治涝标准。反之,对任一希求的治涝标准T治,也可反推出各种年最大暴雨重现期TH(h)和最高水位重现期TM(z)的组合,即可得到T治的等值线。阳山涝区和坦洲涝区不同重现期的年最大暴雨与年最高承泄区水位遭遇的联合概率密度分布图分别见图1的(a)和(b);部分计算成果见表1和表2;阳山涝区和坦洲涝区的治涝标准等值线图分别见图2的(a)和(b)。对比图1的(a)和(b)、表1和表2中的概率密度函数可知,对于阳山涝区不超过十年一遇的暴雨,遭遇同频率的外江水位的机率最大;但对超过十年一遇以上暴雨,最大遭遇机率的外江水位的的重现期低于暴雨重现期;而对于坦洲涝区,各重现期暴雨的最大遭遇机率的外江潮水位的重现期均不超过2年一遇。对比表1和表2可知,对同一标准设计暴雨和同一治涝标准,若遭遇可能性愈大,则组合的设计外水位重现期愈大;若遭遇可能性愈小,则相组合的设计外水位重现期愈小。例如对于20年一遇暴雨,若要治涝标准达15年一遇,则相组合的设计外水位的重现期,阳山涝区需21年一遇,而坦洲涝区为4.5年。而且可见,治涝标准即小于设