孔 凡，李 杰，2

(1.同濟(jì)大學(xué) 土木工程學(xué)院，上海 200092;2.同濟(jì)大學(xué) 土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092)

小波分析發(fā)展歷程充分體現(xiàn)出工程與數(shù)學(xué)領(lǐng)域相互促進(jìn)關(guān)系。小波(wavelet)一詞由Morlet等［1］對地震數(shù)據(jù)分析時提出。Meyer等［2-4］提出連續(xù)小波容許性條件及重構(gòu)公式，并將計算機(jī)視覺領(lǐng)域中多尺度分析思想引入小波分析中，提出了多分辨率分析概念及Mallat算法。Daubechies［5］創(chuàng)立了支持離散小波的二進(jìn)小波理論，并構(gòu)造出緊支的正交小波基。在小波發(fā)展史中，Grossman等［6-7］也做出了杰出貢獻(xiàn)。

在完善和發(fā)展小波的純數(shù)學(xué)理論同時，小波理論在應(yīng)用數(shù)學(xué)及工程領(lǐng)域中也頗具前景。Newland等［8］首先將小波分析應(yīng)用于工程振動信號分析，并提出了諧和小波(Harmonic Wavelet，HW)及廣義諧和小波(Generalized Harmonic Wavelet，GHW)概念。HW 與GHW為一類具有解析形式，在頻率域上緊支但時域衰減較慢(依t-1)的小波［8-14］。HW被應(yīng)用在求解常/偏微分方程［15-17］、積分與積分 - 微分方程［18-19］、時頻分析［20］、估計隨機(jī)過程功率譜密度［21-22］以及計算隨機(jī)激勵下系統(tǒng)二階響應(yīng)［23-24］等方面。

諧和小波良好的工程應(yīng)用前景在于其特殊的諧和成分與實(shí)際中發(fā)生的振動信號結(jié)構(gòu)有相似之處［19］。快速Fourier變換(FFT)技術(shù)的應(yīng)用，使諧和小波變換計算效率更高。HW與GHW均在無窮時域上正交，工程實(shí)際中發(fā)生的振動信號卻不具有無窮持時特征。因此，基于GHW概念，本文提出以有限時間T0為基本周期的周期廣義諧和小波(Periodic Generalized Harmonic Wavelet，PGHW);證明了此類小波在有限區(qū)間內(nèi)的正交性;提出了基于FFT的PGHW快速冗余/非冗余小波變換及重構(gòu)算法，指出了冗余和非冗余小波變換間關(guān)系;利用所提小波變換算法，對某地震動時程進(jìn)行周期廣義諧和小波變換及其重構(gòu)。

1 廣義諧和小波

HW與GHW均為具有解析表達(dá)、頻域緊支且時域衰減很慢的正交復(fù)小波［8-10］。其中母小波及尺度函數(shù)表達(dá)為:

式中:i為虛數(shù)單位(除說明顯外，正體i為虛數(shù)單位，斜體i為下標(biāo)符號);t為時間。經(jīng)平移與尺度變換后的HW為:

式中:j為尺度因子;k為平移變換因子，k=2j，2j+1，…，22j-1。

經(jīng)尺度與平移變換后，HW的Fourier變換為:

式中:Ψ(ω)為母諧和小波的Fourier變換。

尺度函數(shù)、母小波函數(shù)及經(jīng)平移與尺度變換后的小波Fourier變換之模見圖1。

圖1 各級諧和小波的Fourier變換之模Fig.1 Squared modulus of the Fourier transform of HW at different scales

由圖1知，諧和小波為二進(jìn)小波，各級小波頻域互不交疊、且逐級寬度加倍高度減倍。如尺度函數(shù)頻域之模分布于［0，2π)上;0級小波頻域之模分布于［2π，4π)上;1級小波頻域之模分布于［4π，8π)上;第j級小波頻域之模分布于［2j2π，2j+12π)上。每級小波頻域?qū)挾葹槎ㄖ?，即第j級小波頻域?qū)挾葹?j2π。因此不具靈活性。為此 Newland［9］提出音樂小波(Musical Wavelet)，也稱為GHW。各級GHW的Fourier頻寬并不固定為2j2π，而由一對尺度指標(biāo)確定，即(mj-nj)2π。其中mj，nj為第j級尺度指標(biāo)。在Newland的原文中，對信號長度進(jìn)行正規(guī)化，因此無論HW或GHW，其各級小波頻寬均為2π的倍數(shù)(信號長度為T0=1 s)。本文中對GHW采用一般表示方法，即:

式中:i為尺度下標(biāo);k=0，1，…，n-m-1為平移變換因子;G表示廣義;Δω=2π/T0為信號頻域采樣間隔;niΔω與miΔω為第i級尺度下小波頻域上下限。特別地，可考慮各尺度下小波頻帶寬相等，即ni-mi=njmj=n-m且ni-1=mi;T0為考慮信號持時。式(5)的GHW在頻域內(nèi)表達(dá)為:

各階GHW在頻域內(nèi)分布見圖2。

圖2 各階GHW的Fourier變換模Fig.2 Squared modulus of the Fourier transform of GHW at different scales

2 GHW周期化

據(jù)周期化非周期小波的一般過程［5］，周期廣義諧和小波可寫為:

其中:上標(biāo)per代表“周期”。

據(jù)Fourier變換性質(zhì)，式(7)右邊可寫為Fourier變換的逆變換，即:

式(8)右邊定義域?yàn)閙Δω≤ω＜nΔω。據(jù)Poisson公式［10，18-19］，有:

式中:δ(ω-qΔω)為函數(shù)。當(dāng)l=π/T0時，式(9)可化為:

結(jié)合式(8)與式(10)知:

因此，PGHW可寫為:

容易證明，式(12)頻域表達(dá)為:

由式(12)知，周期廣義諧和小波在時域上表現(xiàn)為若干諧和項(xiàng)之和;且平移因子表現(xiàn)為諧和項(xiàng)相位滯后;在頻域上表現(xiàn)為若干δ函數(shù)之和，且其中每個δ函數(shù)坐落于此尺度的頻域區(qū)間內(nèi)每個頻率采樣點(diǎn)處。與GHW相同，PGHW為復(fù)小波，其實(shí)部為偶函數(shù)虛部為奇函數(shù)。式(12)中，平移變換因子k=0，1，…，n-m-1，PGHW平移一次，實(shí)際在時域上平移了T0/(n-m)?？梢宰C明，平移后的PGHW在基本周期［0，T0)內(nèi)是正交的。

如果使PGHW每平移一次，在時域上表現(xiàn)為平移了Δt，Δt為時間取樣間隔。則式(12)可改寫為:

式中:k=0，1，…，N-1為平移變換因子，N為時間取樣總數(shù)。式(14)可稱為冗余PGHW。可證，平移后的PGHW在基本周期上并非兩兩正交。

圖3為在某尺度與平移因子下且基本周期為T0=20.48 s的正交PGHW。由圖3知，PGHW的實(shí)部與虛部分別為關(guān)于其基本周期內(nèi)中點(diǎn)的偶函數(shù)與奇函數(shù)。尺度因子越大，振動頻率越高。

圖3 正交(t)，(t)，(t)(t)小波的實(shí)部(實(shí)線)與虛部Fig.3 Real and imaginary part of orthogonal PGHW，(t)(t)(t)(t)

據(jù)PGHW頻域簡單性，易證式(12)PGHW在其基本周期［0，T0)上是正交的，即:

其中:〈·〉為函數(shù)內(nèi)積符號，且:

由于PGHW的正交性，可保證其完全重構(gòu)性(Perfect Reconstruction)。

3 變換與重構(gòu)

根據(jù)定義，正交PGHW的小波變換可寫為:

結(jié)合式(12)，式(17)可寫為:

據(jù)式(18)，正交PGHW變換見圖4。

圖4 f(t)正交PGHW的小波變換Fig.4 Flow chart of the fast orthogonal PGHW transform

當(dāng)原始信號為復(fù)值信號時，其重構(gòu)公式為:

以工程教育認(rèn)證為指導(dǎo)綱要,在培養(yǎng)社會所需人才的模式和方法上,突出工程應(yīng)用教育的特點(diǎn),強(qiáng)調(diào)“三結(jié)合”(理論與實(shí)踐結(jié)合、課內(nèi)與課外結(jié)合、學(xué)校與企業(yè)結(jié)合),堅(jiān)持系統(tǒng)更新優(yōu)化、穩(wěn)中求新、求進(jìn)的原則,考慮社會發(fā)展、技術(shù)革新、需求引導(dǎo)、市場反饋意見等多方面因素,從而形成持續(xù)改進(jìn)的測控技術(shù)與儀器專業(yè)人才培養(yǎng)模式。

式中:Re表示求實(shí)部。

數(shù)值試驗(yàn)表明，利用式(21)、式(12)進(jìn)行PGHW小波逆變換，由于循環(huán)求和原因，計算效率不佳。利用FFT技術(shù)，可加快PGHW小波變換速度。考察PGHW周期性后，本文提出快速PGHW逆變換算法。式(19)可表達(dá)為:

式中:N為時間取樣總點(diǎn)數(shù)。

式中:t1，t2，…，tN為時間取樣點(diǎn)，)，k為的復(fù)共軛。

同理，式(14)的非正交小波可變換為:

圖5 快速正交PGHW重構(gòu)算法Fig.5 Flow chart of the fast reconstruction algorithm of orthogonal PGHW

將式(14)代入式(28)得:

非正交PGHW的小波變換見圖6。

圖6 f(t)非正交PGHW的小波變換Fig.6 Flow chart of the fast redundant PGHW transform

非正交PGHW的小波逆變換也可表達(dá)為式(19)、式(21)形式，區(qū)別在于該小波變換表達(dá)式為式(28)且:

式中:k=0，1，…，N-1為時間平移因子。

雖然快速非正交PGHW的小波重構(gòu)算法與快速正交PGHW的小波重構(gòu)算法類似，但區(qū)別為無圖5中上插值(Upsampling)環(huán)節(jié)。

顯然，式(12)小波在基本周期上構(gòu)成一完備正交系，小波之間兩兩正交，稱其小波變換為非冗余小波變換。式(14)小波集合包含式(12)小波集合，稱其小波變換為冗余小波變換。非冗余小波變換為與離散信號同維向量;而冗余小波變換為矩陣，其中每個尺度下冗余小波變換為矩陣中的某一向量。即，非冗余小波變換與原始信號之間所包含的信息相等，當(dāng)且僅當(dāng)有該非冗余小波變換時，即可完全重構(gòu)原始信號;而冗余小波變換中所包含信息大于原始信號中信息。因此，非冗余小波變換較冗余小波變換更高效，但重構(gòu)原始信號時，前者較后者對誤差更敏感。由于冗余小波的光滑性(圖7)，在利用PGHW估計隨機(jī)過程功率譜應(yīng)用中，較正交非冗余PGHW更具優(yōu)勢。對比式(18)與式(28)，二者之關(guān)系為:

4 數(shù)值算例

以非平穩(wěn)地震動時程為例，利用所提非冗余及冗余PGHW變換的快速算法，進(jìn)行時程樣本的小波變換及重構(gòu)。平穩(wěn)隨機(jī)地震動功率譜采用Kanai-Tajimi譜，即:

式中:ωs，ζs為場地卓越頻率及等效阻尼;S0為基巖輸入白噪聲幅值。

地震動非平穩(wěn)性由包絡(luò)函數(shù)表達(dá)，即:

式中:g(t)為包絡(luò)函數(shù);S0分別為平穩(wěn)隨機(jī)過程與非平穩(wěn)隨機(jī)過程樣本函數(shù);α=0.25，β=0.5為包絡(luò)函數(shù)參數(shù);k為使gmax(t)=1的正規(guī)化系數(shù)。

目標(biāo)非平穩(wěn)地震動時程樣本由式(32)及平穩(wěn)譜表達(dá)方法［25-26］(Spectral Representation Method)生成。

圖7為樣本地震動在某尺度下的非冗余及冗余小波變換，二者之關(guān)系見式(31)。圖8為目標(biāo)非平穩(wěn)地震動時程與非冗余PGHW快速小波重構(gòu)后的非平穩(wěn)時程對比。利用冗余PGHW快速小波重構(gòu)算法也可得到類似結(jié)果。由圖8可見，小波變換后重構(gòu)信號與原信號吻合良好，說明本文所提PGHW快速小波變換與重構(gòu)算法的可靠性。

圖7 某樣本地震動在某尺度下非冗余及冗余小波變換Fig.7 Non-redundant and redundant PGHW transform of an earthquake excitation at a certain scale

圖8 樣本目標(biāo)地震動與重構(gòu)地震動對比Fig.8 Comparison between the original and the reconstructed signal

5 結(jié)論

(1)本文提出在基本周期［0，T0)內(nèi)正交的諧和小波，即周期廣義諧和小波(PGHW)。PGHW母小波可看作經(jīng)若干諧和分量疊加而得;不僅保留了諧和小波族在頻域內(nèi)的簡單性，且非冗余PGHW在基本周期［0，T0)上正交。

(2)本文給出基于快速Fourier變換的快速正交(非冗余)PGHW變換與冗余PGHW變換及相應(yīng)的快速逆變換，指出冗余及非冗余小波變換之間的關(guān)聯(lián)性。數(shù)值算例表明，由于PGHW在有限區(qū)間內(nèi)的正交性，保證了其在相應(yīng)區(qū)間內(nèi)的完全重構(gòu)性。

［1］ Grossmann A，Morlet J.Decomposition of Hardy function into square intergrable wavelets of constant shape［J］.Journal on Mathematical Analysis，1984，15:723-736.

［2］Meyer Y.Principe d’incertitude，bases hilbertiennes et algebres d’opérateurs［J］.Seminaire Bourbaki，1985，622:1985-1986.

［3］Meyer Y.Ondelettes et operateurs，I.Ondelettes［M］.Paris:Hermann，1990.

［4］Mallat S. Multiresolution approximations and wavelet orthonormal bases of L2(R)［J］.Transactions of the American Mathematical Society，1989，315(1):69-87.

［5］Daubechies I.Ten lectures on wavelets［M］.Philadelphia:Society for Industrial and Applied Mathematics，1992.

［6］Coifman R R.Wavelet analysis and signal processing，signal processing part I［M］.Springer-Verlag，1990:59-68.

［7］Chui C K.An introduction to wavelets［M］.San Diego:Academic Press Professional，Inc，1992.

［8］ Newland D E.An introduction to random vibrations，spectral and wavelet analysis［M］.New York:Longman Scientific＆Technical，1993.

［9］Newland D E.Harmonicandmusicalwavelets ［J］.Proceedings of the Royal Society of London.Series A:Mathematical and Physical Sciences，1994，444(1922):605-620.

［10］ Newland D E.Harmonic wavelet analysis［J］.Proceedings of the Royal Society of London.Series A:Mathematical and Physical Sciences，1993，443(1917):203-225.

［11］ Newland D E.Ridge and phase identification in the frequency analysis of transient signals by harmonic wavelets ［J］.Journal ofVibrationandAcoustics， 1999， 121(2):149-155.

［12］ Newland D E.Practical signal analysis:do wavelets make any difference?［A］.Proceedings of the 16th ASME Biennial Conference on Vibration and Noise.Sacramento，1997.

［13］ Newland D E.Wavelet analysis of vibrations，part II:wavelet maps ［J］.Journal of Vibration and Acoustics，1994，116(4):417-425.

［14］ Newland D E.Wavelet analysis of vibration，part I:theory［J］.Journal of Vibration and Acoustics，1994，116(4):409-416.

［15］Muniandy S V，Moroz I M.Galerkin modelling of the burgers equation using harmonic wavelets［J］.Physics Letters A，1997，235(4):352-356.

［16］ Cattani C.The wavelet-based technique in dispersive wave propagation ［J］.International Applied Mechanics，2003，39(4):493-501.

［17］ Cattani C. Harmonic waveletstowardsthe solution of nonlinear PDE ［J］.Computers and Mathematics with Applications，2005，50(8-9):1191-1210.

［18］ Cattani C，Kudreyko A.Harmonic wavelet method towards solution of the fredholm type integral equations of the second kind ［J］.Applied Mathematics and Computation，2010，215(12):4164-4171.

［19］Cattani C，Kudreyko A.Application of periodized harmonic wavelets towards solution of eigenvalue problems for integral equations ［J］.Mathematical Problems in Engineering，2010，2010:570136-570143.

［20］ Spanos P，Giaralis A，Politis N.Time-frequency representation of earthquake accelerograms and inelastic structural response recordsusing the adaptive chirpletdecomposition and empirical mode decomposition ［J］.Soil Dynamics and Earthquake Engineering，2007，27(7):675-689.

［21］ Spanos P D，Tezcan J，Tratskas P.Stochastic processes evolutionary spectrum estimation via harmonic wavelets［J］.Computer Methods in Applied Mechanics and Engineering，2005，194(12-16):1367-1383.

［22］ Tratskas P. Wavelet-based excitation Representation and response determination of linear and nonlinear systems［D］.Houston:Rice University，2001.

［23］Spanos P D，Kougioumtzoglou I A.Harmonic wavelets based statisticallinearization for response evolutionary power spectrum determination ［J］. Probabilistic Engineering Mechanics，2012，27(1):57-68.

［24］ Tratskas P， Spanos P D. Linear multi-degree-offreedomsystem stochastic response by using the harmonic wavelet transform［J］.Journal of Applied Mechanics，2003，70(5):724-731.

［25］Liang J W，Chaudhuri S R，Shinozuka M.Simulation of nonostationary stochastic process by spectral representation［J］.Journal of Engineering Mechanics，2007，133(6):616-627.

［26］Shinozuka M.Simulation of stochastic processes by spectral representation［J］.Applied Mechanics Reviews，1991，44(4):191-204.