基于EMD的岷江上中游流域流量特性分析

王俊鴻，覃光華,2，童 旭

(1.四川大學水利水電學院，成都 610041；2.水力學與山區(qū)河流開發(fā)保護國家重點實驗室，成都 610041)

山區(qū)和丘陵地帶河流的徑流變化過程往往較平原地區(qū)更為復雜。受人為或自然等眾多因素的影響，這一地區(qū)河流的徑流變化越來越顯著。研究徑流在長時間序列中的變化規(guī)律，有利于進行水文規(guī)律的預測預報，同時指導區(qū)域河流水量調(diào)度。

國內(nèi)外研究長時間序列尺度的徑流變化規(guī)律常用的主要有小波形分析、多維混沌理論、BP激勵函數(shù)神經(jīng)網(wǎng)絡(luò)模型和經(jīng)驗模態(tài)分解(EMD)等方法[1]。小波形分析[2]是運用最為廣泛的一種數(shù)學-水文耦合方法，理論上是一種基于時頻局部化的信號，進入新世紀以來在流量、水位、輸沙率等多個水文要素時間序列研究中均取得了一定的研究成果。而具有優(yōu)良自適應性的BP人工神經(jīng)網(wǎng)絡(luò)方法[3]，被廣泛應用于更復雜的水文變化特征中，例如水文與氣象結(jié)合的二重尺度交互變化規(guī)律[3]。多維混沌理論方法——基于雜亂無序的隨時間變化的大量水文數(shù)據(jù)，分解出非線性特征并對其進行分析，從而進一步創(chuàng)建序列[4]。

經(jīng)驗模態(tài)分解方法(EMD)雖然沒有以上三種方法使用廣泛，但也常被用來分析非線性和非平穩(wěn)時間序列[5]，它能夠把復雜的時間序列組拆分成若干個本征模函數(shù)(Intrinsic Mode Function)[6]，即關(guān)于序列自身形態(tài)模式的函數(shù)。這種分解方式具有很強的自適應性，對數(shù)據(jù)本身特征進行針對性的分解，這在一定程度上提高了分解效率。岳相臣[7]在傳統(tǒng)EMD方法的基礎(chǔ)上利用端點效應和模態(tài)混疊提出新的改進方法，攻破了信號因接受頻率而受限的問題；于本成[8]在進行熵性能研究時用到了改進了的EMD算法。另一方面，EMD方法用于水文預測時，其精度比單一神經(jīng)網(wǎng)絡(luò)更高[9,10]，因而EMD具有一定的推廣價值，目前基于EMD的優(yōu)化算法已經(jīng)在其他各學科得到有效應用[11]，在水文領(lǐng)域也有涉及[12]。

岷江上中游，作為典型的山區(qū)丘陵河流，基于EMD方法的規(guī)律研究較少。本文運用EMD方法對岷江上中游流域6個代表水文站觀測到的1956～2016年平均流量序列分別進行多時間尺度分析，以此獲得徑流序列的時間變化規(guī)律，并對其進行特征分析。

1 研究方法

EMD的基本思想是將一個頻率不規(guī)則的波化為若干個單一頻率的波(IMF)與殘波的和的形式，本質(zhì)上是對數(shù)據(jù)序列或信號的平穩(wěn)化處理。在分解過程中，每個IMF均需符合以下2個性質(zhì)：一是信號在時間序列上任意一個特定長度區(qū)間段的最大值或最小值的數(shù)目和過零點數(shù)目相等，或滿足指定的某一標準(最大不超過1)，這個條件是為適應數(shù)據(jù)在頻率區(qū)間上的密集性需求；二是由特定長度區(qū)間段的若干最大值構(gòu)成的上圍線和由若干最小值構(gòu)成的下圍線的中值為零。其目的是為使序列的波形趨于對稱，從而方便對瞬時頻率進行操作。

根據(jù)上述IMF的兩個性質(zhì)，對原始時間序列進行分解時分別采用由最大值產(chǎn)生的局部上包絡(luò)線和最小值產(chǎn)生的局部下包絡(luò)線確定若干個極值。當極值確定后，利用三次樣條曲線將所有的局部最大值和局部最小值分別擬合成上、下兩條包絡(luò)線，并且需要保證所有的時間序列數(shù)值都介于這兩條線之間，取m1(t)為這兩條包絡(luò)線的中位值。

初始目標值通過x(t)與m1(t)相減得到：

h1(t)=x(t)-m1(t)

(1)

若一系列的h1(t)不平穩(wěn)，則視為不滿足作為IMF分量的條件，因此要對h1(t)重復以上過程p次進行篩分，有：

hp(t)=hp-1(t)-mp(t)

(2)

直到hp(t)符合IMF分量的條件為止，此時將hp(t)作為確定的第一個IMF分量，記為imf1(t)。為簡化篩選程序，可以利用限制標準偏差STDV的值來判斷每次篩分結(jié)果是否為IMF分量，其定義為：

(3)

式中：T表示待研究的時間序列總長度。

霍譯：The words were as oil upon fire.The old lady blazed.“Who gave orders for those coffins…”

這樣得到的imf1(t)應是信號的特征時間尺度最小的高頻IMF分量，然后進行下一步驟：

x(t)-imf1(t)=r1

(4)

因為剩余序列r1中也許會含有比imf1(t)的特征時間尺度還大的IMF分量，所以又要把 當成新的序列進行上述篩選過程，以求得第二個IMF分量imf2(t)。

這便是EMD分解核心步驟——篩分過程，其目的是為消除畸形波并平滑細微的彎曲，保存大的振幅，簡化波形；同時使序列曲線變成圍繞零均值線的局部極大值和極小值基本對稱的波形。

為使得篩分過程不至于持續(xù)進行下去，我們還有這樣的標準[13]，在EMD過程中當滿足0.2

圖1 EMD模型分解步驟Fig.1 EMD model decomposition steps

各個IMF分量分別表示不同的特征時間尺度，用連續(xù)兩個極值之間的時間間隔作為本征模態(tài)的時間尺度，不僅能清晰地分辨出不同尺度的震蕩分量，還可用于數(shù)值全為正或負的沒有穿越零點的序列。根據(jù)以上篩分過程，直到imf1(t)或rn(t)的值小于預先設(shè)計的定值，或者rn(t)不能進一步分解為IMF分量，以單調(diào)函數(shù)的形式截止。這里的rn(t)表示原時間序列的趨勢，記為趨勢項RES[14]。最后，原時間序列可以改寫為：

(5)

2 研究區(qū)概況

岷江是長江上游的重要支流，在樂山城區(qū)以上為岷江的上中游，以大渡河為正源，其流域面積為126 478 km2，河長為1 129 km，落差達到2 795 m，這一范圍大致位于龍門山、岷山等眾多盆地邊緣山脈與成都平原之間組成藏區(qū)高原東緣地形陡變帶內(nèi)，平均海拔2 640 m，河道不僅比降大，且兩岸以砂壤土質(zhì)為主。又因其位于我國沿海季風濕潤區(qū)向高原高寒區(qū)的過渡地帶，氣候條件錯綜復雜，流域年平均降水500～950 mm，多年平均氣溫12～20 ℃。岷江上中游出口斷面處徑流較為充沛，據(jù)有關(guān)水文資料多年平均流量為1 000 m3/s，豐水期流量占全年流量的75%，而枯水期流量僅占全年流量的25%，豐枯季界限分明。并且，岷江上中游流域因其水系發(fā)達，從古至今一直以來也是成都平原的水源補給河流。

3 實例分析

青衣江和大渡河是岷江中游以上重要的一級支流，并且本次選擇的這6個站點有上游的大金站、都江堰站，中游的瀘定站、彭山站，下游的夾江站和沙灣站，既有山區(qū)站點也有丘陵站點，同一河道上最近的兩站點間隔在100 km以上，覆蓋面較廣，同時這些站點都有超過60 a的長時間序列資料，因此都具有較好的代表性，通過這些代表站點的水文資料能較好地反映研究流域的流量特征。各代表站在研究流域上的分布見圖2。

圖2 各代表站在研究流域上的分布Fig.2 The distribution of representative gaging stations on the study basin

針對上述站點1956-2016 a實測徑流資料，運用EMD方法對年平均流量進行分析[15]。經(jīng)過EMD拆分過程，年均流量將得到若干個IMF分量(分別記作imfi)及一個趨勢項RES。見圖3，圖3(a)～圖3(f)中signal表示原始時間序列，imf1/2/3和res.表示原始時間序列經(jīng)過EMD分解得到的3個ifm分量和1個趨勢項[16]，各自的平均周期如表1所示。

表1 各站IMF平均周期Tab.1 Average period of IMF in each station

從表1可以看出，較高頻率波動的 所反映的周期基本相同[17]，各站點imf1顯示周期集中在2.5～3.5 a，imf2顯示周期集中在5.0～7.0 a。而低頻波動的imf3顯示的周期差異則較大，最小的都江堰水文站周期為10.0 a，最大的瀘定水文站周期為17.4 a[18]。同時，岷江上中游各站點年均流量的EMD趨勢均反映出近乎直線型的下降趨勢[19]，其中大金站的年均流量在2007年開始停止下降，此后近十年趨于平穩(wěn)。這是由于大金站地處大渡河上游高原山區(qū)，在全球氣候變化及地區(qū)地理氣候、和人類環(huán)保意識加強的多重影響下，該地區(qū)的流量在經(jīng)歷了前幾十年的削減后，在近十年得到有效緩解[20]，而這一緩解現(xiàn)象在未來將在岷江流域的其他站點重現(xiàn)。

岷江上中游流域徑流量的下降趨勢通過樂山市區(qū)的五通橋水文站1956-2016年觀測得到的數(shù)據(jù)也能夠發(fā)現(xiàn)。五通橋水文站地處岷江中游干流上，是岷江上中游流域的出口斷面，因而該水文站近半個世紀的年均流量觀測數(shù)據(jù)能較好地反映岷江上中游流域徑流量的發(fā)展趨勢。將所研究的6個水文站和五通橋水文站1956-2016年觀測計算得到的年均流量用折線圖表示，見圖3，圖3中虛線為趨勢線，趨勢線斜率均為負值。

圖3 岷江上中游典型水文站及五通橋水文站1956-2016年均流量Fig.3 The average flow of hydrological station of the upper and middle reaches of the Minjiang River and Wutongqiao Hydrological Station from 1956 to 2016

imfi(i=1,2,3)所對應的各站點標準偏差見表2。

表2 imfi所對應各站點的標準偏差Tab.2 Standard deviation of each station corresponding to imf1/2/3

通過計算imfi(i=1,2,3)各站點的STDV，分別為0.38、0.69、2.45，這反映了6個站點在imf1，imf2，imf33種情況下均存在離散情況，并且隨著頻率波動的減小，其離散程度相應增大。這是因為大渡河、青衣江、岷江三條河流的流域下墊面存在差異，大渡河流域主要以山地為主，青衣江和岷江主要以丘陵為主，這決定了河道的匯流方式和流量的年際分布，所以在高頻率波動的imf1，imf2中地勢差異起伏大的站點STDV普遍高于地勢平緩的站點；而低頻率波動的imf3中由于受氣象和下墊面等因素的耦合作用[21]，這種差異被放大，同時受人為干預或高原復雜地下水演變的影響，都江堰站和大金站顯著偏小，導致imf3的離散程度最大。為試圖通過imfi之間的STDV規(guī)律來預測岷江上中游流域流量的趨勢(圖4)，可以發(fā)現(xiàn)STDV(imf2)約為STDV(imf1)的2倍，STDV(imf3)約為STDV(imf1)與STDV(imf2)和的2倍。

圖4 岷江上中游代表水文站年均流量EMD結(jié)果Fig.4 EMD results of annual average flow of representative hydrological stations on the upper and middle reaches of Minjiang River

4 結(jié) 語

采用經(jīng)驗模態(tài)分解(EMD)對岷江上中游流域6個代表水文站的1956-2016年均流量進行分析，并通過樂山市區(qū)五通橋水文站的同時期資料加以復核，得到相應的周期及趨勢和統(tǒng)計規(guī)律，結(jié)論如下。

(1)岷江上中游流域年均流量存在2.5～3.5 a、5.0～7.0 a的短周期和10.0～18.0 a的長周期，且近半個世紀來岷江上中游徑流量持續(xù)減少。

(2)除水文規(guī)律外，還發(fā)現(xiàn)隨著頻率波動的減小(周期增大)，多個站點的年均流量標準偏差相應增大，且低頻周期的離散程度在高頻周期的基礎(chǔ)上并不是逐級增加，而是呈倍數(shù)放大，且倍比與高頻周期的離散程度密切有關(guān)，STDV(imf2)約為STDV(imf1)的2倍，STDV(imf3)約為STDV(imf1)與STDV(imf2)和的2倍。這一規(guī)律過去很少有研究涉及，這表征了水文規(guī)律的預測不僅可以從單一序列自身規(guī)律著手，還可以從一個流域的多個序列為對象，通過他們之間潛在的關(guān)系來進行預測，這符合水文實際情況，同時也使預測結(jié)果更客觀全面。

但由于本研究僅就岷江上中游流域作為研究區(qū)，這一地區(qū)的地形、氣候、植被[22]等變化呈規(guī)律性階梯分布，并且受人為影響明顯，因此這一結(jié)果在未對更多流域進行驗證時僅適用于岷江流域的諸多河流。望今后更多學者在基于標準偏差的多站點水文規(guī)律研究中有所涉足。