基于VMD 的谷物清選篩振動分析與結構優化

龐靖，林毅，王升升，杜哲，謝樂樂，陳新奇,2

（1. 河南科技大學農業裝備工程學院，洛陽 471003；2. 洛陽智能農業裝備研究院有限公司，洛陽 471003）

0 引言

隨著農業機械化進程的不斷推進，農業裝備的作業質量、可靠性等也需不斷提高。振動是任何結構都具有的固有特性，對機械裝備的振動特性被廣泛關注。聯合收獲機作為田間收獲作業的重要裝備，其復雜的振動特性不僅會影響機器的可靠性和壽命，也會造成難以預料的作業損失[1-4]。

目前國內外對于聯合收獲機的相關研究主要包含整機振動測試[1-3]、應力分析[4]、有限元模態分析與模態試驗[5]、清選篩內部氣流場分析[6]、作物特性研究[7-8]以及關鍵結構的優化設計等[9]。張黎驊等[9]針對大豆玉米套作場景下的收獲機清選作業質量問題，設計了一種應用于大豆玉米兼收的收獲機清選裝置，使用流固耦合和臺架試驗探究了振動頻率、上下篩傾角等對含雜率和損失率的影響，并得到了工作指標最優的影響因子組合；劉鵬等[10]以多參數可調的清選系統為分析對象，對大豆聯合收獲作業中的清選參數進行了研究，得到了魚鱗篩開度、風機轉速和作業速度等因素組合對損失率等指標的影響；李永磊等[11]對所提出的批次式種子清選機橡膠球清篩裝置測定了橡膠球的激振力，探究橡膠球激振力對清選性能的影響；張學軍等[12]針對油葵聯合收獲中籽粒含雜率和損失率偏高的問題，以振動頻率為因素之一，探究其對籽粒損失率和含雜率的影響；LI等[13]針對谷子聯合收獲時的清選效率問題，以探究自主設計的雙風機振動篩谷子清選裝置的工作性能為目標，設計了因素試驗，其中，篩體振頻是其考慮的因素之一。EBRAHIMI等[14]構建了描述聯合收獲機割臺振動的有限元模型，使用頻域分解法對割臺在運行狀態下的振型參數進行了估計，發現在50 Hz 附近存在共振，采用質量變化策略估計割臺振動系統的頻響函數矩陣，最終，通過優化割臺結構，使得割臺的第五階振型50 Hz 的固有頻率降低為48 Hz，有效地降低了割臺的共振情形；YILMAZ等[15]針對聯合收獲機作業中的噪聲振動問題，使用兩個加速度計和兩個麥克風采集傳統聯合收獲機脫粒過程的噪聲和振動，得到了發動機是主要的激振源；GALKIN等[16]針對谷物振動氣動分選機清種效率低、損失率高等問題，推導了分選機在工作時種流的理論運動模型，計算了顆粒物料分層的氣流速度和物料表面的流動速度，研究結果表明，使用改進型振動氣動分選機在進行谷物清選時，可降低30%～50%的功耗，減少1.5 倍的種子損失。

綜上，目前針對聯合收獲機的相關研究已取得了較多的成果，但有關清選系統振動特性研究的相關文獻相對較少。清選是聯合收獲機作業流程中的重要一環，由于加工、裝配誤差，以及清選篩在機器作業時受到的復雜外界激勵，會引起清選篩的不平衡振動，進而加劇清選篩關鍵部件的磨損，影響其清選性能，導致含雜率增高等問題。因此，本文將基于動態信號分析理論中的變分模態分解（VMD，variational mode decomposition）結合峭度及模糊熵等指標對谷物清選系統進行振動測試與分析，圍繞其中因轉子系統不平衡而造成的不平衡振動問題，對采集到的原始信號進行VMD 分解和特征提取，然后使用懲罰函數法完成對清選篩轉子系統的結構優化。以期提高聯合收獲機械作業穩定性并為其減振分析提供參考。

1 谷物清選系統組成及工作原理

谷物清選系統結構簡圖如圖1a 所示，主要包括風機1、篩箱2、篩體3 和篩體驅動裝置4。清選系統模型如圖1b，其中，驅動裝置為曲柄滑塊機構，篩體上篩為魚鱗篩、下篩為編織篩。在風機作用下，篩箱內部形成風場，由于谷物脫出物空氣動力學特性的差異，將其中的輕質雜余吹出篩箱；動力通過鏈傳動傳輸至曲柄，曲柄帶動篩體運動，除了驅動裝置，上篩和下篩分別通過螺栓固定在篩體支架上，隨著篩體運動而振動，籽粒落在篩體網面上受到振動落入網面空隙而達到清選的目的。

圖1 清選系統結構示意圖Fig.1 Schematic diagram of cleaning system structure

在清選系統工作時，受自身結構以及一些外部激勵的作用，造成轉子系統不平衡，使得整個裝置振動不平衡，進而影響清選篩的工作效率與壽命。因此，需要對清選篩的振動特性進行研究。

2 驅動機構轉子系統不平衡分析

2.1 驅動機構轉子系統組成

本研究的驅動機構轉子系統示意如圖2 所示，主要由軸承座、偏心輪、中間軸、篩體連接盤、驅動鏈輪和配重盤等組成。其中，配重盤的質量為3 500 g，質心與驅動鏈輪中心的距離為135 mm，中間軸軸徑為41 mm，偏心輪和篩體連接盤以及中間軸之間通過一對軸承連接；中間軸通過軸承固定在軸承座上，實現動力的傳遞。

圖2 驅動機構轉子系統模型Fig.2 Driving mechanism rotor system model

2.2 轉子系統動力學模型

根據有限元分析理論，轉子系統可描述為由輪盤單元、軸段單元和軸承座單元組成的系統整體，其中，各個單元之間通過節點來聯結。根據轉子動力學理論，轉子系統的運動方程可分解為各個單元的運動方程。拉格朗日方程是建立轉子系統動力學模型主要的廣義模型[17]，其數學表達式為

式中T為系統總動能（J）；V為系統總勢能（J）；qi為廣義坐標；Qi為系統所受廣義力（N）。

假設剛性輪盤單元質量為m1、外徑為D1、內徑為D0、極轉動慣量為J1，輪盤單元的慣性力主向量在空間坐標系中x軸和y軸上投影的坐標為(x1,θy)和(y1,-θx)，則該輪盤單元處節點位移可表示為：u1=[x1,θy]T，u2=[y1,-θx]T。

基于拉格朗日方程可進一步求得輪盤單元的動力學方程為

式中Qwdk(k=1,2)為輪盤單元所受的廣義力(N)。

同理可得軸段單元的動力學方程為

式中Ms為軸段單元的運動慣性矩陣，ΩJs為回轉矩陣，Ks為剛度矩陣，Qks為軸段單元所受的廣義力(N)。

對于軸承座單元的動力學分析，不僅要考慮其固體接觸，也要分析油液接觸，基于Hertz 接觸理論，受到載荷的作用時，假設兩構件接觸時一構件的粗糙表面微型凸起與另一構件剛性結合面接觸時的真實接觸半徑為r1，不考慮變形時兩構件接觸半徑為r2。根據Lankarani-Nikravesh 碰撞接觸模型[18]，綜合考慮兩碰撞體的材料屬性、碰撞體間的相對碰撞速度及碰撞體的局部變形等信息可得固體接觸部分碰撞力模型為 δ(-)為兩部件碰撞相對速度（m/s）。

式中K'為接觸剛度系數，δ為碰撞刺穿深度（mm），n為非線性指數（金屬材料取值1.5），ce為碰撞恢復系數，

此外，軸與軸承配合之間存在潤滑劑，因此，需要分析軸與軸承結合部位的液體接觸。

由于間隙往往非常小，即油膜厚度很小，因此，假設油膜上的承載力方向在油膜厚度方向（y）上保持不變，且不考慮動力粘性系數μ隨空間位置的變化。則基于液體平均流動的廣義雷諾方程，軸與軸承之間油膜流動的動力學方程可表示為[19]

式中?2為拉普拉斯算子，div 為散度，Fi(i=x,y,z)為作用在油膜內部微元體內所有質量上的力。

2.3 驅動機構轉子系統不平衡分析

旋轉件的異常振動會影響整機的正常運轉，轉子不平衡是旋轉機械的主要激振源[17]，不僅會影響轉子系統的穩定性，也會加劇磨損等進而造成轉軸疲勞、聯結失效等。旋轉件在運行中即使是很小的質量偏心也會產生較大的離心力而引發振動，對于平衡轉子的常見方法包括靜平衡和動平衡等。無試重動平衡法[20-21]旨在通過理論分析來代替實際測試進而對不平衡量進行估計或計算，常見的如動力學模型求解、影響系數估計和智能算法求解等。為了分析本文驅動機構轉子系統的不平衡特性，將配重塊視為偏心配重圓盤（如圖3 所示），對偏心配重圓盤建立如下的動力學方程：

圖3 偏心配重盤示意圖Fig.3 Schematic diagram of eccentric counterweight plate

式中J2為偏心配重盤的極轉動慣量，e為配重盤安裝偏心距，m2為配重盤質量，? 為轉子不平衡量相位，Ω=?。

最終，將各個單元的動力學方程進行綜合可得轉子系統響應的運動方程為

其中，M為系統各單元集成后的質量矩陣；C為各單元集成后的阻尼矩陣；K為各單元集成后的剛度矩陣；R為系統的位移向量；Q為各單元集成后的廣義力向量。

針對清選篩的不平衡振動問題，假設其不平衡響應解為

c為阻尼系數；k為剛度系數。

由式（9）可求解處系數A1、A2、B1、B2，進而得到不平衡信號幅值特征矩陣。

3 不平衡振動分析

3.1 變分模態分解

變分模態分解于2014 年由DRAGOMIRETSKIY等[22]提出。該方法是一種非遞歸的變分分解方法，和其他信號分解方法相比，可以指定分解的模態數，在對模態函數求解過程中，采用鏡像延拓有效避免了類經驗模態分解在信號分解中的端點效應[23]，但對模態函數的求解前需要指定模態數、二次罰項等參數。作為一種信號分解方法，其將任何信號都假設為由一系列具有特定中心頻率、有限帶寬的子信號組成，基于經典維納濾波、希爾伯特變換和混頻和外差解調等原理，構建描述約束變分問題的數學表達式[22-23]，并對其進行迭代求解，得到中心頻率和帶寬。

3.2 模糊熵與峭度

模糊熵[24-25](FE,fuzzy entropy)是ZADEH 在1968 年提出的來描述一個模糊集模糊程度的概率學概念，該方法用均值代替絕對幅值差，以指數函數作為模糊函數來衡量兩個向量之間的相似度。在信號分析中以待測信號作為樣本集，計算其模糊熵YFE，YFE越大，兩信號相關性越差；YFE越小，兩信號相關性越好。

峭度(Kurtosis)是一種無量綱指標，該值越大，信號中所包含的沖擊成分和特征信息就越豐富[26]，在正常工作狀態下，軸承座處振動信號接近正態分布，峭度指標接近3[27]。而在一些不平衡沖擊下，會使軸承座處振動信號偏離正態分布，不平衡振動沖擊成分會混疊在原始信號中，峭度值越大，不平衡振動沖擊成分特征更易于提取。因此，通過對峭度的計算可以有效的識別特征信號。其定義如下：

式中μi為本征模式分量ui的均值，σi為ui的方差。

3.3 基于VMD 及綜合指標信號分析流程

盡管基于VMD 的振動信號分析在信號時間序列的預測與成分分解方面是可行的，但僅僅依靠該方法無法準確得到關于所要提取特征本身的信息。在一些研究中，模糊熵和峭度被構建為一種綜合指標應用于軸承的故障診斷中[28-29]，綜合指標K的定義如下：

通過綜合指標K可有效的找到振動信號中包含沖擊成分和特征信息最豐富的分量，進一步的，對提取后分量信號的進行包絡譜分析，得到不平衡振動頻率。該方法的主要步驟如下：

1）采集清選篩在不同工況下的振動數據樣本，對其進行時、頻域分析，選取頻率成分豐富的測點采集到的振動信號進行VMD 分解，得到有限個本征模式分量ui。

2）計算各分量的綜合指標，選擇綜合指標最小的分量作為敏感分量，用于后續的特征提取。

3）對選擇的敏感本征模式分量進行希爾伯特變換求其包絡譜，提取不平衡振動。

4 清選篩振動臺架試驗

4.1 清選篩試驗臺工作參數

本文所用清選篩試驗臺為河南科技大學自主設計的谷子切-軸流脫粒清選系統，如圖4a 所示。其結構主要包含兩級脫粒系統、清選系統和電機等，兩級脫粒系統主要由一級切流滾筒和一級軸流滾筒以及凹板篩組成。如圖4b 所示，在裝置運行時，動力由電機輸出，通過動力中間傳動軸傳輸至風機上，并通過一級帶傳動和一級鏈傳動傳輸至清選篩驅動轉軸上，其次，通過帶傳動將動力輸出到脫粒滾筒上，兩級脫粒滾筒之間通過鏈傳動。

圖4 谷物清選篩試驗臺及其傳動示意圖Fig.4 Grain cleaning sieve test-bed and its transmission schematic diagram

激勵源是引起結構振動的核心因素，除此之外，軸承缺陷、結構件的制造或裝配誤差以及旋轉件的配重設計等也是重要因素。其中，外部激勵源包括驅動裝置或作業環境中可能存在的激勵源，如地面起伏等。本文旨在分析轉子系統不平衡引起的振動問題，忽略軸承缺陷、系統的制造和裝配誤差，且不考慮地面的激勵作用，對清選篩試驗臺進行了傳動計算和工作參數分析，其中，電機、脫粒滾筒、中間傳動軸和風機為主要激勵源，分別計算了其理論工作轉速和激勵頻率，如表1 所示。

表1 試驗臺工作參數Table 1 Working parameters of the test bed

4.2 振動試驗測試系統及方法

如圖5 所示，振動試驗測試系統由加速度傳感器、DH5902 動態信號采集儀、計算機、變頻器、異步電機和試驗臺架組成。其中，三向加速度傳感器量程為±50 g，頻響范圍為0.3～6 kHz，x向、y向和z向靈敏度分別為101.1、97.3 和99.3 mV/g；單向加速度傳感器靈敏度為96.1 mV/g；數據采集儀共有16個通道，最高采樣頻率為100 kHz。

圖5 測試系統Fig.5 Testing system

為了準確反映清選篩的振動特性，對試驗臺進行了預試驗，基于試驗臺的動力傳動特點與主要的激勵源，選擇篩體支撐（篩體驅動裝置支撐）軸承座為測點5、篩體尾部側板為測點6，其次，為了反映各個激勵源的實際振動特點，選擇動力中間傳動軸軸承座為測點1、第一級脫粒滾筒轉軸軸承座為測點2、風機轉軸軸承座為測點3、第二級脫粒滾筒轉軸軸承座為測點4，共6個測點（圖6）；其中，測點1 和3 分別布置一個單向加速度傳感器，其余測點布置為三向加速度傳感器，共使用DH5902 采集系統的14個通道；傳感器的安裝選擇膠粘和磁吸兩種方式，在軸承基座上（測點1、測點3）由于基座表面不平整，選擇粘接對傳感器進行固定，其余測點均選擇磁吸方式固定。

圖6 試驗測點布置Fig.6 Layout of test measuring points

此外，為了便于提取和分析各個加速度傳感器檢測的振動信號，從而更加準確的得到每個測點的振動情況，構建如圖6 所示的試驗臺全局坐標系和各個測點的局部坐標系。以物料喂入方向為全局坐標系x方向，垂直大地方向為全局坐標系z向。測點1、3 分別檢測中間傳動軸和風機軸承座的徑向振動，規定其局部坐標系如圖6，測點2 與測點4-6 分別檢測一、二級脫粒滾筒軸承座、篩體支撐軸承座和篩體尾部側板在全局坐標系3個方向上的振動。

4.3 信號采集與時、頻域分析

在完成測試系統的調試和預試驗后，考慮到試驗臺的實際應用場景，即履帶式聯合收獲機的作業場景，基于試驗臺的工作參數，以設計的動力輸入軸額定轉速900 r/min 為基礎，并根據采樣頻率和采樣定理，選擇5 Hz為激勵源的頻率分辨率，分別在電機轉速為300、600、900 r/min 三種工況下對6個測點進行測試。限于篇幅，選取清選篩上的測點5 和測點6 在局部坐標系3個方向上的時域信號及其均方根值如圖7、圖8 所示。

圖7 測點5 在3個方向上的時域信號Fig.7 Time domain signals of measuring point 5 in three directions

圖8 測點6 在3個方向上的時域信號Fig.8 Time domain signals of measuring point 6 in three directions

通過計算振動信號的加速度均方根值（root mean square，RMS）可見，隨著轉速的提高，各個測點在不同方向上的信號的RMS 值越來越大，其中，測點5、6在3個方向上信號的最大RMS 值分別為0.004、5.982、3.617、6.581、24.327 和25.123 m/s2，即振動強度越來越大；其次，三種工況下各個通道采集的信號趨勢相當，隨著整個裝置輸入轉速的增大，脫粒滾筒1 和2（測點2、測點4）、清選裝置尾部側板（測點6）與篩體支撐軸承座（測點5）以及風機轉軸（測點3）的振幅均呈增大趨勢且信號中帶有大量的噪聲和不平衡成分，但中間傳動軸（測點1）和篩體支撐軸承座局部坐標系x方向（測點5-x向）的振幅在0 附近波動且未有較大改變，究其原因，主要是因為中間傳動軸處與軸承座的裝配精度高，從而使得該部位軸在運轉過程中的徑向振動小。而篩體支撐軸承座處局部坐標系x方向的振幅趨近于0 說明清選篩在全局坐標系x方向的振動小，強度趨近于0，這也為下文進一步的分析提供了參考。

為了初步確定測點5、6 處的振動情況，對工況3 下測點5 在其局部坐標系y和z方向與測點6 在其局部坐標系各個方向上的振動信號進行功率譜密度計算如圖9所示。

圖9 轉速為900 r/min 時測點5、6 不同方向功率譜密度Fig.9 Power spectral density in different directions of measurement points 5 and 6 at a rotational speed 900 r/min

其中，由于測點5-x向處信號的振動強度很小，所以不考慮在該處局部坐標系x方向上的振動。除此之外，測點5 處采集到的信號中所包含的特征頻率主要分布在12、13、18、19、20、24、25、30 和36 Hz 附近；測點6 處采集到的信號中所包含的特征頻率主要分布在6、12、13、18、19、20、24、25、26、30、31 和36 Hz 附近。

其中，在12、13、18、20、24、25、30 和36 Hz 附近，兩個測點處均存在對應的振動。

4.4 基于VMD 及綜合指標的不平衡振動提取

對于實測振動信號，其中往往夾雜著噪聲和特征信號等復雜的成分，僅僅通過對實測振動信號的時、頻域分析，不足以明確其準確的振動特點，需要進一步分析。

由于本文主要探討轉子質量偏心所造成的的不平衡振動問題，而測點6 主要測的是篩體尾部側板的振動，考慮其中可能存在因零件的加工、裝配誤差而造成的振動，在轉速為900 r/min 下，選擇篩體驅動裝置支撐軸承座處測點5 局部坐標系y和z方向的振動信號使用變分模態分解（VMD）法對其進行分解（如圖10 所示）。有關模態個數k的選擇，基于對整個裝置運動特點和理論激勵源的分析，初步選擇k為5，7，9 進行預分解，考慮到模態混迭，最終確定k取7。之后分別對測點5局部坐標系y、z方向的振動信號進行VMD 分解。得到測點5 在局部坐標系y和z方向上信號的VMD 分量后，分別計算各分量的模糊熵和峭度，并得到綜合指標K，如表2 所示。

表2 VMD 分量的特征量計算結果Table 2 Calculation results of characteristic quantities of VMD components

圖10 測點5在y和z 向的VMD 結果Fig.10 VMD results of measuring point 5 in y and z directions

由表2 可以看出，在測點5 局部坐標系y和z兩個方向振動信號VMD 后的7 階本征模式分量(u1-u7)中，u1的綜合指標K最小，可初步確定u1與原始信號之間的敏感性最優。

5 清選篩結構優化

5.1 驅動機構轉子系統不平衡模型計算

為了減小因驅動機構因轉子系統不平衡而造成的清選篩不平衡振動，應盡可能的使不平衡量最小，即不平衡振動信號幅值盡可能小。基于清選篩曲柄滑塊機構的設計參數，選擇配重盤質量m2和安裝偏心距e的優化區間分別為[500g,5 000g]和[41mm,200mm],利用約束優化方法建立了驅動機構轉子系統的不平衡優化數學模型如下：

其中：ψ=[m2,e]為優化變量矩陣。

然后通過懲罰函數法[29]利用SUMT 程序在Matlab中完成優化參數的計算。

其運算過程為：首先隨機選取優化變量初始值ψ0、懲罰因子初值r0（根據經驗公式試算選擇）和縮減系數c（經驗取值在0.1～0.7 之間），然后根據收斂條件依次迭代k次，滿足收斂條件后迭代結束，輸出約束最優解ψ*。將驅動機構轉子系統的不平衡優化數學模型導入SUMT 程序中，通過Matlab 進行計算后得到：m2=3 650 g，e=136.7 mm。

5.2 優化效果驗證

在轉速為900 r/min，優化前后測點5 局部坐標系y和z方向上采集的振動信號進行VMD 分解并對其敏感分量進行包絡譜分析如圖11 所示。

圖11 優化前后測點5在y和z 兩個方向敏感分量包絡譜分析結果Fig.11 The envelope spectrum analysis results of the sensitive component of the measuring point 5 in y and z directions beforeand after optimization

從圖11 可以看出，在30.13 Hz 處振動強度最大，分別為0.29 和0.24 dB。此外，基于4.1 節有關激勵源的分析，剔除激勵源附近的頻率點，選擇30.13 Hz 為不平衡振動的主要特征頻率。

相比優化前，優化后在30.13 Hz 的不平衡沖擊振動強度分別為0.12 dB 和0.11 dB，相較于優化前的0.29和0.24 dB，振動能量最大下降了58.62%，證明了對驅動機構轉子系統優化的有效性。

6 結論

本文針對谷物清選篩的振動問題，基于VMD 及綜合指標對清選篩實測振動信號進行分解、特征提取和驅動機構的優化。基于動態信號峭度和模糊熵所構建的綜合指標對VMD 后的模式分量進行計算、選取敏感分量，最后分析敏感分量的包絡譜確定不平衡振動的特征頻率，主要結論如下：

1）基于VMD 信號分解的特征提取方法及由峭度和模糊熵構建的綜合指標，通過選擇最優模式分量并計算其包絡譜，能夠有效的提取信號中的特征成分，在谷物清選篩不平衡振動分析與減振分析中具有良好的實用性。此外，也可為其他結構的振動分析提供參考。

2）通過計算清選篩篩體支撐軸承座處信號VMD 后本征模式分量的綜合指標，得到第1 階本征模式分量的綜合指標最小，分別為0.706 和0.241，選擇第1 階本征模式分量為原始信號敏感分量并計算其包絡譜，得到在電機轉速(即動力輸入軸轉速)為900 r/min 時，清選篩的不平衡振動主要在30.13 Hz 處。

3）采用懲罰函數法建立了清選篩驅動機構轉子系統的結構優化函數并進行計算。優化前30.13 Hz 處的振動強度在0.24～0.29 dB；優化后轉子系統配重塊質量為3 650 g，距輪盤中心的偏心距為136.7 mm，30.13 Hz處的振動強度在0.11～0.12 dB。