漓江最枯徑流量演變分析

羅書文，鄧亞東，覃星銘，史文強

（中國地質科學院 巖溶地質研究所／國土資源部 廣西壯族自治區巖溶動力重點實驗室，廣西 桂林541004）

漓江流域以山清水秀而著名，前來旅游的國內外游客對漓江水的清、潔、靜贊不絕口，但其枯季水量減少的問題不容忽視?？菁荆?月至次年2月）水量嚴重影響桂林漓江旅游觀光，不僅縮短航行里程甚至停航。在20世紀80年代初期，常年都可乘船從桂林市區到陽朔縣，全長約80km漓江風景盡收眼簾?，F在枯水期，只能由楊提鄉到興坪鎮10km河段可行船，且在2003年12月至2004年2月出現了記載以來的第1次斷流，致使游船停航兩月。漓江枯季徑流量減少，河水對流域內的工業廢水和生活污水的凈化稀釋能力降低，從而造成水環境不斷惡化。

近年來隨著城市規模擴大人口劇增，導致生活用水、工業用水及農業灌溉用水供需矛盾與日俱增特別是枯季尤為凸顯。因此，對于漓江最枯徑流量演變研究具有現實意義?？λ固亓饔蚍治銎淇菁緩搅髁浚芏鄬W者在方面這做了大量研究［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內有個小強區表明徑流震蕩頻繁周期短?？傮w來說桂林站年最枯徑流量序列表現出大—小—大—小的交替震蕩。

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.