徑流序列變異診斷的近似熵與樣本熵法對(duì)比研究
——以新疆葉爾羌河為例

薛嵩嵩， 高 凡， 何 兵

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

1 研究背景

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

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

2 研究區(qū)概況

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

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

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

3.1 數(shù)據(jù)來(lái)源

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

3.2 研究方法

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

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

(1)對(duì)時(shí)間序列{u(i),i=1,2,…,N}進(jìn)行空間重構(gòu)，重構(gòu)維數(shù)為m，由此可以構(gòu)造出一組維數(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)計(jì)算任意向量X(i)與其余向量X(j)之間的相對(duì)歐式距離d[X(i),X(j)]；

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

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

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

(1)

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

3.2.2 樣本熵法 樣本熵(SpEn)法由Richman和Moorman于2000年在近似熵法基礎(chǔ)上提出，該方法克服了近似熵在計(jì)算時(shí)的偏差，近年來(lái)廣泛應(yīng)用于機(jī)械故障診斷和醫(yī)學(xué)信號(hào)處理等領(lǐng)域[20]。該方法計(jì)算步驟與近似熵法類(lèi)似，是對(duì)近似熵法的改進(jìn)，詳細(xì)計(jì)算步驟可參考文獻(xiàn)[21]。

4 結(jié)果與分析

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

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

圖2 1957-2016年卡群水文站年徑流量變化過(guò)程線(xiàn)及年徑流累積距平變化過(guò)程線(xiàn)

4.2 Mann-Kendall法突變點(diǎn)診斷

運(yùn)用Mann-Kendall法對(duì)卡群水文站1957-2016年徑流量年序列進(jìn)行突變點(diǎn)檢驗(yàn)，結(jié)果見(jiàn)圖3。

圖3 葉爾羌河卡群水文站年徑流Mann-Kendall突變檢驗(yàn)過(guò)程線(xiàn)

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

4.3 近似熵與樣本熵法變異點(diǎn)識(shí)別與對(duì)比

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

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

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

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

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

圖6 近似熵法與樣本熵法診斷結(jié)果數(shù)據(jù)對(duì)比曲線(xiàn)圖

5 討 論

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

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

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

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

圖7 1960-2016年葉爾羌河平均氣溫變化曲線(xiàn)

圖8 1960-2016年降水量變化曲線(xiàn)

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

注：“*”表示通過(guò)0.05顯著水平檢驗(yàn)。

6 結(jié) 論

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

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

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