任 健，史 紅 玲

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

任 健1，史 紅 玲2

（1.中國水利水電科學研究院，北京 100048;2.國際泥沙研究培訓中心，北京 100048）

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

Hilbert-Huang變換;經驗模態(tài)分解;多時間尺度;年最大洪峰流量

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

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

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

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

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

式中：cj（t）為各IMF分量;Res為趨勢項;t為時間變量。

對每個IMF分量cj（t）進行Hilbert變換，即：

p為Cauchy主值，cj（t）和bj（t）可以構成一個復序列：

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

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

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

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

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

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

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

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

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

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

圖3 基于HHT分析的花園口站年最大洪峰流量序列時幅和時頻譜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可知，對于花園口站年最大洪峰流量序列而言，從第1個分量到第4個分量，其中心頻率分別為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)出中心頻率由高到低的變化特點，而平均周期則由短變長。振幅反映了各分量波動能量的大小。分析各個分量的平均振幅和最大振幅，中短周期波動的振幅一般較大，而長周期波動的振幅一般較小，呈現(xiàn)出頻率越高（或周期越短）振幅越大的特點。方差貢獻代表了各分量對原序列信號波動趨勢的影響程度，方差貢獻越大，說明該分量的波動對原時序變化的影響程度越強。不難發(fā)現(xiàn)，中短周期波動的方差貢獻率遠大于長周期波動，這表明中短周期的波動分量是引起年最大洪峰流量變化的主要原因，而長周期波動分量對年最大洪峰流量變化過程的影響主要體現(xiàn)在對其整體發(fā)展趨勢的控制性作用。

4 討論

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

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

5 結論

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

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

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

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

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

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

［5］張少文，丁晶，廖杰，等.基于小波變換的黃河上游天然徑流變化特性分析［J］.四川大學學報（工程科學版），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］彭永清，王盤興，吳洪寶.大氣低頻變化的分析與應用［M］.北京：氣象出版社，1997.

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

［11］孫林海，趙振國.我國暖冬氣候及其成因分析［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è)專項（200901021）

任?。?982－），男，博士研究生，研究方向為水力學及河流動力學、水文泥沙與河床演變分析。E－mail：ren＿jian＠126.com