信 亮

(遼寧省河庫(kù)管理服務(wù)中心(遼寧省水文局)，沈陽(yáng) 110001)

0 前 言

排列熵作為一種新型水文系統(tǒng)序列的突變檢測(cè)方法，主要特征是計(jì)算便捷、抗噪聲強(qiáng)，快速檢測(cè)等，在氣象學(xué)、信號(hào)學(xué)、生物學(xué)等方面有諸多應(yīng)用[1]。文章在渭河流域的水文時(shí)間序列突變研究中加入排列熵，分析水文系統(tǒng)的動(dòng)力學(xué)結(jié)構(gòu)突變，旨在為水文序列的突變分析提供一種新的研究方法。

1 滑動(dòng)PE基本原理及步驟

1.1 基本原理

排列熵作為解析系統(tǒng)復(fù)雜性的度量方法，其基本原理為：

1)設(shè)有一維水文時(shí)間序列{x(i)，i=1，2，…，n}，空間重構(gòu)后的一維向量為：

X(i)=[x(i),x(i+τ)，…，x(i+m(m-1))τ]

(1)

式中：m為嵌入維數(shù)；τ為滯后時(shí)間。

2)將第m個(gè)重構(gòu)分量[x(i),x(i+τ)，…，x(i+m(m-1))τ]從小到大進(jìn)行次序排序，得到：

[x(i+(k1-1))≤x(i+(k2-1))
τ≤…x(i+(kp-1))τ]

(2)

排序時(shí)，若x(i+(ki1-1))τ=x(i+(ki2-1)τ)，根據(jù)k值決定順序，則任意向量X(i)的符號(hào)序列為：

B(t)=[k1，k2,…,km]
1≤t≤n-m+1

(3)

3)將相同次序的符號(hào)序列B(t)=[k1,k2,…,km]進(jìn)行歸類，根據(jù)出現(xiàn)個(gè)數(shù)計(jì)算在n-m+1組中的發(fā)生概率P1,P2，…，Pl[2]。

4)按照Shannon信息熵的形式定義概率P1,P2，…，Pm-p+1，得到：

(4)

一般采用In(n-m+m)對(duì)Hp(m)標(biāo)準(zhǔn)化處理，表達(dá)式為：

(5)

Hp小，序列規(guī)則，復(fù)雜度低；Hp大，序列接近隨機(jī)，復(fù)雜度高。

1.2 徑流突變分析步驟

1)選擇子序列長(zhǎng)度n。

2)選擇滑動(dòng)步長(zhǎng)h，滑動(dòng)窗口n，在時(shí)間序列中第1個(gè)數(shù)據(jù)以滑動(dòng)步長(zhǎng)h逐步滑動(dòng)，選取子序列，直至原始序列結(jié)束，通常h=1。

3)選取嵌入維數(shù)m和滯后時(shí)間τ，對(duì)各子序列進(jìn)行PE分析，計(jì)算Hp值；根據(jù)步長(zhǎng)h羅列Hp序列。

4)結(jié)合Hp序列確定水文突變點(diǎn)[3]。

2 滑動(dòng)PE在徑流突變分析中的應(yīng)用

2.1 滑動(dòng)PE計(jì)算

根據(jù)渭河流域張家山站1960-2002年(共41年)的日徑流量資料x(chóng)(i),(i=1,2,…14965)，先選取子序列長(zhǎng)度 n=5a(a取365d)，滑動(dòng)步長(zhǎng)h為1天，嵌入維數(shù)m=5,計(jì)算的徑流PE值結(jié)果如圖1所示。

圖1 渭河流域張家山站日徑流PE值

由圖可知，徑流序列曲線在h=3650附近出現(xiàn)明顯變化，隨后曲線逐步上升，表明張家山站日徑流在h=3630(1971年)徑流出現(xiàn)較大變化，分析原因是自1970年以來(lái)該區(qū)域推廣的水土保持措施大范圍應(yīng)用，下墊面變化帶來(lái)的流域徑流量變化[4]。

2.2 滑動(dòng)PE準(zhǔn)確性影響因素

1)滑動(dòng)步長(zhǎng)：

取滑動(dòng)步長(zhǎng)h=a，子序列長(zhǎng)度n=5a，嵌入維數(shù)m=5，計(jì)算結(jié)果見(jiàn)圖2。由圖可知，當(dāng)子序列遠(yuǎn)>滑動(dòng)步長(zhǎng)時(shí)，滑動(dòng)步長(zhǎng)的小幅度擴(kuò)大可使曲線平滑，計(jì)算量變小，結(jié)果便于讀取觀測(cè)，但水文序列的突變點(diǎn)位置不變。

圖2 張家山站日徑流PE值(h=a，n=5a，m=5)

2)子序列長(zhǎng)度：

滑動(dòng)步長(zhǎng)取h=a，嵌入維數(shù)取m=5，子序列不同長(zhǎng)度n=a，3a，4a，6a，10a時(shí)的PE突變檢測(cè)結(jié)果見(jiàn)圖3。

(b)n=4a

(c)n=a

(d) n=6a

(e)n=10a

從圖3(a)、(b)中可知當(dāng)n=3a，4a時(shí)，突變點(diǎn)基本保持一致，突變時(shí)間點(diǎn)均為1971年。為了檢驗(yàn)子序列長(zhǎng)度的影響，分別進(jìn)行極大、極小探討。當(dāng)n=a時(shí)，圖(c)的時(shí)間序列在1972年附近出現(xiàn)突變，曲線彎折多變，上升過(guò)程出現(xiàn)較多拐點(diǎn)，影響序列突變點(diǎn)的有效判斷；當(dāng)n=6a時(shí)，圖(d)的序列突變點(diǎn)向左移動(dòng)，突變點(diǎn)為1970年，與小序列所測(cè)結(jié)果不同；當(dāng)n=10a時(shí)，圖(e)序列PE曲線逐漸平穩(wěn)，無(wú)突變點(diǎn)出現(xiàn)，由此可知子序列長(zhǎng)度的選值對(duì)最終突變點(diǎn)的準(zhǔn)確預(yù)測(cè)具有決定性作用[5]。

3)維數(shù)變化影響：

子序列長(zhǎng)度n=5a，滑動(dòng)步數(shù)取h=a，分別給出了P=3,4,5時(shí)的分析結(jié)果，如圖4所示。

(a)P=3

(b)P=4

(c)P=2

(d)P=15

由圖4可知，在當(dāng)嵌入維數(shù)P=3，4時(shí)，序列的整體變化曲線未出現(xiàn)明顯變化，突變的時(shí)間點(diǎn)基本相同，在1971年處產(chǎn)生徑流突變，當(dāng)P=2時(shí)，PE值Hp較小，且小范圍內(nèi)波動(dòng)，不具備排列統(tǒng)計(jì)意義；當(dāng)P=15時(shí)，排列熵值變幅較大，在多個(gè)區(qū)間明顯轉(zhuǎn)折，突變點(diǎn)結(jié)果檢測(cè)失準(zhǔn)。當(dāng)P=5左右時(shí)，PE值Hp變化范圍適中，曲線較平滑，突變檢測(cè)診斷較為準(zhǔn)確。

3 結(jié) 論

文章基于熵理論的基本原理和解析步驟，結(jié)合渭河張家村站1960-2000年日徑流序列資料，利用排列熵的滑動(dòng)PE方法計(jì)算徑流時(shí)間序列突變點(diǎn)，得出1971年為該流域的徑流水文序列突變點(diǎn)，同時(shí)開(kāi)展了子序列長(zhǎng)度，滑動(dòng)步長(zhǎng)，嵌入維數(shù)的研究，綜合認(rèn)為n=5a，h=a，m=5時(shí)的診斷結(jié)果最準(zhǔn)確，具備一定的實(shí)際應(yīng)用價(jià)值。