基于低秩Hankel矩陣逼近的模態參數識別方法

包興先，李昌良，劉志慧

（中國石油大學（華東）石油工程學院，山東 青島 266580）

模態參數識別技術廣泛應用于工程領域，已成為解決各類工程問題的重要手段。傳統模態參數識別方法建立在已知系統輸入、輸出基礎上，據輸入輸出信息利用實驗模態分析獲得系統模態參數。但工程實際中激勵信號較難獲得，通常只獲得含噪聲的響應信號，因此基于輸出響應的模態參數識別技術頗受重視。受測量噪聲影響，識別結果通常含虛假模態，如何辨識并剔除虛假模態成為參數識別過程關鍵。目前模態參數識別通常為計算用模型階次高于真實模型階次，引入“噪聲模態”［1］。為區分真實、虛假模態大多用穩定圖方法［2］。穩定圖可表示模態參數與模型階次關系，理論上隨模型階次的增加，真實模態參數趨于穩定而虛假模態不會，但實際上識別結果中真實模態參數也存在不穩定性，導致穩定圖法失效。尤其隨模型階次的升高，虛假模態易趨于穩定，用穩定圖較難正確識別結構真實模態參數［3－4］。尤其信噪比較低時區分大量虛假、真實模態更困難。為準確確定模型階次，基于奇異值分解方法［4－7］一般可通過觀察確定奇異值曲線突變點對應的奇異值個數即模型階次。王樹青等［8］提出基于奇異值相對變化率方法，引入量化模型階次指標，以該指標最大值對應的階次作為模型階次。對輸出響應信號降噪，基于奇異值分解方法應用最廣。研究表明，奇異值分解具有理想的去相關特性，基于奇異值分解的信號分析方法可對信號進行重構，能較好從背景噪聲中分離出有用信號［9］；但利用奇異值分解的降噪方法通常只能進行一次截斷奇異值分解（TSVD）［10－12］，不能獲得最佳降噪效果。理論上對響應信號降噪技術，屬線性數學中矩陣低秩逼近（Low Rank Approximation）范疇，通過獲得Frobenius范數意義的最佳逼近矩陣降低噪聲［13］。低秩逼近可解決信號過程及圖像增強中減噪問題。所謂結構低秩逼近指保持矩陣結構情況下對其秩進行約束。本文提出低秩Hankel矩陣逼近方法降低響應信號中噪聲干擾。利用結構測量脈沖響應信號構造Hankel矩陣，并進行低秩逼近計算后獲得降噪信號，并用其進行模態參數識別。通過數值算例研究測量噪聲、矩陣維數等對低秩Hankel矩陣逼近方法影響，利用某懸臂梁實驗數據驗證本方法的有效性。

1 低秩Hankel矩陣逼近方法

1.1 Hankel矩陣及奇異值分解

一個N自由度動力系統脈沖響應函數可表示為

實測脈沖響應函數h（t）中含未知的M階模態，并以采樣間隔Δt表示成離散形式時，式（1）可表示為

式中：M≤N；l＝0，1，2，…。

用該實測脈沖響應信號hl構建m×n維Hankel矩陣，得

式中：m，n為 Hankel矩陣行、列數，m，n≥2M；s＝m＋n－2。

對Hm×n進行奇異值分解，得

式中：上標“T”表示矩陣轉置；U，V為正交矩陣；∑為對角矩陣，對角元素為降序排列的奇異值。而∑可分解為r個非零奇異值子矩陣∑r及幾個零子矩陣，即

式（5）表明矩陣Hm×n的秩為r。理論上超出矩陣秩的奇異值應為零。實測信號因隨機噪聲影響，奇異值并不等于零，但會變得很?。?4］。

1.2 低秩逼近

當實測脈沖響應序列hl受隨機噪聲干擾時可寫為

式中：l，el為真實信號及噪聲。

理論上，由含噪信號hl構建的Hankel矩陣H可分為兩部分，即

式中：為真實信號構建的Hankel矩陣；E為噪聲矩陣。已證明，若信號l包含M階模態，則矩陣的秩等于2M［14］。即矩陣的秩即為模型階次，為信號中包含模態數的兩倍。

低秩Hankel矩陣逼近降噪方法的基本思想為：基于H獲，即通過與H最接近的Hankel矩陣＾（秩為2M）逼近，使矩陣H與之差的 Frobenius范數最小。具體步驟為：① 對Hankel矩陣H進行截斷奇異值分解，即H＝U∑VT，基于確定的模型階次，即矩陣的秩r獲得獲得低秩逼近矩陣。此時＾不滿足Hankel矩陣形式。② 將矩陣A＾中各元素由其所在反對角線上元素的數學平均值代替，獲得滿足 Hankel形式的矩陣 H＾。H＾的秩不為 r。交替迭代步驟①、②，直至滿足收斂標準。

1.3 關鍵點

對含噪信號用低秩Hankel矩陣逼近降噪的兩關鍵點在于如何確定Hankel矩陣維數即m，n及如何保留由真實信號產生的奇異值分離階數，即如何確定模型階次r。本文通過數值算例研究矩陣維數選擇對計算效率、降噪效果影響規律，并基于奇異值分解方法確定模型階次。

2 數值算例

2.1 數據模擬

建立5自由度質量－彈簧－阻尼系統數值模型見圖1。單元質量、剛度、阻尼系數分別為 mn＝50 kg，kn＝2.9×107N／m，cn＝1 000 N·s／m。經特征值分析，得模態頻率理論值為 34.499 Hz，100.70 Hz，158.730 Hz，203.880 Hz，232.520 Hz；模態阻尼比的理論值為 0.003 737 4，0.010 909，0.017 197，0.022 092，0.025 198。

圖1 五自由度質量－彈簧－阻尼系統Fig.1 A 5－DOF massspringdashpot system

在任一自由度輸入單位脈沖激勵可得任一自由度輸出的脈沖響應函數。用Matlab編制程序可得25個脈沖響應函數。本例采用第1自由度輸入、第1自由度輸出的脈沖響應函數，1 024個采樣點，采樣頻率500 Hz，經傅里葉變換可得頻率響應函數。通過對精確（不含噪聲）脈沖響應函數疊加高斯白噪聲模擬含噪的脈沖響應函數。噪聲水平用百分比定量描述，該百分比定義為白噪聲標準差與精確信號標準差之比。精確信號與含5%、20%噪聲信號頻率響應函數見圖2，圖中5個峰值分別對應5階模態頻率，第5階模態（頻率232.520 Hz）相對較弱，易受噪聲干擾。

圖2 精確信號與含5%、20%噪聲信號對比Fig.2 The comparison of noise free signal and signal with 5%and 20%noise

2.2 確定模型階次

用基于奇異值的相對變化率［8］確定模型階次。利用結構脈沖響應信號構造Hankel矩陣，進行奇異值分解后計算相鄰奇異值相對變化率。變化率最大值即為對應模型的階次。不含噪聲情況下對1 024個采樣點構建Hankel矩陣。多次計算發現矩陣維數只要滿足行數、列數大于模型階次即可，即m，n＞10。故構建Hankel矩陣H513×512，并進行奇異值分解，計算奇異值相對變化率－模型階次指標，見圖3（a）。很明顯，模型階次指標最大值對應10階。對含5%、20%噪聲脈沖響應信號，模型階次計算結果見圖3（b）、（c）?？梢钥闯?，隨噪聲增強，模型階次最大值指標減小。5%噪聲時可確定模型階次為10；20%噪聲時模型階次為8，說明信號中某階模態被噪聲湮沒。由于第5階模態對脈沖響應貢獻較小，更易受噪聲影響。噪聲嚴重時第5階模態完全被噪聲湮沒，故模型階次由10降為8。

圖3 不同噪聲的模型階次指標Fig.3 Model order indicators with different noise level

2.3 信號降噪

含噪信號確定模型階次后采用低秩Hankel矩陣逼近方法進行降噪的另個關鍵點在于如何確定矩陣的維數。以含5%噪聲信號為例，理論上對秩約束為10的Hankel矩陣低秩逼近計算，矩陣維數只需滿足行數、列數大于10。因此矩陣最小維數為H1014×11，最大維數為H513×512。當選擇矩陣最小維數n＝11（n為矩陣列數）時，用低秩 Hankel矩陣逼近方法分別經100、1 000、10 000次迭代，獲得降噪效果見圖4。由圖4看出，盡管低頻區噪聲可有效消除掉，但高頻區噪聲即使經10 000次迭代仍明顯存在。而采用矩陣的最大維數n＝512時，僅經2次迭代即達較好降噪效果，見圖5。

圖4 n＝11時經不同迭代次數后降噪效果Fig.4 Performance of noise elimination for theHankel matrix with n＝11

為量化降噪效果，引入作為評價指標。其中F，F＾分別為精確及降噪后頻響函數，2代表Frobenius范數。實際上α在數值上等于已定義的噪聲水平，含5%噪聲信號的α值為0.05。對應不同維數矩陣（n＝11、30、100、512），α與迭代次數關系見圖6。由圖6看出，n＝512時只需2次迭代α即達到水平漸近線值0.007；n＝11時即使經10 000次迭代α為0.01，仍大于n＝512時經1次迭代的α值；n＝512時降噪效果最好，n＝11時降噪效果最差。

用同一臺計算機，n＝11、512的計算用時與迭代次數比較見表1。由表1看出，n＝512的單次迭代用時為n＝11的30倍，達到收斂所需迭代次數n＝11是n＝512的約10000倍。因此，n＝512的計算效率更高。需說明的是，圖6中未列出其它維數矩陣降噪情況。實際上，經多次計算發現，約n＝300時曲線與n＝512的曲線幾乎重合，考慮計算用時，n＝300的計算效率更高，可認為n＝300為更好選擇；但為降噪算法應用簡便，對任意個數的脈沖響應信號，仍建議選方陣或接近方陣的Hankel矩陣形式。

圖5 n＝512時2次迭代后降噪效果Fig.5 Performance of noise elimination for the Hankel matrix with n＝512 after 2 iterations

圖6 不同維數矩陣α與迭代次數關系Fig.6αagainst the number of iterations for Hankel matrices

表1 不同維數Hankel矩陣計算效率比較Tab.1 Comparison of the computational time for Hankel matrices

2.4 模態參數識別

對降噪后脈沖響應信號采用單參考點復指數法進行模態參數識別［1］。識別的模態頻率、阻尼比分別見表2、表3。與理論值、含噪信號識別值對比看出，采用降噪后信號識別的模態頻率、阻尼比與理論值非常接近，模態頻率相對誤差不超過0.1%，阻尼比前4階最大相對誤差為2.725%，第5階阻尼比識別誤差較大，因第5階模態貢獻最小，最易受噪聲影響，含噪信號無法識別第5階模態頻率及阻尼比。前4階識別結果相對誤差亦較大。

表2 用降噪信號與含5%噪聲信號識別模態頻率 （Hz）Tab.2 The estimated modal frequencies from the filtered and 5%corrupted signals（Hz）

表3 用降噪信號與含5%噪聲信號識別阻尼比Tab.3 The estimated damping ratios from the filtered and 5%corrupted signals

2.5 含20%噪聲信號

由以上知，信號中第5階模態完全被噪聲湮沒，確定模型階次為8，構建Hankel矩陣H513×512，通過秩約束為8的矩陣低秩逼近計算獲得降噪信號；但第5階模態經降噪處理后丟失。通過對降噪及含噪信號分別進行模態識別，結果見表4、表5。對比二表發現，采用降噪信號識別的前4階模態頻率與理論值較接近，模態頻率相對誤差范圍為0～0.358%，阻尼比相對誤差范圍1.109%～8.705%。相反，采用含噪信號識別的模態頻率及阻尼比的相對誤差均較大。

表4 用降噪信號及含20%噪聲信號識別的模態頻率 （Hz）Tab.4 The estimated modal frequencies from the filtered and 20%corrupted signals（Hz）

表5 用降噪信號及含20%噪聲信號識別阻尼比Tab.5 The estimated damping ratios from the filtered and 20%corrupted signals

3 懸臂梁模型試驗

為驗證本文所提方法的有效性，用懸臂梁進行沖擊實驗。布置傳感器測量振動響應，見圖7。采樣頻率400 Hz。采用上端兩傳感器測量脈沖響應為例，說明本方法的應用步驟。對兩組響應各取512個數據點分別構建Hankel矩陣H257×256，通過計算模型階次指標得兩組響應信號模型階次均為6，見圖8。分別對基于秩約束為6的兩矩陣H257×256進行低秩逼近計算，獲得降噪信號，見圖9。由圖9看出，實測信號中毛刺（噪聲）已被去除。對降噪信號進行模態參數識別，并與實測信號識別結果對比，見表6、表7。盡管無法獲得實驗模型真實模態參數，但可通過對比不同信號識別結果的一致性評價。由降噪信號識別的模態參數一致性較好，其中 3階模態頻率相對誤差分別為 0.001%，0.002%，0.007%；阻尼比相對誤差分別為 0.018%，0.060%，0.553%，均遠小于實測信號識別值的相對誤差。

圖7 懸臂梁物理模型及傳感器布置Fig.7 A cantilever beam experimental setup and one closeup accelerometer

圖8 兩組信號模型階次指標Fig.8 Model order indicators for measured signals at two locations

圖9 兩組信號降噪前后對比Fig.9 Performance of noise elimination for signals related to sensor 1 and sensor 2

表6 用實測信號與降噪信號識別的模態頻率 （Hz）Tab.6 The modal frequencies estimated from the measured and filtered signals（Hz）

表7 用實測信號與降噪信號識別的阻尼比Tab.7 The damping ratios estimated from the measured and filtered signals

4 結 論

（1）利用結構的實測脈沖響應信號構造Hankel矩陣，提出利用低秩Hankel矩陣逼近降噪，采用降噪信號進行模態參數識別。經數值算例與模型實驗已驗證本方法的有效性。

（2）隨噪聲增強信號中較弱模態易受干擾，模型階次的確定亦會受影響。

（3）矩陣列數最小時，即使經多次迭代仍無法獲得較好降噪效果；矩陣為方陣或接近方陣時，只需較少幾次迭代即可獲得較好降噪效果。建議使用該方法時用方陣或接近方陣的Hankel矩陣形式。

［1］Ewins D J.Modal testing：theory，practice and applications（2nd edition）［M］.England：Baldock， Hertfordshire，Research Studies Press，2000.

［2］Allemang R J，Brown D L.A unified matrix polynomial approach to modal identification［J］.Journal of Sound and Vibration，1998，211（3）：301－322.

［3］常軍，張啟偉，孫利民.穩定圖方法在隨機子空間識別模態參數中的應用［J］.工程力學，2007，24（2）：39－44.CHANG Jun，ZHANG Qiwei，SUN Limin.Application of stabilization diagram for modal paramenter identification using stochastic subspace method［J］.Engineering mechanics，2007，24（2）：39－44.

［4］易偉建，劉翔.動力系統模型階次的確定［J］.振動與沖擊，2008，27（11）：12－16.YI Weijian，LIU Xiang.Order identification of dynamic system model［J］.Journal of Vibration and Shock，2008，27（11）：12－16.

［5］Hu SL J，Bao X X，Li H J.Improved polyreference time domain method for modal identification using local or global noise removal techniques［J］.Sci ChinaPhys Mech Astron，2012，55：1464－1474.

［6］趙學智，葉邦彥，陳統堅.奇異值差分譜理論及其在車床主軸箱故障診斷中的應用［J］.機械工程學報，2010，46（1）：100－108.ZHAO Xuezhi，YE Bangyan，CHEN Tongjian.Difference spectrum theory of singular value and its application to the fault diagnosis of headstock of lathe［J］.2010，46（1）：100－108.

［7］包興先，劉福順，李華軍，等.復指數方法降噪技術及其試驗研究［J］.中國海洋大學學報，2011，41（1／2）：155－160.BAO Xingxian，LIU Fushun，LI Huajun，et al.The complex exponential method based on singalnoise separation for modal analysis［J］.Periodical of ocean university of China，2011，41（1／2）：155－160.

［8］王樹青，林裕裕，孟元棟，等.一種基于奇異值分解技術的模型定階方法［J］.振動與沖擊，2012，31（15）：87－91.WANG Shuqing，LIN Yuyu，MENG Yuangdong，et al.Model order determination based on singular value decomposition［J］.Journal of Vibration and Shock，2012，31（15）：87－91.

［9］齊子元，米東，徐章隧，等.奇異譜分析在機械設備故障診斷中的應用［J］.噪聲與振動控制，2008，28（1）：82－85.QI Ziyuan，MI Dong，XU Zhangsui，et al.Application of singular spectral analysis in mechanical device fault diagnosis［J］.Nosie and Vibration Control，2008，28（1）：82－85.

［10］段向陽，王永生，蘇永生.基于奇異值分解的信號特征提取方法研究［J］.振動與沖擊，2009，28（11）：30－33.DUAN Xiangyang，WANG Yongsheng，SU Yongsheng.Feature extraction methods based on singular value decomposition［J］.Journal of Vibration and Shock，2009，28（11）：30－33.

［11］徐鋒，劉云飛.基于中值濾波－奇異值分解的膠合板拉伸聲發射信號降噪方法研究［J］.振動與沖擊，2011，30（12）：135－140.XU Feng，LIU Yunfei.Noise reduction of acoustic emission signals in a plywood based on median filteringsingular value decomposition［J］.Journal of Vibration and Shock，2011，30（12）：135－140.

［12］錢征文，程禮，李應紅.利用奇異值分解的信號降噪方法［J］振動、測試與診斷，2011，31（4）：459－463.QIAN Zhengwen， CHENG Li， LI Yinghong. Noise reduction using singular value decomposition［J］.2011，31（4）：459－463.

［13］De Moor B.Total least squares for affinely structured matrices and the noisy realization problem［J］.IEEE Trans Signal Process，1994，42（11）：3104－3113.

［14］Hu SL J，Bao X X，Li H J.Model order determination and noise removal for modal parameter estimation［J］.Mechanical Systems and Signal Processing，2010，24（6）：1605－1620.