一種新的固有頻率跟蹤及虛假頻率剔除方法*

周昊天， 趙 婕， 曹為政, 趙 銳， 于開平， 白云鶴

(1.哈爾濱工業(yè)大學(xué)航天學(xué)院 哈爾濱，150001) (2.黑龍江科技大學(xué)計(jì)算機(jī)與信息工程學(xué)院 哈爾濱，150022) (3.中國電子科技集團(tuán)公司第54研究所 石家莊，050081)

引 言

結(jié)構(gòu)模態(tài)參數(shù)辨識是結(jié)構(gòu)動力學(xué)的重要研究內(nèi)容，當(dāng)前已經(jīng)有眾多成熟方法用于解決定常系統(tǒng)的非時(shí)變模態(tài)參數(shù)辨識問題，例如隨機(jī)減量法[1]及隨機(jī)子空間方法[2]等。由于線性時(shí)變系統(tǒng)比定常系統(tǒng)復(fù)雜，同時(shí)，適用于定常非時(shí)變系統(tǒng)的模態(tài)參數(shù)辨識方法并不能直接應(yīng)用于線性時(shí)變系統(tǒng)，而時(shí)變系統(tǒng)在實(shí)際工程應(yīng)用中又是普遍存在的，因此對線性時(shí)變系統(tǒng)的模態(tài)參數(shù)辨識一直是研究的熱點(diǎn)。

在時(shí)變系統(tǒng)中，傳統(tǒng)的模態(tài)參數(shù)的概念不復(fù)存在。Liu[3]提出了“偽模態(tài)參數(shù)”的概念，并將其推廣用于線性時(shí)變系統(tǒng)[4]。目前，已有的用于時(shí)變模態(tài)參數(shù)辨識的方法包括基于時(shí)頻分析的辨識方法[5]、基于盲源分離的辨識方法[6]和基于時(shí)間序列模型的方法[7]。對于工程實(shí)際而言，應(yīng)用最為廣泛的是基于子空間的辨識方法[8]。Liu等[9]提出一種基于子空間的時(shí)變模態(tài)參數(shù)辨識方法，并對一個(gè)移動質(zhì)量梁的偽固有頻率進(jìn)行了識別。于開平等[10]提出一種利用系統(tǒng)輸入輸出整體數(shù)據(jù)進(jìn)行時(shí)變模態(tài)參數(shù)辨識的方法，對一個(gè)移動簡支梁的偽固有頻率進(jìn)行了辨識，并提出改進(jìn)方法[11]。龐世偉等[12]引入投影估計(jì)子空間算法(projection approximation subspace tracking，簡稱PAST)，提出了采用固定長度平移窗的遞推子空間辨識方法，利用該方法對移動質(zhì)量簡支梁進(jìn)行了辨識。楊凱等[13]引入自然冪迭代算法(natural power iteration，簡稱NPI)和逼近冪迭代算法(approximated power iteration，簡稱API)，對移動分布質(zhì)量-懸臂梁和兩自由度彈簧-質(zhì)量系統(tǒng)進(jìn)行了時(shí)變模態(tài)參數(shù)辨識。倪智宇等[14]提出了一種用于時(shí)變模態(tài)參數(shù)辨識的截?cái)啻氨平鼉绲涌臻g方法(truncated window approximated power iteration，簡稱TW-API)，并將該方法用于二連桿機(jī)械臂和在軌航天器[15]的偽固有頻率辨識中。

由于子空間跟蹤算法本身的限制和測量數(shù)據(jù)受噪聲污染的影響，識別結(jié)果中通常混雜了虛假的模態(tài)信息[9]，因此提高最終識別結(jié)果的精度是一項(xiàng)重要的研究內(nèi)容。除了采用精度更高的跟蹤算法以外，另一種研究思路是對識別結(jié)果進(jìn)行處理，采用特定的算法識別出真實(shí)模態(tài)、剔除虛假模態(tài)。Liu等[9]提出了一種針對慢變的偽固有頻率篩選方法，該方法計(jì)算量小，但需要給出若干估計(jì)的頻率值。Zhou等[16]引入fuzzy C-means方法，提出基于模糊聚類的模態(tài)參數(shù)驗(yàn)證方法，即通過歐氏距離檢測不同階次之間的差異，從而實(shí)現(xiàn)對偽固有頻率的篩選。

為了提高最終結(jié)果的精度，筆者引入新信息準(zhǔn)則子空間跟蹤算法[17](novel information criterion，簡稱NIC)，針對時(shí)域類模態(tài)參數(shù)辨識方法常見的虛假頻率問題，提出一種基于滑動數(shù)據(jù)窗的固有頻率篩選方法。兩種不同頻率變化規(guī)律的數(shù)值仿真算例驗(yàn)證了辨識算法的有效性，辨識結(jié)果精度得到了提升。通過對高溫環(huán)境下的熱防護(hù)結(jié)構(gòu)的時(shí)變偽固有頻率進(jìn)行辨識，驗(yàn)證了篩選方法的有效性。

1 NIC子空間跟蹤算法與時(shí)變模態(tài)參數(shù)辨識

1.1 輸入量的迭代更新

由于基于NIC的子空間跟蹤方法是采用遞推最小二乘算法求解其中的子空間優(yōu)化問題，因此振動信號的前處理是必需步驟。考慮線性時(shí)變系統(tǒng)離散時(shí)間狀態(tài)空間模型

(1)

其中：x(k)∈Rn×1為第k個(gè)采樣時(shí)刻的系統(tǒng)狀態(tài)向量；n為系統(tǒng)階次；A(k)∈Rn×n為k～k+1時(shí)刻的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣；B(k)∈Rn×m為第k時(shí)刻的輸入矩陣；u(k)∈Rm×1為第k時(shí)刻的m維輸入向量；y(k)∈Rr×1為第k時(shí)刻的r維響應(yīng)向量；C(k)∈Rr×n為第k時(shí)刻的輸出矩陣。

利用k時(shí)刻附近的輸入輸出數(shù)據(jù)可組成相應(yīng)的廣義Hankel矩陣

(2)

其中：M為廣義Hankel矩陣的塊行數(shù)。

對應(yīng)地，下一時(shí)刻的廣義Hankel矩陣為

(3)

(4)

根據(jù)文獻(xiàn)[8]提出的數(shù)據(jù)迭代更新規(guī)則，令

得到子空間跟蹤輸入量z(k)的迭代更新形式

(7)

其中：上標(biāo)符號“?”代表矩陣偽逆。

式(5)和式(7)中的[U(k-1)UT(k-1)]-1以及[Y(k-1)U?(k)]更新形式分別為

通過以上公式可得到滿足子空間跟蹤算法要求的輸入數(shù)據(jù)z(k)，且有

其中：Γ(k)為k時(shí)刻的能觀矩陣。

通過對Y(k)U⊥(k)進(jìn)行奇異值分解，得到

Y(k)U⊥(k)=R(k)Σ(k)V(k)T

(9)

則R(k)的前n列組成的矩陣即為能觀矩陣Γ(k)，同時(shí)也是矩陣Y(k)U⊥(k)的信號子空間W(k)。通過尋找合適的跟蹤算法可以避免計(jì)算量較大的奇異值計(jì)算步驟，達(dá)到節(jié)省計(jì)算量的目的。

1.2 NIC子空間跟蹤

NIC算法將子空間跟蹤視為無約束優(yōu)化問題，NIC算法選取的目標(biāo)方程為

(10)

其中：W為信號子空間矩陣；tr為矩陣的跡；Rpp=E[ppT]為輸入量的協(xié)方差矩陣。

J(W)對W的梯度可以寫為

(11)

根據(jù)梯度上升規(guī)則，W(k)的迭代更新規(guī)則可寫為

(12)

其中：η為學(xué)習(xí)步長。

(13)

(14)

(15)

(16)

(17)

初始化步驟：

1) 令P(0)=δIr，其中：δ為一個(gè)正小量；Ir為r×r單位矩陣；

3)W(0)通過式(9)進(jìn)行計(jì)算。

迭代更新步驟：

1.3 時(shí)變模態(tài)參數(shù)提取

(18)

對系統(tǒng)矩陣的估計(jì)進(jìn)行特征值分解，得到

(19)

其中：ψ(k)為特征向量矩陣，Λ(k)=diag(λ1(k),…,λn(k))為特征值對角矩陣。

由此可以求得k時(shí)刻第i階偽固有頻率ωi(k)和偽模態(tài)阻尼比ξi(k)分別為

其中：上標(biāo)“R”表示取特征值的實(shí)部；Δt為數(shù)據(jù)采樣間隔。

2 仿真算例

為驗(yàn)證所提出算法的有效性以及抗噪能力，對如圖1所示的兩自由度彈簧-質(zhì)量系統(tǒng)進(jìn)行辨識。

圖1 2自由度彈簧-質(zhì)量系統(tǒng)示意圖Fig.1 Schematic diagram of the 2-DOF spring-mass model

該彈簧-質(zhì)量系統(tǒng)的運(yùn)動控制方程為

(21)

質(zhì)量矩陣M為

(22)

剛度矩陣K為

(23)

阻尼矩陣C為

(24)

為了研究所提出的方法對不同頻率變化模式的適用情況，對剛度為指數(shù)變化和剛度為正弦變化兩種情況進(jìn)行討論。

2.1 剛度為指數(shù)變化的情況

為模擬頻率為指數(shù)變化的模式，式(23)中的各剛度取值為K1=1 000exp(-0.4t)，K2=2 000。仿真計(jì)算采用的時(shí)間步長Δt=0.001 s。分別采用文獻(xiàn)[12]和筆者提出的辨識方法進(jìn)行辨識，為模擬真實(shí)的測量情況，結(jié)構(gòu)響應(yīng)中加入一定量的噪聲。兩種方法均采用相同的計(jì)算參數(shù)，結(jié)構(gòu)響應(yīng)信號信噪比(signal to noise ratio，簡稱SNR)為10時(shí)，辨識結(jié)果如圖2，3所示。可以看出，PAST和NIC方法都能很好地跟蹤頻率變化。在初始時(shí)刻，PAST的誤差較大，而NIC算法的結(jié)果更接近理論參考值。在跟蹤過程中，PAST算法結(jié)果呈現(xiàn)出較大的離散性，而NIC算法的結(jié)果相對比較平滑。為了量化評估在不同信噪比條件下兩種算法的計(jì)算結(jié)果，引入平均絕對誤差率(mean absolute percentage error，簡稱MAPE)

(25)

圖2 PAST辨識結(jié)果與理論參考值對比(第1組)Fig.2 Comparison of the results by PAST algorithm and reference results (the first group)

圖3 NIC辨識結(jié)果與理論參考值對比(第1組)Fig.3 Comparison of the results by NIC algorithm and reference results (the first group)

通過兩種算法5次計(jì)算得到系統(tǒng)固有頻率的MAPE誤差取均值如表1所示。可以看出，當(dāng)結(jié)構(gòu)響應(yīng)信號的信噪比較高時(shí)，兩種算法的MAPE誤差相差不大。隨著噪聲量的增加，PAST算法的MAPE誤差呈現(xiàn)出大幅增加的趨勢，NIC算法辨識誤差雖然也有增加，但是仍然優(yōu)于PAST算法的辨識結(jié)果。

表1 辨識結(jié)果MAPE比較(第1組)

Tab.1MAPE comparison of the identified results(the first group)%

SNRPASTNIC第1階第2階第1階第2階5011.404.2910.243.722011.745.8210.403.921011.896.3410.805.70512.2412.9110.838.67

2.2 剛度為正弦變化的情況

式(23)中，剛度取值為K1=2 000+1 000×sin(πt/5)，K2=2 000，仿真計(jì)算的時(shí)間步長Δt=0.01 s。信噪比為5時(shí)的辨識結(jié)果如圖4，5所示。可以看出，兩種算法對正弦形式的頻率變化都可以很好地跟蹤辨識，兩種算法結(jié)果差異不大，在個(gè)別位置如2 s時(shí)刻左右，NIC算法得到的結(jié)果較PAST算法得到的結(jié)果更接近理論參考值。

圖4 PAST辨識結(jié)果與參考值對比(第2組)Fig.4 Comparison of the results by PAST algorithm and reference results (the second group)

圖5 NIC辨識結(jié)果與理論參考值對比(第2組)Fig.5 Comparison of the results by NIC algorithm and reference results (the second group)

不同信噪比情況下得到的MAPE誤差如表2所示。可以看出，兩種方法的識別結(jié)果誤差都在10%以內(nèi)，在添加了20%的噪聲以后都可以保證很好的識別精度，NIC算法在不同SNR情況下的識別結(jié)果誤差均不同程度地優(yōu)于傳統(tǒng)的PAST算法。

表2 辨識結(jié)果MAPE比較(第2組)

Tab.2MAPE comparison of the identified results(the second group)%

SNRPASTNIC第1階第2階第1階第2階504.821.524.131.35204.931.854.271.70105.211.914.591.8455.512.594.992.67

3 實(shí)驗(yàn)驗(yàn)證與應(yīng)用

飛行器熱防護(hù)結(jié)構(gòu)在溫度變化環(huán)境下的振動特性一直受到研究人員的關(guān)注，這里對一塊陶瓷基復(fù)合材料熱防護(hù)板的偽固有頻率進(jìn)行辨識。熱防護(hù)結(jié)構(gòu)懸吊在剛性結(jié)構(gòu)上，加熱面通過遠(yuǎn)紅外熱輻射加熱至800℃，冷端此時(shí)溫度為80℃，停止加熱使熱防護(hù)結(jié)構(gòu)自然降溫。在冷端通過激振器施加可視為高斯白噪聲的隨機(jī)載荷，并由4個(gè)加速度傳感器進(jìn)行數(shù)據(jù)采集，采樣頻率為1 400 Hz。實(shí)驗(yàn)裝置如圖6所示。

圖6 實(shí)驗(yàn)裝置圖Fig.6 Image of the experiment setup

利用NIC方法對偽固有頻率進(jìn)行辨識，計(jì)算參數(shù)Hankel矩陣的塊行數(shù)M=100，遺忘因子β=0.94，辨識結(jié)果如圖7所示。可見，大部分辨識結(jié)果都聚集在220，300和440 Hz左右，其他點(diǎn)可視為由于測量噪聲等原因?qū)е碌奶摷倌B(tài)。

圖7 NIC方法辨識結(jié)果Fig.7 Identified results using NIC algorithm

采用比較篩除的方法剔除虛假模態(tài)點(diǎn)，定義一個(gè)加權(quán)平均偽固有頻率為

(26)

其中：θ為加權(quán)因子；fs(k)為經(jīng)過篩選的k時(shí)刻的偽固有頻率。

圖8 篩選后的辨識結(jié)果Fig.8 Identified results with selection and sifting procedure

圖8為篩選后的辨識結(jié)果。可以看出，自由狀態(tài)下熱防護(hù)結(jié)構(gòu)從800℃自然降溫的150 s內(nèi)，前3階固有頻率分別由221，295，434 Hz升至230，304，444 Hz。經(jīng)過篩選算法處理的結(jié)果不再包含虛假模態(tài)，結(jié)構(gòu)的時(shí)變趨勢更加容易分辨。

4 結(jié)束語

提出了一種基于NIC子空間跟蹤算法的偽固有頻率辨識方法，通過一個(gè)兩自由度彈簧-質(zhì)量系統(tǒng)兩種不同頻率變化形式的仿真數(shù)值算例，驗(yàn)證了算法的有效性。數(shù)值算例表明，該方法在低SNR情況下仍有較好的辨識結(jié)果。與傳統(tǒng)的PAST方法相比，該方法辨識結(jié)果與理論參考值更接近，具備更高的識別精度。針對時(shí)域辨識算法常見的虛假固有頻率問題，從辨識結(jié)果后處理的角度出發(fā)，提出了一種虛假固有頻率剔除方法。該方法具備計(jì)算量小、計(jì)算所需參數(shù)少的特點(diǎn)。最后，識別了飛行器熱防護(hù)結(jié)構(gòu)降溫過程的偽固有頻率變化過程。實(shí)驗(yàn)結(jié)果表明，該方法處理效果明顯，所提出的辨識方法和后處理方法具備較高的應(yīng)用價(jià)值。