應用VMD與Teager能量算子的結構模態系統辨識*

靳 行， 林建輝

(西南交通大學牽引動力國家重點實驗室 成都，610031)

引 言

對于整體結構如齒輪箱、整備后的車輛等進行模態測試時，由于施加人工激勵較為困難，因此通常直接利用環境激勵下的輸出信號識別模態參數。與傳統的輸入、輸出的模態實驗分析相比，直接利用輸出信號識別模態參數具有測量簡單、對近頻與重頻敏感、更接近真實動力特性等優點[1]。

環境激勵下的模態參數研究除了常用的頻域法與時域法，時頻法也是近些年被國內外廣泛關注的方法。例如Bao等[2]開發了一種自適應的時頻方法來識別電纜的張力變化。Spiridonakos等[3]在實驗室建立了橋梁-車輛系統時變頻率的時變自回歸移動平均模型。經驗模態分解[4](empirical mode decomposition，簡稱EMD)通過將信號分解成多個單一頻率的固有模式函數(intrinsic mode functions，簡稱IMF)，然后通過希爾伯特變換提取瞬時頻率。陳雙喜等[5]應用該方法有效地提取車輛-軌道垂向耦合系統動態特性。李康強等[6]應用EMD與能量算子有效提取行星齒輪箱各階次的模態參數。然而在一些情況下，使用EMD方法提取的IMF由于模態混疊導致一系列頻率，即IMF非單一頻率。由于邊界效應，導致IMF本身可能是虛假的偽分量。

為了克服EMD的種種限制與缺陷，變分模態分解(variational mode decomposition，簡稱VMD)應運而生[7]。 VMD方法繼承了EMD的迭代思想，通過尋找中心頻率，將信號進分解成不同尺度的細節信號。已經有研究表明，VMD在機械振動信號分析[8]、故障診斷[9-10]、電網預測[11]和地震信號分析[12]等領域的性能優于EMD方法。盡管應用VMD在許多方面表現出比EMD 的性能更加優異，但VMD也存在一些缺點，如層數參量K與罰參量α需要提前給定，其取值大小對計算結果影響尚缺乏理論依據[8]。

在此基礎上，筆者提出基于VMD與TEO對輸出信號的結構模態識別方法，并針對VMD參數選擇這一缺陷，提出了通過調整層數參量K尋找穩定極點的作法，實現模態參數準確識別。首先，使用VMD算法對測量結構輸出的加速度信號進行分解；其次，應用TEO算法得到瞬時頻率與瞬時幅值來擬合固有頻率與阻尼比；最后，將層數參量K作為一種最優解參考，通過觀察不同VMD層數參量K時固有頻率與阻尼比組成的極點，判斷真實的模態參數。筆者分析了在沖擊力作用下結構的數值模型，對比現有的方法以驗證所提方法識別模態參數的準確性和有效性。實驗研究了動車組齒輪箱的錘擊法與環境激勵的振動數據，結果表明，所提方法在識別結構模態數據方面是有效而準確的。

1 基本理論

1.1 結構模態系統

具有n階自由度(degrees of freedom，簡稱DOF)的結構的運動方程可以寫成

(1)

每階模態固有頻率都是一個窄帶信號，由模態疊加原理可知，每個自由度均可通過 對n階模態固有頻率進行疊加，得到加速度響應

(2)

其中：Φi為第i階模態振型；qi(t)為第i階模態對應的響應函數。

當在第z個DOF上施加沖擊力，并將式(2)帶入式(1)中解耦，可以得到第i個廣義模態坐標上i的加速度響應

(3)

因此第p個自由度上的加速度響應可以寫成

(4)

將式(4)寫成調幅調制信號的形式

(5)

1.2 變分模態分解

VMD的方法原理與EMD方法十分相似，但是它放棄了EMD方法中循環篩選極值進行濾波的方法，而是將信號分解引入變分模型中來解決濾波問題。為了找到最優解，通過尋求約束變分模型，從而實現信號分解。在VMD分解的過程中，每個BIMF分量的中心頻率及其帶寬不斷地交替迭代，自適應地分解為合適的信號頻帶，得到K個預設尺度下的BIMF 分量。變分約束問題為

(6)

其中：uk為模態函數集；ωk為每個模態函數集的中心頻率；δ(t)為脈沖單位函數。

為了解決上述約束最優化問題，VMD 算法在實施過程中綜合利用了二次懲罰項和拉格朗日乘子法的優勢，引入了增廣拉格朗日函數ζ，如式(7)所示

ζ({uk},{ωk},λ)=

(7)

其中：α為罰參量；λ(t)為拉格朗日乘子。

解決式(7)中最小化問題的方法是：用交替方向乘子法進行一系列的迭代優化找到的極小值點。

VMD算法的分解詳細過程可參考文獻[7]。與EMD方法相似, 原始信號x(t)被分解為K個穩定的BIMF分量uk(t)及一個能量較低的殘差信號ε(t)

(8)

VMD通過在變分問題框架中迭代搜索變分模型最優解來實現信號的分解，其本質是維納濾波，具有很好的噪聲魯棒性，并且解決了EMD 算法存在的模態混疊、偽分量以及易受噪聲干擾等缺點，保留了EMD算法可以有效抑制ln[Ak(t)]交叉項的優點。

1.3 Teager 能量算子

Teager能量算子是一個非線性算子，為了減少離散TEO算法的誤差，筆者對TEO算子的求解進行了改進。采用改進后的TEO算法能準確地解調出原信號的瞬時頻率TEO[13]，信號u(t)的能量算子Ψ[·]運算定義為

(9)

對該無衰減自由振蕩的線性振子，其振動位移為x(t)=Acos(ωt+θ)。其中：A為振幅；ω為角頻率；θ為初相。

對上式求1階微分

(10)

改進的TEO算法是通過對Ψ[·]進行低通濾波，數字濾波器采用切比雪夫Ⅰ型，其單位脈沖響應為δ(n)，則有

Ψd[·]=δ(n)*Ψ[·]

(11)

用Ψd[·]替代Ψ[·]，便可得到改進的離散時間能量分離算法，即改進的離散TEO算法。由此得到瞬時包絡A(t)和瞬時角頻率ω(t)

(12)

1.4 模態固有頻率與阻尼比的識別

對于小阻尼系統，可以根據模態響應的衰減幅值來計算阻尼比與固有頻率[14]，由于模態響應與時間t之間存在以下關系

(13)

為了估計得到固有頻率與阻尼比，使用最小二乘法進行線性擬合，可以得到固有頻率ωn,k與阻尼比ξk。

2 分析步驟與方法

分析步驟與方法如下：

1) 為了減少計算量，避免弱衰減和零衰減部分對結果的影響，對原始信號截取明顯衰減的部分進行VMD 分解，得到多個BIMF 分量；

2) 預設最大分解層數Kmax，將VMD層數參量設為K=1,2,…,Kmax，依次對n個測點的時域信號進行VMD分解，得到n∑Kmax個BIMF分量；

3) 設定低通頻率，應用改進TEO算法得到瞬時幅值與瞬時頻率，根據式(13)擬合得到不同層數參量K下每一個BIMF固有頻率與阻尼比；

4) 根據層數參量K重新排列固有頻率與阻尼比，在穩態圖中觀察極點列，選擇不隨層數參量K變化的穩定極點作為最終識別結果。

其中穩態圖包含3種極點:

1) 首次出現的極點為原始極點;

2) 分解層數增加，原始極點的固有頻率變化小于閾值的極點為頻率極點；

3) 分解層數增加，原始極點或頻率極點的頻率變化與阻尼比變化均小于閾值的極點為穩態極點。

在穩態圖中選擇表現結構物理特性的極點，從一列極點的最低分解層數中找到穩定極點，與下一層極點比較，如果阻尼比較穩定，則選擇該極點的固有頻率與阻尼比作為最終識別結果。

3 仿真信號分析

3.1 仿真信號

圖1 五自由度模型Fig.1 The 5-dof numerical model

質量矩陣為

(14)

剛度矩陣為

(15)

仿真中使用Rayleigh阻尼模型，阻尼為C=0.5M+0.000 4K，Newmark算法用于獲得系統加速度曲線。采樣頻率為128Hz，沖擊力施加在第1個自由度

(16)

圖2 系統的時域自由振動響應和相應的功率譜密度Fig.2 Time-domain free-vibration responses of the system and the corresponding power spectral density

3.2 同類方法比較

首先對振動響應進行VMD分解，最大分解層數Kmax=10 ，罰參量設為α=1 024，統計不同層數參量下模態識別結果。

為了突出本研究方法的創新性，采用兩種同類方法進行對比。圖3所示為采用文獻[15]所述的Hilbert方法得到的穩態結果，圖4所示為采用文獻[16]所述的中心頻率方法得到穩態結果。圖5為本研究提出的TEO方法穩態結果。仿真信號中，頻率變化閾值為1%，阻尼比變化閾值為5%。

圖3 Hilbert模態參數極點穩態圖Fig.3 Modal parameter identification for Hilbert of stabilization diagra

如圖4～6所示，穩態圖中均存在5列明顯的極點列。表1是根據穩態圖中選擇得到的固有頻率。3種方法識別的固有頻率與理論值接近。

表2對比了識別結果與理論值的誤差，可以看出，筆者所提的TEO法識別的結果平均誤差要小于前2種方法，固有頻率識別的平均誤差降了0.14%。

圖4 中心頻率模態參數極點穩態圖Fig.4 Modal parameter identification for center frequency of stabilization diagra

圖5 TEO算法模態參數極點穩態圖Fig.5 Modal parameter identification for TEO of stabilization diagra

表3為3種方法識別的阻尼比結果，結果與理論值存在一定偏差，如表4所示，筆者所提方法平均誤差小于其他兩種方法。

表1 仿真信號固有頻率

表2 固有頻率誤差

表3 仿真信號阻尼比

表4 阻尼比誤差

3.3 VMD參數的影響

EMD是一種自適應的分解算法，但是過程總可能會產生虛假分量[6]。VMD分解依賴參數的選擇，如圖6所示，當罰參量α=128，隨著層數參量增加，出現虛假的極點列，虛假的極點列中不會出現穩定極點，對于真實的模態結果，不會受到參數的影響。

圖6 α=64時TEO算法模態參數極點穩態圖Fig.6 Modal parameter identification for TEO of stabilization diagram at α=64

選擇不同的罰參量α，觀察對結構系統參數識別的結果，如表5所示。5階固有頻率和阻尼比均達到穩定狀態需要的層數，隨著罰參量α的增加而增加。

通過觀察仿真信號，不難發現，罰參量α不是最優參數時，會引入虛假分量，但是虛假分量的阻尼比不會穩定，層數參量K偏小時，無法得到全部模態頻率，因此只需最大分解層數Kmax，真實的極點會形成穩定的極點列。由此可見，應用VMD的層數參量可以有效地幫助工程人員找到真實的模態固有頻率與阻尼比，排除虛假分量。

表5 參量α與穩態層數關系

Tab.5 Relationship between parameterαand steady-state layer number

αα 646512 812851 0241025662 04810

3.4 測量噪聲的影響

模態實驗中通常信噪比在30～40 dB時可以確保結果的準確[17]，因此為了驗證所述方法的魯棒性，對仿真信號添加噪聲，對信噪比10～40 dB的測試信號進行結構模態參數的識別。

圖7 固有頻率的誤差趨勢Fig.7 Error of natural frequency

圖8 阻尼比的誤差趨勢Fig.8 Error of damping ratio

如圖7所示，估計得到固有頻率除第1階以外，固有頻率誤差趨勢與信噪比關系不明顯。圖8所示阻尼與信噪比的誤差趨勢，同樣除第1階以外，與信噪比關系不明顯。由此得出結論，本研究所提基于VMD與TEO的結構模態辨識方法具有較好的魯棒性，計算效果良好。

4 實測算例

實測算例來自某動車組齒輪箱模態實驗，如圖9所示。動車組運行速度越來越高,其安全可靠性的要求就越來越高。齒輪箱作為動車組車輛傳遞扭矩驅動車輛行進的關鍵零部件,其安全穩定性直接影響動車組車輛的運行可靠性[18]。由于低階模態對箱體的振動特性起主要作用，因此在持續工況時，齒輪箱箱體的固有頻率應避開轉動頻率，避免共振現象。

圖9 實驗系統Fig.9 Experimental system

分別通過錘擊法與環境激勵兩種方法進行實驗，測點布置如圖10所示，采樣頻率為5kHz，改進TEO算法中低通濾波閾值設為2kHz，分析頻率分析范圍為0～2 kHz。從所布置傳感器的對立面的水平徑向齒輪箱進行錘擊實驗，并采集敲擊過程中的振動信號。環境激勵時，選取齒輪箱從最高速減速至靜止。

圖10 齒輪箱測點布置Fig.10 Measurement setup of gearbox

從水平徑向( 圖9所示x方向) 對齒輪箱敲擊并采集其振動加速度信號，時域信號如圖11(a)所示，環境激勵信號如圖11(b)所示。利用VMD 對其進行分解，最大分解層數Kmax=20，罰參量設為α=5 000。應用改進的TEO算法得到每一層BIFM的固有頻率與阻尼比，頻率變化誤差設為1%，阻尼變化誤差設為5%，得固有頻率與阻尼比的穩態結果如圖12所示。

圖11 時域波形Fig.11 Time-domain wave

圖12 齒輪箱穩態圖Fig.12 Stability diagram of gearbox

由圖12可看出，第1條明顯的極點列及30 Hz附近，由于阻尼不滿足閾值條件，沒有穩態極點，因此不是齒輪箱固有模態頻率。根據前述作法，找出極點列穩態極點，如表6所示為穩態極點組成的前7階模態固有頻率與阻尼比。

表6 信號前7階固有頻率與阻尼比

Tab.6 Signal In front of the seven order natural frequency and damping ratio

/Hz/Hz1 61.60.048 671.70.025 1 83.10.039 52172.40.029 2178.40.011 63245.20.024 9223.1 0.013 3 4354.50.014 8342.90.009 75437.10.015 7434.50.011 16547.70.011 7547.30.008 47694.20.009 1696.8 0.006 2

根據設計要求齒輪箱在持續工況時的小齒輪轉速為4 174 r/min,車軸的轉頻為19.37 Hz，小于箱體的1階模態頻率(71.7 Hz)，而該工況下的齒輪對嚙合頻率為1 530.47 Hz，遠高于前7階固有頻率。由于低階模態對箱體的振動特性起主要作用，因此在持續工況時，齒輪箱箱體不會發生共振現象。

如表6所示，錘擊法中60～80 Hz附近極點列隨著K的增加出現頻率分離，而環境激勵法中僅存在一個穩態極點71.7 Hz。因此錘擊法中識別的第1階固有頻率不能確定，而環境激勵法中穩態極點十分確定。除第1階固有頻率不做評判，其他各階識別結果偏差均小于10%，錘擊法和環境激勵法的模態分解結果基本一致。環境激勵的阻尼比識別結果均小于力錘法30%～50%。由此可以得出，本研究所提方法適用于環境激勵下的模態系統辨識。

5 結 論

1) 應用變分模態分解與Teager能量算子相結合分解提出了一種新的結構系統辨識方法：應用VMD得到滿足TEO法的單一頻率BIMF分量；通過VMD層數參量K的變化尋找穩態極點，作為模態固有頻率與阻尼比的識別作法。

2) 通過改變層數參量K來尋找穩定的固有頻率與阻尼比，解決了VMD參數選擇這一方法缺陷。通過極點穩定的作法識別模態參數不僅適用于本研究所提的VMD-TEO識別方法，同時適用與VMD希爾伯特法與VMD中心頻率法。通過仿真信號對比，所提方法將固有頻率識別的平均誤差降了0.14%。在不同白噪聲情況下，本研究所提方法具有較好的魯棒性。

3) 數值和實驗結果表明，所提出的方法能夠很好地利用各個模式的有限帶寬進行信號分解，并準確地提取結構模態系統的固有頻率與阻尼比。