岩溶区流域年最大洪峰流量变化特征分析

莫崇勋，班华珍，谢燕平，何嘉奇，阮俞理，孙桂凯

(1. 广西大学土木建筑工程学院，南宁 530004；2. 工程防灾与结构安全教育部重点实验室，南宁 530004；3. 广西防灾减灾与工程安全重点实验室，南宁 530004)

0 引 言

目前，关于洪峰流量的研究成果多聚焦于场次洪峰流量[1-5]，较少涉及洪峰流量序列演变规律方面。纵有研究，也多集中于非岩溶区，如任健等[6]运用HHT分析法研究黄河花园口站年最大洪峰流量序列的多时间尺度特征；唐权辉等[7]基于北江干流4个水文站资料建立年最大洪峰流量序列，利用Mann-Kendall和EMD法对序列进行特征分析；孙思瑞等[8]采用水文变异诊断系统对洞庭湖三口的年最大洪峰流量和年最高洪峰水位序列进行变异诊断，结果表明各站点洪峰流量多呈现下降的趋势变异。这些研究都是针对非岩溶区进行的，而岩溶地貌作为世界上最脆弱的地貌类型之一，对环境变化的响应更为敏感。随着气候变化和人类活动影响的加剧，对岩溶区流域的洪峰流量序列进行特征分析具有一定的研究意义。鉴于此，论文基于岩溶区坝首水文站1963-2014年实测逐日入库流量，运用水文变异诊断系统中的多种统计方法和Morlet小波分析法，研究年最大洪峰流量序列的演变规律，以期为岩溶区流域澄碧河水库制定防洪减灾对策和水环境治理措施提供科学依据。

1 研究区概况与数据来源

澄碧河隶属于珠江流域西江水系右江干流的一级支流，发源于广西百色市凌云县青龙山北麓一带，途经弄林、浩坤、平塘、百色右江区等地，水流呈明伏流交替出流，至那东村为澄碧河水库坝址。水库坝址以上流域集雨面积2 000 km2，流域地貌以弄林为界分为两个部分，浩坤～弄林以上为典型的岩溶峰林区，浩坤以下为非岩溶区。岩溶区流域，由于石灰岩裸露，在长期风化作用下，形成大量的裂隙、溶沟、漏斗和落水洞，雨洪过程中，地面径流通过这些通道大量涌入地下，形成地下径流，短历时内对洪峰起到削峰作用。流域内的浩坤溶洞伏流河形成天然水库，库容约3 亿m3，对澄碧河水库起调节作用。流域形状近似矩形，地势西北高而东南低，流域概况图如图1。

图1 澄碧河流域水系图Fig.1 Chengbi river Basin

澄碧河流域内共有4个水文站，分别为坝首、平塘、浩坤和下甲水文站。从地理位置看，浩坤站和下甲站较能体现岩溶区的特性，但因这两个水文站建站时间较晚，水文序列均不足30年，平塘站水文序列不够完整。坝首水文站位于澄碧河流域出口断面处，站点数据系列足够长，可靠性高，且一定程度上能表征岩溶区特性。故论文基于广西岩溶区流域澄碧河水库坝首水文站1963-2014年共52 a的日均入库流量数据，采用年最大值法对洪峰流量进行选样，得到52 a的年最大洪峰流量序列。

2 研究方法

2.1 趋势和突变检验方法

目前，可用于检验水文时间序列变异情况的方法很多，但因水文时间序列(特别是洪水序列)变异形式相对复杂，单一检验方法往往可靠性较低。故本文采用谢平等[9]提出的水文变异诊断系统中多种统计方法，对年最大洪峰流量序列的趋势成分和突变(跳跃)成分进行检验分析。检验过程分为3个阶段：①初步诊断，选取线性趋势法、滑动平均法和Hurst系数法对序列的趋势和突变进行检测，定性分析序列的变异情况；②详细诊断，选用线性回归法、Spearman秩次相关法和Man-Kendall秩次相关法共3种方法分析序列的趋势显著性，结合Mann-Kendall检验法、Lee-Heghinan检验法、滑动 T 检验法、滑动秩和检验法、有序聚类检验法、累积距平检验法、滑动游程检验法和滑动F检验法共8种变异诊断方法对序列的突变进行分析；③综合诊断，对详细诊断结果进行汇总，趋势检验结果一般较统一，对显著性进行求和即可，若结果大于等于1则显著，反之，则不显著；突变检验结果相对复杂，除了进行显著性综合外，还需进行突变权重综合，论文以所有可能突变点权重最大者作为最可能突变点。上述方法的具体过程见相关文献[10-13, 15]。

2.2 Morlet小波分析法

小波分析是广泛应用于研究水文序列周期性规律的方法，Morlet小波分析的基本原理及计算步骤[14]如下。

小波函数ψ(t)∈L2(R)且满足：

(1)

式中：ψ(t)为小波基函数，其作用是通过尺度伸缩和在时间轴上的平移构成一簇函数系。

(2)

式中：ψa,b(t)为基小波；a为尺度因子；b为时间平移因子。

由于年最大洪峰流量序列是离散的，设函数f(kΔt)(k=1，2，…，N)，Δt为时间间隔。则离散序列小波变化表示为：

(3)

式中：Wf(a,b)为小波变化系数；其余同上。

对时间域上关于a的全部小波系数的平方进行积分，即为小波方差：

(4)

3 结果与分析

3.1 趋势性分析

采用一元线性法和5 a滑动平均法对年最大洪峰流量序列的趋势性进行初步检验，结果见图2。由图2可知，年最大洪峰流量序列整体呈微弱下降趋势，年均下降率为6.04 m3/s；1979年之前，洪峰流量序列大多处于均值上方，其中1967年达到最大值，为3 399 m3/s，1977年次之，为3 075 m3/s；1979年之后，洪峰流量序列在均值附近平稳波动，其中2013年达到最小值，为301.77 m3/s。为进一步检验序列趋势变化的显著性，采用线性回归法、Spearman秩次相关法及Man-Kendall秩次相关法进行详细诊断，显著性α均取为0.05，结果见表1(图略)，经计算，线性回归法和Spearman秩次相关法的统计量T值分别为-1.1和-0.15，绝对值均小于α=0.05显著性水平临界值1.64；Man-Kendall秩次相关法的统计量U为-0.209，绝对值小于α=0.05显著性水平临界值1.96；记不显著性为“-1”，显著性为“1”，则综合诊断趋势显著性为“-3”，表明年最大洪峰流量序列的变化趋势不显著，即序列呈不显著的下降趋势。

图2 年最大洪峰流量序列趋势线Fig.2 The trend of annual maximum peak discharge series

诊断方法显著值(|T|/|U|)检验结果趋势诊断线性回归法T=-1.1不显著(-1)Spearman秩次相关检验法T=-0.15不显著(-1)Kendall秩次相关检验法U=-0.209不显著(-1)综合显著性--不显著(-3)

3.2 突变性分析

采用R/S分析原理计算出年最大洪峰流量序列的Hurst系数[15]h值为0.69，表明序列呈现中变异，初步判定该序列可能存在趋势或突变(跳跃)变异。由3.1节知，趋势变化不显著，故主要对突变(跳跃)成分变异情况进行详细诊断。采用Man-Kendall检验法、Lee-Heghinan法等8种变异检测方法对序列进行分析，显著性α均取为0.05，Man-Kendall检验法结果如图3所示。由图3可知，统计量UF和UB在α=0.05显著性水平临界线内存在多个交点，分别为1966年、1971年、1980年、1994年、1998年和2011年，因1966、1971、1998和2011年过于接近序列两端，可靠性相对较低，故选取1980年和1994年为可能突变点；其余方法的检测结果见表2。由表2可知，各方法检测结果不一致，突变点及其显著性均有差异，故进行综合诊断分析，结果见表3。由表3可知，突变(跳跃)点为1979年，其综合权重达到最大值为0.375，综合显著性为3，表明突变显著。

图3 Mann-Kendall变异检测Fig.3 Mann-Kendall Mutation detection

综上可知，52 a来，年最大洪峰流量序列整体呈不显著的下降趋势，并于1979年发生以跳跃(突变)形式为主的变异。

3.3 周期性分析

利用Morlet小波分析序列的周期，结合Matlab软件绘制出如图4和图5所示。图4(a)为小波系数实部等值线图，小波系数正相位表示洪峰流量处于偏高值，负相位表示处于偏低值；由图4可知，年最大洪峰流量主要存在5～10、11～20和25～30 a的时间尺度周期性变化规律。对应于时间尺度的交替，序列出现明显变化，其中11～20 a尺度震荡最强烈，并以15 a为中心时间尺度，正负交替持续到1985年前后；5～10 a表现较明显，其中心时间尺度为8 a左右，正负交替持续到1975年；25～30 a尺度出现得较微弱，其中心时间尺度为28 a左右，正负交替持续到1985年左右。1985年之后，洪峰流量未出现明显的时间尺度周期性规律，但其相位几乎全为正值，表明洪峰流量于1985年之后几乎全部大于均值。图4(b)为小波系数模平方图，反映不同变化周期对应的能量密度在时间域的分布强弱，对应值越大，能量越强。由图可知，序列存在两个相对显著的能量聚集中心：5～10和11～20 a的时间尺度，表明洪峰流量变化对应于这两个尺度周期性比较显著。其中以11～20 a时间尺度震荡能量变化最强，主要发生在1960-1985年，震荡中心在1967年；5～10 a时间尺度十分突出，主要发生在1960-1975年，震荡中心也为1967年。图5表示波动能量尺度分布的小波方差图，由图可知，年最大洪峰流量序列存在3个比较明显的峰值。其中，最大峰值对应15a的时间尺度，表明15 a的周期震荡最强烈，即为年最大洪峰流量变化过程的主周期；8和28 a对应着第二峰值，即为序列的次周期。

表2 年最大洪峰流量序列详细诊断结果Tab.2 The detail diagnosis of annual maximum peak discharge series

表3 年最大洪峰流量序列综合诊断结果Tab.3 The comprehensive diagnosis of annual maximum peak discharge series

图4 年最大洪峰流量序列Morlet小波变化系数时频分布Fig.4 The time-frequency distribution of the Morlet wavelet coefficients of annual maximum peak discharge

图5 年最大洪峰流量序列小波方差图Fig.5 The wavelet variance of annual maximum peak discharge

4 讨论与小结

4.1 讨 论

影响径流量的因素众多，其中气候变化和人类活动是两个重要因素[16]，洪水是径流的一个极端表现，洪峰是洪水过程的一个重要指标[17]，降雨、径流、气温等水文气象要素变化一定程度上会引起洪峰流量的改变。澄碧河流域是典型的岩溶区与非岩溶区结合的流域，弄林以上岩溶区以峰从洼地为主，其特点为峰林峰从连座、洼地闭塞，土壤覆盖层稀薄，岩溶裂隙、漏斗、溶沟和落水洞很发育，形成具有不规则性、多层性和连通性的地下河网。大部分洼地的消水溶洞与地下暗河、管道连接，形成畅排型的地下河网通道，能快速补给地下径流；少部分洼地虽不存在消水洞，但土壤覆盖层松散薄弱，垂向补给速率快，雨洪过程，洼地蓄水，短历时调节中，起到削减洪峰作用[18, 19]。整个澄碧河流域近50 a年来气候变化较明显，澄碧河流域降雨量、气温和蒸发量分别呈下降、上升和下降趋势，降雨量在20世纪60年代、80年代和 21世纪发生强变异[20, 21]， 这与坝首站年最大洪峰序列呈不显著下降趋势，并于1979年发生显著突变基本一致。此外，20世纪70年代至80年代，澄碧河流域无巨大人类活动影响(如建设大型水利工程、大面积植被砍伐等)。据此，可推断影响年最大洪峰流量序列发生变异的主要原因可能是气候变化。

相关研究指出，太阳黑子活动与全球气候异常，与大气环流和热带气旋，与厄尔尼诺和拉尼娜现象，与各地区的降水异常、洪水灾害等均有良好的响应关系[22]。太阳黑子峰谷值出现时易发大洪水，从变化周期来看，岩溶区流域坝首站年最大洪峰流量存5～10、11～20和25～30 a的震荡周期，与太阳黑子的震荡周期(11、20 a)[23]有一定的吻合度。

4.2 小 结

论文利用水文变异诊断系统中的多种统计方法和Morlet小波分析法对岩溶区澄碧河水库52 a年最大洪峰流量进行变化特征分析，得到了如下结论。

(1)1963-2014年，岩溶区流域年最大洪峰流量整体呈不显著的下降趋势，年均下降率为6.04 m3/s。

(2)水文序列于1979年发生突变，其综合权重为0.375，综合显著性为3，说明52 a来，岩溶区流域年最大洪峰流量存在跳跃形式为主的变异。

(3)1985年之前，年最大洪峰流量存在15 a的主周期，8和28 a的次周期；1985年之后，不存在明显的周期变化规律。

(4)建议对年最大洪峰流量变化特征的影响因素分析时，结合气候因素深入研究，探讨1985年之后序列不存在周期性规律的原因。

(5)建议采用多个站点的数据对岩溶区流域进行系统的研究分析。