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

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

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

引 言

結構模態參數辨識是結構動力學的重要研究內容，當前已經有眾多成熟方法用于解決定常系統的非時變模態參數辨識問題，例如隨機減量法[1]及隨機子空間方法[2]等。由于線性時變系統比定常系統復雜，同時，適用于定常非時變系統的模態參數辨識方法并不能直接應用于線性時變系統，而時變系統在實際工程應用中又是普遍存在的，因此對線性時變系統的模態參數辨識一直是研究的熱點。

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

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

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

1 NIC子空間跟蹤算法與時變模態參數辨識

1.1 輸入量的迭代更新

由于基于NIC的子空間跟蹤方法是采用遞推最小二乘算法求解其中的子空間優化問題，因此振動信號的前處理是必需步驟?？紤]線性時變系統離散時間狀態空間模型

(1)

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

利用k時刻附近的輸入輸出數據可組成相應的廣義Hankel矩陣

(2)

其中：M為廣義Hankel矩陣的塊行數。

對應地，下一時刻的廣義Hankel矩陣為

(3)

(4)

根據文獻[8]提出的數據迭代更新規則，令

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

(7)

其中：上標符號“?”代表矩陣偽逆。

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

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

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

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

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

(9)

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

1.2 NIC子空間跟蹤

NIC算法將子空間跟蹤視為無約束優化問題，NIC算法選取的目標方程為

(10)

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

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

(11)

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

(12)

其中：η為學習步長。

(13)

(14)

(15)

(16)

(17)

初始化步驟：

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

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

迭代更新步驟：

1.3 時變模態參數提取

(18)

對系統矩陣的估計進行特征值分解，得到

(19)

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

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

其中：上標“R”表示取特征值的實部；Δt為數據采樣間隔。

2 仿真算例

為驗證所提出算法的有效性以及抗噪能力，對如圖1所示的兩自由度彈簧-質量系統進行辨識。

圖1 2自由度彈簧-質量系統示意圖Fig.1 Schematic diagram of the 2-DOF spring-mass model

該彈簧-質量系統的運動控制方程為

(21)

質量矩陣M為

(22)

剛度矩陣K為

(23)

阻尼矩陣C為

(24)

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

2.1 剛度為指數變化的情況

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

(25)

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

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

通過兩種算法5次計算得到系統固有頻率的MAPE誤差取均值如表1所示?？梢钥闯?，當結構響應信號的信噪比較高時，兩種算法的MAPE誤差相差不大。隨著噪聲量的增加，PAST算法的MAPE誤差呈現出大幅增加的趨勢，NIC算法辨識誤差雖然也有增加，但是仍然優于PAST算法的辨識結果。

表1 辨識結果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，仿真計算的時間步長Δt=0.01 s。信噪比為5時的辨識結果如圖4，5所示。可以看出，兩種算法對正弦形式的頻率變化都可以很好地跟蹤辨識，兩種算法結果差異不大，在個別位置如2 s時刻左右，NIC算法得到的結果較PAST算法得到的結果更接近理論參考值。

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

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

不同信噪比情況下得到的MAPE誤差如表2所示?？梢钥闯觯瑑煞N方法的識別結果誤差都在10%以內，在添加了20%的噪聲以后都可以保證很好的識別精度，NIC算法在不同SNR情況下的識別結果誤差均不同程度地優于傳統的PAST算法。

表2 辨識結果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 實驗驗證與應用

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

圖6 實驗裝置圖Fig.6 Image of the experiment setup

利用NIC方法對偽固有頻率進行辨識，計算參數Hankel矩陣的塊行數M=100，遺忘因子β=0.94，辨識結果如圖7所示?？梢?，大部分辨識結果都聚集在220，300和440 Hz左右，其他點可視為由于測量噪聲等原因導致的虛假模態。

圖7 NIC方法辨識結果Fig.7 Identified results using NIC algorithm

采用比較篩除的方法剔除虛假模態點，定義一個加權平均偽固有頻率為

(26)

其中：θ為加權因子；fs(k)為經過篩選的k時刻的偽固有頻率。

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

圖8為篩選后的辨識結果?？梢钥闯?，自由狀態下熱防護結構從800℃自然降溫的150 s內，前3階固有頻率分別由221，295，434 Hz升至230，304，444 Hz。經過篩選算法處理的結果不再包含虛假模態，結構的時變趨勢更加容易分辨。

4 結束語

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