淺談氣動阻尼識別方法

饒明航, 谷任奇, 李紅義

(西南交通大學，四川成都 610031)

流體對結構振動造成阻礙并消耗其能量，這一影響結構振動的氣動力參數稱為氣動阻尼[1]。氣動阻尼受結構外形、風場特性、風速等很多因素影響，產生機理十分復雜。氣動阻尼的理論解析和數值風洞研究都有一定難度，而風洞試驗是目前研究氣動阻尼的常用手段。風洞試驗研究的關鍵問題之一就是如何從試驗數據中準確識別氣動阻尼，傳統的阻尼識別方法主要分為兩類：譜分析方法、時間序列法[2]。其中時間序列法里的隨機減量法在氣動阻尼識別上應用較多，但是對于多自由度系統，隨機減量法會混有多階模態和噪聲信號。因此，隨機減量法常與其他方法相結合來識別氣動阻尼，其中Hilbert-Huang變換(Hilbert-Huang transform,簡稱HHT)和小波變換是研究最多的兩種氣動阻尼識別方法。另外，大多的阻尼識別方法包括隨機減量法，都無法確定阻尼比的識別精度，而譜分析方法中的快速貝葉斯FFT則可以通過阻尼的后驗變異系數確認識別精度[3]。

因此，本文首先對譜分析方法和時間序列法這兩類方法進行了介紹和分析；其次，對以隨機減量法為基礎的HHT和小波分析進行了詳細介紹，并進行分析和比較；最后，對譜分析方法中快速貝葉斯FFT方法進行了詳細闡述和分析。

1 氣動阻尼識別方法

氣動阻尼是通過識別結構自身阻尼和總阻尼，再由總阻尼減去自身阻尼得到的。傳統的阻尼識別方法主要分為譜分析方法和時間序列法[2]。譜分析方法包括：譜距法、譜曲線擬合法、半功率帶寬法、對數遞減法、快速貝葉斯FFT等。時間序列法又包括：最大熵法、自回歸法、自相關衰減法、特征系統實現法、隨機子空間法以及隨機減量法等。

1.1 譜分析方法

對有限時長信號進行傅里葉變換的譜分析方法必然會有誤差。選擇較長的數據是減小譜分析方法誤差的有效辦法，但是這又會降低精度。因此，需要通過增加自由度、合并相鄰的譜估計、對幾個譜估計進行整體平均等方法來提高譜分析方法的精度。由于較長數據的不穩定性，運用這些方法時通常需要通過曲線擬合來保證譜估計值的平滑性，如Breukelman[4]和Jones等[5]分別采用最大似然估計和最小二乘法對單自由度系統譜估計進行了擬合。但是Jeary[6]指出，曲線擬合適用于線性系統，卻無法反映非線性系統的特征。

1.2 時間序列法

時間序列法相比較于譜分析方法，不存在短數據的分辨率不足問題和長數據的不穩定性問題。但是文獻[7-9]指出，這一方法對模態階數十分敏感，以自回歸模型為例，模態階數太少可能得不到主峰值，而模態階數太多又會產生偽峰值。因此，選擇合適的模態階數才能通過時間序列法得到比較精確的阻尼。

時間序列法中的隨機減量法是目前評估風振結構阻尼較多的一種方法。Cole[10]在1969年提出了隨機減量法，基本思想是：從零平均隨機激勵產生的響應采集樣本，保證每個樣本的初始條件相同；再對樣本進行整體平均，消除了外界激勵引起的響應，進而得到相同初始條件的自由振動響應；最后對自由振動響應進行結構阻尼的評估。譜分析方法大都要求信號是平穩的和線性的，而隨機減量法不僅可以進行線性結構的阻尼評估，也適用于非線性結構[11]。Jery[11]采用隨機減量法對已有的幾個非線性結構風振響應的信號進行阻尼評估，得到的阻尼與公認正確的結論吻合良好。

對于多自由度系統信號，采用隨機減量法通常會混有多階模態和噪聲信號。另一方面，大多的阻尼識別方法包括隨機減量法，都無法確定阻尼比的識別精度。因此，本文將對目前較為前沿的HHT方法、小波變換以及快速貝葉斯FFT進行闡述分析，其中HHT和小波變換都是與隨機減量法相結合來識別氣動阻尼的。

2 Hilbert-Huang變換

Huang[12]在1998年提出的一種分析非線性非平穩信號、具有自適應性的處理方法，即Hilbert-Huang變換。這一方法由經驗模態分解法(Empirical mode decomposition，簡稱EMD)和希爾伯特變換(Hilbert transform，簡稱HT)兩部分組成。

2.1 經驗模態分解法

對信號進行EMD處理，可以使信號平穩化，并產生一組具有不同特征尺度的序列，即本征模函數(intrinsic mode Function，簡稱IMF)。EMD的具體步驟是：分別對信號X(t)的極大值和極小值進行三次樣條曲線擬合，形成上下包絡線。將上下包絡線進行整體平均，得到均值m1(t)，再用X(t)減去m1(t)得到新的序列h1(t)：

h1(t)=X(t)-m1(t)

(1)

重復以上過程k次，若mk(t)趨近于零，就獲得了第一個IMF分量C1(t)：

C1(t)=hk(t)=hk-1(t)-mk(t)

(2)

用原始序列X(t)減去C1(t)得到剩余項r1(t)，再對r1(t)重復上述平穩化過程，就可以獲得第二個IMF分量C2(t)，如此重復直到余項rn(t)不能再被分解為止。最后原始序列X(t)可以表示為n個IMF分量和余項：

(3)

2.2 希爾伯特變換

通過HT對信號進行處理，可以得到信號幅值和頻率隨時間變化的關系[13]。對信號δ(t)進行希爾伯特變換可得：

(4)

(5)

式中：a(t)為幅值函數，θ(t)為相位角函數：

a(t)=Ae-ξω0t

(6)

θ(t)=ωdt+φ0

(7)

對式(6)進行對數變換可得：

Ina(t)=InA-ξω0t

(8)

HHT識別氣動阻尼的具體步驟為：先對原始數據進行EMD分解；再將關心的IMF分量用隨機減量法處理，得到自由衰減曲線；最后對自由衰減曲線進行希爾伯特變換，并對結果進行最小二乘擬合求得阻尼比。

孫旭峰[14]和樓文娟等[15]采用HHT和隨機減量法相結合的方法，分別識別了肋環型索穹頂結構和大跨度柔性屋蓋結構的氣動阻尼，表明HHT適用于大跨度這種非線性結構，可以較好的分析非線性非平穩信號。但是使用HHT識別結構參數會有以下問題：(1)在低頻段產生不期望的本征模函數[16]、難以分離微弱信號或能量小的模態[17]；(2)EMD分解過程中的端點效應會不斷積累誤差，以致高階本征模函數失真[18]。為了解決第一個問題，潘峰等[19]和徐正[20]在運用HHT識別結構氣動阻尼時引入了帶通濾波。除此以外，還有引入間歇準則、加入偽信號等方法也可以解決第一個問題。對于端點效應問題，Huang提出采用特征波進行端點延拓，但是并沒有給出確定特征波的方法[18]。

3 小波分析

小波分析是一種時頻局部化方法，其窗口大小一定而形狀可以改變。若平方可積的實數域L2(R)內有信號x(t)和基小波g(t)，其中g(t)滿足小波允許條件，則x(t)經過小波變換可得小波系數W(a,t)。

(9)

式中：a為尺度因子，*表示共軛。模態參數識別通常采用復Morlet小波，主要是因為復Morlet小波實部和虛部的幅值函數與結構的自由響應信號有相同的函數形式，且復Morlet小波不存在尺度函數，可以實現任意時頻域的分辨率，對密集模態的識別效果較好[21]。復Morlet小波的時域表達形式如下：

(10)

其中fc為中心頻率，fb為帶寬頻率。對式(10)進行傅里葉變換，可得復Morlet小波頻域表達形式為：

G(af)=e-π2fb(af-fc)2-e-π2fb((af)2+fc2)

(11)

(12)

G(af)=e-π2fb(af-fc)2

(13)

對式(13)取極值可知因子a=fc/f。假設結構為有N個節點的多自由度粘性阻尼系統，其自由響應信號為：

(14)

其中Ak為第k階模態幅值；θk為相位角；fk、fdk和ξk分別為第k階的無阻尼頻率、有阻尼頻率和阻尼比，且滿足：

(15)

將式(12)、式(14)代入式(9)中可得：

(16)

由第k階尺度因子ak=fc/fk可得：

(17)

(18)

ξk=-λ1/(2πfk)

(19)

當采用小波分析識別結構氣動阻尼時，首先通過隨機減量法得到結構的自由衰減曲線，再進行小波分析，按上述方法得到阻尼比。

文獻[22-23]聯合采用隨機減量法和小波變換，對結構的模態參數進行了識別，從而成功的引入了新的時頻分析方法。吳海洋等[24]正是運用這一方法識別了一多自由度高層建筑的氣動阻尼，表明基于隨機減量法的小波分析適用于多自由度結構的氣動阻尼識別。在模態參數分析的過程中，小波分析會產生邊界效應問題。針對這一問題，吳海洋等[24]從實用角度出發，除去端部而只取中間數據進行分析，較好的識別了氣動阻尼。另外，Kijewski等[23]詳細分析了邊界效應問題，提出在端部增補數據來減少邊界效應的影響。

小波變換和HHT變換都是處理非平穩非線性信號的時頻分析方法，與隨機減量法相結合都可用于識別結構的氣動阻尼。HHT變換不需要任何的先驗知識，具有良好的自適應性，但是沒有完善數學基礎，同時存在模態混淆和端點效應問題。小波變換則具有完善的數學基礎，可以對信號進行詳細的時頻分析和多分辨率分析，但存在邊界效應問題。

4 快速貝葉斯FFT方法

Yuen和Katafygiotis[25-27]提出了基于貝葉斯快速傅里葉變換(Bayesian Fast Fourier Transform, Bayesian FFT)的模態參數識別方法。這一方法有以下幾個難點：(1)自由度和模態參數的增加導致求解目標函數時會涉及到一個病態矩陣的逆矩陣；(2)模態參數的最佳估計需要求解的數值優化問題維度過多；(3)三是需要通過有限差分法來確定Hessian矩陣，這些難點使得該方法的應用受到很大的限制[28]。針對這些難點，Au[28]提出了快速貝葉斯FFT(Fast Bayesian FFT)方法。該方法是對某一共振頻率帶上的單一模態進行分析，因此不需要考慮自由度的數量，并且模態參數的最佳估計和Hessian矩陣可以分別通過求解一個四維數值優化問題和對數似然函數的二階導數來確定。Au等[29]運用快速貝葉斯FFT方法識別了一個主-次結構的模態參數。

快速貝葉斯FFT方法將加速度響應表示為：

(20)

(21)

式中：i2=-1，Δt為釆樣時間間隔，Fk和Gk為傅里葉變換的實部及虛部，由于采樣通道的電壓偏置和傅里葉變換共軛等原因，僅取k=2,...,Nq，Nq=int[N/2]+1。

根據貝葉斯理論，后驗概率密度函數PDF可以由對數似然函數L(θ)表示：

p(θ|{Zk})∝exp[-L(θ)]

(22)

L(θ)= -nNfIn2+(n-1)NfInSe

(23)

θ為模態參數，包括自振頻率f、阻尼比ξ、振型Φ、功率譜密度S以及預測誤差的譜密度Se，Nf為計算模態相應頻段上的離散頻率值，k為對應頻帶內傅里葉變換所對應的值，并且有：

(24)

(25)

(26)

通過對式(23)的無約束優化就可以求得阻尼比，而后驗不確定性可通過后驗協方差矩陣來計算，其大小可以根據各自的后驗變異系數來進行衡量。

快速貝葉斯FFT方法已被國內外多次用來識別結構模態參數，并且該方法還可以評估參數識別結果的不確定性。黃銘楓[30]采用快速貝葉斯FFT方法對高層結構氣動阻尼進行了識別并與HHT方法識別的結果進行了對比，兩者識別出的阻尼較為接近。在采用快速貝葉斯FFT方法識別氣動阻尼時，需要注意兩點[3]：一是選擇不同共振頻率帶會得到不同的阻尼比識別結果，因此需要多選擇幾次才能得到理想的結果；二是采樣數據盡量要多一些，識別精度會有所增加。

由于快速貝葉斯FFT是直接對原始數據進行傅里葉變換，并考慮了統計估計誤差，因此其在識別阻尼比時減少了人為因素的影響。另外，大多數模態識別方法都無法確定阻尼比的識別精度，而快速貝葉斯FFT方法可以根據模態參數的后驗變異系數確認各模態參數的識別精度。

5 結論

本文對傳統的氣動阻尼識別方法進行了闡述分析，并進一步詳細地介紹了HHT、小波分析和快速貝葉斯FFT，得出以下結論：

(1)譜分析方法要求信號是平穩的和線性的，而時間序列法中的隨機減量法則可以適用于非線性系統，從而可以更廣泛的應用于結構的氣動阻尼識別。(2)小波變換和HHT變換都是處理非平穩非線性信號的時頻分析方法。HHT變換具有良好的自適應性，但是沒有完善數學基礎，同時存在模態混淆和端點效應問題；而小波變換具有完善的數學基礎，可以對信號進行詳細的時頻分析和多分辨率分析，但存在邊界效應問題。(3)快速貝葉斯FFT方法考慮了統計估計誤差，減少了人為因素對阻尼比識別的影響，并且快速貝葉斯FFT方法可以確認對各模態參數的識別精度。