基于改進MCSSD的時變系統模態參數識別

張志超，史治宇，張杰

(南京航空航天大學 機械結構力學及控制國家重點實驗室，江蘇 南京 210016)

0 引言

伴隨著科技的發展和社會的進步，在航空航天、土木建筑等領域，結構動力學參數的時變特性越來越引起工程師和研究人員的注意，與此同時，有關時變參數影響的動力學工程需求越來越迫切[1]。例如，大型運載火箭或導彈在發射過程中隨燃料消耗質量減小所產生的結構時變特性；列車高速過橋所引起的車橋系統時變特性；飛行器飛行過程中狀態監測所表現出的時變特性等。要解決這類時變結構的動力學問題，時變結構的參數識別方法研究便顯得尤為重要。

隨著研究的深入，國內外學者提出了不同的時變系統參數識別方法。20世紀90年代Feldman[2-3]應用Hilbert變化對單自由度非線性振動自由響應信號和受迫響應信號進行了研究分析，識別了結構的瞬時頻率、阻尼等時變參數。Staszewski和Cooper[4]將連續小波引入到動力學結構分析當中，之后Ruzzene等[5]和Staszewski[6]根據脈沖響應函數的Morlet小波變換與結構模態參數之間的關系識別了結構的動力學信息。Argoul[7]還把Cauchy小波也運用到參數識別當中。1997年Liu[8]將線性時不變系統的子空間方法推廣到時變系統，提出了基于自由響應信號的線性時變系統偽模態參數識別方法。1999年其[9]再次提出了利用任意激勵下輸入輸出信號的時變系統子空間參數識別方法。

上述時變系統參數識別方法在不同方面均存在一定的缺陷。基于小波變換時頻分析，其使用的平行劃分的等面積時頻窗，不能自適應地分解出密集模態成分[10]；而經驗模態分解雖然能自適應地將輸出響應分解為各個對應于各階模態的固有模態函數(Intrinsic mode function，IMF)，但仍然存在端點效應和模態混疊等問題[11]。

Mallat等[12]人于1993年提出了稀疏信號分解的概念，根據信號本身的特點，在過完備原子庫上，自適應的選擇基函數對信號進行投影分解。CANDES等[13]人于2007年提出了線調頻小波路徑追蹤算法，由于此算法僅能對單一分量信號進行分解，目前也只應用于地震引力波分析中。PENG等[14]人在此基礎上提出多尺度線調頻基稀疏分解方法，可以對多分量信號進行多尺度分解，此方法多用于故障診斷。

本文在多尺度線性調頻原子稀疏信號分解及多項式相位信號參數估計理論的基礎上，提出一種時變系統的模態參數識別方法。對于n自由度時變系統的自由響應是n個獨立單模態自由響應的線性疊加，因此通過多尺度線調頻原子稀疏信號分解可得到n個獨立分量，最小二乘瞬時頻率擬合可以保證瞬時頻率的擬合精度。至此，n自由度系統自由響應在時域上達到解耦，將多自由度時變系統參數識別轉換為單自由度時變系統參數識別，根據分解得到的幅值估計和瞬時頻率估計，可以計算得到系統各階模態頻率和模態阻尼比，進而達到多自由度時變系統模態參數識別。

1 改進多尺度線性調頻原子稀疏信號分解

1.1 信號分解

根據信號分析理論，任意信號f(t) 可以展開一組基的線性組合，如：

(1)

式中：hn為基函數，an為基函數所對應的系數。

如果選取的基函數正交，則其展開系數可通過內積來計算，如：

(2)

為使信號分解過程中不出現能量泄漏，使分解信號具有重構性，使分解系數具有唯一性，通常選取的基函數為標準正交基。本文選取多尺度線調頻原子作為基函數[15]，其表達式如：

Dha,b,I=ha,b,It=Ka,b,Iexp-iat+bt2lIt

(3)

根據稀疏信號分解理論，在第I個支撐區上對原信號進行投影，其投影系數按式(4)計算。

(4)

在第I個支撐區上最大投影系數代表的分析信號為：

cIt=2βIexp-iat+bt2-angle(2βI)lIt

(5)

所以在第I個支撐區上有：

fIt=cIt+rIt

(6)

式中：rIt為殘余信號。

根據線性調頻小波路徑追蹤算法，要找到一種動態時間支撐區的鏈接方法，將每個動態時間支撐區上投影系數最大的線性調頻原子按照此種連接算法進行連接，此種方法要保證信號殘余能量最小，即分量信號能量最大。

(7)

基于以上理論，對多尺度線調頻原子進行最優鏈接。此算法可以保證分量信號與原信號中同頻率成分信號分量最接近，即保證了信號分離的精度。

1.2 瞬時頻率估計

在第I個動態時間支撐區上，瞬時頻率如：

(8)

式中：angle2βI在每一個動態時間支撐區上為常數，因此對時間求導為0。

將不同的動態時間支撐區的瞬時頻率按照上述最優鏈接算法進行鏈接，則可得到完整頻率序列。

(9)

但由于噪聲等其他因素的影響，使得在分解時所選取的線調頻原子存在一定誤差，因此直接按照最優鏈接算法鏈接得到的頻率序列曲線不光滑，同時存在較大的誤差，對后續計算帶來很大的影響。為了解決上述存在的問題，本文將各個動態時間支撐區首尾節點頻率提取出來，根據多項式相位信號參數估計理論[16]，利用最小二乘法對提取節點頻率進行多項式擬合，當誤差達到閾值時，擬合多項式即為完整頻率序列，進而頻率估計轉化為多項式參數求解。

由多項式相位信號參數估計理論可知，由多尺度線調頻原子分解得到的分量信號瞬時頻率可由多項式進行逼近。

(10)

最小二乘曲線擬合核心思想就是尋找形如式(10)的n階多項式，使：

(11)

式(11)滿足了所擬合完整頻率序列光滑，但由于分量信號不具有先驗信息，所以要達到較高的精度，需要通過增加多項式階數來實現，仿真數據發現當階數增加到一定程度時，前后擬合頻率序列相差較小，因此定義連續兩次擬合頻率序列相對差值控制擬合精度，即：

(12)

式中：δ為閾值，當前后擬合頻率序列相對誤差小于閾值時，即達到擬合精度。

則具體瞬時頻率估計步驟如下：

1) 根據多尺度線調頻原子稀疏信號分解方法得到信號分量，設定閾值δ，取初始多項式階數為1；

2) 根據最優線調頻原子鏈接算法，提取每個動態時間支撐區首尾節點瞬時頻率；

3) 根據最小二乘多項式相位信號參數估計方法對頻率序列進行項擬合，計算前后兩次相對差值ξδ；

4) 如果ξδ>δ，增加擬合多項式階數，返回步驟3)；如果ξδ<δ，停止擬合，得到頻率序列的最小二乘估計多項式形式。

2 基于改進多尺度線調頻基信號分解的時變模態參數識別

2.1 時變系統振動響應分析

對于n自由度粘性阻尼系統自由振動方程可表示為如下形式：

(13)

式中：m,k,c是n×n階質量，剛度和阻尼矩陣；Xt=x1,x2,…,xnT為響應矩陣。

(14)

式(14)特征方程為：

sA+By=0

(15)

Dy=sy

(16)

采用“凍結法”對式(16)進行特征值求解，第i對復特征值為：

(17)

將式(17)帶入式(15)，可求得特征值對應的特征向量為：

(18)

定義坐標變換：

(19)

至此系統方程可以在模態坐標下進行解耦，系統自由響應可以表示為：

(20)

則第k點的響應可以表示為：

(21)

式中：Ajt=2ηijtgjtexp-αjtt。

2.2 模態參數識別

由式(21)可以得出，n自由度時變系統的第k點自由振動響應可以表示為n個獨立模態自由振動響應的線性疊加。改進多尺度線調頻原子稀疏信號分解可以將響應信號分解為n個獨立模態自由振動響應信號，至此n自由度時變系統實現了時域解耦。

對于具有粘性阻尼的單自由度時變系統，其自由振動微分方程可表達為：

(22)

其位移響應如：

ut=Atcosθt

(23)

構造解析信號zt如：

(24)

將式(24)帶入式(22)得：

(25)

(26)

解方程式(26)得到系統模態頻率和模態阻尼比為：

(27)

式中：At和θt為改進多尺度線調頻原子稀疏信號分解中得到的幅值包絡和頻率序列。

本文所采用的基于改進多尺度線調頻原子稀疏信號分解的時變系統模態參數識別方法具體過程如下：

1) 測量系統自由響應數據；

2) 對位移響應信號進行n(n為系統自由度數)次改進多尺度線調頻原子稀疏信號分解，得到n個信號分量；

3) 對n個分量的各節點瞬時頻率進行最小二乘逼近，得到頻率序列；

4) 將得到的幅值包絡At和頻率序列θt代入式(27)，求得系統的模態頻率和模態阻尼比。

3 仿真信號分析

選取一個3自由度的彈簧阻尼系統作為算例，其結構如圖1所示。

圖1 3自由度系統示意圖

初始時系統參數如下：

不同于其他參數單一變化，本文考慮結構多參數變化，主要包括線性變化和周期變化2種情況。

1) 參數隨時間線性變化

假設彈簧k1隨時間線性變化：k1=50×tN/m;彈簧k2隨時間線性變化：k2=40×tN/m；阻尼c1隨時間線性變化：c1=0.05×tN·s/m；阻尼c2隨時間線性變化：c2=0.04×tN·s/m。利用Newmarkβ算法對結構計算10s內位移響應信號。選取位移響應較為復雜的m3來作為分析信號，m3位移響應信號如圖2所示。

圖2 m3位移響應信號

對m3位移響應信號進行3次改進多尺度線調頻原子稀疏信號分解，取搜索頻率0～20Hz，搜尋最大頻偏2，得到3個分量信號，對其各節點頻率進行最小二乘逼近，最終分量信號如圖3所示。

圖3 3個分量信號

將3個分量信號的幅值包絡及擬合逼近得到的頻率序列帶入式(27)得到模態頻率，最終得到的模態頻率識別結果如圖4所示。

圖4 模態頻率識別結果

從圖4可以看出，改進方法的各階模態頻率的識別精度都很高，第一階和第二階模態頻率比較接近，但也都較好的被識別出來了，識別值相對理論值存在一定的偏差，但曲線趨勢和理論值趨勢基本一致，與理論值偏差存在的原因在于幅值包絡提取不夠精準。

MAPE值越小，代表識別結果總體誤差越小。對模態頻率識別值進行平均絕對百分誤差計算，結果如表1所示。

表1 識別結果平均絕對百分誤差

由表1可以看出，改進線調頻原子稀疏信號分解模態參數識別方法的各階識別精度較高，但第一階相對于其他兩階存在較大的誤差，原因在于第一階能量較小，同時是最后一次分離出來，受前兩次分離殘余信號響應，使得第一階模態頻率識別存在相對于其他兩階較大的誤差，但仍在誤差允許范圍內，因此總體上有很高的精度。

2) 參數隨時間周期變化

假設彈簧k1隨時間周期變化：k1=120cos0.5πt；彈簧k2隨時間周期變化：k2=60cos0.5πt；阻尼c1隨時間周期變化：c1=0.01cos0.4πt；阻尼c2隨時間周期變化：c2=0.02cos0.4πt。利用Newmarkβ算法對結構計算10s內位移響應信號。選取位移響應較為復雜的m3來作為分析信號，m3位移響應信號如圖5所示。

圖5 m3位移響應信號

對m3位移響應信號進行3次改進多尺度線調頻原子稀疏信號分解，取搜索頻率0～20Hz，搜尋最大頻偏1，得到3個分量信號，對其各節點頻率進行最小二乘逼近，最終分量信號如圖6所示。

圖6 3個分量信號

將3個分量信號的幅值包絡及擬合逼近得到的頻率序列帶入式(27)得到模態頻率，最終得到的模態頻率識別結果如圖7所示。

圖7 改進方法模態頻率識別結果

為了對比改進方法的精度提高，同時對m3位移響應信號用原方法進行分解，得到3個信號分量，將信號分量幅值包絡及頻率序列帶入式(27)得到模態頻率，最終得到的模態頻率如圖8所示。

圖8 原方法模態頻率識別結果

分別計算改進方法與原方法所識別的模態頻率MAPE值，如表2所示。

由表2可知，改進線調頻原子稀疏信號分解模態參數識別方法較原方法提高了模態頻率的識別精度，但由于信號分解過程的殘余信號對下次分解的影響，使得原本能量較小的第一階頻率的識別受到了很大的影響，但總體的精度有了很大的提升，各階的識別誤差也在允許范圍內。

表2 誤差對比

4 結語

1)n自由度時變系統的自由振動響應信號是由n個獨立單模態自由振動響應信號的線性疊加，而每一個獨立模態自由振動響應信號又可以表達為多個線調頻信號的線性組合。因而，多尺度線調頻原子稀疏信號分解可以將自由振動響應信號分解為n個獨立分量，這n個獨立分量正是n個單模態自由振動響應。因此n自由度時變系統自由響應信號實現時域上解耦。

2) 由于多尺度線調頻基稀疏信號分解得到的分量是調頻信號，對于瞬時頻率的提取存在一定的誤差，根據多項式相位信號參數估計理論，采用最小二乘法對各節點頻率進行擬合逼近，通過提高多項式階數，保證了各分量瞬時頻率提取的準確度，進而保證了參數識別的準確度。

3) 仿真算例中誤差對比可以看出，本文所提出的改進方法相對于原方法，在識別精度上有了很大提升，但由于此方法還受到幅值提取以及殘余信號等影響，已然存在一定的誤差，但總體精度滿足工程需要。

4) 由于此種方法不需要信號的先驗信息，只需測量結構的位移響應信號，因此具有較好的工程實用價值。

參考文獻:

[1] 于開平, 龐世偉, 趙婕. 時變線性/非線性結構參數識別及系統辨識方法研究進展[J]. 科學通報, 2009, 54(20): 3147-3156.

[2] Feldman M. Non-linear system vibration analysis using Hilbert transform - I: free vibration analysis method FREEVIB[J]. Mechanical Systems and Signal Processing, 1994, 8(2): 119-127.

[3] Feldman M. Non-linear system vibration analysis using Hilbert transform - II: forced vibration analysis method FORCEVIB[J]. Mechanical Systems and Signal Processing, 1994, 8(3): 309-318.

[4] Staszewski W J, Cooper J E. Flutter data analysis using the wavelet transform, Proceedings of international Congress on MV2: New Advances in Modal Synthesis of Large Structures[J]. Non-linear, Damped and Non-Deterministic Cases, Lyon, 1995(11):549-561.

[5] Ruzzene M, Fasana A, Garibaldi L, Piombo B. Natural frequencies and dampings identification using wavelet transform: application to real data[J]. Mechanical Systems and Signal Processing, 1997,11(2):207-218.

[6] Staszewski W J. Identification of non-linear systems using multi-scale ridges and skeletons of the wavelet transform[J]. Journal of Sound and Vibration, 1998,214(4): 639-658.

[7] Argoul P, Le T P. Instantaneous indicators of structural behavior based on the continuous Cauchy wavelet analysis[J]. Mechanical Systems and Signal Processing, 2003, 17: 243-250.

[8] Liu K. Identification of linear time-varying systems[J]. Journal of Sound Vibration, 1997, 204: 487-500.

[9] Liu K. Extension of modal analysis to linear time-varying systems[J]. Journal of Sound Vibration, 1999, 226(1): 149-167.

[10] 許鑫， 史治宇. 狀態空間下基于小波變換的時變系統參數識別 [J]. 振動工程學報， 2010， 23 (4) : 416-419.

[11] 徐晴晴，史治宇. 基于改進EMD分解的時變結構密集模態的瞬時參數識別[J]. 機械科學與技術，2015, 34(8):1161-1165.

[12] Mallat S, Zhang Z. Matching pursuit with time-frequency dictionaries [J]. Signa l Processing, 1993, 41 (12): 3397-3415.

[13] Candès E J, Charlton P R, Helgason H. Detecting highly oscillatory signals by chirplet path pursuit [J]. Applied and Computational Harmonic Analysis, 2008,24(1): 14-40.

[14] 彭富強，于德介，羅思潔，等. 基于多尺度線調頻基稀疏信號分解的軸承故障診斷[J]. 機械工程學報，2010，46(7) : 88-95.

[15] 彭富強, 于德介, 劉堅. 一種基于多尺度線調頻基的稀疏信號分解方法[J]. 振動工程學報, 2010, 23(3):333-338.

[16] 李越雷，張天騏，黃銚，等. 基于曲線擬合的多項式相位信號參數估計[J]. 計算機應用，2009，29:349-352.