基于R/S和Morlet小波分析的北洛河上游徑流變化特征

謝智博, 穆興民,2, 高 鵬,, 邱德勛

(1.西北農林科技大學 水土保持研究所 黃土高原土壤侵蝕與旱地農業國家重點實驗室,陜西 楊凌 712100; 2.中國科學院 水利部 水土保持研究所, 陜西 楊凌 712100; 3.中國科學院大學, 北京 100049)

徑流量是流域氣候及下墊面綜合作用的結果，同時是流域水資源分配的重要依據，其變化通常可以使水資源供需關系發生改變從而導致旱澇災害的發生。黃土高原地區氣候暖干化以及大范圍退耕還林(草)等生態修復工程的實施導致黃河干流及其支流的徑流量發生了顯著的變化，進而對流域生態環境和社會經濟可持續發展造成一定的威脅[1-3]。目前，國內眾多學者對于徑流變化已經做了大量研究。如李二輝等[4]通過分析黃河干流陜縣站和河口鎮站徑流量演變過程，得到自1985年以來，黃河上、中游年徑流量呈不斷減小趨勢，徑流變化具有明顯的階段性，經歷了枯—豐—枯3個時期。在渭河流域徑流年際變化方面，黃晨璐等[5]依據渭河1956—2016年實測徑流資料，得到渭河徑流量呈顯著減少趨勢，且分別在1970年和1994年發生突變。對于徑流量減少的原因，張麗梅[6]、Zhao[7]等對渭河流域多年實測徑流數據進行了系統分析，指出劇烈的人類活動所引起的下墊面變化是徑流量減少的主要因素，而且在不同水土保持措施的實施階段，其對徑流變化的影響程度要遠遠大于氣候變化(如降水變化和潛在蒸散發)。針對渭河流域徑流量未來的變化趨勢，王小杰等[8]研究表明，渭河流域徑流量未來變化趨勢與歷史演變規律具有正相關性，即渭河干流徑流將呈持續下降趨勢。北洛河是黃河中下游多砂粗砂區和暴雨洪水等自然災害的主要源頭，其水文水資源的變化同樣引起了學者們的廣泛關注。如晁智龍[9]和董其廣等[10]以實測水文數據為基礎，分析了北洛河流域徑流的變化，得出志丹、吳旗、劉家河以及交口的年徑流量均表現為極顯著減少趨勢，而且在20世紀60，70年代，其徑流量在年內分配比較集中，豐水期和枯水期徑流量變化趨勢與全年基本保持一致。隨著近幾年北洛河上游年徑流量持續偏枯，研究其徑流變化趨勢和分配特征及周期性變化規律則具有一定的現實意義。

本文基于北洛河上游吳旗水文站1971—2014年逐月徑流量數據，采用M-K(Mann-Kendall)趨勢分析法和R/S法(rescaled range analysis)分析汛期和全年徑流量的變化趨勢及其變化連續性，同時采用Morlet小波變換分析徑流量年周期特征，以期為該流域生態環境建設以及可持續發展提供理論參考。

1 研究區概況

北洛河是渭河的一級支流，黃河二級支流。本文中北洛河上游指的是吳旗水文站以上流域(107°32′39″—108°32′45″E，36°44′57″—37°19′30″N)。流域面積為3 408 km2，河長275.6 km，吳旗水文站位于36.9°N，108.2°E，海拔1 331.4 m。地貌類型屬于黃土高原丘陵溝壑區，地勢總體為西北高，東南低。研究區所屬氣候為溫帶大陸性季風氣候，春季干旱多風，溫差較大，夏季悶熱，秋季霜凍期較早，冬季風多雪少。流域內主要土壤類型為黃綿土，年均氣溫7.7℃，無霜期96～145 d，多年平均降雨量454.9 mm，降水量年內分配不均，主要集中在7—9月，且多以暴雨形式出現[11-13]。

2 數據來源與研究方法

2.1 數據來源

根據《黃河流域水文年鑒》中北洛河流域吳旗水文站1971—2014年逐日徑流量資料整理計算得到月徑流量和年徑流量，數據序列經過嚴格檢查，以確保其準確和可靠。

2.2 研究方法

2.2.1 Mann-Kendall趨勢檢驗法 Mann-Kendall趨勢檢驗法(簡稱M-K檢驗)是一種非參數統計檢驗方法，廣泛應用于水文時間序列的分析研究。因為其能夠對水文氣象要素的整體變化趨勢進行定量刻畫，在時間序列分析上具有獨特的優勢，能夠很大程度上減少非正常值對結果的影響，也不會受到數據分布特征的影響，可以對水文氣象序列的趨勢進行有效的分析研究[14-16]，其統計量的計算與分析可參考文獻[17—19]。

2.2.2 R/S分析法 R/S分析法是一種識別時間序列長程相關(反相關)持續性(反持續性)的簡單而又有效的方法，能夠很好地運用于水文時間序列的分析，為研究時間序列長程相關性提供了理論支撐，其最大優勢在于不必假定時間序列的分布特征，能夠大大簡化對系統的長程相關特性的分析[20-21]。根據R/S分析可以得出Hurst指數，其值大小可作為時間序列是否存在趨勢性成分的依據，并表征其趨勢強度[22]。對于一個非隨機水文序列，首先應該滿足[20,23]：

R(h)/S(h)=(ah)H

(1)

式中：R(h)/S(h)為重標極差；h為增量區間長度；a為常數；H為Hurst指數，其過程為：

其次，計算Dm(h個期間)的累積離差Xk,m：

(2)

接著用Xk,m來算Dm極差Rm以及其準差值Sm，根據以上公式結果來計算Dm的重標極差。

(R/S)m=Rm/Sm

(3)

對于每一個子序列都需要重復以上步驟，那么可以得到整個序列的重標極差.

(4)

不斷增大h，直到滿足h=N/3，最終可以得到{x(t)}的重標極差系列，并且滿足公式(1)，對公式(1)兩邊取對數可得：

lg(R/S)h=Hlga+Hlgh

(5)

解釋變量為lgh，被解釋變量為lg(R/S)，以此來作散點圖，用最小二乘來擬合直線，直線斜率就是Hurst指數。

該方法判斷規則：若0

2.2.3 小波分析法 小波分析(wavelet analysis)在分析水文時間序列方面具有強大的分辨多尺度的功能，而且在時域和頻域上具有良好的局部化功能，能識別出水文序列或高或低不同的頻率成分和不同時間尺度的主要變化周期，從而對未來發展趨勢做出準確預測[17,24]。本文所用的Morlet小波定義如下[25]:

(6)

式中：Wf(m,n)為小波變換系數；m為尺度因子，代表小波的周期尺度；n為時間因子，代表時間上的平移。

Wf(m,n)能同時反映參數m和n的特征，并隨著二者變化而變化。用n當橫坐標，m當縱坐標來繪制關于Wf(m,n)的小波變換系數圖，根據圖中二維等值線閉合中心的正負值和小波系數的零點來判斷徑流的增減性和突變性。以不同時間尺度m為橫坐標，以各尺度下的小波方差為縱坐標繪制小波方差圖，徑流序列的主要變化周期可參照圖中各個峰值來確定[26]。

本文使用Matlab軟件的Complex Continuous Wavelet工具箱計算得到北洛河上游吳旗站1971—2014年年徑流量和汛期徑流量的小波系數后，用Surfer 16軟件繪制小波變換實部分布圖和小波系數模等值線圖，用Excel 2019繪制小波方差圖和第一主周期時間尺度小波系數實部過程線圖。小波變換實部分布圖能夠反映徑流量在不同振蕩周期尺度下的周期性變化；小波系數模等值線圖能夠反映不同振蕩周期尺度下周期性的強弱，其模越大，周期性特征就越明顯；小波方差圖反映不同振蕩周期尺度下波動能量的分布情況，峰值用來確定徑流量的變化主周期[27]。

3 結果與分析

3.1 徑流量年內變化

從圖1可以看出，北洛河上游徑流主要集中在6—9月，峰值一般出現在8月。1971—2014年各年代徑流量的年內分配過程呈現出基本一致的特點，為“雙峰型”，峰值分別出現在3月和8月，形成春汛和夏汛，且峰值在20世紀70—80年代以及20世紀90年代—21世紀20年代均呈現降低趨勢。從全年變化趨勢來看，1—3月徑流略微增加，主要原因是冰雪融水和春季降水較多，同時由于溫度較低，蒸發較弱，所以比較容易產生徑流；3—5月徑流略微下降，說明隨著溫度升高，蒸發隨之變大，導致春旱嚴重；6—8月徑流顯著增加，主要原因是夏季降水充沛，雨量大且歷時短，地表起伏不平，容易形成地表徑流。

圖1 吳旗水文站不同年代徑流量年內分配曲線

由表1可以看出，徑流量最大極值比出現在90年代，而最小的極值比則出現在2011—2014年，在整個44 a內，極值比有減小的趨勢，說明北洛河上游年內分配的相對變化幅度有減小的趨勢。北洛河上游不同年代年內分配變差系數為0.48～1.09，變化幅度較大。20世紀90年代變差系數和完全調節系數均大于多年平均值且為最大值，表明該時期內月徑流序列間的差異極大，徑流年內分配狀況十分不均，可能與這一時期極端天氣導致的多場暴雨有關，尤其是發生于1994年8月的百年一遇特大暴雨，吳旗水文站產生的徑流量是其多年平均徑流量的1.3倍。2001—2014年，偏態系數越來越小，說明北洛河上游年內徑流分配逐漸趨于對稱。整體上，北洛河上游徑流年內分配變差系數與完全調節系數以及偏態系數之間的變化規律呈現出一致性：1971—2014年，三者數值大體呈不斷減小趨勢，年內月徑流相對月均值的偏斜程度不斷減小，年內分配不均勻性逐漸變小，說明隨著時間的推移，徑流量年內分配不均勻狀況有所緩解，即北洛河上游徑流量年內分配逐漸趨于均勻化。

表1 吳旗水文站不同年代徑流量年內分配特征值參數

3.2 徑流量年際變化

北洛河上游全年、汛期的年際徑流量序列變化趨勢具有較高的一致性(圖2)，具體表現為兩個徑流量序列有著相似的“峰值”和“谷值”特性，并呈現出逐漸減小的趨勢。其中，全年徑流量序列與汛期徑流量序列峰值出現的時間基本相同，以1975年和1981年為例，汛期和全年徑流量序列均各出現一個峰值。兩個徑流量序列的最大值分別為汛期1.36億m3(1994年)、全年1.62億m3(1994年)，最小值分別為汛期0.20億m3(2001年)、全年0.42億m3(2008年、2012年)。根據徑流量序列線性回歸方程可知：全年和汛期的徑流量減少速率較為接近，全年徑流量平均每10 a減少0.134億m3，汛期徑流量則平均每10 a減少0.115億m3。由于汛期徑流量在年徑流量中占有較大比重，因此兩個徑流量序列減少速率相近。徑流量減少一方面是因為氣候變化，另一方面是因為人類活動的加劇引起流域下墊面發生變化，且后者對徑流量變化的影響程度要遠大于前者[6]。為了進一步反映北洛河上游徑流變化趨勢特性，采用M-K趨勢檢驗法和R/S分析法對兩個徑流量序列變化趨勢的顯著性和持續性進行分析，結果見表2。兩個徑流序列M-K檢驗值和線性擬合結果較為一致，得到的統計量Z值均小于0且數值遠遠超出臨界值-2.58，表明徑流量序列變化趨勢均在0.01置信水平上顯著減少。根據R/S分析得出汛期徑流量和全年徑流量Hurst指數為分別為0.69，0.53，均介于0.5～1，表明兩個徑流量序列均存在一定的長期記憶性，即遞減趨勢在未來都具有一定的持續性，但全年徑流量持續性的可能性要小于汛期徑流量序列。

圖2 吳旗水文站各徑流量序列線性趨勢

表2 吳旗水文站各徑流量序列趨勢特征

3.3 徑流量周期特征分析

根據多年徑流量小波分析結果(圖3)，同一時間尺度下，小波系數正值等值線代表豐水期，負值等值線代表枯水期，零值等值線則代表平水期。在44 a時間序列中，年徑流量存在2～5 a，16～27 a兩種時間尺度的枯—豐變化振蕩周期和8～14 a時間尺度的豐—枯變化振蕩周期，其中8～14 a豐—枯變化震蕩周期和16～27 a枯—豐變化振蕩周期最為顯著。8～14 a周期性變化具有全域性，共經歷了豐—枯—豐—枯—豐—枯—豐—枯—豐—枯—豐—枯12個交替變化，具體豐水時段為1971—1972年、1977—1980年、1984—1987年、1992—1994年、2000—2002年、2006—2009年、2013—2014年，且2014年后實部正值等值線未閉合，說明豐水期在2014年后還會持續一段時間。16～27 a尺度上，全年徑流也存在豐枯振蕩，經歷了枯—豐—枯—豐—枯—豐—枯—豐8個交替變化，但比較紊亂。在2～5 a尺度上，周期性最強(圖4)，在時域上分布也比較明顯，主要發生在1971—2000年，振蕩中心在1994年左右；10～14 a和16～24 a尺度上周期性也比較強，其中10～14 a尺度具有全域性，16～24 a主要發生在1971—1994年，振蕩中心在1971—1973年以及1975年左右；其他時間尺度上周期性則相對較弱。

從圖5可以看出，3個較為明顯的峰值，可以看出年徑流存在3，11，19 a的周期變化特征。其中3 a時間尺度對應最大峰值，所以3 a的時間尺度對應的周期變化最強烈，是其年徑流變化的第一主周期。11 a和19 a時間尺度分別對應第二主周期和第三主周期。將小波方差分析與小波實部時頻變化相結合，可以推測出年徑流3個主周期變化規律相互對應，并且3個主周期均和其時間尺度中心相重合。其中，在3 a時間尺度上，徑流的平均振蕩周期為1.5 a左右，1971—2014年大約經歷了19個豐水期和18個枯水期，可以作為水文周期預報的參考性周期，同時根據徑流的豐枯周期，可以適當調整工業、農業用水，保證城市生活供水情況(圖6)。另外，年徑流序列還可能存在更長時間尺度的變化周期，但這需要更長的徑流序列來進一步驗證。

圖5 年徑流量小波系數方差

圖6 年徑流量3 a時間尺度的小波系數實部過程線

汛期徑流量序列與年徑流存在相似的豐—枯交替變化周期，其所經歷的豐枯變化過程也極為相似，且包含的主周期相同(圖7)，即均包含3，11，19 a共3個變化主周期，推測可能和汛期徑流量在全年徑流量中占有較大比重有關。已有研究表明：西北地區旱澇存在9 a和2～3 a的短準周期振蕩[28]，而且本文得出的結論在有關降水量的周期性研究中與主汛期降水與年降水的周期變化較接近的結論相一致[29]，由此可以說明北洛河流域年徑流量與西北地區氣候的變化周期基本一致。

圖7 汛期徑流量小波系數方差

4 結 論

(1) 北洛河上游年內徑流量主要集中分布在6—9月，占全年徑流的70%以上，最大徑流量一般出現在8月，徑流量年內分配曲線呈“雙峰型”。2000年以前，徑流量年內分配不均勻程度較大，2000年之后，年內分配不均勻程度相對較小且逐漸減弱。

(2) 1971—2014年北洛河上游全年、汛期兩個徑流量序列的變化均呈顯著減少趨勢，其減少過程在未來均存在一定的持續性。

(3) 北洛河上游徑流量存在豐—枯交替的多尺度變化特征和多主周期變化規律。年徑流量和汛期徑流量序列呈現出相似的變化周期，其3個主周期依次為3，11，19 a左右。其中，當時間尺度為3 a時，徑流量的平均振蕩周期為1.5 a左右。