花園口站年最大洪峰流量時(shí)間序列的HHT分析

任 健，史 紅 玲

花園口站年最大洪峰流量時(shí)間序列的HHT分析

任 健1，史 紅 玲2

（1.中國(guó)水利水電科學(xué)研究院，北京 100048;2.國(guó)際泥沙研究培訓(xùn)中心，北京 100048）

Hilbert-Huang變換能夠定量描述非線性、非平穩(wěn)復(fù)雜時(shí)間序列的時(shí)頻特性，較傳統(tǒng)分析方法更具優(yōu)勢(shì)。通過對(duì)時(shí)間序列進(jìn)行EMD分解，得到變化過程的內(nèi)在模態(tài)函數(shù)和趨勢(shì)項(xiàng)函數(shù)，而后對(duì)各內(nèi)在模態(tài)函數(shù)進(jìn)行Hilbert-Huang變換，從而揭示出時(shí)間序列的多時(shí)間尺度特征。以黃河花園口站1952－2009年的年最大洪峰流量時(shí)間序列為例，對(duì)其進(jìn)行多時(shí)間尺度分析，得到不同波動(dòng)周期的振蕩分量及趨勢(shì)分量，具體分析了各分量的變化特征。結(jié)果表明，花園口年最大洪峰流量變化過程中存在準(zhǔn)3.2 a、準(zhǔn)6.4 a、準(zhǔn)11.8 a和準(zhǔn)31.0 a周期的波動(dòng)，其中準(zhǔn)3.2 a和準(zhǔn)6.4 a的周期波動(dòng)是引起原序列波動(dòng)的主要原因，近60年來花園口年最大洪峰流量變化呈遞減趨勢(shì)，由此揭示了年最大洪峰流量變化過程的多時(shí)間尺度特征。在此基礎(chǔ)上，探討了各波動(dòng)分量變化的影響因素，其變化與大氣低頻振蕩、ENSO、太陽(yáng)活動(dòng)及氣候變遷等因素有關(guān)。

Hilbert-Huang變換;經(jīng)驗(yàn)?zāi)B(tài)分解;多時(shí)間尺度;年最大洪峰流量

黃河花園口站是黃河下游重要的控制性監(jiān)測(cè)斷面，其控制流域面積達(dá)73萬(wàn)km2，約占黃河流域面積的97%，分析花園口站的水沙變化情勢(shì)是研究黃河下游河道演變的前提和基礎(chǔ)。長(zhǎng)期以來，關(guān)于黃河下游水沙變化特征的研究很多［1－4］，但多依據(jù)傳統(tǒng)的數(shù)理統(tǒng)計(jì)或簡(jiǎn)單的數(shù)據(jù)平滑等方法，不能揭示水沙變化的多時(shí)間尺度特性［5］。而Fourier分析和小波分析等方法雖然可以進(jìn)行頻譜分析，但其方法本身存在一些不足和缺點(diǎn)［6，7］。希爾伯特—黃變換［6，7］（Hilbert-Huang Transform，HHT）技術(shù)是近年來發(fā)展的一種最新的處理非線性、非平穩(wěn)時(shí)間序列的方法，它直接由時(shí)間序列本身自適應(yīng)地構(gòu)造出基函數(shù)，進(jìn)而分離得到不同尺度的特征分量，而不用預(yù)先設(shè)定基底函數(shù)再展開，因而該法可更準(zhǔn)確反映序列的內(nèi)在規(guī)律特征。本文運(yùn)用H HT技術(shù)，對(duì)黃河下游花園口站1952－2009年的實(shí)測(cè)年最大洪峰流量序列進(jìn)行多時(shí)間尺度分析，診斷花園口站年最大洪峰變化過程中蘊(yùn)含的多時(shí)間尺度振蕩結(jié)構(gòu)和特征，揭示其在不同尺度和層次上的變化規(guī)律，以期為黃河下游的水資源利用、河床演變、河道整治及防洪調(diào)度等提供技術(shù)支撐和科學(xué)依據(jù)。

1 研究方法和數(shù)據(jù)資料

經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）方法［6，7］是一種全新的處理非線性、非平穩(wěn)數(shù)據(jù)序列的方法，它將時(shí)間信號(hào)f（t）分解成一系列本征模態(tài)函數(shù)（Intrinsic Mode Function，IMF），每個(gè)IMF分量都具有如下特征：從全局特性看，極值點(diǎn)數(shù)必須和過0點(diǎn)數(shù)一致或至多相差一個(gè);在某個(gè)局部點(diǎn)，極大值包絡(luò)和極小值包絡(luò)在該點(diǎn)的值的算術(shù)平均和是0。

EMD分解過程為：找出序列f（t）所有的極大值和極小值點(diǎn)，分別用三次樣條函數(shù)擬合成上下包絡(luò)線，得到平均包絡(luò)線m1，將原序列減去m1得到去掉低頻序列的新序列h1。一般h1并不能滿足IMF的兩個(gè)特征要求，仍是非平穩(wěn)的，因而可以多次重復(fù)上述過程，使平均包絡(luò)線趨近于0，得到第一個(gè)IMF分量c1，c1代表原始序列中最高頻的分量。即：

對(duì)r1（t）繼續(xù)進(jìn)行上述分解，直到剩余部分為單一信號(hào)或其值小于預(yù)先給定的值，分解結(jié)束。原始時(shí)間序列f（t）可以表示為：

式中：cj（t）為各IMF分量;Res為趨勢(shì)項(xiàng);t為時(shí)間變量。

對(duì)每個(gè)IMF分量cj（t）進(jìn)行Hilbert變換，即：

p為Cauchy主值，cj（t）和bj（t）可以構(gòu)成一個(gè)復(fù)序列：

而由瞬時(shí)位相求得相應(yīng)的瞬時(shí)頻率為ωj（t）＝dθj（t）/dt。瞬時(shí)振幅aj（t）具有明確的物理意義，它表示某IMF分量的波動(dòng)振蕩能量，aj（t）的峰值對(duì)應(yīng)于該分量的主要強(qiáng)振蕩，而aj（t）的低谷則對(duì)應(yīng)著該分量的小幅變化，各點(diǎn)瞬時(shí)振幅的平均值即為平均振幅，表示該分量振蕩強(qiáng)弱的平均情況。瞬時(shí)頻率ωj（t）描述的是IMF分量在頻域上的變化特征，它給出了某IMF分量在某個(gè)時(shí)刻出現(xiàn)波動(dòng)的頻率，與該時(shí)刻的波動(dòng)周期互為倒數(shù)關(guān)系，因而其物理意義很明晰。由此可知，中心頻率即為某IMF分量各點(diǎn)瞬時(shí)頻率的平均值，其倒數(shù)即為平均周期。EMD方法不可避免地存在邊界效應(yīng)，本文采用鏡像對(duì)稱延伸方法對(duì)邊界效應(yīng)進(jìn)行了處理［8］。

采用黃河下游花園口水文站1952－2009年的年最大洪峰資料（來源于黃河下游河床演變基本資料匯編及黃河泥沙公報(bào)）獲得近60年來花園口站年最大洪峰流量的變化過程（圖1）。通過圖1只能直觀判斷近60年花園口站年最大洪峰流量變化的趨勢(shì)，并不能深入、細(xì)致地揭示最大洪峰流量的演變周期及多時(shí)間尺度特征，因此有必要運(yùn)用基于EMD的HHT技術(shù)對(duì)其進(jìn)行多時(shí)間尺度分析，從有限的數(shù)據(jù)序列中挖掘盡量多的有用信息，以便更深刻認(rèn)識(shí)洪水變化的內(nèi)在規(guī)律和特性。

圖1 花園口站年最大洪峰流量時(shí)間序列Fig.1 Annual maximum peak flow series at Huayuankou Station

2 年最大洪峰流量變化的EMD分解

運(yùn)用EMD方法對(duì)花園口站1952－2009年的年最大洪峰流量時(shí)間序列進(jìn)行多時(shí)間尺度分解，得到4個(gè)本征模態(tài)函數(shù)（c1～c4分量）和一個(gè)殘余趨勢(shì)項(xiàng)（Res分量）（圖2）。

圖2 花園口站年最大洪峰流量時(shí)間序列的EMD分解Fig.2 Four decomposed IMFs and a trend of maximum peak flow at Huayuankou Station

由圖2可知：1）花園口站實(shí)測(cè)年最大洪峰流量變化過程是非線性和非平穩(wěn)的，是多種波動(dòng)成分共同作用的結(jié)果，其可分解為4個(gè)具有不同波動(dòng)周期的振蕩分量和1個(gè)趨勢(shì)分量，反映出花園口站年最大洪峰流量變化具有復(fù)雜的多時(shí)間尺度特性。2）c1分量具有準(zhǔn)2～4 a的波動(dòng)周期，且不同時(shí)期的振幅也不同，呈現(xiàn)出明顯的階段性特征。1964年前，洪峰流量的變化幅度一般在5 000～6 000 m3/s，振幅遠(yuǎn)大于之后其他時(shí)段;1964－1980年，波動(dòng)幅度急劇降低，變化范圍平均為500～1 000 m3/s;1980－1986年，波動(dòng)幅度有所增大，一般為500～1 000 m3/s;1986－1999年，波動(dòng)周期有所增大，波動(dòng)幅度再次降至500～1 000 m3/s;1999年以來波動(dòng)幅度進(jìn)一步減小，平均在300 m3/s左右。3）c2分量大體具有準(zhǔn)5～7 a波動(dòng)周期，并且不同時(shí)期的振幅也不同，呈現(xiàn)出明顯的階段性特征。其中1964年之前，波動(dòng)幅度呈急劇衰減的趨勢(shì)，波動(dòng)幅度由1950s的7 000 m3/s降至1 000 m3/s;1964－1973年，波動(dòng)幅度保持在500～1 000 m3/s;1974－1986年，波幅有所恢復(fù)，增至2 000～3 000 m3/s;1986－2000年，波動(dòng)幅度再次降至500～1 000 m3/s;2000年以來，波動(dòng)周期和波動(dòng)幅度有所增大，約為2 500 m3/s。4）c3分量大致具有以準(zhǔn)9～12 a為主的波動(dòng)周期。其波動(dòng)幅度在1964年前為2 000～2 500 m3/s;1964－1974年，波動(dòng)幅度降為500～1 000 m3/s;1974－1986年，波動(dòng)幅度為500～800 m3/s;1986－2000年，波動(dòng)幅度平均約為1 000 m3/s;2000年后，波動(dòng)周期變長(zhǎng)，波動(dòng)幅度也有所增大。5）c4分量具有準(zhǔn)30 a的波動(dòng)周期，整個(gè)時(shí)段內(nèi)波動(dòng)幅度基本保持不變，平均約為2 000 m3/s，其衰減程度很小。6）Res分量顯示的是花園口站年最大洪峰流量的整體變化趨勢(shì)，1950s以來該站年最大洪峰流量整體呈急劇衰減的趨勢(shì)。

需要指出的是，花園口站年最大洪峰流量的變化趨勢(shì)項(xiàng)中可能仍含有屬于更長(zhǎng)周期（更小頻率）波動(dòng)的組分，而限于觀測(cè)資料時(shí)段長(zhǎng)度，這種波動(dòng)的周期、頻率和振幅尚不能從趨勢(shì)項(xiàng)Res分量中有效分解出來。

3 年最大洪峰流量變化的HHT分析

將經(jīng)過EMD分解得到的花園口站年最大洪峰流量的各IMF分量進(jìn)行Hilbert變換，得到各IMF分量的時(shí)幅和時(shí)頻譜圖（圖3），同時(shí)，計(jì)算IMF分量的方差貢獻(xiàn)并統(tǒng)計(jì)其Hilbert變換后的特征值（表1）。從圖3可以看出，各IMF分量的頻率并不是一個(gè)恒定值，而是圍繞中心頻率上下波動(dòng)，頻率越高的IMF分量，其波動(dòng)的幅度越大。盡管IMF分量的頻率圍繞中心頻率波動(dòng)，但波動(dòng)的范圍通常有限，相互之間少有交叉重疊現(xiàn)象，保持了一種界限較為清晰的分布特征。由時(shí)頻圖和時(shí)幅圖也可以看出，周期越長(zhǎng)，在該周期上瞬時(shí)頻率的波動(dòng)就越小，振幅的變化也越小，即對(duì)應(yīng)于大尺度周期，年最大洪峰流量的變化相對(duì)穩(wěn)定，但波動(dòng)能量較小。

圖3 基于HHT分析的花園口站年最大洪峰流量序列時(shí)幅和時(shí)頻譜Fig.3 The time-amplitude and time-frequency relationship of annual maximum peak flow at Huayuankou Station after Hilbert-Huang Transform

表1 花園口站年最大洪峰流量的各IMF分量Hilbert變換后的數(shù)字特征值Table 1 Statistics of IMFs of annual maximum peak flow at Huayuankou Station after Hilbert-Huang Transform

從表1可知，對(duì)于花園口站年最大洪峰流量序列而言，從第1個(gè)分量到第4個(gè)分量，其中心頻率分別為0.31 a－1、0.16 a－1、0.08 a－1和0.03 a－1，而平均周期分別為3.2 a、6.4 a、11.8 a和31.0 a，總體上表現(xiàn)出中心頻率由高到低的變化特點(diǎn)，而平均周期則由短變長(zhǎng)。振幅反映了各分量波動(dòng)能量的大小。分析各個(gè)分量的平均振幅和最大振幅，中短周期波動(dòng)的振幅一般較大，而長(zhǎng)周期波動(dòng)的振幅一般較小，呈現(xiàn)出頻率越高（或周期越短）振幅越大的特點(diǎn)。方差貢獻(xiàn)代表了各分量對(duì)原序列信號(hào)波動(dòng)趨勢(shì)的影響程度，方差貢獻(xiàn)越大，說明該分量的波動(dòng)對(duì)原時(shí)序變化的影響程度越強(qiáng)。不難發(fā)現(xiàn)，中短周期波動(dòng)的方差貢獻(xiàn)率遠(yuǎn)大于長(zhǎng)周期波動(dòng)，這表明中短周期的波動(dòng)分量是引起年最大洪峰流量變化的主要原因，而長(zhǎng)周期波動(dòng)分量對(duì)年最大洪峰流量變化過程的影響主要體現(xiàn)在對(duì)其整體發(fā)展趨勢(shì)的控制性作用。

4 討論

黃河下游大的洪峰流量多發(fā)生于伏汛、秋汛期間，是由暴雨引起，因而年最大洪峰流量是由年極端降水事件造成的，二者具有天然的、復(fù)雜的聯(lián)系，年極端降水事件的波動(dòng)影響著年最大洪峰流量的變化過程。短時(shí)期極端降水事件的發(fā)生受到大氣運(yùn)動(dòng)、氣象等因素的影響，但長(zhǎng)時(shí)期極端降水事件發(fā)生的頻率（周期）與強(qiáng)度受到全球及區(qū)域的物理、氣候系統(tǒng)變化的控制與影響。已有的研究表明平流層的大氣運(yùn)動(dòng)存在準(zhǔn)兩年周期的振蕩（QBO），而ENSO事件具有3.5 a和5～6 a周期等大氣低頻變化成分［9］。這些氣象及氣候變化事件的周期與花園口年最大洪峰流量變化過程中存在的準(zhǔn)3.2 a和準(zhǔn)6.4 a分量（c1、c2分量）的周期基本一致;太陽(yáng)黑子活動(dòng)具有準(zhǔn)11 a的周期［10］，這與花園口年最大流量過程中存在的準(zhǔn)11.8 a周期基本一致;我國(guó)氣溫變化存在準(zhǔn)30 a的冷暖周期［11］，這與花園口年最大流量過程中存在的準(zhǔn)31.0 a周期基本一致。因而，花園口年最大洪峰流量變化過程中蘊(yùn)含不同周期的波動(dòng)特征，一定程度上反映了大氣、氣象及氣候等因素對(duì)黃河流域水文要素可能存在的影響作用。

隨著黃河流域內(nèi)大型水庫(kù)調(diào)度運(yùn)營(yíng)和人類活動(dòng)影響程度的加劇，花園口站年最大洪峰流量各波動(dòng)分量的周期和振幅也發(fā)生不同程度的調(diào)整，如圖1中各分量波動(dòng)曲線及圖3中時(shí)幅、時(shí)頻譜圖所示。對(duì)于中短周期分量（c1～c3分量）而言，人類活動(dòng)的影響導(dǎo)致其周期在局部時(shí)段內(nèi)發(fā)生變化，同時(shí)振幅也發(fā)生階段性的變化;并且各分量隨著其中心頻率的遞減，周期和振幅的調(diào)整變化也趨緩。而對(duì)于長(zhǎng)周期分量（c3分量）而言，人類活動(dòng)的影響對(duì)其周期和振幅的影響作用則不太明顯。這在一定程度上說明，人類活動(dòng)對(duì)不同時(shí)間尺度下黃河流域水文過程的影響機(jī)制不同：人類活動(dòng)對(duì)水文過程的中短周期變化的影響作用較大，通常能夠較大程度地影響波動(dòng)幅度，并在有限時(shí)段內(nèi)調(diào)整波動(dòng)周期;而人類活動(dòng)對(duì)水文過程長(zhǎng)周期變化的影響作用較小，對(duì)長(zhǎng)時(shí)間尺度水文過程的影響可能很大程度上受到天然降水波動(dòng)及氣候變化等自然因素的制約。而對(duì)于包括人類活動(dòng)影響在內(nèi)的各種因素對(duì)不同時(shí)間尺度下黃河流域水文變化過程的影響機(jī)制仍需深入探討。

5 結(jié)論

利用基于EMD的HHT方法對(duì)1952—2009年的黃河下游花園口站年最大洪峰流量序列進(jìn)行多時(shí)間尺度分析，研究結(jié)果表明：1）花園口站年最大洪峰流量變化過程包含4個(gè)IMF分量和1個(gè)單調(diào)遞減的趨勢(shì)分量，反映了花園口年最大洪峰流量變化過程存在準(zhǔn)3.2 a、準(zhǔn)6.4 a、準(zhǔn)11.8 a和準(zhǔn)31.0 a周期的波動(dòng)特征，并且精確顯示出近60年來花園口年最大洪峰流量變化整體呈單調(diào)遞減趨勢(shì)。2）在年最大洪峰流量序列的多時(shí)間尺度結(jié)構(gòu)中，各IMF分量的中心頻率呈遞減的趨勢(shì)，平均周期呈遞增的趨勢(shì)，波動(dòng)幅度及方差貢獻(xiàn)率隨著中心頻率的遞減而遞減，揭示出中短周期波動(dòng)分量對(duì)原序列的波動(dòng)振蕩起著主要作用，而長(zhǎng)周期波動(dòng)分量則對(duì)原序列的整體變化趨勢(shì)起著控制性作用。3）從各IMF分量的波動(dòng)周期看，年最大洪峰流量變化過程中蘊(yùn)含的不同尺度的波動(dòng)特征可能與大氣低頻振蕩、ENSO、太陽(yáng)活動(dòng)及氣候變遷等因素有關(guān)，加之人類活動(dòng)的嚴(yán)重影響，使得各波動(dòng)分量的變化規(guī)律更加復(fù)雜。而對(duì)于包括人類活動(dòng)影響在內(nèi)的各種因素對(duì)黃河流域年最大洪峰流量變化過程的影響機(jī)制有待深入研究。

此外，HHT技術(shù)對(duì)于受干擾因素多、變化規(guī)律復(fù)雜、隨機(jī)成分大的非線性、非平穩(wěn)的水文序列具有較好的分析和診斷效果，在水文過程的多時(shí)間尺度分析、建模和預(yù)測(cè)中值得應(yīng)用和推廣。

［1］陳浩，周金星，陸中臣，等.黃河中游流域環(huán)境要素對(duì)水沙變異的影響［J］.地理研究，2002，21（2）：179－187.

［2］劉昌明，鄭紅星.黃河流域水循環(huán)要素變化趨勢(shì)分析［J］.自然資源學(xué)報(bào)，2003，18（2）：129－135.

［3］劉勇勝，陳沈良，陳小英，等.黃河中下游泥沙通量變化規(guī)律［J］.地理與地理信息科學(xué)，2006，22（4）：47－51.

［4］吳保生，張?jiān)h，申冠卿，等.維持黃河主槽不萎縮的水沙條件研究［M］.鄭州：黃河水利出版社，2010.9－15.

［5］張少文，丁晶，廖杰，等.基于小波變換的黃河上游天然徑流變化特性分析［J］.四川大學(xué)學(xué)報(bào)（工程科學(xué)版），2004，36（5）：32－37.

［6］HUANG N E，SHEN Z，LONG S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proc.R.Soc.Lond.A，1998，454：903－995.

［7］HUANG N E，SHEN Z，LONG S R.A new view of nonlinear water waves：The Hilbert spectrum［J］.Ann.Rev.Fluid.Mech.，1999，31：417－457.

［8］FLANDRIN P，RILLING G，GONCALVES G.Empirical mode decomposition as a filter bank［J］.IEEE Signal Processing Letters，2004，11（2）：112－114.

［9］彭永清，王盤興，吳洪寶.大氣低頻變化的分析與應(yīng)用［M］.北京：氣象出版社，1997.

［10］楊建平，丁永建，陳仁升.長(zhǎng)江黃河源區(qū)水文和氣象序列周期變化分析［J］.中國(guó)沙漠，2005，25（3）：351－355.

［11］孫林海，趙振國(guó).我國(guó)暖冬氣候及其成因分析［J］.氣象，2004，30（2）：57－60.

Analysis on Annual Maximum Peak Flow Series at Huayuankou Station Based on Hilbert-Huang Transform

REN Jian1，SHI Hong－ling2
（1.ChinaInstituteofWaterResourcesandHydropowerResearch，Beijing100048;2.InternationalResearchandTrainingCenteronErosionandSedimentation，Beijing100048，China）

Compared with traditional methods for time series，Hilbert-Huang Transform（HHT）has advantages to quantitatively describe time-frequency characteristics for nonlinear，unsteady and complex time series.Empirical Mode Decomposition method is adopted to analyzed time series，and Intrinsic Mode Function（IMF）and trend term are obtained.Then HHT is made for each IMF.Taking annual maximum peak flow series at Huayuankou Station in lower Yellow River from 1952 to 2009 as example，different oscillation components with different fluctuation periods and residual component of annual maximum peak flow series are obtained，and variation of each component is analyzed in detail.The result shows as follows that the trend of annual maximum peak flow series at Huayuankou Station is stepwise decreasing，there are four periodic oscillations of 3.2 a，6.4 a，11.8 a and 31.0 a for the process of annual maximum peak flow variation，and the variation of annual maximum peak flow is caused mainly by the periodic oscillations of 3.2 a and 6.4 a in fact.Then multiple time-scale structure of the variation of annual maximum peak flow series can be revealed.Based on the above，the reasons and infection factors for the variation of each fluctuation component are discussed.The variation of each IMF may be caused by QBO，ENSO，solar activity and climate change etc.

Hilbert-Huang Transform（HHT）;Empirical Mode Decomposition（EMD）;multiple time-scale;annual maximum peak flow

TV142

A

1672－0504（2012）03－0083－04

2011－11- 02;

2012－02－26

水利部公益性行業(yè)專項(xiàng)（200901021）

任?。?982－），男，博士研究生，研究方向?yàn)樗W(xué)及河流動(dòng)力學(xué)、水文泥沙與河床演變分析。E－mail：ren＿jian＠126.com