基于Copula函数的长江流域上中游干旱概率分析

The Drought Probability Analysis in the Upper and Middle Reaches of the Yangtze River Basin Based on Copula Function

摘要: 受气候变化和强人类活动影响，长江流域年际和年内干旱频发，给长江经济带建设带来严峻挑战，干旱概率分析为把握应对极端天气挑战的主动权提供科学依据。本研究基于长江流域上中游气象站点和水文站点的实测数据以及ERA5数据集，通过标准化降水指数和标准化降水蒸散指数识别了气象干旱事件，并运用Copula函数构建气象–水文联合分布模型，系统解析了长江流域上中游地区气象干旱与水文干旱的联合概率分布特征。结果表明：长江流域上中游的气象干旱发生频率较高，干旱概率为25%~40%；气象干旱与水文干旱的联合概率在夏季最为显著，尤其是在岷江、乌江等流域，夏季的气象水文联合干旱概率达到峰值，夏季降水量减少与径流量下降之间相关性最为显著。两种干旱指标在长江流域上中游的表现一致性较好，自2006年后识别出的干旱事件明显增多，SPI和SPEI识别年均干旱事件次数较2006年之前分别增加了16次和17次。本研究揭示了长江流域上中游干旱发生规律，可为流域抗旱减灾策略的科学制定和水资源优化配置提供基础支撑。

Abstract: Due to the impacts of climate change and intense human activities, interannual and seasonal droughts in the Yangtze River basin have become frequent, posing severe challenges to the construction of the Yangtze River Economic Belt. Drought probability analysis provides a scientific support for taking proactive measures to address the challenges of extreme weather. This study identified meteorological drought events based on measured data from meteorological and hydrological stations in the upper and middle reaches of the Yangtze River Basin, as well as the ERA5 dataset, using the standardized precipitation index and standardized precipitation evapotranspiration index. It also employed Copula functions to construct a meteorological-hydrological joint distribution model, systematically analyzing the joint probability distribution characteristics of meteorological and hydrological droughts in the study area. The results indicate that the frequency of meteorological droughts in the sub-basins of the upper and middle reaches of the Yangtze River basin is relatively high, with drought probabilities ranging from 25% to 40%. The joint probability of meteorological and hydrological droughts is most significant in the summer, especially in the Min River and Wu River basins, where the joint drought probability reaches its peak during this season. Additionally, the correlation between reduced summer precipitation and decreased runoff is most pronounced. The two drought indices (SPI and SPEI) exhibit good consistency in the upper and middle Yangtze River basin. Since 2006, the number of identified drought events has significantly increased, with the annual average number of drought events identified by SPI and SPEI increasing by 16 and 17 times, respectively, indicating an intensification of drought conditions. The findings not only contribute to revealing the drought occurrence patterns in the upper and middle reaches of the Yangtze River basin but also provide a hydrological foundation for the scientific formulation of regional drought mitigation strategies and the optimal allocation of water resources.

1. 引言

干旱作为一种全球性频繁发生的自然灾害，对人类社会和自然有着显著的负面影响[1]。2022年，长江流域因长期高温少雨引发了特大干旱，严重危害了当地人民生产生活和粮食安全[2]。根据政府间气候变化专门委员会(IPCC)关于气候变化和土地的特别报告，在全球变暖的背景下，未来干旱事件的规模将更加严重，会出现更多的极端灾害[3]。我国位于亚欧大陆东部，气候复杂，地形多样导致干旱事件频发[4]。因此，全面有效地监测、预测干旱的发生有助于提高旱情的应对能力，对防灾减灾有着重要的意义。

干旱现象根据其成因和影响范围可划分为四个主要类别：由降水异常引发的气象干旱、水资源短缺导致的水文型干旱、土壤水分不足造成的农业型干旱，以及对社会经济活动产生影响的综合型干旱，这些类型之间存在着从自然因素到社会经济影响的递进关系[5]。作为干旱监测的核心量化工具，干旱特征指标在区域干旱研究中具有关键作用。这类指标不仅能够动态反映干旱强度及其时空演变特征，还在构建干旱预警系统、指导防灾减灾决策等方面具有重要的应用价值。在干旱监测研究领域，学界已开发出多个具有代表性的评估指标体系，其中较为典型的有反映水分平衡的帕默尔干旱指数(PDSI)、基于降水概率分布的Z指数、考虑蒸散发影响的标准化降水蒸散指数(SPEI)，以及用于评估地下水状况的标准化水位指数(SWI)等[6] [7]。例如，陶新娥[8]等基于SPI和SPEI指数分析了汉江流域干旱变化特征。

除了气象条件外，干旱事件还有多种特征，基于一种变量的研究很难全面地分析水文事件的变化规律，因此，越来越多的学者使用Copula函数构建多变量之间的关系。Copula函数能够客观地描述多个变量之间的相关性，是进行多元水文分析和预测的有效方法[9]，例如，龙鹏[10]等基于二维Copula分析了抚河流域水文干旱特征。Tang [11]等结合SRI和SPEI，分析了中国西南四个省份干旱的时空演变特征，验证了干旱的相关性和不同干旱指数的适用性。

长江流域作为中国最重要的流域之一，其干旱事件对社会经济和生态环境的影响极为显著。当前干旱研究多聚焦于流域尺度干旱事件的历时–强度联合特征解析或单一指标下的干旱演变规律识别，尚未建立气象–水文多要素耦合的干旱评估体系，特别是缺乏从子流域单元解析长江流域上中游地区气象干旱与水文干旱的联合干旱分布特征。本研究选取标准化降水指数(SPI)和标准化降水蒸散指数(SPEI)作为关键评价指标，构建了基于Copula函数的气象–水文干旱耦合分析模型。通过该模型系统识别气象与水文干旱的协同发生特征，深入探究两者的联合概率分布规律，从而揭示长江流域上中游地区干旱事件的时空演变规律。本研究对进一步认识长江流域上中游气象水文干旱的响应关系有重要意义，为制定有效的抗旱策略和优化水资源管理提供基础支撑。

2. 研究区域和数据

2.1. 研究区域

长江流域(24~36˚N, 90~123˚N)位于中国中部，西起青藏高原，横跨中国地势的三大台阶，流域地貌类型多，成因复杂。流域面积约180万km²，约占全国陆地总面积的18.8%，是我国最为重要的经济区域带之一。

长江流域位于亚热带季风区，雨热同期，年均降水量为1067 mm，年内降水时空分布不均[12]，上游的年降水量小于500 mm，而中下游降雨量充沛，超过1500毫米。如图1所示，本研究将长江流域上中游分为10个子流域。

图1. 长江流域上中游分区和水文站、气象站分布图

2.2. 数据资料

本研究选取了长江流域上中游119个气象站点的逐时降水数据、长江干支流的流域出口水文站点的逐时流量数据以及ERA5数据集的逐时网格降水、气温数据。由于气象站点数据较短，为提高统计分析的可靠性，通过原始降水数据修正ERA5数据，延长气象站数据的序列长度，且修正后数据比原ERA5数据更准确贴合区域实际气象特征。站点分布如图1所示。对于每个子流域，先将各气象站逐时降水数据累加为逐月数据，再通过算术平均法换算成流域月平均降水数据。将流域出口水文站的逐时流量数据转换为逐月径流量数据。

由于每个站点的数据集长度存在差异，为了确保降水与流量数据之间的一致性，本研究统一了不同子流域气象水文数据集的资料长度(表1)。

表1. 长江上中游子流域气象水文数据集的资料长度

流域/区域	水文站	原始资料长度	修正后长度
金沙江流域	向家坝	2001.1~2020.12	1990.1~2022.12
雅砻江流域	桐子林	2001.1~2020.12	1990.1~2022.12
岷江流域	高场	2001.1~2020.12	1990.1~2022.12
嘉陵江流域	北碚	2001.1~2020.12	1990.1~2022.12
乌江流域	武隆	2001.1~2020.12	1990.1~2022.12
长江干流	汉口	2001.1~2020.12	1990.1~2022.12
清江流域	高坝洲	2001.1~2020.12	1990.1~2022.12
洞庭湖流域	城陵矶	2001.1~2020.12	1991.1~2022.12
汉江流域	皇庄	2001.1~2020.12	1990.1~2022.12
鄱阳湖流域	湖口	2005.1~2020.12	2005.1~2022.12

ERA5数据集(https://cds.climate.copernicus.eu/datasets)是ECMMF生产的第五代大气再分析产品，空间分辨率为0.25˚ × 0.25˚，提供了从1950年至今的连续数据，涵盖多种大气变量[13]。本研究以NetCDF格式下载了研究区域内1990~2022年的ERA5逐时“总降水量”和“2 m温度”数据。ERA5逐时气温数据可近似代表研究区域地表气温。基于气象站的实测降水数据修正后的ERA5降水数据，可近似代表研究区域逐时落地降水量。

3. 研究方法

本研究的技术路线主要包含以下计算步骤：1) 基于历史气象数据，分别计算长江上中游子流域(表1)的标准化降水指数(SPI)和标准化降水蒸散指数(SPEI)；2) 利用Copula函数构建气象干旱与水文干旱的联合分布函数，通过最大似然估计法确定Copula函数的最优参数；3) 计算长江干支流气象水文干旱的联合概率分布，对比分析不同子流域的联合分布特征参数，揭示区域干旱演变规律。

3.1. 干旱指标

3.1.1. 标准化降水指数(SPI)

标准化降水指数(SPI)是确定一个地区在一定时间内降水概率的指标，能有效反映干旱季节性和年际变化特征，适用于月尺度及以上的干旱评估[14]。它的输入数据为长期降水序列，计算方便，可用来评价和比较不同地区的干旱强度和范围。SPI的计算分为两步：首先通过Gamma函数拟合给定的降水序列获得其概率分布，再将其转换成标准正态分布。降水序列的概率密度函数定义为：

$f\left(x\right)=\frac{1}{{\beta }^{\gamma }\Gamma \left(\gamma \right)}{x}^{\gamma -1}{\text{e}}^{-x/\beta }$ (1)

式中：x > 0是区域月降水量，γ > 0，β > 0分别为形状参数和尺度参数。

将累计概率$F\left(x\right)$ 近似转换为均值为0且方差为1的标准正态分布结果，即SPI：

$SPI=-\left(t-\frac{c_0+c_{1}t+c_2{t}^2}{1+d_{1}t+d_2{t}^2+d_3{t}^3}\right),t=\sqrt{\ln \left(\frac{1}{F{\left(x\right)}^2}\right)},0<F\left(x\right)\le 0.5$ (2)

$SPI=\left(t-\frac{c_0+c_{1}t+c_2{t}^2}{1+d_{1}t+d_2{t}^2+d_3{t}^3}\right),t=\sqrt{\ln \left(\frac{1}{{\left(1-F\left(x\right)\right)}^2}\right)},0.5<F\left(x\right)<1$ (3)

式中：$c_0=2.515517$ ，$c_1=0.802853$ ，$c_2=0.010328$ ，$d_1=1.432788$ ，$d_2=0.189269$ ，$d_3=0.001308$ 。

3.1.2. 标准化降水蒸散指数(SPEI)

标准化降水蒸散指数(SPEI)作为一种改进的干旱指数，在降水量的基础上，考虑了潜在蒸散发对干旱严重程度的影响[15]。SPEI的计算方法与SPI具有相似的标准化处理框架，其核心计算过程如下：首先构建月尺度降水量与潜在蒸散发的差值序列(D_i)，随后运用三参数对数逻辑分布(Log-Logistic)对该序列进行概率密度函数拟合，最终通过正态标准化转换获得标准化指数序列。D_i序列的概率密度函数定义为：

$f\left(D\right)=\frac{\beta }{\alpha }{\left(\frac{D-\gamma }{\alpha }\right)}^{\beta -1}{\left(1+{\left(\frac{D-\gamma }{\alpha }\right)}^{\beta }\right)}^{-2}$ (4)

式中：α，β，γ分别为尺度、形状和位置参数。本研究采用Thornthwaite方程计算潜在蒸散发，该方法计算简单，仅需月气温数据即可推求出蒸散发。

再按照与SPI相同的方法，把累积分布函数F(D)转换为标准正态分布，即可得到各月的SPEI值。为了反映长江干支流水文干旱变化特征，选用月尺度为1的值进行分析。SPI，SPEI和径流量的干旱等级的划分见表2，径流量按保证率划分干旱情况。

表2. 干旱等级划分

等级	类型	标准化降水指数	标准化降水蒸散指数	类型	径流量
1	无旱	−0.5 < SPI	−0.5 < SPEI	特丰	P ≤ 12.5%
2	轻旱	−1 < SPI ≤ −0.5	−1 < SPEI ≤ −0.5	偏丰	12.5% < P ≤ 37.5%
3	中旱	−1.5 < SPI ≤ −1	−1.5 < SPEI ≤ −1	平水	37.5% < P ≤ 62.5%
4	重旱	−2 < SPI ≤ −1.5	−2 < SPEI ≤ −1.5	偏枯	62.5% < P ≤ 87.5%
5	特旱	SPI ≤ −2	SPEI ≤ −2	特枯	87.5% < P

3.2. 边缘分布函数

建立联合分布函数前需确定每个变量的边缘分布函数。本研究采用广泛使用的P-Ⅲ型曲线描述月径流数据的概率分布特征[16]；由于SPI和SRI有正值和负值，而有些边缘分布函数要求正值随机变量，所以选用Normal和Gumbel分布确定SPI和SRI和概率分布特征，并采用最大似然函数估计参数[17]。

3.3. Copula联合分布函数

Copula函数由Sklar在1959年提出[18]，定义为在[0, 1]上的均匀分布，广泛应用于多元变量的干旱频率分析。二元联合分布表述为：

$F\left(x,y\right)=P\left(X\le x,Y\le y\right)=C\left(F_1\left(x\right),F_2\left(y\right)\right)$ (5)

式中：F₁，F₂表示两个随机变量的边缘分布函数，C为Copula函数，x为SPI或SPEI，y为径流量。

目前常用的Copula函数有四种，本研究选用三种阿基米德型Copula函数(Frank、Gumbel和Clayton)拟合，并采用极大似然法进行参数估计。

3.4. 函数的检验和选择

本研究采用K-S检验法对边缘分布函数和Copula函数进行显著性检验。通过显著性检验的候选函数进一步运用AIC进行优选，其中AIC统计量的最小值对应最优拟合效果。AIC准则法表达式如下：

$AIC=2k+n\ln \left(MSE\right)$ (6)

式中：MSE为均方误差，k为参数个数，n是样本个数。

4. 结果与分析

4.1. 边缘分布函数的确定

使用两种分布函数Normal分布和Gumbel分布，分别拟合子流域的SPI和SPEI值，对各流域边缘分布函数进行K-S检验后，根据AIC准则法选取拟合程度最高的分布函数作为长江干支流的指标边缘分布函数(表3)。

以岷江流域为例，在6月、8月Gumbel函数的AIC值最小，选择其作为最优边缘分布函数，其余月份选择Normal函数。整体上看Normal函数出现次数比Gumbel函数更多，且同一个子流域在同一时段内SPI和SPEI的最优边缘分布函数选择多数是相同的。

表3. 长江干支流1~12月拟合最优边缘分布函数

流域	分布函数	Normal	Gumbel
金沙江流域	SPI	1~2月、4月、9月、12月	3月、5~8月、10~11月
	SPEI	1~2月、4月、9月、12月	3月、5~8月、10~11月
雅砻江流域	SPI	2月、4月、6~12月	1月、3月、5月
	SPEI	1~2月、4月、6~10月、12月	3月、5月、11月
岷江流域	SPI	1~5月、7月、9~12月	6月、8月
	SPEI	1~5月、7月、9~12月	6月、8月
嘉陵江流域	SPI	2~12月	1月
	SPEI	1~12月	无
汉江流域	SPI	1~3月、6~12月	4~5月
	SPEI	1~3月、5~12月	4月
乌江流域	SPI	1~2月、4~8月、10月、12月	3月、9月、11月
	SPEI	1~2月、4~8月、10月、12月	3月、9月、11月
长江干流	SPI	2~3月、7~8月、12月	1月、4~6月、9~11月
	SPEI	1~2月、4月、7~10月、12月	3月、5~6月、11月
清江流域	SPI	1~3月、5~10月、12月	4月、11月
	SPEI	1~3月、5~10月、12月	4月、11月
洞庭湖流域	SPI	1~2月、4~8月、10~12月	3月、9月
	SPEI	1月、3~12月	2月
鄱阳湖流域	SPI	1~2月、5~6月、8~12月	3~4月、7月
	SPEI	1~12月	无

4.2. Copula函数的确定

本研究选用Frank、Gumbel和Clayton三种联合分布函数，构建长江干支流干旱指标与流域出口径流量(Q)间的二维联合分布模型。三种联合分布函数在长江干支流各时间段的p值均大于0.05，即所选联合分布函数通过K-S检验，接着根据AIC最小值原则，选取最佳二维联合分布模型。各流域气象干旱指标与径流量联合分布函数选择结果见表4。以岷江流域为例，SPI-Q联合分布拟合效果最佳的Copula函数在1月、2月、6月为Frank，4月、7月、8月、11月为Gumbel，其余月份为Clayton；SPEI-Q联合分布拟合效果最佳的Copula函数在1月、2月、6月为Frank，4月、7月、8月、11月为Gumbel，其余月份为Clayton。

与边缘分布函数优选结果规律类似，长江干支流SPI-Q和SPEI-Q联合分布函数在相同月份内选择结果多数是相同的，仅出现少数不相符结果，如嘉陵江流域2月SPI-Q最佳联合分布函数为Frank函数、SPEI-Q最佳联合分布函数为Clayton。

表4. 基于AIC准则选择长江干支流气象干旱指标与径流量最佳联合分布函数

流域	分布函数	Frank	Clayton	Gumbel
金沙江流域	SPI-Q	1月、3月、5月、10月、12月	2月、6月、11月	4月、7~9月
	SPEI-Q	1月、3月、5月、10月、12月	2月、6月、11月	4月、7~9月
雅砻江流域	SPI-Q	1~3月、5月、9月、11~12月	无	4月、6~8月、10月
	SPEI-Q	1~3月、5月、9月、11~12月	无	4月、6~8月、10月
岷江流域	SPI-Q	1~2月、6月	3月、7月、10月	4~5月、8~9月、11~12月
	SPEI-Q	1~2月、6月	3月、7月、10月	4~5月、8~9月、11~12月
嘉陵江流域	SPI-Q	1~3月、6~7月	5月	4月、8~12月
	SPEI-Q	1月、3月、6~7月	2月、5月	4月、8~12月
汉江流域	SPI-Q	1月、5月、8月、12月	无	2~3月、6~7月、9~11月
	SPEI-Q	1月、5月、8月、12月	无	2~3月、6~7月、9~11月
乌江流域	SPI-Q	2月、6月、11月	1月、10月	3~5月、7~9月、12月
	SPEI-Q	6月、11月	1~3月、10月	4~5月、7~9月、12月
长江干流	SPI-Q	1~4月、6月、10月	5月、8月	7月、9月、11~12月
	SPEI-Q	1~4月、6月、10月	5月、8月	7月、9月、11~12月
清江流域	SPI-Q	2~3月、5月、7~12月	1月	4月、6月
	SPEI-Q	2~3月、5月、7~12月	1月	4月、6月
洞庭湖流域	SPI-Q	10月	2月、4~7月、9月	1月、3月、8月、11~12月
	SPEI-Q	10月	2月、4~7月、9月	1月、3月、8月、11~12月
鄱阳湖流域	SPI-Q	无	3月、5~6月、9~10月	1~2月、4月、7~8月、11~12月
	SPEI-Q	5月	3月、9~10月	1~2月、4月、6~8月、11~12月

4.3. 干旱概率分析

根据干旱等级划分表，当SPI或SPEI ≤ −0.5时，子流域发生干旱。为保证分析的时段统一，统计整理了长江流域上中游(鄱阳湖流域)从1991年至2022年(从2005年至2022年)的干旱情况。

SPEI共识别出1132次干旱事件，占总时长的32.75%，SPI共识别出1079次干旱事件，占总时长的31.22%；SPI和SPEI识别出的干旱次数相近。SPI识别干旱事件次数最多的是2011年，共有58次，SPEI识别干旱事件次数最多的是2013年，共有62次；清江流域SPI和SPEI识别出的干旱事件最多，分别为129次和133次。长江干支流的干旱发生概率大多在30%~40%之间，表明长江流域上中游的干旱事件发生比较频繁；自2006年后，发生气象干旱的次数明显上升，SPI和SPEI干旱事件年均识别次数较2006年以前分别增加了16次和17次，流域的干旱情况逐渐严重，且SPI和SPEI的干旱发生概率在同一流域相同月份上接近，在干旱识别上具有较好的一致性。

基于所选的Copula函数，进一步分析长江干支流气象水文干旱的联合概率(图2)。图3展示了岷江流域SPI与流域出口水文站径流量(Q)的二维联合分布。在1990年至2022年8月份，岷江流域的气象水文干旱联合概率最大，为26.87%，气象水文联合干旱事件共出现8次，其中有6次出现在2006年以后，最严重的年份为2011年，对应的SPI为−2.23 (特旱)；岷江流域SPEI与径流量干旱联合概率最大也出现在8月，为26.43%，识别出联合干旱事件9次，其中7次出现在2006年以后，最严重的事件也出现在2011年，对应的SPEI为−1.79 (重旱)；同时结果显示，岷江流域在5~8月气象水文干旱联合概率较高，即在夏季岷江流域更容易出现降水量与径流量同枯的情况。

图2. 长江干支流1~12月SPI-Q与SPEI-Q干旱联合概率

对比长江干支流气象水文干旱联合概率，与图3类似，其他流域SPI、SPEI与流域出口径流量都增大时，其联合概率值均出现先慢后快再慢的增加过程。气象水文联合干旱概率越大，表明气象干旱与水文干旱一致性越强，由图2可知，大部分流域呈现出夏季气象水文联合干旱概率偏高(见黄色部分)的情况，同时，各联合概率结果均不超过30%，这也表明气象干旱和水文干旱很难完全同时出现。总体而言，长江流域上中游夏季气象水文干旱相关性更强，其次是春秋季，且流域出现干旱事件的次数呈现增长趋势，说明长江流域上中游的干旱情况将进一步加剧。

图3. 岷江流域1~12月降水指数(SPI)与径流量(Q)的二维联合概率分布图(C)

5. 结论

本研究基于标准化降水指数(SPI)和标准化蒸散降水指数(SPEI)，运用Copula函数构建联合分布模型，对长江流域上中游地区长系列气象干旱指数与径流量的联合分布特征进行了系统性研究。主要结论如下：

1) 月尺度上，金沙江、雅砻江、岷江、嘉陵江、汉江、乌江、清江、洞庭湖、鄱阳湖流域及长江部分干流SPI识别干旱概率最高的月份分别为7月、1月、8月、1月、4月、3月、11月、9月、3月和9月；SPEI识别干旱概率最高的月份分别为10月、3月、8月、4月、4月、9月、11月、2月、9月和3月，干旱概率位于30%~40%，发生概率较高。

2) 气象干旱与水文干旱的联合概率在夏季最为显著，尤其是在岷江、乌江等流域，夏季的气象水文联合干旱概率达到峰值，这表明夏季降水量减少和气温升高与径流量下降之间存在较强的相关性。此外，清江和鄱阳湖流域在春秋季的气象水文联合干旱概率也较高，表明这些流域在非夏季时段同样面临较高的干旱风险。

3) 长江流域上中游两种气象指标表现出较好的一致性，识别出的气象干旱事件次数呈现明显的增长趋势，尤其是在2006年之后，干旱事件的频率和强度均有所增加，表明在气候变化和强人类活动影响下，长江流域的干旱程度会进一步加剧。

本研究深入分析长江中上游干旱规律，有助于提升区域抗旱减灾能力和水资源节约集约利用能力。

基金项目

国家重点研发计划(2021YFC3200303)资助。

NOTES

作者简介：薛凯元(2001-)，男，硕士研究生，主要从事水库抗旱调度研究，Email: [e-mail]

^*通讯作者Email: [e-mail]