基于Copula熵方法的河流之间的相关性研究

Analysis of Correlation between River Flows Using Copula-Entropy Method

摘要: 研究长江干流及其支流之间的总相关性对于长江上游的水力设计、防洪和风险控制非常重要。针对现有的相关性计算方法的不足与缺陷，本文引入Copula熵方法，用来计算多变量之间的相关结构，并推导出Copula熵与互信息的关系与计算方法，采用多重积分法和蒙特卡罗法估计多变量之间的相关性特征。以长江上游的五条主要干支流：金沙江、岷江、沱江、嘉陵江和乌江为研究对象。计算结果表明，当使用不同的Copula函数时，总相关值存在显著差异。长江上游干支流河流之间的总相关性不大，其中岷江和沱江之间的相关性最大，金沙江，岷江和沱江之间也存在一定的相关性，因此，这几条河流有洪水遭遇的可能，对三峡大坝的防洪构成了威胁。

Abstract: Analysis of the dependence between the main stream and its upper tributaries is important for hydraulic design, flood prevention and risk control. In order to solve the disadvantages of the current hydrologic correlation analysis method, the method of copula entropy was introduced to estimate the dependence between Hydrological variables. The relationship between copula entropy and mutual information was discussed and the calculation procedures of copula entropy were deduced, and multiple integration and Monte Carlo methods were used to calculate the copula entropy. The upper Yangtze River was selected for case study. Results show that there is a significant difference in total correlation values, when different copula functions were used. The total correlation among the rivers is not high, and the one between Min and Tuo Rivers is the largest. There are some dependence among Jinsha, Min and Tuo Rivers, which constitutes a threat to flood control by the Three Gorges Dam (TGD). The flows of Jinsha, Jialing, Min and Tuo Rivers significantly influence the flood occurrence in the Yangtze River.

1. 引言

洪水事件一般具有多方面因素的特征属性，而各个变量之间通常具有相关性，研究变量之间的相关性对正确认识水文事件的客观规律，合理地进行水资源开发利用和制定有效的防洪减灾措施，具有十分重要的意义。目前，表达水文变量之间相关性的方法一般采用Pearson线性相关方法。该方法有以下几个缺点：① 只适用于线性相关关系，② 变量必须服从多元正态分布的假设；③ 无法计算多元变量之间的相关性。但在实际应用中，并不是所有的水文变量之间都是线性相关且满足服从正态分布的假设。

近年来基于熵理论的相关性分析方法备受推崇，许多研究者已经成功的利用互信息(Mutual Information)表征随机变量概率分布之间的相关性，该方法的优点是：① 它是一种非参数方法；② 不对边缘分布的函数形式作出假设，变量可以服从任何分布；③ 它可以用于多变量之间的相关性计算。该方法是一种实现多变量之间的相关性计算的有效途径，并在水文学领域中得到了广泛的应用。陈璐等 [1] 利用Copula熵方法对神经网络的径流预报模型中的预报因子进行了选择；李帆等 [2] 利用Copula熵方法对三变量洪水频率进行了分析；韩敏等 [3] 基于Copula熵方法对互信息进行了估计。本文以长江上游干支流为研究对象，采用信息论和Copula函数相结合的非线性技术，即Copula函数的熵，对长江上游干支流之间的相关性进行了分析，为长江上游的洪水遭遇以及三峡的防洪风险分析提供了依据。

2. Copula熵

2.1. Copula熵的定义

令 $x\in R_d$ 为d维随机变量，其边缘分布函数为 $F_i\left(X\right)$ ， $u_i=F_i\left(X\right)$ ， $i=\text{1},\text{2},\cdots ,d$ 。其中， $u_i$ 为服从均匀分布的随机变量，将Copula函数的熵定义为Copula熵，可以表示为：

$H_c\left(u_1,u_2,\cdots ,u_d\right)=-{\displaystyle \underset{0}{\overset{1}{\int }}\cdots {\displaystyle \underset{0}{\overset{1}{\int }}c\left(u_1,u_2,\cdots ,u_d\right)\log \left(c\left(u_1,u_2,\cdots ,u_d\right)\right)\text{d}{u}_1\text{d}{u}_2\cdots \text{d}{u}_d}}$ (1)

式中： $c\left(u_1,u_2,\cdots ,u_d\right)$ 为Copula函数的概率密度函数，可表示为 $\frac{\partial C\left(u_1,u_2,\cdots ,u_d\right)}{\partial u_1\partial u_2\cdots \partial u_d}$ 。下面探讨联合熵和Copula熵之间的关系，变量X的联合概率密度函数可以定义为 [4] ：

$f\left(x_1,x_2,\cdots ,x_d\right)=c\left(u_1,u_2,\cdots ,u_d\right){\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)}$ (2)

结合等式(1)，多维联合熵可表示为：

$\begin{array}{c}H\left(X_1,X_2,\cdots ,X_d\right)=-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x_1,x_2,\cdots ,x_d\right)\log \left[f\left(x_1,x_2,\cdots ,x_d\right)\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}\\ =-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}c\left(u_1,u_2,\cdots ,u_d\right){\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)\log \left[c\left(u_1,u_2,\cdots ,u_d\right){\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)}\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}}\\ =-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}c\left(u_1,u_2,\cdots ,u_d\right){\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)\left\{\log \left[c\left(u_1,u_2,\cdots ,u_d\right)\right]+{\displaystyle \underset{i=1}{\overset{d}{\sum }}\log \left[f\left(x_i\right)\right]}\right\}\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}}\\ =-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}c\left(u_1,u_2,\cdots ,u_d\right)}}{\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)\cdot \log \left[c\left(u_1,u_2,\cdots ,u_d\right)\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}\\ -{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}c\left(u_1,u_2,\cdots ,u_d\right){\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)\cdot {\displaystyle \underset{i=1}{\overset{d}{\sum }}\log \left[f\left(x_i\right)\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}}}\\ =A+B\end{array}$ (3)

其中A可表示为：

$\begin{array}{c}A=-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}c\left(u_1,u_2,\cdots ,u_d\right){\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)}\cdot {\displaystyle \underset{i=1}{\overset{d}{\sum }}\log \left[f\left(x_i\right)\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}}\\ =-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x_1,x_2,\cdots ,x_d\right)\cdot {\displaystyle \underset{i=1}{\overset{d}{\sum }}\log \left[f\left(x_i\right)\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}}\\ =-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x_1,x_2,\cdots ,x_d\right)\cdot \log \left[f\left(x_1\right)\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}\\ -\cdots -{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x_1,x_2,\cdots ,x_d\right)\cdot \log \left[f\left(x_d\right)\right]\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_d}}\\ =-{\displaystyle \underset{0}{\overset{\infty }{\int }}\log \left[f\left(x_1\right)\right]\cdot \left[{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x_1,x_2,\cdots ,x_d\right)\text{d}{x}_2\cdots \text{d}{x}_d}}\right]\text{d}{x}_1}\\ -\cdots -{\displaystyle \underset{0}{\overset{\infty }{\int }}\log \left[f\left(x_d\right)\right]\cdot \left[{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x_1,x_2,\cdots ,x_d\right)\text{d}{x}_1\text{d}{x}_2\cdots \text{d}{x}_{d-1}}}\right]\text{d}{x}_d}\\ =-{\displaystyle \underset{i=1}{\overset{d}{\sum }}{\displaystyle \underset{0}{\overset{\infty }{\int }}f\left(x_i\right)\log \left[f\left(x_i\right)\right]\text{d}{x}_i}}={\displaystyle \underset{i=1}{\overset{d}{\sum }}H\left(X_i\right)}\end{array}$ (4)

根据等式 $\text{d}u=\text{d}x\cdot f\left(x_i\right)$ ，

$\begin{array}{c}B=-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}c\left(u_1,u_2,\cdots ,u_d\right){\displaystyle \underset{i=1}{\overset{d}{\prod }}f\left(x_i\right)\cdot \log \left[c\left(u_1,u_2,\cdots ,u_d\right)\right]\text{d}{u}_1\text{d}{u}_2\cdots \text{d}{u}_d}}}\\ =-{\displaystyle \underset{0}{\overset{\infty }{\int }}\cdots {\displaystyle \underset{0}{\overset{\infty }{\int }}c\left(u_1,u_2,\cdots ,u_d\right)\cdot \log \left[c\left(u_1,u_2,\cdots ,u_d\right)\right]\text{d}{u}_1\text{d}{u}_2\cdots \text{d}{u}_d}}\\ =H_C\left(u\right)\end{array}$ (5)

所以，联合熵可表示为：

$H\left(X_1,X_2,\cdots ,X_d\right)={\displaystyle \underset{i=1}{\overset{d}{\sum }}H\left(X_i\right)+H_C\left(u_1,u_2,\cdots ,u_d\right)}$ (6)

由方程(6)可知，联合熵可分解为两部分：d维边缘分布熵的和与Copula函数的熵。

2.2. Copula熵的计算

2.2.1. 多重积分法

根据等式(1)，Copula熵可以使用多重积分法求得。首先，需要估计Copula函数的参数，从而确定Copula的概率密度函数。多重积分法由Berntsen等提出(1991)，用于计算多重积分，将Copula概率密度函数进行积分，得到相应的Copula熵值 [5] 。

2.2.2. 蒙特卡罗法

对于变量较多的情况，多重积分法较难求得，可采用蒙特卡罗法计算Copula熵。对于在[0,1]中的多元向量，Copula熵可以表示为：

$H_C\left(u_1,u_2,\cdots ,u_d\right)=-{\displaystyle \underset{{\left[0,1\right]}^d}{\int }c\left(U\right)\ln c\left(U\right)\text{d}U}=-E\left[\ln c\left(U\right)\right]$ (7)

Copula熵等于lnc(U)的期望值，可以通过蒙特卡罗法求得，与多重积分法类似，首先需确定Copula函数的相关性系数和参数，从而确定Copula函数，然后计算-lnc(U)的期望值。

3. 总相关

因为计算多变量联合概率时比较复杂，所以有必要计算多个变量之间的相关性。总相关性可以定义为 [6] ：

$I\left(X_1,X_2,\cdots ,X_d\right)={\displaystyle \underset{i=1}{\overset{d}{\sum }}H\left(X_i\right)-H\left(X_1,X_2,\cdots ,X_d\right)}$ (8)

可得，对于d = 2的情况，总相关性等效于两变量公共的互信息。

根据方程(6)，方程(8)可以写为：

$I\left(X_1,X_2,\cdots ,X_d\right)=-H_C\left(u\right)=-B$ (9)

因多变量的边缘熵的和总是大于多变量的联合熵，因此相关性总为正值；当且仅当所有被选择的变量相互独立时，I等于零。然而，如果两个变量有一些相关性，即使其余的变量是独立的，I仍大于零。根据等式(9)，随机变量的总相关性等于它们的负Copula熵。基于公式(9)采用Copula熵方法可直接计算总相关，避免了通过公式(8)计算联合熵和边缘熵带来的误差累积效应，且基于Copula熵计算的总相关性仅依赖于由Copula参数所确定的Copula函数。因此，对于总相关性的估计仅需要确定Copula参数。

4. 应用研究

4.1. 研究区域

长江上游汇集众多支流，主要是左岸的雅砻江、岷江、沱江和嘉陵江，右岸的乌江。雅砻江汇入金沙江，因此，本研究考虑了金沙江，而不考虑雅砻江。金沙江、岷江、沱江、嘉陵江、乌江从上游到下游的五个测站分别是屏山、高场、李家庄、北碚、武隆，采用各站1951年至2007年的日平均流量数据进行计算。采用年最大取样方法，应用最大熵原理方法获得边缘分布。Kendall和Pearson相关系数的计算结果如表1所示。

假设边缘分布是正态分布，通过正态分布和POME方法估计河流的边缘概率密度函数，如图1所示，可以看出，对于岷江和乌江，POME估计的分布所拟合经验分布的效果优于正态分布，两江数据表现出较高的峰度和偏度，所以在这种情况下，正态分布的假设是不合适的，然而，线性相关方法主要取决于边缘分布是正态分布的假设，因此，选择POME方法获得边缘分布是更加合理的。

Table 1. Dependence measures for the upper Yangtze River, China, based on annual maximum data

表1. 上游长江年最大值的相关性

注：右上角元素是Kendall相关系数；左下角元素是Pearson相关系数。

Figure 1. Results of fitting by POME method and normal distribution

图1. 通过POME方法和正态分布拟合的结果

4.2. 拟合优度检验

本文利用AIC准则对所生成的联合分布进行拟合优度检验，结果可以判断Copula能否可以很好的体现出变量之间的相关性，其中最小AIC为Copula函数选取的重要依据。AIC可以通过计算极大似然值或函数的均方误差来得到 [7] ，在一般情况下，AIC可以表示为：

$\text{AIC}=2k-2\ln \left(L\right)$ (10)

式中：k为统计模型中的参数的数量，L为估计模型的极大似然值。

4.3. 双变量模型

首先，确定长江上游任何两条河流的联合分布，根据该区域的5条河流，建立了10个双变量联合分布。Gumbel、Clayton、Frank、Normal和Student t Copula分别用于模拟五个站之间的相关关系。利用伪极大似然法估计这些Copula的参数，AIC值用于选择最合适的Copula函数。

所选的Copula和它们的参数如表2。一般来说，阿基米德Copula比椭圆Copula能更好的拟合联合分布函数。由Genest等人定义的Cramer-von Mises函数S_n (2009)用于拟合优度检验，计算其相应的P值，如表2所示，结果表明所选择的双变量Copula可用于模拟站点之间的相关关系。

分别用多重积分法和蒙特卡罗法计算Copula熵，对于第一种方法，采用Berntsen等人在1991年提出的多重积分法进行计算；对于第二种方法，产生10000对u，并计算 $\ln \left[c\left(u\right)\right]$ 的平均值，这两种方法的计算结果如表3

Table 2. Selections of copulas and determination of the parameters

表2. Copula的选择和参数的确定

注：右上角元素是选择的Copula；左下角元素是相对应的参数；括号中的值为P值。

Table 3. Total correlation values of two tributaries in upstream Yangtze River

表3. 长江上游两条支流的总相关值

注：右上角元素是多重积分法的值；左下角元素是蒙特卡罗法的值。

所示，结果表明两种计算结果相似。从表3中还可看出，总相关值并不大。与Pearson和Kendall相关性系数相比，Copula熵计算出的总相关性都为大于0的值，且一般小于前两种方法估计的值，根据总相关性的定义，Copula熵计算出的值更加可靠。对于Normal Copula，总相关可以表示为 [8] ：

$I_N=-\frac{1}{2}\log \left(1-{\rho }^2\right)$ (11)

式中：ρ是Pearson相关系数。

为了评估所提出方法的合理性，采用等式(11)和提出的两种方法分别计算Normal Copula的总相关性。对于双变量的情况，可以采用Pearson相关系数来估计Normal Copula的参数，这三种方法的计算结果如表4所示，结果表明通过这三种方法计算的总相关性非常接近。因此，对特定情况的分析，可以证明所提出的方法是合理的，并且可以用于相关性的分析计算。

4.4. 多变量模型

阿基米德类的三维Gumbel、Clayton和Frank Copula以及椭圆类的Normal Copula和t Copula可用于模拟三个站之间的相关性特征。由于不对称的阿基米德Copula只能模拟正相关，对有负相关的情况，只能采用Normal Copula和t Copula建立联合分布。利用多重积分法计算得到的总相关性结果如表5所示，表中是具有代表性的组合情况，结果表明，当选择不同的Copula函数时，总相关值存在差异，因此，选择合适的Copula函数对于估计变量之间的总相关性是十分重要的。

建立四条河流的联合分布，由于复杂的相关性结构，采用椭圆Copula，即Normal Copula和t Copula。利用伪极大似然法估计参数，估计的参数如表6所示。从表6可以看出，金沙江、岷江、沱江和嘉陵江的最大总相关性为0.36。由于金沙江、岷江、嘉陵江、沱江河流存在较大相关性，因此，这几天河流有洪水遭遇的可能，

Table 4. Comparisons of total correlation of Gaussian correlated variables in different methods

表4. 不同方法中相关变量之间的总相关性的比较

Table 5. Total correlation analysis of trivariate joint distribution

表5. 三维联合分布的总相关性分析

注：v表示这些模型的最大似然值，下同。

Table 6. Total correlation analysis of four-dimensional joint distribution

表6. 四维联合分布的总相关性分析

Table 7. Total correlation analysis of five-dimensional joint distribution

表7. 五维联合分布的总相关性分析

且对三峡的防洪构成了威胁。

建立了五条河流的联合分布，使用椭圆Copula，估计的参数如表7所示。计算的Normal Copula和t copula的总相关值为0.39，大于四变量总相关值。

5. 结论

本研究分析了长江上游五条河流的相关性。首先在水文领域引入由Copula和熵理论构建的Copula熵。由于变量之间存在复杂的相关性结构，所以采用Copula熵计算总相关性，即首先估计Copula函数的参数从而确定Copula函数，最后计算总相关值。主要结论如下：

1) 应用阿基米德和椭圆Copula建立多元变量的联合分布。一般阿基米德Copula更适合于维数少的情况，

对于维数多的情况，其相关性复杂，只能使用椭圆Copula，对于所建立的联合分布的理论曲线都能很好地拟合经验频率。

2) 基于信息理论和Copula函数，引入Copula熵方法，这是一种非参数方法，可以表示线性和非线性相关性，且不对边缘分布做出假设，可用于更高的维数。此外，该方法仅需要通过计算Copula熵来直接地估计总相关性，而不需要计算边缘熵和联合熵的值，从而避免了偏差的累积效应。

3) 使用多重积分法和蒙特卡罗法计算总相关值，得到的计算结果相近。对于特殊情况，可以采用不同类型的Copula函数，当使用不同的Copula函数时，总相关值存在显著差异，因此，选择合适的Copula函数对于相关性的估计是很重要的。

4) 计算结果表明，河流之间的总相关性不是很大，这与研究区域的气候特征有关。由于岷江和沱江这两条河流之间的距离最短，且属于同一暴雨区，所以它们之间的总相关系数最大为0.33。金沙江，岷江和沱江之间也存在一定的相关性，在平水年，由于长江左岸和右岸一般不会同时发生降雨，所以岷江和乌江的相关性不能忽视。由于金沙江，岷江，嘉陵江，沱江河流存在较大相关性，所以这几条河流有洪水遭遇的可能，且对三峡大坝的防洪构成了威胁。

基金项目

国家自然科学基金面上项目(51679094)；国家自然科学基金重点资助项目(91547208)；中央高校基本科研业务费专项资金资助(2017KFYXJJ194, 2016YXZD048)；中国水利水电科学研究院流域水循环模拟与调控国家重点实验室开放研究基金(IWHR-SKL-201408)资助。