Copula函数在水文多变量分析计算中的问题

宋松柏

（1.西北农林科技大学 水利与建筑工程学院，陕西杨凌712100；2.西北农林科技大学旱区农业水土工程教育部重点实验室，陕西杨凌712100）

水文事件一般具有多种相依性。选择一种特征属性变量，应用单变量频率分析难以为水利工程规划、设计和管理提供全面的水文设计值信息。传统多变量联合分布要求其变量边际分布必须为同一类型的分布，而实际水文现象具有高度的非线性和复杂性，同一水文事件的多种特征属性变量不可能服从同一类型的分布。因而，传统多变量联合分布的这种不足限制了其应用［1-5］。1959年以来，随着 Copula函数理论的发展，基于Copula函数的多变量分析开始受到重视。20世纪 90年代以后，Copula函数被引入洪水［1-14］、暴雨［1-4］、干旱［15-29］、径流预测［30］、径流模拟［31］、枯水概率［32］、水 文 预 报 不 确 定性［33-34］、 水文 丰枯 遭 遇 特征［35-36］、水量水质联合分布［37］、水资源供需风险损失［38］、水沙关系演变［39］、外调水供水补偿特性［40］、缺水风险［41］、用水效率［42］、设计潮位［43-45］和大坝洪水漫顶风险率［46］等的联合概率分析计算。从研究文献来看，目前在水文学中应用的Copula函数主要有：对称Archimedean Copula以及非对称的Archimedean Copula、Plackette Copula、Metaelliptical Copula 和混合Copula 等［2，5］。 Copula 函数参数估计方法有精确极大似然法、边际函数推断法和半参数法。其中：边际函数推断法也称两阶段法（two stage method），是目前水文学中广泛采用的Copula函数参数估计方法。按照边际函数推断法，基于Copula函数的水文多变量联合概率计算步骤主要有［2］：①单变量的边际分布参数估计和拟合度检验；②计算单变量的概率分布值；③变量间的相依性度量指标计算；④选择Copula函数；⑤Copula函数参数估计和拟合度检验；⑥多变量联合概率、条件概率以及相应重现期计算。武汉大学郭生练教授和熊立华教授在国内率先开展了基于Copula函数的水文分析计算［1-4］。2017年8月，国际统计水文学委员会副主席、河海大学陈元芳教授主持举办了河海大学“Copula在水文与环境科学中的应用”国际课程班，培训期间，国际水文科学协会国际统计水文学委员会主席Salvatore Grimaldi教授来华授课。所有这些工作推动了Copula在水文科学领域的应用。目前，国内代表性的专著有《Copulas函数及其在水文中的应用》《Copula函数理论在多变量水文分析计算中的应用研究》《Copulas and Its Application in Hydrology and Water Resources》以及《非一致性水文概率分布估计理论和方法》等。

基于Copula函数的多变量水文频率已有20多a的研究和实践历史，是水文学研究的活跃领域之一，经过国内外水文科学工作者的不懈努力，丰富了Copula函数的理论体系，拓宽了Copula函数的应用领域。基于Copula函数的多变量水文频率计算也是一个崭新的研究领域，涉及高等数学、概率论与数理统计、水文学、数值计算、优化计算等学科的交叉和渗透，面临一系列亟待解决的科学问题。本文不再重复归纳总结水文频率计算研究进展，而是重点探讨Copula函数进行水文多变量联合概率计算的几个问题，以期促进基于Copula函数的多变量水文频率计算理论的进一步发展。

1 Copula函数的主要类型

Copula函数分类很多。按照参数的多少，可以分为单参数Copula函数和多参数Copula函数。按照变量间相依性特性，可以分为对称Copula函数和非对称Copula函数。本节主要介绍水文中应用的几种Copula函数。

1.1 对称Archimedean Copula函数

对称 Archimedean Copula（Symmetric Archimedean Copula）也称为可交换 Copula（Exchangeable Archimedean Copula，EAC）。 对称 Archimedean Copula函数主要有20种。这类函数的特点是Copula函数仅含有一个参数，函数构造简单，相对而言求解容易［2］。其中：Cook-Johnson（Clayton）Copula函数、Gumbel-Hougaard Copula函数、Ali-Mikhail-Haq Copula函数和Frank Archimedean Copula函数是水文多变量频率计算中应用最多的函数。这4种常见的对称Archimedean Copula函数的生成函数φ（t）、参数θ取值范围和Copula分布函数（CDF）表达式见表1（t为生成函数自变量，d为Copula函数维数，uj为第j个变量的边际分布，R为实数）。

表1 常见的4种对称Archimedean Copula函数形式

因为对称Archimedean Copula函数只含有一个参数，仅用一个生成函数描述正的相依性，所以要求变量为对称相依。实际上，这种假定通常是不合理的。如洪水事件可用洪量、洪水历时和洪峰流量来描述，显然，它们任意两两变量间的相依性是不对等的。表1中二维对称Archimedean Copula函数参数与Kendall’s相关系数间存在一定的函数关系，可先计算出样本的Kendall’s相关系数，再推算Copula函数参数。但是，三维以上的对称Archimedean Copula函数需要利用极大似然法，通过迭代法计算获得Copula函数参数。

1.2 非对称Archimedean Copula函数

非对称Archimedean Copula函数通过若干个生成函数描述变量间的相依性，克服了对称Archimedean Copula函数的缺陷，是水文多变量分析计算较为合理的选择函数，这些函数参数的计算需要利用极大似然法。目前，水文分析中最为常用的非对称Archimedean Copula函数有嵌套 Archimedean（Nested Archimedean Copula，简称NAC）构造函数、层次Archimedean Copula（Hierarchical Archimedean Copula，简称 HAC）构造函数和配对 Copula（Pair-Copula Copula，简称 PCC）构造函数等［2］。

1.2.1 嵌套Archimedean构造函数

这类函数主要有M3函数、M4函数、M5函数、M6函数和 M12函数，见式（1）～式（4）。

（1）M3函数。其表达式为

式中：u1、u2、u3分别为边际分布值；θ1、θ2分别为 Copula函数参数， θ2≥ θ1∈ ［0，∞） 。

（2）M4函数。其表达式为

式中： θ2≥ θ1∈ ［0，∞） 。

（3）M5函数。其表达式为

式中： θ2≥ θ1∈ ［1，∞） 。

（4）M6函数。其表达式为

式中： θ2≥ θ1∈ ［1，∞） 。

（5）M12函数。其表达式为

式中： θ2≥ θ1∈ ［1，∞） 。

1.2.2 层次Archimedean Copula构造函数

层次 Archimedean Copula（Hierarchical Archimedean Copulas）构造函数又称为广义嵌套Archimedean构造函数，其一般形式主要基于嵌套多元Archimedean Copula函数的框架。每一级Archimedean Copula函数由前一级聚集构成，顶层级最后以层次Archimedean Copula函数结束，形成d维标准均匀随机变量（U1，…，Ud）的联合分布，联合分布值可在点 u＝（u1，…，ud）∈［0，1］d估算［2］。

1.2.3 配对Pair-Copula构造函数

配对Pair-Copula构造函数是把多元密度分解为d（d-1）/2个二维Copula密度函数，其中前d-1个为无条件Copula函数，其余为条件Copula函数，主要有Canonical vines和 D-vines两类结构［2］。 其中：D-vines结构的表达式为

Canonical vines结构的表达式为

式中：fk为第k个变量的边际概率密度；xk为第k个边际变量；ci，i+j为第i个边际变量与第i+j个边际变量构成的Copula密度函数；xi+j为第i+j个边际变量；xj为第c14｜23表示 c14｜23（F（x1｜x2，x3），F（x4｜x2，x3）） 。

1.2.4 三维Plackett Copula函数

设三变量（u，v，w），其两两交乘比率分别为二维Copula分布 CUV、CVW和 CUW，则三维 Plackett Copula参数 ψUVW定义为［2］

其中

式中：在给定ψUVW（u，v，w）下，ψUV、ψVW和ψUW分别为二维Plackett Copula函数CUV、CVW和CUW的参数。

记z＝ CUVW，则三维Plackett Copulas分布CUVW（u，v，w）为

其中

给定ψUV、ψVW、ψUW和ψUVW，三维Plackett Copula分布 CUVW（ u，v，w）可用式（8）、式（9）表示，但是三维Plackett Copula的分布密度较为复杂，其参数可采用极大似然法进行计算。不难看出，三维Plackett Copula分布有3个参数，属于不对称Copula函数。

1.2.5 Metaelliptical Copula函数

参数为 μ （p×1） 和Σ（p×p）（p为边际变量维数）的d维随机变量z具有elliptical类分布（elliptical distribution，ED），可定义为［2］

式中：r≥0，为随机变量；u为Rd上的均匀分布变量，且独立于r；A为d×d的常数矩阵，且满足A AT＝Σ；符号“＝d”的含义是式（12）两边具有相同的分布。

当r有密度函数时，z的密度函数可以表示为

式中g（·）为一个尺度函数（Scale Function）。

常见的d维对称elliptical类分布见表2。

表2 常见的d维对称elliptical类分布

从表2可以看出，elliptical类分布函数非常复杂，也属于不对称Copula函数，只能采用极大似然法进行参数求解。

2 Copula函数在水文多变量分析中面临的几个问题

根据文献报道和应用实践，Copula函数仍然处在不断发展和完善的阶段，基于Copula函数的水文多变量分析计算中面临的主要问题有以下几方面。

2.1 变量边际分布值计算

按照两阶段法估计参数，Copula函数计算首先需要进行边际变量的分布函数值计算。由于水文事件的概率（频率）分布函数未知，气候变化和高强度人类活动影响使水文数据难以满足独立和同分布条件，以及水文序列长度有限等影响，边际变量的分布参数值计算受到许多挑战［5-6］，因此边际变量函数值计算是影响Copula函数选择和参数计算的关键步骤。

2.2 Copula函数参数计算

选择Pearson古典相关系数rn、Spearman秩相关系数ρn、Kendallτ系数、Chi图和K图方法进行变量间的相依性度量。若相依性存在，则可计算Copula函数参数，否则按变量独立进行多变量联合概率计算。除表2常见的4种二维对称 Archimedean Copula外，Copula函数参数需要应用极大似然法进行估计。一般来说，对称 Archimedean Copulas、非对称 Archimedean copulas和Plackette Copula分布函数已知，但是其密度函数的推求过程比较复杂。因此，正确的样本对数极大似然函数参数的偏导数方程组是保证极大似然法正确求解的关键。Metaelliptical Copula函数的密度函数见表2，以下列出代表性Copula函数的密度函数表达式［2］。

2.2.1 常见的对称三维Copula函数的密度函数

（1）Gumbel-Hougaard Copula。 其密度函数表达式为

式中：w ＝ （-ln u1）θ+ （-ln u2）θ+ （-ln u3）θ。

（2） Clayton （Cook-Johnson） Copula。 其密度函数表达式为

（3）Frank Copula。其密度函数表达式为

2.2.2 常见的非对称Archimedean Copula函数的密度函数

（1）M3 Copula。其密度函数表达式为

式中： w ＝ 1- （1-e-θ2）-1（1-e-θ2u1）（1-e-θ2u2） ； G ＝1-e-θ1。

（2）M4 Copula。其密度函数表达式为

（3）M5 Copula。其密度函数表达式为

（4）M6 Copula。其密度函数表达式为

（5）M12 Copula。 其密度函数表达式为

2.2.3 PlacketteCopula密度函数

对于二维PlacketteCopula函数，其密度函数为

式中：ψUV为二维Plackette Copula参数。

对于三维Plackette Copula函数，其密度函数推导极其复杂，函数表达式较长，本文不再列出，具体表达式见文献［2］。

2.3 Copula函数值计算

Copula函数参数求解后，对称 Archimedean Copula、PlacketteCopula和完全嵌套的非对称Archimedean Copula的分布函数具有显函数的表达式，可直接进行Copula函数值计算。但是，非对称Archimedean Copula函数中的三维 Pair Copula函数以及Metaelliptical Copula函数需要数值积分才能获得Copula 函数值［2］。

2.3.1 三维Pair Copula函数值计算

三维Pair Copula函数构造灵活，但是其分布函数值计算起来比较困难，可采用数值积分计算：

通过积分变量代换，有

显然式（23）为一个二维条件Copula分布的一维积分，可应用高斯数值积分求解出三维概率分布。

2.3.2 Metaelliptical Copula函数值计算

这类Copula函数的密度函数含有特殊类函数，无法表示为显函数形式，只能通过数值积分进行求解计算。Kotz和 Nadarajah（2001）、Nadarajah和 Kotz（2007）推导了二维对称Kotz type分布密度和分布的超几何级数，二维Pearson typeⅡ、Ⅶ类不完全Beta函数的边际分布表达式，详细计算步骤见文献［2］。

2.4 Copula函数的选择

Copula函数实际上是把几个相依边际变量的分布函数值［0，1］通过某一函数连接起来的函数，本身不具有严格的水文物理基础，也就是说Copula函数同样与单变量分布函数一样，它们都不能从物理意义上解释几个相依边际变量的联合发生概率。因此，Copula函数只能根据计算概率与实测数据变量的联合经验概率拟合效果进行定性评估，在此基础上应用分布拟合度法进行假设检验。联合经验频率与相依边际变量的联合分布有关，实际应用中，联合经验概率可采用Gringorten公式。显然，这是一种联合经验频率的近似计算。

Copula拟合度检验通用的方法为CPI Rosenblatt转换法。 设边际分布 FX1（x1），…，FXd（xd） 有联合分布 Copula 函数C（FX1（x1），…，FXd（xd））＝F（x1，x2，…，xd） ，按照 Copula函数模拟， Φ-1（Zi） （其中 i＝ 1，2，…，d）服从标准正态分布 N（0，1） ， Φ 为标准正态分布； Φ-1（Zi） 为标准正态分布 N（0，1） 的逆函数，则Rosenblatt转换检验法的步骤：①提出原假设 H0，即（X1，X2，…，Xd） 具有 C（FX1（x1），…，FXd（xd））＝ C（u1，u2，…，ud）；②选择检验方法，如 Kolmogorov检验、Cramér von Mises检验和Anderson Darling（AD）检验；③根据显著水平α，确定相应的临界值；④根据样本计算统计量的观测值；⑤比较统计量的观测值与临界值，对原假设H0进行判断。

在单变量分布拟合度检验中，当实测数据序列长度较大时，样本统计量分布显著水平α相应的临界值需要提取模拟序列第（1-α）%的分位数来确定。同样，Copula拟合度检验没有给出显著水平α相应的临界值表，也需要进行模拟试验进行确定。

3 基于Copula函数的水文多变量分析计算的建议

根据现有基于Copula函数的水文多变量分析面临的问题，建议深入开展以下研究工作，以期为我国水文频率计算提供参考和支撑。

3.1 非平稳多变量联合概率分布与设计值计算的实用化

基于Copula函数的水文多变量联合概率计算仍然假定边际变量满足独立、同分布条件。若边际变量不满足平稳性，则现有Copula函数的计算方法体系不能用于计算非平稳变量的联合概率与设计值。虽然一些学者提出了时变Copula函数的计算方法，但是其模型求解难度增大。因此，这些方法仍需在实用化方面进一步研究。另外，单变量情况下设计标准（设计重现期）的设计值可以进行清晰的确定并广泛应用于工程实践，但是多变量设计值则无法明确确定。Salvadori G等提出了基于联合概率值定义危险区域的Kendall重现期计算方法［47-50］，他们认为该方法比传统的重现期计算更为合理，其基本原理是将事件发生区域划分为超临界区域、临界面和亚临界区域3类，事件发生在超临界区域上的平均时间间隔长度为多变量事件重现期。这种计算非常复杂，目前仍然停留在研究层面，需要在实用化方面进一步研究。

3.2 提高Copula函数的参数估计精度

如上所述，Copula函数的参数估计方法大多采用极大似然法。按照极大似然法原理，一般需要求解样本似然函数对Copula函数参数偏导数的非线性方程组。若Copula函数参数较少、不含有特殊函数时，则非线性方程组求解相对容易。当Copula函数参数较多、非线性方程组含有特殊函数和边际变量非显式分布函数时，则求解非常复杂。另外，非线性方程组的初值选择有时会影响非线性方程组的求解。因此，当非线性方程组求解困难时，以样本似然函数为最大目标函数，考虑Copula函数参数取值范围，采用遗传算法、粒子群算法、蚂蚁群算法、混合蛙跳算法、头脑风暴算法、蜻蜓算法和水循环算法等智能优化求解方法，可能不失为一种求解途径。

3.3 提高随机变量和、差、积、商分布计算的精度

流域设计断面各部分洪水的地区组成、干旱分析中特征变量的复合运算、水资源评价和配置中研究断面来水与区间耗水量的组合分布分析等都可以归结为水文变量和、差、积、商的分布计算。传统多变量理论上可以通过数学变换和高维数值积分来获得，但是其求解非常复杂，难以保证积分值的可靠性。因此，应用Copula函数原理，研究提高水文变量和、差、积、商的分布计算精度是目前急需解决的科学问题。