水文变异条件下基于Copula函数的非一致性干旱频率分析方法——以鄱阳湖为例

许 斌，袁 喆，孙可可，鄢 波

(1.长江科学院 水资源综合利用研究所，武汉 430010；2.流域水资源与生态环境科学湖北省重点实验室，武汉 430010)

1 研究背景

鄱阳湖流域位于长江中下游南岸，是我国最大的淡水湖泊[1]。鄱阳湖承接赣江、抚河、信江、饶河及修水五大河之水，并与长江保持着天然的联通，湖区平原土壤不仅肥沃，而且复种程度较高，是我国水稻种植的优势区域之一[2]，为我国的粮食生产和粮食安全做出了巨大的贡献[3]。然而近年来，由于气候变化和人类活动的影响，以鄱阳湖蓄水量减少、水位异常偏低为代表的干旱事件日益增多[4]，说明鄱阳湖湖区范围内干旱重现期出现了较为明显的水文变异问题[5]。干旱的频繁出现，不仅会造成湖区取水困难、湿地萎缩，也会对区域粮食生产、社会经济发展和生态环境保护造成严重冲击[6]。

对干旱事件的研究，经历了维度不断扩展的过程，例如最早的干旱事件识别通常是针对区域内干旱的持续时间或者空间上的一维变化[7]。为了克服一维指标无法反映干旱迁移演变规律的问题，根据干旱影响大、范围广等特性而提出的干旱历时、干旱力度、干旱间隔等二维至多维干旱研究取得了较多的研究成果[8-9]。干旱频率分析是评价干旱事件及其影响程度的重要方法，也是干旱风险分析的重要内容[10]；通过干旱频率分析得出的干旱重现期，对干旱条件下的水资源利用非常有价值[11]。Sklar[12]提出的Copula函数，由于其变量分析的灵活性等特性，能克服传统干旱分析方法的一些不足[13]，因此，Copula函数在二维至多维干旱频率分析的研究中得到了充分的应用和发展。Tsakiris等[14]采用Archimedean Copula函数从干旱烈度和干旱范围2个维度对干旱的频率进行了分析；Saghafian等[15]、Ma等[9]、徐翔宇等[7]均对干旱特征分析的三变量Copula函数方法开展了研究。由于Copula函数等方法的不断进步，目前对于干旱事件的刻画，已经能够对时间、空间上连续的干旱变化特征进行描述，而且能够反映干旱的迁移演变规律。

干旱频率分析的本质是从降水、径流等水文气象等要素的随机性入手，以概率的方式对干旱历时、强度等干旱特性进行描述[16]，需要水文气象等时间序列满足一致性的要求。然而，近年来由于气候变化和人类活动的影响，鄱阳湖流域近年来干旱事件出现了较为明显的“非一致性”问题[17]。Copula函数不受相关变量个数和边缘分布函数类型限制等优点，使得其在干旱分析领域具有很好的优势。但Copula函数的发展多集中于不断扩展的干旱特征变量及其频率分析方面，而对变化环境下可能出现的非一致性，考虑并不十分充分。针对变化环境下的非一致性干旱频率分析问题，本文拟结合水文变异诊断系统[18-19]，提出一种基于Copula函数的非一致性干旱频率分析方法，并以鄱阳湖为研究对象，对近年来鄱阳湖出现的干旱情况进行分析，为变化环境下的干旱识别、干旱频率分析提供新的思路，并为鄱阳湖流域旱灾风险分析等工作提供一定的依据。

2 基于Copula函数的非一致性干旱频率分析方法

为了能够适应环境变化对干旱事件分析带来的影响，基于Copula函数的非一致性干旱频率分析方法首先需要对收集到的水文、降水等时间序列进行水文跳跃、趋势变异分析，以判断时间序列是否满足一致性的要求。对于非一致性的时间序列，则需要划分时段以区别不同的干旱形成条件，为下一步干旱频率分析打下基础。基于Copula函数的非一致性干旱频率分析方法流程如图1所示。

图1 基于Copula函数的非一致性干旱频率分析方法流程Fig.1 Flowchart of inconsistent drought frequency analysis based on Copula function

2.1 水文变异诊断系统

水文序列由随机性成分和确定性成分构成，其中确定性成分主要包括趋势、跳跃和周期成分。当水文序列中不包含确定性成分时，则它是平稳的时间序列，即水文序列形成的物理成因相同，在统计上满足一致性要求，水文序列只在均值上下随机波动，其分布形式或分布参数保持不变。当水文序列中含有确定性成分时，表明水文序列形成的物理成因受到影响，其统计规律是非一致的，即时间序列的分布形式或分布参数发生了变化。从统计学的角度，水文序列变异主要是指水文序列的分布形式或(和)分布参数在整个序列时间范围内发生了显著变化[20]。

在水文序列变异诊断方法上，有诊断跳跃成分变异的秩和检验法、有序聚类法等，诊断趋势成分变异的Spearman秩次相关检验法、Kendall秩次相关检验法等。为了从整体上识别和检验时间序列变异及其变异程度，克服单一检验方法有时检验结果不合理、多种检验方法常常检验结果不一致的问题，谢平等[18]提出了用于水文序列变异诊断的水文变异诊断系统，并不断对其进行完善[21-22]。该系统由初步诊断、详细诊断和综合诊断3个部分组成，如图2所示。该系统能够识别与检验时间序列是否发生变异及其变异程度、变异形式，检验指标全面，通过统计试验获取的不同方法权重赋值客观、诊断结果可信度高。因此，本文采用水文变异诊断系统进行水文序列的变异诊断。

图2 水文变异诊断系统流程Fig.2 Flowchart of hydrological alteration diagnosis system

2.2 非一致性时间序列时段划分

以时间序列的变异诊断结果为基础，对非一致性的时间序列进行不同物理成因时段的划分。由于周期成分在年际间相对变化较小，因此，本文主要对跳跃和趋势2种变异形式的时段划分方法进行阐述。

2.2.1 时间序列变异形式为跳跃

当时间序列存在变异且变异形式为跳跃变异时，根据跳跃点的时间(年份)将时间序列进行划分。变异点之前(含变异点)的时段，其时间序列代表没有受到确定性成分影响之前的随机状态，即过去的物理条件下产生的时间序列；变异点之后的时段，其时间序列代表受到环境变化的影响，序列中存在确定性成分的非随机状态，即现状的物理条件下产生的时间序列。

2.2.2 时间序列变异形式为趋势

当时间序列存在变异且变异形式为趋势变异时，将趋势变异视为多级跳跃变异的合成。以变异诊断结果中跳跃变异分析方法的变异点集中程度为依据，近似地把趋势变异看作多级跳跃变异进行处理，即将趋势变异划分为2级及以上的跳跃变异。再参考跳跃变异的分析方法对时间序列进行时段划分。

2.3 基于Copula函数的非一致性干旱频率分析

水文变异条件下基于Copula函数的非一致性干旱频率分析方法的分析过程，可以概况为以下几个步骤。

(1)选取干旱指数：根据不同的基础数据、干旱分析目标等，选用不同的干旱指数进行干旱事件提取，例如降水时间序列可选用标准化降水指数(SPI)、降水距平指数等；考虑了蒸发时间序列的标准化降水蒸发指数(SPEI)，同时还需要前期雨量等时间序列的帕摩尔干旱指数(PDSI)；考虑径流、水位时间序列的标准化径流指数(SRI)、水位指数(SZI)等等。

(2)分时段干旱特征识别：结合选取的干旱指数，根据划分好的不同时段，基于不同时段的干旱判别标准，分别进行不同时段的干旱特征识别。干旱特征变量可以根据研究需要进行选择，包括干旱历时、干旱烈度、干旱范围、干旱间隔时间等。目前最常见的是选取2～3个干旱特征变量，选定合适的干旱识别阈值，对干旱过程、合并以及小干旱事件进行分析和处理，形成干旱特征变量序列。

(3)分时段构建联合分布Copula函数：对分时段的干旱特征变量序列采用合适的分布函数进行拟合，例如分时段采用伽玛分布、广义帕累托分布、对数正态分布等，并采用K-S检验[23]选用拟合优度最好的线型，选择线性矩估计其分布函数的参数[24]作为边缘分布函数。分时段选取适合的Copula函数类型，对Copula函数的参数进行估计和假设检验[25]，并采用拟合优度最好的Copula函数。

(4)变化环境下干旱特征演变规律分析：采用最优Copula函数计算干旱特征变量的联合频率及重现期，对分时段特征频率的重现期及其变化情况进行分析，进而分析变化环境影响下的干旱演变规律，以及对干旱风险产生的影响，从而提出能够适应变化环境的干旱对策。

水文变异条件下基于Copula函数的非一致性干旱频率分析方法，不仅能够为现阶段分析变化环境对干旱演变规律及其造成的影响提供一定的思路，还能为实际干旱政策的制定提供一定的参考。

3 研究区域概况及站点选择

鄱阳湖位于长江中下游南岸，流域跨江西、安徽、浙江、福建、广东和湖南6省，流域面积16.22万km2，其中96.7%的面积位于江西省境内，占长江流域总面积的8.97%[26]。鄱阳湖是我国最大的淡水湖泊，承接赣江、抚河、信江、饶河及修水五大河之水，是一个天然吞吐型、季节性湖泊。近数十年以来，鄱阳湖湖区范围内的干旱事件就出现了较为明显的非一致性特征，包括湖面极端低水位频发、枯期水位出现时间提前等。

鄱阳湖具有一定的调蓄作用，而水位则是其调蓄结果的最终表现，水位的高低可以最直接地表现出鄱阳湖受干旱影响的状态。湖口站是鄱阳湖入长江干流的控制性水文站，其水位是鄱阳湖与长江干流互相影响的最直接体现。因此，本文选取鄱阳湖的湖口站1955—2015年(共61 a)的水位序列，作为干旱分析的时间序列。

4 非一致性干旱频率分析

依据本文提出的变化环境下基于Copula函数的非一致性干旱频率分析方法，以鄱阳湖湖口站水位序列为研究对象，对变化环境下鄱阳湖湖区范围的干旱情况进行分析，以验证本文提出的非一致性干旱频率分析方法的可行性。受文章篇幅限制，本文在进行非一致性干旱频率分析计算的过程中，主要阐述方法的可行性，对干旱特征指数的分布函数选取、Copula函数选择等较为成熟的部分，本文参考相关研究成果进行了适当的简化。

4.1 湖口站水位序列变异分析诊断

在第一信度水平α=0.05，第二信度水平β=0.01的条件下，利用水文变异诊断系统对湖口站的水位序列进行变异诊断。结果如表1所示。

从表1可以看出，湖口站水位序列在2003年出现了跳跃向下的中变异，且为无趋势变异的跳跃变异。跳跃变异结果(图3)显示，湖口站的水位序列均值从变异前(1955—2003年)的12.855 m变化为变异后(2004—2015年)的12.153 m，下降了0.702 m。

图3 湖口站水位序列跳跃变异结果Fig.3 Diagnosis result of water level series alteration at Hukou station

表1 湖口站年均水位序列变异诊断结果Table 1 Diagnosis result of variation in annual average water level series at Hukou station

变异结果表明鄱阳湖湖口站变异形式比较集中，且发生了对湖口站水位影响较大的突发事件。

因此，湖口站的水位序列已经不再满足一致性的要求，不能直接用于干旱指标的计算及干旱频率分析。

4.2 非一致性湖口站水位序列时段划分

以湖口站水位序列的变异诊断结果为基础，根据其跳跃变异诊断的变异年份2003年，将水位序列划分为变异前(1955—2003年)和变异后(2004—2015年)2个时段。从诊断结果可以看出，湖口站水位序列并无趋势变异。因此，变异前的时段代表没有受到确定性跳跃成分影响之前的随机状态，即过去的物理条件下产生的时间序列；变异后的时段，其时间序列代表受到环境变化的影响，序列中存在确定性跳跃成分的非随机状态，即现状的物理条件下产生的时间序列。

4.3 湖口站非一致性水位序列干旱频率分析

4.3.1 干旱指数选取

依据选取的水位时间序列，本文选取水位距平指数(ZA)作为湖口站非一致性水位序列干旱频率分析的干旱指数，即

ZA=(Zi-Z′)/Z′ 。

(1)

式中：Zi为水位序列；Z′为水位序列均值。

湖口站变异前后不同时段ZA指数的计算结果如图4所示。

图4 湖口站水位序列变异前后的ZA指数序列计算结果Fig.4 Calculation result of ZA index series before and after water level series variation at Hukou station

4.3.2 分时段干旱特征识别

干旱特征识别是干旱频率分析的起点和基础，本文选取干旱历时和干旱烈度[27]对湖口站ZA指数进行干旱事件识别。在干旱指标阈值的选取上，本文将鄱阳湖湖相和河相的临界水位(12 m)选取为是否发生干旱的阈值，将鄱阳湖重度干旱时的水位(10 m)选为重度干旱阈值，并内插获得中度干旱的阈值(11 m)。

由于湖口站变异前后水位序列的均值发生了变化，与之相对应的干旱指标阈值R0、R1、R2也有所不同：R0、R1、R2在变异前分别为-0.067、-0.144、-0.222，在变异后分别为-0.013、-0.095、-0.177。依据干旱指标阈值提取出变异前后鄱阳湖ZA指数序列的干旱历时和干旱烈度，如表2所示。

表2 湖口站变异前后干旱历时和干旱烈度计算结果Table 2 Calculation results of drought duration and intensity before and after variation at Hukou station

通过查阅鄱阳湖干旱的相关资料[28-29]，本文识别出的干旱事件具有较好的准确性，可以认为本文的干旱识别结果是合理的。

4.3.3 基于Copula函数的干旱频率计算

考虑到文章篇幅并为简化计算过程，本文参考杨志勇等[25]的研究成果，直接选用伽马分布对变异前后的干旱历时和干旱烈度进行拟合。采用K-S检验对变异前后的干旱历时和干旱烈度序列进行检验，其检验结果如表3所示。

表3 变异前后的干旱历时和干旱烈度序列伽马分布拟合的K-S检验结果Table 3 K-S test results of fitting gamma distribution of drought duration and drought intensity series before and after variation

从表3中的检验结果可以看出，湖口站水位序列变异前后的干旱历时和干旱烈度序列通过了α=0.05的显著性检验(h=0表示服从伽马分布)，其拟合结果如图5所示。

图5 变异前后干旱历时和干旱烈度经验频率曲线Fig.5 Empirical frequency curves of drought duration and drought intensity before and after variation

结合变异前后干旱历时和干旱烈度的分布函数，采用GH Copula函数构造联合分布推求干旱历时和干旱程度的联合重现期和同现重现期的设计值。其计算结果如图6所示。

图6 变异前后干旱历时和干旱烈度联合概率分布Fig.6 Joint probability distribution of drought duration and drought intensity before and after variation

4.3.4 湖口站干旱特征演变规律分析

4.3.4.1 变异前后干旱历时频率演变规律

计算变异前后干旱历时的理论频率如图7所示。从图7可以看出，相较于变异前，同频率变异后的干旱历时整体偏高；且随着频率的降低，干旱历时增加的幅度不断扩大。这说明随着水位变异的影响，同频率的情况下，变异后的湖口站水位序列更容易出现历时较长的干旱事件。

图7 变异前后干旱历时理论频率曲线Fig.7 Theoretical frequency curves of drought duration before and after variation

4.3.4.2 变异前后干旱烈度频率演变规律

计算变异前后干旱烈度的理论频率如图8所示。从图8可以看出，相较于变异前，同频率变异后的干旱烈度整体偏大；且随着频率的降低，干旱烈度增加的幅度不断扩大。说明随着水位变异的影响，同频率的情况下，变异后的湖口站水位序列更容易出现烈度较大的干旱事件。

图8 变异前后干旱烈度理论频率曲线Fig.8 Theoretical frequency curves of drought intensity before and after variation

随着湖口站水位变异的出现，从频率分析的角度，鄱阳湖更容易出现干旱历时更长、干旱烈度更大的干旱事件，分析结果与鄱阳湖近年来的实际干旱事件发生情况相吻合，说明鄱阳湖干旱发生的风险有所增加。

4.3.4.3 变异前后干旱历时与干旱烈度联合频率演变规律

以鄱阳湖2009年干旱事件为研究对象，结合基于Copula函数的干旱历时与干旱烈度联合频率分析结果，2009年干旱历时(3 a)与干旱烈度(2.5)在变异前后的联合频率及重现期如表4所示。

表4 变异前后干旱历时与干旱烈度相应的联合频率及重现期Table 4 Joint frequency and return period of drought duration and drought intensity before and after variation

从表4可知，2009年的干旱事件，变异后比变异前的频率有较大的增长，其重现期从变异前的20 a一遇，演变为变异后的约2.6 a一遇，重现期存在较大幅度的缩减。这表明相同干旱历时和干旱烈度的干旱事件，在变异后将更容易出现，对鄱阳湖区的干旱防控和水资源利用造成较大的影响。变异后的鄱阳湖湖区水位，其针对不同频率的干旱应对措施应尽快作出相应的调整，以适应水文变异带来的挑战。

5 结论与讨论

针对目前干旱频率分析对“非一致性”考虑不足的问题，本文提出了一种基于Copula函数的非一致性干旱频率分析方法，并以鄱阳湖湖口站为研究对象，采用干旱烈度、干旱历时对湖口站水位序列进行干旱识别，并对近年来鄱阳湖出现的干旱情况进行分析，主要结论如下。

(1)提出的基于Copula函数的非一致性干旱频率分析方法主要分为4个步骤：首先是选取干旱指数；其次在水文变异诊断的基础上，分时段确定干旱阈值并提取干旱特征变量；然后确定干旱特征变量的分布函数并采用Copula函数确定其联合分布；最后根据不同时段的频率和重现期计算结果，对水文变异条件下的干旱演变规律进行分析。

(2)以鄱阳湖湖口站水位序列为研究对象，在第一信度水平α=0.05、第二信度水平β=0.01的条件下，湖口站水位在2003年出现了跳跃向下的中变异，水位序列均值从变异前(1955—2003年)的12.855 m变化为变异后(2004—2015年)的12.153 m，下降了0.702 m。

(3)选取水位距平指数作为干旱指数，并结合鄱阳湖水位确定无旱、中旱和重旱的水位分别为12、11、10 m。

(4)频率计算结果表明，相较于变异前，同频率变异后的干旱历时整体偏高，同频率变异后的干旱烈度整体偏大，与鄱阳湖近年来的实际干旱事件发生情况相吻合，说明鄱阳湖干旱发生的风险有所增加。

(5)结合2009年的干旱事件，相同干旱历时和干旱烈度的重现期从变异前的近20 a一遇，演变为变异后的约2.6 a一遇，重现期存在较大幅度的缩减，对鄱阳湖区的干旱防控和水资源利用造成较大影响。

(6)鄱阳湖相同干旱历时和干旱烈度的干旱事件，在变异后将更容易出现，对鄱阳湖区的干旱防控和水资源利用造成较大的影响。鄱阳湖区有关水资源管理部门应尽快对干旱应对措施做出相应的调整，以适应水文变异带来的挑战。

(7)分析结果对鄱阳湖非一致性干旱演变规律具有一定的参考意义，但后期还应结合鄱阳湖湖口站水位影响因素，例如降水、蒸发、长江干流径流条件、湖区湖底高程变化、五河来水等，对鄱阳湖湖口干旱的成因及不同影响因素的影响占比进行进一步的研究。