環境激勵下基于信號降噪的模態參數識別研究

包興先

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

基于環境激勵的結構模態參數識別是當前結構健康監測與安全評估的主要研究內容之一。對于大型結構，如橋梁、海洋平臺等，進行在線實時的模態參數識別往往只能依靠環境荷載的激勵。與傳統的模態參數識別方法相比，它不需要額外的激振設備，不影響結構物的正常使用，并且操作方便，測試費用較低，具有廣闊的發展空間和應用前景。在工程實際中，輸入荷載(環境荷載)是未知的，往往只能獲得含噪聲的響應信號，如何有效獲得結構物的模態參數是研究的關鍵［1－2］。

研究表明，將環境激勵近似當作白噪聲處理是可行的，而自然激勵技術(NExT)對白噪聲的處理是非常有效的。復指數(Prony)法是目前較完善的模態參數識別方法之一。NExT與Prony結合的時域聯合算法為識別結構模態參數提供了可行性［3］。由于測量的信號不可避免地受到噪聲的干擾，識別結果通常包含虛假模態，如何辨識并剔除虛假模態成為參數識別過程的關鍵問題，直接影響了參數識別技術在工程實際中應用的效果。目前模態參數識別的通常做法是計算采用的模型階次高于真實的模型階次，從而引入了“噪聲模態”［3］。為了區分真實模態和虛假模態，大多采用穩定圖方法。穩定圖表示的是模態參數與模型階次的關系，理論上隨模型階次的增加，真實的模態參數會趨于穩定而虛假模態卻不會，但實際上識別結果中真實模態參數也會存在不穩定性，導致穩定圖法有時也會失效［4－5］。尤其是當信號的信噪比低的時候，如何區分大量的虛假模態和真實模態將變得很困難。

與傳統的模態參數識別方法不同，本文提出了環境激勵下基于信號降噪的模態參數識別方法。該方法首先對測量的隨機響應信號采用NExT處理獲得互相關函數，然后對其進行模型定階及結構矩陣低秩逼近(Structured Low Rank Approximation，SLRA)計算，獲得降噪信號，最后利用Prony方法對降噪信號進行模態參數識別。文中將通過一個四自由度的質量—彈簧—阻尼系統數值算例和導管架式海洋平臺物理模型實驗驗證該方法的有效性。

1 NExT

James等［6］指出環境激勵下結構兩點之間響應的互相關函數和脈沖響應函數有相似的表達式，在求得兩點之間響應的互相關函數后，可運用時域模態識別方法識別結構的模態參數。

一個N自由度線性阻尼結構的動力學方程表示如下:

其中，［M］、［C］、［K］分別是質量、阻尼和剛度矩陣;{x··(t)］、{x·(t)］、{x(t)］是加速度、速度和位移向量;{f(t)］是作用在系統上的輸入力向量。對于穩態的互不相關的輸入力，兩個測試的位移之間的互相關函數能夠表示如下:式中Φmr是振型矩陣，Gnr是一個同參考點n和模態階數 r有關的常數項，mr是第 r階模態質量，ξr、ωnr、ωdr是第r階模態的阻尼比、固有頻率和阻尼頻率。而脈沖響應函數能夠表示為:

可以看出，互相關函數Rmn(τ)與系統的脈沖響應函數Xml(t)有相似的形式，都能表示成一系列衰減正弦函數的和，即每個衰減正弦都有一個固有頻率和阻尼比同結構的各階模態相對應，振型向量的位置也相同。因此，可以將兩點間響應的互相關函數Rmn(τ)代替系統的脈沖響應函數Xml(t)。

2 Prony方法

對脈沖響應函數式(3)進行整理，可以得到表達式(4):

如果一個N自由度動力系統的實測脈沖響應函數中包含未知的M階模態，而且，當h(t)以采樣間隔Δt表示成離散形式時，式(4)可表示為:

由于式(5)是2M階線性微分方程齊次解的一般形式，那么一定存在一個與脈沖響應序列hq相關的線性差分方程:

如果式(7)括號內的每一項等于零，那么式(7)恒成立。而此時，Vr(r=1，…，2M)是系數為βm的2M階多項式的根，即

其中，H∈R2M×2M是方陣形式的Hankel矩陣。由方程(10)，我們可得到未知的多項式系數:

如果方程(11)包含多于2M個方程，則用最小二乘法來求解該超定方程的根βm。

3 結構矩陣低秩逼近(SLRA)

低秩逼近是解決信號過程和圖像增強中減噪問題的一種普遍工具。理論上，對于脈沖響應信號的降噪技術，屬于線性數學中的矩陣低秩逼近范疇，通過獲得Frobenius范數意義下的最佳逼近矩陣來降低噪聲［7］。如今，在許多實際應用領域中提出了更高的要求，在保持矩陣結構的情況下，還對其秩作出了約束。從而提出了一類新問題，即結構矩陣低秩逼近［2，8］。所謂“結構”，指矩陣的特定形式，如Hankel、Toeplitz等。

當脈沖響應函數hq中包含未知的M階模態，采用hq構建m×n維的Hankel矩陣，得到:

式中，m，n分別為Hankel矩陣的行數和列數，m，n≥2M，s=m+n－2。

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

式(13)中，上標“T”表示矩陣轉置;矩陣U和V是正交矩陣;Σ是對角矩陣，其對角元素為降序排列的奇異值。而Σ可分解為k個非零奇異值子矩陣Σk和幾個零子矩陣:

這一分解表明矩陣Hm×n的秩是k。理論上，那些超出矩陣的秩的奇異值應當等于零。對于實測信號，由于隨機噪聲的影響，這些奇異值并不等于零，但是會變得很小［9］。

若hq受到隨機噪聲的干擾時，則可以寫成:

式中和eq分別代表真實信號和噪聲。理論上，由式(15)中含噪信號hq構建的Hankel矩陣H可以分為兩部分:

式中，代表真實信號構建的Hankel矩陣，E代表噪聲矩陣。已經證明，如果信號包含有M階模態，則矩陣的秩等于 2M［2］。

結構矩陣低秩逼近降噪算法的基本思想為:基于H得到H—，即通過與H最接近的Hankel矩陣(秩為2M)來逼近，使得矩陣H和之差的Frobenius范數最小。降噪步驟如下:

步驟2，將矩陣中的各元素由其所在反對角線上的元素的數學平均值代替，得到滿足Hankel形式的矩陣。此時的秩不為k。

步驟1和2交替迭代，直到滿足收斂標準。

4 數值算例

4.1 白噪聲激勵

建立一個四自由度的質量－彈簧－阻尼系統的數值模型如圖1所示。單元的質量、剛度和阻尼系數分別為 m1，2，3，4=50 kg、k1，2，3=2 × 107N/m、k4=2 × 106N/m 、c1，2，3，4=500 N·s/m。通過特征值分析，可得到模態頻率的理論值為:26.457 Hz、53.135 Hz、127.05 Hz、181.68 Hz;模態阻尼比的理論值為:0.013 243、0.018 167、0.012 377、0.014 760。

圖1 四自由度質量－彈簧－阻尼系統Fig.1 A 4 － DOF mass-spring-dashpot system

采用Matlab編制程序，以均值為零、方差為1的高斯白噪聲為輸入，將其右向加載在m1上。采樣頻率500 Hz，得到 m1、m2兩處隨機響應信號，約 3 s，見圖 2。以m1處響應信號為參考測點，采用NExT方法計算m2測點與參考測點的互相關函數，取500個數據點，見圖3。

圖2 白噪聲激勵下兩處響應信號Fig.2 Response signals related to two locations with white noise excitation

圖3 由NExT方法得到的互相關函數圖Fig.3 Cross-correlation function obtained by NExT

4.2 噪聲模擬

實測信號不可避免地會受到各種噪聲的影響。通過對精確的(不含噪聲的)隨機響應信號疊加高斯白噪聲來模擬含噪聲的隨機響應信號。噪聲水平通過一個百分比來定量描述，該百分比定義為白噪聲的標準差和精確信號的標準差之比。圖4為精確信號和含10%、40%噪聲信號分別由NExT方法得到互相關函數，再經傅里葉變換后的頻域圖。可以看出，在隨機響應信號的基礎上即使疊加了10%的白噪聲，經過NExT處理后，噪聲信號與精確信號基本吻合;當噪聲水平較高，比如40%時，噪聲會對原始信號形成較大干擾。這表明，NExT方法對隨機白噪聲有一定的抑制作用。

圖4 不同精確度互相關函數頻域對比圖Fig.4 Frequency domain of exact and corrupted cross-correlation functions

4.3 信號降噪

對含10%、40%噪聲的互相關函數分別構建維數為251!250的Hankel矩陣。采用結構矩陣低秩逼近方法，首先需要確定矩陣的秩，即模型階次。對Hankel矩陣進行奇異值分解，采用奇異值歸一化曲線的方法來確定模型階次［2，10］。如圖5所示，曲線突降到漸近線時對應的奇異值個數即為模型階次。可以看出，含10%、40%噪聲的矩陣的秩都為8，也就是說，信號中都包含4階模態。確定模型階次后，對含10%、40%噪聲的Hankel矩陣分別進行結構低秩逼近計算，得到降噪的互相關函數。圖6為含40%噪聲的信號經過降噪后與精確信號的對比。可以看出，噪聲被較好地消除了。

圖5 不同噪聲情況下的模型定階曲線Fig.5 Model order indicators with different noise level

圖6 含40%噪聲信號降噪后與精確的互相關函數頻域對比圖Fig.6 Frequency domain of filtered(40%corrupted)and exact cross-correlation functions

4.4 模態參數識別

采用Prony方法分別對精確、含噪及降噪的互相關函數進行模態參數識別。對含10%噪聲及其降噪信號的識別結果見表1和表2，對含40%噪聲及其降噪信號的識別結果見表3和表4。可以發現，由精確信號識別的模態參數與該系統的理論值較為接近。由于含噪信號是在精確信號的基礎上疊加一定程度的白噪聲，所以通過比較含噪信號、降噪信號與精確信號的識別結果可以說明本文提出方法的有效性。從表1和表2可以發現，由含10%噪聲信號識別的頻率的誤差范圍為0.53% －1.48%，而由降噪信號識別的頻率的誤差均小于1%;由含10%噪聲信號識別的阻尼比的誤差較大，其范圍為13.7% －75.6%，而由降噪信號識別的阻尼比的誤差相對較小，范圍為3.46% －9.83%。同樣，從表3和表4可以看出，由降噪信號識別的頻率和阻尼比的誤差范圍分別為0.41% －1.26%和12.1% －24.4%，均小于由含40%噪聲識別的頻率和阻尼比的誤差。

表1 采用含10%噪聲信號和降噪信號識別的模態頻率(Hz)Tab.1 The estimated modal frequencies from the 10%corrupted and filtered signals(Hz)

表2 采用含10%噪聲信號和降噪信號識別的阻尼比(%)Tab.2 The estimated damping ratios from the 10%corrupted and filtered signals(%)

表3 采用含40%噪聲信號和降噪信號識別的模態頻率(Hz)Tab.3 The estimated modal frequencies from the 40%corrupted and filtered signals(Hz)

表4 采用含40%噪聲信號和降噪信號識別的阻尼比(%)Tab.4 The estimated damping ratios from the 40%corrupted and filtered signals(%)

5 模型實驗

建立一鋼質導管架式海洋平臺物理模型，主腿尺寸為Φ25 mm×2.5 mm，橫撐及斜撐尺寸均為16 mm×1.5 mm，模型示意圖如圖7。將模型置于振動臺上，在各層關鍵節點布置X向加速度傳感器，通過振動臺施加X向的白噪聲激勵，采樣頻率500 Hz，得到各節點處隨機響應信號。以節點17的響應信號為參考測點，采用NExT方法分別計算節點25、27、29、31與參考測點的互相關函數。以節點25為例，圖8為截取其中2 s的隨機響應信號及其互相關函數。

圖7 導管架式海洋平臺模型Fig.7 Jacket offshore platform model

圖8 節點25處的隨機響應信號及互相關函數圖Fig.8 Response signals related to node 25 with white noise excitation and cross－correlation function

采用奇異值歸一化曲線的方法確定模型階次為4(如圖9)，即信號中包含2階模態。對實測信號的互相關函數進行SLRA計算，得到降噪的互相關函數。圖10為實測信號與降噪信號的互相關函數頻域對比圖，可以看出，噪聲被較好地消除了。

圖9 節點25處互相關函數的模型定階曲線Fig.9 Model order indicator related to cross-correlation function of node 25

圖10 實測信號與降噪信號的互相關函數頻域對比圖Fig.10 Frequency domain of filtered and measured cross-correlation functions

采用Prony方法分別對實測和降噪的互相關函數進行模態參數識別。另外，對節點27、29、31的實測信號分別采用相同的處理步驟，不再贅述。由四個不同位置處的實測及降噪信號識別的頻率和阻尼比分別列于表5和表6。盡管無法得到物理模型的真實模態參數，但是可以通過比較不同位置信號的識別結果的一致性來評價識別效果。從表5和6中可以發現，由不同位置處的降噪信號識別的模態頻率和阻尼比都非常一致，而由實測信號的識別結果則有較大差異，特別是阻尼比的識別結果有較大誤差。

表5 采用實測信號和降噪信號識別的模態頻率(Hz)Tab.5 The estimated modal frequencies from the measured and filtered signals(Hz)

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

6 結論

(1)本文提出一種在環境激勵下基于測試信號降噪的模態參數識別方法，即NExT－SLRA－Prony方法。該方法首先對測試的隨機響應數據采用NExT計算而獲得互相關函數，進而基于SLRA方法得到降噪后的信號，最后通過Prony方法識別結構的模態參數。數值算例和模型實驗結果表明，該方法對測量信號有很好的降噪作用，識別精度高。

(2)NExT對隨機白噪聲有一定的抑制作用。含噪的隨機響應信號經NExT處理后得到互相關函數，基于此識別的頻率值受噪聲影響較小，但阻尼比受噪聲影響相對較大。

［1］ Wang Shu-qing，Liu Fu-shun.New accuracy indicator to quanfity the true and false modes for eigensystem realization algorithm［J］.Structural Engineering and Mechanics，2010，34(5):625－634.

［2］ 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.

［3］ Ewins D J.Modal Testing:Theory，Practice and applications［M］.2nd Edition， Baldock， Hertfordshire， England:Research Studies Press，2000.

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

［5］常軍，張啟偉，孫利民.穩定圖方法在隨機子空間識別模態參數中的應用［J］.工程力學，2007，24(2):39－44.CHANG Jun，ZHANG Qi-wei，SUN Li-min.Application of stabilization diagram for modal paramenter identification using stochastic subspace method［J］.Engineering Mechanics，2007，24(2):39－44.

［6］JamesⅢ，G H，Carne T G，Lauffer J P.The natural excitation technique(NExT)for modal parameter extraction from operating structures［J］.The International Journal of Aanalytical and Experimental Modal Analysis，1995，10(4):260－277.

［7］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.

［8］ Hu S L J，Bao X X，Li H J.Improved polyreference time domain method for modal identification using local or global noise removal techniques［J］. Science China:Physics Mechanica Astronomy，2012，55:1464－1474.

［9］包興先，劉福順，李華軍，等.復指數方法降噪技術及其試驗研究［J］.中國海洋大學學報，2011，41(1):155－160.BAO Xing-xian，LIU Fu-shun，LI Hua-jun，et al.The complex exponential method based on singal-noise separation for modal analysis［J］.Periodical of Ocean University of China，2011，41(1):155 －160.

［10］王樹青，林裕裕，孟元棟，等.一種基于奇異值分解技術的模型定階方法［J］.振動與沖擊，2012，31(15):87－91.WANG Shu-qing，LIN Yu-yu，MENG Yuang-dong，et al.Model order determination based on singular value decomposition［J］.Journal of Vibration and Shock，2012，31(15):87－91.