周鋮，羅楊，魏江，曹宏瑞，2，蘭海，張萬昊

(1.西安交通大學 機械工程學院，陜西 西安 710049;2.西安交通大學 機械制造系統(tǒng)工程國家重點實驗室，陜西 西安 710049;3.中國北方車輛研究所，北京 100072)

0 引言

履帶車輛因其良好的機動性能在現(xiàn)代軍事、農(nóng)業(yè)及建筑業(yè)都發(fā)揮著關(guān)鍵性作用［1-2］。隨著我國綜合國力的提升，科學技術(shù)以及軍事力量也有極大的發(fā)展，為保證履帶車輛在運輸和戰(zhàn)斗過程中具有更高的機動性能，履帶車輛發(fā)動機也向高速大功率方向發(fā)展［3］，對于工況復(fù)雜的履帶車輛而言，其制動器的扭振現(xiàn)象愈加顯著。

發(fā)動機是履帶車輛的核心，除提供動力外，也是扭振的激勵源［4-5］。當發(fā)動機產(chǎn)生扭振后經(jīng)過一系列的傳動裝置最終將激勵傳遞到制動器，從而使得制動器發(fā)生扭振。制動器發(fā)生扭振對其可靠性是一個潛在的極大威脅，在很大程度上扭振會造成制動器制動效果降低甚至失效，從而引發(fā)事故。因此對制動器扭振的研究在履帶車輛的研究中尤為重要，對其參數(shù)的測試不僅可以為制動器的內(nèi)部尺寸設(shè)計提供參數(shù)，降低制動器的扭振幅值，并且能為履帶車輛制動器中制動片的更換提供參考，保證履帶車輛制動器在工作過程中的可靠性。

目前，測量扭振的方法還在不斷發(fā)展的過程中，對于扭振信號的獲取精度、有效的數(shù)據(jù)分析方法和扭振監(jiān)測的實時性都需進一步得到提高［6-7］，尤其針對處于復(fù)雜惡劣、空間狹小結(jié)構(gòu)緊湊的環(huán)境中的旋轉(zhuǎn)部件。本文就某型履帶車輛制動器的測試環(huán)境惡劣、工況復(fù)雜、干擾強等特點，展開了對履帶車輛制動器扭振瞬時頻率特征提取的研究，采用脈沖時序計數(shù)法［8-10］的扭振信號獲取方法，研究基于參數(shù)優(yōu)化變分模態(tài)分解(VMD) 的扭振信號瞬時頻率特征提取方法，通過采用優(yōu)化的VMD 算法實現(xiàn)對扭振信號的分解降噪以及重構(gòu)，采用零點插值計算扭振信號兩脈沖之間的時間間隔，完成對扭振信號瞬時轉(zhuǎn)速波動信號的計算，而后采用頻譜分析，提取扭振信號的瞬時頻率特征。

1 扭振信號瞬時轉(zhuǎn)速獲取方法

1.1 扭振信號的測量原理

脈沖時序計數(shù)法扭振測量技術(shù)是利用安裝在軸上的等分齒輪結(jié)構(gòu)，通過非接觸式磁電、光等傳感器產(chǎn)生相應(yīng)的脈沖信號［11］。碼盤或等分齒輪每轉(zhuǎn)過一定角度就會產(chǎn)生一個脈沖，測得兩個脈沖間隔時間，就可獲得兩脈沖之間的平均角速度，進而獲得扭轉(zhuǎn)角和角速度差。當轉(zhuǎn)軸不發(fā)生扭振時，轉(zhuǎn)軸的瞬時速度等于其平均速度，傳感器經(jīng)過等分齒輪時輸出的脈沖為等間隔的［12］;當轉(zhuǎn)軸發(fā)生扭振時，相當于在轉(zhuǎn)軸平均速度上疊加了轉(zhuǎn)速波動，則傳感器輸出的脈沖為不等間隔的。

假設(shè)等分齒輪齒數(shù)為Z 個，則每個齒的分度Δθ 為

或用弧度制(rad) 表示Δθ 為

根據(jù)傳感器輸出波形可采用高頻A/D 采樣法［13］或高頻計數(shù)法［14］求出轉(zhuǎn)過第i 齒時的脈沖信號的周期Ti(s)，則轉(zhuǎn)過第i 齒(i=1，2，3，…，Z) 時的瞬時轉(zhuǎn)速ni為

由此可得轉(zhuǎn)軸轉(zhuǎn)過一周的平均轉(zhuǎn)速為

則等分齒輪轉(zhuǎn)過第i 齒的扭振瞬時轉(zhuǎn)速波動為

1.2 扭振信號瞬時轉(zhuǎn)速的提取

由1.1 節(jié)計算公式可知，扭振測量的關(guān)鍵是獲取扭振信號周期和計算瞬時轉(zhuǎn)速，瞬時轉(zhuǎn)速的提取主要采用零點插值法［15-16］。

零點插值法是指通過線性插值的方式找到原始信號波形與零線的交點序列{p1，p2，…，pZ}，如圖1所示。

圖1 零點法求瞬時轉(zhuǎn)速示意圖Fig.1 Diagram of instantaneous speed by the zero point method

分析圖2 所示的信號片段，可知單個齒經(jīng)過所產(chǎn)生的信號為pi-1至pi+1段。

圖2 線性插值零點法計算瞬時轉(zhuǎn)速Fig.2 Calculation of instantaneous speed by the zero point linear interpolation method

轉(zhuǎn)過第i 齒經(jīng)歷的時間為

式中:T 為采樣卡采樣周期。

第i 齒的瞬時轉(zhuǎn)速為

2 參數(shù)優(yōu)化VMD

2.1 VMD 算法基本原理

VMD 是一種非遞歸、自適應(yīng)分解算法，每個本征模式分量(IMF) 分量的帶寬和中心頻率利用迭代確定［17］，由K 個IMF 組成的模態(tài)函數(shù)uk(t) (k=1，2，…，K) 的具體步驟如下:

1) 通過Hilbert 變換，得到各模態(tài)函數(shù)uk(t) 的解析信號，獲得uk(t)的單邊頻譜uk(t)，δ(t) 為Dirichlet 函數(shù)，* 為卷積符號。

2) 進行頻率混合，對各模態(tài)函數(shù)uk(t) 的解析信號混合一個預(yù)先估計中心頻率{ω1，ω2，…，ωK}表示各IMF 分量的頻率中心。將各模態(tài)的頻譜變換到基頻帶

3) 通過高斯平滑解調(diào)信號獲得每段帶寬，即L2范數(shù)梯度的平方根。構(gòu)造出約束的變分模型為

式中:{uk}={ u1，u2，…，uK} 表示分解得到的K 個IMF 分量;?t為對函數(shù)求時間的偏導(dǎo)數(shù);x(t) 為分解的主信號。

為求解上述約束的變分模型的最優(yōu)解，引入懲罰因子α 構(gòu)造增廣Lagrange 函數(shù)為

式中:λ 為Lagrange 乘子;α 為懲罰參數(shù);〈·〉表示內(nèi)積運算;λ(t) 用來保持約束條件的嚴格性。

將上述Lageange 函數(shù)從時域轉(zhuǎn)換到頻域，并進行極值求解，可獲得迭代n +1 次過程中模態(tài)分量在頻域中的表達式，即

4) 使用交替方向乘子算法(ADMM)［18］進行約束變分模型的最優(yōu)解的求取，將原始信號將被分解為K 個窄帶IMF 分量。

具體算法如下:

3) 根據(jù)如式(13) 更新λn+1:

式中:τ 為正則系數(shù)。

2.2 參數(shù)優(yōu)化VMD

由PSO 算法［19-21］的原理可知，PSO 能夠根據(jù)適應(yīng)度函數(shù)的最值作為優(yōu)化準則來尋找VMD 算法的最優(yōu)參數(shù)組合。利用PSO 搜尋VMD 算法的影響參數(shù)時，需確定一個適應(yīng)度函數(shù)，粒子每次更新位置時計算一次適應(yīng)度值，通過對比新粒子適應(yīng)度值進行更新。本文基于扭振信號的特點采用能量熵［22］作為PSO 的適應(yīng)度函數(shù)搜尋VMD 最優(yōu)參數(shù)組合［Kopt，αopt］。

能量熵Lp為

式中:Lp(k) 表示第k 個IMF 分量的能量熵;pk表示第k 個單分量IMF 能量的歸一化形式;E(k) 表示第k 個單分量IMF 的能量。

基于PSO 優(yōu)化VMD 參數(shù)組合的主要操作步驟如下:

1) 定義粒子的位置和速度范圍。設(shè)粒子群個體的位置對應(yīng)VMD 算法參數(shù)［K，α］的一個候選解;初始化所有粒子位置和速度，并在設(shè)定范圍內(nèi)選取隨機值，通過計算確定目標函數(shù)。

2) 確定PSO 算法的適應(yīng)度函數(shù)，在不同位置對扭振信號進行VMD 運算，并采用式(14) 和式(15) 計算各IMF 的能量熵;

3) 根據(jù)式(16) 計算最小能量熵為局部極小值min Lp(k)，即VMD 的最優(yōu)參數(shù)組合［Kopt，αopt］的一個候選解;

4) 比較不同位置的適應(yīng)度函數(shù)值局部極小值min Lp(k)，并對個體局部極值和群體全局極值進行更新;

5) 對粒子的速度和位置進行更新。轉(zhuǎn)至步驟2，直至達到預(yù)設(shè)迭代次數(shù)，輸出最佳適應(yīng)度函數(shù)值min Lp和粒子位置［Kopt，αopt］。

3 扭振信號處理流程

由扭振信號的測量方法可知，扭振的獲取是一種間接測量法，獲取的扭振信號中包含大量的噪聲，故而必須采用一定的手段或方法進行降噪。本文采用VMD 算法對原始扭振信號進行分解降噪，并利用PSO 算法搜尋VMD 算法最優(yōu)參數(shù)組合，由此提出基于PSO-VMD 的扭振信號處理方法，具體實現(xiàn)步驟如下:

1) 初始化PSO 算法的相關(guān)參數(shù)。PSO 算法設(shè)定的各項參數(shù)如表1 所示。

表1 PSO 各項參數(shù)Table 1 PSO parameters

種群數(shù)量m 反映了計算時需要的粒子的數(shù)目，m 過小，PSO 計算速度快但易于過早收斂，過大，尋優(yōu)效率降低，通常設(shè)置為［10 -100］，Nmax為算法最大迭代次數(shù)。

種群維度D 一般用于反映優(yōu)化的參數(shù)個數(shù)，本文主要優(yōu)化VMD 算法的分解個數(shù)K 和懲罰因子α兩個參數(shù)，故而設(shè)置為2。

學習因子c1和c2反映了粒子間的學習能力，c1表示自我經(jīng)驗、c2表示種群經(jīng)驗。一般c1和c2取固定常數(shù)，并且取c1和c2相等，早期的研究中通常將取值規(guī)定為2［19-21］。

慣性權(quán)重可調(diào)節(jié)PSO 算法局部和全局搜索能力，Shi 發(fā)現(xiàn)，當w 在0.9～0.4 線性變化時PSO 算法優(yōu)化能力更好，即采用線性遞減策略:

2) PSO 搜尋VMD 算法的參數(shù)組合。以VMD算法中的分解個數(shù)K 和懲罰因子α 目標函數(shù)，以能量熵為適應(yīng)度函數(shù)進行優(yōu)化，得到最優(yōu)參數(shù)組合［Kopt，αopt］。

3) 將最優(yōu)參數(shù)組合［Kopt，αopt］重新代入到VMD 算法中完成扭振信號的分解，得到Kopt個IMF分量;

4) 計算各IMF 與原始扭振信號的相關(guān)系數(shù)rK，并取相關(guān)系數(shù)rK最大(或較大) 對應(yīng)的IMF 進行扭振信號的重構(gòu);

相關(guān)系數(shù)用來描述經(jīng)VMD 后得到的各IMF 與原始扭振信號的關(guān)聯(lián)程度，通過計算各IMF 與原始信號自相關(guān)系數(shù)將其歸一化處理［23］，其表達式為

式中: N0為扭振信號y(t) 對應(yīng)的原始信號個數(shù);為第k 個IMF 的自相關(guān)系數(shù);Ry(k) 為扭振信號y(t) 的自相關(guān)系數(shù)。rK越大表示第K 個IMF與原始信號的相關(guān)性越大，含有的轉(zhuǎn)頻信息也越豐富。

5) 將重構(gòu)后的扭振信號進行零點插值，提取扭振瞬時轉(zhuǎn)速波動，最后對扭振瞬時轉(zhuǎn)速波動進行頻譜分析，得到扭振瞬時頻譜特征。

以上扭振信號處理流程如圖3 所示。

圖3 算法流程圖Fig.3 Flow chart of the algorithm

4 實驗分析

為了驗證所提方法的有效性，搭建如圖4 所示的扭振頻率可調(diào)的扭振實驗樣機進行實驗?，F(xiàn)場測試中采用加工有100 齒的等分齒輪，安裝在實驗樣機的合適位置(見圖4(b))，美國Lion 公司的CPL290 電容式位移傳感器采集扭振數(shù)據(jù)。通過億恒8 通道ECON MI-7008 數(shù)據(jù)采集儀記錄，采樣頻率為20 480 Hz，采樣時間為2 s，設(shè)定實驗臺轉(zhuǎn)速為600 r/min、扭振頻率為50 Hz。采集的扭振信號時域波形如圖5 所示。

圖4 扭振實驗樣機Fig.4 Torsional vibration experimental prototype

圖5 扭振信號時域圖Fig.5 Time domain diagram of torsional vibration signals

如圖5 所示為扭振試驗樣機轉(zhuǎn)速為600 r/min、扭振頻率為50 Hz 的扭振信號時域波形圖以及扭振信號的局部放大圖。從時域圖中可明顯觀察到采集到的等分齒輪上信號具有很強的周期性，經(jīng)過計算每一個周期內(nèi)共有100 個脈沖信號(類正弦信號)，對應(yīng)等分齒輪的齒數(shù)，說明采集到的扭振信號是正確的。

圖6 為扭振信號的頻譜分析，從頻譜的低頻部分可得到扭振轉(zhuǎn)速600 r/min 的轉(zhuǎn)頻10 Hz，從高頻部分可得到等分齒輪齒的通過頻率1 000 Hz，且存在10 Hz 的邊頻帶，從而說明采集的扭振數(shù)據(jù)的正確性。但是從頻譜圖中無法得知50 Hz 的來源，說明采用簡單的頻譜分析并不能有效提取出扭振頻率。

圖6 扭振信號頻譜分析圖Fig.6 Spectrum analysis of torsional vibration signals

采用PSO 優(yōu)化VMD 參數(shù)，得到如圖7 所示的適應(yīng)度函數(shù)變化曲線，從圖7 中可得，PSO 迭代1 000 次，粒子群適應(yīng)度函數(shù)值最后達到了收斂，完成收斂時粒子群的迭代次數(shù)為331，對應(yīng)VMD 參數(shù)組合為［K，α］=［5，300］，將優(yōu)化得到的最優(yōu)參數(shù)代入VMD 中對扭振信號進行分解，得到5 個單分量IMF，如圖8 所示。

圖7 適應(yīng)度函數(shù)曲線Fig.7 Curve of fitness function

從圖8 中可以發(fā)現(xiàn)經(jīng)過參數(shù)優(yōu)化后的VMD 能夠?qū)⑴ふ裥盘柾耆蛛x。

圖8 VMD 各IMFFig.8 Decomposition of each IMF by VMD

計算VMD 的5 個IMF 與扭振信號的相關(guān)系數(shù)rK，結(jié)果如表2 所示，選擇相關(guān)系數(shù)最大的IMF2重構(gòu)扭振信號，IMF2時域圖和頻譜分析如圖9 所示。從圖9 中可以得到等分齒輪的工頻以及以轉(zhuǎn)軸轉(zhuǎn)頻10 Hz 的邊頻帶。對重構(gòu)的位移信號采用式(6) 的零點插值法求取單齒經(jīng)過的時間間隔，并用式(7)求解瞬時轉(zhuǎn)速，進而利用式(5) 求解得到如圖10 所示的瞬時轉(zhuǎn)速波動數(shù)據(jù)。對瞬時轉(zhuǎn)速波動信號做頻譜分析得到如圖11 所示的轉(zhuǎn)速波動頻率即為扭振頻率。

圖9 IMF2時域波形及頻譜Fig.9 IMF2 Time domain waveform and spectrum

圖10 瞬時轉(zhuǎn)速波動Fig.10 Instantaneous speed fluctuation

表2 各IMF 分量相關(guān)系數(shù)Table 2 Correlation coefficients of IMF components

從圖11 頻譜分析局部放大圖中可以清楚地看到50.02 Hz 的扭振頻率及其倍頻，同時在扭振倍頻兩側(cè)存在10 Hz 的邊頻帶，從而證明所提算法能夠?qū)崿F(xiàn)對扭振信號頻率的提取，且在考慮頻譜分辨率的情況下可以認為扭振提取誤差應(yīng)小于1%。

圖11 瞬時轉(zhuǎn)速波動的頻譜分析Fig.11 Spectrum analysis of instantaneous speed fluctuation

為了進一步驗證本算法的優(yōu)越性，現(xiàn)對未經(jīng)優(yōu)化的扭振信號和經(jīng)集合經(jīng)驗?zāi)B(tài)分解(EEMD) 分解［24］的扭振信號進行分析，分析結(jié)果分別如圖12和圖13 所示。

從圖12(a) 可看出扭振瞬時轉(zhuǎn)速波動信號毫無規(guī)律，從圖12(b) 扭振瞬時轉(zhuǎn)速波動信號的頻譜局部放大可看出，在頻譜圖中并未找到轉(zhuǎn)軸轉(zhuǎn)頻10 Hz分量，雖然找到了50.02 Hz 成分，但是在其二倍頻處并未找到以10 Hz 為間隔的邊頻帶，并不能說明此50.05 Hz 是否就是所要尋找的扭振瞬時頻率特征。

圖12 未經(jīng)優(yōu)化的扭振信號分析結(jié)果Fig.12 Results of unoptimized torsional vibration signals

采用EEMD 對扭振信號進行分解，經(jīng)過分解發(fā)現(xiàn)，前5 個IMF 基本能反映扭振信號的模態(tài)，而后求解前5 個IMF 與原始信號的相關(guān)系數(shù)，選擇相關(guān)系數(shù)最大的IMF2重構(gòu)扭振信號。

從圖13(d) 扭振瞬時轉(zhuǎn)速波動信號的頻譜局部放大可看出，在頻譜圖中雖然找到了50.024 4 Hz 的扭振頻率，但是并未找到轉(zhuǎn)軸轉(zhuǎn)頻10 Hz 分量，同時也無法找到扭振頻率的倍頻及以10 Hz 為邊頻帶的成分，從而證明了所提的算法的分解效果更佳。同時又分別對經(jīng)VMD 分解得到的其他IMF 進行了同樣的分析，經(jīng)過分析發(fā)現(xiàn)均無法得到扭振頻率特征，充分證明了選擇IMF2進行分析的合理性，同時進一步驗證了所提算法的有效性。而對于各瞬時轉(zhuǎn)速波動圖中所產(chǎn)生的突變成分可能是由于被測轉(zhuǎn)軸質(zhì)量的不平衡或傳感器不嚴格對中所導(dǎo)致，將在后續(xù)的工作中作進一步分析。

圖13 經(jīng)EEMD 的扭振信號分析結(jié)果Fig.13 Analysis results of torsional vibration signal decomposed by EEMD

5 結(jié)論

本文提出了一種使用粒子群參數(shù)優(yōu)化的VMD和零點線性插值法的扭振特征頻率提取方法。得出主要結(jié)論如下:

1) 所提的方法關(guān)鍵在于設(shè)置信號能量熵為粒子群參數(shù)優(yōu)化的適應(yīng)度函數(shù)，為VMD 算法提供分解個數(shù) 和懲罰因子的最優(yōu)參數(shù)組合。在最優(yōu)參數(shù)組合下，VMD 算法能夠準確提取原信號中轉(zhuǎn)頻信息最豐富的分量，因而使所提方法具有獨特優(yōu)越性。

2) 對最優(yōu)分量信號使用零點線性插值法能夠更為準確地還原軸系瞬時轉(zhuǎn)速波動信號，即扭振信號，其頻譜能夠準確體現(xiàn)扭振特征頻率。

3) 慣性載荷扭振實驗結(jié)果表明: 在轉(zhuǎn)速為600 r/min、扭振激勵頻率為50 Hz 的條件下，對采樣頻率為20 480 Hz 的位移信號使用所提方法能精確提取扭振瞬時頻率特征，誤差小于1%。