(1华中科技大学数字化工程与仿真中心,湖北 武汉,430074; 2华中科技大学水电与数字化工程学院,湖北 武汉,430074)
[摘 要]在mexh小波函数的基础上,采用了一种新颖、快速、实用的非正交小波变换方法,并给出了该方法的两种推导过程。文中对丹江口以上地区径流和降水量时间序列进行了小波分析,揭示了它们的多时间尺度结构,分析了不同时间尺度下时间序列变化的周期特征。根据小波系数拟合出周期变化曲线,并与实际数据进行对比。在实际数据受到随机成分影响的情况下,比较结果仍然相当吻合。
[关键词]小波变换;快速算法;周期分析;小波方差
水文时间序列的周期成分主要受地球绕太阳公转以及地球自转的影响而形成,同时也受气候、地理以及人类活动的影响。其中, 确定性成分是主要项,决定了周期性特征,而随机
[1]
成分是干扰项。水文时间序列的周期不是严格意义上的周期,只是概率意义上的周期,并且同一时段中可能包含多种时间尺度的周期变化,即水文时间序列的变化在时域中存在多层次时间尺度结构和局部性特征。水文过程周期规律的研究对探讨水循环的时空变化是很有意义的。水文序列分析的传统方法有滑动平均、滤波、傅氏分析等,它们的不足在于:①在时
[2]
域和频域上不具备局部化性质;②对突变点的诊断缺乏数学上的严谨性。20世纪80年代
[3]
发展起来的小波分析具有多分辨率分析的特点,为刻画非平稳信号提供了有力工具。小波变换表征信号局部特征的能力,可以剖析时间序列内部精细结构,有利于揭示水文序列的多时间尺度变化特征。
[4]
在利用小波计算水文时间序列周期方面,已经涌现出几种方法。刘晓安采用离散小散变换函数和db4小波对流量序列进行小波变换获得周期信息,但是此方法的理论依据不足,
[5]
而且目测图像信息的误差较大。赵利红采用离散小波变换公式和Morlet小波算得周期结论,
[6]
但王红瑞等证明由于Morlet小波不满足容许性条件,使得结论的精确度不如满足容许性条
[1]
件的小波。桑燕芳采用db4小波变换方法进行分解、去噪和重构获得周期信息,但结果显示出序列长度对周期识别存在一定影响, 而且合理确定随机成分对应的小波系数阈值是很
[7]
困难的。汤成友等采用Mallat算法和Daubechies 小波先对水文时间序列去噪,再对确定性成分进行方差分析从而算得周期,该方法实际上是以统计方法为主,小波分析为辅,分解尺度具有主观性,而且工作量较大。
1 mexh小波变换理论
1.1 小波变换介绍
For a given wavelet function aa, hydrological time series of the continuous wavelet transform ab to cd
对于给定的小波函数(t),水文时间序列x(t)L2R的连续小波变换为
收稿日期: 基金项目:“十一五”国家科技支撑计划(2006BAB04A07,2006BAB04A09) 作者简介:杨正祥(1976-),男,汉族,湖北汉川人,博士生,研究方向为计算智能算法与工程应用等。E-mail:yangzhengxiang999@yahoo.com.cn 1
Wx(a,b)a12tbx(t)adt (1) tb为(t)经a式中,a为尺度因子,反映频域特征;b为时间因子,反映时域特征;伸缩和平移后取共轭得到的一簇函数。Wx(a,b)为小波变换系数,是连续小波在尺度a、位移b上与信号的内积,表示信号与该点所代表的小波的相似程度。当x(t)与[4]
a,b(t)相等时,
该点的小波系数值为1。当a较小时,对频域的分辨率低,对时域的分辨率高;当a增大
[8]
时,对频域的分辨率高,对时域的分辨率低。正是在这个意义上小波变换被誉为“数学显微镜”。
Where, aa for the scaling factor, reflecting the characteristics of the frequency domain; bb for the time factor, reflecting the characteristics of time-domain; cc for by the expansion and post-translational conjugation obtained from a cluster of function. dd for the wavelet coefficients, continuous wavelet in the scale is aa, with aa displacement of the inner product signal, said signal with the point represented by the degree of similarity of the wavelet. When the ee and ff, etc., the point of the wavelet coefficient is 1. When gg smaller, the resolution of the low frequency domain of high-resolution time-domain; when jj increases, the resolution of the high frequency domain of low-resolution time-domain. It is in this sense, wavelet transform known as the \"mathematical microscope.\"
1.2 mexh小波系数的计算方法
本文采用的是mexh小波,其函数形式为:t(1t)e212t2。它是高斯函数的二阶导
数(加负号),可知在=0处有二阶零点,所以满足容许条件,而且其小波系数随衰减较快。mexh小波的时、频两域都有很好的局部性,并且满足tdt0。由于它不
R存在尺度函数,所以此小波函数不具有正交性。
This article is based on wavelet mexh, the function of the form: ff. It is a Gaussian function of the second derivative (plus minus), we can see in the dd = 0 Department ff there are second-order zero, so to meet the permit conditions, and the wavelet coefficients decay rapidly with the cc. mexh wavelet, the frequency of the two domains have very good locality, and to meet xx. Because of its scaling function does not exist, this wavelet function does not have orthogonality.
对于正交小波的快速计算,已经发展了相当完善的Mallat算法,而对于非正交的小波变换,文献资料鲜有报道。本文采用了一种快速、实用的算法,流程如图1所示。其中,
2
FT和IFT分别表示傅氏和逆傅氏变换。 For fast orthogonal wavelet basis, has developed a comprehensive Mallat algorithm, while for non-orthogonal wavelet transform, the literature has been reported rarely. In this paper, a fast and practical algorithm, the process shown in Figure 1. Which, FT and IFT, respectively, Fourier and inverse Fourier transform.
水文时间序列x(t) mexh小波:t(1t)e212t2 XFTxt a2e212a2 FTWxa,baXa Wxa,bIFTaXa
图1 mexh小波系数的计算流程图
以下从两个不同的角度----Parseval定理和时域卷积定理推导上述算法。
①从Parseval定理的角度推导:
由Parsevarl定理知对于两函数f1t和f2t有:
The following from two different angles ---- Parseval theorem and time-domain convolution theorem of the above algorithm is derived. ① from the perspective of the Parseval theorem is derived: Known theorem from the two Parsevarl function aa and bb are:
f1tf2tdt=
12F1F2d=
12F1F2d
由小波变换公式可以转换为: Wxa,b=xt,a,bt=
12XFTa,btd=
121tbXFTdaa=
a2Xaejbd=IFTaXa
②从时域卷积定理的角度推导:
小波变换公式变形为: ② from the time-domain point of view convolution theorem is derived: Deformation of wavelet transform formula as follows:
3
Wxa,b=
xt,ta,b=
1atbxtdt=
a1abtxtdt=
a1txt
aa当a看作常数而b看作连续时间变量时,上式为xt和
1at的卷积。由时域a卷积定理即时域中两函数的卷积等效于频域中频谱的乘积得:
When the cc as a constant and continuous-time as a variable, for the aa and ss-type of convolution. Theorem by the time-domain convolution function of two real-time domain convolution is equivalent to the frequency domain spectrum of the product was:
FTWxa,bFTxt1taa1tFTxtFTaa===
Xt1aaa=aXta
可见,从Parsevarl定理的角度和时域卷积的角度推导出的小波变换系数公式是一致的。
由此得非正交小波变换系数为:
Can be seen from the perspective of Parsevarl theorem and time-domain convolution of the point of view derived from the formula of the wavelet transform coefficients are the same. This was non-orthogonal wavelet transform coefficient: aa
Wxa,bIFTaXta (2)
上述算法求解非正交小波变换系数相比常规的离散小波变换算法具有如下明显的优点:
①将小波变换系数的计算转换到纯频域中进行,只需计算时间序列和小波函数的傅氏变换即可,可以充分利用FFT(快速傅氏变换)算法优势,运算速度较快;
②由于式(2)是在时间序列xt和小波函数t的整个时间段上进行,因此小波变换系数中的变量b自动映射到整个时间段上而不用专门设置,所以计算结果只是a的函数,只需
[9]
对a离散化即可,计算过程较为简单,并且易于计算机程序实现。
Algorithm for solving the above-mentioned non-orthogonal wavelet transform coefficients compared to conventional discrete wavelet transform algorithm has obvious advantages as follows:
① to the calculation of wavelet coefficients to switch to pure frequency domain, simply calculating the time series and wavelet function to the Fourier transform, can take full advantage of FFT (Fast
4
Fourier Transform) algorithm advantage of faster;
② as a result of type (2) ss in the time series and wavelet function of the cc on the entire time period, the wavelet coefficients of the variables aa automatically mapped to the entire period of time rather than specialized settings, the results only a function of xx, cc discretization only for you to calculate the process more simple and easy realization of a computer program.
2 应用实例
丹江口水库是汉江的南水北调中线工程水源地。汉江作为长江的一条支流,既要提供本流域的“三生”用水,又要承担向河南、河北、北京、天津供水的重任,具有较大挑战性。本文运用mexh小波分析方法,对丹江口地区的径流和降水的震荡变化规律进行研究,识别周期成分,揭示它在不同时间尺度下的波动特性,为汉江流域水资源开发利用和南水北调中线工程水资源优化配置提供技术支撑。本文所采用的数据资料是1933-2001年的径流量(来自长江水利委员会)和1961-2006年的降水量(来自中国水利水电科学研究院)。 2.1径流周期性分析 2.1.1 初始数据处理
为减少径流序列中季节变化及短期噪声的干扰,首先对年平均流量序列进行标准化处理:
In order to reduce seasonal variation in runoff series and short-term noise, first of all, the average annual flow of the standardization of sequence processing:
RiRiR/ (3) 式中:Ri为第i年的标准径流量;Ri为第i年实测平均径流量;R、分别为多年平均径流量和均方差。标准化处理后,实测序列即转换为标准化序列。
Where: aa for the first years of the standard cc runoff; zz measured for the first year the average zz runoff; xx, cc, respectively, the average runoff for many years and the mean square deviation. Standardization of treatment, measured sequence that is converted to standardized sequence.
丹江口以上地区1933~2001年(69年)的标准径流序列如图7或图8中的无点实线所示。 2.1.2 小波变换
用mexh小波对标准径流序列进行连续小波变换,由计算结果绘制小波变换系数等值线
[11]
图(图2),图中反映出降水量序列的周期变化、突变点分布和位相结构特征。小波系数为正时,表示径流量偏多,图中用实线绘出;为负时表示径流量偏少,图中用虚线绘出。由图3可见,能量中心的频域尺度主要集中在8~10年,22~24年(即在这两个区域内有规
5
[10]
''律地出现丰枯交替过程),代表了较为明显的两个主要周期。
Mexh wavelet with the standard run-off continuous wavelet transform sequences, drawn from the calculation of wavelet transform coefficient contour map (Figure 2), the figure reflects the sequence of the cycle of precipitation change, mutation distribution and phase structural characteristics. Wavelet coefficients for the time being, said runoff on the high side, with solid line drawn map; negative that low runoff, map drawn with dotted lines. Can be seen from Figure 3, the energy center of the frequency-domain measure mainly concentrated in the 8 to 10 years, 22 ~ 24 years (that is, in these two regions appear regularly alternating丰枯process), on behalf of the more obvious the two main cycles .
图2 径流的小波变换系数等值线图
2.1.3 小波方差
小波方差公式为:
Wavelet variance formula is:
Vara它的离散形式为:
Its discrete form as follows: Vara1[12]
Wxa,bdb (4)
2Wa,x (5) n2tt1n式中:W2a,xt为尺度a、时间xt处的小波系数的平方;n为系列的长度。
Where: aa to scale cc, time xx Department of the square of the wavelet coefficients; n is the length of series.
6
根据式(5)计算小波方差,并作小波方差图(图3),该图能反映波动的能量随尺度a[9]
的分布,可以用来辩识时间序列中各种尺度扰动的相对强度和周期特征。小波方差随时间尺度变化的过程中有3个峰值,横坐标分别是3年、9年、23年。其中后两个峰值在相应尺度下信号震荡强烈,所以9年和23年是主要周期,3年则是次要周期。18年处也对应一个峰值,但此峰值不明显,可视为小波方差值的小波动,予以忽略。
According to type (5) calculation of wavelet variance, and wavelet variance map (Figure 3), the plan to reflect the fluctuations in the distribution of energy with the scale, time series can be used to identify various scales and the relative intensity of the disturbance cycle特征. Wavelet variance changes with the time scale of the process of three peaks, respectively, the abscissa is 3 years, 9 years, 23 years. The latter two of which scale the peak signal in the corresponding strong shocks, so 9 years and 23 years is a major cycle, 3-year cycle is of secondary importance. 18 years also correspond to a peak, but this peak was not obvious, can be regarded as wavelet small fluctuations in the margin to be ignored.
图3 径流的小波方差图
2.2降水量周期性分析 2.2.1 初始数据处理
丹江口以上地区1961~2006年(46年)的标准降水量序列如图9中的无点实线所示,对实测数据的处理方法同上。 2.2.2 小波变换
对上述标准降水量序列作小波变换,结果见图4。图5显示,从1961-2006年整个时间序列都可见,能量中心的频域尺度主要集中在22~24年,代表了较为明显的主要周期。在1961-1990年间,在频域尺度8-10年内也有规律地丰枯变化,但后面没有延续这种规律,这与实际情况相符。进入90年代以来,实际降水偏枯。
7
图4 降水量的小波变换系数等值线图
2.2.3小波方差 进一步计算出小波方差,并作小波方差图(见图5)。小波方差随时间尺度变化的过程中有3个峰值,横坐标分别是3年、9年、23年。其中最后一个峰值在相应尺度下信号震荡强烈,所以23年是主要周期,3年和9年则是次要周期。5年处也对应一个峰值,但此峰值不明显,可视为小波方差值的小波动,予以忽略。
图5 降水量的小波方差图
2.3结果分析
选取标准径流序列的两个尺度(9年和23年),分别绘制小波系数曲线,如图6和图7中的带点实线所示,可得它们的周期分别为9年和23年。图中的小波系数曲线实际上是对实际数据进行了一定程度“去噪”、“拟合”的结果,即保留低频成分,滤掉高频成分。其震幅表明水文信号的强度。
Select the standard two-scale runoff series (9 and 23 years), respectively, the wavelet coefficients drawn curve in Figure 6 and Figure 7 as shown in solid line a bit to get their cycle of 9 years and 23 years. Graph curve of the wavelet coefficients is the actual data to a certain extent \"de-noising\\"fit\" the results, that is to retain low-frequency components, high frequency components filtered out. Its amplitude shows that the intensity of the hydrological signal.
8
图6 9年尺度小波系数与标准径流序列对比图
图7 23年尺度小波系数与标准径流序列对比图
选取降水量的23年尺度,绘制该尺度的小波系数曲线,如图8中的带点实线所示,其周期为23年。
图8 23年尺度小波系数与标准降水量对比图
可见,在实际数据受到随机成分影响的情况下,实际曲线与拟合曲线仍然相当吻合。拟合曲线都呈现周期性,但同一拟合曲线各周期间的震幅并不严格相同,恰好说明小波分析在注重整个时间序列的频域信息的同时,也充分考虑了局部的时域信息。
Can be seen in the actual data by the random element of the impact of cases, the actual curve and fitting curve are still very fit. Were the only cyclical curve, but the same fitting curve of the amplitude of the Week is not strictly the same, just that the whole emphasis on wavelet analysis of time series in the frequency domain information, but also take full account of the partial time-domain information.
由前面分析知,对于1933~2001年的径流,9年和23年是主要周期,3年则是次要周期;对于1961~2006年的降水量,23年是主要周期,3年和9年则是次要周期。结果表明降水和径流具有明显相关性。实际数据显示进入90年代以来,降水量和径流都有所减少,9年周期强度减弱,表明人类活动的影响加大。因为1933~2001年径流序列较长,能较大程度地“屏蔽”90年代以来的流量序列,所以结果仍呈现较强的9年周期强度。而1961~2006年的降水量序
9
列则相反。
3 结束语
实际上,历史水文记录,都只是随机来水的样本。当抽样系列相当长时,它的统计特性和随机过程的特性是一致的。这样,根据历史资料所得出的结论便可以作为实际工作的依据。水资源优化配置是一个具有多目标的复杂系统工程,深入研究和把握其中存在的定性规律是十分必要的,即便这些规律只能在一定条件下成立且只能反映这一复杂系统工程的某些侧面。设法使这些侧面问题得到较好解决,将有利于整个系统的解得到逐步改善,最终达到较
[13]
为满意的境界。
[参 考 文 献]
[1]桑燕芳,王栋.水文时间序列周期识别的新思路和两种新方法[J].水科学进展,2008,19(3):412-417. [2]吕翠美,吴泽宁.刘文立,等.伊河流域径流周期变化特征的小波分析[J].人民黄河,2007,29(5):26-28. [3]邓跃红,聂双双.小波分析在水下回波奇异性检测中的应用[J].计算机工程与科学,2008,30 (3):123-125. [4]刘晓安.小波分析在径流分析和预报中的应用研究[D].武汉:华中科技大学,2006. [5]赵利红.水文时间序列周期分析方法的研究[D].南京:河海大学,2007.
[6]王红瑞,叶乐天,刘昌明,等.水文序列小波周期分析中存在的问题及改进方式[J].2006,16(8):1002-1008. [7]汤成友,缈韧.基于小波变换的水文时间序列分解及周期识别[J].人民长江,2006,37(12):32-34. [8]Torrence C, Compo G P. A practical guide to wavelet analysis[J]. Bull. Amer. Meteor. Soc. ,1998, 79:61-78. [9]罗光坤. Morlet小波变换理论与应用研究及软件实现[D].南京:南京航空航天大学,2007.
[10]NAKKEN M. Wavelet analysis of rainfall-runoff variability isolating climatic from anthropogenic patterns[J]. Environmental Modelling&Software,1999,14:283-295.
[11]刘俊萍,田峰巍,黄强,等.基于小波分析的黄河河川径流变化规律研究[J].自然科学进展,2003,13(4):383-387.
[12]Cahill A T.Determination of changes in stream flow varicance by means of a wavelet-based test[J]. Water Resources Research,2002,38(6):1-14.
[13]张勇传.系统辩识及其在水电能源中的应用[M].武汉:湖北科技出版社,2008 ,112-155.
Application of a New Wavelet Algorithm in the Hydrological Periodical Analysis
YANG Zheng-xiang1, HE Xiao-Cong1
(Digital Engineering and Simulation Research Center, HUST , WuHan HuBei,430074, China)
(School of Hydropower and Information Engineering, HUST ,WuHan HuBei,430074, China)
Abstract::On the basis of mexh wavelet function, a novel, rapid and practical method
of non-orthogonal wavelet transform is used, and two types of derivation for the method are given. In this paper, the time series of stream-flow and precipitation of Danjiangkou region above are analyzed by using wavelet transformation. Their structures of multiple-scale time series are revealed and the periodicity at different scales are analyzed. The Periodical curves simulated according to wavelet coefficients are in good agreement with the actual data affected by the random component.
Keywords: wavelet transformation; rapid algorithm; periodical analysis; time scale
10
11
因篇幅问题不能全部显示,请点此查看更多更全内容