薛嵩嵩， 高 凡， 何 兵

(1.新疆農業(yè)大學 水利與土木工程學院， 新疆 烏魯木齊 830052 ；2.西安理工大學 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室， 陜西 西安 710048)

1 研究背景

河川徑流是水文循環(huán)的基本環(huán)節(jié)，可以直接反映氣候變化和人類活動對水文循環(huán)的影響[1]。受變化環(huán)境影響，徑流序列若發(fā)生變異將導致水文資料一致性破壞，不利于水文分析與模擬及水資源規(guī)劃與利用工作。識別徑流序列變異點是進行徑流序列重構與一致性再建立的關鍵技術環(huán)節(jié)，對于認知變化環(huán)境下水文序列形成、演變及響應機理和應對徑流變化進行水資源適應性管理具有重要意義[2]。傳統(tǒng)的水文數(shù)理統(tǒng)計方法，如Mann-Kendall法，滑動T檢驗，滑動F檢驗，貝葉斯變點分析法，pettitt檢驗法等，多傾向于識別徑流序列的突變點。如祁文燕等[3]采用Mann-Kendall法對湟水干流近60年徑流進行突變特征分析；陳立華等[4]采用有序聚類法和滑動T檢驗法對西江干流近60年徑流資料進行突變分析；王琪森等[5]采用時序累計值相關曲線法識別南四湖流域突變點。黃強等[6]認為，突變點不等同于變異點。由于傳統(tǒng)的徑流序列突變點診斷方法對參考數(shù)據(jù)序列長度要求嚴格和對噪聲有較強敏感性等特點，應用于非線性、非穩(wěn)定性的徑流時間序列變異點診斷尚存在一定的局限性。伴隨非線性理論的發(fā)展和不斷完善，Pincus于1991年提出近似熵(ApEn)法，以其較強抗干擾抗噪能力、所需數(shù)據(jù)序列短、能獲得穩(wěn)定的計算成果等優(yōu)勢得到廣泛應用[7-8]。在近似熵法基礎上，2000年Richman和Moorman提出一種新的時間序列復雜性度量方法——樣本熵法，該方法是近似熵的改進算法，能更好地維持計算結果的一致性[9]。

葉爾羌河流域是塔里木河流域四源一干水系格局的重要組成部分，葉爾羌河徑流序列變化直接關系到塔里木河流域水生態(tài)安全和水資源開發(fā)利用。變化環(huán)境下，葉爾羌河流域氣溫、降水序列演變特征與新疆氣候變化趨勢一致，即呈暖濕化發(fā)展[10]，在此背景下葉爾羌河徑流序列演變響應及變異點診斷研究成果較少。本文分別采用近似熵法和樣本熵法對葉爾羌河源流區(qū)卡群水文站1957-2016年60 a的日徑流序列進行變異診斷，并對比分析兩種方法變異診斷效果，以期為氣候變化和人類活動影響下干旱區(qū)內陸河流水文過程演變以及水資源的適應性管理研究提供科學決策支持。

2 研究區(qū)概況

葉爾羌河流域東依塔克拉瑪干大沙漠，西接布古里、托格拉克沙漠，南以喀喇昆侖山為屏障，北與天山南麓余脈相毗鄰，流域總面積8.57×104km2, 位于74°28′～80°54′E、34°50′～40°31′N[11]。葉爾羌河是喀什地區(qū)第一大河流，也是塔里木河的主要源流之一，發(fā)源于喀喇昆侖山北坡的拉斯開木河，全長1 179 km，多年平均年徑流量65.45×108m3[12]。

圖1 研究區(qū)示意圖

3 數(shù)據(jù)來源與研究方法

3.1 數(shù)據(jù)來源

卡群水文站(海拔1 450 m)位于76°54′E、37°59′N，是葉爾羌河總水量控制站，多年平均徑流量為68.75×108m3，其中，冰雪融水占70%左右[13]。本文采用葉爾羌河源流區(qū)卡群水文站1957-2016年60 a實測日徑流數(shù)據(jù)序列，并收集同期氣溫和降水數(shù)據(jù)序列。數(shù)據(jù)資料均來源于新疆水文年鑒及葉爾羌河流域水文水資源局。

3.2 研究方法

本文采用累積距平法[14]分析研究區(qū)年徑流序列的趨勢性變化特征，采用突變檢驗法——Mann-Kendall法[15]識別徑流序列的突變點。分別運用近似熵法和樣本熵法診斷長系列日徑流序列的變異點并進行方法對比，限于篇幅關系，本文僅介紹近似熵法與樣本熵法，累積距平法、Mann-Kendall檢驗法詳見參考文獻[14-16]。

3.2.1 近似熵方法 近似熵(ApEn)法由Pincus于1991年提出，該方法基于時間序列復雜性角度度量信號中產生新模式概率的大小，即時間序列在模式上的自相似程度[17]。其值越大，表明產生新模式的概率越大，即時間序列的復雜程度越高[18]，復雜度的高低決定了序列是否存在變異情況，即最值點極可能成為序列的變異點。計算步驟如下：

(1)對時間序列{u(i),i=1,2,…,N}進行空間重構，重構維數(shù)為m，由此可以構造出一組維數(shù)為m的新向量X(1),X(2),…,X(N-m+1)，其中X(i)={u(i),u(i+1),…,u(i+m-1)},i=1,2,…,N-m+1；

(2)計算任意向量X(i)與其余向量X(j)之間的相對歐式距離d[X(i),X(j)]；

(5)維數(shù)m增加 1，重復步驟(1)～(4)，得Cim+1(i)和φm+1(r);

(6)近似熵的估計值定義為：

ApEn(m,r)=φm(r)-φm+1(r)

(1)

式中：ApEn的取值與m和r有關，Pincus建議取m=2，r=(0.1～0.25)σ，其中σ為原始數(shù)據(jù)的標準差[19]。

3.2.2 樣本熵法 樣本熵(SpEn)法由Richman和Moorman于2000年在近似熵法基礎上提出，該方法克服了近似熵在計算時的偏差，近年來廣泛應用于機械故障診斷和醫(yī)學信號處理等領域[20]。該方法計算步驟與近似熵法類似，是對近似熵法的改進，詳細計算步驟可參考文獻[21]。

4 結果與分析

4.1 葉爾羌河源流區(qū)1957-2016年徑流序列變化特征

1957-2016年卡群水文站年徑流量變化過程線見圖2(a)。從圖2(a)可以看出，卡群水文站多年平均年徑流量為67.08×108m3，年徑流量序列呈鋸齒狀波動增加趨勢，Mann-Kendall法檢驗統(tǒng)計量Z值為1.062，年徑流量序列增加趨勢不顯著。由1957-2016年卡群水文站年徑流量累積距平變化過程線(見圖2(b))可以看出，1957-20160年年徑流序列呈現(xiàn)出先下降后上升的“V”型階段性變化特征。其中，1957-1963、2002-2016年呈顯著上升，表明該階段年徑流量豐水期多于枯水期；1964-1976、1990-2001年呈顯著下降，表明該時段年徑流量豐水期少于枯水期；1967-1989年徑流序列呈頻繁波動，表明該時段年徑流量處豐枯交替階段。

圖2 1957-2016年卡群水文站年徑流量變化過程線及年徑流累積距平變化過程線

4.2 Mann-Kendall法突變點診斷

運用Mann-Kendall法對卡群水文站1957-2016年徑流量年序列進行突變點檢驗，結果見圖3。

圖3 葉爾羌河卡群水文站年徑流Mann-Kendall突變檢驗過程線

由圖3可以看出，在檢驗線之間(±1.96)UF和UB曲線共產生9個交點，分別為1958、1959、1960、1982、1987、1990、1994、1996、2000年，即該序列存在9個突變點。由于1958、1959、1960年這3個年份距離該系列的起始時間較近，靠近邊界的交點尚需更長的水文資料進行驗證，且短時間內不會出現(xiàn)頻繁突變現(xiàn)象，故將該3個年份視為一個突變點，即1959年。此外，UF和UB曲線在20世紀80年代至90年代末產生多次交點，表明該時段內年徑流呈波動趨勢，豐枯突變較頻繁，突變點較多，與實際情況不相符，且均未超過95%顯著水平，突變不顯著。通過以上分析表明，Mann-Kendall法診斷得出的突變點無法擁有足夠論據(jù)判定變異點。

4.3 近似熵與樣本熵法變異點識別與對比

4.3.1 近似熵法變異點診斷 1957-2016年卡群水文站日徑流量序列近似熵值(ApEn)變化過程線見圖4。依據(jù)近似熵計算步驟及文獻[22]，本文取m=2，r=0.15σ，(σ為原始數(shù)據(jù)的標準差)。由圖4可以看出，研究區(qū)日徑流序列近似熵值變化過程可以分為4個變化階段，分別為1957-1959、1960-1987、1988-2004、2005-2016年。其中，1957-1959年近似熵值呈持續(xù)增加；1960-1987年間近似熵值呈波動遞增變化，1987年近似熵值達極大值，其徑流復雜程度極高；1988-2004年近似熵值呈波動遞減變化，2004年近似熵值達極小值，其徑流復雜程度最低；2005-2016年近似熵值均高于2004年，可預測性降低。結合曲線趨勢及復雜度可以判斷，1959、1987、2004年卡群水文站徑流序列發(fā)生了變異。

圖4 1957-2016年葉爾羌河卡群水文站年徑流量序列近似熵變化過程

4.3.2 樣本熵法變異點識別 1957-2016年葉爾羌河卡群水文站日徑流序列樣本熵(SpEn)變化過程線見圖5。依據(jù)樣本熵計算步驟及文獻[23]，本文取m=2，r=0.15σ(其中σ為原始數(shù)據(jù)的標準差)。由圖5可以看出，研究區(qū)日徑流序列樣本熵值變化過程可分為3個階段，分別為1957-1987、1988-2004、2005-2016年，整體上呈先上升后下降再上升的階段性變化特征。1987年樣本熵出現(xiàn)極大值，其徑流復雜度最高；2004年樣本熵出現(xiàn)極小值，其徑流復雜度最低，可預測性最高。結合曲線趨勢及復雜度可以判斷，1987、2004年卡群水文站徑流序列發(fā)生了變異。

4.3.3 近似熵法與樣本熵法變異點識別效果對比 數(shù)據(jù)量越大，近似熵值越小[24]，當數(shù)據(jù)量小于一定程度時，近似熵不能反映序列真實的波動狀態(tài)；同理，樣本熵也存在相同情況。因此，兩種方法在維數(shù)一定，相同數(shù)據(jù)量的情況下，統(tǒng)計值越大，表明波動越接近真實狀態(tài)[25]。近似熵與樣本熵法診斷結果對比曲線見圖6，由圖6可以看出，樣本熵值大于同期近似熵值，且樣本熵值變化幅度較近似熵值明顯，表明樣本熵具有靈活的敏感度，可更好反映徑流量時間序列的波動情況，即更加適應徑流量時間序列的變異點診斷。

圖5 葉爾羌河卡群水文站徑流量序列樣本熵變化過程

圖6 近似熵法與樣本熵法診斷結果數(shù)據(jù)對比曲線圖

5 討 論

為了進一步驗證樣本熵法識別研究區(qū)徑流序列變異點的可靠程度，對圖5中徑流量樣本熵變化過程3個階段內的均值、方差及變差系數(shù)進行比較分析，計算結果見表1。

表1 卡群水文站1957-2016年徑流量序列分階段均值、方差、變差系數(shù)參數(shù)

由表1可以看出：(1)1957-1987、1988-2004、2005-2016年3個階段的年徑流量均值呈先減少后增加趨勢，1988-2004年均值最小，1988-2004年徑流量均值相對于1957-1987年減少了14.88%，而2005-2016年徑流均值相對于1988-2004年增加了20.39%；(2)1957-1987、1988-2004、2005-2016年3階段的方差差異性較大，1988-2004年間的年徑流量序列方差最小，說明該階段離散程度越小，且1988-2004年徑流量序列的方差值較1957-1987年減少了103.71%，2005-2016年徑流量方差相對于1987-2004增加了122.89%；(3)1988-2004年間變差系數(shù)最小，相對于1957-1987減少了80%，2005-2016年徑流變差系數(shù)較1988-2004增加了70%。從整體來看，無論是均值、方差，還是變差系數(shù)，均具有相似的變化規(guī)律，即呈先減少后增加的趨勢，進一步驗證樣本熵識別的變異點可靠程度。

葉爾羌河屬于融雪與雨水混合補給型河流，夏季融雪徑流大，冬季徑流量穩(wěn)定[26]。徑流量變異是不同驅動因素(人為因素、自然因素)共同作用的結果，由于卡群站以上為高山分布區(qū)，人為干擾較小，自然因素成為影響研究區(qū)徑流量變化的主要因素[27]。已有相關文獻顯示[28-29]，1986-1987年葉爾羌河流域氣候轉型顯著，2001-2005年暖濕化趨勢進一步加強。氣溫升高，水循環(huán)加快，使得上游冰川物質平衡處于波動狀態(tài)，造成冰川融雪徑流量增加及冰川阻塞潰壩型突發(fā)洪水。本文通過對1960-2016年葉爾羌河卡群站徑流量、年平均氣溫和年降水量相關性分析以及數(shù)據(jù)統(tǒng)計(圖7、8和表2)，所得結論與高鑫等[29]結果一致，即氣溫和降水量呈增長趨勢，氣溫升高對徑流量的影響遠大于降水量的影響。

圖7 1960-2016年葉爾羌河平均氣溫變化曲線

圖8 1960-2016年降水量變化曲線

表2 1960-2016年卡群站徑流量與降水量、氣溫因素相關系數(shù)

注：“*”表示通過0.05顯著水平檢驗。

6 結 論

(1)采用傳統(tǒng)的突變檢驗方法Mann-Kendall法對葉爾羌河出山口卡群水文站1957-2016年年徑流量序列進行突變檢驗，發(fā)現(xiàn)出現(xiàn)多個虛假變異點，無法準確斷定突變年份。

(2)采用近似熵法對徑流量序列進行變異點診斷，發(fā)現(xiàn)1959、1987和2004年3個變異點，采用樣本熵法對徑流序列進行變異點診斷，發(fā)現(xiàn)1987和2004年兩個變異點。通過兩種方法的對比分析發(fā)現(xiàn)，在同一維度下，樣本熵法較近似熵法能更好地反映徑流量序列的波動情況，并且樣本熵法具有更好的靈敏度。

(3)卡群站以上為高山分布區(qū)，人為干擾較小，自然因素成為影響研究區(qū)徑流變化的主要因素。由于氣候轉型，氣溫升高，水循環(huán)加快，上游冰川物質平衡處于波動狀態(tài)，造成冰川融雪徑流增加及冰川阻塞潰壩型突發(fā)洪水。相關性分析結果表明，氣溫和降水量呈增長趨勢，氣溫升高對徑流量的影響遠大于降水量增加的影響。