非平稳时间序列

非平稳时间序列（精选7篇）

非平稳时间序列 第1篇

时间序列分析有时域和频域分析两类方法, 由于频域法的计算复杂且结果不够直观。小波分析是以多分辨率分析是一种时/频综合描述方法, 相对于加窗傅立叶变换, 小波基函数具有适应频率变化的可变窗宽, 集成了时域和频域分析方法的理论精华, 尤其对非平稳信号中短时高频现象和低频轮廓信息具有很好的“显微”和“概括”效果。

1 小波变换

对于能量有限信号f (t) ∈L2 (R) , 有

其中Wf为小波系数, ψa, b为积分为零且平方可积的母小波函数。a, b是将连续小波函数中相应参数离散化的结果, 其中a为尺度参数, b为平移参数或时间参数。随a, b变换, ψa, b产生一系列小波基函数, 对ψi用不同尺度的滤波器进行滤波。相对于加窗傅立叶变换, ψa, b具有适应频率变化的可变窗宽, 即, 高频时ψa, b的时窗变窄, 低频时则时窗变宽。傅立叶分析是将信号分解为一系列不同频率的正弦波的加权叠加, 而小波分析是将其分解为一系列小波基函数的加权叠加, 这些小波基函数由一个母函数经过平移和尺度两个指标的协同伸缩而获得的, 其权值序列即为小波系数和逼近系数序列, 表征了信号与该尺度下的小波函数的相关程度。

2 小波分析在时间序列分析中的应用

小波变换将时间序列按不同频率特征分解为逼近系数序列和细节系数序列, 事实上是将信号的某些特定信息用新的方法呈现出来, 信号和它的小波变换是同一数学实体的两个不同表现。传统时间序列统计分析的中心目的是寻找将相关序列表示成不相关序列的组合的方法。通过对小波变换得到的细节系数序列进行自相关检验, 发现其相关性已大大降低, 即小波分析能有效地进行去相关处理。通过快速算法, 使得分解时间序列的过程更为简便。

与传统方法相比较, 逼近系数序列与传统方法中趋势项和季节项类似。传统方法剔除趋势后, 再根据固定周期内序列的方差分析, 提取出一个固定季节分量, 则认为余下残量即为随机分量。而小波分析是将第一次低通滤波处理后的逼近系数序列进一步进行多次不同频宽的低通滤波处理, 从而得到多组不同频率特征的细节分量。因此, 小波分解得到的逼近系数序列不仅表现出固定的季节成分, 同时将每个季节周期内低频段的其他特征也表现出来, 这比传统方法更为细致。

3 乘性时间序列的小波方差分析

小波变换在时间序列分析方面的一个重要应用就是在不同尺度上对样本方差进行分析。小波方差就是时间序列分解后各层小波系数的方差。

对于乘性季节分量时间序列, 小波分解得到的小波系数序列仍随季节变化, 这样的结果不满足方差齐性的要求。其原因可由小波分析的快速算法可知, 小波系数序列也可以表示为上一层逼近系数序列与当前层逼近系数序列之差, 如果逼近系数序列代表季节分量, 则小波系数代表的随机分量, 两者是加性关系。因此小波系数不能准确地描述乘性序列的特性以及小波系数方差。

对于这类乘性数据, 本文提出使用当前层的逼近系数序列视作乘性季节分量, 作为乘积因子对小波系数序列作方差齐性处理, 即以逼近系数序列对应的点为分母, 将小波系数序列的季节性成分剔除, 则得到原序列的成分模型为：

其中, 函数表示求序列的均值。

4 试验结果

由小波分解的各层系数序列图可以看出, 细节系数仍随季节变化。日均水量图d1层小波系数序列作去乘性季节分量的处理, 得到图1。

处理后的系数序列相关图 (图2) 表明, d1层周期性基本消除。偏相关的截尾阶数也降低很多。故使用逼近系数序列对小波系数序列进行去乘性季节分量具有一定的效果。

结束语

本文提出利用尺度系数作为季节性乘积分量, 将小波系数序列处理为方差平稳的平稳序列, 以验证原序列的乘性特征, 并将此方法用于样本数据异常点的识别。实验证明, 使用修正了异常点的样本数据所建立的预测模型, 在预测精度上有所提升。

摘要：本文提出了基于小波变换的非平稳时间序列分解方法。通过小波分解, 将原时间序列依尺度分解成不同层次, 利用尺度系数作为季节性乘积分量, 将小波系数序列处理为方差平稳的平稳序列, 并应用于样本数据异常点的识别和修正, 使模型的预测精度得到提升。实例验证该方法是有效的。

关键词：非平稳时间序列,小波变换,小波方差

参考文献

[1]赵黎明, 吴文清, 刘嘉焜.基于小波分析的游客流量神经网络预测研究[J].系统工程学报, 2006.

[2]郭其一, 路向阳, 李维刚, 王钰.基于小波分析和模糊神经网络的水文预测[J].同济大学学报 (自然科学版) , 2005 (1) :130-134.

非平稳时间序列 第2篇

非平稳环境下混沌信号的小波去噪方法

针对传统的去噪方法对混沌信号除噪的盲目性,及往往仅对平稳环境或缓慢变化的噪声有效的局限性,本文提出根据信号与噪声在小波域的分布特性及信号和噪声的.模极大在细尺度下收敛的横坐标点来检测信号的奇异性,以分形维树理论为依据决定阈值,得到噪声在小波域中的位置以及小波系数大小实现去噪.实验结果证明此法有效可行.

作 者：赵颖 孙鹏勇 ZHAO Ying SUN Peng-yong 作者单位：辽宁工程技术大学通信与信息系,辽宁,阜新,123000刊 名：电光与控制 ISTIC PKU英文刊名：ELECTRONICS OPTICS & CONTROL年，卷(期)：13(6)分类号：V271.4 TN911.4关键词：非平稳环境 混沌信号 离散小波变换 分形维数

非平稳条件下的图像去噪 第3篇

1 信号噪声的小波分析

在实践中,可以遵循以下3个步骤:

(1)信号的小波变换。选择小波和德组合数M,然后做M小波构成并得到水平系数记录wj,k。

(2)高频系数的量化。对于每一个级别,选择一个阈值并且得到估计系数记录

(3)信号的重新配置。此方法是基于M水平的低频系数和每个级别上高频系数的估计。

在这3个步骤中,关键的一步是如何选择一个合适的阈值,这关系到去噪效果的好坏。

2 阈值的选择

有许多选择阈值的方法。一个较简单的方法就是使用固定的λ并使它始终保持不变。通常它可以计算如下:注意:σ是比例系数;Ν是采样点的数目。根据以上规则可以得到如下算法:

称为硬阈值算法。

称为软阈值算法。

前面提到的硬或软的阈值方法在国内外已经用于实践中,然而效果却不尽如人意。因此此方法在实际应用中受到很大的局限性。对此本文引入了一种基于梯度的选择阈值的新的改进算法。

3 基于梯度的阈值选择

目前原始信号可以用来估计均方误差。本文引入了一个基于梯度的新的改进法方法,它可以实现自适应的均方误差估计以及选择合适的阈值。

定义为:且g=[g0,g1,...,gN-1],如果g(y)是可微的,则均方误差表示为:

且:

式中y的无偏估计可以定义为:

式中t是阈值匹配式(4)和式(5)可以发现式(5)是无偏估计的均方误差,他的梯度可以表

且:

式(3)中的阈值函数η(y,t)是微分无穷大。因此可以得到t的步长:

α表示学习效率,应该找一个合适的只来保证其实现:。则可以得到:

由于在此算法中必须给阈值赋初始值,选择:

式中:σlev是噪声水平估计;N是信号长度。如果满足Δtktk<ε→0此算法停止运行。此算法的流程如图1所示。

这一新的阈值选取方法在图像去噪中取得了很好的效果。该算法的优点是视线里在最小均方误差下获取阈值。

4 结果分析与结论

选择一个不平稳的方波信号并使用两种不同的方法来加以处理,通过Matlab仿真得到如下结果。图2(a)是原始信号;(b)是含有SNR=3的噪声信号;(c)是使用新方法后的效果图;(d)是使用固定阈值方法所得效果图。

通过信噪比结果,可以看出小波方法比传统傅里叶方法在非平稳信号去噪方面有更好的效果。新的自适应方法也要比之前的固定阈值法要好,因为它可以获得和原始图像更相似的信号。由于小波变换可以提供信号的频率及时间域信息,因此它被广泛的应用在非平稳信号的处理中。

摘要：介绍一种基于小波变换的图像去噪方法。这一新的方法是基于阈值的小波变换。并且此方法更适用于非平稳信号下的图像去噪。该算法根据图像与噪声在小波域的分布特性以及它们的小波变换模极大值随尺度的变化大小不同,运用迭代算法得到不同尺度小波域中噪声的具体位置以及小波系数大小,完成了图像去噪。实验结果表明,对峰值信噪比较低的图像,该方法去噪后峰值信噪比比传统方法的高,并且保留了较多的图像细节,同时对平稳和非平稳的噪声都能进行较好地去噪。

关键词：图像去噪,小波变换,非平稳性噪声去噪,阈值选择

参考文献

[1]BOLL S.Suppression of acoustic noise in speech using spectralsubtraction[J].IEEE Transactions on Acoustics,Speech andSignal Processing,1979,4(2):113-120.

[2]何炜.基于Daubechies小波的水下图像去噪方法研究[D].哈尔滨:哈尔滨工程大学,2006.

[3]LOCKWOOD P,BOUDY J.Experiments with a nonlinear spec tral subtractor(NNS),hidden Markov models and projectionfor robust recognition in cars[J].Speech Communication,1992,11(2/3):215-228.

[4]HYVRINEN A.Survey on independent component analysis[J].Neural Computing Surveys,1999,2(4):94-128.

[5]COHEN I,BERDUGO B.Speech enhancement for non stationarynoise environments[J].Image Processing,2001,81(11):2403-2418.

[6]SIM B L,TONGY C,TAN C T,et al.A parametric formulationof the generalized spectral subtraction method[J].IEEETransactions on Speech and Audio Processing,1998,6(4):328-337.

[7]ZHANG Z,ZHA H.Principal manifolds and nonlinear dimen sionality reduction via tangent space alignment[J].SIAM Jour nal of Scientific Computing,2005,6(1):313-338.

非平稳时间序列 第4篇

1 时频方法简介

设s (t) 是一实的非平稳信号, 在进行时频分析之前, 需要先将实信号s (t) 转变为复信号z (t) 的形式[3], 常用的方法是利用Hilbert变换定义实信号s (t) 的复信号。这种转换的定义式为:

z (t) =s (t) +jH[s (t) ] (1)

式中, H[s (t) ]是s (t) 的Hilbert变换, z (t) 称为s (t) 的解析信号。解析信号的优点是它剔除了实信号中的负频率成分, 同时不会造成任何信息损失, 也不会带来虚假信息。

1.1 短时傅里叶变换

短时傅里叶分析方法是Gabor于1946年提出的[4], 其基本思想是:在信号傅里叶变换前乘上一个时间有限的窗函数, 并假定非平稳信号在分析窗内的短时间隔内是平稳的, 通过在时间轴上移动窗从而使信号逐段进入被分析状态, 这样就可以得到信号的一组“局部”频谱, 从不同时刻“局部”频谱的差异上, 便可以得到信号的时变特性。

给定一个时间宽度很短的窗函数η (t) , 让窗滑动, 则信号z (t) 的短时傅里叶变换 (STFT) 定义为:

STFTz (t, f) =∫∞-∞z (τ) η (τ-t) e-j2πfτdτ (2)

STFT方法分析的信号必须具有分段平稳的特性。对于时间尺度不同的平稳小段, 则需要和自己相适应的窗的长度和类型, 才能获得最佳结果。然而由STFT的定义可知, 只能选定一个固定的窗函数, 即时频分辨率是固定不变的, 这是短时傅里叶分析方法的主要缺陷。

1.2 Wigner-Ville分布

Wigner-Ville分布可被看作信号能量在时域和频域中的分布, 它属于Cohen类时频分布[5]。信号s (t) 的Wigner-Ville分布定义为:

undefined

式中, z (t) 是s (t) 的解析信号, *表示取共轭。

虽然Wigner-Ville分布具有好的时频聚集性, 但是对于多分量信号, 根据卷积定理, 其Wigner-Ville分布会出现交叉项, 产生“虚假信号”, 这是应用中存在的主要缺陷。交叉项是二次型时频分布的固有结果, 它来自于多分量信号中不同信号分量之间的交叉作用, 而交叉项的抑制又主要通过核函数的设计来实现。常用的加核函数后的Wigner-Ville分布有伪Wigner-Ville分布 (PWD) :

undefined

式中, h (τ) 是窗函数。

为了提高信号分量的时频聚集性, 同时减小交叉项, 一种有效的方法是对信号进行重排, 下面介绍Cohen类时频分布的重排。信号s (t) 的Cohen类时频分布可以看成Wigner-Ville分布的二维卷积, 即:

Cz (t, f;∏) =∫∞-∞∫∞-∞∏ (s-t, ζ-f)

Wz (s, ζ) dsdζ (5)

式中,

∏ (t, f) =∫∞-∞∫∞-∞Φ (ζ, τ) e-j2π (fτ+ζt) dτdζ (6)

Φ (τ, v) 为核函数。∏可以看成是一个平滑核函数, Cohen类时频分布可以理解为平滑后的Wigner-Ville分布, 当Φ (τ, υ) =1, 就是Wigner-Ville分布。核函数∏定义了Cohen类时频分布的重排方法, 即:

undefined

则, 重排后的Cohen类时频分布为:

Cundefined (t′, f′;∏) =∫∞-∞∫∞-∞Cz (t, f;

undefined

2 合成模拟的非平稳信号

为了证明二次型时频表示分析非平稳信号的有效性, 并比较STFT和时频表示分析非平稳信号的效果, 模拟生成一非平稳信号 (未加噪) 。该信号包含:两个频率分别为1000Hz和555 Hz的高斯调幅信号, 一个抛物线调频信号、一个正弦调频信号和一个线性调频信号:

undefined

式中,

undefined

采用频率设置为10kHz, 采样点数为N=512, 信号s (n) 和其相应的频率波形如图1所示。

3 仿真结果与讨论

利用本文第二部分所述内容, 求得合成信号s (n) 的时频分布如图2所示。其中STFT使用的分析窗函数是Hamming窗, 仿真时发现, 当窗的点数较小时, 信号的前三个高斯调幅信号的时频分辨率较好, 但后两个调频信号的时频分辨率较差;随着窗点数的增大, 调幅信号的分辨率降低, 而调频信号的分辨率增大。因此综合考虑整个信号段的时间分辨率和频率分辨率, 窗点数取51点较好。从图2 (b) 可以看到, 信号的WV分布存在比较严重的交叉项问题, 抛开交叉项问题可以看到, WV分布的分辨率较STFT要好。为了抑制交叉项, 本文求出信号的伪WV分布如图2 (c) 所示, 同时还计算量信号伪WV分布的重排如图2 (d) 所示。在计算信号伪WV分布的重排时, 分别选用了29点和53点Kaiser窗作为时间平滑窗和频率平滑窗。

另外, 因为实际信号都含有噪声, 所以本文也仿真处理了含噪信号。模拟产生一高斯噪声, 与信号s (n) 叠加, 所得含噪信号sn (n) 的信噪比为1dB。含噪信号sn (n) 的波形和其伪WV分布重排波形如图3所示, 由图可知, 二次型时频分析本身具有一定的抗噪性能, 基本上能提取含噪信号的时频分布。

4 结束语

由本文的仿真结果可知, 时频分析是处理非平稳信号的一种有效方法, 比较线性时频分析结果和二次型时频分析结果, 可知二次型时频分析方法效果较好。另外二次型时频分析具有一定的抗噪性能, 若结合有效的去噪方法, 时频分析方法可广泛用于雷达、生物、医学和故障诊断等领域。图2 (d) 与图1 (b) 比较可知, 所求得的伪WV分布没能完整的表示前两个常频率高斯调幅信号, 信号始末部分的信息丢失, 不同频率信号相叠加的部分也没检测到。这是目前二次型时频分析方法的缺陷, 也是今后研究所要解决的问题。

摘要：使用时频分析方法处理非平稳信号, 提取非平稳信号的时频分布。分别采用线性时频分析和二次型时频分析处理模拟仿真信号, 仿真结果显示用二次型时频分析处理所得的时频分布的时频分辨率较高。仿真结果证明, 二次型时频分析方法是处理非平稳信号的有效方法。

关键词：时频分析,非平稳信号,信号处理

参考文献

[1]Rajagopalan S, Restrepo J A, Aller J M, et al.Nonstationary MotorFault Detection Using Recent Quadratic Time-Frequency Representa-tions[J].IEEE Transactions on Industry Applications, 2008, 44 (3) :735-744.

[2]Saad N M, Abdullah A R, Yin Fen Low.Detection of Heart Blocks inECG Signals by Spectrum and Time-Frequency Analysis[C].4th Stu-dent Conference on Research and Development, June 2006:61-65.

[3]葛哲学, 陈仲生.Matlab时频分析技术及其应用[M].人民邮电出版社, 2006.

[4]Gabor D.Theory of communication[J].IEEE, 1946, 93:429-457.

非平稳时间序列 第5篇

在地震勘探中, 地震波在地下介质中传播时, 由于散射和大地滤波的作用, 使地震波非平稳、非单一频率的特征, 即实际地震资料是频率随时间变化的非平稳信号。时频分析方法是非平稳信号常用的处理方法。本文选取其中的连续小波变换、广义S变换、子波匹配追踪这三种方法进行了理论说明, 并在地震非平稳信号上的应用及分析优缺点。

二、三种时频分析方法原理简介

2.1连续小波变换

20世纪80年代末首先由法国地球物理学家J.Morlet提出, 之后再由理论物理学家A.Grossmann将小波变换理论完善。小波变换是以一个形状可变但是面积固定的窗口在信号上滑动的一种局部时频域分析方法。

信号x (t) 的连续小波变换 (CWT) 表示为

式 (1) 中是ω的复共轭, ω是满足容许条件的母小波函数, a是伸缩因子, b是平移因子。

常用的母小波函数是Morlet小波, 其表达式为

式 (2) 中f0是Morlet小波的中心频率。

2.2广义S变换

1996年首先由R.G.Stockwell等人提出S变换。S变换是连续小波变换方法的一种改进。陈学华在Gabor的短时傅里叶基础上提出的变窗函数广义S变换。

1946年, Gabor提出得短时傅里叶变换表示为

其中, x (t) 为原始信号, ω (t) 为高斯函数, τ为时间参数。窗函数的选择决定了短时傅里叶变换的时频分辨率。

现定义上式中的窗函数为

其中δ是确定窗函数时间宽度的尺度因子

令时窗尺度因子为

其中λ、p为调节因子 (λ>0, p>0)

上式 (4) 和 (5) 代入式 (3) 中, 可以得广义S变换的表达式为

2.3子波匹配追踪算法

1993年首先由Mallat等人提出的。其算法是基于Gabor函数组成的超完备子波库贪婪算法, 并以子波库中各匹配子波的Wigner-Ville分布之和来对原始信号进行时频表征。

设输入信号为s (t) 。

匹配追踪算法每迭代一次即提出一个最优子波wγn, n为迭代次数, 经过N次迭代之后有

其中, an是第n个匹配子波wγn的振幅, R (N) s是迭代N次之后的残差。

其中WVD|wγn (t, f) |为匹配子波wγn的Wigner-Ville分布。

三、地震非平稳信号上的实际应用

在实际地震资料中取出一道地震记录来进行时频分析, 有下面几张图。

图1即为实际地震资料中一道非平稳的地震记录。图2是基于Morlet小波变换方法做出的时频谱, 小波变换在多尺度方面效果突出, 它能在给出比较好的时间精度时同时保持好频率精度, 然而小波变换在时间方向上仅仅是瞬时振幅的体现, 体现不出相位信息。图3是广义S变换的时频谱, 广义S变换也能分辨多频率信号, 但是广义S变换窗口的选择使时频谱在频率间断的地方上频率刻画不准确。图4是Ricker子波匹配追踪算法做出的时频谱, 不仅解决了小波变换时频方法没有相位信息的物理问题, 同时相比广义S变换有更高的时间精度和频率精度, 但计算量大是现今面临的问题。

四、结束语

通过分析比较以上三种时频分析方法, 可以得出以下结论:

(1) 连续小波变换本身存在的无相位信息这个问题不可避免。

(2) 广义S变换尽管解决了小波变换不具有物理意义的这个问题, 但是广义S变换的窗口选择也很重要。

非平稳时间序列 第6篇

由于降雨过程属于非平稳随机过程, 降雨量的大小受到多种复杂因素的影响, 因此要对降雨量进行准确的预测是非常困难的。平稳时间序列预测在降雨量的预测中应用较为广泛, 但预测精度仍然不够理想。近年来发展起来的小波分析方法在水文水资源预测等领域有了广泛的应用, 并取得了良好的效果[1,2,3], 但在降雨量预测中的应用则较少。因此, 本文尝试在基于小波消噪的基础上, 应用平稳时间序列方法对降雨量进行预测, 以提高预测的精度。

1 平稳时间序列预测模型

1.1 平稳时间序列预测方法基本原理

平稳时间序列是平稳随机过程的一个现实[4], 平稳随机过程是随机过程的一种特殊形式。如果一个时间序列的概率分布与时间无关, 则称该序列为严格的 (狭义) 平稳时间序列。平稳时间序列预测是建立在对历史数据分析之上, 历史数据越多, 越准确, 预测也越可靠。平稳时间序列的主要特点是过程的统计性质不随时间的平移而变化[5]。通常, 将时间序列描述成一个有序的数列:y1, y2, …, yt。其中下标t表示时间序号。

1.2 平稳时间序列预测方法建模步骤

(1) 由{yt}的一组观测值y1, y2, …, yN, 先按状态分类, 然后计算样本平均数和各个数据的距平值。即:

$\{\begin{array}{l}\overline{y}=\frac{1}{{\rm N}}{\displaystyle \sum _{i=1}^{{\rm N}}y}{}_i\\ y{}_i^{´}=y{}_i-\overline{y}\end{array}(1)$

(2) 由相关函数的性质及计算误差精度确定预报方程阶数m, 实际应用时预报阶数可取

$m=\min \{{\rm N}-1,({\rm N}/4)\}(2)$

(3) 在最小二乘法原则的基础上, 确定线性最优预报方程:

$\widehat{y}{}_{{\rm N}+\tau }^{´}=\alpha {}_{1}y{}_{{\rm N}}^{´}+\alpha {}_{2}y{}_{{\rm N}-1}^{´}+\cdots +\alpha {}_{m}y{}_{{\rm N}-m-1}^{´}(3)$

式中:τ取正整数 (外推步长) 。

(4) 由下式:

$B(k)=\frac{1}{{\rm N}-k}{\displaystyle \sum _{i=1}^{{\rm N}-k}y}{}_{i+k}y{}_i^{´}$

计算B (0) , B (1) , …, B (m) 。

(5) 解线性方程组:

$\{\begin{array}{l}B(0)\alpha {}_1+B(-1)\alpha {}_2+\cdots +B(1-m)\alpha {}_m=B(\tau )\\ B(1)\alpha {}_1+B(0)\alpha {}_2+\cdots +B(2-m)\alpha {}_m=B(\tau +1)\\ ⋮\\ B(m-1)\alpha {}_1+B(m-2)\alpha {}_2+\cdots +\\ B(0)\alpha {}_m=B(\tau +m-1)\end{array}(4)$

(6) 取定外推步长τ后, 得出最优估计[6,7]。

2 小波消噪

2.1 小波消噪原理

小波是一种建立在泛函分析、调和分析、傅里叶分析基础上的视频原子, 它在时域和频域同时具有良好的局部化特性和多分辨率特性, 常被誉为信号分析的“数学显微镜”。近10多年来, 小波分析的理论和方法在许多领域得到了蓬勃发展和广泛应用[8]。

一般的, 有用信号通常为低频信号或是一些平稳的信号, 而噪声信号表现为高频信号。所以消噪过程主要就是通过小波分析, 将高频成分和低频成分有效分离出来, 根据不同信号在小波变换后表现出的不同特性, 对小波分解序列进行处理, 将处理后的序列加以重构, 即可实现对原始信号的消噪处理。

2.2 小波消噪的步骤

目前小波消噪主要采用阈值消噪, 主要步骤分为3步[9]:

第1步:序列的小波分解。选择一个合适的小波函数并确定小波分解的层次N, 然后对原始序列进行分解计算, 得到高频系数和低频系数。

第2步:小波分解高频系数的阈值量化。对各个分解尺度下的高频系数选择一种阈值进行阈值量化处理。

第3步:小波重构。根据小波分解的最低层低频系数和各层高频系数进行小波重构, 得到消噪后的序列。

3 基于小波消噪的平稳时间序列预测模型

基于小波消噪的灰色拓扑预测方法基本思想是:首先利用小波分解和重构消除降雨量序列中的噪声, 然后对重构的降雨量序列应用平稳时间序列预测方法进行预测。

小波消噪的关键环节是选择适合的小波函数和分解层数, 本文选用Daubechies (dbN) 小波函数。Daubechies小波是由世界著名的小波分析学者Inrid Daubechie构造的小波函数, 具有非常重要的性质, 他不仅是连续的和正交的, 而且在分解与重构算法中所需要的计算量少, 这在信号处理中非常重要。因此, 本文采用db3小波对降雨量序列分别进行一层分解, 分解过程中阈值的确定采用Stein (斯坦) 的无偏估计阈值方法, 整个分解过程应用MatLAB软件来实现[10]。

本文采用平稳时间序列的线性外推法对丹东地区降雨量进行预测, 依此来预测旱涝灾害。以1971-2000年的降雨量为时间序列, 应用1.2节的方法建立自回归预测模型, 以2001-2006年的降雨量为检验值进行拟合。

4 实例分析

4.1 项目区概况

丹东地区位于辽宁省东部, 地处辽东半岛经济开发区东南部的鸭绿江畔, 属辽东山地南坡, 呈扇形, 地势北高南低, 依次为中低山、丘陵、平原, 呈阶梯状分布, 地理坐标为东经123°22′30″～125°42′, 北纬39°43′31″～41°09′21″, 东西长195.85 km, 南北宽160.2 km。丹东地区属暖温带季风型大陆性气候, 夏无酷暑, 冬无严寒, 四季分明, 雨热同季, 位于欧亚大陆东岸中纬度地带。全市降雨受地形等因素的影响, 年降雨量的分布自西南向东北递增, 是辽宁省年降雨量的最高值, 也是东北地区的暴雨中心。同时降雨时空分布很不均匀, 80%以上集中在6-9月份, 导致丹东地区旱涝灾害频繁。

4.2 降雨量预测

丹东地区1971-2006年的降雨量序列如表1所示, 以1971-2000年的降雨量为原始序列建模, 以2001-2006年的降雨量为检验值进行拟合。

mm

首先通过做自相关图 (ACF) 对丹东地区1971-2000年的实测降雨量序列的平稳性进行分析, 为了提高精度, 对明显没有落入置信区间的异常值弃掉不用。去掉异常值后的降雨量时间序列的ACF如图1所示。图1中纵坐标ACF为自相关系数, 横坐标Lag Number是滞后值的相关系数。由图1可见, 去掉异常值的样本序列数据的自相关系数很快落入置信区间, 序列基本平稳。

然后应用db3小波对降雨量进行分解, 重构降雨量序列。由于已具备一定的平稳性, 噪声污染较轻, 因此只对其进行一层分解, 分解过程如图2所示。

最后对小波消噪后重构的降雨量序列再采用平稳时间序列的线性外推法进行预测。采用stain无偏估计法分别对分解过程的高频系数进行软阈值量化处理, 经计算一层分解的系数阈值为T1=2.512 8, 结合低频系数进行小波重构, 即得到消噪后的降雨量序列。由新的序列建立平稳时间序列的回归模型, 对降雨量进行预测。预测结果及对比情况如表2所示。

5 结 论

(1) 通过预测结果对比分析可以看出, 应用小波消噪后的平稳时间序列预测方法, 其预测结果除了2001年的相对误差有所增加外, 其他的年份预测值的相对误差有了很明显的降低, 因此基于小波消噪的平稳时间序列预测方法对丹东地区年降雨量进行预测是可行的。

(2) 小波消噪的关键环节是小波函数的选择和分解层数的确定, 本文只应用了一种小波函数进行了一层的尝试, 对于应用其他小波函数的效果有待进一步的研究。另外, 由于不同地区的影响降雨量大小和时空分布的因素不同, 所以各种预测模型在不同地区的适用性也不尽相同, 因此在应用时, 要根据不同地区的降雨量规律选择适合当地的小波函数。

(3) 由于丹东地区降雨量的年内分配极不均匀, 所以往往在年降雨量趋于正常的年份中, 也时有旱涝灾害发生。针对这种现象, 可以应用上述方法对丹东地区的降雨量进行逐月的预测, 并与年降雨量预测值综合考虑。这样不仅能对预测年份的降雨量有个定量的认识, 还能对降雨量的年内分配有定性的分析。特别是在降雨量对农作物生长影响较大的月份, 降雨量的预测结果还可以与该地区其他气象因素的多年平均值综合考虑, 以便更详细地分析旱涝灾害发生的等级状况, 为科学决策和防灾减灾提供依据。

摘要：小波消噪与时间序列分析方法在预测领域中应用十分广泛, 但是在降雨量的预测中应用不多。在基于小波消噪的基础上应用时间序列中平稳时间学列方法对降雨量进行预测, 结果显示, 应用该方法有效地提高了降雨量的预测精度。用丹东地区1971-2006年的降雨量作为历史数据, 建立降雨量预测模型, 结果表明新模型算法简单、精度较高, 比传统的拓扑预测模型效果更好, 为降雨量预测提供了一种行之有效的方法。

关键词：小波消噪,平稳时间序列,降雨量预测

参考文献

[1]汤成友, 王瑞, 缈韧.离散小波变换在水文序列分解中的应用[J].中国农村水利水电, 2007, (2) .

[2]曹梦成, 刘小玲.小波去噪的GPS算法在振动变形监测中的应用[J].中国农村水利水电, 2008, (4) :88-91.

[3]刘东, 付强.小波最近邻抽样回归耦合模型在三江平原年降雨预测中的应用[J].灌溉排水学报, 2007, (4) :82-85.

[4]温品人.时间序列预测法的实际应用分析[J].江苏广播电视大学学报, 2001, (12) :64-65.

[5]刘乐平.时间序列分析在灾变预测中的应用[J].华东地质学院学报, 1997. (6) :191-195.

[6]张宁宁.朝阳地区干旱特征分析及抗旱对策研究[D].沈阳:沈阳农业大学, 2006.

[7]张世英, 陆晓春, 李胜朋.时间序列在城市交通预测中的应用[J].天津大学学报 (社会科学版) , 2006, 8 (5) :370-372.

[8]葛哲学, 沙威.小波分析理论与MATLABR2007实现[M].北京:电子工业出版社, 2007.

[9]葛哲学, 陈仲生.Matlab时频分析技术及其应用[M].北京:人民邮电出版社, 2006.

[10]张善文, 雷英杰, 冯有前.MATLAB在时间序列分析中的应用[M].西安:西安电子科技大学出版社, 2007.

[11]李毅, 周慧, 李秦.基于R/S分析的参考作物腾发量时间序列分析特征[J].节水灌溉, 2009, (4) :4-7.

非平稳时间序列 第7篇

在临床医学领域中, CFI能够实时准确无创地探测血流速度, 是诊断心血管疾病的重要技术[1]。在CFI中, 探头接收到的回波信号包括血流信号、杂波信号 (由血管搏动和组织慢速移动引起) , 通常杂波信号强度比血流信号强度高出40~80 d B[2], 这给准确地估计血流速度带来了极大困难。因此, 为了得到真实可靠的血流速度, 必须对回波信号中的杂波进行充分地抑制。

目前常用的杂波抑制器有:低阶FIR滤波器[3]、基于投影初始化的IIR滤波器[4]和回归滤波器[5]。其中低阶FIR滤波器、投影初始化IIR滤波器和回归滤波器均属于静态滤波器, 它们的通带截止频率、阻带衰减等特性较为固定。当杂波信号是平稳信号时, 这些滤波器能够获得较为理想的效果。但在实际临床诊断中, 由于人体呼吸、脉搏等因素造成组织加速运动, 导致杂波信号属于非平稳随机信号, 静态滤波器对此类信号的处理效果不理想。后来有学者提出非平稳杂波抑制法[1], 该方法对杂波抑制效果较好, 但是在杂波信号较弱的区域会造成误消除, 导致血管内血流信息的损失, 对血流速度的估计造成极大误差。

本文在非平稳杂波抑制法的基础上提出一种基于动态区域划分的非平稳杂波抑制方法。该方法首先利用能量特性对回波信号进行动态区域划分, 然后结合非平稳杂波抑制法和多项式回归法对不同区域的信号分别进行处理。实验结果表明, 该方法的算法复杂度低, 满足超声诊断设备的实时性要求;同时, 该算法能较好地抑制杂波并保证流速度剖面的完整性, 是一种快速有效的杂波抑制方法。

1 原理与方法

1.1 动态区域划分和局部平滑

多普勒回波信号是在同一探测位置发射K次脉冲波得到的一维向量矩阵。将轴向采样容积为M的K个复解调信号排成二维矩阵, 则对回波信号处理的过程即是对二维矩阵进行运算。其中, 信号构成如图1所示。通常, 信号矩阵的行向量被称为慢时信号, 列向量被称为快时信号[2]。

复解调回波信号可表示为矩阵X:

其中K是慢时信号方向的采样容积 (脉冲重复数) , M是快时信号方向的采样容积。

回波信号慢时信号方向的平均能量E_X表示为:

同理, 快时信号方向上不同深度的慢时信号能量E_Xs (m) 可由式 (3) 表示:

其中j是慢时信号方向上的第j个采样点, m是快时信号方向上的第m个采样点。由于不同区域的组织产生的回波信号强度不同, 以式 (4) 的方式对信号进行划分, 可以区分不同组织区的信号。

通常在动态组织区, 杂波信号的强度比血流信号强度高出40~80 d B, 因此可以由式 (5) 再次对动态组织区的信号进行划分:

其中E_M是动态组织区快时信号上慢时方向的平均能量。

静态组织区远离血管, 整体趋于静止, 所以回波信号趋于平稳, 且不携带任何速度信息;血流区分布在血管内部, 加之血管中的血流速度呈抛物线分布且相对较稳定, 因此血流区的信号相对平稳;杂波区受呼吸、脉搏等影响, 造成血管及血管附近组织具有加速度, 因此杂波区信号属于非平稳信号。

为保证杂波滤波器自适应于组织运动且不造成血流信息的损失, 本文算法采用空间平滑法处理动态组织区, 如式 (6) :

其中N是快时信号方向上参与空间平滑的区间半径。

1.2 杂波频率估计[1]

复解调回波信号x (m, k) 用复指数形式表示为:

其中A (m, k) 是信号的幅值, φ (m, k) 是信号的相位。通常在多普勒回波信号中, 杂波信号的相位可以用一个低阶多项式表示:

对式 (8) 求导即可求得杂波的瞬时频率w (m, k) :

由此可见, 对杂波信号的瞬时频率进行估计, 就是求出多项式系数ad的过程。

引入矩阵R和向量φ, 对拟合多项式 (9) 求解得多项式系数向量a:

又由式 (9) 可得杂波瞬时频率矩阵W=[w (m, 1) w (m, 2) …w (m, K) ]T的表达式:

将式 (10) 得到的系数向量a带入式 (11) 就可以进一步求得杂波的瞬时频率矩阵W。对于D的取值, 参考文献[5]通过分析低阶多项式拟合的原理并进行大量的实验后得出:当D≤4时, 多项式能够有效地拟合低频信号。在超声彩色血流成像中, 当D=3时, 拟合多项式能够精确地描述杂波的瞬时频率。

1.3 局部非平稳杂波抑制

多普勒回波信号经过动态区域划分和杂波瞬时频率估计后就可进行混频处理。其中混频是为了将非平稳杂波信号的频率搬移到零频附近。为保证血流信号点的完整性, 本文算法只对杂波区信号进行混频处理。利用前文得到的杂波瞬时频率, 混频信号可表示为:

混频过程可用式 (14) 表示:

通过混频处理, 杂波信号的非平稳性大大降低, 此时通过常用的静态滤波器就可以对杂波信号进行充分抑制。由于血流区信号可能含有较弱的杂波, 因此本文算法选用能够避免数据点损失的多项式回归滤波器对杂波信号进行抑制, 这样既保证了杂波信号的有效去除, 又保证了血流信息的完整性。算法流程如图2所示。

2 仿真实验

2.1 仿真模型与参数

为了分析和验证所提出的杂波抑制方法的有效性, 本文在计算机上用MATLAB R2012a软件平台进行仿真, 运行环境为Centrino2 (1.66 GHz) , 采用丹麦技术大学学者Jenson等研制的Field ii软件包仿真超声多普勒回波信号。仿真所用的超声血流模型如图3所示, 仿真参数如表1所述。图4是通过仿真模型得到的一条回波信号波形图。

2.2 仿真实验结果

为验证本文算法对杂波的抑制效果, 将其与投影初始化滤波器、多项式回归滤波器、SVD滤波器和非平稳滤波器进行比较。图5所示的是一组 (16条) 回波信号经过各滤波器后自相关估计得到的相应的血流速度剖面波形图, 图6所示的是218组相邻回波信号通过各滤波器后自相关估计、编码映射得到的血流速度图。两图直观地展示了不同滤波器对杂波的抑制性能。

由图5 (e) 可以看出, 采用本文算法进行杂波抑制后, 自相关估计得出的血流速度信息几乎仅存在于血管内部 (约在扫描深度24 mm~34 mm之间, 血管理论直径为8 mm) , 且速度分布相对比较完整, 这与实际情况基本相符。由图5 (a) ~5 (d) 可以看出, 经过其他滤波器处理后的回波信号经自相关估计后得到的血流速度信息不仅存在于血管内部, 也存在于组织区, 这显然与实际情况不符。

从图6 (e) 可以看出, 本文算法得到的血管壁清晰、完整, 血管截面速度分布均匀;从图6 (b) 和6 (d) 可以看出, 经过回归滤波器和非平稳滤波器得到的血流速度分布不均匀, 且组织区出现伪血流信息;从图6 (c) 中可以看出, 经过SVD滤波器后得到的血流速度分布不完整, 血管内速度出现零点, 且血管直径明显小于理论值, 成像效果较差;从图6 (a) 可以看出, 投影初始化IIR滤波器由于采用固定频率作为截止频率而造成数据点损失较大, 使得血管直径明显小于理论值 (理论血管直径为8 mm) 。

此外, 为了从客观上对本文算法进行评价, 文章从最大速度估计、杂波血流比和算法运行时间三方面对算法进行比较。比较结果如表2所示。

从表2中可以看出, 本文算法估计的最大速度值接近理论最大速度值0.972 m/s (本文算法的最大估计速度与多项式回归滤波器相同, 这是因为该算法采用多项式回归法对血管中心区域信号进行处理) , 表明采用本文算法对杂波信号进行抑制时, 原始信号中的血流信号损失较小;经本文算法滤波后, 信号的杂波血流比明显低于其他几种滤波算法, 这表明该算法能够有效地对原始信号中所含的杂波信号进行抑制;在运行时间上, 本文算法的运行时间远小于其他几种算法, 这是由于该算法对静态组织区的信号不做处理导致的。综上所述, 本文算法在杂波抑制效果、血流信息的保留和算法时间复杂度上都能取得较为满意的结果, 由此证明此算法较其他杂波抑制算法更有效。

3 结论

本文提出的基于动态区域划分的非平稳杂波抑制方法一改以往采用单一滤波方法进行杂波抑制的模式, 通过分析回波信号的特点, 首先采用动态区域划分法将信号分割为不同部分, 再根据各部分信号的特点做出相应处理。本文算法采取对静态组织区信号完全保留的方法, 既保证滤波器不会在此区域引入噪声, 又保证了运行速度。同时, 对杂波区和血流区采取不同的滤波方法, 既保证了杂波的有效去除, 又避免了血流信号的损失。此外, 该算法采用矩阵进行推理运算, 易于工程实现。

另外, 杂波瞬时频率的估计是基于低阶多项式来实现的, 当多项式阶数取值适当时, 能较好地描述杂波频率, 但当阶数取值欠佳时, 不能很好地对频率进行描述。因此, 采取一种能精确描述杂波瞬时频率的方法将是下一步研究的目标。

参考文献

[1]王沛东, 沈毅, 王艳.超声彩色血流成像中非平稳杂波的抑制[J].中国生物医学工程学报, 2008, 27 (2) :240-243.

[2]尤伟, 汪源源.超声彩色血流成像中基于动态区域划分抑制杂波的方法[J].仪器仪表学报, 2010, 31 (3) :594-599.

[3]JENSEN J A.Estimation of blood velocities using ultrasound[C].Cambridge, UK:Cambridge University Press, 1999:195-226.

[4]CHORNOBOY E S.Initialization for improved IIR filter performance[J].IEEE Trans.Signal Processing, 1992, 40 (5) :543-550.