漓江最枯径流量演变分析

罗书文，邓亚东，覃星铭，史文强

（中国地质科学院 岩溶地质研究所／国土资源部 广西壮族自治区岩溶动力重点实验室，广西 桂林541004）

漓江流域以山清水秀而著名，前来旅游的国内外游客对漓江水的清、洁、静赞不绝口，但其枯季水量减少的问题不容忽视。枯季（9月至次年2月）水量严重影响桂林漓江旅游观光，不仅缩短航行里程甚至停航。在20世纪80年代初期，常年都可乘船从桂林市区到阳朔县，全长约80km漓江风景尽收眼帘。现在枯水期，只能由杨提乡到兴坪镇10km河段可行船，且在2003年12月至2004年2月出现了记载以来的第1次断流，致使游船停航两月。漓江枯季径流量减少，河水对流域内的工业废水和生活污水的净化稀释能力降低，从而造成水环境不断恶化。

近年来随着城市规模扩大人口剧增，导致生活用水、工业用水及农业灌溉用水供需矛盾与日俱增特别是枯季尤为凸显。因此，对于漓江最枯径流量演变研究具有现实意义。喀斯特流域分析其枯季径流量，很多学者在方面这做了大量研究［1－11］，主要是对枯季径流影响要素、承载力、径流中长期预报和径流多时间尺度序列周期等方面研究，很少对最枯径流序列内发生的频率及各周期内的径流演变进行研究。

本研究利用Morlet小波变换法、R／S分析和P－Ⅲ型曲线对漓江流域年最枯径流量序列进行分析，旨在探讨漓江年最枯径流量在时间序列上的演变特性，为漓江生态环境保护、治理、旅游观光等社会经济活动的规划和实施提供理论依据。

1 研究区概况

研究区涉及桂林水文站（桂林市瓦窑渡头村）以上的漓江流域范围，面积2 762km2。该区主要由喀斯特和非喀斯特共同组成，主源发育于兴安县华江乡老山界南侧，地势为北高南低，主干长105km由北向南流。

分水岭地带主要为非喀斯特深山区，植被良好森林密茂，河床坡度上游极陡；中下游为丘陵、峰丛、峰林等喀斯特地貌，其河床坡度较缓。流域内多年平均降雨1 853.7mm，雨季（3—7月）占全年的67.52%，其中5月多年平均降水334.2mm占全年18.03%；枯季（9月至次年2月）降水仅占全年39.43%，其中12月多年平均降水46.9mm仅占全年2.53%。多年平均径流量1 588.7m3／s，丰水期总径流量1 094.8m3／s占全年径流量81.65%，其中6月份多年平均径流量318.1m3／s，占全年20.02%；枯水期径流总量291.5m3／s，其中12月份多年平均径流量33.6m3／s，占全年2.12%。

本研究收集了漓江流域桂林站（1942—1996年）53a最枯径流数据，数据缺失1961和1962年两年水文数据，为了保证其具有连续性，首先对缺失的原始进行拉格朗日插值法获取缺失数据（图1）。

图1 漓江桂林站1942－1996年最枯径流量变化

2 最枯径流量分析

2.1 震荡周期性分析

2.1.1 周期分析 采用小波分析方法研究了漓江最枯径流的变化特征。自从法国人Morlet提出小波分析法后很快成为国际研究热点［12－13］，其不仅能反映信

式中：t——时间参数，反映时间上相对于t的平移；α——小波周期；θ（t）——小波母函数，从公式（1）中可以看出其在实数域上的积分为零，即具有波动性。而G（φ）是θ（t）信号 的小波变换为：号在时频域上的总体特征而且能提出时域和频域的部化信息，克服了窗口傅叶氏变换带来的噪声的弱点，认为满足如下条件的任意函数为小波函数：

式中：Wω（τ，α）——ω（t）在时域τ和频域α 通过单位脉冲响应函数的输出，其基本函数有：Shannon，Gaussan、Mexican hat，Wave等函数（小波）。通过Morlet函数对漓江流域最枯年径流量序列进行分析，Morlet函数为复数函数，表示为：

式中：c——常数；i——虚部。在实数域，公式（3）的离散表达式：

式中：Δt—样本间隔，当τ较小时，对频域的分辨率低，对时域的分辨率高；当τ增大时，对频域的分辨率高，对时域的分辨率低。所以，小波变换实现了窗口的大小固定，形状可变的时频局部化［14－16］。

漓江最枯径流量序列存在显著的周期性震荡变化如图2所示。

图2 桂林站Morlet小波变换系数时频结构及方差分析

由图2可以看出，以10～15a尺度和2～35a尺度的信号最为显著，其中在以12a为中的闭合值2～－2具有较强的信号特征，以30a尺度为中心的表现次之，同时在以4和6a为中心也表现处一定的信号，所以认为漓江流域桂林站最枯径流量序列不仅存在12，30a尺度的两个主周期，同时还存在4，6a的两个次周期。由图2还可看出，1942—1950年径流量大，1950—1957年径流量小，1957—1967年径流量大，1967—1970年径流量小，1970—1972年径流量大，1972—1974年径流量小，1974—1985年径流量大；1985—1987年径流量小、而在1987年以后小波信号在1～5a内有个小强区表明径流震荡频繁周期短。总体来说桂林站年最枯径流量序列表现出大—小—大—小的交替震荡。

2.1.2 周期分析检验 为了保证小波分析结果是否具有有效性，可以运用小波方差来进一步进行检验，即对时间域上的不同尺度的所有小波系数的平方后积分如公式（6）。

其函数var（t）随时间尺度 的连续变化过程反映了信号函数中在各种尺度上的信号波动和强弱随t的变化特征。波峰处函数信号显著，说明最枯径流量序列在时间尺度上具有一个震荡的主要周期。因此，将桂林年最枯径流量1942—1996年的序列计算出来的小波系数带入公式（6）并绘制成曲线图（图2）。由图2中可以看出，其小波方差变化曲线表现出两个较为明显，位于12a左右的波峰较为明显且高而窄，而位于30a附近的峰值宽缓；同时在4，6a附近也有两个较小的波峰，且4a附近的波峰较6a对应的波峰尖。说明桂林站年最枯径流量序列存在12和30a两个主周期以及4和6a的两次周期。因此，通过小波分析桂林站年最枯径流量存在周期性变化具有一定的可靠性。

2.2 最枯径流变化趋势分析

为了进一步了解其演变趋势，运用R／S分析法研究桂林站年最枯径流量时间序列及期周期变化的趋势。通过提取各个周期径流量组成不同的径流量时间序列，建立径流序列持续性，使其反映径流时间序列前后数据之间的相互关联作用与径流序列变化趋势是具有持续性还是反持续性［17－19］，其基本模型［20］为：

对于时间序列｛X（t）｝（t＝1，2，…，n）对于任意正整数τ≥1

定义均值序列：

累积离差：

极差序列：

标准差序列：

对于比值R（τ）／S（τ）≡R／S如果存在如下关系：R／S∝τH则说明时间序列｛X（t）｝（t＝1，2，…，n），存在Hurst现象，H称为Hurst指数，H值可根据计算出的（τ，R／S）的值，在双对数坐标系〔ln（τ），ln（R／S）〕中用最小二乘法拟合，H 对应于拟合直线的斜率。根据H的大小可以判断时间序列趋势成分是表现为持续性，还是反持续性。Hurst等人证明，如果｛X（t）｝是相互独立、方差有限的随机序列，则有H＝0.5。对于不同的Hurst指数H（0＜H＜1），存在3种情况：

（1）H＝0.5时，表明时间序列变化是随机的；

（2）0＜H＜0.5时，表明时间序列具有长期相关性，但将来的总体趋势与过去的相反，过程具有反持续性。H值越接近于0，反持续性越强；

（3）0.5＜H＜1时，表明时间序列过程具有持续性，H越接近1，持续性越强。

运用Excel中的VBA编辑器工具，根据R／S分析理论编程计算桂林站1942—1996年最枯径流量时间序列以及4，6，12和30a径流量序列的Hurt系数值（分别为：桂林站 H＝0.647 209，H4a＝0.245 899，H6a＝0.654 19，H12a＝0.660 992及 H30a＝0.650 547）绘制成图（图3）。研究结果表明：（1）漓江流域桂林站最枯径流量持续减小的趋势；（2）最枯径流量周期变化中6，12，30a仍然为主要周期且12，6a逐渐加强；（3）4a变化和过去具有反持续性其信号较为强烈，由图2可知4a周期在1980年以前很少出现，所以在将来的时间里4a也将成为主要周期。

图3 桂林站最枯径流R／S分析

2.3 最枯径流量的发生频率分析

P-Ⅲ型曲线对漓江流域桂林站1942—1996年最枯径流量拟合程度较好（图4），该曲线是用于洪峰流量设计较为广泛同时也用于最枯径流量预测和计算，所以运用曲线计算每个周期最低少径流量发生频率。

图4 漓江最枯径流P－Ⅲ型拟合曲线

由图2看出，4a周期最枯径流发生在1990年（8.09m3／s），6a周 期 最 枯 径 流 发 生 在 1977 年（5.3m3／s），12a周期最枯径流发生在1966年（10.2 m3／s），30a 周 期 最 枯 径 流 发 生 在 1950 年 （3.8 m3／s）。根据P-Ⅲ型曲线预测枯季径流方法可以得出：4a周期最枯径流量8.09m3／s发生频率35%～40%，6a周期最枯径流量5.3m3／s发生频率10%～15%，12a周期最枯径流量10.2m3／s发生频率55%～60%，30a周期最枯径流量3.8m3／s发生频率1%～5%。

3 结果讨论

研究表明，在过去的几十年里，随着社会的飞速发展，漓江流域下垫面性质发生了很大改变（表1），漓江径流量的变化主要是受到下垫面性质变化所影响［15］，年最枯径流演变也不例外。漓江最枯径流量减小的直接原因有两个：（1）下垫面性质的改变是最直接最时效的原因之一；（2）研究区内由喀斯特和非喀斯特地貌单元共同组成，喀斯特地貌单元具有地表地下立体水文场系统，同时具有独有的地貌、岩性、土壤和植被等构建系统，构成了脆弱的生态环境，地表水随着地壳抬升沿着地质构造薄弱带改道或渗漏潜入地下，此过程虽然是漫长的但不能忽视。最枯径流量的周期性波动可能主要是受到气候周期性波动所制。从最枯径流量发生频率来看12a为主要周期其最枯径流量多数在10.2m3／s左右。但随着气候、环境等因素变化，使其它周期在将来也可能成为主导周期，最枯径流量变为更小，如4a周期，无论是从径流量演变趋势分析还是从周期和发生频率分析其结果都表现出较强信号，说明将成为主要周期且最枯径流量将会从10.2m3／s变为8.09m3／s为主导甚至更小。

表1 漓江上游不同生态区分类统计［15］

4 结论

（1）漓江年最枯径流存在12，30a两个主周期和4，6a两个次周期，分析了漓江年最枯径流和各周期内径流变化趋势，并可靠地计算了各周期最枯径流量发生频率值，结果可为漓江流域水资源开发利用、河流生态环境保护治理、旅游观光和抗旱救灾等提供重要的参考数据。

（2）通过运用小波分析方法使用方差检验结合频率计算发现，小波方差曲线波峰尖而窄表示最枯径流量发生的频率较高；波峰宽而缓则表示发生的频率较低。最枯径流量发生频率与振幅无关。

（3）小波和R／S分析方法对漓江径流研究具有一定的适用性，但并不能说明对于所有的喀斯特流域径流研究均适用。喀斯特水文过程与非喀斯特流域更是大相径庭［18］，表现出极强的非线性和随机性特性［11］使得其水系发育、水文动态上表现出其独有性，加上人类活动的干预下变得尤为复杂。所以研究喀斯特流域径流演变情况还需长期观测收集更多数据，在研究方法上还有待于进一步提高。

［1］ 王在高，梁虹.基于GIS分析喀斯特流域下垫面因素对枯季径流的影响：以贵州省河流为例［J］.中国岩溶，2002，21（1）：55-60.

［2］ 孔兰，梁虹，黄法苏，等.基于喀斯特流域径流量多时间尺度小波分析［J］.人民长江，2008，39（5）：17-26.

［3］ 罗书文，梁虹，杨桃，等.南明河流域枯水径流量的长期预报［J］.水土保持通报，2008，28（5）：44-47.

［4］ 郝庆庆，陈喜，马建良.南方喀斯特流域枯季退水影响因子分析［J］.水土保持研究，2009，16（6）：22-29.

［5］ 孔兰，梁虹，戴洪刚.基于灰色关联法的喀斯特流域枯水影响因素分析［J］.水科学与工程技术，2007（4）：1-3.

［6］ 谢永玉，石朋，瞿思敏，等.岩溶流域枯季径流的区域频率分析［J］.水电能源科学，2012，30（6）：24-27.

［7］ 罗书文，梁虹，杨桃，等.基于分形理论的喀斯特流域枯水径流影响因素分析［J］.水科学与工程技术，2008（5）：44-46.

［8］ 梁虹，王在高.喀斯特流域枯水径流频率分析：以贵州省河流为例［J］.中国岩溶，2002，21（2）：106-113.

［9］ 孔兰，梁虹，黄法苏.喀斯特流域径流量时序演变特征分析：以贵州省为例［J］.中国岩溶，2007，26（4）：341-346.

［10］ 梁虹，王剑.喀斯特地区流域岩性差异与洪、枯水特征值相关分析：以贵州河流为例［J］.中国岩溶，1998，17（1）：67-73.

［11］ 梁虹.喀斯特流域尺度与枯水流量初步研究：以贵州为例［J］.贵州师范大学学报：自然科学版，1997，15（3）：1-5.

［12］ 王文圣，丁晶，向红莲.小波分析在水文学中的应用研究及展望［J］.水科学进展，2002，13（4）：515-517.

［13］ 王文圣，丁晶，向红莲.水文时间序列多时间尺度分析的小波变换法［J］.四川大学学报：工程科学版，2002，34（6）：14-17.

［14］ 张学真，刘燕.灞河出山径流序列变化的小波分析［J］.水资源保护，2006，22（3）：12-15.

［15］ 郭纯青，方荣杰，代俊峰，等.漓江流域上游区水资源与水环境演变及预测［M］.北京：中国水利水电出版社，2011.

［16］ 张超，杨秉赓.计量地理学［M］.北京：高等教育出版社，2004.

［17］ 张少文，丁晶，廖杰，等.基于小波的黄河上游天然年径流变化特性分析［J］.四川大学学报：工程科学版，2004，36（3）：32-37.

［18］ 杨明德，谭明，梁虹.喀斯特流域水文地貌系统［M］.北京：地质出版社，1998.

［19］ 黄勇，周志芳，王锦国，等.R／S分析法在地下水动态分析中的应用［J］.河海大学学报，2002，30（1）：83-87.

［20］ 徐建华.现代地理学中的数学方法［M］.北京：高等教育出版社，2002.