黄河唐乃亥站径流干旱多变量联合特征分析

马明卫，程玉佳,2，崔惠娟，臧红飞，王梓宇

(1.华北水利水电大学 水资源学院,河南 郑州 450046; 2.河海大学 水文水资源学院,江苏 南京 210098;3.中国科学院地理科学与资源研究所 陆地表层格局与模拟重点实验室,北京 100101；4.华北水利水电大学 水利学院,河南 郑州 450046)

干旱是一种成因复杂的自然现象，不仅与自然因素息息相关，还受到经济、社会、农业条件等的影响。由于干旱具有形成慢、历时长的特点，所以经常容易被人忽视，但其潜在破坏性不亚于地震、洪涝等灾害。同时，也正是因为干旱具有上述特性，在信息较充分的情况下，可以对其进行预防，进而减轻其影响[1-2]。基于不同的角度，干旱大致被划分为气象、农业、水文及社会经济干旱等4类[3]。其中，水文干旱是指因降水缺乏导致地表水减少[4]，如河川径流、水库和湖泊蓄水等减少，可能造成严重的后果[5]。径流干旱是水文干旱的主要研究对象，通常可以基于游程理论识别径流干旱事件，在识别过程中所选取的阈值可以是固定的，也可以是随时间变化的。但考虑到丰水期与枯水期流量相差较大，若选取固定阈值来识别干旱事件显然是不合理的。此外，在识别干旱事件的过程中经常会出现相邻干旱事件时间间隔很短，以及大量历时很短或烈度很小的干旱事件的情况，它们都会影响干旱特征统计分析的客观性。因此，有必要选择多阈值来识别径流干旱事件并对它们进行融合剔除。

干旱事件具有多变量的特征属性，且各特征变量常服从不同的边缘分布，这给干旱多变量研究带来较大的困难，而Copula函数可以联结服从任意边缘分布的多个变量来构造它们的联合分布函数[6]。姚蕊等[7]基于Joe Copula函数分析淮河流域干旱历时和干旱烈度的联合频率特征，并对引起水文干旱特征频率变化的原因及影响因素进行分析。马建琴等[8]为了进一步研究沙颍河流域水文干旱特征及其演变规律，选取多种Archimedean Copula函数拟合干旱特征变量的联合分布，并计算干旱事件的联合重现期和同现重现期。涂新军等[9]通过构建水文干旱特征两变量联合分布Copula模型，来分析变化环境下东江流域水文干旱特征及缺水响应机制。黄河流域绝大部分属于干旱半干旱区，人口集中、地多水少、资源环境问题突出。习近平总书记指出，黄河流域在我国经济社会发展和生态安全方面具有十分重要的地位，保护黄河是事关中华民族伟大复兴的千秋大计。然而，黄河中上游流域，特别是黄河源区的自然生态十分脆弱，对气候变化极为敏感。近年来，在全球变暖的背景下，该区域气温上升明显，地表蒸散发量增大，径流量显著减少，生态环境遭到不同程度破坏。作为黄河流域的“水塔”，黄河源区的多年平均径流量约占整个黄河流域的35%。由于地处东亚季风区的边缘，气候的轻微改变都可能引起该区域较为强烈的反应，特别是西太平洋副热带高压强度和位置的变化，导致东南暖湿气流水分通量有很大的不确定性，因而，其降水在年际和季节尺度上振荡较为明显，造成黄河中下游历史上旱涝交替剧烈，干旱和洪水灾害频发。本文针对黄河源区，根据唐乃亥站的历史逐日流量资料，基于游程理论和可变日阈值识别径流干旱事件，并对干旱事件进行融合剔除，并采用多种Copula函数构建径流干旱多变量联合分布模型，以期为黄河源区干旱风险评估和水资源利用管理提供科学依据与数据支持。

1 研究区与数据

黄河源区位于青藏高原东北部，其水资源与生态环境状况关乎整个黄河流域的水安全和区域发展。作为黄河源区的门户，唐乃亥站(100°09′E，35°30′N)位于青海省兴海县唐乃亥乡的黄河干流上，该站可以提供反映黄河上游水文水资源时空变化和水沙变化的重要基础资料。唐乃亥站断面以上集水区域地理位置如图1所示。笔者选用该站1956—2018年的逐日平均流量实测资料，对黄河源区1956—2018年的历史径流干旱特征及其变化进行分析。

2 研究方法

2.1 径流干旱事件的识别与融合剔除

文中基于游程理论来识别一系列的径流干旱事件，通常将设定的流量作为阈值，当逐日流量序列在一定时段内持续小于给定阈值时则出现负游程，此时，则认为发生了径流干旱事件，如图2所示。

图2 径流干旱事件的识别与提取结果

一次径流干旱事件的特征变量通常包括干旱历时D(d)、干旱烈度S(亿m3)、干旱烈度峰值P(m3/s)等。其中，干旱历时反映了干旱从开始到结束所持续的时间，干旱烈度反映了干旱历时内的总缺水量，干旱烈度峰值反映了干旱历时内的最大缺水流量，如图3所示。根据以往的研究成果[10-11]，通常选取Q70～Q95(即日流量70%～95%的分位数)作为阈值，其中日阈值每年有365个(闰年有366个)值，它们可从所要研究年份的全部日流量数据中提取。

图3 径流干旱事件的融合剔除

对于逐日流量序列，在相邻的两个历时较长、烈度较大的径流干旱事件之间，经常会出现流量在短时间内超越设定阈值的情况。在这种情况下，通常认为这两个相邻干旱事件之间存在联系，应将它们融合为一个干旱事件，如图3所示。TALLAKSEN L M等[12]通过判断相邻干旱间隔时间与超过设定阈值流量的情况来对它们进行融合。假设两个相邻的径流干旱事件为Ri与Ri+1，它们的干旱历时、干旱烈度与干旱烈度峰值分别为Di、Si、Pi和Di+1、Si+1、Pi+1，若它们满足式(1)，则需要对它们进行融合。式(1)为:

(1)

设融合后的径流干旱事件为Rp，融合方法如式(2)所示:

(2)

式中Dp、Sp和Pp分别为Rp的干旱历时、干旱烈度和干旱烈度峰值。

若融合后的径流干旱事件与下一个相邻径流干旱事件之间超过设定阈值的流量过程的历时与流量仍同时满足上述判别标准，则继续进行融合，直至相邻径流干旱事件的间隔时间或超过设定阈值流量不满足上述判别标准为止。图3中左边的3个干旱事件需经过2次融合。

进行干旱事件融合之后，一般仍存在大量干旱历时较短或干旱烈度较小的干旱事件，这些小干旱事件对干旱特征变量的统计分析意义不大，但会影响分析结果的客观性，因此需要将它们剔除，如图3中第4个干旱事件。设径流干旱事件系列的干旱历时均值为E{D}、干旱烈度均值为E{S}，对于某一个径流干旱事件，若其历时小于rDE{D}或烈度小于rSE{S}，则该径流干旱事件被剔除，通常rD=rS=0.3[12]。

2.2 干旱特征变量边缘分布

首先，选择指数分布、对数正态分布、威布尔分布、伽马分布、广义帕累托分布和广义极值分布分别拟合干旱特征变量的边缘分布，并采用极大似然法进行参数估计。然后，经拟合度检验，确定对数正态分布作为干旱历时(D)、干旱烈度(S)的理论边缘分布，伽马分布作为干旱烈度峰值(P)的理论边缘分布。

对数正态分布函数为:

(3)

式中：X为随机变量，x为其样本值；μ、σ分别为期望和标准差；erf(·)为误差函数。

伽马分布函数为:

(4)

式中：α、β分别为形状参数和尺度参数；Γ(·)为伽马函数。

采用Kolmogorov-Smirnov(K-S)检验、Cramer-von Mises(C-M)检验和Anderson-Darling(A-D)检验，进行上述边缘分布的拟合度检验，若检验统计量小于临界检验值，说明通过检验。

2.3 干旱特征多变量联合分布模型

Copula函数可以联结任意类型的边缘分布函数在[0,1]区间上形成服从均匀分布的多元分布函数[13]。假设F1、F2、…、Fn为连续随机变量的边缘分布函数，则对于任意的x∈Rn，多元分布函数F(x1,x2,…,xn)与n维Copula函数C满足以下关系：

F(x1,x2,…,xn)=P{X1≤x1,X2≤x2,…,Xn≤xn}

=C[F1(x1),F2(x2),…,Fn(xn)]

=C(u1,u2,…,un) 。

(5)

式中：x1、x2、…、xn分别为观测样本；ui=Fi(xi),为边缘分布函数。

水文领域较为常用的Copula函数有Meta-elliptical Copula(椭圆型)函数和Archimedean Copula(阿基米德型)函数，其中,Archimedean Copula函数的计算简便、结构简单，且构成的多变量联合分布函数的形式多样。因此，二维Archimedean Copula函数在水文学领域得到广泛应用。多维Archimedean Copula(n≥3)函数一般适用于描述变量之间具有正相关的情况，但用它们来描述变量间相关性结构时存在诸多的局限和不足。Meta-elliptical Copula函数虽结构相对复杂，但在用于模拟多维随机变量时，不受变量间相关性范围或相关性结构类型等的限制[14]。

采用Meta-elliptical Copula函数中的Gaussian Copula函数以及Archimedean Copula函数中的GH Copula函数与Frank Copula函数来构建两变量和三变量联合分布模型，利用极大似然估计法进行参数估计，并运用赤池信息量准则(Akaike Information Criterion，AIC)、贝叶斯信息准则(Bayesian Information Criterion，BIC)和均方根误差准则(Root Mean Square Error Criterion，RMSE)进行拟合优度评价，以AIC、BIC和RMSE值最小作为拟合程度最优的判别依据,相关计算公式如下：

(6)

式中：MSE为均方根误差(Mean Square Error,MSE)；n为样本容量；k为Copula函数参数个数；Pei为联合分布的经验频率；Pi为联合分布的理论概率。

2.4 重现期计算

单变量重现期的计算公式为：

(7)

设干旱历时(D)、干旱烈度(S)和干旱烈度峰值(P)的边缘分布函数分别为u=FD(D)，v=FS(S)，w=FP(P)。D、S的同现重现期(Ta)和联合重现期(To)的计算公式分别如式(8)和式(9)所示：

(8)

(9)

式中d、s分别为某一单变量重现期下的干旱历时和干旱烈度值。

D、P和S、P的情况类似，受篇幅限制，此处不再列出。D、S、P的同现重现期(Ta)和联合重现期(To)的计算公式分别如式(10)和式(11)所示：

(10)

(11)

式中p为某一单变量重现期下的干旱烈度峰值。

3 结果与讨论

3.1 径流干旱事件的提取

为研究不同日阈值对径流干旱事件识别的影响，并确定最适合用于识别唐乃亥站径流干旱事件的阈值，文中选取6个日阈值，分别为Q70、Q75、Q80、Q85、Q90和Q95来提取干旱历时(D)、干旱烈度(S)与干旱烈度峰值(P)，并对相关结果进行统计分析，具体见表1。

表1 不同日阈值时的径流干旱事件统计分析

由表1可以看出，随着日阈值逐渐减小(从Q70到Q95)，干旱次数、短历时干旱次数、平均干旱历时和平均干旱烈度大致呈下降趋势。然而，阈值过小时，如日阈值为Q95，虽然短历时干旱所占比例依旧较大，但平均干旱历时和平均干旱烈度都比日阈值为Q90的要大，说明日阈值过小会使相邻短历时干旱事件交织连接，导致历时长、烈度大的干旱事件大幅增加；除此以外，虽然短历时干旱次数随日阈值的减小而降低，但是短历时干旱所占比例整体呈上升趋势，说明用小日阈值来识别径流干旱事件时，短历时干旱所占比例增加，这些短历时干旱对干旱特征变量的统计分析意义不大，反而会给计算带来干扰。所以，在选择日阈值时要尽量避免选择太小的日阈值。从表1中还可以看出，日阈值为Q80时，短历时干旱所占比例最小，所以下文重点分析Q80作为日阈值时所识别的唐乃亥站的径流干旱事件。

对初步识别的径流干旱事件进一步融合剔除，并对该过程中的干旱特征变量进行统计分析，结果见表2。

表2 干旱事件的融合剔除过程分析

从表2中可以看出：在1956—2018年间，唐乃亥站共发生429次干旱，即平均1.76个月发生1次干旱；融合剔除之后，该站共发生93次干旱，即平均8.13个月发生1次干旱，平均干旱历时从10.3 d增加至41.2 d。可以认为，经过融合剔除之后，径流干旱事件更符合干旱发生的规律和干旱持续性的特点。因此，非常有必要对识别的径流干旱事件进行融合剔除处理。

3.2 干旱特征单变量边缘分布

经计算，干旱历时(D)的对数正态分布参数σ=0.89、μ=3.29；干旱烈度(S)的对数正态分布参数σ=1.11、μ=0.31；干旱烈度峰值(P)的伽马分布参数α=2.33、β=69.11。采用K-S检验、C-M检验和A-D检验进行拟合度检验，结果见表3。

从表3中可以看出，α分别为0.01、0.05、0.10时，检验统计量均小于临界检验值，表明唐乃亥站干旱特征变量D、S、P的理论边缘分布均通过了拟合度检验。

表3 干旱特征变量边缘分布拟合检验

3.3 干旱特征多变量联合分布

采用Copula函数描述干旱特征变量的相关结构之前需进行相关性分析。本文利用Pearson相关系数γn、Spearman秩相关系数ρn和Kendall秩相关系数τn分别度量其相关性。干旱历时(D)、干旱烈度(S)和干旱烈度峰值(P)3个特征变量两两间相关系数的计算结果见表4。

从表4中可以看出，S-P的相关程度最高，D-S的次之，D-P的最低，即D、S、P之间存在一定相关性，可利用Copula函数构建干旱特征的多变量联合分布模型。

表4 干旱特征变量间的相关性计算结果

3.3.1 干旱特征变量二维联合分布

二维Copula函数的参数估计和拟合优度评价结果见表5。干旱特征变量二维联合分布如图4所示。

表5 二维copula函数参数估计和拟合优度评价结果

图4 干旱特征变量二维联合分布的P-P图

由表5可知，Gaussian Copula、GH Copula、Frank Copula函数分别是拟合D-S、D-P、S-P二维联合分布的最优模型。

图4中将相应二维联合分布的理论概率与经验频率进行了比较。从图4中可以看出，D-S、D-P和S-P的经验点据大致分布在45°线附近，表明所采用的最优Copula函数模型均能很好地拟合相应的经验点据。需要指出的是，尽管Frank Copula函数拟合S-P二维联合分布的效果最好，但在进行同现重现期分析计算时发现其结果并不合理。故文中采用AIC、BIC和RMSE值次小的Gaussian Copula函数来拟合S-P二维联合分布，其理论概率与经验频率的对比结果如图4(d)所示。从图4(d)中可以看出，经验点据大致分布在45°线附近，且与图4(c)中Frank Copula函数的拟合结果相差不大。

3.3.2 干旱特征变量三维联合分布

三维Copula函数的参数估计和拟合优度评价结果见表6。将相应三维联合分布的理论概率与经验频率进行比较，结果如图5所示。

表6 三维copula函数参数估计和拟合优度评价结果

从表6中可以看出，Gaussian Copula函数是拟合干旱历时(D)、干旱烈度(S)和干旱烈度峰值(P)三维联合分布的最优模型。

图5 干旱特征变量D-S-P三维联合分布的P-P图

从图5中可以看出，D-S-P的经验点据大致分布在45°线附近，表明所采用的最优Copula函数模型能够很好地拟合相应的经验点据。

3.4 单变量与多变量的干旱事件重现期

在单变量干旱事件的重现期分别为2、5、10、20、50和100 a条件下，根据式(7)可分别反求相应干旱烈度(S)和干旱烈度峰值(P)，根据式(8)和式(9)可得采用Frank Copula函数拟合S-P下的同现重现期(Ta)与联合重现期(To)，结果见表7。

表7 S、P边缘分位数和基于Frank Copula函数得到的联合分布重现期

计算D-S、D-P、S-P的To与Ta，并据此绘制干旱特征两变量Ta与To的等值线图，如图6所示。

图6 干旱特征两变量同现重现期与联合重现期等值线图

从图6中可以看出，相同条件下Ta大于To，表明干旱特征两变量同时大于等于某一特定值的重现期比其中任意一个大于等于某一特定值的重现期长。

由式(8)至式(11)可计算得到不同单变量重现期条件下干旱特征单变量分位数及两变量与三变量同现重现期(Ta)与联合重现期(To)，结果见表8。

表8 干旱特征变量边缘分布及其对应的联合分布重现期

由表8可以看出，单变量重现期都介于两变量或三变量重现期Ta与To之间，即可以把多变量重现期Ta与To分别视为单变量重现期的上、下临界值，可利用它们估计实际干旱重现期的变化区间，对干旱预测有一定参考价值。同时，也可将它们作为单变量重现期分析的补充，以便更加全面地反映干旱的真实特征，有利于对实际旱情的分析。此外，单变量重现期对应的Ta大于To，且Ta与To之间的差距会随着单变量重现期的增大变得越来越大。例如，单变量重现期为2 a时，D≥38.73 d∩S≥2.17亿m3的同现重现期为3.13 a，D≥38.73 d∪S≥2.17亿m3的联合重现期为1.47 a，而单变量重现期为100 a时，D≥239.81 d∩S≥21.29亿m3的同现重现期为464.23 a，D≥239.81 d∪S≥21.29亿m3的联合重现期为56.04 a。

同时可以发现，若采用Frank Copula函数作为S-P的联合分布模型，单变量重现期为100 a时，S-P对应的同现重现期为2 060.04 a(见表7)；采用Gaussian Copula函数作为D-S的联合分布模型，单变量重现期为100 a时，D-S对应的同现重现期为464.23 a(见表8)。由表4可知，S-P的相关程度比D-S的高，根据以往研究结果，相关性高表明某一干旱特征值出现时，另一干旱特征出现的概率也高。因而，在单变量重现期相同的情况下，S-P的同现重现期应小于D-S的同现重现期。鉴于由Frank Copula函数所得结果与此矛盾，故它不适合作为S-P的重现期的最佳拟合模型。而运用Gaussian Copula函数作为S-P的联合分布模型，单变量重现期为100 a时，S-P对应的同现重现期为288.11 a(见表8)，小于D-S对应的同现重现期464.23 a，所得结果较Frank Copula函数的更为合理。

4 结语

本文利用黄河源区唐乃亥站1956—2018年的逐日流量资料，基于游程理论、可变日阈值和径流干旱事件融合剔除，获得了干旱历时(D)、干旱烈度(S)、干旱烈度峰值(P)的时间序列，采用Copula函数构建了D、S、P之间的二维和三维联合概率分布，计算并分析了相应单变量和多变量的联合重现期，主要研究成果如下：

1)短历时干旱所占比例随着日阈值的减小总体呈上升趋势，在选择日阈值时要尽量避免出现过多短历时干旱事件。

2)经过融合剔除后的径流干旱事件更符合干旱发生的规律和干旱持续性的特点，所以，非常有必要对识别的径流干旱事件进行融合剔除处理。

3)Gaussian Copula函数是拟合唐乃亥站径流干旱D-S、S-P和D-S-P的最优分布模型，GH Copula函数是拟合D-P的最优分布模型。虽然用Frank Copula函数拟合S-P时AIC、BIC和RMSE值最小，但是存在单变量重现期较大时相应同现重现期不合理的问题。因此，在选择拟合性最好的Copula函数时，应将AIC、BIC和RMSE值最小的Copula函数作为首选，同时也应根据实际情况多方面考虑。

4)相比单变量重现期，径流干旱多变量联合重现期(To)更短。换言之，在同时考虑干旱事件多个特征属性的情况下，出现某一个或两个超过设定阈值干旱特征的可能性将增大、期望时间间隔缩短。文中计算结果表明，多变量视角下黄河唐乃亥区域径流干旱的总体程度和风险均较高，多变量视角下的分析能够有效弥补单变量干旱频率分析的缺陷与不足。

5)径流干旱多变量的同现重现期(Ta)与联合重现期(To)可以分别看作单变量重现期的上、下临界值，可据此开展干旱事件概率预测研究。