張 勇,吳 福,翟國軍,鄧洪星
(廣西壯族自治區地質環境監測站,南寧530029)
對于水文要素的序列研究,主要集中在時域和頻域兩方向上,然而水文要素受到多種自然因素及其人為活動的影響,其序列在發展過程中不僅具有一定的趨勢性和周期性,還受各種因素的影響伴有隨機性和突變性,其本身并非平穩序列。傳統的統計方法在分析水文要素的時間序列時域和頻域特征上具有很大的局限性,不能揭示水文要素在變化上多尺度和多層次的結構[1]。而小波分析是一種對數據進行調和分析的方法,不僅沿襲了傳統譜分析方法中時頻聯合分析的優越性,能有效解決了時域和頻域在傳統分析中所出現的局部化問題,更能解決水文要素長時間序列所涉及的突變點及多時間尺度特征等問題[2-4]。
珠江流域流經滇、黔、桂等63個地(州)市,關于其水文要素的研究成果較為豐富:戴仕寶等通過分析珠江流域主要水文站的水文要素資料探討輸沙量與氣候變化的關系[5];陳立華選取西江中上游的天峨站、遷江站和梧州站等3 個水文站近30年的逐月徑流資料研究了徑流演變趨勢[6];吳創收利用小波分析和Mann-Kendall 非參數秩次檢驗兩種方法對珠江流域入海水沙序列通量變化特點進行分析[7]。目前珠江流域自然要素和人類活動加劇對氣候和流域下墊面環境產生直接影響,而氣候變化和下墊面改變進而影響流域的產匯流和侵蝕產沙,最終可能致使原有的水文過程發生在其顯著的變化。其主要表現為流域水文過程存在發生異變的點(跳躍點),這個點指的是時間和空間上的節點。水文時間序列在過程和特性上已經發生了巨大的改變,以上工作從不同角度得到了很多有意義的結論,但主要集中在較大的流域尺度上;并且以小流域為單元的小波分析時間序列變化研究,能為綜合研究大流域的水沙演變提供參考,仍然具有一定的科學意義。而小流域的水沙變化特性的耦合關系研究可以探尋兩者的發展態勢,從而能較好地反映兩系統的相互脅迫作用,為地區發展起一定的指導作用。
20世紀50年代以來,受到人為活動(主要為大煉鋼鐵運動和修建水庫)、氣候波動及水土保持等因素的影響,西江水沙序列變化過程更加復雜。基于此,本次將選取西江中游缺乏水沙序列耦合分析的遷江站作為研究對象,采用Mann-kendall 檢驗法進行水沙時間序列的趨勢性變化分析,可得西江中游水沙序列的變化趨勢,進一步采用小波分析方法對遷江站水沙序列進行周期特征分析,獲取較長的時序長度進行周期識別,得到水沙序列的變化規律,并采用耦合的方法得到徑流和輸沙序列的發展特性,可以獲得西江中游水沙序列的演化特征,為西江流域水沙序列的多時間尺度特征提供工程參考依據。
西江是珠江流域的主要干流,位于22°15′~26°30′N 和106°42′~112°30′E 之間,流域面積33 萬km2。其發源于云南省曲靖市,地處低緯度地區,北回歸線橫貫中部,屬于亞熱帶季風氣候區域,多年平均氣溫介于19.7~21.3 ℃之間,年際變化不大;多年平均年降水量一般在1 537.1 mm 左右。西江自西向東流經滇、黔、桂和粵等省區,下游與珠江水網相接,一直流到港澳,是構成泛珠江三角區域經濟體系的重要海動脈。
遷江站位于西江中游,是西江干流紅水河的主要控制站,其流域控制面積137 760 km2,汛期在4-9月,占年徑流量的78.2%。該站建站以來,為紅水河的防汛抗旱,給排水和水利建設等部門提供各種寶貴資料。本研究利用《中華人民共和國水利部刊發的中國河流泥沙公報》和梧州市水利局所提供的1957-2016年遷江站水沙序列數據進行研究分析,該數據資料具有較好的代表性。
1.2.1 Mann-Kendall趨勢檢驗
Mann-Kendall 趨勢檢驗是研究長時間序列趨勢變化的有效工具,其優點是所受到少數異常數據的干擾較小,不需要要求樣本自身具有一定的分布規律,被廣泛用于分析水文序列趨勢變化特征的研究[8-10]。假設某一水文序列X(x1,x2,…,xn),所涉及的時間序列長度為n,水文要素統計變量S的具體計算方法可由下方方程表示:

式中:sgn為sign函數;xm為序列第m年的數值;xj為序列第j年的數值,且j>m。
指定Z為標準統計值,計算方法如下:

當水文要素統計量Z值是正的,說明該水文序列具有增加趨勢,若Z是負值,說明該序列具有下降的趨勢。經查表可知,當水文要素的顯著性水平為0.05 時,所對應的統計檢驗臨界值為Za/2=±1.96,若|Z|大于|Za/2|=1.96,則表現為明顯著,反之則為不顯著;同理,當顯著性水平是0.01 時,水文要素的統計檢驗臨界值變為Za/2=±2.58,|Z|大于|Za/2|=2.58,則表現為明顯著,反之則為不顯著。
1.2.2 小波分析
水文系統中水沙變化受到氣候、地形和人類活動等綜合影響,其在演化過程中具備多時間尺度和多空間尺度重要特征。小波分析采用一簇小波函數來表述某一信號的基本思想,能較好地解決多尺度所帶來的多層次問題。因此結合水文序列的特性,本文采用Morlet連續復小波函數[11-13],分析探討遷江站的水沙演化規律,其函數為:

式中:w0為常數;t為虛部。
對于已知并滿足一定條件的小波函數φ(t),時間序列f(t) ∈L2(R)的連續小波變換為:

式中:?φ(t)為φ(t)的復共軛函數;φ(t)為基小波;a為對應的尺度因子;b為平移因子;Wf(a,b)為對應于不同尺度a下的不同位置b的小波系數。一般時間序列以離散的形式存在,則式(2)可表示為:

式中:N為函數的離散數;?t為抽樣的時間間隔;k為抽樣時間間隔的數目;Wf(a,b)則能夠同時反饋時頻參數b和頻域參數a所具有的特性。
小波方差是反饋波動的能量隨著尺度的分布,用來確定水沙序列的主要周期,其公式為:

1.2.3 水沙耦合模型
耦合是指多個(含兩個)系統相互影響和相互作用的過程,通過借鑒物理學的演化思想[14-16],可構建某個系統的非線性演化方程式,對該系統的發展進行描述:

式中:f為Xi的非線性函數。
在原點處對非線性函數進行泰勒級數展開,利用利亞普斯諾夫第一近似定理[17],保證改系統穩定性,從而得到近似線性系統:

按上述方法建立徑流量(R)與輸沙量(S)這兩個水文要素系統的一般函數為:

式中:x為徑流系統的元素(以時間為變量的函數);y為輸沙系統的元素(以時間為變量的函數);ai為徑流系統中元素xi對應權重;bi為輸沙系統中元素yi對應權重。
由于西江中游的徑流量與輸沙量系統相互影響與作用,將兩者視作一個復合系統,f(R)與f(S)在系統中占處于主要地位,其演化方程可表征為:

式中:A為徑流系統本身和多種外界因素影響下的演化狀態;B為輸沙系統本身和多種外界因素影響下的演化狀態;VA為徑流系統的演化速度;VB為輸沙系統的演化速度。
由于整個系統只有f(R)與f(S)這兩個元素,所以當f(R)與f(S)協調時,代表水沙序列這一系統的運行也是協調的,其演化速度V為關于徑流量(VA)和輸沙量(VB)的非線性函數,所以V=f(VA,VB)。通過控制VA,VB的時間變量,對V進行分析,可以獲知整個系統及f(R)與f(S)的協調關系。
首先建立滿足S 型發展機制的V 模型,通過將VA,VB的演化軌跡投影在直角坐標系中進行分析。將代表徑流系統和輸沙系統演化速度的VA,VB的比值反正切函數,a為年徑流量與年輸沙量系統的耦合度,也是復合系整體統的主要演變過程(圖1),大致可分為3個階段:

圖1 系統演化軌跡
當0°≤α<30°時,系統處于初級發展時期,流域徑流量和輸沙量受到人類活動影響較小,氣候環境穩定,徑流量與輸沙量兩個系統之間的影響與限制關系不強,徑流量與輸沙量協調發展。
當30°≤α<60°時,系統進入極限發展時期。此時,人類活動對流域的影響加劇,加之氣候環境的變遷,原本和諧的徑流量和輸沙量系統的限制加劇。人類活動對自然環境的改造會直接影響流域土地利用方式的改變,導致輸沙量系統改變,加之氣候因子的變化,徑流系統也隨之改變,兩個系統之間的矛盾日益嚴重。此時,為保證兩系統的不斷進步,人們將采取一系列的措施對出現的矛盾進行緩解,科學技術的進步促使兩系統的耦合朝著良性階段進行發展,從而進入一個更為的高級關系階段。
當60°≤α<90°時,原有劇烈矛盾的系統在相對正確的措施作用下,會進入再生時期。此時,科學技術的運用會緩解乃至解決日益嚴重的矛盾,兩個系統之間的關系會逐步和諧,最終達到較高的水平,但新的問題也將會隨之而來,新的矛盾將會出現。
1.2.4 數據處理

式中:xij為第i年的第j個指標的初始值;Zij為經過標準化處理后相對應的數據[18]。
根據遷江站(1957-2016年)實測年徑流量的數據,做出該水文要素的年徑流量變化過程線、均值線和一元方程趨勢線,如圖2。近60年來遷江站的多年徑流量平均值為646.6 億m3,最大年徑流量和最小年徑流量分別出現于1979年和2013年,與之相對應的流量分別為1 043 億m3和377.3 億m3,大部分年份的徑流量圍繞于趨勢線上下波動,多數年份的徑流量均高于多年平均值,多年徑流量總體上呈現下降趨勢。

圖2 遷江站年流量序列波動變化
而根據趨勢線判斷,因此結合表1,即Mann-Kendall趨勢檢驗法的計算結果來判斷遷江站的變化顯著情況,經計算遷江站的標準統計量值Z為正值,且|Z|大于|Za/2|=1.96,表明西江中游遷江站的年徑流量下降的變化趨勢顯著。

表1 1957-2016西江流域干流年徑流變化趨勢的Mann-kendall檢驗結果
為了解遷江站輸沙量的演變特征,繪制遷江站年輸沙量序列過程線,圖中輸沙量過程曲線在20世紀90年代后期顯著下降(圖3)。對于年徑流量的變化趨勢,20世紀90年代之前,年輸沙量與徑流量呈波動變化具有一定的相關性,總體上均為豐水多沙,少水少沙的規律;而20世紀90年代后,輸沙量的波動與徑流量的波動并無相關性。

圖3 西江流域遷江站年輸沙序列波動變化

表2 西江流域干流年輸沙量變化趨勢的Mann-kendall檢驗結果
用Mann-kendall趨勢檢驗法計算出1957-2016年遷江站的年輸沙量統計量值Z(見表4),其年輸沙量統計值|Z|均大于置信度α等于0.05 時的臨界值,故西江流域干流遷江站的年輸沙量變化趨勢均為顯著;且統計量值Z均為負數,因此年輸沙量的下降趨勢。
根據1957-2016年遷江站水沙資料,繪制水沙序列的小波系數實部等值線圖(圖4),規避系數虛部所可能帶來的誤差,進一步知道年徑流量在不同尺度的豐、枯交替情況,以及突變點的分布。圖4(a)中在流域徑流序列演變過程中,有25~30 a,13~24 a 和3~12 a 的3 類周期變化。其中,20~30 a 與13~20 a 是較明顯的兩個主要周期,在20~30 a 尺度上出現了3 次震蕩,在13~20 a時間尺度上徑流序列存在著6次明顯的豐枯交替震蕩;圖4(b)中可知,遷江站輸沙序列演變過程中,主要存在20~30 a,9~20 a 和3~9 a 的三類周期的明顯變化。其中,20~30 a 是遷江站比較顯著的主周期,在1957-2016年期間經歷了4 次輸沙量震蕩,變化較為穩定,具有全局性。9~20 a 和3~9 a 時間尺度是遷江站中輸沙量變化相對復雜的周期,存在一定的輸沙增減交替變化,在9~20 a 的時間尺度下,1957-1995年內發生了4 次震蕩,而1995年后,輸沙量的變化趨于穩定;在3~9 a 時間尺度下,1957-2000年內,發生了10次震蕩,在2003年后輸沙量的變化趨于穩定。圖4中兩幅在正負相位的交替上出現相似的地方,在20~30 a時間尺度均發生3次震蕩,震蕩時頻幾乎一致,但其相位交替出現相反的規律。

圖4 遷江站水沙序列的小波系數實部等值線圖對比圖
小波系數的模方圖作為小波能量譜(圖5),可以很好地反映震蕩能量[10],知道水沙在各時間尺度的具有不同強弱分布。圖5(a)中遷江站徑流序列的20~30 a 時間尺度變化最強,在整個時域中處于主導地位,震蕩中心位于2013年;13~20 a 和3~12 a時間尺度震蕩強度雖然弱于20~30 a時間尺度的震蕩,但仍占據了整個研究時頻,它們的震蕩中心分別為1960年和2015年。圖5(b)中在遷江站年輸沙序列20~30 a時間尺度震蕩能量較強,其周期占據整個研究時域,震蕩中心在1978年左右;而其在9~20 a 和3~9 a 時間尺度能量較弱,輸沙序列的震蕩中心分別在1967年和1983年左右,震蕩能量分別于1995年左右和2003年左右開始衰弱。

圖5 遷江站水沙序列模方等值線圖
為進一步了解遷江站水沙序列的周期顯著性隨時間尺度的變化特征,采用式(4)繪制水沙序列小波方差檢驗圖(圖6),方差曲線中每一峰值可以對應每一時間尺度下存在的顯著周期。從圖6可以看出,遷江站水沙序列的小波峰的出現呈現一定的相似規律,但主波峰的出現存在一定的差異。根據小波方差圖,遷江站徑流序列的第一主周期是28 a,而輸沙序列的主峰值變成了32 a;水沙序列的第二主周期均為16 a;徑流序列第三主周期出現在5 a,輸沙序列的則出現在了6 a上,大致相同。造成水沙序列不同主周期的差異,可能是由于受到人類活動或外界因素的影響。

圖6 遷江站徑水沙序列的小波方差圖
小波系數實部變化過程(圖7)能很好地反映水文要素在該時間尺度下的小波系數實部變化特征。其中,正值小波系數實部代表水沙偏多期,負值代表水沙偏少期,而系數為0時則為突變點。

圖7 遷江站水沙序列在不同時間尺度下的小波系數實部變化曲線圖
在5 a時間尺度上,遷江站徑流序列大約發生了36個枯-豐周期轉換變化,平均一個周期的變化在4 a 左右;輸沙序列大約經歷了30 個偏多期-偏少期,平均變化周期在4.5 a,2003年以后輸沙序列的變化趨于平穩。由于水沙序列的小波系數實部曲線的振幅大小與其量值變化程度有關,則遷江站在這一段時間徑流與輸沙這兩個水文要素的量值波動性較大。20世紀50-60年代,小尺度上,水沙序列變化趨勢基本同步,輸沙序列振幅要小于徑流序列。70年代期間,由于流域自身條件及不合理的開發導致了水沙序列發生紊亂,期間,農業和經濟的改革導致毀林開荒造田情況嚴重,水土流失情況進一步惡化[19];70年代中期以后,兩條曲線的變化規律越發不一樣,輸沙曲線振幅已經超過徑流曲線,大規模開發所帶來的生態環境等問題逐步引起了人們的重視[20]。1991年《中華人民共和國水土保持法》的通過,人工造林、封山育林等一系列造林措施的實施有效抑制了中國生態環境的惡化,在法案實施后,截止2011年末,廣西壯族自治區在水土流失面積從30 600 km2減少至28 122 km2[21];1999年所啟動的“退耕還林工程”,宣告紅水河流域的大規模開墾時代結束[20]。這些政策和措施的實施都有助于輸沙量的減少。
在16 a時間尺度上,遷江站水沙序列大約經歷了11個豐枯的周期轉換變化,平均一個周期的變化在13 a 左右。20世紀50-60年代期間,紅水河流域沒有水利設施的興建,水沙變化較為一致。在70-90年代期間,龍灘等水利樞紐的先后建成,對水沙起到了調配的作用,使得輸沙序列曲線相位變化逐漸滯后于徑流序列的變化。1982-2003年西江上游的植被覆蓋指數NDVI出現輕微下降的趨勢[20],由此推測當時所實施的水土保持等措施并沒有顯著增加植被的覆蓋率,而是降低了植被覆蓋度的減少趨勢,而輸沙量產生變異可能主要受龍灘樞紐等大型水利設施投入使用的影響。從這種趨勢的發展來看,輸沙序列的豐枯突變時期要遠短于徑流序列,會相較于徑流長期處于偏少狀態。
在28 a 時間尺度上,遷江站水沙序列大約經歷了6 個豐枯的周期轉換變化,但出現了明顯不同,徑流序列平均變化周期為19 a 左右,而輸沙序列為16 a 左右,并且輸沙序列相位變化明顯滯后于徑流序列。徑流在1957-1965年、1975-1985年、1992-2003年和2010-2016年的各時段內處于正位相,表明這幾個時段內徑流為偏豐狀態;1966-1974年、1986-1991年和2004-2009年徑流位相為負,表明處于偏少狀態。輸沙量1962-1971年、1981-1989年和2001-2009年處于偏多狀態,1957-1961年、1972-1980年、1990-2000年和2010-2016年處于偏少狀態。
通過小波分析,可知遷江站的水沙序列主要以5、16 和28 a左右的周期作用為主導,這3 個周期在不同程度上共同決定遷江站水沙序列在整個時域上的周期變化特性。大尺度的未來趨勢可以對小尺度未來進行預測[22],從而為系統中、長期預測提供依據。基于此,遷江站的輸沙序列由于水土保持等相關措施的實施和大型水利設施的投入使用,已滯后于徑流序列的變化,預計長時間內輸沙序列相對于徑流序列將處于偏少狀態;在5a 時間尺度下,徑流將處于偏枯期,之后由于16 a 和22 a 時間尺度的影響會從偏枯逐漸進入偏豐狀態。總體來說,徑流序列相較于輸沙序列振幅穩定,輸沙序列對于環境的變更反應劇烈但逐漸趨于平穩。
根據所處理后的數據,對西江中游的水沙系統進行非線性擬合,結合公式(8)~(12),可以求得以t來表示的VA和VB:

根據上述方程,可計算出1957-2016年間西江中游遷江站徑流系統演變速度VA,輸沙系統演變速度VB,正切值tanα以及兩系統耦合度的α值。徑流序列的R2為0.41,輸沙的為0.71,輸沙序列的顯著較徑流序列高,這是由于20世紀80年代后降雨量有所波動,導致徑流變化更為復雜[6]。為了更直觀地反映出西江中游遷江站徑流系統與輸沙系統的耦合演化過程,依據結果做出了兩系統的耦合曲線,見圖8。

圖8 西江中游徑流與輸沙系統耦合演化曲線
1957-2016年間,徑流與輸沙系統絕大部分處于極限發展時期和再生時期。伴隨著科學技術的進步,人為活動對于徑流和輸沙系統的影響日益提高,徑流系統與輸沙系統的耦合度逼近臨界值(90°),兩個系統達到和諧,但新的矛盾將會出現,預計耦合系統會逐步步入以新矛盾為主體的新初級發展時期,主要集中表現在:
(1)20世紀50年代,大煉鋼鐵運動和農業發展活動導致大規模的植被被砍伐,水土流失嚴重,是造成河道輸沙量增加的主要原因;到1983年后,以植樹造林為核心的水土保持措施開始實施,西江中游的輸沙量開始下降,且徑流與輸沙系統有回歸初級發展時期的趨勢。然而,這并不意味著人為活動引起的矛盾得以解決,徑流輸沙系統回歸至原有的協調發展。
(2)《中華人民共和國水土保持法》實施以來,截止1998年末廣西壯族自治區植樹造林面積達3 591 km2;1992年巖灘水電站開始運行,西江中游的輸沙量明顯減少,大型水利工程對輸沙量下降的影響不可忽視。此外,在1964-2016年間水庫的修建致使紅水河流域的輸沙量減少約80%,植被覆蓋的增長僅僅貢獻輸沙量變化的20%,大型水利工程的出現對于西江中游輸沙量減少起主導作用[23]。在1992年后,科學技術不斷發展,人類活動通過修建大型水利工程逐步影響西江中游的輸沙情況,徑流與輸沙系統逐步邁入再生時期。隨著水利工程的使用年限推后,以及人類活動的影響,西江中游的水沙序列不可避免的進入新的初級發展時期,新的水沙系統矛盾將會出現。
通過Mann-Kendall 趨勢檢驗可以得出遷江站年水沙序列具有下降的趨勢,并結合小波分析,利用耦合系統方程可以得到水沙序列在演化過程中有以下特征。
(1)西江中游年徑流量與年輸沙量整體存在水多沙多、水少沙少的情況;通過Mann-Kendall 趨勢檢驗,年徑流量和年輸沙量均存在顯著下降的趨勢。
(2)遷江站水沙序列在演化過程中存在相似周期變化特征。在25~30 a 的時間尺度上水沙序列周期震蕩十分顯著,乃至波動影響整個時間序列;周期分析結果存在較為明顯的峰值,且峰值出現所對應的時間尺度較為相近,表明作用水沙序列的主周期基本一致。
(3)研究表明,近些年水沙序列所演化出的局部特征,如5 a、16 a和28 a時間尺度周期。這種局部特征產生的主要受紅水河大型水利設施投入影響,輸沙序列的變化相較于徑流序列更為明顯。在5 a 時間尺度下,通過對于水沙序列的周期變換分析得知,徑流序列會在第一主周期28 a 及第二主周期16 a 時間尺度的影響下由偏枯狀態進入偏豐狀態;輸沙序列由于水土保持等相關措施的實施和大型水利設施的作用,將長期處于平穩狀態,其量相較于徑流序列要偏少。
(4)耦合系統中,徑流序列耦合方程的R2為0.43,輸沙序列的為0.71,預計西江中游的水沙系統從極限發展時期進入了再生時期。隨著人類環保意識的增強,加之科學技術的進步,大型水利工程及水土保持措施的應用導致輸沙量產生巨大變化,將之前水沙系統的矛盾得以有效解決,但水沙系統將會面臨新的矛盾,進入新的初級發展時期。在進入新的初級發展時期后,原有的水沙耦合系統可能無法準確反應新的矛盾特征,可能需要進一步加入其他系統進行補充。 □