基于啟發式分割和近似熵法的徑流序列變異診斷*

黃生志，黃 強，王義民，陳昱潼

(西安理工大學西北旱區生態水利工程國家重點實驗室培育基地，陜西 西安 710048)

隨著全球氣候的變化以及日益加劇的人類活動的影響，流域水文要素以及水文情勢諸如河川徑流、湖泊蓄水以及流域降雨等無論是在時間上還是空間上均有可能發生變異，致使水文資料的一致性遭到破壞[1-2]。水文資料一致性的破壞將會給水文分析、水文模擬、水文計算以及流域的水資源規劃等帶來巨大的影響[3-5]，進而導致變異前的成果不能再繼續使用，否則將產生極大的偏差。對流域水文變異進行診斷，不僅有助于掌握流域水文要素的演變規律，還能避免因使用一致性被破壞的水文資料進行設計或分析而產生重大失誤或者損失。

基于此，國內外不少學者對水文變異進行了研究并取得了大量的研究成果。諸如Lee等[6]運用了貝葉斯推斷理論對水文變點進行了研究。然而，該方法對資料的要求過高且對樣本分布的假設過多，因此在實際的應用過程中操作性不強[7]；陳廣才等[8]針對潮白河水資源分區的45 a年徑流量序列，采用基于平穩與線性的檢驗假定的滑動F檢驗法、滑動T檢驗法以及秩和檢驗法等統計的方法，進行了徑流序列的變異診斷。事實上，水文序列受到氣候變化以及人類活動的影響，往往呈現出高度的非線性和非平穩性，使用傳統的統計診斷方法對水文序列的變異點進行檢驗則往往產生一定的偏差[9]。

由Bernaola-Galván[10]提出的啟發式分割法主要用于處理非線性及非平穩時間序列的均值變異點。同傳統的變異檢驗方法相比，啟發式分割法能夠將一個非平穩的時間序列分割成多個平穩的子序列，各子序列的均值互不相同并表征不同的物理背景，能較好地彌補傳統檢驗方法的不足。由于該方法在分割時采取了一分為二的迭代算法，因此，極大地減少了計算量，具有較好的實用性[9]。此外，由Pincus于20世紀90年代提出的近似熵(approximate entropy)，因其具有較好的抗干擾與處理非穩態序列能力且所需數據量較少等優點，故在機械設備故障的診斷及醫學等領域得到廣泛地應用[11]。鑒于此，本文基于啟發式分割算法與近似熵理論，對渭河流域的年徑流量進行變異診斷，并對其變異進行歸因分析。

1 研究區域與數據

1.1 研究區域概況

渭河作為黃河流域最大的一級支流，發源于甘肅省內渭源縣的烏鼠山，自西向東分別流經甘肅的隴西、天水，陜西的寶雞、西安、渭南，最終于潼關匯入黃河[12]。渭河的干流全長約為818 km, 總流域面積約為13.5萬km2(圖1)。渭河流域位于34°-38°N，104°-110°E，屬于大陸性季風氣候，春季溫暖少雨，夏季雨熱同期且有伏旱，秋季涼爽濕潤，冬季寒冷少雨。流域的最冷月平均氣溫大約是-3～-1℃，最熱月平均氣溫大約為23～26℃[12]。流域降雨空間分布不均，由東南向西北遞減，秦嶺南麓雨水充沛，最大年降水量超過1 000 mm，而平原河套地區年降水量僅為500 mm左右[12]。

圖1 渭河流域及相關水文站位置Fig.1 The location of the Wei River Basin and the related hydrological stations

流域降水年內分配不均且年際變化大，汛期降水量約為年降水量的65%。自從國家實行西部大開發這一國策以及成立了關中-天水國家級經濟開發區后，渭河流域的經濟社會高速發展，工業用水、農業用水量劇增。因此，人類大量開采與利用地下水及地表水，致使渭河流域的河川徑流顯著減少。此外受全球氣候變化的影響，近年來渭河流域的降雨呈現下降的趨勢，如此一來渭河流域水文情勢將發生巨大的變化，日益減少的水資源將難以滿足渭河流域經濟社會的可持續發展。鑒于渭河流域用水安全的重要性，為了進一步掌握渭河流域河川徑流的演變規律，本文對該流域的年徑流量進行變異診斷，并對變異結果進行歸因分析。本文將研究區域劃分為三個部分，分別為林家村水文站以上流域，該區域代表渭河流域的中上游；張家山水文站以上流域，該區域代表渭河最大的支流涇河流域；華縣水文站以上流域，該區域代表整個渭河流域。

1.2 數據

本研究所用數據為渭河流域21個氣象站的1960-2005年的日氣象數據，通過求和運算得到各個氣象站的降雨、氣溫、風速、日照變化率等年氣象序列以及林家村、張家山以及華縣水文站的1960-2005年的日徑流序列，通過求和而得到各水文站的年徑流量，渭河流域及各站點所在的位置如圖1所示。以上氣象數據均來源于中國國家氣象中心，水文數據摘自水文手冊。所有數據在使用之前均進行了三性審查，審查結果表明本文所采用的水文氣象數據具有較好的可靠性、一致性以及代表性，可以用作相應分析。潛在蒸發量是通過Penman-Monteith 公式計算而成。

2 研究方法

2.1 啟發式分割算法

某一時間序列x(t)由N個點組成，其中一個分割點i從序列的左邊沿著該序列依次向右邊滑動，分割點左邊和右邊部分的平均值分別為μ1(i)和μ2(i)，其標準差分別是s1(i)和s2(i)，i點的合并偏差SD(i)可表示為[9,13-14]：

(1)

式中，N1、N2分別表示i點左邊與右邊部分的點數。其中，i點左右兩邊兩子序列均值的差異能夠使用t檢驗統計值T(i)來進行度量[9]：

(2)

T值越大，則說明該點左右兩邊兩子序列的差異越明顯。計算T(t)中最大值Tmax所對應的統計顯著性P(Tmax)，其計算公式如下[9]：

(3)

通過蒙特卡洛模擬可得：η=4.19lnN-11.54且δ=0.40。其中，N表示該序列的長度，v=N-2,Ix(a,b) 是不完全β函數。預先設定一個臨界值P0，當P(Tmax)≥P0，則在該點處將此序列分割成左右兩個均值差異較大的子序列，否則不進行分割[9]。對得到的新序列不斷進行迭代并重復以上操作，直到子序列的長度小于0(0是最小分割尺度)時便停止對其分割。由此，便將原序列分割成幾個均值不同的子序列，而分割點則為該序列的均值變異點。

通常情況下，P0取值范圍為0.5～0.95，而0的取值則不應小于25[15]。通過不斷調整的大小，能夠實現對該時間序列不同尺度上的變異檢測[9]。

2.2 近似熵方法

近似熵方法實質上是一種基于邊緣概率分布統計的定量描述某一時間序列的復雜程度的方法，該算法的具體步驟簡要介紹如下[11]。

假定某一長度為N的時間序列u(1),u(2), …,u(N)，在該時間序列中構造一組m維數的向量X(1),X(2), …,X(N+m-1),其中向量為[11]：

X(i)={u(i),u(i+1),…,u(i+m-1)}

i=1,2,...,N-m+1

(4)

定義兩個向量X(i)與X(j)之間的距離d[X(i),X(j)]是兩個向量中對應的元素的最大值[11]：

d[X(i),X(j)]=max[|u(i+k)-u(j+k)|]

k=0,1,...,m-1

(5)

對于每一個i，定義兩個向量X(i)與X(j)之間的關聯程度：

(j)]≤r的個數}/

(N-m+1)

(6)

(7)

(8)

然而，實際上N不可能為無窮大。當N取有限值時，便可得到近似熵的估計值，表示為：

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

(9)

3 結果與討論

3.1 基于啟發式分割算法的年徑流變異結果分析

根據上述2.1的具體步驟，將P0設定為0.95，0的取值則設為25，分別得到林家村以上流域、張家山以上流域以及整個渭河流域的年徑流變異分析結果，如圖2所示。

由圖2(a)中的第一次分割可知，1971年對應的T值最大，并且其相應的P(Tmax)等于1要大于臨界值0.95。因此，1971年是林家村以上流域年徑流序列的第一個突變點；從第二次分割可知，1994年所對應的T值最大，且其相應的P(Tmax)等于0.998 5，也大于臨界值0.95，故1994年為該流域年徑流序列的第二個突變點。

由圖2(b)可得，張家山以上流域的年徑流序列并無突變點，因其最大T值的P(Tmax)小于0.95。對于整個渭河流域的年徑流序列，由圖2(c)中第一次分割可知，1969年的T值為最大值，其對應的P(Tmax)等于0.992 9，大于臨界值0.95，因此1969年為第一個變異點；由其第二次分割可得，1993年的T值為最大值，其對應的P(Tmax)等于0.971 8，大于臨界值0.95，故1993年為第2個變異。

3.2 基于近似熵的年徑流變異分析結果

根據上述2.2求近似熵的具體步驟，采用滑動窗口的方法，將窗口寬度設為5，滑動步長設為1，m和r分別設為2和0.15σ。其中，σ為原始序列的標準差，分別得到林家村以上流域、張家山以上流域以及整個渭河流域的基于近似熵的年徑流變異分析結果于圖3(a)-(c)中。

這時零售商的和供應商的期望收益均是供應鏈整體期望收益函數的仿射函數，同時各企業的最優定貨量q*相同，這表明在無突發事件下基準契約能實現供應鏈協調.

由圖3(a)可知，林家村以上流域年徑流1960-1971年間的近似熵較大，而1971-1994年間的近似熵較小，1971年明顯是一個變異點，而1994-2005年間的近似熵較大，易知1994年為另一變異點。由圖3(b)可知，張家山以上流域年徑流序列的近似熵沒有顯著的變異情況。因此，該流域的年徑流序列無變異點。而由圖3(c)中整個渭河流域年徑流的近似熵明顯可知，渭河流域年徑流序列明顯存在兩個變異點，分別為1969與1993年。該結果與基于啟發式分割算法的結果一致，進一步論證了以上3個流域年徑流序列變異點的正確性，也充分展示了基于啟發式分割算法及近似熵理論在水文變異診斷上的優越性。

3.3 流域年徑流變異點的印證

3.3.1 累積年徑流與年份關系圖 雖然通過啟發式分割算法及近似熵方法已經發現林家村以上流域年徑流序列存在2個變異點(1971和1994年)，張家山以上流域年徑流序列無變異點存在，整個渭河流域年徑流序列也存在兩個變異點(1969和1993年)。

圖2 年徑流啟發式分割檢驗圖Fig.2 The test chart of annual runoff based on heuristic segmentation algorithm

為了進一步證明這兩個算法的在變異診斷上的準確性，繪制3個流域的累積年徑流量與年份的關系圖，其結果如圖4所示。

由圖4實心圓部分可以明顯看出林家村以上流域累積年徑流被2個變異點分成了3部分，且每部分的線性擬合度都很高。其中，1960-1971與1972-1994年2個部分的線性擬合相關系數高達0.99，而1995-2005的線性擬合系數為0.96。而這兩個變異點恰好與上述3.2節中基于啟發式分割法及近似熵理論所找到的關于林家村以上流域年徑流序列的變異點相一致。圖4空心五角星部分所表示的是張家山以上流域累積年徑流量與年份之間的關系，在整個研究階段1960-2005年間，所有的點均擬合于一條直線，且擬合精度高，其相關系數高達0.99。因此，在1960-2005年間，張家山以上流域年徑流序列無明顯變異點存在。由圖4的空心圓部分可知，整個渭河流域的累積年徑流量被1969與1993年這2個變異點分成了3個部分，每個部分點的線性擬合程度都很高，3個部分的相關系數均高達0.99。因此，通過累積年徑流量與年份的關系圖可知，整個渭河流域的年徑流序列存在1969和1993年這2個變異點，而這又與上述利用啟發式分割法及近似熵理論所發現的整個渭河流域年徑流量的變異的十分吻合。由此進一步論證了基于啟發式分割法及近似熵理論在水文變異診斷中的準確性；與此同時，進一步印證了渭河流域年徑流序列存在變異的事實。其中，林家村以上流域與整個渭河流域的年徑流序列的確存在兩個變異點，而張家山以上流域年徑流序列并無變異點存在。此外，對于林家村以上流域和整個渭河流域被兩個變異點分割成的3個部分而言，其擬合直線的傾斜度均隨著時間推移而不斷減小，由此可知這兩個流域的徑流量有著明顯的下降趨勢，而實際上也確實如此。該研究結果與拜有存[16]的基于差異信息理論的渭河徑流序列變異診斷結果基本一致，他發現渭河流域的年徑流存在這兩個變異點。

3.3.2 變異前后各參數對比 分別計算存在變異點的林家村以上流域以及整個渭河流域變異前后徑流序列的均值、方差以及變差系數等參數，進一步對比分析變異點前后序列的差異性，具體的參數值如表1所示。由表1可知：① 林家村以上流域由兩個變異點所分割成的3部分徑流序列無論是均值還是方差均有較大差別，就均值而言，1972-1994年間的年徑流均值較1960-1971年減少了42.2%，而年徑流均值1995-2005年比1960-1972年減少了77.3%，可以說3個部分的均值差別比較明顯。② 而就方差而言，1972-1994年間的年徑流序列的方差較1960-1971年減少了71.4%，而1995-2005年的年徑流序列的方差同1960-1972年相比減少了86.9%，由此可知，這3部分的方差差異性較大。但是，對于變差系數而言，1972-1994年間的變差系數值較1960-1972年變化不大，而1995-2005的年徑流序列的變差系數與1960-1972年相比減少了40.9%。③ 整個渭河流域的年徑流序列變異前后參數的變化規律與林家村以上流域類似，3個部分的均值與方差的差異性較大，1970-1994年間的變差系數值同1960-1970相比變化不大，而1994-2005的年徑流序列的變差系數較1960-1969年減少了25%。

圖3 年徑流的近似熵檢驗圖Fig.3 The test chart of annual runoff based on approximate entropy method

圖4 三個流域累積年徑流量與年份的關系圖Fig.4 The relationship between accumulative runoff and year in the three basin

總的來說，被兩變異點分割成的3個年徑流序列在均值與方差上差異明顯，并且均值、方差及變差系數有明顯的減小趨勢。由此可側面證明變異點的正確性。

表1 兩流域年徑流序列變異前后參數對比表Table 1 The comparison of parameters of the annual runoff before and after the change points in the two basin

3.4 徑流變異的歸因分析

3.4.1 氣候變化的影響 采用Mann-Kendall法求算變點前后降雨和潛在蒸發的變化趨勢，其結果如表2所示。

由表2可知，林家村以上流域與渭河流域的降雨和潛在蒸發變點前后的變化趨勢相似。降雨先減小后變大，但在整個時間序列表現為不顯著地下降趨勢。潛在蒸發量先增大后減小再增大，在整個時間序列表現為不顯著的上升趨勢。從常規分析上可知渭河流域與徑流聯系密切的降雨與潛在蒸發等氣象因子均未發生顯著變化，人類活動對徑流變異的影響或許更強。然而，渭河流域地處中國大陸性季風區，受西太平洋副熱帶高壓系統的影響。1971年時，西太平洋副熱帶高壓系統有所減弱并發生一個El Nio事件致使中國大部分地區的大氣環流及氣候發生變異，從而導致渭河流域當年降雨量劇減[17]，而此時恰好是林家村以上流域年徑流序列的第一個變異點所對應的時間。而1994年時，另外一個ENSO事件發生，此時西太平洋副熱帶高壓系統異常強大，導致渭河流域當年的降雨量較年平均降雨量減少了14.9%[18]，而此時恰好對應于林家村以上流域及整個渭河流域年徑流序列的第2個變異點。

3.4.2 人類活動的影響 20世紀70年代初，渭河流域修建了不少水利工程，如1970年修建的羊毛灣水庫和馮村水庫，1972年修建的大峪水庫和石門水庫。這些水庫的修建，極大地改變河川徑流的時空分布規律，破壞天然的徑流規律，從而導致徑流序列變異。除了興建水利工程外，流域下墊面的變化也改變河川徑流的演變規律。50年代末60年代初，由于經濟困難與政治運動的緣故，大量的樹木遭到破壞，尤其在1958年的“大煉鋼鐵”時期，集中砍伐了大量的林木，從而造成了嚴重的水土流失，由此也使得1960-1970年間的徑流量較大[19]。但是，隨著水土保持措施的開展，尤其是70-80年代，渭河流域對水土流失開展了一系列有組織、有計劃、有步驟的防治工作，并取得了很好的成效，從而在一定程度上減小了渭河流域的水土流失量。因此，在1973-1994年間河川徑流較1960-1970年有一定的減少。此外，隨著人口的迅速增加以及經濟的高速發展，渭河流域的工業用水和農業用水增長迅速。其中，1990-2000年間，國民經濟各行業用水增加了17.1億m3，年均增長率為3.2%，遠高于全國平均1.1%和黃淮海流域的0.3%。90年代后的渭河流域國民經濟平均用水量為42.62億m3，比90年代前增加了52.6%。大量地開采與開發地表水和地下水，是渭河流域徑流序列變異的又一重要因素。

表2 變點前后降雨和潛在蒸發的變化趨勢Table 2 The trend of precipitation and potential evaporation before and after change points

4 結 論

為了探究渭河流域河川徑流變異情況，本文采用1960-2005年林家村、張家山以及華縣3個水文站的日徑流數據，將研究區域劃分為林家村以上流域、張家山以上流域以及整個渭河流域3個部分。基于啟發式分割算法與近似熵方法對各流域的年徑流序列進行了變異診斷，并利用累積年徑流量與年份的關系圖對變異結果進行判定分析，得到主要結論如下：

1)啟發式分割算法與近似熵方法均具有較好的檢驗均值突變點的能力，具有較強的抗噪抗干擾的能力，十分適應于非線性非穩態的徑流序列的變異診斷，能夠準確地找到其變異點。

2)林家村以上流域存在2個變異點，分別是1971和1994年；張家山以上流域無變異點存在；整個渭河流域存在兩個變異點，分別是1969和1993年。

3)對于林家村以上流域和整個渭河流域而言，被兩個變異點分割成3個子序列的線性趨勢線的斜率隨著年份的增大而不斷減小，由此說明在這兩個流域，其徑流量有明顯的減小趨勢。

4)造成渭河流域河川徑流變異的主要原因是氣候變化與日益增長的人類活動。其中，氣候變化主要是1971年的El Nio事件和1994年的ENSO事件的發生致使當時渭河流域降雨量劇減，從而引起徑流量的顯著減少。人類活動則具體包括水利工程建設、流域下墊面的改變以及由于人口和經濟的增長而大量開采與開發地下水和地表水，由此而使得河川徑流發生變異。

[1] 鄧建偉，宋松柏，盧書超. 石羊河流域年徑流序列的變異診斷[J]. 西北農林科技大學學報:自然科學版，2006，34(4) ：121-124.

[2] 李艷，陳曉宏，張鵬飛. 北江流域徑流序列年內分配特征及其趨勢分析[J]. 中山大學學報:自然科學版， 2007， 46(5) ：113-116.

[3] 丁晶，劉權授. 隨機水文學[M]. 北京：中國水利水電出版社，1997.

[4] 王孝禮，胡寶清，夏軍. 水文序列趨勢與變異點的R/S分析法[J]. 武漢大學學報：工學版，2002，35(2)：10-12.

[5] 詹道江，葉守澤. 工程水文學[M]. 北京：中國水利水電出版社，2000.

[6] AUSTIN F S L, SYLVA M H. A shift of the mean level in a sequence of independent normal random variable-A Bayesian Approach [J]. Technometrics, 1997, 19(4): 503-506.

[7] 肖宜，夏軍，申明亮，等. 差異信息理論在水文時間序列變異點診斷中的應用[J]. 中國農村水利水電，2001(11)：28-30.

[8] 陳廣才，謝平. 水文變異的滑動F識別與檢驗方法[J]. 水文，2006，26(2)：57-60.

[9] 陳廣才，謝平. 基于啟發式分割算法的水文變異分析研究[J]. 中山大學學報:自然科學版，2008，47(5)：122-125.

[10] BERNAOLA-GALVN P, IVANOV P C, AMARAL L A N,et al. Scale invariance in the nonstationarity of human heart rate[J].Physical Review Letters, 2001, 87:168105.

[11] 王啟光，張增平. 近似熵檢測氣候突變的研究[J]. 物理學報，2008，57(3)：1976-1982.

[12] 張宏利，陳豫，任廣鑫，等. 近50年來渭河流域降水變化特征分析[ J].干旱地區農業研究，2008，26(4)：236- 241.

[13] 李海彬，張小峰，胡春宏,等. 基于B-G分割算法的河川年輸沙量突變分析[J]. 水利學報，2010，41(12)：1387-1391.

[14] 封國林，龔志強，董文杰. 基于啟發式分割算法的氣候突變檢測研究[J]. 物理學報，2005，54(11)：5494-5499.

[15] 龔志強. 基于非線性時間序列分析方法的氣候突變檢測研究[D]. 揚州：揚州大學，2006.

[16] 拜有存. 基于差異信息理論的渭河年徑流變異點診斷[J]. 人民黃河，2009，31(7)：29-30.

[17] BAI H，MU X M，WANG F，et al. Analysis on evolution law of meteorological and hydrological drought and wetting[J]. Agricultural Research in the Arid Areas，2012，30:237-241.

[18] ZHANG H L，CHEN Y， REN G X. The characteristics of precipitation variation of Weihe River Basin in Shaanxi Province during recent 50 years[J]. Agricultural Research in the Arid Areas，2008，26:236-240.

[19] 張淑蘭，王彥輝，于澎濤,等.定量區分人類活動和降水量變化對涇河上游徑流變化的影響[J]. 水土保持學報，2010，24(4)：53-58.