任 健，史 红 玲

花园口站年最大洪峰流量时间序列的HHT分析

任 健1，史 红 玲2

（1.中国水利水电科学研究院，北京 100048;2.国际泥沙研究培训中心，北京 100048）

Hilbert-Huang变换能够定量描述非线性、非平稳复杂时间序列的时频特性，较传统分析方法更具优势。通过对时间序列进行EMD分解，得到变化过程的内在模态函数和趋势项函数，而后对各内在模态函数进行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变换;经验模态分解;多时间尺度;年最大洪峰流量

黄河花园口站是黄河下游重要的控制性监测断面，其控制流域面积达73万km2，约占黄河流域面积的97%，分析花园口站的水沙变化情势是研究黄河下游河道演变的前提和基础。长期以来，关于黄河下游水沙变化特征的研究很多［1－4］，但多依据传统的数理统计或简单的数据平滑等方法，不能揭示水沙变化的多时间尺度特性［5］。而Fourier分析和小波分析等方法虽然可以进行频谱分析，但其方法本身存在一些不足和缺点［6，7］。希尔伯特—黄变换［6，7］（Hilbert-Huang Transform，HHT）技术是近年来发展的一种最新的处理非线性、非平稳时间序列的方法，它直接由时间序列本身自适应地构造出基函数，进而分离得到不同尺度的特征分量，而不用预先设定基底函数再展开，因而该法可更准确反映序列的内在规律特征。本文运用H HT技术，对黄河下游花园口站1952－2009年的实测年最大洪峰流量序列进行多时间尺度分析，诊断花园口站年最大洪峰变化过程中蕴含的多时间尺度振荡结构和特征，揭示其在不同尺度和层次上的变化规律，以期为黄河下游的水资源利用、河床演变、河道整治及防洪调度等提供技术支撑和科学依据。

1 研究方法和数据资料

经验模态分解（Empirical Mode Decomposition，EMD）方法［6，7］是一种全新的处理非线性、非平稳数据序列的方法，它将时间信号f（t）分解成一系列本征模态函数（Intrinsic Mode Function，IMF），每个IMF分量都具有如下特征：从全局特性看，极值点数必须和过0点数一致或至多相差一个;在某个局部点，极大值包络和极小值包络在该点的值的算术平均和是0。

EMD分解过程为：找出序列f（t）所有的极大值和极小值点，分别用三次样条函数拟合成上下包络线，得到平均包络线m1，将原序列减去m1得到去掉低频序列的新序列h1。一般h1并不能满足IMF的两个特征要求，仍是非平稳的，因而可以多次重复上述过程，使平均包络线趋近于0，得到第一个IMF分量c1，c1代表原始序列中最高频的分量。即：

对r1（t）继续进行上述分解，直到剩余部分为单一信号或其值小于预先给定的值，分解结束。原始时间序列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分量在某个时刻出现波动的频率，与该时刻的波动周期互为倒数关系，因而其物理意义很明晰。由此可知，中心频率即为某IMF分量各点瞬时频率的平均值，其倒数即为平均周期。EMD方法不可避免地存在边界效应，本文采用镜像对称延伸方法对边界效应进行了处理［8］。

采用黄河下游花园口水文站1952－2009年的年最大洪峰资料（来源于黄河下游河床演变基本资料汇编及黄河泥沙公报）获得近60年来花园口站年最大洪峰流量的变化过程（图1）。通过图1只能直观判断近60年花园口站年最大洪峰流量变化的趋势，并不能深入、细致地揭示最大洪峰流量的演变周期及多时间尺度特征，因此有必要运用基于EMD的HHT技术对其进行多时间尺度分析，从有限的数据序列中挖掘尽量多的有用信息，以便更深刻认识洪水变化的内在规律和特性。

图1 花园口站年最大洪峰流量时间序列Fig.1 Annual maximum peak flow series at Huayuankou Station

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

运用EMD方法对花园口站1952－2009年的年最大洪峰流量时间序列进行多时间尺度分解，得到4个本征模态函数（c1～c4分量）和一个残余趋势项（Res分量）（图2）。

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

由图2可知：1）花园口站实测年最大洪峰流量变化过程是非线性和非平稳的，是多种波动成分共同作用的结果，其可分解为4个具有不同波动周期的振荡分量和1个趋势分量，反映出花园口站年最大洪峰流量变化具有复杂的多时间尺度特性。2）c1分量具有准2～4 a的波动周期，且不同时期的振幅也不同，呈现出明显的阶段性特征。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波动周期，并且不同时期的振幅也不同，呈现出明显的阶段性特征。其中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分量的方差贡献并统计其Hilbert变换后的特征值（表1）。从图3可以看出，各IMF分量的频率并不是一个恒定值，而是围绕中心频率上下波动，频率越高的IMF分量，其波动的幅度越大。尽管IMF分量的频率围绕中心频率波动，但波动的范围通常有限，相互之间少有交叉重叠现象，保持了一种界限较为清晰的分布特征。由时频图和时幅图也可以看出，周期越长，在该周期上瞬时频率的波动就越小，振幅的变化也越小，即对应于大尺度周期，年最大洪峰流量的变化相对稳定，但波动能量较小。

图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变换后的数字特征值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，总体上表现出中心频率由高到低的变化特点，而平均周期则由短变长。振幅反映了各分量波动能量的大小。分析各个分量的平均振幅和最大振幅，中短周期波动的振幅一般较大，而长周期波动的振幅一般较小，呈现出频率越高（或周期越短）振幅越大的特点。方差贡献代表了各分量对原序列信号波动趋势的影响程度，方差贡献越大，说明该分量的波动对原时序变化的影响程度越强。不难发现，中短周期波动的方差贡献率远大于长周期波动，这表明中短周期的波动分量是引起年最大洪峰流量变化的主要原因，而长周期波动分量对年最大洪峰流量变化过程的影响主要体现在对其整体发展趋势的控制性作用。

4 讨论

黄河下游大的洪峰流量多发生于伏汛、秋汛期间，是由暴雨引起，因而年最大洪峰流量是由年极端降水事件造成的，二者具有天然的、复杂的联系，年极端降水事件的波动影响着年最大洪峰流量的变化过程。短时期极端降水事件的发生受到大气运动、气象等因素的影响，但长时期极端降水事件发生的频率（周期）与强度受到全球及区域的物理、气候系统变化的控制与影响。已有的研究表明平流层的大气运动存在准两年周期的振荡（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周期基本一致。因而，花园口年最大洪峰流量变化过程中蕴含不同周期的波动特征，一定程度上反映了大气、气象及气候等因素对黄河流域水文要素可能存在的影响作用。

随着黄河流域内大型水库调度运营和人类活动影响程度的加剧，花园口站年最大洪峰流量各波动分量的周期和振幅也发生不同程度的调整，如图1中各分量波动曲线及图3中时幅、时频谱图所示。对于中短周期分量（c1～c3分量）而言，人类活动的影响导致其周期在局部时段内发生变化，同时振幅也发生阶段性的变化;并且各分量随着其中心频率的递减，周期和振幅的调整变化也趋缓。而对于长周期分量（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、太阳活动及气候变迁等因素有关，加之人类活动的严重影响，使得各波动分量的变化规律更加复杂。而对于包括人类活动影响在内的各种因素对黄河流域年最大洪峰流量变化过程的影响机制有待深入研究。

此外，HHT技术对于受干扰因素多、变化规律复杂、随机成分大的非线性、非平稳的水文序列具有较好的分析和诊断效果，在水文过程的多时间尺度分析、建模和预测中值得应用和推广。

［1］陈浩，周金星，陆中臣，等.黄河中游流域环境要素对水沙变异的影响［J］.地理研究，2002，21（2）：179－187.

［2］刘昌明，郑红星.黄河流域水循环要素变化趋势分析［J］.自然资源学报，2003，18（2）：129－135.

［3］刘勇胜，陈沈良，陈小英，等.黄河中下游泥沙通量变化规律［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］杨建平，丁永建，陈仁升.长江黄河源区水文和气象序列周期变化分析［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

水利部公益性行业专项（200901021）

任健（1982－），男，博士研究生，研究方向为水力学及河流动力学、水文泥沙与河床演变分析。E－mail：ren＿jian＠126.com