气候变化对岩溶区流域极端径流频率分析的影响研究

2021-05-31王大洋王大刚

中国农村水利水电 2021年5期

王大洋，王大刚

（中山大学地理科学与规划学院，广州510275）

当今环境下，气候变化正时刻扰动着地表流域的水文循环过程，使许多流域的水文系统发生了明显的改变。极端气候的频繁发生导致极端水文事件频发，量级也明显增大，这其中最为突出的当属降水和洪水［1，2］。传统的洪水频率分析多是基于径流序列满足一致性的前提，而气候变化和人类活动的影响，很可能会使原始洪水序列的一致性遭到破坏。因此，研究气候变化条件下的洪水频率分析是非常重要且有现实意义的［3-5］。此外，岩溶区作为中国西南地区常见的地貌类型，其水文循环过程本身就较普通流域复杂，因此，有必要针对岩溶区流域开展气候变化对其极端径流频率分析的影响研究。

1 研究方法

1.1 趋势分析和突变检验

研究分别采用线性趋势分析和Mann-Kendall 检验方法对研究流域的径流序列进行趋势分析和突变点确定，从而将径流序列划分为两个时段。突变前时段被认为是未受到气候变化影响的时段，称为天然时段；突变之后时段被认为是受到气候变化影响的时段，称为气候变化影响时段。由于线性趋势分析和Mann-Kendall 检验方法较为常用，故本文不再赘述相关理论，具体可参考相关文献［6-7］。

1.2 径流频率分析模型优选

针对极端径流情况，研究采用广义极值分布（GEV）和广义帕累托（GPD）分布。选择两者的原因有两点：①与皮尔逊III（P-III）型曲线在中国的广泛应用相比，GEV分布和GPD分布在中国的研究起步相对较晚，但近些年的关注度在持续上升；②多数分布都是以两者的统一形式为基础进行衍生，如Gumbel分布、Fréchet 分布和Weibull 分布是GEV 分布的形状参数选取不同符号（正、零、负）的具体表现。Pareto 分布、指数分布和Beta分布分别是GPD 形状参数为（正、零、负）的具体表现［5］。GEV分布和GPD分布模型的表达式如下。

GEV的概率分布表达式为：

GPD的概率分布表达式为：

式中：μ、σ、ξ分别为模型的位置参数，尺度参数和形状参数，满足μ∈R，σ＞0，ξ∈R和1+ξ(x-μ)/σ＞0的条件。

1.3 样本抽取

样本抽选是频率分析的基础，针对GEV和GPD两种不同的分布模型，其抽样方法也是有所差异的。对于GEV，通常选择序列中的“年最大值（AM）”作为抽样依据；而对于GPD，通常选择超过某个定量值（POT）的所有数据来建模，适用于超门限值等系列的频率分析。本研究中，参考我国学者进行POT 模型分析的经验，以样本年限为基数，以年均超过定量值发生次数2倍为原则进行样本抽取，此外，对于POT 抽样方法，还应当注意抽取样本的独立性，从而避免因样本依赖导致的抽样误差［5］。

1.4 参数估计

由于极大似然估计法（MLE）结果能够在总体上较好地反映样本的统计信息，具有良好的统计性质［8，9］。因此本研究采用极大似然估计法进行分布模型的参数估计。用MLE 对GEV 和GPD进行参数估计的似然函数如下。

对于GEV，计算极大似然估计参数需要求解方程组：

对于GPD，其极大似然方程为：

式中：xi表示样本值；ε表示门限阈值。

为了便于求解，通常将其转化为对数形式的方程：

根据牛顿迭代法，可以计算得到相关参数。

1.5 模型优选检验

模型建立后，通过对观测样本进行拟合和参数估计，可得到各模型的拟合结果。结合模型检验方法可以优选出最合适的线型。本研究采用Kolmogorov-Smirnov 检验［10］（简称K-S 检验），PP概率图检验，QQ 图检验3种检验方法，根据GEV 和GPD模型的拟合结果进行检验，从而确定最优分布模型。

（1）Kolmogorov-Smirnov 检验。该方法以经验分布函数和理论分布函数之差为基础，进行模型拟合效果的检验，其检验统计量为：

式中：yi是经验频率去线上的点；F(yi)是累积分布函数，N为样本大小。

（2）PP 图和QQ 图。PP 累积分布概率（CDF）图是根据变量的累积概率对应所指定的理论分布概率来绘制散点图，它可以较为直观地反映样本数据与拟合模型分布的一致性。当样本数据符合某一分布时，样本点会较为接近45°对角线。QQ 图则和PP图比较相似，只是采用变量数据分布的分位数与所指定分布的分位数的关系进行比较。两种方法的基本原理如下。

若样本为{x1，x2，…，xN}，若经验分布函数为(x) ={x1，x2，…，xN}，模型估计的拟合分布函数为(x)，则可以通过(x)和(x)的差异反映模型拟合效果。若用概率进行表示（PP图），则为：

若用样本点的分位数对比，即QQ图：

2 研究区域概况

澄碧河发源于广西壮族自治区西部的百色市，属于珠江流域的西江水系，河流长度为127 km，流域集雨面积2 087 km2。高程由西北向东南逐渐走低，流域内分布着较多的典型石灰岩岩溶峰林，伏流河，落水洞较为常见，是世界上最大的连续岩溶地区之一。流域多年平均降水量为1 560 mm，因地形原因和岩溶发育较广，汛期降水多集中且易引发洪涝灾害，造成较为严重的社会经济损失和负担。

流域内部设有4个水文站，其中部分站点由于迁移、改测和撤销等原因使水文序列时间较短，且不能满足一致性要求。研究选取平塘站的1963-2017年的实测径流数据，该站点以上集雨面积占据整个流域的70%左右，且观测资料较为完整，因此可以作为研究数据。

3 研究结果及分析

3.1 趋势分析及突变检测

气候变化影响流域水文循环，水文循环中水分经过降水产流汇流过程，最终在河道的径流环节得以体现。研究选取澄碧河水库流域的年最大1日径流量、年最大3日径流量、年最大5日径流量和年最大洪峰流量作为极端径流的研究对象。首先采用线性趋势分析，对4个极端径流指标进行趋势分析，然后采用Mann-Kendall 方法对其进行突变检验，从而确定极端径流指标的突变时间点。

采用线性趋势分析时，研究发现年最大1日径流量、年最大3日径流量、年最大5日径流量和年最大洪峰流量均表现为下降的趋势，将其趋势线方程汇总在表1。从结果可知，最大5日径流的下降率最为明显，年最大3日次之，年最大洪峰流量相比之下降率最小。但4 个极端径流指标总体均表现为下降趋势，较为一致。

表1 各个水文特征值趋势线方程

采用Mann-Kendall 检验方法对澄碧河水库流域的上述4个极端径流指标进行趋势分析和突变检验，结果见图1。从结果来看，年最大1日、年最大3日、年最大5日径流的趋势上具有较强的一致性。UF曲线的正负显示了特征值的变化趋势，图1表明，年最大1日、年最大3日、年最大5日径流在1990年之后，UF一直小于0，说明三者均表现出明显的下降趋势，且这一下降趋势在1998年左右超出了α= 0.05 的显著水平，表明下降趋势显著。相比之下，年最大洪峰流量，虽然在1983年之后呈现下降趋势（UF＜0），但未超出置信区间，因此其表现为不显著下降。此外，由图可知UF和UB的曲线出现交点，且该交点位于置信区间范围内部。从交点位置来看，年最大1日径流量、年最大3日径流量、年最大5日径流量的交点位置均出现在1991年，而年最大洪峰流量的交点较多，为1973、1976、2002 和2013年。综合各检测结果，确定澄碧河流域的气候突变发生在1991年。

3.2 序列分段与样本抽取

根据检测出的突变点，将整个样本分为两个时段：天然时段和气候变化影响时段。研究中天然时段为1963-1991年，气候变化影响时段为1992-2017年。因此，分别以两段时间为基础进行样本抽选。

样本抽选是频率分析的基础，针对GEV和GPD两种不同的分布模型，研究采用不同的适应性抽样方法。对于GEV，结合“最大值”原则，本研究中选取最大1日径流量、年最大3日径流量、年最大5日径流量和年最大洪峰流量为GEV 分布拟合的基础数据；对于GPD，本研究，以年均超过定量发生次数2 倍为原则，针对突变前和突变后分别以58 和52 为样本容量进行样本抽取，针对澄碧河水库的汇流时间，以3 d 作为抽样间隔，以保证每个样本的独立性。

3.3 分布模型拟合及优选

研究选取了GEV和GPD分布模型进行拟合与对比分析，其突变前、突变后时段拟合的PP 图和QQ 图为图2和图3，同时将拟合模型的各个参数值和K-S 模型检验结果汇总在表2。图表中am 和pot 分别表示年最大取样法和超定量取样法，pre 和pos分别表示气候突变前、后时段，1 d、3 d、5 d 和peak 分别对应最大1日径流量（万m³）、年最大3日径流量（万m³）、年最大5日径流量（万m³）和年最大洪峰流量（m³/s）。

从GEV分布和GPD分布模型的PP图和QQ图可知，两个分布均能较好地反映经验观测值与模拟值之间的关系，从散点图与45°对角线的关系可以看出，两者也都紧密地环绕在对角线周围。通过K-S 检验量化的统计结果（表2）发现：①所有情况下的统计量值均小于对应的显著性标准值。由此表明，GEV 分布和GPD 分布均能通过α=0.05 的显著性检验标准。②除了天然时段的年最大5日径流量（pre-5d）和气候变化影响时段的年最大1日径流量（pos-1d）之外，其他所有情况下的K-S 检验结果中GPD 分布对应的统计量均小于（或等于）GEV 分布对应的统计量。因此综合上述比较分析可知，虽然GEV分布和GPD分布均能很好地反映澄碧河流域的极端径流序列，但相比之下，GPD分布在拟合澄碧河流域的极端径流方面更具优势。

表2 GEV和GPD分布模型参数及K-S检验结果

从GPD 分布的参数变化来看，突变后时段相比突变前时段，年最大1日径流量和年最大洪峰的形状参数均有所减小，而年最大3日径流和年最大5日径流则表现为形状参数增大。尺度参数方面，除最大1日径流增大外，其余3 个指标均有所减小。而位置参数方面，只有年最大洪峰流量表现为增加，其余则都是减小的趋势。由此表明，分布参数的变化在一定程度上受到气候变化的影响。

3.3 气候变化对各重现期下极端径流的影响

通过分布模型的对比分析，确定澄碧河流域的极端径流频率分析模型为GPD 模型。因此，以POT 抽取的样本观测数据为基础，采用GPD 分布进行数据拟合，采用极大似然估计进行参数估计得到的线型可以作为重现期计算的依据。根据重现期与累积概率分布函数的关系可以计算天然时段和气候变化影响时段的各重现期及其对应的极端径流，进而研究在气候变化条件下各重现期对应的极端径流变化情况。

对计算的两个时段的重现期进行分析，结果见图4。

从图4和表3变化可知，突变后时段的各个极端径流指标在不同重现期水平下，较突变前均有不同程度的减小。其中，年最大洪峰变化最为明显，10年一遇重现期水平下，其值由550.54 m3/s 减小至448.09 m3/s，减小比例为18.61%；100年一遇重现期水平下，其值由1 141.02 m3/s 减少至602.63 m3/s，减小比例为47.18%；1 000年一遇重现期水平下，其值由2 130.04 m3/s减少至718.52 m3/s，减小比例为66.27%。其次是最大1日径流量，其减小比例为23.03%～65.51%。相比之下，年最大5日径流量变化最小，其减小比例为12.95%～21.94%；由此可见，气候变化在一定程度上确实影响了极端径流的频率，使其发生了变化。

表3 澄碧河流域气候变化前后极端径流重现期

4 讨论与结论

4.1 讨论

受气候变化的影响，处于岩溶区的澄碧河流域的极端径流序列平稳性发生了变化，作为极端径流指标的洪峰和洪量序列的平稳性遭到了破坏，都由原本的“平稳序列”转变为了“非平稳序列”。综合分析，序列的突变年发生在20 世纪90年代初。研究表明亚洲区域，尤其是中国范围内，气候变化在90年代表现最为剧烈，这在一定程度上佐证了本研究的正确性。此外，据广西百色市《凌云县县志》考证，澄碧河流域自20 世纪70年代之后就没有再有过大规模的人类活动，因此造成极端径流频率变化的原因应该还是来自于气候变化。气候变化影响着流域的气温、降水、蒸发等水文循环过程，最终表现在极端径流上［11］。

从结果来看，气候变化影响最大的是年最大洪峰流量和年最大1日径流，这可能与澄碧河流域的山区地形和岩溶区地貌有关。澄碧河流域弄塘以上流域多为岩溶区，岩石裂隙、洼地、溶沟、溶槽和地下暗河等地质形态发育而成了交错纵横的地下河网体系。加上该地区土壤相对疏松，水流垂向运动频繁，暴雨转化为径流的效率较高，一旦发生极端降水形成大洪水，也能很快被削峰坦化。从流域的汇流面积及汇流时间分析，即便是大洪水，基本也能在3 d 左右完成汇流过程。所以该流域对于短时降水引发的极端径流，表现是较为明显的。这与研究得到的年最大洪峰流量和年最大1日径流变化率较大，而年最大5日径流量变化率最小的结论相符。通过岩溶区释水过程及径流组分的相关研究［12，13］可知，由于岩溶洼地、溶沟、溶槽和地下暗河等地质构造的存在，使得岩溶区流域释水过程要明显快于非岩溶区流域，其衰减系数也相对较大。此外，岩溶区径流组分随洪峰的变化也较非岩溶区剧烈，这也是在气候变化影响下，极端径流频率变化明显的重要原因。