漓江最枯徑流量演變分析

羅書文，鄧亞東，覃星銘，史文強(qiáng)

（中國地質(zhì)科學(xué)院 巖溶地質(zhì)研究所／國土資源部 廣西壯族自治區(qū)巖溶動(dòng)力重點(diǎn)實(shí)驗(yàn)室，廣西 桂林541004）

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

近年來隨著城市規(guī)模擴(kuò)大人口劇增，導(dǎo)致生活用水、工業(yè)用水及農(nóng)業(yè)灌溉用水供需矛盾與日俱增特別是枯季尤為凸顯。因此，對(duì)于漓江最枯徑流量演變研究具有現(xiàn)實(shí)意義?？λ固亓饔蚍治銎淇菁緩搅髁?，很多學(xué)者在方面這做了大量研究［1－11］，主要是對(duì)枯季徑流影響要素、承載力、徑流中長期預(yù)報(bào)和徑流多時(shí)間尺度序列周期等方面研究，很少對(duì)最枯徑流序列內(nèi)發(fā)生的頻率及各周期內(nèi)的徑流演變進(jìn)行研究。

本研究利用Morlet小波變換法、R／S分析和P－Ⅲ型曲線對(duì)漓江流域年最枯徑流量序列進(jìn)行分析，旨在探討漓江年最枯徑流量在時(shí)間序列上的演變特性，為漓江生態(tài)環(huán)境保護(hù)、治理、旅游觀光等社會(huì)經(jīng)濟(jì)活動(dòng)的規(guī)劃和實(shí)施提供理論依據(jù)。

1 研究區(qū)概況

研究區(qū)涉及桂林水文站（桂林市瓦窯渡頭村）以上的漓江流域范圍，面積2 762km2。該區(qū)主要由喀斯特和非喀斯特共同組成，主源發(fā)育于興安縣華江鄉(xiāng)老山界南側(cè)，地勢為北高南低，主干長105km由北向南流。

分水嶺地帶主要為非喀斯特深山區(qū)，植被良好森林密茂，河床坡度上游極陡；中下游為丘陵、峰叢、峰林等喀斯特地貌，其河床坡度較緩。流域內(nèi)多年平均降雨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最枯徑流數(shù)據(jù)，數(shù)據(jù)缺失1961和1962年兩年水文數(shù)據(jù)，為了保證其具有連續(xù)性，首先對(duì)缺失的原始進(jìn)行拉格朗日插值法獲取缺失數(shù)據(jù)（圖1）。

圖1 漓江桂林站1942－1996年最枯徑流量變化

2 最枯徑流量分析

2.1 震蕩周期性分析

2.1.1 周期分析 采用小波分析方法研究了漓江最枯徑流的變化特征。自從法國人Morlet提出小波分析法后很快成為國際研究熱點(diǎn)［12－13］，其不僅能反映信

式中：t——時(shí)間參數(shù)，反映時(shí)間上相對(duì)于t的平移；α——小波周期；θ（t）——小波母函數(shù)，從公式（1）中可以看出其在實(shí)數(shù)域上的積分為零，即具有波動(dòng)性。而G（φ）是θ（t）信號(hào) 的小波變換為：號(hào)在時(shí)頻域上的總體特征而且能提出時(shí)域和頻域的部化信息，克服了窗口傅葉氏變換帶來的噪聲的弱點(diǎn)，認(rèn)為滿足如下條件的任意函數(shù)為小波函數(shù)：

式中：Wω（τ，α）——ω（t）在時(shí)域τ和頻域α 通過單位脈沖響應(yīng)函數(shù)的輸出，其基本函數(shù)有：Shannon，Gaussan、Mexican hat，Wave等函數(shù)（小波）。通過Morlet函數(shù)對(duì)漓江流域最枯年徑流量序列進(jìn)行分析，Morlet函數(shù)為復(fù)數(shù)函數(shù)，表示為：

式中：c——常數(shù)；i——虛部。在實(shí)數(shù)域，公式（3）的離散表達(dá)式：

式中：Δt—樣本間隔，當(dāng)τ較小時(shí)，對(duì)頻域的分辨率低，對(duì)時(shí)域的分辨率高；當(dāng)τ增大時(shí)，對(duì)頻域的分辨率高，對(duì)時(shí)域的分辨率低。所以，小波變換實(shí)現(xiàn)了窗口的大小固定，形狀可變的時(shí)頻局部化［14－16］。

漓江最枯徑流量序列存在顯著的周期性震蕩變化如圖2所示。

圖2 桂林站Morlet小波變換系數(shù)時(shí)頻結(jié)構(gòu)及方差分析

由圖2可以看出，以10～15a尺度和2～35a尺度的信號(hào)最為顯著，其中在以12a為中的閉合值2～－2具有較強(qiáng)的信號(hào)特征，以30a尺度為中心的表現(xiàn)次之，同時(shí)在以4和6a為中心也表現(xiàn)處一定的信號(hào)，所以認(rèn)為漓江流域桂林站最枯徑流量序列不僅存在12，30a尺度的兩個(gè)主周期，同時(shí)還存在4，6a的兩個(gè)次周期。由圖2還可看出，1942—1950年徑流量大，1950—1957年徑流量小，1957—1967年徑流量大，1967—1970年徑流量小，1970—1972年徑流量大，1972—1974年徑流量小，1974—1985年徑流量大；1985—1987年徑流量小、而在1987年以后小波信號(hào)在1～5a內(nèi)有個(gè)小強(qiáng)區(qū)表明徑流震蕩頻繁周期短?？傮w來說桂林站年最枯徑流量序列表現(xiàn)出大—小—大—小的交替震蕩。

2.1.2 周期分析檢驗(yàn) 為了保證小波分析結(jié)果是否具有有效性，可以運(yùn)用小波方差來進(jìn)一步進(jìn)行檢驗(yàn)，即對(duì)時(shí)間域上的不同尺度的所有小波系數(shù)的平方后積分如公式（6）。

其函數(shù)var（t）隨時(shí)間尺度 的連續(xù)變化過程反映了信號(hào)函數(shù)中在各種尺度上的信號(hào)波動(dòng)和強(qiáng)弱隨t的變化特征。波峰處函數(shù)信號(hào)顯著，說明最枯徑流量序列在時(shí)間尺度上具有一個(gè)震蕩的主要周期。因此，將桂林年最枯徑流量1942—1996年的序列計(jì)算出來的小波系數(shù)帶入公式（6）并繪制成曲線圖（圖2）。由圖2中可以看出，其小波方差變化曲線表現(xiàn)出兩個(gè)較為明顯，位于12a左右的波峰較為明顯且高而窄，而位于30a附近的峰值寬緩；同時(shí)在4，6a附近也有兩個(gè)較小的波峰，且4a附近的波峰較6a對(duì)應(yīng)的波峰尖。說明桂林站年最枯徑流量序列存在12和30a兩個(gè)主周期以及4和6a的兩次周期。因此，通過小波分析桂林站年最枯徑流量存在周期性變化具有一定的可靠性。

2.2 最枯徑流變化趨勢分析

為了進(jìn)一步了解其演變趨勢，運(yùn)用R／S分析法研究桂林站年最枯徑流量時(shí)間序列及期周期變化的趨勢。通過提取各個(gè)周期徑流量組成不同的徑流量時(shí)間序列，建立徑流序列持續(xù)性，使其反映徑流時(shí)間序列前后數(shù)據(jù)之間的相互關(guān)聯(lián)作用與徑流序列變化趨勢是具有持續(xù)性還是反持續(xù)性［17－19］，其基本模型［20］為：

對(duì)于時(shí)間序列｛X（t）｝（t＝1，2，…，n）對(duì)于任意正整數(shù)τ≥1

定義均值序列：

累積離差：

極差序列：

標(biāo)準(zhǔn)差序列：

對(duì)于比值R（τ）／S（τ）≡R／S如果存在如下關(guān)系：R／S∝τH則說明時(shí)間序列｛X（t）｝（t＝1，2，…，n），存在Hurst現(xiàn)象，H稱為Hurst指數(shù)，H值可根據(jù)計(jì)算出的（τ，R／S）的值，在雙對(duì)數(shù)坐標(biāo)系〔ln（τ），ln（R／S）〕中用最小二乘法擬合，H 對(duì)應(yīng)于擬合直線的斜率。根據(jù)H的大小可以判斷時(shí)間序列趨勢成分是表現(xiàn)為持續(xù)性，還是反持續(xù)性。Hurst等人證明，如果｛X（t）｝是相互獨(dú)立、方差有限的隨機(jī)序列，則有H＝0.5。對(duì)于不同的Hurst指數(shù)H（0＜H＜1），存在3種情況：

（1）H＝0.5時(shí)，表明時(shí)間序列變化是隨機(jī)的；

（2）0＜H＜0.5時(shí)，表明時(shí)間序列具有長期相關(guān)性，但將來的總體趨勢與過去的相反，過程具有反持續(xù)性。H值越接近于0，反持續(xù)性越強(qiáng)；

（3）0.5＜H＜1時(shí)，表明時(shí)間序列過程具有持續(xù)性，H越接近1，持續(xù)性越強(qiáng)。

運(yùn)用Excel中的VBA編輯器工具，根據(jù)R／S分析理論編程計(jì)算桂林站1942—1996年最枯徑流量時(shí)間序列以及4，6，12和30a徑流量序列的Hurt系數(shù)值（分別為：桂林站 H＝0.647 209，H4a＝0.245 899，H6a＝0.654 19，H12a＝0.660 992及 H30a＝0.650 547）繪制成圖（圖3）。研究結(jié)果表明：（1）漓江流域桂林站最枯徑流量持續(xù)減小的趨勢；（2）最枯徑流量周期變化中6，12，30a仍然為主要周期且12，6a逐漸加強(qiáng)；（3）4a變化和過去具有反持續(xù)性其信號(hào)較為強(qiáng)烈，由圖2可知4a周期在1980年以前很少出現(xiàn)，所以在將來的時(shí)間里4a也將成為主要周期。

圖3 桂林站最枯徑流R／S分析

2.3 最枯徑流量的發(fā)生頻率分析

P-Ⅲ型曲線對(duì)漓江流域桂林站1942—1996年最枯徑流量擬合程度較好（圖4），該曲線是用于洪峰流量設(shè)計(jì)較為廣泛同時(shí)也用于最枯徑流量預(yù)測和計(jì)算，所以運(yùn)用曲線計(jì)算每個(gè)周期最低少徑流量發(fā)生頻率。

圖4 漓江最枯徑流P－Ⅲ型擬合曲線

由圖2看出，4a周期最枯徑流發(fā)生在1990年（8.09m3／s），6a周 期 最 枯 徑 流 發(fā) 生 在 1977 年（5.3m3／s），12a周期最枯徑流發(fā)生在1966年（10.2 m3／s），30a 周 期 最 枯 徑 流 發(fā) 生 在 1950 年 （3.8 m3／s）。根據(jù)P-Ⅲ型曲線預(yù)測枯季徑流方法可以得出：4a周期最枯徑流量8.09m3／s發(fā)生頻率35%～40%，6a周期最枯徑流量5.3m3／s發(fā)生頻率10%～15%，12a周期最枯徑流量10.2m3／s發(fā)生頻率55%～60%，30a周期最枯徑流量3.8m3／s發(fā)生頻率1%～5%。

3 結(jié)果討論

研究表明，在過去的幾十年里，隨著社會(huì)的飛速發(fā)展，漓江流域下墊面性質(zhì)發(fā)生了很大改變（表1），漓江徑流量的變化主要是受到下墊面性質(zhì)變化所影響［15］，年最枯徑流演變也不例外。漓江最枯徑流量減小的直接原因有兩個(gè)：（1）下墊面性質(zhì)的改變是最直接最時(shí)效的原因之一；（2）研究區(qū)內(nèi)由喀斯特和非喀斯特地貌單元共同組成，喀斯特地貌單元具有地表地下立體水文場系統(tǒng)，同時(shí)具有獨(dú)有的地貌、巖性、土壤和植被等構(gòu)建系統(tǒng)，構(gòu)成了脆弱的生態(tài)環(huán)境，地表水隨著地殼抬升沿著地質(zhì)構(gòu)造薄弱帶改道或滲漏潛入地下，此過程雖然是漫長的但不能忽視。最枯徑流量的周期性波動(dòng)可能主要是受到氣候周期性波動(dòng)所制。從最枯徑流量發(fā)生頻率來看12a為主要周期其最枯徑流量多數(shù)在10.2m3／s左右。但隨著氣候、環(huán)境等因素變化，使其它周期在將來也可能成為主導(dǎo)周期，最枯徑流量變?yōu)楦?，?a周期，無論是從徑流量演變趨勢分析還是從周期和發(fā)生頻率分析其結(jié)果都表現(xiàn)出較強(qiáng)信號(hào)，說明將成為主要周期且最枯徑流量將會(huì)從10.2m3／s變?yōu)?.09m3／s為主導(dǎo)甚至更小。

表1 漓江上游不同生態(tài)區(qū)分類統(tǒng)計(jì)［15］

4 結(jié)論

（1）漓江年最枯徑流存在12，30a兩個(gè)主周期和4，6a兩個(gè)次周期，分析了漓江年最枯徑流和各周期內(nèi)徑流變化趨勢，并可靠地計(jì)算了各周期最枯徑流量發(fā)生頻率值，結(jié)果可為漓江流域水資源開發(fā)利用、河流生態(tài)環(huán)境保護(hù)治理、旅游觀光和抗旱救災(zāi)等提供重要的參考數(shù)據(jù)。

（2）通過運(yùn)用小波分析方法使用方差檢驗(yàn)結(jié)合頻率計(jì)算發(fā)現(xiàn)，小波方差曲線波峰尖而窄表示最枯徑流量發(fā)生的頻率較高；波峰寬而緩則表示發(fā)生的頻率較低。最枯徑流量發(fā)生頻率與振幅無關(guān)。

（3）小波和R／S分析方法對(duì)漓江徑流研究具有一定的適用性，但并不能說明對(duì)于所有的喀斯特流域徑流研究均適用?？λ固厮倪^程與非喀斯特流域更是大相徑庭［18］，表現(xiàn)出極強(qiáng)的非線性和隨機(jī)性特性［11］使得其水系發(fā)育、水文動(dòng)態(tài)上表現(xiàn)出其獨(dú)有性，加上人類活動(dòng)的干預(yù)下變得尤為復(fù)雜。所以研究喀斯特流域徑流演變情況還需長期觀測收集更多數(shù)據(jù)，在研究方法上還有待于進(jìn)一步提高。

［1］ 王在高，梁虹.基于GIS分析喀斯特流域下墊面因素對(duì)枯季徑流的影響：以貴州省河流為例［J］.中國巖溶，2002，21（1）：55-60.

［2］ 孔蘭，梁虹，黃法蘇，等.基于喀斯特流域徑流量多時(shí)間尺度小波分析［J］.人民長江，2008，39（5）：17-26.

［3］ 羅書文，梁虹，楊桃，等.南明河流域枯水徑流量的長期預(yù)報(bào)［J］.水土保持通報(bào)，2008，28（5）：44-47.

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

［5］ 孔蘭，梁虹，戴洪剛.基于灰色關(guān)聯(lián)法的喀斯特流域枯水影響因素分析［J］.水科學(xué)與工程技術(shù)，2007（4）：1-3.

［6］ 謝永玉，石朋，瞿思敏，等.巖溶流域枯季徑流的區(qū)域頻率分析［J］.水電能源科學(xué)，2012，30（6）：24-27.

［7］ 羅書文，梁虹，楊桃，等.基于分形理論的喀斯特流域枯水徑流影響因素分析［J］.水科學(xué)與工程技術(shù)，2008（5）：44-46.

［8］ 梁虹，王在高.喀斯特流域枯水徑流頻率分析：以貴州省河流為例［J］.中國巖溶，2002，21（2）：106-113.

［9］ 孔蘭，梁虹，黃法蘇.喀斯特流域徑流量時(shí)序演變特征分析：以貴州省為例［J］.中國巖溶，2007，26（4）：341-346.

［10］ 梁虹，王劍.喀斯特地區(qū)流域巖性差異與洪、枯水特征值相關(guān)分析：以貴州河流為例［J］.中國巖溶，1998，17（1）：67-73.

［11］ 梁虹.喀斯特流域尺度與枯水流量初步研究：以貴州為例［J］.貴州師范大學(xué)學(xué)報(bào)：自然科學(xué)版，1997，15（3）：1-5.

［12］ 王文圣，丁晶，向紅蓮.小波分析在水文學(xué)中的應(yīng)用研究及展望［J］.水科學(xué)進(jìn)展，2002，13（4）：515-517.

［13］ 王文圣，丁晶，向紅蓮.水文時(shí)間序列多時(shí)間尺度分析的小波變換法［J］.四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2002，34（6）：14-17.

［14］ 張學(xué)真，劉燕.灞河出山徑流序列變化的小波分析［J］.水資源保護(hù)，2006，22（3）：12-15.

［15］ 郭純青，方榮杰，代俊峰，等.漓江流域上游區(qū)水資源與水環(huán)境演變及預(yù)測［M］.北京：中國水利水電出版社，2011.

［16］ 張超，楊秉賡.計(jì)量地理學(xué)［M］.北京：高等教育出版社，2004.

［17］ 張少文，丁晶，廖杰，等.基于小波的黃河上游天然年徑流變化特性分析［J］.四川大學(xué)學(xué)報(bào)：工程科學(xué)版，2004，36（3）：32-37.

［18］ 楊明德，譚明，梁虹.喀斯特流域水文地貌系統(tǒng)［M］.北京：地質(zhì)出版社，1998.

［19］ 黃勇，周志芳，王錦國，等.R／S分析法在地下水動(dòng)態(tài)分析中的應(yīng)用［J］.河海大學(xué)學(xué)報(bào)，2002，30（1）：83-87.

［20］ 徐建華.現(xiàn)代地理學(xué)中的數(shù)學(xué)方法［M］.北京：高等教育出版社，2002.