盧 戀， 任偉新， 王世東

（1. 合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，合肥 230009；2. 深圳大學(xué) 濱海城市韌性基礎(chǔ)設(shè)施教育部重點(diǎn)實(shí)驗(yàn)室（籌），廣東 深圳 518060）

引 言

工程結(jié)構(gòu)運(yùn)營時(shí)所處的復(fù)雜環(huán)境荷載和激勵(lì)，其本質(zhì)不僅是隨機(jī)而且是非穩(wěn)態(tài)的.列車通過橋梁發(fā)生振動，結(jié)構(gòu)長期運(yùn)營發(fā)生累積損傷導(dǎo)致剛度退化，結(jié)構(gòu)在地震、強(qiáng)風(fēng)等災(zāi)變荷載作用下，在一定程度上表現(xiàn)出較強(qiáng)非線性行為等，均會導(dǎo)致結(jié)構(gòu)特性參數(shù)隨時(shí)間變化，因而，識別結(jié)構(gòu)的時(shí)變和非線性參數(shù)更符合實(shí)際情況，這不僅是工程實(shí)際和結(jié)構(gòu)健康監(jiān)測的重大需求，同時(shí)也是工程結(jié)構(gòu)參數(shù)識別理論發(fā)展的趨勢[1].

工程結(jié)構(gòu)的時(shí)變參數(shù)識別是近年來的一項(xiàng)前沿研究課題，以往基于線性系統(tǒng)和動力響應(yīng)信號平穩(wěn)性假定的分析方法，如傳統(tǒng)的Fourier變換（FT）是一種全局變換，其頻域上不包含時(shí)間信息，無法反映信號頻率隨時(shí)間變化的趨勢.時(shí)頻分析方法能同時(shí)在時(shí)域和頻域內(nèi)分析信號的時(shí)頻特征，能更好地反映信號的本質(zhì)特性，是分析非平穩(wěn)信號的有力工具.其中包括：Hilbert變換（HT）[2-3]、Wigner變換[4]、短時(shí)Fourier變換（STFT）[5]、連續(xù)小波變換（CWT）[6]、分?jǐn)?shù)階Fourier變換（FRFT）[7]等.在Hilbert變換中，F(xiàn)eldman[8]提出了一種新的非平穩(wěn)信號分解公式——Hilbert振動分解方程（HVD），實(shí)現(xiàn)了對待測信號進(jìn)行適當(dāng)?shù)姆纸猓⑹狗蛛x出的信號具有良好的Hilbert變換特性.為了研究Wigner變換在實(shí)際工程中的適用性，續(xù)秀忠等[9]利用Wigner變換和短時(shí)Fourier變換兩種方法識別時(shí)變剛度系統(tǒng)的模態(tài)參數(shù)，可以得到較為準(zhǔn)確的時(shí)間和頻率分辨率的識別結(jié)果.從理論上講，Wigner變換為非平穩(wěn)單分量線性調(diào)頻（FM）信號提供了很好的時(shí)頻分布，但是在多分量和非線性調(diào)頻單分量信號的情況下引入了交叉項(xiàng).為了克服這一缺點(diǎn)，提出了基于Fourier-Bessel級數(shù)展開的方法[10].2012年，Chen和Wang[11]提出了一種新的信號分解方法——解析模態(tài)分解（AMD），通過對該方法進(jìn)行擴(kuò)展，可以提取非平穩(wěn)信號的振動特征.

在短時(shí)Fourier變換中，利用間隔較短的時(shí)間窗函數(shù)將信號分為平穩(wěn)的若干小段，并對所有時(shí)間間隔進(jìn)行Fourier變換得到窗內(nèi)信號的頻率成分，隨著時(shí)間窗函數(shù)在信號的時(shí)間軸上滑動，最終可以獲得信號頻率與時(shí)間的對應(yīng)關(guān)系.然而，固定形狀的窗函數(shù)導(dǎo)致了時(shí)間和頻率分辨率無法同時(shí)提高，使得合理選擇窗函數(shù)和控制窗長度成為短時(shí)Fourier變換的關(guān)鍵.廖慶良[12]以一大跨度懸索橋?yàn)檠芯繉ο?，利用試?yàn)數(shù)據(jù)進(jìn)行了驗(yàn)證，探討了如何利用短時(shí)Fourier變換提高結(jié)構(gòu)參數(shù)識別的準(zhǔn)確度.連續(xù)小波變換是對短時(shí)Fourier變換的進(jìn)一步改進(jìn)，其本質(zhì)上是一種基于可調(diào)窗口的信號分析工具.與單一分析窗的短時(shí)Fourier變換相比，連續(xù)小波變換在高頻處使用短窗，低頻處使用高窗，有效克服了短時(shí)Fourier變換的分辨率限制.研究者提出了很多方法來提高小波脊線的時(shí)頻分辨率，便于識別結(jié)構(gòu)的時(shí)變參數(shù)[13-14].時(shí)變結(jié)構(gòu)的瞬時(shí)頻率的識別與小波脊線有直接關(guān)系，研究者提出了許多方法以獲取清晰的小波脊線.王超和任偉新等[15-16]基于復(fù)Morlet小波，利用小波系數(shù)模極大值方法提取小波脊線，從而識別出時(shí)變結(jié)構(gòu)的瞬時(shí)頻率.劉景良和任偉新等[17]通過最大坡度法識別出時(shí)變結(jié)構(gòu)的瞬時(shí)頻率，其識別精度高于小波系數(shù)模極大值方法，并通過數(shù)值模擬和試驗(yàn)證明了該方法對瞬時(shí)頻率提取的正確性.

從以傳統(tǒng)Fourier變換為基礎(chǔ)的時(shí)頻分析方法發(fā)展歷程可知，小波變換改進(jìn)了短時(shí)Fourier變換的分辨率問題，但在信號處理過程中也面臨著基函數(shù)合理選擇的難題.與此同時(shí)，分?jǐn)?shù)階Fourier變換同樣作為傳統(tǒng)Fourier變換的推廣，相比較于小波變換的雙變量，僅有單一變量的分?jǐn)?shù)階Fourier變換引起了學(xué)者們的關(guān)注，并且分?jǐn)?shù)階Fourier變換的核函數(shù)為chirp函數(shù)，避免了類似小波變換在運(yùn)用過程中基函數(shù)選取的問題.分?jǐn)?shù)階Fourier變換這一概念最早由Wiener在1929年提出，國內(nèi)陶然、鄧兵等[18-19]對分?jǐn)?shù)階Fourier變換進(jìn)行了系列研究.分?jǐn)?shù)階Fourier變換是由Fourier變換推廣而來的一類廣義線性變換，它可以被認(rèn)為是n次的Fourier變換，n不需要是整數(shù).因此，可以把一個(gè)函數(shù)變換到時(shí)間和頻率之間的任意中間域，為分?jǐn)?shù)階Fourier域（FRFD），其同時(shí)反映了信號在時(shí)域和頻域的信息.與常用二次型時(shí)頻分布不同的是，用單一分量（旋轉(zhuǎn)參數(shù)α）來表示時(shí)頻信息[20]，沒有交叉項(xiàng)困擾，在信號分析領(lǐng)域有更強(qiáng)的實(shí)用性.

分?jǐn)?shù)階Fourier變換作為一種新型的信號處理工具，其發(fā)展是建立在Fourier變換與短時(shí)Fourier變換的基礎(chǔ)上，所以其原理有諸多相似性.本文首先從理論分析出發(fā)，考慮到連續(xù)分?jǐn)?shù)階Fourier變換的采樣型離散分?jǐn)?shù)階Fourier變換與快速Fourier變換的計(jì)算效率相同，充分利用了Fourier變換計(jì)算效率高的優(yōu)勢，接著借鑒改變短時(shí)Fourier變換窗函數(shù)的思路，由分?jǐn)?shù)階Fourier變換定義式推導(dǎo)出一般信號的頻率與分?jǐn)?shù)階Fourier變換單一變量α的關(guān)系式，得到非平穩(wěn)信號瞬時(shí)頻率的理論表達(dá)式，實(shí)現(xiàn)一般時(shí)變信號的時(shí)頻分析，從而達(dá)到時(shí)變結(jié)構(gòu)瞬時(shí)頻率的識別的目的，通過非線性調(diào)頻信號和時(shí)變系統(tǒng)數(shù)值模擬算例驗(yàn)證方法的可行性和有效性.

1 分?jǐn)?shù)階Fourier變換基本特性

1.1 分?jǐn)?shù)階Fourier變換

相應(yīng)的逆變換為

信號x(t)的Fourier變換及對應(yīng)的逆變換，使用算子符號的表示如下：

分?jǐn)?shù)階Fourier變換作為Fourier變換的推廣，在時(shí)頻平面中，分?jǐn)?shù)階Fourier變換實(shí)際上可以被視為坐標(biāo)軸的逆時(shí)針旋轉(zhuǎn)而得到分?jǐn)?shù)階Fourier域中的表示，如圖1所示，等量關(guān)系為

式中u為 分?jǐn)?shù)階Fourier軸，可以看出，傳統(tǒng)的Fourier變換可視為當(dāng)旋轉(zhuǎn)角度為時(shí)的分?jǐn)?shù)階Fourier變換.信號x(t)的分?jǐn)?shù)階Fourier變換的定義為

Xp(u)的逆變換為

式（5）中的角度α僅在三角函數(shù)中出現(xiàn)，因此只考慮范圍 α∈(-π,π]（p ∈(-2,2]）.當(dāng) α= 0時(shí)， X0(u) = x(t)，為原信號；當(dāng)時(shí)， X1(u) =為傳統(tǒng)Fourier變換；當(dāng)時(shí)， X-1(u)為信號 x(t)的逆Fourier變換.連續(xù)分?jǐn)?shù)階Fourier變換的一些重要屬性如下[21]：

② 酉性：(Fp)-1= (Fp)H；

③ 階數(shù)、疊加性： Fp1Fp2= Fp1+p2.

圖 1 時(shí)頻面旋轉(zhuǎn)Fig. 1 Rotation of the time-frequency plane

1.2 離散分?jǐn)?shù)階Fourier變換（DFRFT）

針對數(shù)字信號，DFRFT目前主要有采樣型、線型加權(quán)型以及特征分解型三種形式.其中以采樣型DFRFT的應(yīng)用最為廣泛，該算法由Ozaktas提出，與快速Fourier變換相比，其計(jì)算復(fù)雜度同為O(NlgN).在數(shù)字信號處理系統(tǒng)中，所使用的信號是從模擬信號采樣的數(shù)字信號，即離散信號，為了實(shí)現(xiàn)采樣型DFRFT的實(shí)際工程應(yīng)用，需要將離散信號量綱歸一化處理[22]，其過程如下.

假定原始信號的時(shí)域表示限定在區(qū)間 [-Δt/2,Δt/2]，頻域表示限定在區(qū)間 [-Δf/2,Δf/2]，Δt和Δf分別為信號的時(shí)寬和帶寬.為了解決量綱不同的問題，引入尺度因子S，將新的尺度化坐標(biāo)定義為

經(jīng)量綱歸一化后的新坐標(biāo)系（x,v）量綱為1，信號被限定在 [-Δx/2,Δx/2]和 [-Δv/2,Δv/2]中：

令時(shí)寬Δt＝T，頻寬Δf＝fs，其中T為觀測時(shí)間，fs為采樣頻率，則有

① 首先信號x(t)被chirp信號調(diào)制，以為采樣間隔并利用Shannon內(nèi)插公式可得

② 將式（11）代入定義式（5）化簡可得

③ 從式（12）可以看出時(shí)域變量已經(jīng)離散化，而分?jǐn)?shù)階Fourier域變量仍是連續(xù)的，類似于離散時(shí)間Fourier變換的定義，可以認(rèn)為式（12）為離散時(shí)間FRFT.下面令u ＝m/(2Δx)，將分?jǐn)?shù)階Fourier域變量離散化，代入式（12）得

進(jìn)一步將式（13）化簡得

式中－N≤m≤N.求和部分為離散卷積形式，從而可以借助快速Fourier變換進(jìn)行快速計(jì)算.

2 分?jǐn)?shù)階Fourier變換瞬時(shí)頻率識別

由分?jǐn)?shù)階Fourier變換特性可知，分?jǐn)?shù)階Fourier變換是一個(gè)時(shí)頻算子，分?jǐn)?shù)階Fourier變換域是一個(gè)時(shí)頻域，而旋轉(zhuǎn)角度α是連接信號在時(shí)（頻）域與分?jǐn)?shù)階Fourier域的橋梁，因此討論參數(shù)α在信號瞬時(shí)頻率求解過程中的作用很有必要.式（5）中，當(dāng)信號的幅值最大時(shí)，其對應(yīng)的最佳分?jǐn)?shù)域坐標(biāo)u與 旋轉(zhuǎn)角度α有一定的關(guān)系，求解過程分為如下步驟.

由式（5）可知核函數(shù) Kα(t,u)與分?jǐn)?shù)域變量u和 旋轉(zhuǎn)角度α有關(guān)，則 Kα(t,u)對α求偏導(dǎo)數(shù)為

因此可以得到

已知

將式（17）代入式（16），可得

進(jìn)一步化簡得

式（20）為信號分?jǐn)?shù)階Fourier變換系數(shù)最大時(shí)，對應(yīng)的分?jǐn)?shù)域最佳變量與旋轉(zhuǎn)角度的相互關(guān)系.為了信號瞬時(shí)頻率的求解，令為量綱歸一化后每段時(shí)間間隔t 的分?jǐn)?shù)域表示，由式（10）可得，對分?jǐn)?shù)階Fourier變換的定義式表現(xiàn)形式進(jìn)行改變，當(dāng) α≠nπ時(shí)，由式（5）得

令

將式（22）代入式（21），可得

其中， gα(t)表示可伸縮變換的窗函數(shù)，伸縮系數(shù)為 (tan α)-1/2， gα(t-)表示窗函數(shù)在時(shí)間軸上的平移變換，平移Fourier變換算法本質(zhì)上是一種普通Fourier變換結(jié)合伸縮平移窗的算法.由傳統(tǒng)Fourier變換定義可知式（23）中 w′是與信號瞬時(shí)頻率有關(guān)的變量，且一般信號的瞬時(shí)頻率表達(dá)式為

式中，f 表示信號的瞬時(shí)頻率，u表 示信號分?jǐn)?shù)階Fourier變換幅值最大時(shí)對應(yīng)的最佳分?jǐn)?shù)域坐標(biāo).因此可得信號 x(t)的分?jǐn)?shù)階Fourier變換的表達(dá)式為

其中

經(jīng)上述分析，得到基于分?jǐn)?shù)階Fourier變換算法對一般任意信號瞬時(shí)頻率的識別步驟如下：

① 將一般信號通過式（4）旋轉(zhuǎn)到分?jǐn)?shù)階Fourier域內(nèi)；

② 由式（10）可以得到信號的分?jǐn)?shù)階域和時(shí)域表示的相互關(guān)系；

③ 通過式（19）得出信號的旋轉(zhuǎn)角度和最佳分?jǐn)?shù)域坐標(biāo)的相互關(guān)系；

④ 由式（26）得到信號最終的時(shí)頻分布.

因此，結(jié)合式（19）、（25）可知分?jǐn)?shù)階Fourier變換識別信號頻率過程中參數(shù)間的理論關(guān)系為

其識別信號作用原理如圖2所示.

圖 2 FRFT識別一般信號的原理示意圖Fig. 2 Schematic diagram of general signal identification based on the FRFT

任意非平穩(wěn)信號經(jīng)分?jǐn)?shù)階Fourier變換后，頻率f 關(guān)于旋轉(zhuǎn)角度及時(shí)間t的三維曲面表示如圖3所示，其時(shí)頻變化曲線，即圖3的投影平面如圖4所示.

圖 3 任意非平穩(wěn)信號關(guān)于頻率的三維曲面圖Fig. 3 The 3D surface of any non-stationary signal with respect to the frequency

圖 4 任意非平穩(wěn)信號時(shí)頻曲線Fig. 4 The time-frequency distribution of arbitrary nonstationary signals

3 數(shù)值驗(yàn)證分析

3.1 非線性調(diào)頻模擬信號

假設(shè)一非線性調(diào)頻模擬信號，信號表達(dá)式為

參數(shù)取值分別為：振幅 A = 1.5 m，初始頻率 f0= 15 Hz ，參數(shù) K1= 3.5， K2= 1，觀測時(shí)間為 [0,10] s，采樣頻率為fs= 800 Hz. 圖5為分別加了噪聲（信噪比（ SNR））-10 dB，10 dB，20 dB的非線性調(diào)頻模擬信號的時(shí)域波形.

圖 5 加噪信號時(shí)域波形：（a）SNR為-10 dB；（b）SNR為10 dB；（c）SNR為20 dBFig. 5 Time domain waveforms of noised signals: (a) SNR is -10 dB; (b) SNR is 10 dB; (c) SNR is 20 dB

圖 6 FRFT對信號瞬時(shí)頻率的識別過程：（a）α與最佳關(guān)系曲線；（b）時(shí)間與α關(guān)系曲線；（c）時(shí)間與最佳關(guān)系曲線Fig. 6 The FRFT identification process of signal instantaneous frequencies: (a) the relationship curve between α and; (b) the relationship curve between t and α;(c) the relationship curve between t and

圖 7 FRFT與STFT對信號瞬時(shí)頻率的識別與理論值的對比：（a）SNR為-10 dB； （b）SNR為10 dB； （c）SNR為20 dBFig. 7 The comparison between FRFT and STFT identification values of signal instantaneous frequencies and theoretical values: (a) SNR is -10 dB;(b) SNR is 10 dB; (c) SNR is 20 dB

圖6為分?jǐn)?shù)階Fourier變換識別瞬時(shí)頻率過程，分別展示了在不同信噪比條件下，與信號的瞬時(shí)頻率有關(guān)的各個(gè)參數(shù)間的相互關(guān)系.其中，為了凸顯本文提出的分?jǐn)?shù)階Fourier變換算法的優(yōu)越性，得到基于短時(shí)Fourier變換與分?jǐn)?shù)階Fourier變換的信號瞬時(shí)頻率識別結(jié)果與理論值對比圖，如圖7所示.從圖中可知，在不同的信噪比條件下，F(xiàn)RFT（虛線）對比于STFT（實(shí)線）而言，更接近于信號的理論值（點(diǎn)線），表明了分?jǐn)?shù)階Fourier變換能更好地識別出模擬的非線性調(diào)頻信號的頻率變化規(guī)律，結(jié)果分析可知，本文方法識別出的瞬時(shí)頻率與理論值吻合良好，且方法具有一定的抗噪能力，相對比于短時(shí)Fourier變換而言在識別精度上有一定的優(yōu)越性.

3.2 時(shí)變結(jié)構(gòu)數(shù)值算例

如圖8所示的三自由度質(zhì)量-彈簧-阻尼線性時(shí)變系統(tǒng)， m1= m2= m3= 1 kg， k1= k3= 500 N/m，c1= c2= c3=c4= 0.15 N·s·m-1，假定剛度 k2，k4隨時(shí)間線性變化，具體變化如下式所示：

根據(jù)“凍結(jié)法”假定結(jié)構(gòu)參數(shù)在小的時(shí)間間隔內(nèi)保持不變，為時(shí)不變系統(tǒng)，可計(jì)算出不同時(shí)刻系統(tǒng)的三階固有頻率，如圖9所示，作為系統(tǒng)頻率的理論值.

圖 8 三自由度線性系統(tǒng)Fig. 8 The 3DOF linear system

圖 9 剛度線性變化時(shí)結(jié)構(gòu)固有頻率Fig. 9 Natural frequencies of the structure with linearly changing stiffness

為模擬出圖8所示系統(tǒng)的響應(yīng)信號，假定系統(tǒng)初始條件 x1(0) = 0.3 m， x2(0) = -0.06 m， x3(0) = 0.01 m，為提高計(jì)算效率，采用自適應(yīng)變步長的四階Runge-Kutta法求解結(jié)構(gòu)的位移響應(yīng)，即使用較大的步長求解變化較慢的位移響應(yīng)以提高計(jì)算速度，反之使用較小的步長求解變化較快的位移響應(yīng)以提高計(jì)算精度.由于計(jì)算過程中時(shí)間步長的不同，可先采用樣條插值法擬合得到位移響應(yīng)，再以200 Hz的采樣頻率得到等步長的結(jié)構(gòu)位移自由響應(yīng)，采樣時(shí)間為20 s.

為模擬實(shí)際噪聲的影響，對離散的響應(yīng)信號添加信噪比為20 dB的Gauss白噪聲，加噪后質(zhì)量1的位移響應(yīng)如圖10所示.對加噪后的位移響應(yīng)信號進(jìn)行分?jǐn)?shù)階Fourier變換，識別出結(jié)構(gòu)的三階頻率與其對應(yīng)的理論值對比圖分別如圖11~13所示.

結(jié)果表明，本文建立的基于分?jǐn)?shù)階Fourier變換的瞬時(shí)頻率識別方法，用時(shí)變結(jié)構(gòu)位移響應(yīng)識別出的瞬時(shí)頻率結(jié)果和理論值吻合良好，且該方法具有一定的抗噪性能.

圖 10 質(zhì)量1的位移響應(yīng)Fig. 10 The displacement response of mass 1

圖 11 FRFT識別及理論結(jié)構(gòu)第一階固有頻率Fig. 11 The FRFT-based identification and the 1st-order theoretical natural

圖 12 FRFT識別及理論結(jié)構(gòu)第二階固有頻率Fig. 12 The FRFT-based identification and the 2nd-order theoretical natural frequencies of the structure

圖 13 FRFT識別及理論結(jié)構(gòu)第三階固有頻率Fig. 13 The FRFT-based identification and the 3rd-order theoretical natural frequencies of the structure

4 結(jié) 論

與傳統(tǒng)的Fourier變換相比，分?jǐn)?shù)階Fourier變換增加了一個(gè)旋轉(zhuǎn)角度α參數(shù)，適用于非平穩(wěn)信號的處理，具有較好的時(shí)頻分辨率和局部化特性.本文基于變量α對于分?jǐn)?shù)階Fourier變換分析信號核心作用進(jìn)行分析，從理論上解釋了分?jǐn)?shù)階Fourier變換本質(zhì)上是一種普通Fourier變換結(jié)合伸縮平移窗的算法，進(jìn)而在分?jǐn)?shù)階Fourier域建立了非平穩(wěn)信號瞬時(shí)頻率的一般表達(dá)式，實(shí)現(xiàn)了結(jié)構(gòu)瞬時(shí)頻率的識別.采用任意非線性調(diào)頻信號仿真算例和三自由度有阻尼時(shí)變結(jié)構(gòu)系統(tǒng)的數(shù)值算例對提出的方法進(jìn)行了驗(yàn)證，得到如下結(jié)論：

1）分?jǐn)?shù)階Fourier變換在信號處理中的作用相當(dāng)于普通Fourier變換在分析信號時(shí)結(jié)合了伸縮平移窗函數(shù)，隨著窗函數(shù)的平移以及彈性伸縮進(jìn)而得到信號的時(shí)頻分布.其中起著至關(guān)重要的作用是旋轉(zhuǎn)參數(shù)p(α= pπ/2)，首先通過旋轉(zhuǎn)參數(shù)旋轉(zhuǎn)坐標(biāo)軸將時(shí)域內(nèi)信號變換到分?jǐn)?shù)階域內(nèi)，繼而在分?jǐn)?shù)階Fourier域中，對加窗信號進(jìn)行常規(guī)Fourier計(jì)算，窗函數(shù)是由旋轉(zhuǎn)參數(shù)的正切函數(shù)進(jìn)行伸縮變換的chirp函數(shù).

2）推導(dǎo)出了一般信號瞬時(shí)頻率的表達(dá)式，論證了信號瞬時(shí)頻率和各個(gè)因素的關(guān)系.通過模擬和數(shù)值算例表明該方法能有一定的抗噪性，識別出的瞬時(shí)頻率值與理論值能較好吻合，在時(shí)變結(jié)構(gòu)的模態(tài)參數(shù)識別方面具有一定的應(yīng)用前景.