基于Tikhonov正則化及奇異值分解的載荷識別方法

郭 榮, 房懷慶, 裘 剡, 于欽林, 朱偉偉

(1.同濟大學(xué) 中德學(xué)院，上海 201804; 2.同濟大學(xué) 汽車學(xué)院,上海 201804; 3.同濟大學(xué) 汽車新能源汽車工程中心,上海 201804)

傳遞路徑分析方法[1-3](Transfer Path Analysis, TPA)為機械系統(tǒng)振動噪聲常用診斷手段，廣泛用于汽車、航空航天、潛艇等領(lǐng)域。TPA主要通過估計噪聲源及各傳遞路徑的貢獻率，找出對響應(yīng)影響較大的聲源及傳遞路徑，其中載荷識別精度對各傳遞路徑貢獻量影響顯著。

載荷識別為結(jié)構(gòu)動力學(xué)難點之一，而頻響函數(shù)求逆法在工程實際中使用較廣。Schevenels等[4-6]利用結(jié)構(gòu)動力學(xué)關(guān)系構(gòu)建載荷識別動力學(xué)方程，并提出通過頻響函數(shù)求逆法進行載荷識別，用奇異值分解法(Singular Value Decomposition，SVD)改善聲源識別中頻響函數(shù)求逆的病態(tài)問題。Choi等[7-8]運用Tikhonov正則化方法，提高載荷識別結(jié)果穩(wěn)定性，并分析對比不同正則化參數(shù)選擇方法對載荷識別精度影響。目前研究采用SVD法、Tikhonov正則化法雖可提高載荷識別精度，并未比較兩方法間的優(yōu)缺點及系統(tǒng)病態(tài)性不同時兩種方法的適用性。對此，本文提出基于Tikhonov正則化與奇異值分解的載荷識別方法，考慮系統(tǒng)病態(tài)性對載荷識別影響，用頻響函數(shù)法矩陣條件數(shù)評價系統(tǒng)的病態(tài)性，通過平板模型仿真分析，在頻域運用奇異值分解法與基于普通交叉驗證法(Ordinary Cross Validation，OCV)、廣義交叉驗證法(Generalized Cross Validation，GCV)、L曲線法(L-curve)三種不同正則化參數(shù)選擇方法的Tikhonov正則化法進行載荷識別，分析總結(jié)系統(tǒng)病態(tài)性不同時，不同方法載荷識別精度及適用性。通過簡支梁臺架載荷識別實驗，驗證該方法的可行性。

1 結(jié)構(gòu)載荷識別方法

TPA方法假設(shè)機械系統(tǒng)某目標響應(yīng)位置的結(jié)構(gòu)噪聲可由各傳遞路徑產(chǎn)生的貢獻分量相加[1]：

(1)

式中：p為目標響應(yīng)點處總聲壓；Hij為傳遞路徑i在j方向頻率響應(yīng)函數(shù)(Frequency Response Function，F(xiàn)RF)；Fj為系統(tǒng)中第j位置施加的激勵。由式(1)知，若求出目標點總聲壓，則須先獲得各條路徑的頻響函數(shù)及工作載荷，其中結(jié)構(gòu)載荷識別精度對各路徑聲壓貢獻量影響較大。

工程中常采用頻響函數(shù)求逆法識別工作載荷。設(shè)待識別工作載荷數(shù)目為s，實測響應(yīng)點位置數(shù)為n，響應(yīng)點數(shù)n與工作載荷數(shù)目s通常不相等，故頻率響應(yīng)函數(shù)的逆為廣義逆[9]：

(2)

其中：a為響應(yīng)點響應(yīng);H+為頻響函數(shù)H的廣義逆。

由式(2)，頻響函數(shù)求逆法應(yīng)用的關(guān)鍵在于解決其可能存在的病態(tài)逆問題。若系統(tǒng)病態(tài)性嚴重，即使較小測量誤差也會被放大。該病態(tài)問題可用方法以改善，其中奇異值分解法因計算簡便應(yīng)用最廣[4-6]；而Tikhonov正則化法求解穩(wěn)定性良好，亦獲得應(yīng)用[7-8]。

2 奇異值分解法

(3)

將式(3)代入式(2)可得奇異值分解法對載荷識別的表達式[5]：

F=VS-1UHa

(4)

3 正則化方法

3.1 Tikhonov正則化理論

測量方法所限，所測響應(yīng)a及頻率響應(yīng)函數(shù)H均存在誤差。為控制及度量誤差，設(shè)整系統(tǒng)絕對誤差為：

e=a-HF

(5)

為使誤差e盡量小，Tikhonov引入罰函數(shù)[10]：

J=(eHe)+λ(FHF)

(6)

當(dāng)J對F的一階導(dǎo)為零時，誤差e存在最小值，則可得F：

F=(HHH+λI)-1HHa

(7)

其中：I為單位陣；λ為正則化參數(shù)。λ等于零時，上述解即為最小二乘解；λ趨向無窮大時，解為零向量。如何選擇正則化參數(shù)λ是正則化算法的關(guān)鍵。正則化參數(shù)λ選擇方法主要有普通交叉驗證法、廣義交叉驗證法及L曲線法。

3.2 正則化參數(shù)選擇方法

3.2.1 OCV法

在OCV方法中，通過式(7)求得激振力F，并代入式a=HF中求得響應(yīng)，比較所求響應(yīng)與實際響應(yīng)間差別，OCV可定義為[8]：

(8)

3.2.2 GCV法

定義對角化循環(huán)酉矩陣W：

(9)

用循環(huán)酉矩陣W與某矩陣相乘，即對該矩陣進行離散傅里葉變換，可用此解決FRF對角化問題。由絕對誤差：

e=a-HF=a-USVHF

(10)

式(10)兩邊同乘WUH，得：

et=at-HtF

(11)

其中：et=WUHe；at=WUHa；Ht=WSUH。對式(11)用OCV法理論可得GCV法表達式為[6,11]：

(12)

3.2.3L曲線法

(13)

4 平板模型仿真分析

4.1 仿真模型建立

研究在系統(tǒng)病態(tài)性情況不同時，采用基于OCV法、GCV法、L曲線法不同正則化參數(shù)選擇的Tikhonov方法及SVD法的載荷識別精度及此方法具體的適用性，利用Patran、Nastran進行仿真分析。仿真中，用600×500×1.5 mm鋼板，彈性模量E=2.07×1011N/m2，泊松比μ=0.3，結(jié)構(gòu)阻尼0.06。在板上隨機選擇激勵及響應(yīng)點位置，盡量避免邊界位置[4]。由于鋼板厚度相對長寬非常小，故可將問題弱化為平面問題[13]，采用殼單元，單元厚度1.5 mm。簡支矩形板約束為：[Ux,Uy,Uz]=[0,0,0], [Rx,Ry,Rz]=[0,0,1](“0”為約束該自由度；“1”為不約束該自由度，在Patran中設(shè)置為空格)，網(wǎng)格大小控制為0.01 m，共有節(jié)點3 111個、單元3 000個。對有預(yù)約束的平板模型進行模態(tài)分析，獲得平板前10階固有頻率，見表1。通過4次在不同位置施加單位激勵進行頻響分析，可確定激勵點及響應(yīng)點間頻響函數(shù)。第5次在不同位置同時施加一定大小激勵，進行頻響分析，獲得響應(yīng)點加速度響應(yīng)用于識別所加激勵。仿真中施加激振力及響應(yīng)點位置見表2?？紤]結(jié)構(gòu)載荷頻率范圍，頻響分析范圍設(shè)為10～250 Hz。

表1 平板前10階固有頻率

表2 激勵幅值、坐標及響應(yīng)位置坐標

4.2 噪聲模型

響應(yīng)點加速度測量及激振力施加會存在噪聲干擾，因此需在仿真數(shù)據(jù)中引入不同信噪比(SNR)噪聲。噪聲用可加性高斯白噪聲，其模型為：

N(ω)=Nndej2πNud

(14)

其中：Nnd為正態(tài)分布隨機數(shù)，均值0，方差A(yù)(ω)10-(B/20)(A(ω)為加速度響應(yīng)，B為信噪比)；Nud為0～1均勻分布隨機數(shù)?？紤]實際情況，在頻響函數(shù)中分別加入40 dB，22 dB信噪比等級噪聲，在響應(yīng)中分別加入40 dB，22 dB，15 dB，10 dB信噪比噪聲，探討不同噪聲等級下基于不同正則化參數(shù)選擇的Tikhonov方法與SVD法的載荷識別精度。

4.3 不同噪聲等級下系統(tǒng)病態(tài)性

頻響函數(shù)測量誤差會直接影響系統(tǒng)的病態(tài)性，而在數(shù)據(jù)處理中病態(tài)性危害較大，會造成結(jié)果失真，若系統(tǒng)嚴重病態(tài)，即使測量中極小誤差也會在數(shù)據(jù)處理過程中被放大。一般采用頻響函數(shù)法矩陣條件數(shù)度量系統(tǒng)的病態(tài)性。設(shè)頻率響應(yīng)函數(shù)矩陣法矩陣HHH的條件數(shù)為K，統(tǒng)計經(jīng)驗表明[14]：01 000，系統(tǒng)嚴重病態(tài)。計算在不同噪聲等級下，系統(tǒng)病態(tài)情況，結(jié)果見圖1。由圖1看出，① 無噪聲干擾時頻響函數(shù)條件數(shù)波峰恰好與結(jié)構(gòu)固有頻率對應(yīng)，最大條件數(shù)位置為系統(tǒng)一階固有頻率；② 隨頻響函數(shù)中噪聲等級增大，條件數(shù)波動越劇烈，尤其一階固有頻率處。

4.4 采用不同逆問題解法載荷識別結(jié)果

為討論系統(tǒng)病態(tài)性不同時采用不同逆問題解法的載荷識別結(jié)果差異，以頻響函數(shù)中含22 dB噪聲、響應(yīng)中含15 dB噪聲為例，用不同逆問題解法的載荷識別結(jié)果見圖2。

圖1 系統(tǒng)在不同噪聲等級下病態(tài)情況

由圖1、圖2看出：

(1) 誤差較大處均在條件數(shù)>1 000區(qū)域，即系統(tǒng)病態(tài)區(qū)域，且均在系統(tǒng)固有頻率處，其中第一階固有頻率誤差最大(圖1)。

圖2 頻響函數(shù)中含22 dB噪聲、響應(yīng)中含15 dB噪聲時載荷識別結(jié)果

(2) 在一階固有頻率附近，用不同正則化參數(shù)選擇的Tikhonov正則化法載荷識別精度普遍優(yōu)于SVD法(圖2)。

(3) 在100～250 Hz頻率范圍內(nèi)，條件數(shù)相對較小，與實際激勵相比，SVD法較Tikhonov正則化方法識別的載荷精度更好。

為定量分析不同載荷識別方法誤差，定義誤差計算表達式為：

(15)

其中：FRf為識別激勵；FTf為施加激勵；N為采集數(shù)據(jù)點數(shù)；n為待識別激勵數(shù)目；εforce單位dB。以頻響函數(shù)法矩陣條件數(shù)1 000為基準，條件數(shù)>1 000時，SVD法及基于不同正則化參數(shù)選擇的Tikhonov方法在不同噪聲等級下載荷識別誤差見圖3；條件數(shù)≤1 000時，SVD法與基于不同正則化參數(shù)選擇的Tikhonov方法在不同噪聲等級下載荷識別誤差見圖4。

由圖3、圖4看出：

(1) 條件數(shù)>1 000時，Tikhonov正則化方法相較SVD法可有效降低載荷識別誤差，其中基于OCV法與L曲線法正則化參數(shù)選擇的Tikhonov方法載荷識別誤差相差不大，而GCV法誤差最小。

圖4 條件數(shù)小于等于1 000時SVD法和基于不同正則化參數(shù)選擇的Tikhonov方法的識別誤差

(2) 條件數(shù)≤1 000時，基于三種正則化參數(shù)選擇的Tikhonov方法載荷識別誤差相差不大，而SVD法除響應(yīng)中含10 dB高噪聲外，載荷識別誤差最小。響應(yīng)中10 dB信噪比噪聲在實際中較高，高噪聲時SVD法的識別精度無Tikhonov正則化方法好。

5 基于Tikhonov法及SVD法的載荷識別方法

由以上分析知，系統(tǒng)嚴重病態(tài)(條件數(shù)>1 000)時，Tikhonov正則化法識別載荷精度較好，其中GCV選擇正則化參數(shù)精度最高；條件數(shù)≤1 000時，SVD法識別載荷誤差較小。為此，提出基于Tikhonov正則化方法與SVD法的載荷識別新方法。該方法以頻響函數(shù)法矩陣條件數(shù)為依據(jù)，條件數(shù)≤1 000時，采用SVD法識別載荷；條件數(shù)>1 000時，采用Tikhonov正則化方法識別載荷，以GCV法選擇正則化參數(shù)為例，對基于Tikhonov法與SVD法的載荷識別流程進行說明，見圖5。

圖5 方法流程圖

采用新的基于Tikhonov正則化方法與SVD法的載荷識別方法與單獨用SVD法和Tikhonov正則化識別誤差對比見圖6。由圖6看出：①響應(yīng)中含40 dB、22 dB、15 dB噪聲時，新方法明顯較單獨用SVD法與Tikhonov正則化方法識別載荷精度高；②響應(yīng)中含10 dB高噪聲時，新方法載荷識別精度明顯優(yōu)于SVD法，但較Tikhonov方法略差。原因在于高噪聲時SVD識別誤差較大，導(dǎo)致新方法識別精度下降。與單獨用SVD法及Tikhonov方法相比，新方法在不同等級噪聲下均可取得較好載荷識別精度。

圖6 SVD法、Tikhonov正則化方法與新方法載荷識別誤差

6 載荷識別實驗驗證

實驗用兩端簡支矩形梁，采用YE15400電動式激振器激勵同時三點激勵，激勵信號為隨機觸發(fā)，測量四點加速度響應(yīng)，數(shù)采LMS SCADAS Ⅲ 316，分析范圍25～800 Hz。簡支梁長、寬、高為：640 mm、56 mm、8 mm，模態(tài)分析獲得簡支梁前四階固有模態(tài)見表3。共進行三組實驗，實驗框圖見圖7，實驗臺架見圖8。

表1 簡支梁前四階固有頻率

圖7 實驗框圖

圖8 實驗臺架

實驗獲得頻響函數(shù)法矩陣條件數(shù)見圖9。由圖9看出，第一階固有頻率處條件數(shù)最大。

圖9 頻響函數(shù)法矩陣條件數(shù)

圖10 SVD法、Tikhonov正則化方法與新方法載荷識別誤差

采用新的基于Tikhonov正則化方法與SVD法載荷識別方法與單獨用SVD法及Tikhonov正則化識別誤差對比見圖10。由圖10看出，新方法可提高載荷識別精度。

7 結(jié) 論

本文采用Tikhonov正則化與SVD法解決病態(tài)逆問題。通過仿真討論系統(tǒng)病態(tài)性情況不同時，采用基于OCV法、GCV法、L曲線法不同正則化參數(shù)選擇的Tikhonov方法及SVD法載荷識別精度具體適用性，據(jù)頻響函數(shù)法矩陣條件數(shù)，提出綜合使用Tikhonov正則化法及SVD法的新載荷識別方法，結(jié)論如下：

(1) 載荷識別誤差較大處均在條件數(shù)>1 000處，即系統(tǒng)病態(tài)區(qū)域，且均在系統(tǒng)固有頻率處，其中第一階固有頻率誤差最大；

(2) 條件數(shù)>1 000，相比SVD法，采用Tikhonov方法的載荷識別誤差較小，其中采用GCV法選擇正則化參數(shù)精度最高；當(dāng)條件數(shù)≤1 000時，SVD法載荷識別精度較好；

(3) 由仿真、實驗研究發(fā)現(xiàn)，與單獨采用SVD法及Tikhonov方法相比，綜合使用Tikhonov正則化法與SVD法可有效提高載荷識別精度。

本文研究結(jié)論以期為動力學(xué)病態(tài)逆問題解決提供參考。

參 考 文 獻

[1]郭 榮,萬 鋼,趙燕男,等.車內(nèi)噪聲傳遞路徑分析方法探討[J].振動、測試與診斷,2007(3):199-203.

GUO Rong, WAN Gang, ZHAO Yan-nan,et al. Study on transfer path analysis method of automobile interior noise[J]. Journal of Vibration, Measurement and Diagnosis, 2007(3): 199-203.

[2]Vecchio A, Valent L, Carluccio L,et al. Experimental noise transfer path analysis on the helicopter A109[C]. National Aerospace Laboratory NLR-32nd European Rotorcraft Forum, ERF,Aerospace Societies,2007:197-210.

[3]Pieter S,Van der Linden P J G. An inverse force measurement method to determine the injected structure-borne sound power from an installation into a building element[J]. Building Acoustics, 2010, 17(3):199-219.

[4]Kim Y, Nelson P A. Optimal regularization for acoustic source reconstruction by inverse methods[J]. Journal of Sound and Vibration , 2004, 275(4-5):463-487.

[5]Thite A N, Thompson D J. The quantification of structure-borne transmission paths by inverse methods[J]. Journal of Sound and Vibration, 2003, 264(3):411-451.

[6]Nelson P A,Yoon S H. Estimation of acoustic source strength by inverse methods[J]. Journal of Sound and Vibration, 2000, 233(4): 643-705.

[7]Choi H G,Thite A N, Thompson D J. A threshold for the use of Tikhonov regularization in inverse force determination[J]. Applied Acoustics, 2006, 67(7):700-719.

[8]Choi H G, Thite A N,Thompson D J. Comparison of methods for parameter selection in Tikhonov regularization with application to inverse force determination[J]. Journal of Sound and Vibration,2007, 304(3-5):894-917.

[9]Janssens M H A, Verheij J W. A pseudo-forces methodology to be used in characterization of structure-borne sound sources [J].Applied Acoustics, 2000, 61(3):285-308.

[10]Tikhonov A N, Arsenin V Y. Solutions of Ill-posed problems [M]. New York: Wiley, 1977.

[11]Golub G H, Heath M, Wahba G. Generalized cross-validation as a method for choosing a good ridge parameter[J]. Technometrics, 1979,21(2): 215-223.

[12]Hansen P C. Analysis of discrete ill-posed problems by means of the L-curve [J].SIAM Review, 1992, 34(4):561-580.

[13]張 方,秦遠田.工程結(jié)構(gòu)動載荷識別方法[M].北京：國防工業(yè)出版社,2011.

[14]王振杰.測量中不適定問題的正則化解法[M].北京：科學(xué)出版社,2006.