基于聚類和馬爾可夫鏈的貨車實際工況構(gòu)建方法研究

周英超,王保林,李連強,田 瑛,蒙心蕊

(1. 山東理工大學(xué) 交通與車輛工程學(xué)院,山東 淄博 255022; 2. 一汽解放青島汽車有限公司,山東 青島 266000)

0 引 言

車輛行駛工況是車輛在特定區(qū)域內(nèi)的行駛特性,用以判斷汽車的排放和油耗,對于新車型的技術(shù)發(fā)展及性能評價有一定的參考價值[1]。目前,國內(nèi)貨車油耗評價采用C-WTVC工況,由于不同地區(qū)之間地理環(huán)境、交通道路、車輛類型和城市規(guī)模的差異,導(dǎo)致法規(guī)工況與實際行駛工況相關(guān)性較差,使得二者在行駛特征和車輛能耗方面存在顯著差異,不利于車輛的精細化設(shè)計與性能評價[2]。因此,有必要根據(jù)車輛實際行駛數(shù)據(jù),構(gòu)建貨車實際工況,為車輛的精細化設(shè)計和實際道路整車性能評價提供條件和依據(jù)。

目前,行駛工況構(gòu)建研究主要采用短行程法、聚類法以及馬爾可夫鏈法(Markov chain,MC)等方法。短行程法將車輛行駛數(shù)據(jù)拆分為特征片段,然后將不同類型的特征片段連接形成短行程,選取其中代表性的部分,拼接合成最后工況[3]。短行程法采用時間隨機的單個片段,合成工況后時長不定。聚類法通過對比主要車輛參數(shù)貢獻率,對主要因素進行聚類分析,以使所構(gòu)建工況最大程度反映實際特征[4-6]。其中,K-均值聚類和模糊聚類方法容易受到局部最優(yōu)解的影響,而層次聚類算法具有較高的時間復(fù)雜度,其結(jié)果還受到聚類合并點和分裂點選擇的影響,這些因素都可能對行駛工況的構(gòu)建精度產(chǎn)生影響。馬爾可夫鏈法將工況看作一個馬爾可夫過程,以加速、減速、勻速及怠速4種狀態(tài)或車速區(qū)間作為狀態(tài)劃分標準,計算狀態(tài)轉(zhuǎn)移概率矩陣合成行駛工況[7-9]。馬爾可夫鏈法可以產(chǎn)生特定時間的運行工況,具有較高的構(gòu)建精度,但僅能對全局的數(shù)據(jù)進行分析,并不能對各種路況進行工況構(gòu)建。

因此,筆者提出一種基于高斯混合模型(Gaussian mixture model,GMM)聚類和馬爾可夫鏈法的工況構(gòu)建方法,這種方法不僅能夠劃分不同道路類型的運動學(xué)片段,還可以在每個類內(nèi)建立馬爾可夫工況,從而顯著提高了工況構(gòu)建的精度。筆者以目標貨車為研究對象,采用車聯(lián)網(wǎng)平臺獲取車輛實際行駛數(shù)據(jù),以速度劃分運動學(xué)片段,利用主成分分析對特征參數(shù)進行降維,并采用高斯混合模型對主成分進行分類,最后運用馬爾可夫鏈法來構(gòu)建目標貨車的實際行駛工況。

1 準備工作

1.1 數(shù)據(jù)采集

為構(gòu)建車輛實際行駛工況,選取目標貨車作為數(shù)據(jù)采集對象,采用自主駕駛法收集車輛行駛數(shù)據(jù),充分考慮不同駕駛員的駕駛行為和道路等級對行駛特征的影響。

采用車載終端及車聯(lián)網(wǎng)平臺進行行車數(shù)據(jù)采集,并將實時數(shù)據(jù)存儲于車內(nèi)服務(wù)器,并將所獲得的數(shù)據(jù)通過車載總線進行實時上傳至車聯(lián)網(wǎng)平臺,以便對數(shù)據(jù)進行調(diào)取和分析。為真實反映實際道路行駛特征,采樣頻率為1 Hz,連續(xù)采集時長為30 d,具體時間段選定為每天07:00-11:30、14:00-18:30、19:00-23:30,共收集試驗車輛共計1 592 280條行駛數(shù)據(jù)。

1.2 數(shù)據(jù)前期處理

由于通信網(wǎng)絡(luò)質(zhì)量、硬件條件等因素的影響,獲取的行駛數(shù)據(jù)存在信號波動、噪聲明顯等情況,需要進行降噪處理。小波去噪能夠?qū)崿F(xiàn)在各個尺度上的降噪,得到各個尺度的低、高頻成分,使其能夠最大限度地保持原有特征[10]。

小波降噪的具體過程(圖1)如下:

圖1 小波去噪原理Fig. 1 Principle of wavelet denoising

1)選擇適當(dāng)?shù)男〔ɑ头纸獬叨?對包含噪聲的車速和加速度信號進行小波變換,生成小波系數(shù)。

2)通過選擇合適的閾值函數(shù),對小波系數(shù)量化,進而更新小波系數(shù)。

3)利用新的小波系數(shù)進行小波逆變換,以重構(gòu)信號,從而得到去噪后的車速和加速度信號。

分別對部分原始數(shù)據(jù)進行3、5、7尺度下的小波去噪,濾波結(jié)果如圖2～圖4。

由圖3、圖4可知,7尺度下濾波加速度數(shù)據(jù)效果較好,濾波加速度結(jié)果與原始數(shù)據(jù)變化趨勢保持一致,并且能有效覆蓋原始數(shù)據(jù)的局部極值,因此選擇7尺度小波濾波結(jié)果作為參考車速。

圖3 各尺度加速度去噪結(jié)果Fig. 3 Acceleration denoising results of each scale

圖4 宏觀時間下7尺度分解結(jié)果Fig. 4 Seven scale decomposition results in macroscopic time

1.3 運動學(xué)片段劃分

以速度特征劃分運動學(xué)片段,取行駛數(shù)據(jù)中怠速到怠速為一個運動學(xué)片段,圖5為一個運動學(xué)片段示意。

圖5 運動學(xué)片段示意Fig. 5 Schematic diagram of kinematic segments

為反映運動學(xué)片段內(nèi)的瞬時特征,對速度v和加速度a進行狀態(tài)劃分:加速狀態(tài)(v>1 km/h且a>0.1 m/s2)、減速狀態(tài)(v>1 km/h且a<-0.1 m/s2)、勻速狀態(tài)(v>1 km/h且|a|≤0.1 m/s2)及怠速狀態(tài)(v≤1 km/h且|a|≤0.1 m/s2)。

1.4 特征參數(shù)選取與計算

表1 運動學(xué)片段主要特征參數(shù)值Table 1 Main characteristic parameter values of kinematic segments

2 聚類分析

2.1 主成分分析

表2 主要特征參數(shù)主成分分析結(jié)果Table 2 Principal component analysis results of main characteristic parameters

經(jīng)過主成分分析降維,十一維的運動學(xué)片段特征參數(shù)矩陣被降低至三維(表3)。

表3 主成分得分矩陣Table 3 Principal component score matrix

2.2 高斯混合模型

2.2.1 高斯混合模型

假設(shè)d維隨機變量x=(x1,x2,…,xd)T,包含K個子模型的高斯混合模型為:

(1)

(2)

(3)

式中:N(x|μg,Fg)為第g個子模型的高斯概率密度函數(shù);wg、μg、Fg分別為混合模型中第g個子模型的權(quán)重、均值、協(xié)方差矩陣;p(x)為高斯混合模型的概率密度函數(shù)。

2.2.2 EM算法

對于參數(shù)wg、μg、Fg,采用期望最大化算法(expectation maximum,EM)進行求解,其計算步驟分為E-step和M-step[11]:

E-step得到的數(shù)據(jù)xi(i=1,2,…,N)來自第g個子模型的概率為:

(4)

M-step計算新一輪迭代的模型參數(shù)見式(5):

(5)

重復(fù)計算式(4)和式(5)直至收斂。

2.2.3 最佳聚類數(shù)確定

利用貝葉斯信息準則判定最佳聚類數(shù)(Bayesian information criterion,BIC),即具有最低值WBIC模型對應(yīng)的K值為GMM聚類最佳聚類數(shù),計算公式為:

(6)

式中:k為模型中不同類型數(shù)據(jù)的個數(shù),取k=3;n為樣本中運動學(xué)片段數(shù)量,取n=2 516。

2.3 聚類結(jié)果分析

在MATLAB中基于Gmdistribution函數(shù)設(shè)置不同的K值,進行最佳聚類數(shù)即最小WBIC值的計算,結(jié)果如圖6。當(dāng)K=2時,WBIC值最大隨著K值變大,WBIC值呈現(xiàn)下降狀態(tài);當(dāng)K=9時,WBIC值達到最小;隨后WBIC值呈現(xiàn)緩慢上升狀態(tài)。因此確定的最佳聚類數(shù)應(yīng)為9。

圖6 不同聚類數(shù)下的WBIC值Fig. 6 WBIC values under different clustering numbers

在確定最佳聚類數(shù)后,對主成分得分進行高斯混合模型聚類,聚類效果如圖7。

圖7 高斯混合模型聚類效果Fig. 7 Clustering effect of Gaussian mixture model

將主成分得分聚類為9類,集合對應(yīng)的各類運動學(xué)片段,主要特征參數(shù)如表4。

表4 高斯混合模型聚類結(jié)果Table 4 Clustering results of Gaussian mixture model

3 馬爾可夫鏈法

3.1 馬爾可夫性驗證

在1、3、10、20 s的時間間隔下,對類4中車速數(shù)據(jù)進行相關(guān)性分析,如圖8。從圖8可以看出,車速數(shù)據(jù)相關(guān)性會隨著時間間隔增加逐漸減弱。由此可以證明,在短時間間隔下,實際道路行駛數(shù)據(jù)具備馬爾可夫特性,可采用馬爾可夫鏈法構(gòu)建貨車實際工況。

圖8 不同時間間隔下數(shù)據(jù)相關(guān)性Fig. 8 Data correlation under different time intervals

3.2 狀態(tài)劃分

以3 km/h的速度區(qū)間作為行駛狀態(tài)劃分標準,將各類速度區(qū)間分別劃分為9、11、13、23、22、26、30、34、34個狀態(tài),即每類看作是一個以{X1,X2,X3,…,Xn}為狀態(tài)空間的馬爾可夫鏈,馬爾可夫性可以表示為:

P(Xn+1=xn+1|X1=x1,…,Xn=xn)=P(Xn+1=xn+1|Xn=xn)

(7)

式中:Xn為狀態(tài)空間;P(Xn+1=xn+1|Xn=xn)為系統(tǒng)從i到j(luò)的狀態(tài)轉(zhuǎn)移概率。

3.3 狀態(tài)轉(zhuǎn)移概率

狀態(tài)轉(zhuǎn)移概率pij為系統(tǒng)當(dāng)前階段處于第i類狀態(tài)時下一階段轉(zhuǎn)移到第j類狀態(tài)的概率。常用的狀態(tài)轉(zhuǎn)移概率計算方法為最大似然估計法。該方法在樣本數(shù)據(jù)量較大的情況,計算效果較好,接近真實情況,但樣本數(shù)很小時,易導(dǎo)致概率為0,結(jié)果偏離實際,影響工況構(gòu)建精度。因此,筆者采用Laplace Smoothing方法[12]來估計狀態(tài)轉(zhuǎn)移概率,避免概率值為0的情況,如式(8):

(8)

在完成運動學(xué)片段分類后,獲得了各類運動學(xué)片段的集合,可以計算各類原始數(shù)據(jù)的狀態(tài)轉(zhuǎn)移概率矩陣,類3的結(jié)果如圖9。

圖9 狀態(tài)轉(zhuǎn)移概率矩陣直方圖Fig. 9 State transition probability matrix histogram

3.4 工況合成

結(jié)合狀態(tài)轉(zhuǎn)移概率矩陣的計算,按流程(圖10)構(gòu)建基于馬爾可夫鏈法的實際行駛工況。

圖10 基于馬爾可夫鏈法的實際行駛工況構(gòu)建流程Fig. 10 The construction process of actual driving condition based on Markov chain method

3.4.1 初始狀態(tài)確定

首先,將轉(zhuǎn)移概率最大的狀態(tài)作為馬爾可夫鏈的起始狀態(tài),在[0,1)區(qū)間內(nèi)產(chǎn)生服從于平均分配的隨機數(shù)r。然后,選擇最靠近這個隨機數(shù)的狀態(tài)轉(zhuǎn)換概率,也就是最小的差值絕對值,從而決定下一馬爾可夫狀態(tài)b:

(9)

式中:Pac為從狀態(tài)a轉(zhuǎn)移到狀態(tài)c的概率;c為下一狀態(tài)的序號;M為狀態(tài)總數(shù)。

3.4.2 構(gòu)建流程

重復(fù)3.4.1節(jié)步驟,直到滿足目標時長,獲得行駛狀態(tài)序列。在此基礎(chǔ)上,將車輛的運行狀態(tài)轉(zhuǎn)換成相應(yīng)的車速-時間曲線。因為該狀態(tài)只表示該速度值的范圍,所以需要根據(jù)式(10)來產(chǎn)生該速度值,即:

vt=[(ct-1)+rt] ×ds

(10)

式中:vt為t時刻狀態(tài)ct對應(yīng)的速度值;ct為t時刻的狀態(tài);r為t時刻MATLAB軟件均勻分布的隨機數(shù);ds為狀態(tài)對應(yīng)的速度區(qū)間長度,取ds=3。

3.4.3 代表工況篩選

由于隨機性的影響,馬爾可夫鏈法每次構(gòu)建行駛工況的結(jié)果均不相同,因此需要生成多條候選工況,并計算候選工況和原始數(shù)據(jù)主要特征參數(shù)的相對誤差ER和平均相對誤差EMR,并選取數(shù)值最小的候選行駛工況作為實際行駛工況。

(11)

式中:θcal為候選行駛工況的評估參數(shù)值;θe為原始試驗數(shù)據(jù)的評估參數(shù)值。

(12)

式中:ER,i為第i個評估參數(shù)的相對誤差;n為評價參數(shù)的數(shù)量。

建立2 000 s時長的工況,構(gòu)建工況由9類運動學(xué)片段組成。每類片段數(shù)目均為1,按照類1～類9的順序合成貨車實際行駛工況,如圖11。

圖11 貨車實際行駛工況Fig. 11 Actual driving conditions of trucks

4 結(jié)果分析

4.1 特征參數(shù)分析

將所構(gòu)建工況的特征參數(shù)分別與原始數(shù)據(jù)和貨車整車測試標準工況(C-WTVC)進行對比,特征參數(shù)結(jié)果如表5。

表5 特征參數(shù)對比Table 5 Comparison of characteristic parameters

所建工況同原始數(shù)據(jù)的相關(guān)特征參數(shù)進行對比,兩者的平均相對誤差為4.20%,所有特征參數(shù)的相對誤差均在10%范圍內(nèi),誤差微小且穩(wěn)定。因此,所構(gòu)工況能夠真實反映車輛的實際運動學(xué)行駛特征。

與C-WTVC工況相比,所建工況勻速比例增加了200%,加速比例減小了30.21%。由此可見,在車輛實際行駛過程中,法規(guī)工況與實際工況部分特征存在明顯差異,對于車輛的精細化控制與管理不利,開發(fā)實際行駛工況具有重要意義。

4.2 速度-加速度聯(lián)合概率密度分布分析

行駛工況中任一狀態(tài)均可通過速度和加速度進行描述,如圖12。由圖12可知,所建工況和原始數(shù)據(jù)的速度-加速度聯(lián)合概率分布(speed acceleration probability distribution , SAPD)整體趨勢比較一致,差異位于 -2～2 m/s2加速度區(qū)間內(nèi),分布誤差在 ±2%范圍內(nèi),二者差異較小。

圖12 速度-加速度聯(lián)合概率密度及誤差分布Fig. 12 Velocity-acceleration joint probability density and error distribution

4.3 仿真驗證

為了進一步驗證所建工況的有效性,作者選用Cruise軟件進行了仿真。試驗車輛的主要參數(shù)如表6。

表6 試驗車輛主要參數(shù)Table 6 Main parameters of the test vehicle

以法規(guī)工況和所建工況作為循環(huán)工況,進行仿真驗證,油耗如圖13。實際運行中,車輛的平均百公里油耗為22.86 L/100 km,而法規(guī)工況經(jīng)過轉(zhuǎn)換后得到的百公里油耗為16.84 L/100 km,兩者的相對誤差為26.33%。由此表明,法規(guī)工況不符合貨車的實際行駛特點。而所建工況經(jīng)過換算后可得百公里油耗為21.75 L/100 km,相對誤差為4.86%,符合車輛實際行駛特征。

圖13 仿真工況燃油消耗Fig. 13 Fuel consumption in simulation conditions

5 結(jié) 語

1)高斯混合模型聚類依據(jù)概率分配聚類成員,通過BIC信息準則確定最佳聚類數(shù),得到了較優(yōu)的聚類效果;馬爾可夫鏈法采用Laplace Smoothing方法來估計狀態(tài)轉(zhuǎn)移概率,建立了基于聚類分析的類內(nèi)工況,并充分考慮了車輛的瞬時運行特性。

2)通過主要的特征參數(shù)和Cruise平臺油耗仿真對比分析,國內(nèi)貨車整車循環(huán)測試所采用的C-WTVC法規(guī)工況并不能完全反映實際交通狀況的特征。因此,有必要建立能夠反映實際行駛特征的行駛工況。

3)構(gòu)建的貨車實際行駛工況是基于車聯(lián)網(wǎng)平臺的實際道路行駛數(shù)據(jù),能夠基本反映貨車實際道路行駛特征,可以為車輛的精細化設(shè)計和實際道路整車性能評價提供條件和依據(jù)。