基于新型综合干旱指数的珠江流域未来干旱变化特征研究

李 军，吴旭树，王兆礼，2，尹家波，陈晓宏

（1.华南理工大学 土木与交通学院 亚热带建筑科学国家重点实验室，广东 广州 510640；2.广东省水利工程安全与绿色水利工程技术研究中心，广东 广州 510641；3.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072；4.中山大学 水资源与环境研究中心，广东 广州 510275）

1 研究背景

干旱是一种自然灾害，具有发生频率高、影响范围广和危害程度深等特点，易对社会经济发展造成不良影响［1-3］。近10年来，中国经历频繁的干旱灾害，造成的损失占气象灾害的50 %以上［4］。随着全球气候变暖，干旱影响将逐渐增大［5-6］，但往往却难以捕捉干旱的开始、持续时间和结束等特征［7-8］。为此，前人做出许多努力来研究干旱［1，9-11］。其中，干旱指数构建是干旱研究的基础，其可对干旱特征进行客观的时空分析和监测，从而达到度量、对比和综合的研究目的［12-14］。故干旱指数的研究具有重要的现实意义。

干旱可分别从气象、水文、农业等角度进行描述和刻画，相应地有多种干旱指数［13］。例如，用于描述气象干旱的标准化降水蒸散指数（SPEI）［15］，用来描述水文干旱的地表水供应指数（SWSI）［16］和标准化径流指数（SRI）［17］，以及用于描述农业干旱的标准化土壤水指数（SMDI）［18］、作物水分指数（CMI）［19］和修正帕默尔干旱指数（scPDSI）［20］等；然而单一指数往往不能综合多种水文气象要素，从而无法综合地从气象、水文、农业等角度全方面地监测干旱，且干旱一般在气象、水文、农业等角度无清晰的界限。因此，需要探索一种综合考虑多种水文气象要素的干旱指数。近年来，为了创建一种融合多种干旱相关水文气象变量的综合干旱指数，国内学者进行了诸多研究。如任怡等［21］采用熵权法和模糊综合法综合PDSI、SPI 和SPEI 指数，建立一种模糊综合评价指数。常文娟等［22］采用主成分分析法将降雨、径流及土壤含水量等要素融合为一种综合干旱指标。然而，模糊综合干旱指数的赋权存在一定主观性，不能有效地客观比较跨时空尺度干旱。主成分分析法的综合干旱指数往往基于降维后的第一主成分构建，不能全面地保持原始变量全部信息，且无法反映水文气象要素信息的非线性影响特征［9］。Copula 函数能构造不同边缘分布的多变量联合分布函数，能巧妙地避免上述问题，且在构建综合干旱指数方面已有较广应用［23-25］。张迎等［23］基于Copula 函数，联合降雨和径流要素，构建了一种能综合表征渭河流域气象和水文干旱的综合干旱指数。Hao 等［24］基于降水、径流两要素，构建美国加利福尼亚州和北卡莱罗纳州的综合干旱指数，并将该指数与标准化降水指数及标准化径流指数对比，论证该综合干旱指数在该地区具有一定的可靠性，具备监测水文、气象干旱的能力。已有的综合干旱指数一般仅考虑径流、降水等要素［23-25］，忽略了蒸散发和土壤水等要素，故难以从水文、气象、农业等多角度综合地刻画干旱。

鉴于此，本文拟基于降雨、蒸散发、径流和土壤水等水文气象要素，运用Copula 函数，构建适宜湿润地区的一种新型综合干旱指数CSDI（Comprehensive Standardized Drought Index），综合表征气象-水文-农业干旱，并以珠江流域为例，分析未来气候变化情景下综合干旱的变化特征。该指数能基于水文，气象，农业等多角度全面地描述干旱，提高综合干旱指数的准确性和适用性，以期为研究区防旱抗旱提供新的参考依据。

2 研究内容和方法

2.1 研究区概况及数据资料珠江流域（图1）面积45.4 万km2，由西江、北江、东江等诸流域组成。珠江流域地处热带和亚热带季风气候区，多年平均气温15～22°C，多年平均降水量1500 mm，年内降水80%集中在4—9月；同时珠江三角洲地区也是中国经济最发达地区之一。近年来，珠江流域干旱现象不断出现，尤其是2003年以来几乎年年出现干旱（2009、2010 和2011年大干旱事件），连续不断的干旱使珠江流域遭受了严重的经济社会损失；干旱灾害频次约占总自然灾害的1/5，均为各项灾害之首；粮食因旱灾减产占总产量的5%以上［26-27］。在气候变化的影响下，珠江流域干旱风险将进一步增加［28］，因此分析研究区未来综合干旱特征变化研究具有重要的意义。

图1 珠江流域地理位置及站点分布

为收集构建综合干旱指数的基础数据，本文基于VIC 水文模型模拟珠江流域蒸散发、土壤水和径流等水文要素。为运行VIC 水文模型，收集了珠江流域及其周边91 个气象站1960—2015年的气象数据，包括降水、气温、风速、相对湿度、气压和日照时数等。水文数据包括高要、石角、博罗水文站点（西江、北江和东江的控制站点）的1960—2000年的逐月流量数据，来源于广东省水文局。地表高程数据来自于90 m 分辨率的STRM 数据。土壤参数来自于FAO 的全球1 km 土壤数据集。植被参数来自于中国科学院1 km 分辨率数据。因CanESM2，CNRM-CM5，CSIRO，MIROC-ESM，MPI-ESM-LR，MRI-CGCM3 等气候模式对中国气候因子均有较好的模拟效果［29］，因此，本文收集这些模式RCP2.6、RCP4.5 和RCP8.5 情景下（2050—2100年）的降水、气温、气压、风速、相对湿度、入射短波及长波辐射等数据。采用Delta 方法［30］对气候模式数据进行降尺度，后用贝叶斯模型平均法进行多模式集成。Li 等［28］研究已对珠江流域VIC 模型进行率定与验证，表明该模型能较为准确地模拟流域的水文过程，可用于流域的干旱评估。同时对多模式集成结果进行评估，表明能较合理地再现珠江流域的气候特征，适合用于珠江流域气候变化的研究（VIC 模型率定与验证及气候模式评估详见Li 等［28］）。VIC 模型在世界范围内已被广泛用于模拟土壤水分和径流，这些研究证明VIC 模型的土壤水数据是可靠且合理的［31-33］。此外，该模型已成功应用于珠江流域土壤水和径流模拟［34-35］。

2.2 综合干旱指数构建现有干旱指数一般分别从农业、水文、气象等角度描述干旱［1，13］，但影响干旱的水文气象因素复杂多样，很难使用一个干旱指数从多种角度完全描述干旱［23］，需要探索能综合从多种角度（如水文、气象、农业等）刻画干旱的综合指数。Copula 函数不受单变量服从何种边缘分布的限制，可客观反映多种水文要素低于某一阈值时的状态，构造不同边缘分布下的多源水文气象变量之间的联合分布。故本文运用Copula 函数构建综合干旱指数。在Copula 函数中，对于多维的随机变量X1，…，Xd，它们的联合分布函数为H(x1，…，xd)，边缘分布函数为F(x1)，…，F(xd)，则一定存在一个Copula 函数C，满足H(x1，…，xd)=C(F(x1)，…，F(xd))［9，23］，该函数将不同边缘分布的相关变量连接组合，不是基于各变量的权重来构建指数［9，23］。为能同时表征气象、水文及农业干旱，本文基于Copula 函数联合降水、蒸散发、土壤水及径流等水文气象要素，构建综合干旱指数（CSDI）。所有要素都可在预定的滑动窗口（例如，3月尺度）累积，以表示不同时间尺度的水分亏缺状况，这使得CSDI 在时间尺度上具有灵活性。前人研究表明3月时间尺度的干旱指数可较好表征中国干旱［36］，故本文以构建3月时间尺度的CSDI 为例。而其他时间尺度的（1、6、9 和12）CSDI 指数具体构建方法跟3月时间尺度基本一样，但降水、蒸散发、土壤水及径流等因子要考虑对应时间尺度的具体值，如6月时间尺度CSDI 指数，相应的水文气象因子要考虑前期6 个月的累积变化。

2.2.1 边缘分布函数 边缘分布考虑广义极值分布（GEV）、皮尔逊Ⅲ型（P-Ⅲ）、韦布尔分布（WBL）、对数正态分布（LOG）、伽马分布（GAM）、和对数逻蒂斯迪克分布（LOL）共6 种常用模型。边缘分布采用Kolmogorov-Smirnov（K-S）方法和均方根误差法（RMSE）进行拟合优度检验，确定最适边缘分布函数。K-S 和RMSE 结果如表1所示，结果表明降水和土壤水的最优模型均为GEV 模型，蒸散发和土壤水则分别为P-III 和WBL 模型。

表1 边缘分布模型拟合优选

2.2.2 最优Copula 函数选择 Copula 函数构造多样，例如Archimedean、Meta-elliptic、Extreme Value 和Miscellaneous 等。Archimedean Copulas 结构简单、代表性广，已在珠江流域水文气象多变量联合分布模拟中取得了良好效果［37］。由于其对高维问题有一定受限，于是更为灵活的嵌套Copulas 应运而生。嵌套Archimedean Copula 函数将二元Archimedean Copula 函数扩展到多维问题框架，用来研究多维变量的相关性，简单和实用［9］。故用嵌套Archimedean Copulas，选取其家族的Gumbel、Clayton、AliMikhail-Haq（AMH）和Frank Copula 函数作为备选函数来选择最优Copula 函数。

利用拟合得到的4 种Copula 函数与经验Copula 的赤池信息量准则（AIC）、RMSE 和相关系数（R）来进行拟合优度评价。AIC 是衡量统计模型拟合优良性的一种标准；RMSE表示模拟值与真实值偏差的平方与模拟次数比值的平方根；R表征两变量之间的线性相关程度；RMSE和AIC值越小，表示模拟拟合效果越优，而R值越大，表示模型效果越优。AIC、RMSE和R的结果如表2所示，表明Archimedean Copulas 的AMH、Clayton、Frank 和Gumbel 模型，均能较好地模拟珠江流域的降水、蒸散发、土壤水和径流等联合分布。而珠江流域的水文气象多变量联合分布的最优模型为Frank Copula 模型。因此，本文选用Frank Copula 函数作为连接函数。

表2 Copula 函数联合分布模型拟合优选

2.2.3 标准化指数 以降水、蒸散发、径流和土壤水序列为随机变量，进行多变量联合分布函数构建，X1，X2，X3和X4为降水、蒸散发、径流和土壤水的某一数字，假设随机变量相应的边缘分布分为F(x1)，F(x2)，F(x3)和F(x4)，则他们的联合分布P可用累积联合概率q和Copula 函数C表示：

上述得到累积联合概率是由所有边际函数得到联合水分亏缺状况，累积概率p值较小表征干燥状况，反之，则表示湿润状况。而累积概率值服从给定的边缘分布，对于不同季节或地区的值，它们通常反映不同水分亏缺状况，不能用于区域水分亏缺比较［38］。因此，用Kendall 分布函数进行转化，可使CSDI 指数用于区域不同干旱的比较。Kendall 分布函数为：

于是得到综合干旱指数：

式中：φ为标准正态分布。

如表3所示，基于CSDI 累积概率值，按照综合严重程度，将干旱等级划分为无旱、轻旱、中旱、重旱和特旱等不同的等级。

表3 CSDI 指数干旱等级划分

2.3 其他方法基于游程理论［37］识别和提取综合干旱事件和特点，包括历时、严重性、峰值等典型特征。为反应CSDI 表征干旱的能力，scPD⁃SI 指数，3月时间尺度SPEI 指数和SRI 指数等作为参考指数，scPDSI 指数计算基于国标GBT 20481—2017 气象干旱等级标准。此外，运用Copula 函数，基于同频放大法，反算给定联合重现期对应的设计组合值。有

式中：干旱历时X、严重性Y和峰值Z为连续随机变量，各自的边缘分布函数分别为Fx(x)、Fy(y)和Fz(z)；P1、P2和P3分别为干旱历时、严重性和峰值的超阈概率；C1、C2和C3分别为历时、严重性与峰值的超阈联合概率；C(u，v)、C(u，w)和C(v，w)均为边缘分布为[0 ，1] 区间的Copu⁃la 联结函数。

3 结果

3.1 综合干旱指数评价图2展示了月尺度下珠江流域1961—2015年CSDI 指数与SPEI 指数、SRI 指数、scPDSI 指数的相关系数空间分布图。由图可知，CSDI 指数与SPEI 指数、SRI 指数、scPDSI 指数的相关系数在全流域基本大于0.7（P＜0.01），呈极强相关性。CSDI 指数与SPEI 指数和SRI 指数在珠江流域上游和下游相关性最大，而与scPDSI 在流域上游端部地区相关性最大。

图2 1961—2015年月CSDI 指数与其他干旱指数相关系数空间分布

对珠江流域1961—2015年CSDI 指数、SPEI 指数、SRI 指数、scPDSI 指数序列进行月干旱频次统计。统计了CSDI 指数至少与SPEI 指数、SRI 指数、scPDSI 指数中的一种同时出现干旱的频次占CSDI指数干旱总频次的比例。图3显示了CSDI 指数干旱监测“准确率”。 CSDI 指数在全流域的干旱监测“准确率”基本在80%以上。这些分析表明CSDI 指数监测干旱具有一定的可靠性。

图3 CSDI 指数干旱监测“准确率”

查阅关于珠江流域干旱研究的相关文献发现［27，39］，珠江流域在1963年遭受特大干旱。因此，以珠江流域1963年的干旱为例，比对CSDI 指数、SPEI 指数、SRI 指数和scPDSI 指数对本次干旱的监测能力。图4显示了基于CSDI 指数、SPEI 指数、SRI 指数和scPDSI 指数的珠江流域1963年干旱的空间演变。1963年1月流域上游开始出现干旱。同年3月，气象干旱发展迅速，基本蔓延至整个流域，流域下游局部地区出现严重气象干旱；同时流域上下游地区出现农业和水文干旱，而中部地区呈无农业和水文干旱，因农业干旱和水文干旱对气象干旱的响应存在一定的延迟［23］。CSDI 指数显示干旱空间分布与SPEI 指数相似，表明CSDI 指数捕捉干旱开始的能力与SPEI 指数相当。1963年7月，整个流域遭受特大气象干旱和水文干旱，CSDI 指数空间分布与SPEI 指数及SRI 指数相似，表明CSDI 指数可监测气象干旱和水文干旱。而流域中上游地区出现特大农业干旱，其余地区以严重农业干旱为主，CSDI 指数可为农业干旱提供及时预警。1963年12月，气象干旱和水文干旱基本在珠江流域结束，而仍有农业干旱，CSDI 指数显示相似的干旱空间分布，这表明CSDI 指数捕捉干旱结束的能力与scPDSI指数相当。这些分析表明综合降雨、蒸散发、径流及土壤水等多种水文气象要素的CSDI 指数能够灵敏、有效地捕捉干旱的开始、持续时和结束等重要特征。

图4 1963年珠江流域干旱空间演变

综上分析可知，本文构建的CSDI 综合干旱指数在干旱监测方面具有一定的可靠性和综合性等优势，兼具SPEI 指数、SRI 指数和scPDSI 指数的能力，如SPEI 指数一样较好地捕捉干旱的开始，也能如scPDSI 指数一样较好地监测干旱的持续时间和结束，并监测水文干旱，能够更加综合、全面地监测水文、气象和农业干旱，提前预警，可为干旱的监测和预警预报提供新的工具。这意味着，当CS⁃DI 指数开始低于干旱阈值时，气象干旱可能开始，应制定气象干旱传播导致农业干旱的应对方案，并发布预警；CSDI 指数继续减小，则是气象干旱开始传播，导致农业和水文干旱，根据应对方案采取及时有效的措施；而CSDI 指数持续低值时，表明3 种干旱可能在该区域存在，必须采取积极的抗旱政策。

3.2 综合干旱时空变化分析既往（1961—2015）和未来情景（2050—2100）综合干旱特点（历时、严重性和峰值）。图5展示了未来情景下干旱严重性、历时和峰值较既往时段变化百分比。RCP2.6 情景下，干旱严重性在珠江流域下游明显低于既往时段，减少20%～40%；流域中上游干旱严重性变大，局部地区增加20%～40%。RCP4.5 和RCP8.5 情景下，整个流域的干旱严重性基本变大，上游严重性明显增加（图5（a））。未来情景的干旱历时变化与干旱严重性变化基本一致（图5（b））。RCP2.6 情景干旱严重性在珠江流域平均减少4.2%，而在RCP4.5 和RCP8.5 情景下，则分别平均增加14.8%和15.3%。 RCP2.6 情景下干旱历时在珠江流域平均减少3.5%，而在RCP4.5 和RCP8.5 情景下，则分别平均增加16.8%和11.5%。RCP2.6 和RCP4.5 情景下，干旱峰值空间变化无明显一致性，上中下游地区均有峰值增加或减少的区域；RCP8.5 情景下，流域中上游大部分地区峰值变大，局部地区甚至增加20%～40%，整个流域平均增加5.1%（图5（c））。

图6为未来情景下干旱最大严重性和历时较既往时段的变化。RCP2.6 和RCP 4.5 情景下，最大严重性和最长历时在流域上下游地区较既往时段降低，局部地区最大严重性和最长历时甚至减少40%以上，而流域中部大部分地区最大严重性和最长历时均变大，局部地区甚至增加60%以上。RCP8.5情景下，最大严重性和最长历时在珠江流域大部分地区均呈增加，分别平均增加11.3%和20.9%。

图6 未来情境下干旱最大严重性、最长历时较既往时段变化

3.3 综合干旱特征联合分布变化图7展示了珠江流域既往和未来情景的干旱特征指标之间的相关性。干旱特征指标之间呈现为正相关性。除了干旱历时与峰值相依程度略微小一些，干旱特征指标之间的相关系数基本上均大于0.5（P＜0.05），相关性程度较高。既往和未来情景的干旱特征严重性和历时相关系数差别都较小。总体上来看，严重性与历时、严重性与峰值的相关性要明显高于历时与峰值。

图7 所有格网既往和未来干旱特点相互关系

为构建干旱特征联合分布，故确认严重性、历时和峰值最优拟合模型。严重性和历时的最优拟合模型均为GEV 模型，峰值的最优模型为WBL 模型。此外，为模拟干旱特征的两变量联合分布，基于Archimedean Copulas 家族的Gumbel、Clayton、AMH 和Frank Copula 函数选取最优拟合模型。干旱特征两变量联合分布最优模型均为Frank Copula。干旱特征边缘分布和联合分布最优模型选取方法采用2.2 节所述方法，详细结果此处不再一一赘述。图8展示了未来情景下干旱特征两变量超阈联合重现期100年的同频设计值较既往时段的变化百分比。虽然严重性与历时、历时与峰值和严重性与峰值的相依性特征有所区别，各自的联合分布模型与参数也存在差异，但相同联合重现期所对应的单变量设计值空间分布变化基本一致。RCP2.6 情景下，珠江流域下游地区超阈联合重现期100年的同频设计严重性和历时值均变小，减少20%～60%；中游地区一百年同频设计严重性和历时值均变大，整个流域分别平均减少13.5%和6.5%；珠江流域北部一百年同频设计峰值均变小，减少20%～40%。RCP4.5 情景下，珠江流域上下游两端地区同频设计严重性和历时值减少40%～60%，而中游地区同频设计严重性和历时值均变大，局部区域甚至增加80%以上；中上游地区一百年同频设计峰值均变小，减少20%～40%。RCP8.5 情景下，珠江流域大部分地区一百年同频设计严重性和历时值均变大，整个流域分别平均增加20.5%和22.1%。

图8 未来情景干旱双变量联合分布100年重现期同频设计值较既往时段变化

4 讨论

在RCP8.5 情景下，珠江流域综合干旱严重性、历时、峰值等均呈明显增加现象，特别中上游地区。影响干旱变化的气象因素一般有降水、气温、蒸散发等因素［1］。降水变化异常一直认为是珠江流域干旱变化的主要因素，特别是中上游地区［26-27］。然而，气候变化背景下，气温显著增加可能导致蒸散发明显增加，而主导干旱变化［40］。在高二氧化碳排放情景下（即RCP8.5 情景），珠江流域未来降水未呈明显增加或减少趋势，而由于气温变化导致，蒸散发呈明显增加趋势，是珠江流域未来干旱严重性增加的主要因素［29］。中上游地区喀斯特地貌的广泛分布也可能加剧干旱。广泛分布的喀斯特地貌改变了地表水文循环过程，特别是降雨径流过程。喀斯特地貌中，植被稀疏且地表覆盖土壤较薄，增加了地表水了下渗，减少地表对降水的蓄留能力。在气候异常情况下（如降水不足伴随极端高温），容易导致频繁的干旱［41］。此外，珠江流域干旱的变化也与大气环流变化有密切关系。就既往时段而言，过去40年，中国西南向东北有一条干旱趋势带［41］，而流域西北部显著变干地区正好与这条趋势带重叠。干旱的发生往往与特殊的大气环流型密不可分。大气风向改变导致水汽输送通量的变化，在很大程度影响珠江流域的降水格局。当西太平洋副高系统东退或减弱时，赤道辐合带东风减弱，流域内输入的水汽量不及输出的水汽多，造成珠江流域西北部上空水汽辐散，削弱云贵地区降水形成条件致使局地干旱发生［26，42］。

干旱作为一种对人类社会影响严重的自然灾害，一直受到学术界和有关部门的持续关注和重视。在全球气候变化背景下，水文、气象或农业干旱频发且交叉，形成综合性干旱，而目前对综合干旱指数构建及其规律的认识尚需更多科学、有效研究，对综合干旱的监测和预报能力尚未达到生产实践的要求［23-25，40］。本文基于降水、径流、蒸散发、及土壤水等水文气象要素，构建适宜于珠江流域的综合干旱指数，该指数与SPEI、SRI 及scPDSI 指数对比，比较了它们相关性及对干旱监测能力，表明该指数兼具表征气象、水文和农业干旱的能力。本文提出的综合干旱指数（CS⁃DI），可对现有的干旱指数进行补充，可为综合干旱的监测及预警预报提供可靠、有力的新工具。但本文尚存一些研究局限。CSDI 指数虽从气象、水文、农业角度等描述和刻画干旱，但未考虑社会经济因素。本文仅验证CSDI 在南方典型湿润气候区的适用性，尚未验证该指数在其他气候区的适用性。此外，本文尚未分析由于模型模拟输入数据的误差，引起干旱特征变化的不确定性，如水文模型模拟水文要素之间可能有重叠部分，而造成干旱严重性和历时增加。未来研究将针对这些不足展开工作，以求为有关部门提供更加可靠、适用的综合性干旱指数，为防灾减灾提供科学依据。

5 结论

本文运用珠江流域降水、蒸散发、土壤水和径流等水文气象数据，基于参数嵌套Frank Copula 函数构建了一种新型综合干旱指数（CSDI），并基于该指数分析了珠江流域未来综合干旱的变化特征。主要结论如下：（1）CSDI指数与SPEI指数、SRI指数和scPDSI指数呈现良好的相关性（P＜0.01），CSDI 指数能够如SPEI 指数一样很好地监测到干旱的开始，也能如scPDSI 指数一样很好地监测到干旱的结束，能够综合表征气象、水文和农业三种干旱。（2）RCP2.6 情景下（2050—2100年），综合干旱严重性和历时在珠江流域均变小，分别平均减少4.2%和3.5%。RCP4.5 和8.5 情景下，整个流域的干旱严重性均变大，分别平均增加14.8%和15.3%，而历时则分别平均增加16.8%和11.5%。（3）综合干旱既往和未来情景下的干旱特征两变量之间具有较高的正相关性（P＜0.05）。RCP2.6 情景下，珠江流域下游地区超阈联合重现期100年的同频设计严重性和历时值均变小，分别平均减少13.5%和6.5%。RCP 8.5 情景下，珠江流域大部分地区100年同频设计严重性和历时值均变大，分别平均增加20.5%和22.1%。