滦河流域气象干旱向水文干旱传播特征及风险分析

张 璇，许 杨，郝芳华，郝增超

（1.北京师范大学 水科学研究院城市水循环与海绵城市技术北京市重点实验室，北京 100875；2.北京师范大学 中国绿色发展协同创新中心，北京 100875）

1 研究背景

我国降水时空分布不均，水资源总体短缺，气候变暖和极端气候事件的增加对我国水资源的可持续利用带来了更加严峻的挑战［1-2］。近年来，我国大范围干旱事件频发，如2000年长江流域干旱，2006年华南、川渝地区高温干旱，2014年华北、河南等地特大干旱等，对农业生产、生态环境和社会经济产生严重的影响［3］。据中国气象灾害年鉴统计，2012—2016年，干旱对我国农业生产和社会经济造成的损失超过2000 亿元。研究表明，在未来自然和人类活动驱动的气候变化背景下，水循环要素间的转化关系将更为复杂［4-5］，给干旱灾害的成因机制、预警预防和应对措施研究带了新的挑战。

国内外学者从水循环不同要素演变过程带来的影响角度进行干旱评价，提出四大类指标来评价干旱：以降水指标划分为主的气象干旱，以地表径流和地下水指标划分为主的水文干旱，以土壤水分和作物指标划分为主的农业干旱，以供水和人类需水指标划分为主的社会经济干旱［6-7］。大部分干旱事件的发生起源于降水的偏少（气象干旱）；降水不足时对土壤的水分的补充减少，而随着蒸散发的能力不断加强，导致土壤水分缺失从而胁迫作物和植物的生长（农业干旱）；当气象干旱持续一定时间并扩展到一定范围时，区域河川径流减少，河流与湖泊的水位下降（水文干旱）；最终影响社会经济系统中的水资源供给（社会经济干旱）。通常气象干旱发生发展及结束较快，而农业/水文干旱滞后于气象干旱且持续时间较长［4］。一次完整的干旱过程包括了干旱的发生、发展、持续、缓和、解除五个阶段［8-10］，在整个干旱事件的发展过程中，干旱的发生（onset）、干旱烈度（severity）、干旱的解除（recovery）等是评价干旱的主要特征因子。研究表明，不同类别干旱的发生时间、强度等特征之间存在明显的相关性，水分缺失信号在不同类型干旱之间的传递被定义为干旱的传播［11］。作为水循环过程的重要组成部分，干旱的传播对水资源管理与干旱事件的预警具有非常重要的意义［12］。

研究农业/水文干旱对气象干旱响应时往往可以采用相关分析［13］、游程理论［10］、函数拟合［9］等分析方法，甄别气象干旱向农业/水文干旱的传播特征。吴杰峰等［9］通过建立Logarithm 函数研究晋江流域气象干旱演进为水文干旱的临界条件，结果表明为气象干旱临界历时约为1.45 个月。Xu 等［13］应用相关分析方法评价滦河流域气象干旱向水文干旱的传播时间约是1 ～ 6 个月。刘永佳等［14］基于相关性分析方法讨论了不同季节气象干旱向水文干旱的传播动态变化受降水、气温、煤炭开采和水库建造等因素影响。目前的研究多基于干旱事件的匹配与相关性分析干旱的传播特征，对于不同类型干旱的特征传播规律研究相对较少，干旱发生时间、影响范围以及结束时间等因素的传播关系的研究仍待厘清。

本研究以滦河流域山区（滦县站以上）为研究对象，采用标准化降水指数（Standardized Precipita⁃tion Index，SPI）评价流域的气象干旱，采用标准化径流指数（Standardized Runoff Index，SRI）评价流域的水文干旱，采用游程理论方法从各站点SPI 和SRI 序列中分离出气象干旱与水文干旱的发生次数、干旱历时和干旱峰值强度等特征；其次对气象干旱向水文干旱的传播事件进行筛选和匹配，讨论干旱传播过程中各要素的演变特征和滞后性；在此基础上建立气象干旱和水文干旱特征关系模型，最后构建气象干旱和水文干旱特征的最可能组合，进一步揭示了干旱特征变量的演变特征和传播风险，为干旱的预警提供支撑。

2 研究区域概况及数据

2.1 研究区概况滦河山区地处39.66°N—42.73°N，115.47°E—119.73°E 之间，属于我国典型的半湿润半干旱区的过渡带。流域总面积为43 940 km2，地势起伏明显，海拔高度从最高的2241 m，到平原地区的3 m 左右。根据流域的地势特征，将研究区域进一步划分成为上游和下游两个区域，上游区域面积2.4 万km2，主要集中在内蒙古草原、坝上高原的南部和河北的山区，下游面积2.0 万km2，主要集中在燕山山脉区域。滦河流域是我国北方重要的生态屏障区，不仅有世界闻名的塞罕坝森林公园，而且还是京津冀重要的水源涵养地区，每年向天津和唐山供水超过21 亿m3。滦河流域年平均温度为1 ～ 11℃，多年平均降水为400 ～ 800 mm，且由于地处大陆季风性气候区，降水主要集中在6—8月。流域地表径流的时空分布差异也非常明显，1960年代上下游的年均径流深分别为44.04 mm 和93.03 mm，2000年代分别减少为19.32 mm 和22.68 mm。据年鉴记载，滦河流域地区发生了多次旱灾，如1972年发生历史上罕见的春夏连旱，承德地区受灾面积达1907 km2，约占流域内耕地面积的80%。1980—1984年、1997—2005年也均有重大干旱灾害事件发生［15］。

2.2 数据来源与处理本研究降水数据来源于中国气象共享服务网站（http：//data.cma.cn/）提供的流域内部及周边的9 个气象站点的1960—2017年的月均降水数据。采用标准正态检验（Standard Normal Homogeneity Test，SNHT）方法对降水数据进行均一性检验［16］，结果显示9 个气象站的降水数据均通过均一性检验，数据能够用于流域气象干旱评价。而后采用反距离加权（Inverse Distance Weight，IDW）插值法［17］得到的上、下游中心处降水量作为上、下游地区的降水量。同时从河北省水务局收集到滦河流域干流三道河子站和滦县站、以及滦河上最大支流伊逊河上控制站点韩家营站的水文资料。根据各控制站的控制面积，将水文站监测到的月径流量转换为上、下游平均月径流深。流域地理位置及站点空间分布如图1所示。

图1 研究区地理位置

3 研究方法

3.1 气象干旱与水文干旱评价由于在不同时间尺度的降水量与径流量的变化幅度较大，直接采用降水量和径流量在不同时空尺度上进行相互比较可能存在偏差，因此，分别采用标准化降水指数（SPI）［18］和标准化径流指数（SRI）［19］来表征某一时段内流域气象干旱和水文干旱状况，两个指数均具有计算简便、意义明确等特点，已在多个国家和地区作为业务化干旱监测指标得到广泛的应用［20］。SPI 计算原理是以某月的降水量为随机变量x，将降水序列拟合成一个合适的概率分布函数。与SPI 计算原理相似，SRI 可以通过拟合径流深度序列获得。已有研究已经证实Gamma 分布能较好地拟合流域降水和径流序列［21］，其概率密度函数可表示为：

式中，α＞0，β＞0，α和β分别为的形状和尺度参数，可用极大似然估计法求得。

已有研究表明，3 个月尺度的SPI 和SRI 对年内干旱事件的爆发与持续影响更为敏感，能够较好的反应研究区干旱的季节性特征［22］，同时与流域实际的干旱情况更为接近［13，23］，可以反映短期的气象和水文的干旱特征。因此，本文选用3 个月尺度的SPI 与SRI（即SPI3 和SRI3）分别对气象干旱和水文干旱事件进行识别与分析。

3.2 干旱事件的识别与特征要素的提取基于SPI3 和SRI3 计算结果，采用游程理论方法［24］分别提取干旱事件，即设定阈值，截取离散的时间序列值，当序列值（x）在一个或多个时间段内连续大于阈值时，判定为正游程；反之则判定为负游程。如图2所示，本研究基于阈值法识别干旱过程设定了3 个阈值，分别为X0，X1，X2，将阈值分别设定为0，0.5，-0.5。当SPI3 或SRI3 小于X0（0）且干旱历时大于1 个月（如干旱事件1），此次干旱就被识别为一次干旱事件；当某时刻X大于X2（-0.5）但干旱历时仅为1 个月时，即图中t3-t4时段内的干旱信号则不被识别为干旱事件；同时，当两次干旱事件间隔小于1 个月时，且其间x小于X1（0.5），如图2 中干旱事件2所示，可以将相邻两次干旱事件合并为一次干旱事件，即图中t5-t8时段内的干旱信号被合并为一次干旱事件。干旱事件从开始到结束的时间也就是负游程的长度定义为干旱历时。干旱事件发生过程中的绝对值最大的干旱强度值定义为干旱峰值强度（Peak intensity，PI），通常达到PI 也被认为是干旱从发展进入到解除阶段的分界点［25］。研究期内年均干旱次数定义为干旱发生的频率。从SPI3 或SRI3 序列中识别出干旱的次数、干旱历时、干旱峰值强度以及干旱不同阶段的时间节点（干旱开始、峰值、结束时间）6 个干旱的特征。

图2 干旱事件识别方法

3.3 干旱的传播规律与特征分析由于一个固定的时段内流域的降雨和径流不可能表现出完全对应的关系，并不是每一次气象干旱都会导致水文干旱的发生，且干旱在传播过程中会发生滞后、减缓、扩大、减弱等现象［10］，已有的相关性分析方法［13，26］可能会导致干旱传播识别结果的不确定性。因此，为更加全面的反映气象干旱向水文干旱传播的特征，本研究将识别出的气象干旱事件与水文干旱事件逐一进行识别和匹配，当某一场水文干旱发生前3 个月到该干旱结束时间内发生与之对应的气象干旱事件，则认为该气象干旱传播形成水文干旱，即发生干旱传播事件。通过分析两类干旱特征之前的演变关系，能够系统的甄别从气象干旱到水文干旱的特征传播规律。由于干旱的历时、峰值强度是干旱事件预测和管理的关键特征，本文重点关注干旱历时和峰值强度的传播规律。

3.4 干旱的传播分析评价由于干旱变量之间传播关系复杂，很难用线性关系直接刻画特征变量之间的联系，Copula 函数能够将若干个随机变量的边缘分布连接起来，构造成为联合分布函数。对于边缘分布为的n维函数，存在n维Copula 函数C满足：

本文关注气象干旱和水文干旱特征变量间的关系，对于水文干旱特征的概率分布函数FX（x）和气象干旱变量概率分布函数FY（y）及它们的联合分布函数F（x，y），联合重现期为：

对于Copula 函数的单变量，在两变量联合重现期等值线上包含无数种组合，Salvadori 等［27］定义了最大可能权函数方法，能够基于边缘分布及联合分布函数推出最大可能发生的组合，从而量化风险。本文采用该方法定义两变量联合重现期下水文干旱和气象干旱特征（干旱历时/峰值强度）的组合情景：

最后，通过上述方法构建不同联合重现期下气象干旱和水文干旱特征的最可能组合，建立回归模型，研究表明［5］气象干旱与水文干旱存在线性相关关系，回归模型如下：

式中E和R分别为回归模型的弹性系数和余项。

4 结果与分析

4.1 气象干旱与水文干旱特征运用游程理论对滦河流域气象和水文干旱进行识别，结果表明（图3），滦河流域在1960—2017年的近60年时间内，气象干旱事件的频率远大于水文干旱事件，识别出的气象干旱事件数量约为水文干旱事件的2 ～ 3 倍。其中上游地区气象干旱发生频率（1.26 次/年）要高于下游地区（1.11 次/年）；而水文干旱发生频率一致，约0.39 次/年。气象干旱比水文干旱的历时短，气象干旱的平均历时在上下游地区分别为4.7 个月和5.3 个月，最长的一次气象干旱持续时间为18 个月，最短的为1 个月，而水文干旱平均历时高于15 个月，在1998 至2006年间滦河流域上下游地区均发生持续的水文干旱。干旱峰值强度方面，气象干旱的平均峰值强度为1.3，高于水文干旱的平均峰值强度（1.0），其中下游地区气象干旱和水文干旱峰值强度略高于上游地区。

图3 滦河流域气象干旱与水文干旱特征

4.2 干旱的传播规律当气象干旱发生后水文循环系统中的水量也发生相应变化，自然条件下，随着气象干旱的不断累积，形成水文干旱，而在不同因素的影响下可能会加剧或缓解水文干旱的形成［28］。对时间相邻的气象干旱和水文干旱进行对应关系识别，发现存在以下4 种场景：S0：气象干旱和水文干旱间无明显的响应关系；S1：气象干旱和水文干旱一一对应；S2：多场气象干旱综合作用引起单次水文干旱；S3：多场气象干旱形成多场水文干旱。例如，图4所示的滦河上游地区1971年至1978年干旱事件中，第16、17、18、24 和25 场气象干旱未发生传播，为S0 场景；第15 场气象干旱传播至水文干旱形成第4 场水文干旱事件，即气象干旱和水文干旱事件一一对应关系（S1）；第22、23 场气象干旱共同形成了第7 场水文干旱，为S2 场景；第19 ～ 21 场气象干旱共同形成了第5、6 场水文干旱即图4 中第20 场干旱事件，为S3 场景。

图4 干旱传播事件识别和划分方法

本文基于3.3 节所述方法，对滦河流域气象干旱和水文干旱事件进行匹配，在滦河流域上、下游地区依次识别出40、36 次干旱事件，各种场景的传播关系占比如图5所示，未发生传播的干旱事件（S0）在滦河流域各区域都占有较大比例，在上下游地区分别占60%和59%，以气象干旱没有发生传播的类型为主。一方面，如图6所示，S0 情景下气象和水文干旱的平均历时以及峰值强度均小于其他情景可知，当气象干旱历时较短或峰值强度较弱时往往不引起水文干旱［9，29］，由于出现短期的降水不足时，土壤水或地下水的补给能够维持流域的产流过程，使径流量没有明显减少［30］。另一方面，人类活动对干旱传播关系的干扰，当气象部门预测或监测到气象干旱发生时，通过采取涵闸引水、水库放水、泵站提水等一系列措施补充河道水量，能够维持生态流量［31-32］，使气象干旱不传播为水文干旱。在识别出的干旱传播事件中，气象干旱和水文干旱事件一一对应的干旱事件（S1）占比最大，在上、下游地区分别占20%和22%。比较S1 情景下的干旱特征（图6）发现此情景下水文干旱比气象干旱的历时短、峰值低，即当一次独立的气象干旱事件传递形成一次水文干旱的过程中，干旱峰值强度会被削弱，而相应的干旱时间也较短。这一过程也符合流域降水产流过程的基本特征。多次气象干旱形成一次水文干旱的情景（S2）也较为常见，上、下游地区分别有12%、11%的干旱事件为此类情景，此时，水文干旱比气象干旱的历时短、峰值高。当气象干旱发生后水文干旱还未恢复，又与新的气象干旱叠加，容易形成持续时间长且难以恢复的水文干旱事件。多次气象干旱引起多次水文干旱的组合情景（S3）发生次数较少，在上、下游地区均占8%，这类干旱事件形成过程复杂，在所有类型中平均历时最长，峰值强度最大，发生后会造成重大影响。

图5 滦河流域上游和下游地区干旱传播关系

图6 干旱传播特征

4.3 干旱传播的滞后效应进一步分析干旱传播事件特征发现，滦河流域上、下游地区干旱传播事件平均历时分别为25.1、27.2 个月。其中，气象干旱约占传播事件的60%。水文干旱约占传播事件的80%，气象和水文干旱在干旱传播事件中的时间占比不同，体现出气象干旱向水文干旱传播过程中存在时间上的交错。本文将干旱的开始、峰值、结束作为关键时间节点，来探究气象干旱向水文干旱传播过程中的滞后效应（图7）。结果表明，干旱开始和结束时间的滞后特征比较稳定，上游较下游地区开始时间的滞后短，分别为1.8 个月和2.9 个月；而上游结束时间的滞后较上游长，分别为2.3 个月和1.9 个月。干旱峰值的滞后时间波动范围较大，在上游地区峰值强度在气象干旱和水文干旱的峰值强度出现先后顺序不定，平均气象干旱和水文干旱同时到达峰值，在下游地区，水文干旱峰值滞后于气象干旱峰值出现时间平均滞后约2.2 个月。由于本文所选择的上、下游面积相当，流域面积对传播过程的影响几乎可以忽略。考虑到滦河上游地区，而承德市主城区在下游地区，同时受到潘家口水库等大型水库的影响下，人工取用水和水库的调蓄等对流域产汇流过程的影响大［33］，这些人类活动的作用可能延长了气象干旱转换为水文干旱的时间，同时也缩短了气象干旱结束后水文干旱的持续时间。

图7 气象干旱向水文干旱传播的时间特征

4.4 干旱传播风险分析由于识别出的气象干旱向水文干旱的传播规律复杂，简单的线性关系难以评估干旱的传播关系。本节采用Copula 函数建立气象干旱和水文干旱变量间的联合分布模型，分别在滦河流域上、中、下游地区评估气象干旱向水文干旱的传播风险。首先对干旱历时和峰值强度进行边缘分布函数的拟合，再建立气象-水文干旱联合模型。如表1所示，为通过K-S 检验从正态（Nor⁃mal）、指数（Exponential）、伽马（Gamma）、对数正态（Lognormal）、对数逻辑斯特（Loglogistic）、广义极值（GEV）和韦布尔（Weibull）等7 种概率分布函数中优选的干旱历时和峰值强度的分布函数。对边缘分布的联合分布进行优选，使用AIC 检验优选出Gaussian、t、Clayton、Frank 和Gumbel 分布中最适合的气象-水文干旱联合分布模型，结果显示Gumbel 分布在构建各干旱特征间的相关关系时均为最适分布。进一步检验联合分布和经验分布的r2和RMSE均大于0.9，模型模拟精度高，适用于干旱传播风险评价。

表1 变量分布函数类型

基于不同情景下的联合分布模型，分别推算出上、下游地区气象干旱特征和水文干旱特征的传播联合重现期。如图8所示，随着气象和水文干旱历时和峰值的不断增大，气象水文干旱联合重现期也随之变化。在已经发生的干旱事件中，整体上干旱传播事件的重现期大于5年小于50年。上游地区干旱历时联合重现期比下游高，说明下游气象干旱与水文干旱历时的联系更加紧密。上游地区干旱峰值强度联合重现期比下游高，上游地区发生干旱传播的风险低于下游地区。

图8 干旱联合重现期和最优组合

图8所示为5 ～ 50年一遇联合重现期下，气象干旱与水文干旱的最可能组合情景设计值的线性回归模型参数。即建立了最可能的气象干旱历时和峰值强度向水文干旱传播的对应关系。表2 为模型相关参数和边界值，其中弹性系数（E）能够反应水文干旱特征随气象干旱特征变化的敏感程度。对干旱历时的传播而言，弹性系数能够反应气象干旱发生后水文干旱的恢复能力，弹性系数越大，水文干旱恢复速度越慢。滦河流域上、下游地区的气象和水文干旱历时间的弹性系数分别为2.26 和2.59，表明流域上游地区水文干旱的恢复速度高于下游地区。当考虑干旱峰值强度时，弹性系数表明水文系统对气象干旱峰值强度的响应情况，弹性系数小于1 能够体现出气象干旱向水文干旱的传播过程中对气象干旱峰值的削减，反之则反。滦河流域上、下游地区的气象和水文干旱峰值强度间的弹性系数均大于1，分别为1.17 和1.03，表明滦河流域气象干旱峰值强度在水文系统的传递过程中将进一步扩大，上游地区扩大速度高于下游地区。方程余项R表明触发干旱传播的风险阈值，当气象干旱发展到一定程度使得水文干旱特征值大于0 时，水文干旱才可能发生。用最可能组合情景回归模型来推求气象干旱传播到水文干旱的临界条件，结果表明，滦河流域上、下游地区气象干旱发展成水文干旱的临界历时分别为0.93 个月和1.26 个月，临界峰值强度分别为0.66 和0.44。即发生历时较短、峰值强度较低的气象干旱时，流域水循环中土壤水或地下水能仍够补给径流；或人为的水库放水、跨流域调水等系列抗旱措施能够补充河道水量［32-33］，使气象干旱不传播为水文干旱。而当气象干旱的历时和峰值强度超过临界条件后，发生水文干旱的风险高，以此为水文干旱的预警标准，能够进一步根据推求的干旱滞后时间预测下一场干旱发生的时间和到达干旱峰值的时间节点，提高水文干旱的预测能力。

表2 最可能组合情景回归模型相关参数和临界值

5 结论与展望

本文应用游程原理识别出滦河流域山区的气象和水文干旱事件，对气象和水文干旱事件的传播关系进行匹配，全面评价干旱传播事件中气象和水文干旱历时、峰值强度、滞后时间和传播风险等特征，得出以下结论：（1）滦河流域气象干旱的发生频率和峰值强度高于水文干旱。而水文干旱历时较长于气象干旱。（2）可将气象干旱和水文干旱的关系划分为：无明显的响应、一一对应、多对一和多对多4 种情景，前两种为主要关系，无明显关系情景在上、下游各占60%和59%，一一对应的情景占20%和22%；后两种关系所占比例较小，多对一情景在上、下游分别占12%和11%，多对多情景在上、下游均占8%。（3）从上游和下游地区气象干旱和水文干旱开始时间分别存在1.8 个月和2.9 个月不等的滞后；干旱峰值的滞后时间波动范围大；干旱结束时间滞后2.3 个月和1.9 个月。（4）各类Copula 函数中，Gumbel Copula 函数最适用于评估干旱传播历时和峰值强度的传播风险，滦河流域下游气象干旱与水文干旱历时的联系更加紧密，上游地区发生干旱传播的风险低于下游地区。由干旱特征联合重现期和最优组合模型得到的滦河流域上、下游地区气象干旱发展成水文干旱的临界历时分别为0.93 个月和1.26 个月，临界峰值强度分别为0.66 和0.44。

通过以上研究以能够更加全面了解流域干旱发生机制，为流域管理者提供气象灾害发生后的风险决策管理依据。本研究在传播时间阈值的选取等方面仍存在不确定性，且干旱传播过程涉及复杂的水文循环过程尚未不明确，未能量化不同影响因素（如气温、蒸散发、地下水以及流域前期土壤水、水体状况以及流域产流汇流等气象水文因素；水库、调水、灌溉等人类活动因素）对干旱传播特征的影响，我们将在未来的研究中通过模型模拟等方式进一步定量阐述。