基于s變換的presley非平穩(wěn)隨機過程演變功率譜密度

基于s變換的presley非平穩(wěn)隨機過程演變功率譜密度

1時間-尺度分析工具之二:st果變換隨機噪聲的研究通常涉及對隨機過程的描述。對于非平穩(wěn)隨機過程而言,其非平穩(wěn)性不僅依賴于時間的幅值,而且其頻率成分也是隨著時間變化的。觀察典型的地震動過程,其幅值表現(xiàn)為上升-平穩(wěn)-衰減的趨勢,此現(xiàn)象稱之為時域非平穩(wěn);而由于地震波在傳播過程中的行波效應(yīng)或高頻成分不斷被吸收導(dǎo)致地震動時程越零率不斷減少,此現(xiàn)象稱之為頻域非平穩(wěn)。獲取非平穩(wěn)隨機過程的時間-頻率特性,成為描述和模擬隨機過程的重要環(huán)節(jié)。最早的平穩(wěn)隨機過程頻域描述為Fourier變換。因其基函數(shù)為諧和函數(shù),決定了Fourier系數(shù)在物理意義上能解釋為頻率。由于借助諧和函數(shù)可為諸多物理和數(shù)學(xué)方程提供解析解,決定了Fourier變換在自然科學(xué)領(lǐng)域數(shù)百年時間內(nèi)的重要地位。但Fourier系數(shù)卻不具有時間定位特征,或更具體地講,其時間特性湮滅在相位信息之中。1946年,Gabor提出了首個時間-頻率分析工具,即Gabor展開(Gaborexpansion),Gabor展開的系數(shù)可以通過短時Fourier變換STFT(ShortTimeFourierTransform)得到。其基本思想是,為待分析信號加窗,然后,對每個加窗的子信號進(jìn)行Fourier變換。這樣,所得到的加窗分段Fourier變換便具有了時間信息。其缺點在于,當(dāng)窗函數(shù)給定后,對于每個時間點和頻帶,其時間和頻率分辨率是無法改變的。Wigner和Ville將經(jīng)典的平穩(wěn)隨機過程譜理論延伸于非平穩(wěn)隨機過程方面,提出了Wigner-Ville分布WVD(Wigner-VilleDistribution)。在這一時間-頻率分析工具中,時變譜是通過時變相關(guān)函數(shù)的Fourier變換得到的。一般來講,由于WVD的變換核仍為時間無限長的諧和函數(shù),因此,從本質(zhì)上無法合理反映局部的時間-頻率信息;而且,在某些情況下會有負(fù)值出現(xiàn),因此很難在物理上給出合理的解釋。小波變換WT(WaveletTransform)的提出,克服了上述時間-頻率分析工具的缺點。小波變換是一種時間-尺度分析工具,依賴于尺度因子,每個小波相對于“母小波”進(jìn)行拉伸或壓縮;依賴于平移因子,每個小波相對于“母小波”進(jìn)行平移。因此,尺度因子提供了局部的頻率信息,而平移因子提供了局部的時間信息。但是,在特定尺度下,并非所有小波基函數(shù)都能與頻率顯示聯(lián)系,比如正交低階Daubechies小波,其在頻域內(nèi)表現(xiàn)為非緊支且多峰的情況。因此,雖然小波分析能很自然地進(jìn)行時間-尺度分析,但在利用其進(jìn)行時間-頻率分析時,必須注意尺度與頻率之間的聯(lián)系。最近提出的一種時間-頻率分析工具,即Stockwell變換ST(StockwellTransform)可以視為WT與STFT的結(jié)合。其優(yōu)點在于:(1)提供了類似于WT的多尺度的恒Q(ConstantQuality)分析;(2)同時提供了STFT才具有的絕對相位信息。近年來,S變換已廣泛應(yīng)用于地球物理,地震和醫(yī)學(xué)信號處理,電子和機械工程以及頻譜分析。關(guān)于S變換的研究,Gibson等分析了S變換與小波變換之間的區(qū)別及聯(lián)系;Pinnegar等對S變換提出了一些改進(jìn);Wong等提出了廣義S變換并對S變換的相位信息進(jìn)行了進(jìn)一步考察。在非平穩(wěn)隨機過程的演變功率譜估計方面,Iyama等和Basu等建立了不同尺度小波變換均方值與隨機過程時變譜之間的關(guān)系。Spanos等得到了諧和小波變換與非平穩(wěn)隨機過程演變功率譜密度之間的聯(lián)系;與此同時,Spanos等結(jié)合Priestley非平穩(wěn)隨機過程,利用一般非正交小波估計了非平穩(wěn)隨機過程的演變功率密度,此文的根本在于,相對于一般小波基函數(shù),Priestley非平穩(wěn)隨機過程的調(diào)制函數(shù)為相對慢變函數(shù)。2009年,Huang等將這種方法推廣到多變量非平穩(wěn)隨機過程的情況。受文獻(xiàn)啟發(fā),本文提出了一種利用S變換估計Priestley非平穩(wěn)隨機過程演變功率譜密度的方法。方法的根本在于,相對于S變換的“變換核”,Priestley非平穩(wěn)隨機過程的調(diào)制函數(shù)為相對慢變函數(shù)。因此,非平穩(wěn)隨機過程的S變換可視為相位修正后的另一非平穩(wěn)隨機過程。進(jìn)而推導(dǎo)出了以非平穩(wěn)隨機過程演變功率譜密度表達(dá)的在特定頻率點處的S變換瞬時均方值?？紤]所有頻率點的表達(dá)式,且通過將前者寫為有限個頻率點的級數(shù)展開,則以S變換瞬時均方值顯式地表達(dá)非平穩(wěn)隨機過程演變功率譜密度的問題,轉(zhuǎn)化為了求解一組代數(shù)方程。此代數(shù)方程組的解為級數(shù)展開中的依賴于時間的未知項,非平穩(wěn)隨機過程的演變功率譜密度隨之而得。2速度s變換與stft的關(guān)系為了與S變換的原始文獻(xiàn)保持一致,本文中所述推導(dǎo)式中的頻率除非明顯標(biāo)出,均為自然頻率。如果確定性過程為x(t)∈L2(R),則其S變換可以表達(dá)為式中f為頻率?？梢?S變換的核函數(shù)是經(jīng)過高斯函數(shù)調(diào)制后的諧和函數(shù),即與小波變換的相同之處在于,其高斯窗函數(shù)的寬度是隨著頻率而變化的,因此導(dǎo)致了其具有多尺度分辨率;與STFT的相同之處在于,其核函數(shù)是隨著平移因子τ而變化的,因而使其具有絕對相位信息。由于S變換可以看作為卷積形式,利用Fourier變換的卷積定理,可以得到S變換的頻域表達(dá)式為式中X(α)為確定性過程x(t)的Fourier變換?；谑?3),其離散S變換可表達(dá)為式中T為確定性過程x(t)的時間采樣間隔,N為其采樣總點數(shù),Sj[T,n/NT]為變換的采樣點。研究表明,如采用快速Fourier變換(FFT)技術(shù),能很大程度上加快S變換的計算速度。可以證明,如果Sxφ(τ,f)為確定性過程x(t)的S變換,則其平移過程x(t-r)的S變換可以表達(dá)為S變換還具有其他重要性質(zhì),限于篇幅限制,這里只給出對于估計隨機過程功率譜重要的性質(zhì),其余性質(zhì)可參見文獻(xiàn)。3子問題的估計考慮一個形如的非平穩(wěn)隨機過程,其中A(ω,t)為依賴于時間的慢變函數(shù),Z(ω)為復(fù)值正交增量的隨機過程,即式中E為期望算子,SX-X-(ω)為相應(yīng)的平穩(wěn)隨機過程的雙邊演變功率譜密度函數(shù)。因此,當(dāng)調(diào)制項A(ω,t)為時間慢變函數(shù)時,非平穩(wěn)隨機過程X(t)的演變功率譜密度函數(shù)可以近似表達(dá)為下文中常把非平穩(wěn)隨機過程的功率譜近似理論表達(dá)(式(9))稱為目標(biāo)功率譜,并與估計功率譜對比以驗證所提方法的合理性?？紤]某特定頻率fj處的S變換,并將式(6)代入式(1)可得根據(jù)Priestley非平穩(wěn)隨機過程模型(式(6))的定義,A(ω,t)相對于S變換核為時間慢變函數(shù)。具體而言,圖1所示為典型地震工程中使用的某調(diào)制函數(shù)與S變換核對比?？梢钥闯?當(dāng)S變換頻率點選定后,在相應(yīng)的S變換核有效積分區(qū)間內(nèi),調(diào)制函數(shù)可視為常值,因此可近似將此慢變函數(shù)移至積分號外,同時其時間變量取為高斯調(diào)制函數(shù)的對稱軸t=τ,即如將式(11)內(nèi)部積分單獨考慮,顯見,式(12)即為h(t)=exp[i2πω(t-τ)]的S變換。利用S變換的平移特性(式(5))可知,如果u(t)=exp(i2πωt)的S變換已知,則h(t)的S變換亦可求得。利用S變換的頻域表達(dá),u(t)的S變換可寫為進(jìn)而,h(t)的S變換可以表示為將式(14)代入式(11)可得式中保留了dZ(ω)的正交特性,亦為正交隨機過程,即若令W(τ,fj)=S(τ,fj)exp(i2πfjπ),則W(τ,fj)亦為Priestley非平穩(wěn)隨機過程。因此,原非平穩(wěn)隨機過程在某特定頻率處的S變換的相位修正可認(rèn)為是連續(xù)變量的非平穩(wěn)隨機過程,且此隨機過程的功率譜密度函數(shù)為因此,這一隨機過程的瞬時均方值為4基于時間條件的變分迭代式(18)建立了在某頻率點fj處原非平穩(wěn)隨機過程的S變換的均方值與該隨機過程演變功率譜密度的關(guān)系。考慮若干頻率點上的S變換的相位修正項,式(18)可以寫為式中N為S變換的頻率點取樣總點數(shù)。如果將原非平穩(wěn)隨機過程寫為的級數(shù)展開。其中cj(τ)為某頻率點的時變系數(shù)。結(jié)合式(20)與式(19),可得式中E[|S(τ,f)|2]為N維列向量,c(τ)為N維列向量,Q為N×N維矩陣,即式中每個元素Qi,j為由此可知,通過解矩陣方程(21),可得時變系數(shù)c(τ)。因此,可由式(20)給出的非平穩(wěn)隨機過程的演變功率譜估計。在矩陣方程(21)中,形狀函數(shù)矩陣Q是不依賴時間τ的。因此,對于所有時間取樣點,此矩陣只需求取一次逆,且由于基于FFT技術(shù)的S變換的引入,故本文建議方法具有較高的計算效率。值得特別指出的是,上述求解過程并不需要隨機過程的時-頻功率譜函數(shù)的具體形式,因此,特別適合于利用實際觀測記錄估計非平穩(wěn)隨機過程的時-頻譜。當(dāng)然,若預(yù)設(shè)某類時-頻譜形式,利用上述譜估計方法,可以通過求解反問題方式識別譜參數(shù)。有關(guān)研究工作,將另文發(fā)表。5機過程估計算法本節(jié)以兩種具有演變功率譜的非平穩(wěn)隨機過程,即均勻調(diào)制和非均勻調(diào)制的非平穩(wěn)隨機過程為例,利用本文所建議的算法,估計其演變功率譜密度。5.1平穩(wěn)過程/功率譜密度的定義考慮如下均勻調(diào)制的非平穩(wěn)隨機過程,即式中g(shù)(t)為只依賴于時間t的調(diào)制函數(shù),為平穩(wěn)隨機過程,其雙邊功率譜密度為Kanai-Tajimi譜式中S0=1.0cm2s-3rad-1為基巖白噪聲輸入強度,ζs=0.4和ωs=20rads-1分別為場地土阻尼與卓越頻率系數(shù)。調(diào)制函數(shù)g(t)定義為如圖1所示,其中a=0.25,b=0.5,k為使gmax=1的正規(guī)化系數(shù)。據(jù)式(9)可知,均勻調(diào)制非平穩(wěn)隨機過程的目標(biāo)演變功率譜密度可以表達(dá)為可以看出,不論在樣本層次上,還是在功率譜意義上,均勻調(diào)制非平穩(wěn)隨機過程的樣本/功率譜都能由相應(yīng)平穩(wěn)過程的樣本/功率譜乘以僅依賴于時間的調(diào)制函數(shù)(/的平方)而得。因此,文獻(xiàn)[20,23,24]也將調(diào)制項A(ω,t)=g(t)的均勻調(diào)制隨機過程稱為可分離(Separable)非平穩(wěn)隨機過程或“準(zhǔn)平穩(wěn)”(Quasi-stationary)隨機過程。在本文建議方法中,確定性樣本函數(shù)的生成為基于S變換的譜估計的基礎(chǔ)。在此,平穩(wěn)過程的樣本函數(shù)用譜表達(dá)方法生成,即式中Δω=ωu/n為頻率取樣間隔,ωu為上限截止頻率,φn為[0,2π]之間均勻分布的隨機變量。在譜表達(dá)方法中,樣本函數(shù)數(shù)目為500個,時間取樣點數(shù)目為1024,取樣間隔為Δt=0.02s,截止頻率取為Nyquist頻率ωu=π/Δtrad/s,諧和項數(shù)目為512。據(jù)第3節(jié)所述估計功率譜密度的方法(式(21)和式(20)),圖2所示為該非平穩(wěn)隨機過程在6s處的瞬時目標(biāo)功率譜密度與瞬時估計功率譜密度對比,可以看出,二者吻合良好。為了更清楚地表現(xiàn)在其他時刻處的瞬時功率譜密度,圖3和圖4分別給出了目標(biāo)演變功率譜密度和估計演變功率譜密度曲面圖。5.2非平穩(wěn)隨機過程的演化過程算例1為均勻調(diào)制的非平穩(wěn)隨機過程的功率譜密度函數(shù),其頻譜特征并不隨時間變化,只是幅值隨著調(diào)制函數(shù)變化。因此,這種均勻調(diào)制的隨機過程并不符合實際地震動的時變頻譜特征。正如引言中所述,由于行波效應(yīng)或傳播路徑中高頻成分不斷被吸收等因素,隨機地震動的瞬時頻譜特征是依賴于時間的。因此,利用本文建議方法估計具有演變譜隨機過程的功率譜密度函數(shù),更具有實際意義。為此,考慮具有目標(biāo)演變頻譜密度的非平穩(wěn)隨機過程為其演變功率譜密度如圖5所示。可以看出,在不同的時間段,其功率譜密度不僅在幅值上,且其卓越頻率所在區(qū)間也隨時間變化,式(29)較為充分地反應(yīng)了地震動的幅值及頻譜特點,即非平穩(wěn)不僅存在于幅值方面,而且存在于頻譜方面。由式(29)亦知其并不能表示為如式(27)所示的時間-頻率分離形式,因此,這種同時反應(yīng)時間-頻率非平穩(wěn)特征的隨機過程,在諸多文獻(xiàn)中稱之為不可分離(Non-separable)非平穩(wěn)隨機過程。非平穩(wěn)樣本函數(shù)由譜表達(dá)方法合成,即式中SXX(t,kΔω)為時變的目標(biāo)演變功率譜密度函數(shù),其他諸參數(shù)與式(28)相同。圖6為利用本文建議方法所估計的演變功率譜密度的三維視圖;圖7為在4s和12s時刻的目標(biāo)瞬時功率譜密度與估計瞬時功率譜密度對比?？梢钥闯?估計功率譜密度與目標(biāo)功率譜密度吻合良好。5.3非平穩(wěn)隨機過程隨時間的變化顯然,調(diào)制函數(shù)相對于S變換核的慢變特性對于演變功率譜估計值有較大影響。但本文所建議方法,對于大多數(shù)工程實際中幅值非快變的非平穩(wěn)隨機過程均適用。另一方面,合適地選擇S變換頻率點使其變換核較之調(diào)制函數(shù)為快變亦至關(guān)重要,如圖1所示。如何定量描述由S變換核和調(diào)制函數(shù)相對變化快慢而導(dǎo)致的估計值偏差,可由后續(xù)研究考察。其次,由式(20)可知,非平穩(wěn)隨機過程瞬時功率譜密度展開為若干高斯函數(shù)平方的加權(quán)和。參與計算的頻率點fj決定了高斯函數(shù)的對稱軸位置及寬度。數(shù)值計算表明,合理選擇頻率點fj決定了估計瞬時功率譜的光滑程度。圖2所示估計功率譜中出現(xiàn)的若干小幅抖動實為式(20)中不同高斯函數(shù)峰值;圖7中4s瞬時功率譜峰值附近處,目標(biāo)與估計功率譜的差別亦由此導(dǎo)致。最后值得注意的是,理論上,式(6)之非平穩(wěn)隨機過程的精確時變譜密度并非為式(9)之演變譜密度?？疾霵riestley非平穩(wěn)隨機過程可知,只有當(dāng)式(6)中的調(diào)制項A(ω,t)為時間慢變函數(shù)時,其時變功率譜密度才能近似以式(9)表示。因此,Priestley定義的非平穩(wěn)隨機過程為一類幅值慢變的振蕩(Oscillatory)非平穩(wěn)隨機過程。只有當(dāng)調(diào)制項A(ω,t)依時間慢變時,才不會對非平穩(wěn)隨機過程功率譜近似值(式(9))帶來時間-頻率域內(nèi)的偏移(Shift)。換言之,當(dāng)調(diào)制項A(ω,t)依時間非慢變時(其自身帶有頻率成分),須考慮由其導(dǎo)致的時間-頻率域內(nèi)的偏移。對演變功率譜分析感興趣的讀者可參見文獻(xiàn)。6非平穩(wěn)隨機過程的s變換本文首先總結(jié)和回顧了各種時間-頻率分析工具,包括Fourier變換、短時Fourier變換、WignerVille分布以及小波變換。針對最近出現(xiàn)的S變換,做了詳細(xì)描述。結(jié)合非平穩(wěn)隨機過程的Priestley模型,提出