劉巧斌,史文庫,高承明,陳志勇,閔海濤,陳龍
(吉林大學汽車仿真與控制國家重點實驗室,130022,長春)
橡膠材料是一種優良的隔振材料,而遲滯非線性建模是對橡膠材料的彈性和阻尼等力學特性進行準確而全面表述的有效數學方法。研究橡膠遲滯非線性,對于汽車和機械裝置上橡膠件的力學性能,以及減振、抗沖擊和降噪等直接影響人的感官評價的性能,乃至可靠耐久性能的預測、優化和控制都具有重要的意義。橡膠懸置、橡膠襯套和推力桿橡膠球鉸是汽車上最常用的3種橡膠零件。推力桿作為重型商用車懸架系統的關鍵導向元件,一般通過襯套式橡膠球鉸與車架和車橋相連接,對橡膠球鉸的結構、力學性能、疲勞進行深入分析,并通過優化設計提升性能,可為提高車輛的舒適性、平順性和耐久性提供科學的理論指導,具有重要的工程實用價值。
由于橡膠材料的使用,導致了推力桿縱向力-位移曲線具備了黏彈性非線性材料的遲滯特性。國內外已有大量關于遲滯非線性的研究報道,但對于不同的結構和材料,并不存在一種普遍適用的遲滯非線性模型,需要根據具體情況選取合適的模型。常用的遲滯非線性模型有雙線性模型、Bouc-Wen模型、Davidenkov模型、多項式模型等[1],而對于本文所研究的推力桿,其遲滯特性具有非對稱特性,直接采用以上方法無法對非對稱遲滯特性進行建模,故在對現有對稱遲滯模型進行相應改進的基礎上,提出可以描述推力桿非對稱遲滯特性的模型是本文的一個研究重點。
為獲取模型的參數,有必要采用參數識別方法從試驗測量的數據中提取有效信息,并在保證模型精度的前提下,最大程度地提高建模效率。由于測試環境的影響、傳感器的限制和數據采集記錄過程中出現的種種誤差,使試驗采集到的信號不可避免地包含了噪聲信號的干擾,因此對采集到的信號進行有效地濾波是進一步建模和參數識別的基礎。常用的時域數據平滑濾波方法有卡爾曼濾波、均值濾波、中值濾波、小波濾波、Savitzky-Golay濾波等[2],而Savitzky-Golay平滑濾波方法具有占用內存小、計算效率高和原理簡單等優點,故廣泛應用于數據平滑濾波處理領域[3-5]。各種濾波方法在使用時若參數選取不同,則濾波效果存在較大的差異,如濾波不足就不能有效去除信號中的噪聲,但濾波過度又會引起信號的失真。因此,提出一個實用的濾波指標,并采用優化方法獲取最佳指標下的最優濾波器參數,以提高濾波效果,是本文的另一個研究重點。
當前行業內對于推力桿及橡膠球鉸的研究,大多集中于結構有限元分析、剛度模態預測、動靜態試驗、疲勞分析和優化等方面[6-8],而對其非對稱遲滯特性和試驗數據濾波處理方面的研究,鮮有文獻報道。建立準確高效的推力桿遲滯特性模型是進一步對推力桿進行性能預測、結構優化和整車多體動力學高精度仿真的基礎。
本文在進行推力桿縱向遲滯特性試驗的基礎上,提出了布谷鳥算法優化的灰色Savitzky-Golay濾波方法對試驗數據進行去噪,以準確預測去噪后的遲滯特性為目標,分別對線性阻尼模型、分段線性阻尼模型、非對稱Bouc-Wen模型和高階多項式模型的參數進行了識別,通過與試驗數據的對比,比較了不同模型的精度和效率。
橡膠球鉸在外力的作用下產生非線性恢復力,在加載和卸載過程中,恢復力軌跡不重合,卸載過程與加載過程相比,恢復力存在一定的損失。從能量的觀點看,加載過程存儲的能量,在卸載過程中有一部分通過阻尼產生熱量而耗散。加載和卸載過程中的恢復力-位移曲線,通常稱為遲滯回線或滯回環,非線性的恢復力也稱為滯回力。由于激勵振幅、頻率和預載荷的不同,滯回環的形狀有所不同。本文主要研究準靜態加載下的橡膠球鉸非線性滯回特性,不考慮振幅、頻率和預載等因素的變化對滯回特性產生的影響。
為了研究的方便,將滯回力分解為廣義彈性力和廣義阻尼力,如下式所示
(1)
對于線性剛度、線性阻尼系統,其滯回特性曲線為如圖1所示的橢圓。

圖1 線剛度、線阻尼系統滯回環曲線
引入復剛度,表達式如下
K*=K′+jK″=K′+jωC=K′(1+jη)
(2)
式中:K*為復剛度;K′為存儲剛度;K″為正交剛度;ω為圓頻率;C為阻尼系數;η為損耗因子,且η=tanφ,φ為復剛度和存儲剛度的夾角,如圖2所示。

圖2 復剛度示意圖
橢圓的斜率為剛度,每個循環的能量消耗即為橢圓的面積,即
(3)
在加載過程中,存儲的彈性能為三角形面積,即
(4)
損耗能量和加載的彈性能的比值為
(5)
故等效阻尼為
(6)
由式(6)可知,系統的等效阻尼系數與能量比、等效剛度和激勵頻率3個因素有關,對于任意非線性剛度、非線性阻尼系統,可通過數值積分方法求得滯回曲線包圍坐標軸的面積,從而計算能量比,再通過滯回環的斜率計算等效線剛度,就可以求得系統的等效阻尼系數。
在本文的研究中,引入等效線性阻尼模型,將式(1)變形為等效線性阻尼滯回特性模型
(7)
由于橡膠球鉸的滯回特性具有拉壓非對稱特性,由試驗可知,壓縮過程與拉伸過程的滯回環面積并不對稱。為提高模型精度,采用拉壓分段線性阻尼模型
(8)
式中:Ceq1為拉伸過程等效線性阻尼系數;Ceq2為壓縮過程等效線性阻尼系數。拉壓過程等效線性阻尼系數可由式(6)計算。
經典Bouc-Wen模型為
Fh=Fk(x)+z(x)
(9)
式中:z(x)為Bouc-Wen滯回力,通過一階微分方程求解;α、β、γ和n為Bouc-Wen模型的4個參數,α控制滯回環的斜率,β和γ控制滯回環的形狀,指數n控制滯回環的光滑度。
經典的Bouc-Wen模型只能描述對稱滯回曲線,利用分段擬合的思想,建立拉壓分段擬合的非對稱Bouc-Wen模型[9-10]如下
(10)
當x>0時,橡膠球鉸處于受拉狀態,i=1,Bouc-Wen模型的參數取值為α1、β1、γ1和n1;當x≤0時,橡膠球鉸處于受壓狀態,i=2,Bouc-Wen模型的參數取值為α2、β2、γ2和n2。
對于任意形狀的曲線,可通過高階多項式逼近,故采用多項式滯回力模型描述橡膠球鉸的非線性滯回特性,即
(11)
式中:n為多項式階次;ai、bi和ci分別為擬合的第i階位移項、速度項和交叉項的系數。
Savitzky-Golay濾波方法是一種局部多項式最小二次卷積擬合的數據平滑方法,通過卷積實現局部區間的多項式擬合,可以有效去除信號的高頻噪聲,是一種高效的信號低通濾波器[11]。Savitzky-Golay方法的原理如圖3所示。將采集到的含噪信號x[n]分為若干個區間,每個區間的信號點數N=2m+1,分別對每個區間內的點進行k階多項式擬合,擬合時通過多項式卷積快速計算多項式的系數,高效地獲得濾波后各點的數值,最后將濾波后的數據重新拼接在一起,形成去噪信號s[n]。

圖3 Savitiky-Golay濾波原理
應用Savitiky-Golay濾波,需要確定局部區間內的點數N和局部多項式擬合階數k這兩個重要參數。選取合適的N、k,不僅可以提高濾波去噪質量,減少信號的失真,同時還可以保證濾波效率,減少計算時間。因此,對于一組給定的信號,存在一組最優參數N、k,使得Savitzky-Golay濾波器的性能最優。
選定一個合適的去噪評價指標是進行濾波器參數優化的前提,現有的信號去噪評價指標有信噪比、信號熵、均方根等。這些評價指標在推力桿橡膠球鉸恢復力、位移和速度的去噪評價中并不適合,因為不同濾波器參數下的區分度不明顯。為此,本文引入灰色關聯度評價指標,對去噪前后的信號進行評價。灰色關聯度分析是我國學者鄧聚龍教授提出的灰色系統理論的重要組成部分[12-13],灰色關聯度描述的是序列之間的幾何相似程度,因此采用灰色關聯度評價去噪前后信號的接近程度,可以實現對濾波前后信號的宏觀幾何相似程度的判斷,從而保證信號的較小失真和相對較好的平滑效果。
去噪前后信號的灰色關聯度計算如下。
步驟1求原信號和濾波后信號的歸一化序列

(12)
式中:x′[i]和s′[i]分別為歸一化后的含噪信號和去噪信號;n為離散信號的總數。
步驟2求歸一化后的兩個信號序列之間的絕對差序列
Δ(i)=|s′[i]-x′[i]|, i=1,2,…,n
(13)
步驟3求絕對差序列的最值

(14)
步驟4計算關聯系數
(15)
式中δ為分辨系數,取δ=0.5。
步驟5求灰色關聯度
(16)

圖4 灰色Savitzky-Golay濾波方法流程圖
本文提出的灰色Savitzky-Golay濾波方法流程圖如圖4所示。對于給定的一組含噪信號x[n],通過調用布谷鳥算法優化程序,以去噪前后信號的灰色關聯度最大化為優化目標,應用Savitzky-Golay濾波算法對信號進行濾波,通過布谷鳥算法的迭代運算獲得一組最優的濾波器參數N、k對信號進行濾波,并輸出濾波后的信號s[n]。
在本文的研究中有兩個核心,即濾波器參數的調試和非線性遲滯模型的參數識別,這兩個問題的本質都是非線性優化問題。傳統的優化方法在效率和精度上有較大的缺陷,而智能算法在這類優化問題上具有明顯的優勢,因此本文選取布谷鳥算法這種智能算法進行濾波器的最優參數確定和非線性遲滯模型的參數識別。
布谷鳥算法由Yang等在2009年提出,是一種模擬布谷鳥尋找鳥巢放置鳥蛋以繁衍下一代的生物學行為的智能算法,算法采用Levy概率分布實現尋優軌跡的迭代更新[14-16]。布谷鳥算法構建在以下3個理想化條件的基礎上:①每只布谷鳥一次只產一個蛋,并隨機選擇鳥巢放置;②最好的鳥巢(解)能夠保留到下一迭代過程;③可用的鳥巢數是固定的,且以pa為巢主鳥發現外來鳥蛋的概率。通常,取pa=0.25。在鳥蛋被巢主鳥發現后,有丟棄寄生鳥蛋或者新建鳥巢兩種選擇。
算法的尋優路徑位置更新公式如下
?L(λ)
(17)

L~u=t-λ,λ∈(1,3]
(18)
其中u∈[0,1]且服從均勻分布。在更新位置坐標后,若u>pa,則重新隨機生成一組新位置,并計算適應度值,與上一次迭代的適應度值比較,保留適應度大的位置,通過不斷地迭代實現適應度函數的尋優,直至滿足收斂條件或達到最大迭代次數,停止迭代,輸出最優解。
在本文的研究中,推力桿試驗數據最優濾波參數優化的適應度函數為灰色關聯度的倒數,而參數識別時適應度函數為模型計算結果和試驗結果誤差的平方和,優化問題最終都轉化為極小尋優問題。
為研究商用車推力桿的縱向遲滯特性,試驗在電液伺服試驗臺上進行,試驗原理圖如圖5所示,試驗現場圖如圖6所示。試驗通過夾具將V型推力桿固定在試驗臺架上,通過控制電液伺服系統進行拉壓載荷的加載,并通過位移傳感器和力傳感器測量試驗過程中的位移和力信號。試驗時加載的載荷為準靜態斜坡力信號,載荷變化可分為4個階段,即0→150→0→-150→0 kN,完成一個拉壓周期的加載。一個周期共254 s,采樣頻率為10 Hz,故一個周期內采集到的力信號和位移信號的樣本點數為2 540。

圖5 V型推力桿滯回特性試驗原理圖

圖6 V型推力桿滯回特性試驗現場圖
試驗采集到的信號是含噪的力信號和位移信號,為了研究推力桿的縱向滯回特性,需要提供準確的力、位移和速度信號。因此,有必要對試驗采集到的力信號和位移信號進行濾波,去除信號的毛刺、尖峰和異常點,并通過位移信號的差分計算獲得速度信號。
采用本文提出的灰色Savitzky-Golay濾波方法對試驗采集到的位移信號進行濾波去噪。圖7所示是不同的濾波器參數N、k下去噪結果與原信號的灰色關聯度三維曲面,由圖可知在區間點數較大而局部擬合維度較小時,去噪前后信號的灰色關聯度最大。圖8所示是采用最優參數的灰色Savitzky-Golay濾波方法去噪后的位移信號,由圖可知,濾波器實現了對采集到的位移信號的有效濾波,在去除了信號中毛刺的同時,最大程度地保留了位移信號的宏觀曲線形狀。圖9所示是由原信號與去噪后信號相減獲得的位移噪聲信號,由圖可知,位移噪聲除了在載荷的拐點有尖峰外,其余部分為隨機噪聲。為進一步研究噪聲的分布情況,采用統計學方法進行分析。圖10是位移噪聲信號的平滑誤差頻數圖,由圖可知,噪聲信號呈現參數為N(0,σ2)的高斯正態分布,圖中曲線是位移噪聲的正態分布密度的擬合曲線。

圖7 不同參數對位移信號濾波效果的影響

圖8 濾波前后的位移信號對比圖

圖9 位移信號噪聲時域圖

圖10 位移信號噪聲統計特性
圖11所示是不同濾波器參數對去噪前后速度信號的灰色關聯度影響曲面圖,與位移信號濾波器參數有類似的規律。圖12所示是濾波后速度信號和直接對含噪位移信號差分獲得的速度信號對比圖,其中,一次濾波為對濾波后的位移信號進行差分計算的速度信號,二次濾波是在一次濾波獲得的速度信號的基礎上再次進行濾波后獲得的速度信號。由圖可知,一次濾波濾除了大部分的毛刺和尖峰,而二次濾波在一次濾波基礎上,進一步去除了位移信號中的隨機干擾噪聲。圖13和圖14所示分別是速度噪聲時域圖和頻數圖,由圖可知速度噪聲也服從高斯正態分布。

圖11 不同參數對速度信號濾波效果的影響

圖12 濾波前后的速度信號對比圖

圖13 速度信號噪聲時域圖

圖14 速度信號噪聲統計特性

圖15 不同參數對力信號濾波效果的影響

圖16 濾波前后的滯回力對比圖

圖17 滯回力信號噪聲時域圖
圖15所示是力信號濾波器參數與去噪前后信號灰色關聯度曲面,圖16是力信號濾波前后對比圖,圖17和圖18是力噪聲時域圖和頻數圖。力信號最優濾波器參數的規律和噪聲的統計分布特性與位移、速度噪聲類似,不再贅述。

圖18 滯回力信號噪聲統計特性
表1所示是采用布谷鳥算法獲得的最優化濾波器參數和噪聲信號分布特性的統計量。由表可知,對于不同的信號,最優濾波器參數不同,這是由噪聲信號和濾波方法共同決定的,采用本文所提出的布谷鳥算法優化后的灰色Savitzky-Golay濾波器能夠根據不同的含噪信號特征,采取最優的濾波參數,在去除噪聲的同時,保留信號的宏觀幾何特征,減少信號的失真。

表1 布谷鳥優化的濾波器參數
圖19所示是濾波前后的推力桿力-位移滯回特性曲線,由圖可知,濾波后的滯回環不存在明顯的毛刺和尖峰,而濾波前后滯回環的宏觀幾何特征很好地保留了下來。圖20所示是采用式(1)的廣義非線性彈性力和廣義非線性阻尼力分解出的彈性力和實測滯回力的對比圖,在分解出了廣義彈性力后,后文分別采用不同的模型對剩余的廣義非線性阻尼力進行擬合。

圖19 濾波前后的總體滯回特性曲線

圖20 分解出的非線性彈性力與實測滯回力的對比

圖21 等效線性模型對阻尼力的擬合效果
圖21是采用等效線性模型式(7)對阻尼力的擬合效果,由圖可知,等效線性模型使用一個等效的橢圓對不規則的遲滯阻尼環進行擬合,橢圓包圍的面積與阻尼滯回環包圍的面積相等。圖22和圖23分別是阻尼力和總體力-位移滯回環的擬合效果,由圖可知,等效線性模型可以對阻尼力和滯回環的趨勢進行擬合,但精度較差,且在等效計算過程中,將非對稱的遲滯特性平均等效而抹除了其非對稱特性。

圖22 等效線性模型仿真的阻尼力時域曲線

圖23 等效線性模型仿真的總體遲滯環

圖24 分段線性模型對阻尼力的擬合效果

圖25 分段線性模型仿真的阻尼力時域曲線
圖24所示是根據式(8)所述的分段線性模型進行的阻尼力的擬合效果,由圖可知,分段線性模型使用兩個等效的半橢圓對不規則的遲滯阻尼環進行擬合。兩個半橢圓包圍的面積分別與拉壓過程的阻尼滯回環包圍的面積相等。圖25和圖26所示分別是阻尼力和總體力-位移滯回環的擬合效果,由圖可知,分段線性模型可以對阻尼力和滯回環的趨勢進行較好地擬合,其精度較等效線性模型有所提高,且可以用于描述推力桿縱向拉壓的非對稱遲滯特性。

圖26 分段線性模型仿真的總體遲滯環

圖27 非對稱Bouc-Wen模型對阻尼力的擬合效果
采用布谷鳥算法,對非對稱Bouc-Wen模型的參數進行識別。圖27所示是根據式(10)所述的非對稱Bouc-Wen模型進行的阻尼力的擬合效果,由圖可知,非對稱Bouc-Wen模型使用兩個等效的半曲線對不規則的遲滯阻尼環進行擬合,兩條曲線包圍的面積分別與拉壓過程的阻尼滯回環包圍的面積相等。圖28和圖29所示分別是阻尼力和總體力-位移滯回環的擬合效果,由圖可知,非對稱Bouc-Wen模型可以對阻尼力和滯回環的趨勢和數值進行較好地擬合,其精度高于分段線性模型,且可以用于描述推力桿縱向拉壓的非對稱遲滯特性。

圖28 非對稱Bouc-Wen模型仿真的阻尼力時域曲線

圖29 非對稱Bouc-Wen模型仿真的總體遲滯環

圖30 多項式模型對阻尼力的擬合效果
應用式(11)的高階多項式模型對濾波后的數據進行模型參數識別,發現采用12階多項式可以對阻尼力進行很好的擬合。圖30所示是高階多項式模型擬合的阻尼力效果圖,由圖可知,高階多項式模型使用高階的多項式生成的不規則曲線對不規則的遲滯阻尼環進行高度逼近,兩條曲線包圍的面積分別與拉壓過程的阻尼滯回環包圍的面積相等。圖31和圖32所示分別是阻尼力和總體力-位移滯回環的擬合效果,由圖可知,高階多項式模型可以對阻尼力和滯回環的趨勢和數值進行很好地擬合,其精度高于其他3個模型,且可以用于描述推力桿縱向拉壓的非對稱遲滯特性,但由于模型是高階多項式模型,參數較多,參數識別難度大,增加了模型的復雜度。

圖31 多項式模型仿真的阻尼力時域曲線

圖32 多項式模型仿真的總體遲滯環
表2所示是4種不同模型的精度和復雜度對比情況。由表可知,隨著模型參數的增加,模型的精度有所提高。等效線性模型所需的參數數量最少,復雜度最小,而精度也較低。等效線性模型適合于定性地分析遲滯特性,不適合精度要求較高的仿真計算。在滿足一般工程計算精度要求的情況下,可以選取分段線性模型或者非對稱Bouc-Wen模型進行非對稱遲滯特性的仿真。若要求進一步提高計算精度,可采用高階多項式模型進行擬合。

表2 不同模型精度和復雜度對比
本文在分析了商用車推力桿橡膠球鉸實測縱向位移-力滯回特性的基礎上,對實測數據進行降噪處理,濾除數據中的毛刺和高頻噪聲等干擾項,以濾波去噪后的試驗數據和模型計算結果之間的誤差平方和最小為目標,應用布谷鳥算法對各模型的參數進行識別,主要結論如下。
(1)由試驗結果可知,商用車推力桿球鉸存在非對稱滯回特性,且具有高度的非線性。縱向拉壓過程中,橡膠球鉸剛度對稱,而阻尼不對稱,且壓縮阻尼大于拉伸阻尼。
(2)試驗數據存在較多的毛刺和尖峰,提出的灰色Savitzky-Golay濾波方法以原信號和去噪后信號的灰色關聯度最大化為目標,采用布谷鳥算法對濾波器的參數進行優化,在濾除了數據中的毛刺和尖峰等噪聲信號的同時,最大程度地保留了原信號的趨勢,減小了信號的失真。
(3)采用布谷鳥算法對非線性遲滯模型的參數進行識別,在滿足參數識別結果精度的同時,提高了參數識別的效率。
(4)等效線性模型和分段線性模型對非對稱遲滯特性的擬合精度較差,平均誤差大于15%,而非對稱Bouc-Wen模型和高階多項式模型可以很好地逼近試驗數據,平均誤差均小于10%,且高階多項式模型的誤差小于5%。