基于模糊C均值聚類小波包重構(gòu)的車輛頻域非平穩(wěn)振動(dòng)譜模擬方法研究

叢 楠， 任焱晞， 陳俊達(dá)， 王彬星

(1.工程裝備系統(tǒng)工程研究所，北京 100093； 2.北京市標(biāo)準(zhǔn)化研究院，北京 100013)

采集車輛在實(shí)際道路上行駛時(shí)的振動(dòng)響應(yīng)，是準(zhǔn)確開(kāi)展道路模擬試驗(yàn)的重要前提。由于受到采集場(chǎng)地與設(shè)備的限制，在有限里程內(nèi)采集的實(shí)車響應(yīng)譜往往無(wú)法滿足在試驗(yàn)臺(tái)架上或虛擬仿真環(huán)境下模擬數(shù)百乃至數(shù)十萬(wàn)公里行程的需求。因此，在開(kāi)展模擬試驗(yàn)時(shí)，實(shí)際輸入振動(dòng)臺(tái)架作為激勵(lì)的往往是根據(jù)采集譜的某些特征進(jìn)行“重構(gòu)”處理后所得到的重構(gòu)譜。

按照重構(gòu)目標(biāo)特征的不同，對(duì)于車輛振動(dòng)譜的重構(gòu)方法可分為時(shí)域循環(huán)、頻域(功率譜)重構(gòu)和幅值概率密度重構(gòu)三個(gè)大類[1]。其中，以功率譜為目標(biāo)的重構(gòu)方法應(yīng)用最為廣泛，并形成了以諧波疊加法、線性濾波法、ARMA法為代表的一系列成熟方法[2-5]。然而，該類方法均以采集譜的總體功率譜為重構(gòu)目標(biāo)，因而無(wú)法重構(gòu)出非平穩(wěn)振動(dòng)譜。

隨著車輛設(shè)計(jì)與試驗(yàn)水平的不斷提高，對(duì)于如何準(zhǔn)確模擬車輛的各種非平穩(wěn)隨機(jī)振動(dòng)得到了越來(lái)越多的關(guān)注。王巖松等[6]由時(shí)變狀態(tài)方程建立一種在車速變化條件下的非平穩(wěn)振動(dòng)輸入模型；張立軍等[7]基于等效協(xié)方差法建立了一種通用的車輛非平穩(wěn)振動(dòng)輸入模型；萇亮等[8]由轉(zhuǎn)移函數(shù)法與狀態(tài)空間法對(duì)車輛對(duì)路面的非平穩(wěn)動(dòng)載荷進(jìn)行了模擬。上述方法均須對(duì)車輛振動(dòng)系統(tǒng)建模，并對(duì)該系統(tǒng)做線性假設(shè)，僅考慮由于車速變化造成路面激勵(lì)關(guān)于時(shí)間的非平穩(wěn)性(假定路面激勵(lì)關(guān)于空間是平穩(wěn)的)所導(dǎo)致的車輛振動(dòng)的在振動(dòng)強(qiáng)度上的非平穩(wěn)性。而事實(shí)上，由于車輛振動(dòng)系統(tǒng)是個(gè)復(fù)雜的非線性系統(tǒng)，車輛行駛時(shí)，即使路面激勵(lì)輸入是(關(guān)于時(shí)間)平穩(wěn)的，其上某一部位的振動(dòng)也會(huì)受到各種結(jié)構(gòu)非線性因素以及振動(dòng)傳遞路線多樣性的影響，在頻域上呈現(xiàn)出一定的非平穩(wěn)特征，即振動(dòng)的瞬時(shí)功率譜具有時(shí)變性。在時(shí)變瞬時(shí)功率譜信號(hào)的模擬方面，包括小波(包)、經(jīng)驗(yàn)?zāi)Ｊ椒纸狻⑾柌刈V分析在內(nèi)多種時(shí)頻分析方法，已被應(yīng)用于包括地震、風(fēng)載、水輪機(jī)以及車輛振動(dòng)的模擬，但這些方法大多要求功率譜的變化須滿足特定的隨機(jī)分布[9-11]，因此上述方法得到的重構(gòu)譜并不能認(rèn)為是嚴(yán)格意義上頻域非平穩(wěn)的，不能很好的與實(shí)際車輛振動(dòng)瞬時(shí)功率譜的變化特點(diǎn)相適應(yīng)。

本文通過(guò)對(duì)實(shí)車采集譜進(jìn)行小波包變換，將一維振動(dòng)譜的重構(gòu)問(wèn)題轉(zhuǎn)化為二維小波包系數(shù)矩陣(圖像)的重構(gòu)問(wèn)題。基于對(duì)小波包系數(shù)歸一化列向量的模糊C均值(Fuzzy C-Means，F(xiàn)CM)聚類，提取包括振動(dòng)譜所包含的主要振動(dòng)類型、各類型發(fā)生的總概率與相鄰條件概率、各類型中圍繞聚類中心的分散程度等采集譜在頻域中的非平穩(wěn)特征，并依據(jù)該特征實(shí)現(xiàn)對(duì)采集譜的頻域非平穩(wěn)化模擬。通過(guò)將該方法應(yīng)用于某試驗(yàn)道路響應(yīng)譜的重構(gòu)工作，所得到的重構(gòu)譜相比于原采集譜，在具有相似的頻域非平穩(wěn)特性的同時(shí)，其在時(shí)域波形連續(xù)性與平滑性、總體功率譜以及幅值概率密度等方面均具有良好的相似程度。本文為車輛振動(dòng)譜的重構(gòu)與模擬工作，提供了一種新的思路與方法。

1 振動(dòng)譜的小波包系數(shù)矩陣及其行列向量

1.1 小波包變換與小波包系數(shù)矩陣

小波變換是目前廣泛應(yīng)用于非平穩(wěn)信號(hào)分析與處理的時(shí)頻分析方法。在信號(hào)合成與重構(gòu)方面，其已在包括地震波、脈動(dòng)風(fēng)載、軌道不平度等諸多領(lǐng)域取得了成功的應(yīng)用[12-14]。

小波包變換由小波變換發(fā)展而來(lái)，其基本思想是在對(duì)小波分解后的高頻部分與低頻部分等同對(duì)待，均做進(jìn)一步分解，從而在全頻帶上具有一致的時(shí)頻分辨率。相比小波變換，更適合于對(duì)已經(jīng)過(guò)低通濾波等前處理的在全頻帶上均須做出準(zhǔn)確分析的車輛振動(dòng)譜進(jìn)行分解與合成。

對(duì)于長(zhǎng)度為N的離散信號(hào)，工程中通常使用Mallat快速(遞推)算法[15]來(lái)進(jìn)行小波包分解

(1)

相應(yīng)地，小波包重構(gòu)(反變換)的遞推算法為

(2)

(3)

通常，依據(jù)上述系數(shù)矩陣中系數(shù)取值的不同，以圖像的形式繪制出該矩陣可以更清晰的顯示出信號(hào)的時(shí)頻特征(由此又被稱為時(shí)頻圖像)。根據(jù)小波包變換性質(zhì)，對(duì)于采樣頻率Fs的信號(hào)進(jìn)行j層小波包分解，其時(shí)頻分辨率分別為

(4)

即隨著層數(shù)的增加，小波包分解的時(shí)間分辨率下降而頻率分辨率上升。

圖1(a)所示為經(jīng)過(guò)預(yù)處理的某型車駛過(guò)試驗(yàn)場(chǎng)塊石路時(shí)的單一輪軸振動(dòng)譜，采樣頻率為128 Hz，使用db6小波包對(duì)其進(jìn)行6級(jí)分解的系數(shù)圖像如圖1(b)所示，由式(4)可知，該圖像的時(shí)間/頻率分辨率分別為0.5 s/1 Hz。

(a) 振動(dòng)譜

(b) 小波包圖像圖1 某試驗(yàn)場(chǎng)塊石路的車輛振動(dòng)譜及其小波包系數(shù)圖像Fig.1 Vehicle vibration spectrum and its wavelet packet coefficients image of a proving ground gravel road

1.2 小波包變換與小波包系數(shù)矩陣

在式(3)所示的小波系數(shù)矩陣中，其第i行系數(shù)組成的行向量

ri={d[i,1],d[i,2],…，d[i,n]}

(5)

即為某一小波包系數(shù)，根據(jù)小波包分解的性質(zhì)，其代表了在頻帶[(i-1)Δf,iΔf]內(nèi)的信號(hào)能量隨時(shí)間的變化，且其上各系數(shù)的平方和

(6)

代表了信號(hào)在該頻帶內(nèi)的廣義能量。由帕賽瓦爾定理，其與信號(hào)的功率譜之間有如下對(duì)應(yīng)關(guān)系

(7)

式中：Cφ為一個(gè)與小波包基的類型有關(guān)的常數(shù)；P(f)為原信號(hào)的功率譜密度函數(shù)。在全頻帶上，亦有

(8)

式(7)、式(8)說(shuō)明了各小波包系數(shù)的平方和直接反映了該頻帶內(nèi)信號(hào)能量的強(qiáng)弱。同理，對(duì)于系數(shù)矩陣中的第k列系數(shù)向量及其平方和

ck={d[1,k],d[2,k],…，d[2j,k]}T

(9)

(10)

其分別代表了在[(k-1)Δt,kΔt]時(shí)段內(nèi)，信號(hào)的廣義瞬時(shí)功率譜與廣義瞬時(shí)能量。

此外，為便于進(jìn)行平穩(wěn)性分析，對(duì)于如式(5)、式(9)所示的小波包系數(shù)行/列系數(shù)向量，分別定義其歸一化版本

(11)

圖2給出了如圖1所示小波包系數(shù)進(jìn)行歸一化后的行向量與列向量圖像。

(a) 歸一化行向量

(b) 歸一化列向量圖2 小波包歸一化系數(shù)向量圖像Fig.2 Normalized wavelet packet coefficients image

2 基于FCM聚類小波包系數(shù)矩陣重構(gòu)的頻域非平穩(wěn)振動(dòng)模擬

正如前文所述，當(dāng)前討論關(guān)于車輛振動(dòng)信號(hào)的平穩(wěn)性的討論，大多單指其振動(dòng)強(qiáng)度關(guān)于時(shí)間的平穩(wěn)性。事實(shí)上，由圖2(a)中所示的歸一化行向量圖像即可看出，當(dāng)路面條件單一時(shí)，車輛響應(yīng)在所有頻帶內(nèi)的強(qiáng)度可認(rèn)為是平穩(wěn)的(此時(shí)振動(dòng)譜可順利通過(guò)平穩(wěn)性檢驗(yàn))。而圖2(b)所示的歸一化列向量圖像則表現(xiàn)出明顯的非平穩(wěn)特性，由歸一化列向量的含義，說(shuō)明該振動(dòng)譜的瞬時(shí)功率譜具有明顯的時(shí)變性。

2.1 歸一化列向量的FCM聚類

(12)

步驟2計(jì)算隸屬度

(13)

步驟3根據(jù)隸屬度修正聚類中心

(14)

由于FCM聚類方法的靈活性，其無(wú)需預(yù)先設(shè)定聚類中心的確切數(shù)量及數(shù)值(即瞬時(shí)功率譜的樣式)，僅需任意指定一個(gè)較大的聚類數(shù)量M(該值由受試車輛的主要激振源的數(shù)量決定，一般取10～20即可滿足要求)作為初始值，即可滿足對(duì)車輛瞬時(shí)功率譜的聚類要求。

2.2 基于FCM聚類結(jié)果的歸一化列向量非平穩(wěn)重構(gòu)

經(jīng)過(guò)FCM聚類，原有n組列向量將被歸入m(m≤M)個(gè)類別里，按照類別內(nèi)所包含列向量數(shù)量由多到少依次排序，各類別包含的數(shù)量為nci|i=1,2,…,m，其對(duì)應(yīng)的聚類中心(曲線)ccenter-i|i=1,2,…,m。為了進(jìn)行歸一化列向量重構(gòu)，須對(duì)聚類結(jié)果做進(jìn)一步統(tǒng)計(jì)：

(2) 統(tǒng)計(jì)每一類別中所屬的各歸一化列向量相對(duì)聚類中心的擴(kuò)散程度。由于各列向量及聚類中心均可表示為曲線，且類別中所包含的原列向量數(shù)量往往比較有限，因此，在不對(duì)擴(kuò)散的分布特征做前提假設(shè)時(shí)，本文使用百分位對(duì)曲線上的各點(diǎn)做單獨(dú)統(tǒng)計(jì)，即對(duì)于第i|1,2,…,m個(gè)分類上所有歸一化列向量上的第k|1,2,…,2j個(gè)值(對(duì)應(yīng)于相應(yīng)的列歸一化小波包系數(shù))，統(tǒng)計(jì)出若干個(gè)百分位點(diǎn)：{Qik-q1,Qik-q2,…,Qik-qn}|0

在完成上述對(duì)聚類結(jié)果的統(tǒng)計(jì)與分析后，即可對(duì)歸一化列向量進(jìn)行重構(gòu)。

重構(gòu)過(guò)程分為兩步：

(1) 指定各重構(gòu)歸一化列向量的類別序列，并使相鄰的概率為pij。具體可由如下遞推過(guò)程實(shí)現(xiàn)：

(i) 隨機(jī)指定一個(gè)類別i∈{1,2,…,m}；

(ii) 產(chǎn)生一個(gè)在[0,1]范圍內(nèi)的均勻分布的隨機(jī)數(shù)u；

(iii) 令i+1←min(j|pcum-ij≥u)；

(iv) 令i←i+1并重復(fù)第(ii)步，直至確定所有重構(gòu)列向量的分類序列。

(2) 按照既定的分類序列生成重構(gòu)歸一化列向量曲線，并使各重構(gòu)曲線滿足與原聚類結(jié)果中各類別中列向量圍繞其聚類中心的散布樣式相一致。對(duì)于基于百分位統(tǒng)計(jì)出的原分散結(jié)果，可由如下過(guò)程生成各重構(gòu)曲線：

(ii) 對(duì)應(yīng)于同一個(gè)分類中的第k|1,2,…,2j個(gè)值，在其各百分位點(diǎn)qn+1形成的qn-1個(gè)分段中按均勻分布隨機(jī)插值，同時(shí)產(chǎn)生Nci|i=1,2,…,m個(gè)隨機(jī)數(shù)；

(iii) 將所有2j組隨機(jī)數(shù)合并，則構(gòu)成了第i類的所有Nci|i=1,2，…，m個(gè)重構(gòu)歸一化列向量；

(iv) 令i←i+1并重復(fù)第(ii)步，直至生成所有類別所需的所有重構(gòu)列向量；

(v) 對(duì)所有最大值不為1的重構(gòu)列向量進(jìn)行校正，使其歸一化；

(vi) 按照已生成的重構(gòu)列向量的分類序號(hào)，從所生成的各類重構(gòu)歸一化列向量中依次抽取，生成重構(gòu)歸一化列向量矩陣。

2.3 小波包系數(shù)矩陣重構(gòu)

對(duì)于已生成的重構(gòu)歸一化列向量陣，僅需補(bǔ)充表示強(qiáng)度的信息，即可將其恢復(fù)為重構(gòu)小波包系數(shù)矩陣。由于由式(10)所示的廣義瞬時(shí)能量被認(rèn)為是平穩(wěn)的，因此利用信號(hào)小波包系數(shù)矩陣的該廣義能量信息，即可完成重構(gòu)小波包系數(shù)矩陣的工作，具體步驟如下：

步驟2將重構(gòu)廣義瞬時(shí)能量賦予各重構(gòu)的歸一化列向量

(15)

3 塊石路響應(yīng)譜重構(gòu)實(shí)例

對(duì)于如圖1所示的塊石路采集譜及其小波包系數(shù)圖像，將圖2(b)所示的歸一化列向量圖像用曲線集合的形式表示，如圖3所示。由小波包分解的含義，圖3即表示了原信號(hào)的時(shí)間分辨率為0.5 s、頻率分辨率為1 Hz的歸一化瞬時(shí)功率譜譜線的集合。

圖3 歸一化小波包列系數(shù)向量(局部)Fig.3 Normalized wavelet packet column coefficient vectors(partial)

通過(guò)FCM對(duì)如圖3所示的歸一化列向量曲線集合進(jìn)行聚類，形成了9個(gè)有效分類(有少量?jī)H包含1條列向量曲線的分類合并另作為一個(gè)類別，在重構(gòu)時(shí)可根據(jù)需要從中隨機(jī)抽取或忽略)。按照聚類中心峰值頻率、峰值頻率變化范圍以及中高頻成分強(qiáng)弱的不同，該9個(gè)有效分類又進(jìn)一步可劃分為4個(gè)大類，如圖4所示。

(a) 第I大類(第2、第3有效分類)

(b) 第II大類(第4、第6有效分類)

(c) 第III大類(第1、第7有效分類)

(d) 第IV大類(第5、第8、第9有效分類)圖4 聚類中心Fig.4 Centers of clustering

由圖4可知：第I類中心線具有尖銳的峰值，代表了較為理想的車輛振動(dòng)響應(yīng)，此時(shí)各種外部非線性因素影響較小，車輛運(yùn)行平穩(wěn)，中心曲線峰值頻率的輕微變化表明了車輛振動(dòng)系統(tǒng)結(jié)構(gòu)動(dòng)力學(xué)參數(shù)具有輕微的非線性變化；第II類中心曲線具有一強(qiáng)一弱兩處峰值，代表了由采集處振動(dòng)受到某一外部(此處為對(duì)側(cè)車輪)激勵(lì)的相對(duì)強(qiáng)烈的影響；第III類中心曲線具有較寬的主峰，其所包含的曲線的峰值頻率具有較大的變化范圍，表明采集處振動(dòng)受到多處外部激勵(lì)的影響以及振動(dòng)系統(tǒng)結(jié)構(gòu)動(dòng)力學(xué)參數(shù)的較強(qiáng)的非線性變化；第IV類中心曲線在中高頻段上的具有多處峰值，表明振動(dòng)中包含的間歇性中高頻或沖擊性激勵(lì)成分。

對(duì)各有效分類進(jìn)行統(tǒng)計(jì)，可得到如圖5所示的所屬曲線集合圍繞中心線的散布情況(以第3有效類為例，并以箱線圖的形式表示)，以及如表1所示的各有效類別之間的一階相鄰條件概率(以第1、第3、第5、第7個(gè)有效分類為例)。

圖5 第3有效分類中所屬各曲線圍繞聚類中心的散布情況Fig.5 Scatter of curves in valid category 3

表1 有效分類間的一階相鄰條件概率(部分)Tab.1 First-order adjacent conditional probabilities among valid categories (partial)

由表1可知，各有效分類之間并非均勻分布的隨機(jī)出現(xiàn)，而是呈較明顯的成組及連續(xù)出現(xiàn)的趨勢(shì)，從而顯現(xiàn)出明顯的頻域非平穩(wěn)特性。根據(jù)上文所述的重構(gòu)方法，首先由上述一階相鄰概率依次指定重構(gòu)列向量的類別，由全概率公式，當(dāng)重構(gòu)列向量數(shù)量足夠大時(shí)，重構(gòu)列向量的各類別的總體概率將與原信號(hào)保持一致。圖6展示了使用本文方法重構(gòu)出10倍于原信號(hào)長(zhǎng)度的重構(gòu)信號(hào)時(shí)，重構(gòu)信號(hào)及原信號(hào)中歸一化列向量各有效分類所占的比例。圖7對(duì)比了原信號(hào)的歸一化列向量與重構(gòu)歸一化列向量圖像。由圖6、圖7可知，對(duì)于歸一化列向量的按類別及散布特征重構(gòu)，可對(duì)車輛振動(dòng)中的頻域非平穩(wěn)特征進(jìn)行有效模擬。

圖6 原信號(hào)與重構(gòu)信號(hào)中歸一化列向量各有效分類的概率Fig.6 Probabilities of valid clustering categories of original vs reconstructed signal

圖7 原信號(hào)與重構(gòu)信號(hào)的歸一化列向量圖像Fig.7 Normalized column vector images of original vs reconstructed signal

原始采集譜及重構(gòu)的廣義能量序列如圖8(a)所示。當(dāng)原信號(hào)為強(qiáng)度平穩(wěn)時(shí)，可認(rèn)為該廣義能量序列為滿足某一分布(此處為對(duì)數(shù)正態(tài)分布)的隨機(jī)序列，根據(jù)其生成的重構(gòu)廣義能量序列及兩信號(hào)廣義能量序列的概率密度分布分別如圖8(b)、圖9所示。將該重構(gòu)廣義能量序列賦予圖7所示的重構(gòu)歸一化向量，得到的重構(gòu)小波包系數(shù)圖像如圖10所示。

(a) 原信號(hào)(b) 重構(gòu)信號(hào)圖8 廣義能量序列Fig.8 Generalized energy sequence圖9 廣義能量序列的幅值概率密度Fig.9 Probability density of generalized energy se-quences 圖10 重構(gòu)小波包系數(shù)圖像Fig.10 Wavelet packet coeffi-cients image of recon-structed signal

基于式(2)所示的小波包重構(gòu)，將重構(gòu)小波系數(shù)矩陣進(jìn)行小波包重構(gòu)，即可得到重構(gòu)振動(dòng)譜。該重構(gòu)譜與原始譜在時(shí)域波形、總體功率譜以及幅值概率密度的對(duì)比情況分別如圖11～圖13所示。

(a) 總體波形

(b) 局部放大圖11 不同尺度下的時(shí)域波形Fig.11 Time domain waveform in different scales

圖12 總體功率譜密度Fig.12 Gross PSDs

圖13 幅值概率密度分布Fig.13 Amplitude PDFs

通過(guò)圖11～圖13可知，采用本文方法得到的重構(gòu)譜在具有與原采集譜相似的非平穩(wěn)特性的基礎(chǔ)上，同時(shí)在時(shí)域波形、總體功率譜及幅值概率密度上與原始采集譜具有較高的一致性。并且，根據(jù)dB小波包自身性質(zhì)，通過(guò)合理選擇小波包消失矩的階數(shù)，可使得重構(gòu)譜在保持原采集譜的平滑性連續(xù)性與功率譜一致性之間，進(jìn)行靈活的取舍——對(duì)于包含更多沖擊樣高幅值激勵(lì)且平滑性差的振動(dòng)譜，可選用消失矩階數(shù)較低的小波包；相反，對(duì)于較為平緩的振動(dòng)譜，應(yīng)選用消失矩較高的小波包。

相比于傳統(tǒng)振動(dòng)譜的重構(gòu)方法，本文方法得到的重構(gòu)振動(dòng)譜能夠更真實(shí)地反映出車輛的實(shí)際振動(dòng)情況，為開(kāi)展各種車輛振動(dòng)模擬試驗(yàn)提供更為精確的振動(dòng)輸入。

4 結(jié) 論

受到地面-車輛系統(tǒng)中多種因素的影響，實(shí)測(cè)車輛振動(dòng)譜在頻域上往往是非平穩(wěn)的；該非平穩(wěn)特征體現(xiàn)為若干種類別振動(dòng)模式的間歇性交替出現(xiàn)并相互交疊；通過(guò)對(duì)振動(dòng)譜小波包系數(shù)歸一化列向量進(jìn)行FCM聚類，可有效地對(duì)各振動(dòng)模式(及其主要交疊模式)進(jìn)行提取與分析；通過(guò)重構(gòu)與聚類結(jié)果相一致的歸一化列向量進(jìn)而生成重構(gòu)系數(shù)矩陣，經(jīng)由小波包逆變換，可得到與原采集譜具有相似振動(dòng)強(qiáng)度、總體功率譜及幅值概率密度，且同時(shí)擁有相似頻域非平穩(wěn)特性的重構(gòu)振動(dòng)譜；通過(guò)選取合適的小波包類型及分解層數(shù)，原采集譜的平滑性與連續(xù)性也可得到有效的模擬。

通過(guò)將該方法應(yīng)用于某試驗(yàn)場(chǎng)塊石路采集譜的重構(gòu)工作，驗(yàn)證了本文方法在車輛振動(dòng)譜模擬方面的正確性和有效性。此外，由于小波包分析及FCM聚類方法本身對(duì)振動(dòng)譜分析的普遍適用性，本文所述方法亦可應(yīng)用于其它非平穩(wěn)振動(dòng)譜的模擬工作，為各種振動(dòng)模擬試驗(yàn)提供與實(shí)際更為吻合的振動(dòng)輸入。