淮河流域极端气候事件非平稳特征研究

王怀军，曹 蕾，肖明贤，冯 如

（1.淮阴师范学院城市与环境学院，江苏淮安223300；2.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，南京210029；3.南京水利科学研究院水利部应对气候变化研究中心，南京210029）

气候变化导致极端气候事件的发生频率、强度、空间范围及持续时间发生改变，并可能引发前所未有的极端事件［1-3］。中国是世界上受自然灾害影响最为严重的国家之一，也是极端气候事件发生频率与强度较高的国家之一［4］。随着全球气候变化以及中国社会经济快速发展对资源环境和生态产生的压力加剧，使得极端气候事件防范应对形势更加严峻复杂［5］。据民政部和国家减灾委员会等部门统计，1990-2010年气象灾害造成的直接经济损失高达2 000 亿元/a，约相当于当年GDP 的1%～3%［6］。并且近年来气象灾害损失有日益加重趋势，气象灾害的频繁发生已成为制约我国国民经济持续稳定发展的主要因素之一。

目前，在极端气候事件时空变化研究方面，主要从以下方面进行研究：①极端气候事件的表征。基于单一台站监测技术的定义，采用“气候极值”量化异常天气气候现象，当气候要素达到定义的阈值时，便可认为极端事件发生［1，7，8］。显然，阈值的确定受人为因素的影响较大；②极端气候事件时空变化特征。主要利用趋势检验、突变分析、小波变换以及地统计学技术分析某流域或者地区的极端气候事件时空分布规律［9，10］；例如，尹义星等［11］对1951-2013年江苏省极端最高和最低气温变化趋势进行研究，结果表明冰冻日数和冷昼日数呈下降趋势，夏季日数和暖昼日数呈上升趋势；③基于平稳极值函数的时空分布特征。极端气候事件是一个随机变量，利用广义极值分布（GEV）和广义帕累托分布（GPD）对其进行拟合以及重现期水平提取［12，13］。如赵玲玲等［14］利用GEV 对华南强降雨特性进行分析，余锦华等［15］采用GPD 对中国东部夏季极端降水频率特征进行分析，结果均表明GEV和GPD模型可以很好的拟合极端气候事件。然而，气候变化和人类活动导致极端气候事件呈非平稳性，传统极值函数分析中，极端气候事件的模拟基于平稳的边际分布函数，这往往忽视气候变化背景下水文变量非平稳特性的影响［16］。温庆志等［17］指出淮河流域年最高气温表现出非平稳性，增幅达1.5 ℃。Kalai 等［18］研究结果表明由于土地使用方式、人类干预和气候变化的迅速变化，非平稳区域洪水频率分析（RFFA）能够更有效地捕获洪水的时变行为。Gu等［19］指出城市化导致局部地区（例如中国北方）的极端降水呈现非平稳，城市化地区气候发生严重的非平稳性的可能性更高。Yang 等［20］则指出在平稳和非平稳条件下检测到的城市洪水存在差异，需要将更多的精力进一步研究非平稳条件下的城市洪水变化。因此，有必要将数据的非平稳性引入到GEV 和GPD 函数中，以便更好地揭示极端气候事件的时空分布特征。

淮河流域地处我国南北气候过渡带，是中国重要的地理生态分界线和生态环境脆弱区，对全球气候变化十分敏感［21］。由于淮河流域本身特殊的地理环境条件，气候系统、水循环系统的不稳定性增大，导致极端气候事件发生的频率和强度都有所增大，这对人民生产生活产生了不可忽视的影响［22］。本研究选择位于南北气候过渡带这一特殊地理环境下的淮河流域为研究区，借助非平稳GPD 和GEV 模型拟合极端气候事件，揭示非平稳条件下极端气候事件时空演变规律，以期为淮河流域减灾防灾提供决策参考。

1 研究区概况、数据与方法

1.1 研究区概况

淮河流域地处我国东部，介于中国两大流域-长江与黄河之间，其地理位置位于30°55′～36°36′N，111°55′～121°25′E 之间（见图1），流域面积27 万km2。淮河流域包括鄂、豫、皖、鲁、苏五省40 个地（市）181 个县（市），流域内总人口1.65 亿人，平均人口密度居各大江大河流域人口密度之首，为611 人/km2。流域地处我国南北气候过渡地带，年平均气温在11～16 ℃之间。最高月平均气温出现在7月份，约为25 ℃，最低月平均气温出现在1月份，约为0 ℃。极端气温最高达44.5 ℃，最低达-24.1 ℃，是中国极端气温事件的高发区之一。多年平均降水量约为920 mm，其分布情况大致是由南向北递减，平原少，山区多，内陆少于沿海。

图1 研究区和气象站点分布图Fig.1 Study area and distributed weather stations

1.2 数 据

1960-2018年淮河流域内70 个站点的每日降水和气温数据来自于中国气象科学数据共享服务网（http：//www.cdc.cma.gov.cn），数据经过重重检查，质量良好。本文首先借助R 语言提取年最大值序列（AM）和超门限序列（POT）。AM 序列包括1日最低气温（日最低气温的年最低值，Tnn）、1日最高气温（日最高气温的年最大值，Txx）和1日最大降水量（日降水量的年最大值，Rx1day）。POT 序列阈值取第95%百分位的日最高气温、日最低气温和日降水量，分别命名为超门限最低气温（Tmingpd）、超门限最高气温（Tmaxgpd）和超门限降水（Pregpd）。

1.3 方 法

1.3.1 平稳性检验

时间序列的平稳性可以表示如下：①均值E（X）=μ是与时间t无关的常数；②方差Var（Xt）=σ2是与时间t无关的常数；③协方差Cov（Xt，Xt+k）=γ k是只与时期间隔k有关，与时间t无关的常数。常用的判断方法包括图示法、单位根检验法、自相关函数（ACF）判断法。图示法通过时间序列的路径图来粗略地判断它是否是平稳的，平稳序列表现出一种围绕其均值不断波动的过程，非平稳序列则往往表现出一定的上升或者下降趋势。为方便进行处理，我们采用趋势和突变检验来进行检验，如果时间序列存在显著趋势和突变点，则说明该序列非平稳。趋势检验采用Man-Kendall 检验，均值突变检验采用Pettitt 检验，方差突变采用单变异点检测方法（AMOC）［23，24］。Man-Kendall 检验通过计算统计量Zs来进行趋势统计的显著性检验［25，26］。当Zs＞0 时，表示时间序列呈增加趋势，Zs＜0 时，表示时间序列呈减少趋势，如果Zs绝对值大于1.96，表明时间序列趋势显著（0.05 显著水平）。Man-Kendall 检验中还计算了倾斜度Sen，表示趋势的方向和大小。计算公式分别为：

式中：1＜i＜j＜n，n为数据资料时间长度，如果Sen的值为正，表示呈现增加趋势，如果Sen的值为负则表示呈现减小趋势。

Pettitt 非参数变点检验假设一个水文时间序列y1，y2，…，yn中的突变点出现在t时刻，将t作为分界点，将该水文时间序列分为两部分：y1，y2，…，yt和yt+1，yt+2，…，yn分别遵循分布Fy1（y）和Fy2（y），然后检验在显著性水平α下突变点的显著度，具体计算过程参照文献［27，28］。

另外文章采用ADF 检验进行单位根检验，对于时间序列AR（p）模型，有

式中：c为漂移项；β、φi为参数；εt为随机误差项，是服从独立分布的白噪声过程。

ADF 检验的原假设H0∶β=0，存在单位根，为非平稳序列，从而可以根据ADF单位根检验结果判断序列是否平稳。

1.3.2 平稳和非平稳GEV模型

GEV模型适合用来处理AM序列（如年最大降水量、年极端最高和最低气温），其累积分布函数为：

式中：μ，σ，ξ分别为位置、尺度和形状参数。

GEV模型T年重现水平为：

式中：P为不同重现期对应的概率。

由于在方差突变检验中，所有的极值时间序列均不显著，且形状参数一般不易改变，因此非平稳GEV 分布模型仅设置为位置参数随时间变化［μ（t）=ρ0+ρ1t］，形式如下：

1.3.3 平稳和非平稳GPD模型

POT序列采用GPD分布进行模拟，其累计分布函数是：

式中：u，σ，ξ分别为分布的阈值、尺度和形状参数。

假定某一极端事件的重现期为1/（1-p），则其对应的重现水平为：

由于阈值为95%分位数，形状参数一般不易改变，因此非平稳GPD 分布模型定义如下［σ（t）=sin（2 π×（t-1960）/365.25］+cos［2π×（t-1960）/365.25）］：

对平稳和非平稳模型进行似然比检验（LR 检验），统计量为：

式中：L1为复杂模型最大似然值；L2为简单模型最大似然值；LR近似符合卡方分布。

对于GEV 模型，分别构造了GEV0（无协变量）和GEV1（位置参数随时间线性变化）；对于GPD 构造GPD0和GPD1（尺度参数随时间线性变化）。通过构造的复杂模型与简单模型比较来检验加入非平稳性后是否能够显著地提高模拟效率。

最后利用选择的平稳或非平稳GPD和GEV模型，通过设置不同的重现期，得到每个站点2年、10年、20年、30年、50年、100年一遇的极端气候事件重现水平，利用克里金插值得到重现水平空间分布图。

2 结果分析

2.1 极端气候事件时空变化趋势及平稳性分析

关于最高气温，淮河流域Tmaxgpd和Txx变化趋势为0.07 d/a和0.01 ℃/a，且变化走势相对一致，都为1960-1985呈减少趋势，而在1985年以后呈增加趋势［图2（a），图2（d）］。从空间变化趋势来看，上述两个指标大部分站点表现为增加趋势，但表现为显著变化的站点相对较少，分别为14%和21%［表1，图3（a），图3（d）］。对于最低气温，Tmingpd和Tnn均为显著增加趋势，变化幅度分别为0.16 d/a和0.04℃/a［图2（b），图2（e）］，其变化幅度均高于最高气温，这也说明与最低气温相关的极值要高于与最高气温相关的极值，这与其他地区的研究结果一致［29］。空间变化趋势具有与区域时间一致的变化趋势，显著增加的站点分别达到了59%和76%。Pregpd和Rx1day的区域变化趋势分别为0.01 d/a 好0.11 mm/a，但不管是从区域上还是空间上均为不显著变化。

图2 淮河流域极端气候事件区域时间变化趋势Fig.2 Regional times series of climate extremes in the Huaihe River Basin

表1 极端气候事件空间变化趋势百分比 %Tab.1 Percentage showed different trends of extreme climate events

为进一步分析极端气候指数的平稳性，采用了ADF 检验和Pettitt 检验进行分析（图4）。79%站点的Tmaxgpd和87%站点的Tmingpd的ADF 检验P值大于0.05，说明具有单位根，表现为非平稳；Txx和Tnn具有单位根的百分比为46%和36%，说明在此检验下表现为非平稳的站点百分比不到一半；Pregpd和Rx1day仅有少部分站点具有单位根，说明在淮河流域降水极值表现为平稳性变化。在Pettitt 检验下［图4（b）］，Tmingpd和Tnn变现为显著变化的站点分别为71%和74%，说明大部分站点在此统计量下表现为非平稳。其他指数表现为显著变化的站点比例很少，均小于30%，且突变年份分布很广［图4（c）］，说明大部分站点表现为非平稳变化。从趋势分析和相关平稳性结果可以看出，Tmingpd和Tnn在淮河流域表现为非平稳性，Pregpd和Rx1day表现为平稳性，Tmaxgpd和Txx表现相对较为复杂，平稳性和非平稳性同时存在。

2.2 基于GEV模型的极端气候事件变化特征

对AM 序列进行GEV 分析。首先利用平稳性GEV 模型对所有站点的GEV 序列进行KS 拟合优度检验［图5（a）］，发现几乎所有的站点的P值均大于0.05，这说明在KS 统计量下所有AM序列均满足平稳性GEV模型。该研究结果表明在进行非平稳性分析时，不能通过拟合优度检验来探查数据的平稳性。通过LR 检验来分析非平稳模型是否对平稳模型有改进［图5（b）］，对所有站点的AM 序列进行平稳和非平稳GEV 函数拟合（图6）。结果表明Tnn80%的站点P值小于0.05，这说明大部分站点的非平稳模型对平稳模型有改进。而Txx和Rx1day，非平稳模型改进的比率仅有30%和1%，这说明大部分站点满足平稳性。所以在后续GEV 分析中，Tnn采用非平稳模型，而Txx和Rx1day采用平稳性模型。

图3 淮河流域极端气候事件空间变化趋势Fig.3 Spatial trend of extreme climate events in the Huaihe River Basin

图4 极端气候事件平稳性检验及突变年份Fig.4 Stationary test and the year of abrupt change of extreme climate events

对AM 序列的2年（RTL2）、10年（RTL10）、20年（RTL20）、30年（RTL30）、50年（RTL50）和100年（RTL100）重现期重现水平进行求取，发现不用重现期下重现水平具有较高的相关系数（图7），这说明不同重现期下极端气候事件重现水平具有一致的空间分布特征，因此，后续研究仅分析30年重现期重现水平的空间分布（图8）。淮河流域Rx1day30年重现水平西北向东南递增［图8（a）］，这可能与东南季风的逐渐增强有关。Txx从西往东减少［图8（b）］，由于海洋对流域东部地区气温的调节作用，离海洋越近受海洋的影响越大，导致Txx自东往西逐渐升高。由于将Tnn视为非平稳模型，重现水平会随着时间而变化。为方便表示，将重现水平分为两个时间段，1960-1989［图8（c）］和1990-2018［图8（d）］，分别对两个时间段重现水平求均值。结果显示，Tnn-pre和Tnn-post具有一致的空间分布趋势，重现水平从南向北递增，这与纬度地带性规律一致。但值得注意的是两个时期的Tnn相差了近乎1～2 ℃，这说明在淮河流域Tnn的增温趋势非常明显。

图5 淮河流域极端气候事件AM序列拟合优度检验Fig.5 Goodness-of-fit test of extreme climate events of AM series in the Huaihe River Basin

图6 基于平稳和非平稳模型的泗洪站Tnn（取绝对值）拟合优度图Fig.6 Goodness-of-fit test of absolute Tnn of Sihong station based on stationary and non-stationary models

2.3 基于GPD模型的极端气候事件变化特征

采用平稳和非平稳性GPD 模型对所有站点POT 序列进行拟合，并利用LR 检验验证模型是否有改进。结果表明，7%、78%和93%的非平稳GPD 模型对降水、最高气温和最低气温有改进（图9、图10），这说明Pregpd表现为平稳变化，Tmaxgpd和Tmingpd表现为非平稳变化。因此，采用平稳GPD 模型对Pregpd模型进行拟合，非平稳模型对Tmaxgpd和Tmingpd进行拟合。同Rx1day空间分布特征一致，Pregpd空间部分从西北向东南递降［图11（a）］，同样反映了东亚季风从东南向西北逐步减弱的态势。将Tmaxgpd和Tmingpd重现水平分为两个时间段，1960-1989（pre）和1990-2008（post），分别对两个时段重现水平求均值，分析其空间分布状态。Tmaxgpd-post和Tmaxgpdpre反映了同样的分布趋势，从东部向西部递增［图11（b）、图11（c）］，这反映了海洋的热力调节作用。前后两个时期的重现水平相对变化较少，且有正有负［图10（f）］，这说明Tmaxgpd30年重现水平并不一定为增加趋势，这可能与部分站点并不为非平稳性有关。Tmingpd-pre和Tmingpd-post空间分布趋势具有一致性，从南向北递降［图11（d）、图11（e）］。相较于Tmaxgpd，Tmingpd-post比Tmingpd-pre数值都大，且大部分地区超过0.5 ℃以上［图11（g）］，这说明最低气温重现水平呈显著增加趋势。

图7 30年重现水平与2、10、20、50、100年重现水平相关系数图Fig.7 Correlation coefficients between the 30-year return level and the 2，10，20，50，and 100-year return levels

图8 AM序列30年重现水平空间分布图（其中Tnn取绝对值）Fig.8 Spatial distribution of the 30-year return level of AM extremes

3 结 论

基于淮河流域最新气象站点日值气温和降水数据，计算目前国际上常用的极端气候事件指数，分析极端气候事件的时空分布趋势。基于平稳和非平稳GEV和GPD极值模型，分析极端气候事件在不同极值模型下的重现水平空间变化特征，主要结果如下。

（1）对于最低气温，超门限最低气温（Tmingpd）和1日最低气温（Tnn）呈显著增加趋势，且空间上亦表现为显著增加趋势。超门限最高气温（Tmaxgpd）、1日最高气温（Txx）在区域上表现为不显著增加，且发生显著变化的站点百分比相对较低。超门限降水（Pregpd）和1日最大降水量（Rx1day）不管是从区域上还是空间上均表现为不显著变化。

图9 淮河流域极端气候事件POT序列LR检验Fig.9 LR test of extreme climate events of POT series in the Huaihe River Basin

（2）不同的平稳性检验方法会得出一定差异的平稳性检验结果，单位根检验下Tmaxgpd和Tmingpd表现为非平稳，Pregpd和Rx1day表现为平稳性变化。在Pettitt检验下，Tmingpd和Tnn表现为非平稳。

（3）AM 序列中，Tnn大部分站点非平稳模型对平稳模型有改进，Txx和Rx1day非平稳模型改进的比率相对较低，说明在AM 序列中，Tnn宜采用非平稳模型，而Txx和Rx1day采用平稳模型。淮河流域Rx1day30年重现水平西北向东南递增，Txx从西往东减少，Tnn绝对值从南向北递增，且不同时期重现水平变化显著。

图10 基于平稳和非平稳GPD模型的泗洪Tmingpd拟合优度图Fig.10 Goodness-of-fit test of Tmingpd of Sihong station based on stationary and non-stationary models

图11 POT 序列30年重现水平空间分布Fig.11 Spatial distribution of the 30-year return level of POT extremes

（4）POT 序列中，Tmaxgpd和Tmingpd非平稳模型改进的比率超过78%，这两个指数宜采用非平稳模型，Pregpd改进比率较低，宜采用平稳GPD 模型。Pregpd重现水平空间分布从西北向东南递减，Tmaxgpd从东部向西部递增，Tmingpd从南向北递降。 □