徐 婷，崔世超，劉 明，賀玉龍，石 鋒

(1. 長安大學(xué) 運(yùn)輸工程學(xué)院，陜西 西安 710064；2. 長安大學(xué) 汽車學(xué)院，陜西 西安 710064；3. 北京工業(yè)大學(xué) 城市交通學(xué)院，北京 100124; 4. 中國汽車工程研究院股份有限公司，重慶 401122)

0 引言

汽車行駛工況是描述特定交通環(huán)境下車輛行駛特征的速度-時(shí)間曲線[1]，它能夠代表被測區(qū)域內(nèi)車速變化規(guī)律[2]，對汽車節(jié)能減排有重大意義。中國國家標(biāo)準(zhǔn)化管理委員會于2019年10月下旬發(fā)布了中國標(biāo)準(zhǔn)行駛循環(huán)工況(CATC)，包括乘用車、輕型商用車、重型商用車的8條整車測試工況曲線及發(fā)動機(jī)工況曲線。此外，很多學(xué)者致力于汽車行駛工況研究:王楠楠[3]通過短行程法，采用模糊C均值聚類方法構(gòu)建了合肥市的行駛工況。張宏等[4]通過短行程法，對不同車速區(qū)間賦予權(quán)重構(gòu)建了呼和浩特市典型汽車行駛工況。姜平等[5]通過馬爾可夫鏈方法，利用最大似然估計(jì)分類方法構(gòu)建了合肥市的行駛工況。劉應(yīng)吉等[6]提出用組合主成分分析進(jìn)行汽車行駛工況的構(gòu)建。李耀華等[7]通過馬爾可夫蒙特卡洛方法構(gòu)建了西安市城市公交行駛工況。國外[8-11]主要針對城市道路進(jìn)行工況構(gòu)建。查閱文獻(xiàn)發(fā)現(xiàn)，目前的工況曲線無法涵蓋高海拔地區(qū)。

高海拔地區(qū)生態(tài)薄弱，且隨著海拔升高，汽車動力性能下降，汽車排放引起的環(huán)境問題尤為嚴(yán)重，多山的地形條件讓汽車行駛工況與城市工況有很大差異。

本文以高原山區(qū)車輛的車況為研究對象，基于青藏高原川藏線(得達(dá)鄉(xiāng)—海子山段)試驗(yàn)道路采集數(shù)據(jù)，根據(jù)汽車運(yùn)行特點(diǎn)，提出了改進(jìn)的短行程法，進(jìn)行高原山區(qū)汽車行駛工況的構(gòu)建。

1 數(shù)據(jù)采集及預(yù)處理

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

通過對我國高原地區(qū)道路線形復(fù)雜程度、地理環(huán)境特點(diǎn)對比，選取青藏高原G318川藏公路南段(得達(dá)鄉(xiāng)—海子山段)道路為試驗(yàn)道路。試驗(yàn)路段最低海拔3 610.4 m，最高海拔4 591.6 m。

本次試驗(yàn)選用多次參加道路試驗(yàn)的3名駕駛員，試驗(yàn)車輛為乘用車，車輛外廓尺寸等參數(shù)均在典型乘用車設(shè)計(jì)范圍內(nèi)，具有代表性，車輛狀況良好。通過自主行駛法駕駛試驗(yàn)車輛在試驗(yàn)道路上行駛。車輛信息如表1所示：

表1 試驗(yàn)車輛技術(shù)參數(shù)Tab.1 Technical parameters of test vehicle

試驗(yàn)車輛上裝有OBD(On Board Diagnostics)車載終端，通過OBD對車輛行駛過程中行駛車速、行駛里程、時(shí)間、發(fā)動機(jī)輸出扭矩、地理位置等信息以1 Hz的頻率進(jìn)行采集記錄。

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

在數(shù)據(jù)采集過程中可能會出現(xiàn)外界引號干擾、車輛顛簸導(dǎo)致數(shù)據(jù)采集出現(xiàn)噪聲或者傳輸錯(cuò)誤。因此，將采集到的數(shù)據(jù)進(jìn)行預(yù)處理，再對過濾后的數(shù)據(jù)進(jìn)行工況構(gòu)建。通過調(diào)查，高原地區(qū)道路車流量不大，服務(wù)水平高。結(jié)合高原地區(qū)汽車行駛特點(diǎn)，設(shè)定以下車速數(shù)據(jù)篩選條件：

(1)車速長時(shí)間(超過200 s)為0 km/h認(rèn)為數(shù)據(jù)采集裝置與試驗(yàn)車輛脫離，采集數(shù)據(jù)無效(試驗(yàn)路段為公路，擁堵情況較少發(fā)生)。

(2)加速度超過5 m/s2認(rèn)為加速度異常(平原地區(qū)汽車百公里加速時(shí)間為7 s左右，高原地區(qū)車輛動力性能會降低)。

(3)速度超過90 km/h為超速值(試驗(yàn)道路設(shè)計(jì)速度80 km/h。但根據(jù)相關(guān)研究，路段限速80 km/h，道路服務(wù)水平較高時(shí)，駕駛員平均傾向速度為90 km/h左右[12])，進(jìn)行剔除。

過濾后的部分速度與時(shí)間工況數(shù)據(jù)片段如圖1所示。

圖1 采集數(shù)據(jù)片段Fig.1 Collected data fragments

2 車速-時(shí)間工況構(gòu)建

常用的汽車行駛工況構(gòu)建方法有馬爾可夫鏈法、短行程法。

馬爾可夫鏈法：將采集數(shù)據(jù)按照統(tǒng)計(jì)學(xué)特征劃分為若干狀態(tài)，計(jì)算出各狀態(tài)的轉(zhuǎn)移概率矩陣，通過隨機(jī)數(shù)和轉(zhuǎn)移矩陣進(jìn)行行駛工況構(gòu)建[13-15]。

短行程法：將采集數(shù)據(jù)以怠速點(diǎn)劃分成若干運(yùn)動學(xué)片段，然后再選取特征值，通過主成分分析，聚類方法進(jìn)行典型行駛工況的構(gòu)建[15-17]。

高原地區(qū)有高海拔、多山的地形特點(diǎn)，且交通量較少，交通處于自由流。汽車行駛過程中為了適應(yīng)道路線形，需要不斷變化車速，導(dǎo)致怠速工況少，采用常規(guī)的短行程法不能構(gòu)建出理想的行駛工況。

首先采用馬爾可夫鏈法進(jìn)行高原山區(qū)行駛工況提取，并基于實(shí)際行駛工況特點(diǎn)，提出改進(jìn)的短行程法：將采集數(shù)據(jù)按照固定時(shí)長劃分為若干片段，計(jì)算各個(gè)片段特征值，再進(jìn)行主成分分析，通過聚類方法進(jìn)行工況構(gòu)建。

2.1 馬爾可夫鏈法

2.1.1 馬爾可夫特性的證明

“無后效應(yīng)”是馬爾可夫鏈的典型特征。下一狀態(tài)的條件概率分布僅依賴于當(dāng)前狀態(tài)，而與歷史狀態(tài)無關(guān)。假設(shè)汽車行駛過程中車速變化具有馬爾可夫特性，根據(jù)試驗(yàn)數(shù)據(jù)分別構(gòu)建不同時(shí)間間隔的子序列，理論上時(shí)間間隔越短，相關(guān)性越強(qiáng)[18]。分別采取1 s,5 s和10 s間隔進(jìn)行車速馬爾可夫特性分析。

相關(guān)性圖可以直觀反映兩個(gè)變量相互影響的強(qiáng)弱。由圖2可知，當(dāng)時(shí)間間隔為1 s時(shí)，采集數(shù)據(jù)呈現(xiàn)出強(qiáng)相關(guān)性，各點(diǎn)分布接近于一條直線；隨著時(shí)間間隔的增加，采樣點(diǎn)松散分布，變量之間的相關(guān)性逐漸減弱。為了量化第t秒車速與第t+s秒車速兩變量相關(guān)性，使用公式(1)皮爾遜相關(guān)系數(shù)進(jìn)行計(jì)算,結(jié)果如表2所示。

(1)

式中，X為第t秒的車速；Y為第t+s秒的車速；s為時(shí)間間隔。

圖2 不同時(shí)間間隔下的車速相關(guān)性圖Fig.2 Correlation diagram of vehicle speed at different time intervals

由圖2和表2可知，車速的時(shí)間間隔越小，第t秒與第t+s秒的相關(guān)性越強(qiáng)。其中，時(shí)間間隔為1 s時(shí)第t秒的車速與第t+1秒的車速相關(guān)性系數(shù)達(dá)到0.990 3，且隨著時(shí)間間隔增加，逐漸失去相關(guān)性。綜上兩點(diǎn)，汽車行駛速度具有明顯的馬爾可夫特性，在研究中能夠采用馬爾可夫鏈法構(gòu)建高原山區(qū)汽車行駛工況。

表2 不同時(shí)間間隔相關(guān)性系數(shù)Tab.2 Correlation coefficient of different time intervals

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

車速從某一時(shí)刻某一狀態(tài)變化到下一時(shí)刻某一狀態(tài)，叫做狀態(tài)轉(zhuǎn)移，所有的狀態(tài)轉(zhuǎn)移組合在一起組成狀態(tài)轉(zhuǎn)移矩陣P，如公式(2)所示。

(2)

且狀態(tài)轉(zhuǎn)移矩陣具有如下特征：

(1)矩陣中每個(gè)元素均非負(fù)；

(2)矩陣中每行和為1。

先將采集數(shù)據(jù)進(jìn)行狀態(tài)劃分，然后再根據(jù)相鄰兩點(diǎn)之間關(guān)系進(jìn)行狀態(tài)模型事件頻率的統(tǒng)計(jì)，用頻率代替概率，如公式(3)所示。

(3)

式中，Nij為狀態(tài)i轉(zhuǎn)移到狀態(tài)j發(fā)生的頻數(shù)；pij為從狀態(tài)i轉(zhuǎn)移到狀態(tài)j的概率。

對預(yù)處理后的數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，得出汽車行駛車速主要分布在40～80 km/h區(qū)間內(nèi)。將片段劃分為10個(gè)狀態(tài)(狀態(tài)1:v≤40 km/h,狀態(tài)2: 40 km/h80 km/h)。根據(jù)狀態(tài)轉(zhuǎn)移矩陣公式及狀態(tài)模型事件頻率得到狀態(tài)轉(zhuǎn)移矩陣P：

2.1.3馬爾可夫鏈法工況構(gòu)建

(4)

圖3為Δ最小的片段，選擇其最后1 s的車速作為馬爾可夫初始狀態(tài)，v初始=53.254 km/h。

圖3 Δ值最小的片段Fig.3 Segment with smallest Δ value

將確定好的馬爾可夫初始狀態(tài)按照蒙特卡洛模擬法，通過狀態(tài)轉(zhuǎn)移矩陣進(jìn)行工況構(gòu)建。因?yàn)轳R爾可夫鏈蒙特卡洛算法采用隨機(jī)數(shù)進(jìn)行工況構(gòu)建，即使初始狀態(tài)相同，也會構(gòu)建出完全不同的行駛工況。本研究構(gòu)建了10條行駛工況曲線，通過計(jì)算加速比例(Pa)、勻速比例(Pc)、減速比例(Pd)、平均速度(vm)、平均加速度(am)、速度標(biāo)準(zhǔn)偏差(vsd)、加速度標(biāo)準(zhǔn)偏差(asd)，選出與實(shí)際工況所有指標(biāo)相對誤差之和最小的一條，如圖4所示。

圖4 馬爾可夫鏈法構(gòu)建工況Fig.4 Driving cycle constructed by Markov chain method

2.2 改進(jìn)的短行程法

2.2.1 特征值選取

根據(jù)高海拔地區(qū)道路線形設(shè)計(jì)規(guī)范和汽車行駛特點(diǎn)，將預(yù)處理后的數(shù)據(jù)按照50 s的時(shí)長劃分為若干運(yùn)動學(xué)片段。選取加速比例(Pa)、勻速比例(Pc)、減速比例(Pd)、平均速度(vm)、平均加速度(am)、速度標(biāo)準(zhǔn)偏差(vsd)、加速度標(biāo)準(zhǔn)偏差(asd)為特征值。這些片段的特征值能夠詳細(xì)地描述汽車在高原地區(qū)的行駛特性，是構(gòu)建典型工況的依據(jù)[17-19]。

2.2.2主成分分析

每個(gè)特征值都能夠提供一定的車輛道路行駛狀態(tài)信息，這些變量有的會包含部分重疊信息，即變量之間具有一定的相關(guān)性，相互之間不獨(dú)立。因此，可以通過主成分分析法簡化變量個(gè)數(shù)，以便于后續(xù)處理。

主成分分析法的基本思想是嘗試組合更多的原始變量，相互獨(dú)立地構(gòu)建一些新的變量，并從實(shí)際的研究需求中選擇一些綜合變量[20]。通常，選擇80%或更多原始變量信息的線性組合可以簡化問題而不會忽略太多實(shí)際信息。

通過主成分分析法對特征值進(jìn)行處理，得到各主成分貢獻(xiàn)率和累積貢獻(xiàn)率，如表3所示。

表3 各主成分貢獻(xiàn)率和累積貢獻(xiàn)率Tab.3 Contribution rate and cumulative contribution rate of each principal component

表3中Mi為各個(gè)主成分序號，特征值的大小反映主成分方差貢獻(xiàn)大小。前3個(gè)主成分累計(jì)貢獻(xiàn)率達(dá)到85.14%，且M1,M2,M3特征值均大于1，可以全面反映片段信息，達(dá)到了主成分分析的目的。

通過主成分分析結(jié)果繪制雙標(biāo)圖如圖5所示。以提取出的前3個(gè)主成分為坐標(biāo)軸，以中心為原點(diǎn)，將各變量在前3個(gè)主成分上的載荷用矢量線段表示，各矢量夾角的余弦值是它們的相關(guān)系數(shù)，余弦值越大，相關(guān)性越強(qiáng)。各片段在前3個(gè)主成分的得分情況可在雙標(biāo)圖上顯示。

圖5 雙標(biāo)圖Fig.5 Biplot

由圖5可知，第1主成分主要反映加速比例、勻速比例、速度標(biāo)準(zhǔn)差、加速度標(biāo)準(zhǔn)差等信息；第2主成分主要反映減速比例、平均加速度等信息；第3主成分主要反映平均速度等信息。

2.2.3k-means聚類

k-means聚類是一種簡潔高效的聚類算法，通過計(jì)算不同樣本間的距離判斷它們之間的相近關(guān)系，將具有相同特性的數(shù)據(jù)對象劃到同一類別中，與相關(guān)度較小的數(shù)據(jù)做以區(qū)分。基于各片段距各聚類中心的距離進(jìn)行聚類[21-22]。

k-means聚類算法具體步驟如下：

(1)在主成分分析后的得分矩陣中選擇k個(gè)樣本點(diǎn)，將k個(gè)樣本點(diǎn)值分別作為k個(gè)初始聚類中心。

(2)依次計(jì)算各個(gè)樣本數(shù)據(jù)集中的所有數(shù)據(jù)點(diǎn)到各個(gè)初始聚類中心點(diǎn)的歐式距離，如公式(5)所示。

dij=|xi-μj|,

(5)

(3)迭代收斂

計(jì)算每個(gè)聚類的平均值，并作為新的中心點(diǎn)，重復(fù)上述過程，直到這k個(gè)中線點(diǎn)收斂了，輸出劃分結(jié)果。

圖6是聚類后的結(jié)果，將 138個(gè)樣本聚為3簇。

圖6 聚類結(jié)果Fig.6 Clustering result

第1簇有53個(gè)片段，第2簇有42個(gè)片段，第3簇有43個(gè)片段。將各運(yùn)動學(xué)片段按照類別重新編號，并根據(jù)典型特征值進(jìn)行分析，如圖7所示。第1簇加速工況占比最大，為加速路段；第2簇勻速工況占比最大，為勻速路段；第3簇的3種工況占比差別不大，分析為試驗(yàn)車輛行駛路段線形復(fù)雜，駕駛員需要不斷變化車速以適應(yīng)道路變化。

圖7 聚類結(jié)果分析Fig.7 Cluster result analysis

根據(jù)各片段到聚類中心的距離大小，按照運(yùn)動學(xué)片段原始序號，第1簇選擇片段63，36，50，68，127，23；第2簇選擇片段129，77，102，108，97；第3簇選擇片段74，31，76，25，41。將所選片段按照時(shí)間順序進(jìn)行工況構(gòu)建。

在工況構(gòu)建過程中，可能會出現(xiàn)車速較大的點(diǎn)與車速較小的點(diǎn)相連接，導(dǎo)致加減速強(qiáng)度過大或加速度數(shù)值異常。對這些點(diǎn)采取中位值平均濾波法進(jìn)行濾波：以加速度異常值處為中心，選擇14個(gè)前后相鄰點(diǎn)，去掉這些點(diǎn)中的速度最大值和速度最小值，用剩余點(diǎn)的平均速度代替加速度數(shù)值異常處的點(diǎn)，并控制連續(xù)加速時(shí)間。構(gòu)建工況如圖8所示。

圖8 改進(jìn)的短行程法構(gòu)建工況Fig.8 Driving cycle constructed by improved short stroke method

3 對比分析

將馬爾可夫鏈法構(gòu)建高原山區(qū)典型行駛工況1和改進(jìn)的短行程法構(gòu)建的高原山區(qū)典型行駛工況2與實(shí)際工況進(jìn)行對比，如表4、表5所示：

表4 構(gòu)建行駛工況1與實(shí)際工況對比Tab.4 Comparison between constructed driving cycle 1 and actual driving cycle

表5 構(gòu)建行駛工況2與實(shí)際工況對比Tab.5 Comparison between constructed driving cycle 2 and actual driving cycle

由表4、表5可知，馬爾可夫鏈法構(gòu)建工況的特征值與實(shí)際工況偏差較大，其中，勻速工況誤差最大。分析原因?yàn)闋顟B(tài)轉(zhuǎn)移矩陣主對角線數(shù)值最大，即通過馬爾可夫鏈法構(gòu)建的行駛工況偏向于車輛穩(wěn)定行駛工況，這一特點(diǎn)與城市工況較為符合；改進(jìn)的短行程法構(gòu)建工況車速變化復(fù)雜，特征值與實(shí)際工況特征值誤差較小，能夠較好地?cái)M合高原山區(qū)汽車行駛工況。

將馬爾可夫鏈法與改進(jìn)的短行程法構(gòu)建的工況進(jìn)行對比，如圖9所示。馬爾可夫鏈方法構(gòu)建的高原地區(qū)行駛工況較為平穩(wěn)，車速變化幅度及頻率較?。桓倪M(jìn)的短行程法構(gòu)建的行駛工況車速變化頻繁，與高海拔地區(qū)多山的道路環(huán)境特征相符。

圖9 構(gòu)建工況對比Fig.9 Comparison of constructed driving cycles

4 結(jié)論

本研究以川藏線道路采集工況數(shù)據(jù)為研究對象提出了構(gòu)建高原山區(qū)乘用車行駛工況的方法。

(1)通過不同時(shí)間間隔證明汽車行駛工況具有馬爾可夫特性。根據(jù)馬爾可夫蒙特卡洛法構(gòu)建的高原山區(qū)行駛工況較為平穩(wěn)，以勻速工況為主，車速變化幅度及頻率較小。

(2)用改進(jìn)的短行程法按照50 s的固定時(shí)長將汽車行駛工況劃分為138個(gè)運(yùn)動學(xué)片段，通過主成分分析，k-means聚類構(gòu)建行駛工況，并通過中位值平均濾波法進(jìn)行異常點(diǎn)的濾波。改進(jìn)的短行程法構(gòu)建的行駛工況加速，減速的波動變化較為頻繁，

(3)通過比較指標(biāo)(加速比例、勻速比例、減速比例、平均速度、平均加速度、速度標(biāo)準(zhǔn)偏差、加速度標(biāo)準(zhǔn)偏差)，將兩種方法構(gòu)建的行駛工況與實(shí)際工況對比，發(fā)現(xiàn)用改進(jìn)的短行程法構(gòu)建的行駛工況誤差較小，可以反映高原山區(qū)車輛行駛特點(diǎn)，這為研究行駛工況提供了新思路。