徐宗煌, 林 瑤
(1.福州理工學(xué)院 應(yīng)用科學(xué)與工程學(xué)院,福建 福州 350506; 2.福州大學(xué) 環(huán)境與資源學(xué)院,福建 福州 350108)
本文的問題和數(shù)據(jù)均源于2019年“華為杯”第十六屆中國研究生數(shù)學(xué)建模競賽D題.汽車行駛工況(Driving Cycle)又稱車輛測試循環(huán)工況,是指汽車在行駛過程中具有代表性的速度-時(shí)間關(guān)系曲線.由于該曲線能夠反映汽車道路行駛的運(yùn)動(dòng)學(xué)特征[1],是汽車行業(yè)的一項(xiàng)重要的、共性基礎(chǔ)技術(shù),也是汽車各項(xiàng)性能指標(biāo)標(biāo)定優(yōu)化時(shí)的主要基準(zhǔn).因此,基于國內(nèi)汽車行駛的工況特點(diǎn)和規(guī)律,構(gòu)建出適合我國國情的汽車行駛工況曲線具有的重要的經(jīng)濟(jì)意義和實(shí)用價(jià)值[2].
目前,許多汽車發(fā)達(dá)國家均已開發(fā)了適合實(shí)際情況的汽車行駛工況,早在20世紀(jì)50年代,美國就制定了世界上第一個(gè)汽車行駛工況標(biāo)準(zhǔn),之后歐洲部分國家、日本等,陸續(xù)頒布了具有代表性的典型汽車行駛工況[3].現(xiàn)階段,國內(nèi)外構(gòu)建汽車行駛工況的方法主要有短行程法[4—5]、小波變換法[6]、聚類分析法[7—9]及馬爾可夫分析法[10].本文基于城市汽車實(shí)際道路行駛采集的數(shù)據(jù),通過對時(shí)間不連續(xù)數(shù)據(jù)、尖點(diǎn)數(shù)據(jù)、毛刺噪聲數(shù)據(jù)以及怠速異常數(shù)據(jù)進(jìn)行一系列預(yù)處理,結(jié)合主成分和K-means聚類分析法,對汽車行駛工況進(jìn)行構(gòu)建并分析誤差.
原始數(shù)據(jù)中不良數(shù)據(jù)的存在,會(huì)嚴(yán)重影響后續(xù)工作的結(jié)果,而數(shù)據(jù)預(yù)處理的方法有很多種,對不同類型的數(shù)據(jù),其處理方法也不盡相同.因此,設(shè)計(jì)合理的方法對不良數(shù)據(jù)依次進(jìn)行預(yù)處理顯得尤為重要,從而使得預(yù)處理后的數(shù)據(jù)能夠反映汽車的真實(shí)行駛情況,為后續(xù)運(yùn)動(dòng)學(xué)片段的提取以及汽車行駛工況的構(gòu)建做準(zhǔn)備.汽車道路行駛數(shù)據(jù)預(yù)處理的流程圖見圖1.
圖1 數(shù)據(jù)預(yù)處理流程圖
1.2.1時(shí)間不連續(xù)數(shù)據(jù)處理 對于時(shí)間不連續(xù)數(shù)據(jù)的處理,關(guān)鍵在于如何定義短時(shí)間丟失和長時(shí)間丟失兩種情況的臨界時(shí)間,考慮一般汽車經(jīng)過隧道的情況,不妨定義臨界時(shí)間為我國公路隧道平均長度與城市汽車平均行駛速度之比,即
(1)
式中:τ為臨界時(shí)間;L為我國公路隧道平均長度,由我國公路隧道長度L′與隧道數(shù)量N的比值確定,即
(2)
查閱資料有L′=14 039.7 km,N=15 181;V為城市汽車平均行駛速度,查閱資料可知V=40 km/h.
由 (1)式和(2) 式,可得
因此,不妨取臨界時(shí)間為τ=85 s,則認(rèn)為τ>85 s屬于“車速凍滯帶”,不進(jìn)行插值補(bǔ)全;而2≤τ≤85為數(shù)據(jù)短時(shí)間丟失,采用線性插值法進(jìn)行平滑處理,即
(3)
式中:(t1,v1)和(t2,v2)為短時(shí)間缺失的兩個(gè)相鄰時(shí)間GPS車速對應(yīng)的坐標(biāo)點(diǎn);t為缺失的時(shí)刻;v為缺失的時(shí)刻要插入的速度值.
分別對文件1~文件3的原始數(shù)據(jù)進(jìn)行時(shí)間不連續(xù)數(shù)據(jù)處理,可得到3個(gè)文件經(jīng)過處理后增加的記錄數(shù)(表1).
表1 時(shí)間不連續(xù)數(shù)據(jù)處理過程記錄
1.2.2尖點(diǎn)數(shù)據(jù)處理 由已知普通轎車一般情況下0~100 km/h的加速時(shí)間大于7 s,可知其最大加速度
而緊急剎車最大減速度為7.5~8 m/s2,不妨令最大減速度為
dmax=7.5 m/s2.
隨后,經(jīng)多次采用線性插值法進(jìn)行平滑處理,使得汽車加、減速度處于下式的正常范圍:
(4)
由于采集數(shù)據(jù)的采樣頻率為1 Hz,則線性插值公式為
(5)
式中:vxt為第t個(gè)時(shí)刻的GPS車速數(shù)據(jù)修正后的值;vt-1和vt+1分別為第(t-1)和第(t+1)個(gè)時(shí)刻的GPS車速數(shù)據(jù).
分別對文件1~文件3的原始數(shù)據(jù)進(jìn)行尖點(diǎn)數(shù)據(jù)處理,可得到3個(gè)文件處理的過程記錄(表2).
表2 尖點(diǎn)數(shù)據(jù)處理過程記錄
值得注意的是:文件3中會(huì)出現(xiàn)異常的GPS速度,而導(dǎo)致汽車加、減速度異常的情況,即GPS速度已超出實(shí)際情況的車速數(shù)據(jù).本文以高速公路上行駛最高限速120 km/h為標(biāo)準(zhǔn),舍棄GPS車速大于120 km/h的數(shù)據(jù).由此可得文件3中刪除的數(shù)據(jù)記錄數(shù)為345.
1.2.3毛刺噪聲數(shù)據(jù)處理 由于可按怠速情況處理毛刺噪聲數(shù)據(jù),因此需要對GPS車速數(shù)據(jù)進(jìn)行修正,以反映汽車的真實(shí)行駛工況.查閱資料[11]可知,最終的怠速修正上限為v0=10 km/h,線性修正上限為15 km/h,即可得修正公式
(6)
通過(6)式對原始數(shù)據(jù)進(jìn)行修正,可得毛刺噪聲數(shù)據(jù)處理過程記錄(表3).
表3 尖點(diǎn)數(shù)據(jù)處理過程記錄
1.2.4怠速異常數(shù)據(jù)處理 對于怠速異常數(shù)據(jù)情況,保留180 s的數(shù)據(jù),直接將怠速時(shí)間大于180 s的數(shù)據(jù)刪除即可.經(jīng)過統(tǒng)計(jì),可得怠速異常數(shù)據(jù)處理的過程記錄(表4).
經(jīng)過上述對原始數(shù)據(jù)預(yù)處理,可得到文件1~文件3各文件數(shù)據(jù)經(jīng)預(yù)處理后的記錄數(shù)(表5).
表5 數(shù)據(jù)預(yù)處理后的記錄數(shù)
按照圖2的研究思路,首先,基于預(yù)處理后的數(shù)據(jù),對運(yùn)動(dòng)學(xué)片段進(jìn)行劃分[12],考慮到所選取汽車運(yùn)動(dòng)特征評(píng)估體系指標(biāo)合理性,同時(shí)引入14個(gè)相關(guān)的運(yùn)動(dòng)特征指標(biāo);其次,由于選取的指標(biāo)較多,可能某些指標(biāo)之間相互存在信息重疊的現(xiàn)象,因此采用主成分分析對運(yùn)動(dòng)特征指標(biāo)進(jìn)行降維處理,通過將較多的指標(biāo)轉(zhuǎn)化為幾個(gè)主要成分,可以大大降低分析和計(jì)算的難度;再次,選取具有代表性的5個(gè)主成分,基于K-means聚類分析算法將運(yùn)動(dòng)學(xué)片段,按照擁堵路況、一般擁堵路況和通暢路況分為3類,再根據(jù)聚類結(jié)果選取最優(yōu)工況片段,同時(shí)以總時(shí)間1 200~1 300 s為要求進(jìn)行組合,從而得到最終的工況曲線;最后,對汽車行駛工況與該城市所采集數(shù)據(jù)源(經(jīng)處理后的數(shù)據(jù))的各指標(biāo)(運(yùn)動(dòng)特征)值進(jìn)行誤差分析,分析所構(gòu)建的汽車行駛工況的合理性[13].
圖2 汽車行駛工況構(gòu)建的主要思路圖
2.2.1劃分條件 由于一個(gè)完整的運(yùn)動(dòng)學(xué)片段包含了怠速、加速、減速以及勻速4個(gè)狀態(tài),而如果持續(xù)時(shí)間小于20 s,那么認(rèn)為該片段汽車處于頻繁啟動(dòng)的狀態(tài),則該片段可能不能刻畫出該時(shí)間段汽車的行駛工況.因此,為了避免這種情況發(fā)生,認(rèn)為運(yùn)動(dòng)學(xué)片段的持續(xù)時(shí)間應(yīng)當(dāng)大于等于20 s;且將小于20 s的片段與其后一個(gè)片段合并為一個(gè)新的片段,直至合并成新片段的持續(xù)時(shí)間大于等于20 s,即認(rèn)為構(gòu)成了一個(gè)完整的運(yùn)動(dòng)學(xué)片段[14].
2.2.2劃分結(jié)果 利用MATLAB編程,根據(jù)劃分條件對預(yù)處理后的數(shù)據(jù)進(jìn)行進(jìn)一步的劃分,可得到文件1~文件3各文件的運(yùn)動(dòng)學(xué)片段數(shù)量記錄(表6).
表6 各文件的運(yùn)動(dòng)學(xué)片段數(shù)量記錄
顯然GPS車速原始數(shù)據(jù)中,最直觀的兩個(gè)汽車運(yùn)動(dòng)特征的指標(biāo)為時(shí)間和速度,而加速度和減速度則可由時(shí)間和速度計(jì)算得到;如果只采用這4個(gè)指標(biāo)描述汽車行駛工況,則可能會(huì)造成汽車行駛特征信息的丟失,不能真實(shí)、全面地反映其行駛工況.因此,需要引入以下運(yùn)動(dòng)特征指標(biāo)[15].
1)運(yùn)行時(shí)間
t=tt
,
(7)
式中:tt為運(yùn)動(dòng)學(xué)片段的總時(shí)間,即數(shù)據(jù)的總個(gè)數(shù).
2)最大速度vmax
vmax=max{vi,i=1,2,…,tt}.
(8)
3)平均速度vm
(9)
式中:s表示運(yùn)動(dòng)學(xué)片段的運(yùn)行距離,即該片段的距離之和.
4)平均行駛速度vmr.平均行駛速度為除了怠速外的平均速度,即
(10)
5)速度標(biāo)準(zhǔn)差Sv
(11)
6)加速度ai,i+1
(12)
式中:ai,i+1為第i秒到第(i+1) 秒的瞬時(shí)加速度;vi為第i秒的GPS車速;ti為第i秒.
7)最大加速度amax
amax=max{ai,i=1,2,…,tt},
(13)
式中:tt表示運(yùn)動(dòng)學(xué)片段的總時(shí)間,即數(shù)據(jù)的總個(gè)數(shù).
8) 平均加速度am
(14)
9) 加速度標(biāo)準(zhǔn)差Sa
(15)
10) 最大減速度vd,max
vd,max=max{vi,i=1,2,…,tt}.
(16)
11) 平均減速度vd,m
(17)
12) 怠速時(shí)間比例Pi
(18)
式中:ti表示運(yùn)動(dòng)學(xué)片段中ai>0.1 m/s2的時(shí)間.
13) 加速時(shí)間比例Pa
(19)
式中:ta表示運(yùn)動(dòng)學(xué)片段中vi=0的時(shí)間.
14) 減速時(shí)間比例Pd
(20)
式中:td表示運(yùn)動(dòng)學(xué)片段中ai<-0.1 m/s2的時(shí)間.
15) 勻速時(shí)間比例Pe
(21)
式中:te=t-ti-ta-td為運(yùn)動(dòng)學(xué)片段中勻速的時(shí)間.
結(jié)合以上分析,由于文件1~文件3為同一輛車在不同時(shí)間段所采集的數(shù)據(jù),因此,先將文件1~文件3整合成新的數(shù)據(jù),利用MATLAB編程,可得到各運(yùn)動(dòng)學(xué)片段運(yùn)動(dòng)特征評(píng)估體系指標(biāo)值(表7).
表7 各運(yùn)動(dòng)學(xué)片段運(yùn)動(dòng)特征評(píng)估體系指標(biāo)值
上述選取的指標(biāo)較多,由于某些指標(biāo)之間可能存在信息重疊的現(xiàn)象,因此采用主成分分析對運(yùn)動(dòng)特征指標(biāo)進(jìn)行降維處理.主成分分析[16—17]是在處理多個(gè)指標(biāo)時(shí)對指標(biāo)進(jìn)行降維的一種方法,其基本思想是將相互之間具有一定關(guān)聯(lián)性的多個(gè)指標(biāo),轉(zhuǎn)化為幾個(gè)互相無關(guān)的主要成分(即綜合指標(biāo)),可以大大降低分析和計(jì)算的難度.采用主成分分析對運(yùn)動(dòng)特征評(píng)估體系指標(biāo)進(jìn)行降維分析的具體算法描述見表8[18].
表8 各文件的運(yùn)動(dòng)學(xué)片段數(shù)量記錄
利用MATLAB編程進(jìn)行主成分分析,計(jì)算可得各運(yùn)動(dòng)特征指標(biāo)按從大到小排列的貢獻(xiàn)率,見表9.由表9可以看出,前5個(gè)特征根的累計(jì)貢獻(xiàn)率已經(jīng)接近90%,表明主成分分析效果很好.因此,可以認(rèn)為,前5個(gè)主成分已經(jīng)包含了14個(gè)運(yùn)動(dòng)特征指標(biāo)的大部分信息,即可選取前5個(gè)主成分進(jìn)行綜合評(píng)價(jià).前5個(gè)主成分對應(yīng)的特征向量見表10.
表9 主成分分析各指標(biāo)貢獻(xiàn)率
表10 標(biāo)準(zhǔn)化變量的前5個(gè)主成分對應(yīng)的特征向量
考慮到一個(gè)路段往往包含了不同的路況類型,將每個(gè)運(yùn)動(dòng)學(xué)片段按照擁堵路況、一般擁堵路況和通暢路況分為3類路況類型[19].將具有相同路況的運(yùn)動(dòng)學(xué)片段分為一類,同時(shí)可將選取的3類路況類型的片段組合成能體現(xiàn)汽車行駛特征的汽車工況曲線.
基于上述主成分分析所選取的5個(gè)主成分?jǐn)?shù)據(jù),采用K-means聚類分析算法對運(yùn)動(dòng)學(xué)片段進(jìn)行分類,該類算法是著名的劃分聚類算法[20—22].該算法由于其簡潔和高效,成為所有聚類算法中最廣泛使用的聚類算法,其核心思想是把n個(gè)數(shù)據(jù)對象劃分為K個(gè)聚類,使每個(gè)聚類中的數(shù)據(jù)點(diǎn)到該聚類中心的平方和最?。@然,可以根據(jù)路況類型采用K-means聚類將運(yùn)動(dòng)學(xué)片段分為3類,其具體的算法描述見表11.
表11 K-means聚類算法描述
經(jīng)過上述分析,利用MATLAB 編程,對5個(gè)主成分進(jìn)行K-means聚類分析可得到分類結(jié)果(表12).
表12 各文件的運(yùn)動(dòng)學(xué)片段數(shù)量記錄
由表12可知,3類工況運(yùn)動(dòng)學(xué)片段數(shù)比例約為
擁堵工況∶一般擁堵工況∶通暢工況=3∶6∶1.
因此,從第1類選取聚類偏差度最小的前3個(gè)運(yùn)動(dòng)學(xué)片段,從第2類選取聚類偏差度最小的前6個(gè)運(yùn)動(dòng)學(xué)片段,從第3類選取聚類偏差度最小的1個(gè)運(yùn)動(dòng)學(xué)片段,同時(shí)以總時(shí)間1 200~1 300 s進(jìn)行組合,總共選取了10個(gè)運(yùn)動(dòng)學(xué)片段(表13).
表13 各文件的運(yùn)動(dòng)學(xué)片段數(shù)量記錄
最終得到的工況曲線見圖3.
圖3 汽車行駛工況曲線
為了分析上述所構(gòu)建的汽車行駛工況的合理性,需要對汽車行駛工況與該城市所采集數(shù)據(jù)源(經(jīng)處理后的數(shù)據(jù))的各指標(biāo)(運(yùn)動(dòng)特征)值進(jìn)行誤差分析[19].不妨令
(22)
由 (22) 式計(jì)算可得到誤差分析統(tǒng)計(jì)表(表14).
由表14可以看出,汽車行駛工況與該城市所采集數(shù)據(jù)源(經(jīng)處理后的數(shù)據(jù))各指標(biāo)的最大相對誤差為平均減速度的誤差σi max=17.20%.最終構(gòu)建的汽車行駛工況,除了勻速時(shí)間比例、平均加速度以及怠速時(shí)間比例大于10%,其余運(yùn)動(dòng)特征指標(biāo)與該城市所采集數(shù)據(jù)的相對誤差均小于10%.這表明,上述所構(gòu)建的汽車行駛工況曲線是有效的,完全能夠體現(xiàn)參與數(shù)據(jù)采集汽車行駛特征的汽車行駛工況曲線,這驗(yàn)證了所構(gòu)建的汽車行駛工況的合理性.
表14 運(yùn)動(dòng)特征指標(biāo)誤差分析統(tǒng)計(jì)表
本文給出了汽車行駛的工況特點(diǎn)和規(guī)律,基于城市汽車實(shí)際道路行駛采集的數(shù)據(jù),通過對數(shù)據(jù)進(jìn)行時(shí)間不連續(xù)數(shù)據(jù)、尖點(diǎn)數(shù)據(jù)、毛刺噪聲數(shù)據(jù)以及怠速異常數(shù)據(jù)一系列預(yù)處理.采用主成分分析對運(yùn)動(dòng)特征指標(biāo)進(jìn)行降維處理,通過將較多的指標(biāo)轉(zhuǎn)化為幾個(gè)主要成分,可以大大降低分析和計(jì)算的難度.最終,結(jié)合主成分和K-means聚類分析,構(gòu)建了能客觀真實(shí)地反映汽車行駛特征的工況曲線.經(jīng)過汽車行駛工況與該城市所采集數(shù)據(jù)源(經(jīng)處理后的數(shù)據(jù))的各指標(biāo)(運(yùn)動(dòng)特征)值的誤差分析,驗(yàn)證了所構(gòu)建的汽車行駛工況的合理性.