包興先，李昌良，劉志慧

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

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

1 低秩Hankel矩陣逼近方法

1.1 Hankel矩陣及奇異值分解

一個N自由度動力系統(tǒng)脈沖響應(yīng)函數(shù)可表示為

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

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

用該實測脈沖響應(yīng)信號hl構(gòu)建m×n維Hankel矩陣，得

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

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

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

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

1.2 低秩逼近

當(dāng)實測脈沖響應(yīng)序列hl受隨機噪聲干擾時可寫為

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

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

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

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

1.3 關(guān)鍵點

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

2 數(shù)值算例

2.1 數(shù)據(jù)模擬

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

圖1 五自由度質(zhì)量－彈簧－阻尼系統(tǒng)Fig.1 A 5－DOF massspringdashpot system

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

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

2.2 確定模型階次

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

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

2.3 信號降噪

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

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

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

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

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

圖6 不同維數(shù)矩陣α與迭代次數(shù)關(guān)系Fig.6αagainst the number of iterations for Hankel matrices

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

2.4 模態(tài)參數(shù)識別

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

表2 用降噪信號與含5%噪聲信號識別模態(tài)頻率 （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階模態(tài)完全被噪聲湮沒，確定模型階次為8，構(gòu)建Hankel矩陣H513×512，通過秩約束為8的矩陣低秩逼近計算獲得降噪信號；但第5階模態(tài)經(jīng)降噪處理后丟失。通過對降噪及含噪信號分別進行模態(tài)識別，結(jié)果見表4、表5。對比二表發(fā)現(xiàn)，采用降噪信號識別的前4階模態(tài)頻率與理論值較接近，模態(tài)頻率相對誤差范圍為0～0.358%，阻尼比相對誤差范圍1.109%～8.705%。相反，采用含噪信號識別的模態(tài)頻率及阻尼比的相對誤差均較大。

表4 用降噪信號及含20%噪聲信號識別的模態(tài)頻率 （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 懸臂梁模型試驗

為驗證本文所提方法的有效性，用懸臂梁進行沖擊實驗。布置傳感器測量振動響應(yīng)，見圖7。采樣頻率400 Hz。采用上端兩傳感器測量脈沖響應(yīng)為例，說明本方法的應(yīng)用步驟。對兩組響應(yīng)各取512個數(shù)據(jù)點分別構(gòu)建Hankel矩陣H257×256，通過計算模型階次指標(biāo)得兩組響應(yīng)信號模型階次均為6，見圖8。分別對基于秩約束為6的兩矩陣H257×256進行低秩逼近計算，獲得降噪信號，見圖9。由圖9看出，實測信號中毛刺（噪聲）已被去除。對降噪信號進行模態(tài)參數(shù)識別，并與實測信號識別結(jié)果對比，見表6、表7。盡管無法獲得實驗?zāi)Ｐ驼鎸嵞B(tài)參數(shù)，但可通過對比不同信號識別結(jié)果的一致性評價。由降噪信號識別的模態(tài)參數(shù)一致性較好，其中 3階模態(tài)頻率相對誤差分別為 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 兩組信號模型階次指標(biāo)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 用實測信號與降噪信號識別的模態(tài)頻率 （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 結(jié) 論

（1）利用結(jié)構(gòu)的實測脈沖響應(yīng)信號構(gòu)造Hankel矩陣，提出利用低秩Hankel矩陣逼近降噪，采用降噪信號進行模態(tài)參數(shù)識別。經(jīng)數(shù)值算例與模型實驗已驗證本方法的有效性。

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

（3）矩陣列數(shù)最小時，即使經(jīng)多次迭代仍無法獲得較好降噪效果；矩陣為方陣或接近方陣時，只需較少幾次迭代即可獲得較好降噪效果。建議使用該方法時用方陣或接近方陣的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］常軍，張啟偉，孫利民.穩(wěn)定圖方法在隨機子空間識別模態(tài)參數(shù)中的應(yīng)用［J］.工程力學(xué)，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］易偉建，劉翔.動力系統(tǒng)模型階次的確定［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］趙學(xué)智，葉邦彥，陳統(tǒng)堅.奇異值差分譜理論及其在車床主軸箱故障診斷中的應(yīng)用［J］.機械工程學(xué)報，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］包興先，劉福順，李華軍，等.復(fù)指數(shù)方法降噪技術(shù)及其試驗研究［J］.中國海洋大學(xué)學(xué)報，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］王樹青，林裕裕，孟元棟，等.一種基于奇異值分解技術(shù)的模型定階方法［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］齊子元，米東，徐章隧，等.奇異譜分析在機械設(shè)備故障診斷中的應(yīng)用［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］徐鋒，劉云飛.基于中值濾波－奇異值分解的膠合板拉伸聲發(fā)射信號降噪方法研究［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］錢征文，程禮，李應(yīng)紅.利用奇異值分解的信號降噪方法［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.