汪雯琦，高廣闊，王子鑒，梁易鑫

(上海理工大學(xué) 管理學(xué)院，上海 200093)

0 引言

目前，我國汽車運行工況的相關(guān)標準主要是在上個世紀引進的德國和歐洲標準。這些標準與我國汽車運行工況的差異很大。近年來，隨著汽車保有量的快速增長和城市化率的不斷增加，我國道路交通狀況變化很大,且不同地區(qū)汽車運行工況差異很大，目前缺乏統(tǒng)一適用的標準。當前，制訂出適合我國各地國情的油耗標準和排放規(guī)范是當務(wù)之急。

運動學(xué)片段是指汽車從怠速狀態(tài)開始至相鄰的下一個怠速狀態(tài)開始之間的車速變化狀況的連續(xù)過程，如圖1所示。本研究設(shè)計了一種合理的方法，將預(yù)處理后的數(shù)據(jù)劃分為多個運動學(xué)片段，并給出通過數(shù)據(jù)文件處理計算得到的運動學(xué)片段的最終數(shù)量[1]。

圖1 運動學(xué)片段的定義Fig.1 Definition of kinematic segment

為了在國際市場上占據(jù)新一代能源汽車性能的領(lǐng)先地位，歐洲、日本、美國加快了汽車運行工況的研發(fā)[2-7]。近年來，國內(nèi)學(xué)者也開始對車輛的運行狀況進行深入研究[8-9]。多數(shù)學(xué)者采用結(jié)合主成分分析(PCA)和K-means聚類法來進行汽車運行工況的擬合，但傳統(tǒng)K-means聚類算法的初始聚類中心會導(dǎo)致聚類結(jié)果下降，陷入局部最大值。胡宸等[10]、胡志遠等[11]使用PCA、K-means聚類分別構(gòu)建了哈爾濱和上海的車輛運行工況，但均未在初始聚類中心的基礎(chǔ)上尋優(yōu)。石琴等[12]改進了初始聚類中心，并使用神經(jīng)網(wǎng)絡(luò)算法構(gòu)建了汽車運行工況，但神經(jīng)網(wǎng)絡(luò)需要設(shè)置初始權(quán)值，導(dǎo)致聚類質(zhì)量下降。

本研究將在公共交通系統(tǒng)、城市發(fā)展及汽車排放標準方面起到推動作用，其意義主要體現(xiàn)在：

(1)有助于公共交通系統(tǒng)的建設(shè)與發(fā)展，減少汽車怠速時間和緩解城市擁堵，從而提高出行效率。

(2)能規(guī)范城市汽車市場準入標準，解決汽車采購中的亂象，有助于有關(guān)部門掌握車輛的實際使用情況，提升政府決策能力。(3)車輛運行工況的變化對汽車燃料消耗和污染物排放有重要影響，在節(jié)能減排和環(huán)境保護以及在國際市場上搶占新能源汽車性能的領(lǐng)先地位有重要作用。

1 模型

1.1 原始數(shù)據(jù)

用行車記錄儀采集記錄上海市某輛輕型汽車在實際普通道路上行駛的數(shù)據(jù)，采集3次包括3個不同時間段的信息，采樣頻率為1 Hz。先判斷各個指標所代表的含義?！八俣取碧峁┝塑囕v的速度信息，由于采樣時間的間隔為1 s，因此也可計算出車輛行駛加速度?！叭S加速度”提供了車輛3個軸向的加速度信息，對后續(xù)速度與加速度的計算與分析意義不大。“經(jīng)緯度信息”提供了計算行駛距離的可能性，可根據(jù)具體公式計算汽車行駛距離?！鞍l(fā)動機轉(zhuǎn)速”和“節(jié)氣門開度”[13]有助于判斷汽車是否處于熄火狀態(tài)。其他變量也會為研究汽車運行狀況提供幫助。車輛主要結(jié)構(gòu)、配置和性能參數(shù)情況見表1。

表1 試驗車輛情況Tab.1 Conditions of experimental vehicles

1.2 數(shù)據(jù)預(yù)處理

1.2.1缺失數(shù)據(jù)填充

由于高層建筑遮擋或經(jīng)過隧道等，行車記錄儀所記錄的信息會出現(xiàn)GPS信號過弱或丟失，造成數(shù)據(jù)記錄不連續(xù)。首先判斷數(shù)據(jù)丟失時汽車可能的行駛狀態(tài)，以便采用不同方法對缺失數(shù)據(jù)進行填充。數(shù)據(jù)丟失時汽車可能處于以下3類行駛狀態(tài)。

第1類：缺失數(shù)據(jù)時間前后的行駛速度均大于10 km/h，認為在信號丟失的這段時間內(nèi)汽車處于正常行駛狀態(tài)，可能由于高層建筑覆蓋或過隧道等原因造成了信號丟失。根據(jù)信號丟失前后的GPS速度，利用線性插值法對丟失的速度進行填充。

第2類：缺失數(shù)據(jù)時間前一個狀態(tài)GPS速度為0且油門踏板開度也為0，則認為發(fā)動機未處于運行狀態(tài)，故車輛處于熄火狀態(tài)，GPS可能也處于待機狀態(tài)，未檢測到數(shù)據(jù)信息，視為汽車怠速進行處理，將速度填充為0。對于間隔時長超過180 s的數(shù)據(jù)，按照怠速長度180 s進行處理。

第3類：缺失數(shù)據(jù)時間前一個狀態(tài)GPS車速為0或小于10 km/h，而加速踏板位移不為0，可能由于進隧道長時間堵車等原因造成，認為這一階段處于怠速，按照車速為0進行數(shù)據(jù)填充。

對于時間不連續(xù)問題的處理，將后一個數(shù)據(jù)的秒數(shù)減去前一個數(shù)據(jù)的秒數(shù)，若差值等于1，則繼續(xù)進行；若差值不等于1，需要將其標記為0并返回該值所在的位置，進而提取缺失時間數(shù)據(jù)前后的速度和加速踏板位移值，以判斷車輛在信號丟失時間段的可能行駛狀態(tài)。對于不同的行駛狀態(tài)，采用上述方法進行相應(yīng)處理。

利用Matlab進行計算，發(fā)現(xiàn)724處時間不連續(xù)的情況，符合上述第1類情況的時間間隔有429處，符合第2類的有189處，符合第3類的有106處，分別進行相應(yīng)的處理與填充。

1.2.2平滑速度異常狀態(tài)數(shù)據(jù)

行車記錄儀記錄的速度單位為km/h，加速度的單位為m/s2，所以先對速度進行單位轉(zhuǎn)化，使其單位變?yōu)閙/s，計算公式為：

i=1, 2，…，k-1,

(1)

式中，ai,i+1為第is到第(i+1)s的加速度；vi+1，vi分別為第is和第(i+1)s的速度；ti+1，ti分別為第is和第(i+1)s的時刻；k為運動學(xué)片段的總時間，也是總速度的個數(shù)。

根據(jù)式(1)算出各個時間點上的加速度值。通常認為，一般情況下，普通轎車0～100 km/h的加速時間大于7 s，緊急制動最大減速度為7.5～8 m/s2。所以可認為當加速度大于3.97 m/s2或減速度小于-8 m/s2時為異常值。采用移動平均濾波器Smooth函數(shù)進行平滑處理以消除這些異常值，默認窗寬為5。

1.2.3長期停車數(shù)據(jù)處理

對長期停車(如停車不熄火等候人、停車熄火但采集設(shè)備仍在工作等)所采集的異常數(shù)據(jù)，首先明確長時間停車的類型，分為是/否熄火兩種情況進行考慮。如果發(fā)動機轉(zhuǎn)速為0，表示發(fā)動機已處于熄火狀態(tài)，將該數(shù)據(jù)做刪除處理；如果發(fā)動機轉(zhuǎn)速不為0，但車速為0，按怠速處理。

1.2.4怠速情況分析

將長時間堵車、斷續(xù)低速行駛狀況(車速小于10 km/h)按照怠速工況處理。首先找出車輛行駛速度長時間小于10 km/h的連續(xù)時段，選擇車速小于10 km/h的時間點進行劃分。根據(jù)其連續(xù)性的判斷與選擇，合并連續(xù)的時間點，劃分為若干個連續(xù)片段。如果某片段的持續(xù)時間長度大于180 s，則將這段時間視為怠速工況狀態(tài)；如果持續(xù)長度未達到180 s，為了保證數(shù)據(jù)的準確性，則不作處理而保留原始速度信息。

1.2.5怠速數(shù)據(jù)處理

找到怠速時長超過180 s的片段，利用滑動窗口法進行解決。首先定義一個長度為180 s的窗口，并用這一窗口按照從上到下的順序去截取速度為0的數(shù)據(jù)。當截取的這段數(shù)據(jù)中車速均為0時，需刪除1個數(shù)據(jù)再繼續(xù)截取，直至遍及所有數(shù)據(jù)。為了提高滑動速度，若截取到的1段數(shù)據(jù)的最后一個車速值不為0，就直接跳過這一片段，繼續(xù)進行滑動。最后保證數(shù)據(jù)集中沒有出現(xiàn)超過180 s的怠速工況。

數(shù)據(jù)的預(yù)處理均利用Matlab軟件編程解決。在補齊缺失數(shù)據(jù)、處理或剔除異常值數(shù)據(jù)后，行車記錄儀在3個不同時間段采集的數(shù)據(jù)最終記錄數(shù)分別為205 990，173 368，180 786。預(yù)處理前后的數(shù)據(jù)對比如圖2所示。

圖2 三個片段運動學(xué)數(shù)據(jù)預(yù)處理前后對比Fig.2 Comparison of kinematic data of 3 segments before and after preprocessing

2 運動學(xué)片段的提取

2.1 運動學(xué)片段相關(guān)理論

在構(gòu)建道路運動工況的過程中存在著多種方法，其中應(yīng)用最為廣泛的是將全部試驗數(shù)據(jù)按照某種既定原則劃分為若干個片段，并對所有小片段進行研究。選擇適當?shù)姆椒?gòu)建運行工況，而其中這些片段即為運動學(xué)片段(從怠速狀態(tài)開始至相鄰的下一個怠速狀態(tài)開始之間車速變化狀況的連續(xù)過程)。用此方法構(gòu)建的運行工況更符合實際。

將上述運動學(xué)片段中的幾種運動學(xué)狀態(tài)進行劃分如下。

(1)怠速工況：發(fā)動機工作，且車輛速度V為0的狀態(tài)。

(2)勻速工況：車輛加速度a小于0.10 m/s2并大于-0.10 m/s2，且V不為0的狀態(tài)。

(3)加速工況：車輛加速度a大于0.10 m/s2，且V不為0的狀態(tài)。

(4)減速工況：車輛減速度a小于-0.10 m/s2，且V不為0的狀態(tài)。

2.2 運動學(xué)片段的劃分

為構(gòu)建運行工況，需要對運動片段進行劃分。劃分標準是：從怠速狀態(tài)開始至相鄰的下一個怠速狀態(tài)開始之間車速變化狀況的連續(xù)過程。由于汽車行駛狀況數(shù)據(jù)的試驗數(shù)據(jù)量較大，且城市道路行駛過程中啟停較為頻繁，需要進行行程劃分的數(shù)據(jù)數(shù)量較大，必須利用軟件進行處理，選用Matlab軟件輔助處理，分別對3個行駛數(shù)據(jù)文件的數(shù)據(jù)進行切分，形成若干運動學(xué)片段，分別劃分成1 716，2 184，1 681個運動學(xué)片段。

2.3 運動學(xué)片段的特征提取與計算

表2 特征參數(shù)Tab.2 Characteristic parameters

由上述選取的特征參數(shù)中，發(fā)現(xiàn)加速度值較多為5個，因其不能直接通過觀測得到，故需要根據(jù)GPS速度，使用相關(guān)公式進行計算。

13個特征值的計算公式如下。

加速時間百分比：

Pa=Ta÷T，

(2)

式中，Ta為運動學(xué)片段中加速度大于0.1 m/s2的時間數(shù)；T為運動學(xué)片段的總時間數(shù)。

減速時間百分比：

Pd=Td÷T，

(3)

式中，Td為運動學(xué)片段中加速度小于-0.1 m/s2的時間數(shù)。

怠速時間百分比：

P0=T0÷T，

(4)

式中T0為運動學(xué)片段中的怠速的時間數(shù)。

勻速時間百分比：

Pc=(T-Ta-Td-T0)/T。

(5)

最大速度：

vmax=max{vi,i=1, 2，…,k}。

(6)

平均速度：

(7)

式中,S為每個運動學(xué)片段所代表的汽車行駛距離。由于該試驗的采樣頻數(shù)設(shè)置為1 Hz，即每兩個點的時間間隔為1 s，所以S為該片段的速度之和，計算式為：

(8)

平均行駛速度：

Vm=S/T。

(9)

速度標準差：

(10)

式中,Vi為汽車行駛數(shù)據(jù)中按時間順序的第i個速度;k為行駛數(shù)據(jù)的總條數(shù)。

加速度相關(guān)特征值：

(11)

通過編程對上述特征參數(shù)進行了計算，并為后續(xù)運動學(xué)片段的篩選及汽車運行工況的構(gòu)建提供了重要信息。分別對3個文件中的數(shù)據(jù)進行了特征值的提取與計算，每份汽車行駛數(shù)據(jù)文件中分別選取了1張運動學(xué)片段進行前后對比。

圖3 三個行駛數(shù)據(jù)文件中預(yù)處理前后的汽車行駛工況對比Fig.3 Comparison of vehicle driving conditions in 3 driving data files before and after preprocessing

由圖3可明確地識別該車在特定時間段的行駛狀態(tài)。接下來對上述所有分割出來的運動學(xué)片段進行有效篩選，最終選取最具代表性且符合真實道路行駛的片段以構(gòu)建較為準確的汽車運行工況。

3 汽車運行工況的構(gòu)建

3.1 構(gòu)建流程

因特征參數(shù)較多，選用主成分分析模型對特征參數(shù)進行降維處理。將13個特征參數(shù)變量轉(zhuǎn)化為數(shù)量較少的幾個主成分，接著以這幾個主成分作為變量，運用K-means聚類模型對運動學(xué)片段進行分類，主要構(gòu)建流程如圖4所示。

圖4 汽車運行工況構(gòu)建流程Fig.4 Flowchart of constructing vehicle operation condition

3.2 基于主成分分析的數(shù)據(jù)降維

按照主成分分析步驟[14]，對標準化后的特征值進行分析。計算13個特征參數(shù)的方差貢獻率和累積貢獻率，得到結(jié)果列在表3中。

選取累積貢獻率達到80%作為基準，達到的為第4項的84.09%，且前4項主成分的特征值均大于1，基本涵蓋全部信息，所以提取前4項主成分表征原始特征參數(shù)。計算因子載荷矩陣，如表4所示。

主成分1描述了平均速度、最大速度、平均行駛速度、速度標準差、最大加速度、怠速時間百分比和等速時間百分比。主成分2描述了平均減速度、平均加速度和加速度標準差。主成分3描述了最大減速度和減速時間百分比。主成分4描述了加速時間百分比。所以主成分1～4就是降維后的主成分。

根據(jù)：

(12)

式中，Mi為主成分得分序列；bij為對應(yīng)的成分系數(shù)；Ym為標準化后的特征值。

從而得到第1～第4主成分的得分。經(jīng)計算后式(12)變?yōu)椋?/p>

(13)

主成分載荷(無量綱)見表5。

所以可得4個主成分為:

M1=0.37Y1+0.36Y2+0.35Y3+0.34Y4+0.31Y5+

0.30Y6+0.29Y7+0.01Y8-0.02Y9-0.06Y10-

0.04Y11-0.03Y12+0.1Y13，

(14)

M2=-0.02Y1-0.13Y2-0.07Y3+0.26Y4+0.31Y5+

0.30Y6+0.34Y7+0.51Y8+0.45Y9-0.44Y10-

0.08Y11-0.16Y12-0.03Y13，

(15)

M3=-0.09Y1+0.15Y2-0.07Y3-0.09Y4-0.14Y5-

0.10Y6-0.22Y7+0.01Y8-0.32Y9+0.17Y10+

0.80Y11+0.78Y12+0.02Y13，

(16)

M4=0.15Y1+0.10Y2+0.05Y3+0.18Y4+0.00Y5+

0.33Y6+0.15Y7-0.14Y8+0.03Y9-0.57Y10+

0.09Y11-0.13Y12+0.94Y13。

(17)

由上列4個表達式，通過計算變量Y1到Y(jié)13得到M1,M2,M3,M4的值。

然后進入SPSS的標準化處理，其中Y1為怠速

表3 運動學(xué)片段的主成分貢獻率矩陣Tab.3 Principal component contribution matrix of kinematic segments

時間百分比，Y2為加速時間百分比，Y3為減速時間百分比，以此類推計算出主成分M1,M2,M3,M4的結(jié)果。

3.3 K-means聚類模型

降維處理后，運用SPSS軟件以及Jupyter Notebook工具進行降維得到4個主成分分別進行分類處理。具體通過聚類分析，即經(jīng)降維處理，對運動學(xué)片段數(shù)據(jù)按各自特性和一定的規(guī)則分成3類，找到了運動學(xué)片段樣本數(shù)據(jù)中的相似和區(qū)別，大量減少了研究的數(shù)量，在車輛研究領(lǐng)域有重要意義。

3.3.1分類數(shù)量的確定

(1)初步設(shè)定

據(jù)文獻[15-16]，汽車運動狀態(tài)通常被分為3種：走走停停、低速行駛、高速行駛。因此，在SPSS運行時將分類數(shù)初步設(shè)定為3。

(2)輪廓系數(shù)法

在具體決定分類法個數(shù)時，采用了輪廓系數(shù)進行聚類評估。其計算公式為：

表4 運動學(xué)片段的因子載荷矩陣Tab.4 Factor load matrix of kinematic segment

表5 主成分載荷Tab.5 Principal component loads

(18)

式中，a(i)為點i到所有它所屬的簇中其他點距離的平均值;b(i)為點i到與它相鄰最近的一簇內(nèi)所有點平均距離的最小值。

通過Jupyter Notebook編程，得出了圖5所示的初始化群集數(shù)折線圖。

圖5 初始化群集數(shù)折線圖Fig.5 Line chart of initialization cluster number

由圖5可知，當分類數(shù)為3時，s(i)更接近1，此時的聚類合理。這再次驗證了分類為3是正確合理的。

3.3.2直觀驗證

根據(jù)Jupyter Notebook動態(tài)編程進行2類、3類分類嘗試。經(jīng)比較，分為3類合理。

3.3.3運動學(xué)片段的聚類分析

結(jié)合以上3種方法，最終確定分類數(shù)為3時最為合適。經(jīng)過100次迭代，主成分被成功地分成3類，聚類中心[17]見表6。

表6 聚類中心Tab.6 Cluster centers

運動學(xué)片段聚類結(jié)果見圖6。由于數(shù)據(jù)眾多，為方便觀測，作隨機失活處理。第1類三角形圖標含有3 112個運動學(xué)片段數(shù)，占總時間比例的55.76%。第2類正方形圖標含有2 282個運動學(xué)片段數(shù)，占總時間比例的40.89%。第3類圓形圖標含有187個運動學(xué)片段數(shù)，占總時間比例的3.35%。如圖6所示，從左至右圖標設(shè)置不同，分別為正方形圖標、三角形圖標、圓形圖標。

圖6 該輛車所有數(shù)據(jù)運動學(xué)片段的聚類Fig.6 Clustering of all data kinematic segments of vehicle

根據(jù)Jupyter Notebook編程將3類數(shù)據(jù)進行分析。

分析發(fā)現(xiàn)，第3類圓形圖標：走走停停；第2類正方形圖標：高速行駛；第1類三角形圖標：低速行駛。這與現(xiàn)實情況基本吻合，再次驗證了選擇3種分類情況是正確的。綜上所述，對于此輛汽車的運行工況數(shù)據(jù)進行分析，得出了車輛運行被分為3類，且相互之間具有明顯的差異。

表7 部分片段分類結(jié)果Tab.7 Partial fragment classification result

3.4 構(gòu)建汽車運行工況

與數(shù)據(jù)預(yù)處理后的原始數(shù)據(jù)進行匹配，得到分類后的數(shù)據(jù)特征值，然后對第1類～第3類的數(shù)據(jù)進行總體描述，得到分類后的總體數(shù)據(jù)的特征值，見表8。

接著，運用Matlab軟件編程計算出與總體特征矩陣的相關(guān)系數(shù)，匯總后得出表9。

表9表示隨機選出的每類與總數(shù)特征值相關(guān)性較高的計算結(jié)果，沒有按照降序排列，且在生成的結(jié)果中所有的相關(guān)性系數(shù)均大于0.5，這說明聚類分析較為正確，再次驗證了分為3類是合理正確的。

3.5 繪制汽車運行工況

由于第1類～第3類運動學(xué)片段的占比分別為55.76%，40.89%，3.35%。故選取了15個運動學(xué)片段，其中第1類情況中選取了6個運動學(xué)片段，第2類情況中選取了4個運動學(xué)片段，第3類情況中選取了5個運動學(xué)片段[18]。之所以選擇較多片段，是因為這些運動學(xué)片段時間較短。選取的15個汽車運動學(xué)片段的情況見表10。

表8 總體分類后的數(shù)據(jù)特征值Tab.8 Characteristic values of data after overall classification

表9 與總體特征值相關(guān)性較高的片段Tab.9 Fragments with high correlation with the overall eigenvalues

表10 選取的運動學(xué)片段的相關(guān)性系數(shù)Tab.10 correlation coefficient of selected kinematic segments

選取的片段與總數(shù)特征值矩陣的相關(guān)性系數(shù)均為0.9以上，可以代表總體特征，相關(guān)性很強，以此繪制汽車運行工況可以起到很好的說明效果。

使用SPSS軟件進行散點圖的繪制，由于時間的間隔為1 s，將各個散點用直線連接起來就是車輛運行工況。圖13很好地反映了汽車走走停停、加速、減速、勻速、怠速的過程，可以有效地反映車輛運行工況。

通過計算誤差率，形成擬合誤差分析表(表11)，擬合運行工況的絕大多數(shù)特征參數(shù)都與原始數(shù)據(jù)的參數(shù)十分接近，在10%以內(nèi)，只有個別參數(shù)的誤差大于10%，但都小于20%?？傮w說明，構(gòu)建出的汽車運行工況與實際工況很相近。

圖7 車輛運行工況Fig.7 vehicle operation condition

表11 擬合誤差分析表Tab.11 Analysis of fitting error

4 結(jié)論

本研究選取不同的方法對采集到的車輛的3段行駛數(shù)據(jù)進行預(yù)處理并進行運動學(xué)片段劃分，從中選取13個特征參數(shù)以代表運動學(xué)片段的信息，并建立PCA和K-means模型對特征參數(shù)進行降維，將數(shù)據(jù)進行分組聚類。不斷提取離聚類中心最近的(即相似度最高的)運動學(xué)片段以實現(xiàn)汽車運行工況的構(gòu)建。最后將擬合工況的各個特征參數(shù)值與原始數(shù)據(jù)的工況參數(shù)值進行對比，誤差較小，工況擬合效果很好。進而得到以下結(jié)論：

(1)本研究的汽車運行工況構(gòu)建方法可以有效地提取汽車行駛信息，有助于構(gòu)建汽車運行工況，得到的汽車工況誤差很小，優(yōu)于學(xué)界大部分方法，對建立我國符合特定城市的交通狀況及實際運乘條件的汽車性能與燃油排放控制的研究提供了一定的參考。

(2)從最終擬合的運行工況不難看出，我國運行工況的怠速階段時間明顯高于歐美國家的工況，歐洲工況具有較少的怠速運行時間[19]。所以若用歐美國家[20-21]的汽車油耗量與污染物排放量的標準來制訂我國的污染控制策略并不符合實際國情。建議結(jié)合自身情況，建立符合我國道路的汽車運行工況。

但由于時間的關(guān)系，沒有嘗試使用其他模型去擬合運行工況[22]，后續(xù)可以嘗試馬爾可夫法、小波變化法、SOM神經(jīng)網(wǎng)絡(luò)模型等方法去改進現(xiàn)有的算法，這樣有可能更有效地構(gòu)建汽車運行工況，而目前采用的方法也可較為準確地構(gòu)建汽車運行工況。