(青海省水利水電勘測設(shè)計研究院，青海 西寧 810001)

徑流和泥沙是造成水土流失的最主要的兩個因素，研究流域水沙特性演變是避免水土流失、防止洪旱災(zāi)害的關(guān)鍵環(huán)節(jié)。近些年，隨著氣候變化和人類活動的加劇，河流水沙的特性也受到影響。高原山區(qū)作為我國大多數(shù)河流的發(fā)源地，河流的演變將會直接影響到下游河流的演變，進而影響人類社會與水資源的關(guān)系。因此，對高原山區(qū)河流的水沙特征演變規(guī)律進行分析具有重要的理論意義和現(xiàn)實價值。本文以格爾木河為例，對高原山區(qū)河流水沙演變特征進行研究。

1 河流概況

格爾木河位于東昆侖山脈與柴達木盆地邊界，是青藏高原環(huán)境演化的關(guān)鍵區(qū)域之一，也是流域內(nèi)最大的河流，發(fā)源于昆侖山北麓，是柴達木盆地的第二大河流[1-2]。格爾木河全長468km，流域面積2.19萬km2，屬于典型的內(nèi)陸河水系。

2 研究方法

2.1 線性回歸方程

當隨機變量y和x的相關(guān)系數(shù)為線性時，稱為兩變量線性相關(guān)，關(guān)系為

y=a+bx+c

(1)

式中：x為自變量；y為因變量；a和b為回歸系數(shù)；c為殘差。

從而兩變量理論回歸方程為

y=a+bx

(2)

采用最小二乘法估計參數(shù)a和b，使觀測點盡量靠近回歸直線，以達到回歸方程擬合效果最好。若假設(shè)離差

(3)

2.2 滑動平均模型法

運用滑動平均模型法可以減小數(shù)據(jù)的隨機性，從而顯現(xiàn)出數(shù)據(jù)的光滑變化趨勢。其數(shù)學(xué)表達式為

(4)

式中：t為經(jīng)過滑動平均后的序列項；x為序列值；m為整體本中的一部分；i為m的相反數(shù)；2m+1為滑動尺度。

2.3 Mann-Kendall秩次相關(guān)檢驗法

Mann-Kendall秩次相關(guān)檢驗法[3]是一種非參數(shù)統(tǒng)計檢驗方法。對于具有N個樣本量的時間序列n，構(gòu)造一個秩序列：

(5)

其中

(6)

式中：秩序列Sk為第i時刻數(shù)值大于j時刻數(shù)值個數(shù)的累計數(shù)。

在時間序列隨機獨立的假定下，定義統(tǒng)計量：

(7)

式中：UFk、E(Sk)和var(Sk)分別為累計數(shù)Sk的均值和方差，在x1,x2,···,xn相互獨立，且有相同連續(xù)的分布情況下：

(8)

(9)

UFk為標準正態(tài)分布，它是按時間順序x1,x2,···,xn計算出的統(tǒng)計量序列，UFk的絕對值在大于等于1.96時分別通過了置信度95%顯著檢驗，若|UFk|>Uα，則表明序列存在明顯的趨勢變化[4]。

按時間序列x逆序，把此方法引用到逆序排列中，再重復(fù)上述過程，并使計算值乘以-1，得出UBk。所有的UFk組成一條曲線，通過置信度檢驗可知其是否具有趨勢。若UFk超出置信區(qū)間，表明有明顯的變化趨勢。若UFk、UBk兩曲線的交點在置信區(qū)間內(nèi)，則該點可能是突變點。

2.4 滑動t檢驗

(10)

其中

(11)

式(10)遵從自由度γ=n1+n2-2的t分布。

2.5 Hurst Exponent檢驗及改進

Hurst Exponent是1951年由英國水文學(xué)家Hurst在研究尼羅河洪水時發(fā)現(xiàn)的數(shù)據(jù)不服從布朗運動及正態(tài)分布的特性，隨后提出了Hurst系數(shù)。通過Hurst系數(shù)可以定量表征時間序列的持續(xù)性或長期相關(guān)性[6]，其Hurst系數(shù)h代表的意義為：當h=0.5時,標志著樣本序列是隨機的，未來的變化趨勢不受現(xiàn)在變化趨勢影響；當h>0.5時，表明序列有正持續(xù)性，未來變化趨勢與現(xiàn)在變化趨勢相同；當h<0.5時，表明序列有反持續(xù)性，未來變化趨勢與現(xiàn)在變化趨勢相反。h越接近于1，表明序列正持續(xù)性越強；反之，h越接近0，表明序列反持續(xù)性越強。

一般可以通過R/S分析計算Hurst指數(shù)值，定義極差與標準差之間的比值為R/S，則

R/S=(cτ)h

(12)

根據(jù)樣本實測資料，可以采用最小二乘法求得c和Hurst系數(shù)b。

雷紅富等[5]在此基礎(chǔ)上，通過計算給定顯著性水平下的Hurst系數(shù)值的相關(guān)函數(shù)，進一步將變異程度細化為5個區(qū)間：無變異、弱變異、中變異、強變異和巨變異5個等級。

3 水沙演變特征分析研究

3.1 年徑流演變特征分析

分別采用線性回歸、滑動平均模型、Mann-Kendall檢驗、滑動t檢驗和Hurst Exponent檢驗等方法對格爾木河流域格爾木站1956—2017年徑流數(shù)據(jù)進行計算分析。

圖1 格爾木站年徑流線性趨勢及滑動平均曲線

由圖1可知，1956—2017年間，格爾木河年徑流整體呈下降趨勢，平均下降為252.49萬m3/a。結(jié)合5年滑動平均的曲線可以看出，1956—1961年徑流量有小幅度的減少；但隨后的1962—1979年間，徑流量表現(xiàn)出連續(xù)18年的上升趨勢；之后，1980—1993年徑流量又逐年減少至1956年以來的最低值；1994—2017年，年徑流開始有緩慢的上升和平穩(wěn)趨勢，但上升幅度較小。總體情況來看，上升經(jīng)歷時間較長，但下降時間集中，且下降幅度較大，年徑流量整體依舊是減少趨勢。

圖2 格爾木站年徑流Mann-Kendall曲線

通過分析圖2可知，UFk值在1956—1959年雖有正負波動，但均不顯著，1960—1975年，UFk值為負數(shù)，表明該段時間內(nèi)年徑流出現(xiàn)減少趨勢,且在1967年、1968年和1969年均超出顯著范圍，表明從1967年開始呈現(xiàn)顯著下降趨勢；直至1976年，UFk出現(xiàn)正值，并在1981年超出顯著范圍，表明從1981年開始年徑流呈現(xiàn)顯著上升趨勢，該趨勢一直持續(xù)到1995年；隨后，1996年開始UFk值再度回落到0值以下，表明年徑流再次呈現(xiàn)下降趨勢；總體而言，年徑流量雖在1976—1995年有所上升，但就總體趨勢而言，年徑流量仍呈現(xiàn)下降趨勢。

從UFk和UBk曲線的交點情況可知，總共出現(xiàn)過6次交點，對應(yīng)的年份分別是1957年、1959年、1970年、1972年、1989年和2009年。它們均可能是突變點出現(xiàn)的年份。

圖3 格爾木站年徑流滑動t檢驗曲線

滑動t檢驗中，采用步長為5年的滑動距離。從圖3可以看出，年徑流在1957年、1989年、1999年、2005年、2008年和2009年均超出了0.05顯著水平區(qū)間臨界范圍，這些點對應(yīng)的年份均可能出現(xiàn)突變現(xiàn)象。

通過將Mann-Kendall檢驗和滑動t檢驗結(jié)果相結(jié)合，兩種方法出現(xiàn)的突變年份中重疊的有1957年、1989年和2009年，兩種方法相互驗證，可以確定這3個年份為年徑流序列的突變年。結(jié)合UFk曲線可知，1957年表現(xiàn)為年徑流下降突變年，1989年為上升突變年，2009年為下降突變年。1956—2017年這52年總體呈下降趨勢。

通過Hurst指數(shù)方法檢驗未來變化趨勢。通過計算可知，年徑流序列的Hurst指數(shù)為0.8094。結(jié)合水文時間序列未來趨勢特征表(表1)可知，Hurst指數(shù)大于0.5時，通過Mann-Kendall檢驗計算，得到U=-0.19，由此表明年徑流序列具有持續(xù)性，年徑流將出現(xiàn)持續(xù)下降的狀態(tài)。

表1 序列未來變化趨勢特征

3.2 年輸沙量演變特征分析

分別采用線性回歸、滑動平均模型、Mann-Kendall檢驗、滑動t檢驗和Hurst Exponent檢驗等方法對格爾木河流域格爾木站1956—2017年輸沙量數(shù)據(jù)進行計算分析。

圖4 格爾木站年輸沙線性趨勢及滑動平均曲線

由圖4可知，1956—2017年間，格爾木河年輸沙量整體呈現(xiàn)出下降的趨勢，平均下降為0.6162萬t/a。結(jié)合5年滑動平均曲線可以看出， 1956—1960年輸沙量有下降趨勢；1961—1981年，輸沙量呈現(xiàn)緩慢上升的趨勢，在1972—1975年期間存在較小的波動；之后在1973—1990年的18年內(nèi)，年輸沙量經(jīng)歷了較大幅度的下降，直至1990年下降到最低點；1990—2000年的10年間，年輸沙量有所回升，但上升幅度不大；隨后的2001—2017年，基本呈現(xiàn)持平的趨勢?？傮w而言，輸沙量總體表現(xiàn)出小幅度的下降趨勢。

圖5 格爾木站年輸沙Mann-Kendall曲線

通過圖5可知，UFk值在1956—1960年雖然有波動，但幅度比較??；在1961—1973年，UFk值均小于0值，在1969年，UFk值出現(xiàn)最小值，達到臨界范圍邊緣，在1961—1973年的13年內(nèi)，輸沙量呈現(xiàn)明顯的下降趨勢；隨后在1974—1996年，UFk值表現(xiàn)為正值，并在1989年達到最大值，表明此階段輸沙量呈現(xiàn)上升趨勢；1997年開始UFk值回落到0點以下，直到2017年UFk值均保持在負值狀態(tài)，表明該階段年輸沙量呈現(xiàn)下降狀態(tài)?？傮w而言，年輸沙量雖然在1974—1995年期間呈上升趨勢，但整體呈下降趨勢。

從UFk和UBk曲線的交點情況可知，總共出現(xiàn)過4次交點，對應(yīng)的年份分別是1958年、1960年、1970年和1989年。它們均可能是突變點出現(xiàn)的年份。

圖6 格爾木站年輸沙量滑動t檢驗曲線

同樣，采用步長為5年的滑動t檢驗進行計算分析。從圖6可以看出，年輸沙量在1958年、1969年、1976年、1977年、1979年、1980年、1989年和1999年均超出了0.05顯著水平區(qū)間臨界范圍，這些點對應(yīng)的年份均可能出現(xiàn)突變現(xiàn)象。

通過將Mann-Kendall檢驗和滑動t檢驗結(jié)果相結(jié)合，兩種方法出現(xiàn)的突變年份中重疊的有1958年和1989年，兩種方法相互驗證，可以確定這2個年份為年輸沙序列的突變年。結(jié)合UFk曲線可知，1958年表現(xiàn)為年輸沙量下降突變年，1989年為上升突變年。在1956—2017年這52年總體呈現(xiàn)下降趨勢。

結(jié)合Hurst指數(shù)方法檢驗未來變化趨勢。計算可知，年輸沙量序列的Hurst指數(shù)為0.9049。根據(jù)未來趨勢特征表(表1)，Hurst指數(shù)大于0.5時，通過Mann-Kendall檢驗計算得到U=-0.23，由此表明年輸沙序列具有持續(xù)性，年輸沙量也將呈現(xiàn)持續(xù)下降的狀態(tài)，這與年徑流得到的下降結(jié)果表現(xiàn)基本一致。

3.3 水沙相關(guān)分析

從趨勢分析可知，格爾木站的年徑流量和年輸沙量均表現(xiàn)出下降的趨勢，通過Hurst指數(shù)可以判斷，在未來一段時間內(nèi)，水沙都將依舊保持下降的趨勢。從突變點出現(xiàn)的情況可知，年徑流量和輸沙量均在1989年出現(xiàn)上升的突變現(xiàn)象，在1989年之前和1989年之后兩個序列都表現(xiàn)出持續(xù)的下降狀態(tài)。故以1989年為突變年份，將序列劃分為突變前和突變后兩部分，分別研究水沙相關(guān)系數(shù)的變化。

圖7 突變前水沙關(guān)系曲線

圖8 突變后水沙關(guān)系曲線

從突變前后的水沙關(guān)系曲線(見圖7、圖8)可知，在突變前水沙相關(guān)系數(shù)為0.486，突變后相關(guān)系數(shù)增大為0.798。說明隨著突變的出現(xiàn)，水沙相關(guān)性逐漸增強。徑流量減少，說明流量也相應(yīng)減少，其攜帶泥沙的能力也將下降，因此輸沙量也出現(xiàn)下降的情況。

4 結(jié) 論

本研究采用5種方法對格爾木站1956—2017年共52年的徑流量和輸沙量演變特征進行了分析。研究表明：格爾木河年徑流量近52年整體表現(xiàn)出減少趨勢，減少率為252.49萬m3/a；相應(yīng)地，格爾木河年輸沙量近52年同樣表現(xiàn)出減少趨勢，減少率為0.6162萬t/a。水沙變化趨勢基本一致，且在未來一段時間內(nèi)減少趨勢將會持續(xù)。此外，根據(jù)兩種方法得到的突變檢測可知，徑流量和輸沙量都在1989年發(fā)生突變，通過對突變前后水沙相關(guān)系數(shù)分析可知，隨著突變的產(chǎn)生，徑流量和輸沙量相關(guān)系數(shù)由0.486增加到0.798。表明在突變后，徑流量和輸沙量之間的相關(guān)性有所增加。