基于s变换的presley非平稳随机过程演变功率谱密度

基于s变换的presley非平稳随机过程演变功率谱密度

1时间-尺度分析工具之二:st果变换随机噪声的研究通常涉及对随机过程的描述。对于非平稳随机过程而言,其非平稳性不仅依赖于时间的幅值,而且其频率成分也是随着时间变化的。观察典型的地震动过程,其幅值表现为上升-平稳-衰减的趋势,此现象称之为时域非平稳;而由于地震波在传播过程中的行波效应或高频成分不断被吸收导致地震动时程越零率不断减少,此现象称之为频域非平稳。获取非平稳随机过程的时间-频率特性,成为描述和模拟随机过程的重要环节。最早的平稳随机过程频域描述为Fourier变换。因其基函数为谐和函数,决定了Fourier系数在物理意义上能解释为频率。由于借助谐和函数可为诸多物理和数学方程提供解析解,决定了Fourier变换在自然科学领域数百年时间内的重要地位。但Fourier系数却不具有时间定位特征,或更具体地讲,其时间特性湮灭在相位信息之中。1946年,Gabor提出了首个时间-频率分析工具,即Gabor展开(Gaborexpansion),Gabor展开的系数可以通过短时Fourier变换STFT(ShortTimeFourierTransform)得到。其基本思想是,为待分析信号加窗,然后,对每个加窗的子信号进行Fourier变换。这样,所得到的加窗分段Fourier变换便具有了时间信息。其缺点在于,当窗函数给定后,对于每个时间点和频带,其时间和频率分辨率是无法改变的。Wigner和Ville将经典的平稳随机过程谱理论延伸于非平稳随机过程方面,提出了Wigner-Ville分布WVD(Wigner-VilleDistribution)。在这一时间-频率分析工具中,时变谱是通过时变相关函数的Fourier变换得到的。一般来讲,由于WVD的变换核仍为时间无限长的谐和函数,因此,从本质上无法合理反映局部的时间-频率信息;而且,在某些情况下会有负值出现,因此很难在物理上给出合理的解释。小波变换WT(WaveletTransform)的提出,克服了上述时间-频率分析工具的缺点。小波变换是一种时间-尺度分析工具,依赖于尺度因子,每个小波相对于“母小波”进行拉伸或压缩;依赖于平移因子,每个小波相对于“母小波”进行平移。因此,尺度因子提供了局部的频率信息,而平移因子提供了局部的时间信息。但是,在特定尺度下,并非所有小波基函数都能与频率显示联系,比如正交低阶Daubechies小波,其在频域内表现为非紧支且多峰的情况。因此,虽然小波分析能很自然地进行时间-尺度分析,但在利用其进行时间-频率分析时,必须注意尺度与频率之间的联系。最近提出的一种时间-频率分析工具,即Stockwell变换ST(StockwellTransform)可以视为WT与STFT的结合。其优点在于:(1)提供了类似于WT的多尺度的恒Q(ConstantQuality)分析;(2)同时提供了STFT才具有的绝对相位信息。近年来,S变换已广泛应用于地球物理,地震和医学信号处理,电子和机械工程以及频谱分析。关于S变换的研究,Gibson等分析了S变换与小波变换之间的区别及联系;Pinnegar等对S变换提出了一些改进;Wong等提出了广义S变换并对S变换的相位信息进行了进一步考察。在非平稳随机过程的演变功率谱估计方面,Iyama等和Basu等建立了不同尺度小波变换均方值与随机过程时变谱之间的关系。Spanos等得到了谐和小波变换与非平稳随机过程演变功率谱密度之间的联系;与此同时,Spanos等结合Priestley非平稳随机过程,利用一般非正交小波估计了非平稳随机过程的演变功率密度,此文的根本在于,相对于一般小波基函数,Priestley非平稳随机过程的调制函数为相对慢变函数。2009年,Huang等将这种方法推广到多变量非平稳随机过程的情况。受文献启发,本文提出了一种利用S变换估计Priestley非平稳随机过程演变功率谱密度的方法。方法的根本在于,相对于S变换的“变换核”,Priestley非平稳随机过程的调制函数为相对慢变函数。因此,非平稳随机过程的S变换可视为相位修正后的另一非平稳随机过程。进而推导出了以非平稳随机过程演变功率谱密度表达的在特定频率点处的S变换瞬时均方值。考虑所有频率点的表达式,且通过将前者写为有限个频率点的级数展开,则以S变换瞬时均方值显式地表达非平稳随机过程演变功率谱密度的问题,转化为了求解一组代数方程。此代数方程组的解为级数展开中的依赖于时间的未知项,非平稳随机过程的演变功率谱密度随之而得。2速度s变换与stft的关系为了与S变换的原始文献保持一致,本文中所述推导式中的频率除非明显标出,均为自然频率。如果确定性过程为x(t)∈L2(R),则其S变换可以表达为式中f为频率。可见,S变换的核函数是经过高斯函数调制后的谐和函数,即与小波变换的相同之处在于,其高斯窗函数的宽度是随着频率而变化的,因此导致了其具有多尺度分辨率;与STFT的相同之处在于,其核函数是随着平移因子τ而变化的,因而使其具有绝对相位信息。由于S变换可以看作为卷积形式,利用Fourier变换的卷积定理,可以得到S变换的频域表达式为式中X(α)为确定性过程x(t)的Fourier变换。基于式(3),其离散S变换可表达为式中T为确定性过程x(t)的时间采样间隔,N为其采样总点数,Sj[T,n/NT]为变换的采样点。研究表明,如采用快速Fourier变换(FFT)技术,能很大程度上加快S变换的计算速度。可以证明,如果Sxφ(τ,f)为确定性过程x(t)的S变换,则其平移过程x(t-r)的S变换可以表达为S变换还具有其他重要性质,限于篇幅限制,这里只给出对于估计随机过程功率谱重要的性质,其余性质可参见文献。3子问题的估计考虑一个形如的非平稳随机过程,其中A(ω,t)为依赖于时间的慢变函数,Z(ω)为复值正交增量的随机过程,即式中E为期望算子,SX-X-(ω)为相应的平稳随机过程的双边演变功率谱密度函数。因此,当调制项A(ω,t)为时间慢变函数时,非平稳随机过程X(t)的演变功率谱密度函数可以近似表达为下文中常把非平稳随机过程的功率谱近似理论表达(式(9))称为目标功率谱,并与估计功率谱对比以验证所提方法的合理性。考虑某特定频率fj处的S变换,并将式(6)代入式(1)可得根据Priestley非平稳随机过程模型(式(6))的定义,A(ω,t)相对于S变换核为时间慢变函数。具体而言,图1所示为典型地震工程中使用的某调制函数与S变换核对比。可以看出,当S变换频率点选定后,在相应的S变换核有效积分区间内,调制函数可视为常值,因此可近似将此慢变函数移至积分号外,同时其时间变量取为高斯调制函数的对称轴t=τ,即如将式(11)内部积分单独考虑,显见,式(12)即为h(t)=exp[i2πω(t-τ)]的S变换。利用S变换的平移特性(式(5))可知,如果u(t)=exp(i2πωt)的S变换已知,则h(t)的S变换亦可求得。利用S变换的频域表达,u(t)的S变换可写为进而,h(t)的S变换可以表示为将式(14)代入式(11)可得式中保留了dZ(ω)的正交特性,亦为正交随机过程,即若令W(τ,fj)=S(τ,fj)exp(i2πfjπ),则W(τ,fj)亦为Priestley非平稳随机过程。因此,原非平稳随机过程在某特定频率处的S变换的相位修正可认为是连续变量的非平稳随机过程,且此随机过程的功率谱密度函数为因此,这一随机过程的瞬时均方值为4基于时间条件的变分迭代式(18)建立了在某频率点fj处原非平稳随机过程的S变换的均方值与该随机过程演变功率谱密度的关系。考虑若干频率点上的S变换的相位修正项,式(18)可以写为式中N为S变换的频率点取样总点数。如果将原非平稳随机过程写为的级数展开。其中cj(τ)为某频率点的时变系数。结合式(20)与式(19),可得式中E[|S(τ,f)|2]为N维列向量,c(τ)为N维列向量,Q为N×N维矩阵,即式中每个元素Qi,j为由此可知,通过解矩阵方程(21),可得时变系数c(τ)。因此,可由式(20)给出的非平稳随机过程的演变功率谱估计。在矩阵方程(21)中,形状函数矩阵Q是不依赖时间τ的。因此,对于所有时间取样点,此矩阵只需求取一次逆,且由于基于FFT技术的S变换的引入,故本文建议方法具有较高的计算效率。值得特别指出的是,上述求解过程并不需要随机过程的时-频功率谱函数的具体形式,因此,特别适合于利用实际观测记录估计非平稳随机过程的时-频谱。当然,若预设某类时-频谱形式,利用上述谱估计方法,可以通过求解反问题方式识别谱参数。有关研究工作,将另文发表。5机过程估计算法本节以两种具有演变功率谱的非平稳随机过程,即均匀调制和非均匀调制的非平稳随机过程为例,利用本文所建议的算法,估计其演变功率谱密度。5.1平稳过程/功率谱密度的定义考虑如下均匀调制的非平稳随机过程,即式中g(t)为只依赖于时间t的调制函数,为平稳随机过程,其双边功率谱密度为Kanai-Tajimi谱式中S0=1.0cm2s-3rad-1为基岩白噪声输入强度,ζs=0.4和ωs=20rads-1分别为场地土阻尼与卓越频率系数。调制函数g(t)定义为如图1所示,其中a=0.25,b=0.5,k为使gmax=1的正规化系数。据式(9)可知,均匀调制非平稳随机过程的目标演变功率谱密度可以表达为可以看出,不论在样本层次上,还是在功率谱意义上,均匀调制非平稳随机过程的样本/功率谱都能由相应平稳过程的样本/功率谱乘以仅依赖于时间的调制函数(/的平方)而得。因此,文献[20,23,24]也将调制项A(ω,t)=g(t)的均匀调制随机过程称为可分离(Separable)非平稳随机过程或“准平稳”(Quasi-stationary)随机过程。在本文建议方法中,确定性样本函数的生成为基于S变换的谱估计的基础。在此,平稳过程的样本函数用谱表达方法生成,即式中Δω=ωu/n为频率取样间隔,ωu为上限截止频率,φn为[0,2π]之间均匀分布的随机变量。在谱表达方法中,样本函数数目为500个,时间取样点数目为1024,取样间隔为Δt=0.02s,截止频率取为Nyquist频率ωu=π/Δtrad/s,谐和项数目为512。据第3节所述估计功率谱密度的方法(式(21)和式(20)),图2所示为该非平稳随机过程在6s处的瞬时目标功率谱密度与瞬时估计功率谱密度对比,可以看出,二者吻合良好。为了更清楚地表现在其他时刻处的瞬时功率谱密度,图3和图4分别给出了目标演变功率谱密度和估计演变功率谱密度曲面图。5.2非平稳随机过程的演化过程算例1为均匀调制的非平稳随机过程的功率谱密度函数,其频谱特征并不随时间变化,只是幅值随着调制函数变化。因此,这种均匀调制的随机过程并不符合实际地震动的时变频谱特征。正如引言中所述,由于行波效应或传播路径中高频成分不断被吸收等因素,随机地震动的瞬时频谱特征是依赖于时间的。因此,利用本文建议方法估计具有演变谱随机过程的功率谱密度函数,更具有实际意义。为此,考虑具有目标演变频谱密度的非平稳随机过程为其演变功率谱密度如图5所示。可以看出,在不同的时间段,其功率谱密度不仅在幅值上,且其卓越频率所在区间也随时间变化,式(29)较为充分地反应了地震动的幅值及频谱特点,即非平稳不仅存在于幅值方面,而且存在于频谱方面。由式(29)亦知其并不能表示为如式(27)所示的时间-频率分离形式,因此,这种同时反应时间-频率非平稳特征的随机过程,在诸多文献中称之为不可分离(Non-separable)非平稳随机过程。非平稳样本函数由谱表达方法合成,即式中SXX(t,kΔω)为时变的目标演变功率谱密度函数,其他诸参数与式(28)相同。图6为利用本文建议方法所估计的演变功率谱密度的三维视图;图7为在4s和12s时刻的目标瞬时功率谱密度与估计瞬时功率谱密度对比。可以看出,估计功率谱密度与目标功率谱密度吻合良好。5.3非平稳随机过程随时间的变化显然,调制函数相对于S变换核的慢变特性对于演变功率谱估计值有较大影响。但本文所建议方法,对于大多数工程实际中幅值非快变的非平稳随机过程均适用。另一方面,合适地选择S变换频率点使其变换核较之调制函数为快变亦至关重要,如图1所示。如何定量描述由S变换核和调制函数相对变化快慢而导致的估计值偏差,可由后续研究考察。其次,由式(20)可知,非平稳随机过程瞬时功率谱密度展开为若干高斯函数平方的加权和。参与计算的频率点fj决定了高斯函数的对称轴位置及宽度。数值计算表明,合理选择频率点fj决定了估计瞬时功率谱的光滑程度。图2所示估计功率谱中出现的若干小幅抖动实为式(20)中不同高斯函数峰值;图7中4s瞬时功率谱峰值附近处,目标与估计功率谱的差别亦由此导致。最后值得注意的是,理论上,式(6)之非平稳随机过程的精确时变谱密度并非为式(9)之演变谱密度。考察Priestley非平稳随机过程可知,只有当式(6)中的调制项A(ω,t)为时间慢变函数时,其时变功率谱密度才能近似以式(9)表示。因此,Priestley定义的非平稳随机过程为一类幅值慢变的振荡(Oscillatory)非平稳随机过程。只有当调制项A(ω,t)依时间慢变时,才不会对非平稳随机过程功率谱近似值(式(9))带来时间-频率域内的偏移(Shift)。换言之,当调制项A(ω,t)依时间非慢变时(其自身带有频率成分),须考虑由其导致的时间-频率域内的偏移。对演变功率谱分析感兴趣的读者可参见文献。6非平稳随机过程的s变换本文首先总结和回顾了各种时间-频率分析工具,包括Fourier变换、短时Fourier变换、WignerVille分布以及小波变换。针对最近出现的S变换,做了详细描述。结合非平稳随机过程的Priestley模型,提出