Tikhonov正則化在運行工況傳遞路徑分析的應用*

成 瑋， 盧英英， 陸建濤， 張周鎖

(1.西安交通大學機械制造系統工程國家重點實驗室 西安，710049)(2.西安交通大學機械工程學院 西安，710049)

Tikhonov正則化在運行工況傳遞路徑分析的應用*

成 瑋1,2， 盧英英2， 陸建濤2， 張周鎖1,2

(1.西安交通大學機械制造系統工程國家重點實驗室 西安，710049)(2.西安交通大學機械工程學院 西安，710049)

針對傳統運行工況傳遞路徑分析(operational transfer path analysis，簡稱OTPA)存在的不足，通過理論和試驗分析，提出基于Tikhonov正則化方法的OTPA反問題模型。首先，分析Tikhonov正則化方法的理論優勢，給出Tikhonov正則化參數選擇的依據，同時調節電機轉速獲得不同運行工況數據，利用奇異值分解方法研究殼體結構的振動傳遞路徑，分析傳統OTPA算法總貢獻量誤差及路徑貢獻量估計精度；其次，分析運行工況數據是否滿足Picard條件，提出基于Tikhonov正則化方法的OTPA算法，并分析Tikhonov正則化參數對所提出算法的影響。分析結果表明，所提出的方法顯著減小了總貢獻量和路徑貢獻量誤差以及路徑誤判現象。該研究可為振動噪聲監控與減振降噪提供理論依據。

運行工況傳遞路徑分析； Tikhonov正則化； 奇異值分解； Picard條件； 減振降噪

引 言

水下航行器的聲隱身性能是衡量其安全性和作戰能力的重要指標[1]，在低、中速航行時，動力機械設備振動是產生輻射噪聲的主要來源[2]。轎車和高速列車的振動噪聲是評價車輛性能的重要指標[3]。因此，機械設備振動噪聲的有效監測與控制對于提高裝備性能具有重要工程意義。

運行工況傳遞路徑分析是一種用于識別和定量分析機械系統振動噪聲傳播路徑的方法[4]。由于建?？焖俑咝?，OTPA廣泛應用于實際工程[5]。Roozen等[6]基于奇異值分解，建立了齒輪箱的OTPA模型并研究了振動傳遞路徑。Putner等[7]利用OTPA成功估計了不同聲源對車外噪聲的貢獻量。Yan等[8]基于局部奇異值分解，研究了板輻射噪聲對車內噪聲的影響。傳統OTPA采用奇異值分解[9]，獲取求解傳遞率函數的主成分，并建立傳遞路徑模型[10]。由于系統模態會產生路徑交叉耦合使系統病態，外加測試噪聲影響，傳統OTPA不可避免出現貢獻量過估計現象和誤判現象[11]。Tikhonov正則化方法[12]充分考慮了不適定問題[13]求解的不穩定性和噪聲污染等因素，廣泛應用于反問題求解[14-15]。OTPA模型屬于反問題，即通過運行工況數據識別系統的傳遞率函數矩陣。當前，利用Tikhonov正則化的OTPA方法研究很少。

筆者提出了基于Tikhonov正則化的OTPA反問題模型。首先，給出了傳統OTPA和基于Tikhonov正則化OTPA的基本理論；其次，分析了加筋圓柱殼試驗臺(用于模擬水下航行器的動態特性)運行工況數據的相關性，建立了試驗臺的傳統OTPA模型并研究路徑貢獻量，分析了奇異值和矩陣條件數對貢獻量準確度的影響；最后，檢驗了Picard條件，分析了正則化參數對基于Tikhonov正則化OTPA模型精度的影響。通過與傳統OTPA模型比較，驗證了基于Tikhonov正則化OTPA模型在路徑貢獻量識別中具有更高的精度。

1 OTPA基本理論

1.1 傳統OTPA

OTPA[4]線性系統為

(1)

其中：T為未知傳遞率函數矩陣；X為已知激勵源參考點試驗工況矩陣；Y為已知目標點試驗工況矩陣；m為參考點個數；n為目標點個數；r為試驗工況個數。

由于OTPA通過試驗工況響應信號估計傳遞率函數矩陣，因此對于式(1)：a.為了保證參考點試驗工況矩陣的可逆性，試驗工況個數必須等于或大于參考點個數，即r≥m;b.結構模態引起的激勵在其他路徑上產生響應(路徑耦合)，容易出現路徑誤判。

1.2 傳統OTPA算法

傳統OTPA方法通過求解最小二乘問題估計傳遞率函數矩陣，即

(2)

具體算法如下：

1) 對X奇異值分解

(3)

2) 確定X截斷奇異值

(4)

3) 估計傳遞率函數矩陣

(5)

4) 分析傳遞路徑

Y′=X′T

(6)

其中：X′為參考點實際工況數據；Y′為傳統OTPA識別的目標點振動響應。

1.3 基于Tikhonov正則化的OTPA

針對傳統OTPA不足，式(2)僅考慮解的擬合程度。式(3)貢獻量為工程經驗值(τ≥85%)。Tikhonov正則化廣泛應用于反問題研究[14-15]，對于線性系統

Ax=b

(7)

其中：Am×n(m≥n)為已知系數矩陣；bm×1為已知系數向量；xn×1為未知向量。

在實際工程中，噪聲干擾會導致式(7)的解失去物理意義。Tikhonov正則化通過修改系數矩陣或約束未知向量獲得最優解。

(8)

其中：‖‖為Euclidean范數；λ>0為正則化參數。

對比式(2)和式(8)，式(8)同時考慮解的擬合程度和穩定性，因此基于Tikhonov正則化的OTPA方法具有理論優勢。

式(8)的解可表達為

(9)

結合式(3)可得

(10a)

(10b)

離散傅里葉系數和奇異值是Picard條件的主要構成要素，且離散傅里葉系數比奇異值趨于零的速度快。若式(7)滿足離散Picard條件[16]，則能夠求其具有物理意義的正則逼近解。式(8)的關鍵是選擇正則化參數λ，即通過λ平衡解的擬合程度和穩定性，Tikhonov正則化參數的選擇方法主要有廣義交叉驗證(generalized cross validation，簡稱GCV)方法[17]和L曲線方法[18]。

GCV通過求解式(11)極小值確定參數λ，即

(11)

(12)

根據式(11)和式(12)，GCV方法和L曲線方法的基本思想都是確定的參數λ使殘量Axλ-b和解xλ的范數同時保持在較小的水平上，即平衡解xλ的擬合程度和穩定性。

1.4 基于Tikhonov正則化的OTPA算法

算法如下：

1) 融合式(7)和式(1)，建立OTPA模型；

2) 融合式(3)和式(9)，推導出式(10)；

4) 分析傳遞路徑

Y′=X′Tλ

(13)

其中：X′為參考點實際工況數據；Y′為基于Tikhonov正則化OTPA識別的目標點振動響應。

2 傳統OTPA試驗研究

2.1 試驗臺及試驗工況簡介

根據某型號艦艇殼體結構，搭建了圓柱加筋殼體試驗臺，主要由薄壁圓柱殼體結構及7個均勻分布加強筋組成。遠離加強筋一端布置有大小兩臺偏心振動電機模擬振源，通過殼體結構不同位置傳感器獲取結構振動加速度響應信號。試驗臺主要設備實物圖與測點布置圖分別如圖1和圖2所示。

圖1 試驗臺及主要設備實物圖Fig.1 The physical diagram of test bed and major equipments

圓圈內數字為傳感器編號；傳感器1，2，3分別為小電機的豎直方向、水平方向和基座附近；傳感器4，5，6分別為大電機的豎直方向、水平方向和基座附近圖2 傳感器空間位置示意圖(單位：mm)Fig.2 Positions of the sensors (unit: mm)

大、小電機參考點分別為傳感器4與傳感器1，目標點為傳感器13，電機不同轉速下的運行工況數據如表1所示。

當試驗運行工況數據相關性較大時，式(5)中小奇異值較多，導致傳遞率函數矩陣估計不準確，影響OTPA模型的精度。傳感器1和4不同運行工況數據的相關系數如圖3所示。不同工況下傳感器1，4和13的相關系數如圖4所示。

表1 不同電機轉速的運行工況

圖3 不同工況數據相關性Fig.3 The correlation of data in different operating conditions

圖4 不同傳感器數據相關性Fig.4 The correlation of data between different sensors

由圖3可得，傳感器4不同工況的相關系數小于0.1，表明參考點4不同工況的相關性很弱，可保證參考點試驗運行工況矩陣求逆良態。

由圖4可得，不同工況下傳感器1和4的相關系數約為0.15～0.4，表明大小電機耦合較弱，參考點1和4分別表示小電機與大電機的振源特性。不同工況下傳感器1和13的相關系數大部分高于0.80，傳感器4和13的相關系數約為0.5～0.7，表明傳感器13的振動來源于小電機和大電機，且與小電機相關性更高。

綜上分析，參考點4不同工況的相關性很弱，可保證矩陣良態；參考點1和4的相關系數較小，表明振源耦合較弱。綜合分析可得參考點和運行工況的選擇較為合理，可用于研究運行工況振動路徑傳遞。

2.2 傳統OTPA

OTPA基本思路：利用已知運行工況數據識別傳遞率函數矩陣，分析任意運行工況下振動噪聲傳遞路徑。利用表1的運行工況數據，根據式(1)可得試驗臺OTPA模型為

(14)

用傳統OTPA算法計算的式(14)參考點矩陣奇異值和條件數分別如圖5和圖6所示。

圖5 試驗臺OTPA矩陣奇異值Fig.5 Matrix singular values of test bench OTPA

圖5表明試驗臺OTPA模型有兩個奇異值，最大奇異值和最小奇異值分別為0.02和3.64×10-19，前者代表小電機振動，后者代表大電機振動，即小電機振動強度較大，與圖4結論吻合。

圖6 試驗臺OTPA矩陣條件數Fig.6 Matrix condition number of test bench OTPA

圖6表明試驗臺OTPA模型矩陣條件數很大，導致不同頻率下試驗臺OTPA模型不穩定。導致試驗臺OTPA模型不穩定因素有：測試噪聲和圖3所示參考點1不同工況相關性(相關系數約為0.7～0.9，即相關性很大)。

傳統OTPA計算傳遞率函數如圖7所示。T11為參考點1與傳感器13的傳遞率函數，T12為參考點4與傳感器13的傳遞率函數。

圖7 傳統OTPA傳遞率函數Fig.7 Traditional OTPA transmissibility functions of the test bed

傳統OTPA算法識別的大小電機與目標點傳遞率函數波動較大，小電機對應的傳遞率函數相比波動稍小。分析認為，參考點不同工況的相關性和測試噪聲導致矩陣條件數很大，使試驗臺OTPA模型不穩定，最終導致傳遞率函數波動較大；小電機對應的奇異值較大，受噪聲干擾較小，因此小電機對應的傳遞率函數較穩定。因此，圖7與圖3～6的結論一致。

假設實際工況為大小電機分別以工作頻率11.8 Hz和15.2 Hz運行。傳統OTPA算法計算式(14)傳感器13的總貢獻量、大電機路徑貢獻量和小電機路徑貢獻量分別如圖8～10所示。

圖8 總貢獻量(傳統OTPA)Fig.8 The total contributions (traditional OTPA)

由圖8可知，15.2 Hz時傳統OTPA算法總貢獻量與試驗值一致，而11.8 Hz時相差較大。

圖9 大電機路徑貢獻量(傳統OTPA)Fig.9 The big motor path contribution (traditional OTPA)

由圖9可得，大電機貢獻量試驗值約為0.000 3g，而傳統OTPA值約為0.002g，即傳統OTPA算法大電機路徑貢獻量誤差很大。

圖10 小電機路徑貢獻量(傳統OTPA)Fig.10 The small motor path contribution (traditional OTPA)

小電機工作頻率為15.2 Hz，由圖10可得，11.8 Hz時傳統OTPA算法小電機路徑貢獻量較大，出現嚴重誤判。由于小電機工作頻率為15.2 Hz，且小電機對應傳遞率函數較穩定，所以15.2 Hz時傳感器13總貢獻量和小電機路徑貢獻量均與試驗值吻合較好。

通過以上分析，傳統OTPA可定性識別振動傳遞路徑，即傳感器13振動能量源于工作頻率分別為11.8 Hz和15.2 Hz振源。穩定性較差的大電機傳遞率函數導致11.8 Hz時出現大電機路徑貢獻量誤差很大(圖9)和小電機路徑貢獻量誤判(圖10)。

傳遞率函數識別方法是OTPA模型的關鍵技術，但是容易受到參考點不同工況相關性和測試噪聲影響。大小電機分別以工作頻率11.8 Hz和15.2 Hz運行時，參考點和目標點信號相關系數如表2所示。傳感器1和4、1和13、4和13的相關系數分別為0.328，0.908和0.571，表明參考點1和4的相關性較弱，大小電機耦合較小，可分別表示振源特性；傳感器13振動主要源于小電機；分析結論與圖3和圖4結論一致。

表2 參考點和目標點的相關系數

Tab.2 Coefficients between reference points and target points

相關系數傳感器1傳感器4傳感器13傳感器110.3280.908傳感器40.32810.571傳感器130.9080.5711

綜上分析，噪聲干擾和參考點不同工況的相關性導致傳統OTPA算法精度較低，在與目標點相關性較弱振源對應頻率下，總貢獻量識別誤差較大，路徑貢獻量失真較嚴重。

3 Tikhonov正則化OTPA模型試驗研究

3.1 Picard條件檢驗

Picard條件是檢驗試驗臺OTPA模型是否具有物理意義Tikhonov正則化解的重要條件，試驗臺OTPA模型Picard條件如圖11所示。奇異值1和2曲線分別位于對應離散傅里葉系數曲線的上方，表明試驗臺OTPA模型滿足Picard條件，即可通過Tikhonov正則化方法求解其有物理意義的正則解。

圖11 試驗臺OTPA模型的Picard條件Fig.11 Picard conditions of OTPA model for the test bed

3.2 基于Tikhonov正則化的OTPA

通過式(10)中不同Tikhonov正則化參數λ識別試驗臺OTPA傳遞率函數如圖12和圖13所示。圖12和圖13表明傳遞率函數受正則化參數λ的影響較大。同時λ越小，傳遞率函數越接近傳統OTPA算法識別的傳遞率函數；λ越大，傳遞率函數出現明顯峰值，且過濾掉表1大小電機工作頻率以外其他頻段的干擾噪聲。

圖12 λ對小電機相應傳遞率函數的影響Fig.12 λ effects on transmissibility functions of small motor

圖13 λ對大電機傳遞率函數的影響Fig.13 λ effects on transmissibility functions of big motor

基于Tikhonov正則化OTPA算法的分析結果如圖14～16所示。

圖14 λ對總貢獻量的影響Fig.14 λ effects on total contributions

圖15 λ對大電機路徑貢獻量的影響Fig.15 λ effects on path contributions of big motor

圖14為不同λ條件下總貢獻量曲線。顯然，λ對總貢獻量影響很大。圖15為不同λ條件下大電機路徑貢獻量曲線。當λ=1×10-3和λ=1×10-4時，大電機路徑貢獻量出現誤判，但當λ=3×10-5時，路徑誤判消除。

圖16 λ對小電機路徑貢獻量的影響Fig.16 λ effects on path contributions of small motor

圖16為不同λ條件下小電機路徑貢獻量曲線，當λ=1×10-3和λ=1×10-4時，小電機路徑貢獻量出現誤判，但當λ=3×10-5時，路徑誤判消失。

綜上分析，應用Tikhonov正則化方法建立OTPA模型具有實用性和有效性；Tikhonov正則化參數λ嚴重影響總貢獻量和路徑貢獻量的分析結果；通過調節Tikhonov正則化參數λ可消除路徑誤判。

結合圖8～10，圖14～16，傳統OTPA算法和基于Tikhonov正則化方法OTPA算法分析結果的相對誤差對比如表3所示。

表3 傳統OTPA與Tikhonov正則化誤差

Tab.3 Errors of traditional OTPA and Tikhonov regularization

算法相對誤差/%f/Hz11.8f/Hz15.2大電機貢獻量小電機貢獻量λ=1×10-373.621.8482.54誤判24.3誤判λ=1×10-47.141.364.72誤判20.92誤判λ=3×10-581.381.36207.78.13傳統OTPA227.211.46635.753.03誤判

表3中，當λ=1×10-4時，與傳統OTPA算法分析結果比較，Tikhonov正則化方法將總貢獻量相對誤差分別從227.3%和1.46%減小到7.14%和1.36%。因此，Tikhonov正則化方法顯著提高了試驗臺OTPA模型總貢獻量識別精度(其中λ=1×10-4最優)。當λ=3×10-5時，與傳統OTPA算法分析結果相比，Tikhonov正則化方法將大電機路徑貢獻量相對誤差從635.75%減小到207.7%，同時消除小電機路徑貢獻量嚴重誤判現象，且小電機路徑貢獻量相對誤差為8.13% 。因此，Tikhonov正則化方法很好地克服了傳統OTPA路徑貢獻量誤差過大缺點，與此同時很好地彌補了傳統OTPA路徑貢獻量誤判缺陷。

綜上分析，Tikhonov正則化方法同時考慮了試驗臺OTPA模型傳遞率函數矩陣的擬合程度和穩定性，顯著提高了OTPA總貢獻量和路徑貢獻量識別準確度，且抑制了OTPA路徑貢獻量誤判缺陷。

4 結 論

1) 傳統OTPA算法僅考慮擬合程度，工程經驗值不足和測試噪聲等因素易導致路徑貢獻量和總貢獻量計算誤差較大，難以定量識別振動傳遞路徑。

2) 基于Tikhonov正則化方法OTPA同時考慮擬合程度和穩定性，顯著提高了路徑貢獻量和總貢獻量計算精度，同時有效克服了傳統OTPA路徑貢獻量誤判的不足。

3) 理論與試驗分析表明，基于Tikhonov正則化OTPA優于傳統OTPA，可實現振動傳遞路徑定量識別，為結構優化設計和減振降噪提供理論依據。

[1] 張磊, 曹躍云, 楊自春，等. 運行工況下艦船振動-聲傳遞路徑分析及試驗[J]. 華中科技大學學報:自然科學版, 2013, 41(2): 42-47.

Zhang Lei, Cao Yueyun, Yang Zichun, et al. Theoretical analysis and experimental validation of vibration-acoustic transferring paths of ships under operational conditions[J]. Journal of Huazhong University of Science & Technology: Natural Science Edition, 2013, 41(2):42-47. (in Chinese)

[2] 成瑋, 張周鎖, 何正嘉，等. 基于獨立分量分析的潛艇振源貢獻量定量計算方法[J]. 機械工程學報, 2010, 46(7): 83-87.

Cheng Wei, Zhang Zhousuo, He Zhengjia, et al. Quantitative calculation of vibration source contributions of submarines based on independent component analysis[J]. Journal of Mechanical Engineering, 2010, 46(7):83-87. (in Chinese)

[3] Li W, Wang D F, Chen S M, et al. Transfer path analysis of power train vibration on driver's noise[C]∥7th International Conference on Natural Computation (ICNC). Shanghai: IEEE, 2011: 2353-2357.

[4] Yoshida J, Noumura K. Method of transfer path analysis for interior vehicle sound using principal component regression method[J]. Honda R&D Technical Review, 2006,18(1):136-141.

[5] Janssens K, Mas P, Gajdatsy P, et al. A novel path contribution analysis method for test-based NVH trouble shooting[C]∥Proceedings of Acoustics 2008. Geelong, Victoria, Australia:[s.n.],2008: 3673-3683.

[6] Roozen N B, Leclere Q, Sandier C. Operational transfer path analysis applied to a small gearbox test set-up[C]∥Proceedings of Acoustics 2012. Nantes, France: [s.n.], 2012:3467-3473.

[7] Putner J，Lohrmann M，Fastl H. Contribution analysis of vehicle exterior noise with operational transfer path analysis[C]∥Proceedings of Meetings on Acoustics. Montreal, Canada:[s.n.], 2013,19,040035(1-9).

[8] Yan Li, Jiang Weikang. Research on the procedure for analyzing the sound quality contribution of sound sources and its application[J]. Applied Acoustics, 2014, 79:75-80.

[9] 胥永剛，謝志聰，孟志鵬，等. 基于奇異值分解的磁記憶信號特征提取方法[J]. 振動、測試與診斷,2014,34(6)：1106-1109.

Xu Yonggang, Xie Zhicong, Meng Zhipeng, et al. The signal feature extraction method of magnetic memory based on singular value decomposition[J]. Journal of Vibration, Measurement & Diagnosis, 2014, 34(6):1106-1109. (in Chinese)

[10]Klerk D, Ossipov A. Operational transfer path analysis: theory, guidelines and tire noise application[J]. Mechanical Systems and Signal Processing, 2010, 24(7): 1950-1962.

[11]Peter G, Karl J, Ludo G. Critical assessment of operational path analysis: mathematical problems of transmissibility estimation[J]. The Journal of the Acoustical Society of America, 2008, 123(5): 9359-9364.

[12]Hadamard J. Sur les problems aux derivees partielles et leur signification physique[R]. Bulletin: Princeton University, 1902.

[13]Tikhonov A N, Arsenin V Y. Solution of ill-posed problems [J]. Mathematics of Computation, 1978,32(14):491.

[14]Tikhonov A. Solution of incorrectly formulated problems and the regularization method[J]. Doklady Akademii Nauk SSSR, 1963,151(3):501-504.

[15]馬超，華宏星 . 正則化技術在狀態空間載荷識別中的應用[J]. 振動、測試與診斷, 2014, 34(6): 1154-1158.

Ma Chao, Hua Hongxing. The application of regularization technique for load identification in state space[J]. Journal of Vibration, Measurement & Diagnosis, 2014, 34(6): 1154-1158. (in Chinese)

[16]Hansen P C. The discrete picard condition for ill-posed problems[J]. BIT Numerical Mathematics, 1990,30(4):658-672.

[17]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.

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

10.16450/j.cnki.issn.1004-6801.2017.01.009

國家自然科學基金資助項目(51305329)；中國博士后科學基金資助項目(2014T70911，2013M532032)；高等學校博士學科點專項科研基金資助項目(20130201120040)；陜西省博士后基金資助項目

2015-02-09；

2015-04-20

TH128； TB535

成瑋，男，1983年6月生，副教授。主要研究方向為機械系統狀態監測與故障診斷。曾發表《A comprehensive study of vibration signals for a thin shell structure using enhanced independent component analysis and experimental validation》(《Journal of Vibration and Acoustics-Transactions of The ASME》2014, Vol.136, No.4)等論文。 E-mail: chengw@mail.xjtu.edu.cn