組合優(yōu)化聚類與馬爾科夫鏈的城市環(huán)衛(wèi)車輛行駛工況構(gòu)建方法

彭育輝，莊 源

(1. 福州大學(xué)機(jī)械工程及自動(dòng)化學(xué)院，福建 福州 350108; 2. 廈門理工學(xué)院福建省客車及特種車輛研發(fā)協(xié)同創(chuàng)新中心，福建 廈門 361024)

0 引言

汽車行駛工況是用于汽車污染物排放、 燃油消耗等性能評(píng)估的基礎(chǔ)數(shù)據(jù)，是影響汽車動(dòng)力匹配與標(biāo)定的關(guān)鍵因素. 截至2016年底，我國環(huán)衛(wèi)車輛保有量達(dá)到34.75萬輛，作為專用作業(yè)車輛，其工作范圍覆蓋絕大部分城市區(qū)域，工作時(shí)間長(zhǎng)，且行駛工況特點(diǎn)明顯區(qū)別于其它類別車輛. 美國能源部對(duì)城市垃圾轉(zhuǎn)運(yùn)車的運(yùn)行工況進(jìn)行研究，并開發(fā)了測(cè)試工況(neighborhood refuse truck cycle, NRTC)，模擬車輛在各垃圾站點(diǎn)間的頻繁啟停工況； Clark等[1-3]根據(jù)紐約市垃圾車的行駛特征開發(fā)的New York garbage truck cycle(NYGTC)也被應(yīng)用于紐約市環(huán)衛(wèi)車的開發(fā)測(cè)試. 但是目前絕大部分國家并沒有建立針對(duì)城市環(huán)衛(wèi)車的專用測(cè)試行駛工況，而是普遍采用商用車的測(cè)試工況. 中國現(xiàn)階段對(duì)城市環(huán)衛(wèi)車的測(cè)試工況是根據(jù)商用車的底盤屬性，采用由世界重型商用車輛瞬態(tài)循環(huán)(WTVC，world transient vehicle cycle)為基礎(chǔ)修訂的C-WTVC工況進(jìn)行測(cè)試[4]. 當(dāng)車輛的最大設(shè)計(jì)總質(zhì)量GVW大于5.5 t，用于測(cè)試油耗的特征里程分配中，市區(qū)比例10%，公路比例占30%～60%，高速占比工況30%～60%，顯然，這與城市環(huán)衛(wèi)車的實(shí)際運(yùn)行狀況差距較大.

國內(nèi)外研究學(xué)者對(duì)專用汽車行駛工況的構(gòu)建方法做了大量研究. 文[5]利用多參數(shù)統(tǒng)計(jì)理論構(gòu)建自卸車循環(huán)工況； 文[6]以聚類分析方法構(gòu)建重型載貨汽車行駛工況； 文[7-8]利用大樣本數(shù)據(jù)，通過短行程的劃分，采用“最佳增量法”分別構(gòu)建了重型客車與城市公交車的行駛工況. 文[9]最早提出馬爾科夫法并應(yīng)用于美國加利福尼亞州的LA01工況構(gòu)建上，文[9-11]進(jìn)一步將車輛行進(jìn)過程抽象為離散隨機(jī)的馬爾科夫過程，基于狀態(tài)轉(zhuǎn)移概率合成車輛行駛工況. 相對(duì)而言，目前國內(nèi)對(duì)于環(huán)衛(wèi)車行駛工況的研究涉及極少，文[12]是以一輛環(huán)衛(wèi)車行駛數(shù)據(jù)為研究對(duì)象，采用兩階段聚類法構(gòu)建了西安市環(huán)衛(wèi)車行駛工況，但是其采集區(qū)域與路線固定，試驗(yàn)車輛單一，采集時(shí)段覆蓋較少，使得構(gòu)建結(jié)果的代表性不夠理想.

為研究城市環(huán)衛(wèi)車輛的實(shí)際工況特征和合理的工況構(gòu)建方法，以福州市內(nèi)3輛清潔環(huán)衛(wèi)車持續(xù)一周各時(shí)段的實(shí)際道路行駛數(shù)據(jù)為研究對(duì)象，提出一種適用于城市環(huán)衛(wèi)車輛的工況構(gòu)建方法. 通過均值化處理改進(jìn)傳統(tǒng)主成分分析 (principal component analysis，PCA)，提高各主成分提取到的車輛行駛信息，對(duì)特征參數(shù)降維后構(gòu)造得到的各主成分進(jìn)行k均值聚類，將原始數(shù)據(jù)樣本劃分為兩類，有效區(qū)分環(huán)衛(wèi)車在清潔作業(yè)時(shí)(以穩(wěn)定中低車速行駛，加減速平緩)和作業(yè)結(jié)束返程時(shí)存在的兩種明顯不同的行駛特征； 再將每一類行駛過程視為一個(gè)隨時(shí)間變化的馬爾科夫過程，基于馬爾科夫模型構(gòu)建其代表工況，充分體現(xiàn)其瞬時(shí)工況的隨機(jī)性. 所構(gòu)建的福州市環(huán)衛(wèi)車代表工況相較于國家現(xiàn)行標(biāo)準(zhǔn)工況能更好地反映福州市環(huán)衛(wèi)車實(shí)際運(yùn)行狀況，更適用于福州市環(huán)衛(wèi)車輛的開發(fā)測(cè)試.

1 數(shù)據(jù)采集與分析

1.1 實(shí)測(cè)數(shù)據(jù)采集

將車載數(shù)據(jù)采集終端接入車輛OBD接口，以1 Hz的頻率實(shí)時(shí)采集測(cè)試車輛在行駛過程中的速度、 GPS位置、 瞬時(shí)油耗等參數(shù)，并與數(shù)據(jù)監(jiān)控平臺(tái)進(jìn)行通訊，傳輸采集數(shù)據(jù). 選取3輛福州市清潔環(huán)衛(wèi)車輛，采用自主駕駛法，按照其正常工作路線行駛采集數(shù)據(jù)，而非規(guī)劃的固定試驗(yàn)路線，保證采集的數(shù)據(jù)能真實(shí)有效地反映福州市環(huán)衛(wèi)車的實(shí)際運(yùn)行工況特點(diǎn). 車輛行駛區(qū)域覆蓋福州市各主要干線道路，數(shù)據(jù)采集周期為連續(xù)一周，含工作日和周末非工作日的各個(gè)時(shí)段. 采集到總計(jì)約180 ks有效的車輛實(shí)際道路行駛數(shù)據(jù).

1.2 實(shí)測(cè)數(shù)據(jù)預(yù)處理及分析

圖1 實(shí)測(cè)數(shù)據(jù)速度-加速度頻率分布圖Fig.1 Frequency of speed-acceleration distribution

車輛行駛數(shù)據(jù)在采集過程中，由于采集設(shè)備故障、 環(huán)境條件等客觀因素導(dǎo)致采集到的數(shù)據(jù)存在異常值，需對(duì)原始采集數(shù)據(jù)做濾波處理，消除影響. 在實(shí)際行駛過程中，受環(huán)衛(wèi)車輛動(dòng)力與制動(dòng)性能、 道路交通條件和環(huán)衛(wèi)作業(yè)規(guī)定等限制，對(duì)加速度超過-4～4 m·s-2閾值區(qū)間、 存在缺失值、 加速或減速時(shí)間比例為0的環(huán)衛(wèi)車輛行駛數(shù)據(jù)，可判定為是由于GPS振動(dòng)漂移或信號(hào)屏蔽等原因?qū)е碌臄?shù)據(jù)異常，此類數(shù)據(jù)片段做剔除處理； 剔除異常值后車輛實(shí)測(cè)數(shù)據(jù)的速度-加速度頻率分布情況如圖1所示，行駛數(shù)據(jù)主要分布在0～40 km·h-1車速，-1～1 m·s-2加速度的區(qū)間內(nèi)，具有車速較低，加減速平緩等專用作業(yè)車輛的行駛特點(diǎn).

2 環(huán)衛(wèi)車代表工況構(gòu)建

2.1 運(yùn)動(dòng)學(xué)片段劃分

圖2 運(yùn)動(dòng)學(xué)片段Fig.2 Micro-trips

車輛行進(jìn)過程中，因路況和交通因素，存在多次怠速、 加速、 減速與勻速的狀態(tài). 如圖2所示，運(yùn)動(dòng)學(xué)片段為某一怠速開始至下一怠速開始前的區(qū)間[13]. 汽車行進(jìn)過程可視為大量運(yùn)動(dòng)學(xué)片段的銜接，按這一原則對(duì)原始工況數(shù)據(jù)進(jìn)行劃分，共得到2 070個(gè)運(yùn)動(dòng)學(xué)片段. 為最大限度反映每個(gè)運(yùn)動(dòng)學(xué)片段的行駛特征，選取如表1所示的15個(gè)特征參數(shù)來描述運(yùn)動(dòng)學(xué)片段，并計(jì)算出每一片段的特征參數(shù)值，構(gòu)造為特征參數(shù)矩陣.

表1 用于描述運(yùn)動(dòng)學(xué)片段的15個(gè)特征參數(shù)

2.2 基于主成分分析法的特征參數(shù)降維

各特征參數(shù)因量綱不同將導(dǎo)致各自取值較為離散，使得大方差變量在后續(xù)分析時(shí)會(huì)優(yōu)先照顧，降低計(jì)算準(zhǔn)確性，因此特征參數(shù)降維前需先做標(biāo)準(zhǔn)化處理，消除量綱不同造成的影響[14]. 構(gòu)造包含15個(gè)特征參數(shù)的矩陣A為:

(1)

式中：aij(i=1, 2, …,m；j=1, 2, …,n)是第i個(gè)運(yùn)動(dòng)學(xué)片段的第j個(gè)參數(shù).

在主成分分析前先求出矩陣A的協(xié)方差矩陣A′，對(duì)A′做均值化處理，避免各參數(shù)變異程度信息丟失. 令：

(2)

即用協(xié)方差矩陣A′中的每一元素aij除以所在列的所有元素均值，得到均值化后的特征參數(shù)矩陣Bm×n. 統(tǒng)計(jì)學(xué)上已證明[15]，經(jīng)過均值化處理后的矩陣B將不改變各參數(shù)變量的相關(guān)系數(shù)，且能全面地包含原始數(shù)據(jù)信息.

對(duì)均值化處理后的特征參數(shù)矩陣做主成分分析，利用線性變換構(gòu)造出新變量(即主成分)，使得各主成分間兩兩獨(dú)立且線性無關(guān)，消除了所選取的各特征參數(shù)在描述運(yùn)動(dòng)學(xué)片段時(shí)存在的信息重疊，同時(shí)保留了原變量的絕大部分特征信息. 研究表明，若主成分的特征值大于1，且累積貢獻(xiàn)率在80%以上，即滿足行駛工況開發(fā)要求[16]. 特征參數(shù)矩陣在主成分分析后得到的各主成分相應(yīng)特征值、 貢獻(xiàn)率及累計(jì)貢獻(xiàn)率如表2所示，可以看到，前三個(gè)主成分的特征值均大于1，累積貢獻(xiàn)率為85.332%. 因此選取前三個(gè)主成分來表征15個(gè)特征參數(shù)所包含的行駛信息，能達(dá)到降維目的.

表2 主成分貢獻(xiàn)率及累積貢獻(xiàn)率

2.3 k均值聚類分析

聚類的過程就是把樣本集分類成不同組別的過程.k均值聚類是非監(jiān)督快速聚類，能夠使樣本分類后達(dá)到類內(nèi)相似度高，類間相似度低的目的[17]. 具體過程為：

1) 針對(duì)樣本數(shù)據(jù)特點(diǎn)確定k個(gè)初始聚類中心.

2)計(jì)算每一樣本xi中的第k個(gè)變量值到最近聚類中心zj的歐式距離，確定樣本類別； 歐式距離能很好表征類別間不相似性. 歐式距離的計(jì)算公式如下：

(3)

3) 重新計(jì)算求出每一聚類的均值，將這一均值作為新聚類中心.

4) 迭代步驟2)～3)，直至聚類中心收斂為定值.

按上述步驟，對(duì)主成分進(jìn)行k均值聚類，把行駛特征相似的運(yùn)動(dòng)學(xué)片段聚為同一類，獲得了兩類片段樣本集，兩類片段樣本各特征參數(shù)值的對(duì)比如表3所示. 可以看到，兩類片段樣本集的特征參數(shù)差異明顯，分別代表兩類不同的行駛特征. 聚類1、 聚類2的樣本時(shí)間長(zhǎng)度比例分別為69.60%和30.40%，其中聚類1的平均車速較低，僅為16.45 km·h-1，同時(shí)勻速比例Pc達(dá)到40.3%，加速比例Pa、 怠速比例Pi比例和平均加速度都不高，表明聚類1樣本為環(huán)衛(wèi)車在進(jìn)行清潔作業(yè)時(shí)，緩慢穩(wěn)定行駛的片段樣本集； 聚類2的加速狀態(tài)比例達(dá)到37.3%，平均車速也較高，而勻速時(shí)間比例僅為18.3%，表明聚類2樣本為環(huán)衛(wèi)車完成清潔作業(yè)后，跟隨城市車流駛回停車場(chǎng)的片段樣本集.

表3 每一聚類的特征參數(shù)

2.4 基于馬爾科夫鏈法合成代表工況

環(huán)衛(wèi)車在行進(jìn)過程中，車速是隨城市車流根據(jù)交通情況隨機(jī)變化的，因此將其行駛過程視為一個(gè)隨時(shí)間變化的離散馬爾科夫過程，能夠充分體現(xiàn)其瞬時(shí)工況的隨機(jī)性.

2.4.1狀態(tài)劃分

由于行車過程是在加速、 減速、 勻速、 怠速工況間不斷輪換，基于馬爾科夫模型進(jìn)行汽車行駛工況構(gòu)建時(shí)，采用這四種確定的工況對(duì)實(shí)測(cè)數(shù)據(jù)的速度-時(shí)間序列進(jìn)行劃分，每個(gè)被劃分出的片段定義為模型事件，記為Zτ={Z1,Z2, …,Zτ| (τ=1, 2, …,T)}，再將運(yùn)動(dòng)學(xué)特征接近的片段整合為模型事件集，形成馬爾科夫模型中的“狀態(tài)”[18]，記為S={S1，S2，…，Sτ| (τ=1, 2, …,k)}. 狀態(tài)劃分后，在所構(gòu)建馬爾科夫模型中，即能以狀態(tài)(隱性信息)表征車速(顯性信息)，因此狀態(tài)劃分是工況構(gòu)建的首要.

2.4.2狀態(tài)轉(zhuǎn)移矩陣的計(jì)算

對(duì)于所有的模型狀態(tài)S1,S2, …,Sτ，狀態(tài)Sτ的概率僅和前一狀態(tài)Sτ-1相關(guān)聯(lián)，離散的馬爾科夫鏈在固定時(shí)域下具有穩(wěn)定性和齊次性，未來時(shí)刻的取值只與當(dāng)前時(shí)刻取值有關(guān)，而與過去時(shí)刻取值無關(guān)[19].

對(duì)于從任意τ-1時(shí)刻到τ時(shí)刻，其條件概率有：

P={Zτ=s|Zτ-1=r}=P{Zτ+τ′=s|Zτ+τ′-1=r}(τ′=-(τ-1),-(τ-2), …, -1, 0, 1, …)

(4)

利用最大似然估計(jì)法計(jì)算各個(gè)狀態(tài)間的轉(zhuǎn)移概率Prs，全部轉(zhuǎn)移概率構(gòu)造為轉(zhuǎn)移概率矩陣.

(5)

其中：Nrs指時(shí)間從τ-1到τ，狀態(tài)從r轉(zhuǎn)移至s的事件數(shù).

基于馬爾科夫模型計(jì)算時(shí)，采用隨機(jī)選擇連續(xù)片段的傳統(tǒng)方法不能確保每個(gè)模型事件在工況中的最佳代表性. 因此本文利用最大似然估計(jì)法，以估計(jì)狀態(tài)轉(zhuǎn)移概率的方法取代傳統(tǒng)的隨機(jī)片段法，在車速隱馬爾科夫鏈中，通過轉(zhuǎn)移概率來不斷預(yù)測(cè)下一狀態(tài)，保證選取模型事件的代表性，各狀態(tài)間轉(zhuǎn)移概率矩陣計(jì)算結(jié)果見表4.

表4 各狀態(tài)間轉(zhuǎn)移概率矩陣

2.4.3環(huán)衛(wèi)車代表工況合成

按照上文提及的馬爾科夫法基礎(chǔ)理論，以下述原則合成候選工況.

1) 根據(jù)聚類1、 聚類2占原始工況數(shù)據(jù)的時(shí)間比例，計(jì)算出環(huán)衛(wèi)車在兩類工作模式下的候選工況長(zhǎng)度；

2) 為便于行駛工況在臺(tái)架測(cè)試中的應(yīng)用，以怠速工況狀態(tài)作為初始狀態(tài)，而后的每一個(gè)行駛狀態(tài)是根據(jù)轉(zhuǎn)移概率矩陣和當(dāng)前時(shí)刻行駛狀態(tài)確定；

3) 狀態(tài)轉(zhuǎn)移矩陣各個(gè)狀態(tài)間的轉(zhuǎn)移概率≥0；

圖3 清潔環(huán)衛(wèi)車代表行駛工況 Fig.3 Representative driving cycle of cleaning trucks

4) 由最大似然估計(jì)法不斷進(jìn)行下一狀態(tài)篩選，直至合成候選工況達(dá)到預(yù)定長(zhǎng)度.

將合成的聚類1、 聚類2候選工況與實(shí)測(cè)數(shù)據(jù)每一個(gè)狀態(tài)的行轉(zhuǎn)移概率做K-S雙尾檢驗(yàn)，篩選出相似性水平不小于0.9的候選工況，以確保二者的狀態(tài)轉(zhuǎn)移矩陣來自同一樣本分布. 將兩類候選工況隨機(jī)排列組合，篩選出特征參數(shù)平均誤差最小的合成工況作為最終代表工況. 最終代表工況中兩類工況的時(shí)間占比與兩種聚類樣本在實(shí)測(cè)數(shù)據(jù)中的時(shí)間占比相同. 按時(shí)間比例合成的福州市清潔環(huán)衛(wèi)車代表工況如圖3所示.

3 行駛工況有效性檢驗(yàn)

3.1 基于運(yùn)動(dòng)學(xué)參數(shù)的有效性檢驗(yàn)

為驗(yàn)證代表工況的準(zhǔn)確性，將所構(gòu)建的福州市環(huán)衛(wèi)車代表工況的主要特征參數(shù)與實(shí)際試驗(yàn)采集數(shù)據(jù)的主要特征參數(shù)進(jìn)行對(duì)比，其相對(duì)誤差結(jié)果如表5所示. 代表工況與實(shí)際試驗(yàn)數(shù)據(jù)相比，各特征參數(shù)的相對(duì)誤差均在7%以下，相對(duì)誤差絕對(duì)值的算術(shù)平均值為4.46%. 同時(shí)，代表工況中的片段運(yùn)動(dòng)學(xué)特征和車輛實(shí)際行駛狀況基本一致，符合環(huán)衛(wèi)車在實(shí)際行進(jìn)過程中的行駛特點(diǎn).

圖4 試驗(yàn)數(shù)據(jù)與代表工況速度-加速度聯(lián)合分布差異率 Fig.4 SAFD difference between the sampling data and the representative driving cycle

其次，從速度和加速度兩個(gè)角度驗(yàn)證代表工況和試驗(yàn)數(shù)據(jù)之間的差異率，采用速度-加速度聯(lián)合分布圖(speed-acceleration frequency distribution, SAFD)來驗(yàn)證代表工況的準(zhǔn)確性, 如圖4所示. 從圖4可見，試驗(yàn)數(shù)據(jù)與代表工況數(shù)據(jù)在速度加速度各個(gè)分布區(qū)間上的差異率都維持在±5%以內(nèi)，說明所構(gòu)建代表工況能很好地反映試驗(yàn)數(shù)據(jù)樣本，表征車輛的實(shí)際行駛特征.

3.2 工況對(duì)比分析

將構(gòu)建的福州市清潔環(huán)衛(wèi)車代表工況(Fuzhou sanitation truck cycle, FZSTC)與現(xiàn)階段我國對(duì)城市環(huán)衛(wèi)車進(jìn)行測(cè)試的C-WTVC標(biāo)準(zhǔn)工況進(jìn)行對(duì)比，如表6所示，其中FZSTC-1為清潔作業(yè)工況，F(xiàn)ZSTC-2為回程工況. 可以看到，環(huán)衛(wèi)車在清潔作業(yè)時(shí)，勻速時(shí)間比例Pc較大，達(dá)到39.19%，而在回程時(shí)加減速較頻繁，加速時(shí)間比例Pa能達(dá)到40.29%. FZSTC與C-WTVC在平均車速va和最高車速vmax上差異明顯，主要原因是環(huán)衛(wèi)車在清潔作業(yè)時(shí)的特殊性要求，必須保持穩(wěn)定的較低車速，同時(shí)為了專用車輛的行車安全，即使在回程時(shí)，環(huán)衛(wèi)部門也規(guī)定環(huán)衛(wèi)車不能以較高的車速行駛，因此FZSTC整個(gè)循環(huán)不存在高速段，最高車速僅有58.12 km·h-1. 考慮到實(shí)際行駛工況特征，將C-WTVC測(cè)試循環(huán)去除432 s高速段后的工況命名為C-WTVC*，考察其與FZSTC的差異. 結(jié)果表明，C-WTVC和C-WTVC*都與FZSTC差異較大，尤其在勻速時(shí)間比例Pc和平均車速va差異明顯. 這說明，現(xiàn)行對(duì)環(huán)衛(wèi)車測(cè)試采用的商用車標(biāo)準(zhǔn)測(cè)試循環(huán)C-WTVC并不能有效反映出城市環(huán)衛(wèi)車輛的實(shí)際運(yùn)行狀況，不適用于城市環(huán)衛(wèi)車輛的工況測(cè)試.

表6 所構(gòu)建代表工況與標(biāo)準(zhǔn)工況的比較

4 結(jié)語

提出一種改進(jìn)PCA降維和k均值聚類的組合優(yōu)化方法，提高各主成分所提取到的車輛行駛特征信息量，同時(shí)對(duì)環(huán)衛(wèi)車不同行駛特征做有效分類； 以估計(jì)狀態(tài)轉(zhuǎn)移概率法取代傳統(tǒng)的隨機(jī)片段法，應(yīng)用馬爾科夫法按時(shí)間比例合成福州市環(huán)衛(wèi)清潔車的代表工況.

1) 所提方法構(gòu)建的代表工況與實(shí)際采集的行駛數(shù)據(jù)相比，各特征參數(shù)的相對(duì)誤差均在7%以下，且SAFD差異率都維持在5%以內(nèi). 說明基于所提方法開發(fā)的代表工況能很好地反映試驗(yàn)數(shù)據(jù)樣本，表征福州市環(huán)衛(wèi)清潔車輛的實(shí)際行駛特征，證明了所提工況構(gòu)建方法的有效性和合理性.

2) 現(xiàn)行對(duì)環(huán)衛(wèi)車測(cè)試采用的商用車標(biāo)準(zhǔn)測(cè)試循環(huán)C-WTVC與福州市環(huán)衛(wèi)車輛的實(shí)際行駛狀況有較大差異，并不能有效反映出城市環(huán)衛(wèi)車輛的實(shí)際運(yùn)行狀況，不適用于城市環(huán)衛(wèi)車輛的工況測(cè)試. 根據(jù)所提方法而所構(gòu)建的環(huán)衛(wèi)車代表工況能更真實(shí)表征城市環(huán)衛(wèi)車輛的行駛工況，為我國環(huán)衛(wèi)車輛行駛工況的開發(fā)提供一種有效的借鑒.