聚類和馬爾科夫方法結(jié)合的城市汽車行駛工況構(gòu)建

姜 平 石 琴 陳無畏

合肥工業(yè)大學(xué),合肥,230009

0 引言

汽車行駛工況分為標(biāo)準(zhǔn)工況和非標(biāo)準(zhǔn)工況。當(dāng)今世界上很多國家都以標(biāo)準(zhǔn)等形式提出了不同車型在不同應(yīng)用時(shí)的標(biāo)準(zhǔn)工況,如日本的10-15工況循環(huán)、歐洲經(jīng)濟(jì)委員會(huì)的ECE-R15工況循環(huán)、美國聯(lián)邦城市及高速公路循環(huán)(CVS-C/H)。我國汽車排放測試規(guī)程采用歐洲的排放標(biāo)準(zhǔn),往往帶來很不滿意的結(jié)果,因此,需要對待開發(fā)的目標(biāo)車型進(jìn)行工況調(diào)查和測試,并得到合適的行駛工況。所以,開發(fā)和設(shè)計(jì)汽車、考核汽車排放污染水平都需要有相應(yīng)的行駛工況。

目前,國內(nèi)學(xué)者大多采用單純的短行程和聚類分析法對汽車行駛工況進(jìn)行研究[1-3]。國外很多學(xué)者采用“微路徑”以及隨機(jī)過程選擇方法來構(gòu)建行駛工況[4-5]。本文在構(gòu)建汽車行駛工況時(shí),不同于國內(nèi)傳統(tǒng)方法。傳統(tǒng)方法未考慮汽車燃油消耗與平均速度、發(fā)動(dòng)機(jī)工作狀態(tài)分布等與行駛工況有關(guān)的因素,只是籠統(tǒng)地對整個(gè)汽車行駛工況進(jìn)行構(gòu)建。本文將聚類分析法、因子分析法與馬爾科夫過程相結(jié)合,按照平均速度和速度標(biāo)準(zhǔn)偏差對汽車行駛工況進(jìn)行類的劃分,將每類行駛工況分別看作一個(gè)馬爾科夫過程,計(jì)算每類行駛工況的片段轉(zhuǎn)移概率和發(fā)動(dòng)機(jī)比功率,提出構(gòu)建候選工況的新方法。

1 因子分析

因子分析的基本目的就是用少數(shù)幾個(gè)因子描述許多指標(biāo)或因素之間的關(guān)系,把變量表示成各因子的線性組合,即將關(guān)系比較密切的幾個(gè)變量歸在同一類中,每一類變量就成為一個(gè)因子,從而以較少的幾個(gè)因子反映原資料的大部分信息。

發(fā)動(dòng)機(jī)工作狀態(tài)分布一般根據(jù)發(fā)動(dòng)機(jī)比功率、汽車行駛速度、瞬時(shí)加速度決定,而且它又與汽車排放、燃油消耗密切相關(guān),因此在構(gòu)建汽車行駛工況時(shí),應(yīng)考慮發(fā)動(dòng)機(jī)工作狀態(tài)分布。比功率綜合了速度、加速度、坡度以及風(fēng)阻等參數(shù),其物理意義為發(fā)動(dòng)機(jī)輸出功率與機(jī)動(dòng)車質(zhì)量的比值,其擬合公式如下[6]：

式中,v為車輛行駛速度,m/s;a為車輛行駛瞬態(tài)加速度,m/s2;θ為道路坡度,(°);Pa為發(fā)動(dòng)機(jī)比功率,k W/t。

參考國外發(fā)動(dòng)機(jī)工作狀態(tài)分布[7],考慮國內(nèi)混合交通實(shí)際情況,將我國城市典型道路行駛工況燃油效率劃分為表1所示的分布。表1的工作狀態(tài)分布是一個(gè)14維變量(B1～B14),這些變量分別由速度分布、比功率來確定。每個(gè)變量類似于一個(gè)單元,每一個(gè)單元類似于一個(gè)“小抽屜”,在這個(gè)“小抽屜”中,存放著該工況下測點(diǎn)的平均燃油消耗。行駛工況是逐秒采集的,每秒所觀察的數(shù)據(jù)被分配在表1中。如B1表示該工況下汽車制動(dòng)時(shí)平均燃油消耗量,B2表示該工況下汽車怠速時(shí)平均燃油消耗量,B3表示該工況下汽車比功率小于0且速度大于0、小于40km/h時(shí)的平均燃油消耗量,依此類推,其他變量表示在該工況下不同比功率和速度區(qū)間的平均燃油消耗量。

表1 行駛工況燃油消耗分布的劃分

本文應(yīng)用統(tǒng)計(jì)分析軟件SPSS15.0中的因子分析方法來減少變量,采用主成分分析法提取公因子,提取特征值大于1的成分為主成分。因子模型的旋轉(zhuǎn)方法選取方差最大變換法。對合肥市5條典型道路工作狀態(tài)分布進(jìn)行主成分分析,特征碎石圖見圖 1。前 2個(gè)主成分的特征值為1.47,累計(jì)方差貢獻(xiàn)率為84.3%,即前2個(gè)主成分包含了原有變量84.3%的信息,所以可以取這2個(gè)因子來代替原有的變量。

圖1 特征值碎石圖

通過對因子模型進(jìn)行旋轉(zhuǎn),得到旋轉(zhuǎn)后的因子載荷矩陣如表2所示。根據(jù)旋轉(zhuǎn)后的因子載荷大小來判別因子主要解釋的變量,本文選取載荷大于0.8的變量為因子的主要解釋變量。因子1主要解釋了變量 B9、B10、B11、B12、B13 、B14,由表1可知,這些變量主要是高速狀態(tài)下的變量。低速狀態(tài)變量因子載荷小,如 B0、B1、B2、B3、B4、B5、B6、B7,說明因子 1與平均速度有關(guān),可以解釋為平均速度。因子2主要解釋了變量B1、B3、B4、B5,對怠速狀態(tài)變量B2影響最小,載荷在低速時(shí)先增大,中間速度時(shí)減小,高速時(shí)先減小后增大,這表明因子2與車速偏離平均速度的程度有關(guān),可以解釋為速度的標(biāo)準(zhǔn)偏差。

表2 旋轉(zhuǎn)后的因子載荷矩陣

2 兩步聚類分析

兩步聚類分析是一種新型的分層聚類方法,目前一般應(yīng)用在數(shù)據(jù)挖掘與多元統(tǒng)計(jì)的交叉領(lǐng)域——模式分類中,其算法適合任何尺度的變量[8]。采用上述與燃油消耗相關(guān)的兩個(gè)主要因子(平均速度和速度標(biāo)準(zhǔn)偏差)進(jìn)行兩步聚類分析。聚類準(zhǔn)則采用貝葉斯信息準(zhǔn)則 Akaike(AIC)。AIC利用N-L平均信息量最小的原則來確定系統(tǒng)的階次和估計(jì)相應(yīng)參數(shù)值,其定階準(zhǔn)則為[9]：

式中,N為模型獨(dú)立未知參數(shù)的個(gè)數(shù);θML為相應(yīng)參數(shù)的估計(jì)值;L(θML)為在θML條件下的似然函數(shù)值。

在工程應(yīng)用中,對于不同階次 N,首先用極大似然估計(jì)法估計(jì)系統(tǒng)的參數(shù)值,然后計(jì)算似然函數(shù)值及其AIC(N)值,以使AIC(N)值達(dá)到最小時(shí)的階次作為該模型的階次。本文利用SPSS15.0統(tǒng)計(jì)軟件的兩步聚類模型,采用AIC準(zhǔn)則對合肥市典型道路的行駛工況進(jìn)行聚類分析,結(jié)果可以聚類為3類：第1類行駛工況(類1)平均速度偏低,為0～20km/h,第2類行駛工況(類2)平均速度中等,為 20～40km/h,第3類行駛工況(類3)平均速度偏高,大于 40km/h。

3 各類行駛工況的構(gòu)建

3.1 樣本的采集

本研究中試驗(yàn)車為轎車,試驗(yàn)車盡可能跟隨測試路段上的車流正常行駛。速度信號(hào)采樣頻率為1Hz。城市道路行駛工況的研究一般是從多條道路中篩選出有代表性的典型道路進(jìn)行研究。典型道路的特征是能綜合反映整個(gè)城市的主要交通狀況,即反映整個(gè)城市道路的道路等級(快速道、主干道、次干道及支路,以及車道數(shù)、機(jī)動(dòng)車非機(jī)動(dòng)車混合等),交通強(qiáng)度(交通流量、飽和度等),交叉口密度(路段內(nèi)交叉形式和數(shù)量等)和時(shí)間(早晚交通流量變化明顯程度等)。本文選取了合肥市勝利路、屯溪路、明光路等5條道路進(jìn)行樣本采集,使采集的數(shù)據(jù)能反映合肥市整體道路的交通情況,統(tǒng)計(jì)的汽車交通流量具有代表性。

3.2 數(shù)據(jù)的分離

數(shù)據(jù)分離就是把行駛工況的速度時(shí)間曲線劃分成許多反映加速度變化的模型事件,包括加速模型事件、減速模型事件、怠速模型事件和勻速模型事件,這些模型事件是構(gòu)建行駛工況的基礎(chǔ)。模型事件使用(a1,a2,δ,n)準(zhǔn)則進(jìn)行定義。為了確定分類的最優(yōu)值,需要考慮以下目標(biāo)函數(shù)和約束條件。

目標(biāo)函數(shù)CAP是一個(gè)加權(quán)函數(shù)：

式中,k1,k2,k3,k4,k5為目標(biāo)函數(shù)的權(quán)值。

為實(shí)現(xiàn)上述目標(biāo),要滿足以下5個(gè)約束條件：①加速度模型事件含有連續(xù)的加速度,因此,加速標(biāo)準(zhǔn)差θ1應(yīng)該很小;②對于任何勻速事件,最大速度和最低速度差θ2應(yīng)該小;③對于任何加速或減速事件,最大速度和最低速度差θ3應(yīng)大;④任何勻速事件的平均加速度θ4接近0;⑤因?yàn)橐粋€(gè)單獨(dú)點(diǎn)的特征依賴于相鄰點(diǎn)特征,所以一個(gè)好的分類標(biāo)準(zhǔn)應(yīng)該包括一小部分單秒數(shù)據(jù)θ5。

通過迭代法使目標(biāo)函數(shù)CAP最小,得到最優(yōu)的分類準(zhǔn)則——(5,1,3,2)分類原則[10],即加速事件定義為單一瞬時(shí)加速度不小于2.22m/s2或者連續(xù)觀測加速度大于0.71 m/s2,或持續(xù)3s或以上時(shí)間,而且要實(shí)現(xiàn)3.2km/h以上的速度遞增;減速片段定義為單一瞬時(shí)加速度不大于-2.22m/s2,或任何連續(xù)觀測加速度不大于-0.71m/s2,持續(xù) 3s或以上時(shí)間,而且實(shí)現(xiàn)3.2km/h以上的速度遞減。其余的被列為怠速或勻速事件。這種分類原則考慮了現(xiàn)實(shí)世界的行駛工況,勻速模式事件可以包括一些小的加速和減速模型事件。按照這個(gè)分類原則,行駛工況片段劃分結(jié)果如圖2所示。怠速模型事件指速度為0的模型事件。

圖2 速度時(shí)間曲線片段劃分結(jié)果

3.3 轉(zhuǎn)移矩陣的計(jì)算

馬爾科夫過程是一個(gè)隨機(jī)過程,記為Zτ,τ=1,2,…,T,Zτ代表每個(gè)狀態(tài)的模型事件。當(dāng)前時(shí)刻τ狀態(tài)的概率只與前一個(gè)時(shí)刻τ-1狀態(tài)有關(guān)。換句話說,汽車下個(gè)時(shí)間τ+1的狀態(tài)(如加速、減速、勻速和怠速片段)的概率只與當(dāng)前時(shí)間τ的加速、減速、勻速和怠速片段有關(guān)。

在馬爾科夫鏈中,假設(shè)X是方程X=pX時(shí)(p為轉(zhuǎn)移矩陣,X為p的特征向量)的唯一解,則X為馬爾科夫鏈的靜態(tài)分布。對于一個(gè)靜態(tài)的馬爾科夫鏈,隨著馬爾科夫鏈的增長,整個(gè)模型事件分布保持不變。應(yīng)用這個(gè)性質(zhì),根據(jù)模型事件轉(zhuǎn)移矩陣構(gòu)建的代表性行駛工況能代表整個(gè)試驗(yàn)數(shù)據(jù)的行駛工況。

為了估計(jì)每類行駛工況的轉(zhuǎn)移矩陣,運(yùn)用Lee等[11]的最大似然估計(jì)法進(jìn)行估計(jì)。對于一個(gè)固定的馬爾科夫過程,應(yīng)用貝葉斯公式得到連續(xù)的有序狀態(tài)的概率為

假定每個(gè)過程都是重復(fù)觀察的,N rs為τ-1時(shí)刻狀態(tài)r轉(zhuǎn)移到τ時(shí)刻狀態(tài)s的事件數(shù),則N次重復(fù)觀察的公式為

根據(jù)極大似然函數(shù),得到狀態(tài)轉(zhuǎn)移概率方程

依據(jù)上述原理,得到合肥市典型道路不同類別的行駛工況的片段轉(zhuǎn)移概率,如表3所示。例如,在類1工況中,若當(dāng)前狀態(tài)是減速片段,則下個(gè)狀態(tài)為加速片段的概率是53%,為勻速片段的概率是38%,為怠速片段的概率是9%。類3工況常為高速工況,通常沒有怠速,故只有3個(gè)狀態(tài)。

表3 片段轉(zhuǎn)移概率 %

3.4 候選工況的構(gòu)建

構(gòu)建行駛工況的基本思路是截取若干段實(shí)測的數(shù)據(jù)片段,拼成一個(gè)完整的速度—時(shí)間歷程。隨機(jī)選擇一個(gè)小于100s的片段作為開始片段,即初始狀態(tài)。初始狀態(tài)的起始速度可以為0。下個(gè)狀態(tài)的選擇,即片段的選擇一般根據(jù)速度、發(fā)動(dòng)機(jī)比功率和轉(zhuǎn)移概率來選擇。依照以下步驟構(gòu)建候選工況：

(1)初始狀態(tài)一旦確定,它的終點(diǎn)速度和比功率就確定了,假設(shè)為當(dāng)前狀態(tài)。下個(gè)狀態(tài)的初始速度與上個(gè)片段最終速度絕對值誤差小于0.5km/h,而且下個(gè)狀態(tài)初始發(fā)動(dòng)機(jī)比功率與當(dāng)前狀態(tài)的終點(diǎn)比功率相比不能超過15k W/t。

(2)根據(jù)轉(zhuǎn)移矩陣和上面提到的速度與比功率原則,選擇下個(gè)狀態(tài)。

(3)依據(jù)上述原則,不斷選擇下一個(gè)狀態(tài),直到候選工況達(dá)到理想的長度為止。

3.5 代表性工況的選取

為了從眾多候選工況中選取合適的代表性行駛工況,本文考慮2個(gè)參數(shù),即候選工況與試驗(yàn)數(shù)據(jù)的速度—加速度頻率差 d 1,4種模型事件(加速、減速、怠速、勻速)時(shí)間比例差值的絕對值d2：

式中,fr(a,v)為候選工況速度v和加速度a的比值;p(a,v)為試驗(yàn)數(shù)據(jù)中速度v和加速度為a的比值;fm為候選工況中第m種模型事件的時(shí)間比值;pm為試驗(yàn)數(shù)據(jù)中第m種模型事件的時(shí)間比值。

速度 —時(shí)間頻率按照速度間隔5km/h、加速度間隔0.278m/s2進(jìn)行計(jì)算。

綜合考慮 d1與 d2,選擇綜合加權(quán)測試法(composite performance measure,CPM)進(jìn)行評價(jià),即取 CPM最小的候選工況為代表性行駛工況：

式中,w1為d1的標(biāo)準(zhǔn)偏差的倒數(shù);w2為d2的標(biāo)準(zhǔn)偏差的倒數(shù)。

4 結(jié)果分析

以合肥市5條典型道路為例,采用因子分析和兩步聚類分析方法把汽車行駛工況分成3類,再分別對每類行駛工況按照馬爾科夫方法構(gòu)建代表性行駛工況,并與試驗(yàn)數(shù)據(jù)進(jìn)行比較。特征參數(shù)包括平均速度v m、速度標(biāo)準(zhǔn)方差σv、正加速度的平均值av+、負(fù)加速度的平均值av-、怠速時(shí)間比例p i、加速時(shí)間比例 p a、減速時(shí)間比例 p d、勻速時(shí)間比例p c、正的發(fā)動(dòng)機(jī)比功率平均值P+、負(fù)的發(fā)動(dòng)機(jī)比功率平均值P-、發(fā)動(dòng)機(jī)比功率平均值Pa、發(fā)動(dòng)機(jī)比功率方差σvsp。

每類行駛工況所構(gòu)建的代表性行駛工況與試驗(yàn)數(shù)據(jù)的特征參數(shù)比較見表4。可以看出,低速行駛工況怠速比例高,高速行駛工況怠速比例為0。類1、類2和類3所構(gòu)建的代表性行駛工況與試驗(yàn)數(shù)據(jù)特征參數(shù)平均誤差分別為3.47%、1.26%和1.28%。表明構(gòu)建的行駛工況與試驗(yàn)數(shù)據(jù)相符,比較合理。

表4 試驗(yàn)數(shù)據(jù)與代表性行駛工況的特征參數(shù)比較

為了檢驗(yàn)各類行駛工況模型事件分布是否與試驗(yàn)數(shù)據(jù)分布一致,將3類代表性行駛工況的片段轉(zhuǎn)移概率與試驗(yàn)數(shù)據(jù)的片段轉(zhuǎn)移概率進(jìn)行χ2檢驗(yàn),3類檢驗(yàn)概率P值都大于0.05,說明構(gòu)建行駛工況與試驗(yàn)數(shù)據(jù)來自同一樣本。構(gòu)建后的代表性行駛工況的χ2檢驗(yàn)統(tǒng)計(jì)值、d 1、d 2和CPM值如表5所示。

表5 代表性行駛工況周期、評價(jià)參數(shù)和χ2分布檢驗(yàn)

類1、類2、類3三類工況構(gòu)建的代表性行駛工況分別如圖3、圖4和圖5所示。

圖3 類1代表性行駛工況

圖4 類2代表性行駛工況

圖5 類3代表性行駛工況

5 結(jié)論

(1)本文先采用兩步聚類分析方法把行駛工況分成低速、中速、高速3類行駛工況,然后采用馬爾科夫方法構(gòu)建行駛工況。構(gòu)建行駛工況時(shí)考慮了與燃油消耗相關(guān)的因素(國內(nèi)發(fā)動(dòng)機(jī)工作狀態(tài)分布),因此構(gòu)成的代表性行駛工況能較好地估計(jì)汽車的實(shí)際燃油消耗。

(2)利用新的分類準(zhǔn)則將行駛片段聯(lián)系為特殊的模型事件來定義行駛變量,將速度—時(shí)間曲線劃分為加速、減速、勻速和怠速片段,根據(jù)速度、比功率和轉(zhuǎn)移矩陣求取各片段間的轉(zhuǎn)移概率,建立了新的特征參數(shù)評價(jià)準(zhǔn)則和代表性行駛工況的選取原則。

(3)利用組合方法構(gòu)建的代表性行駛工況與試驗(yàn)數(shù)據(jù)的特征參數(shù)相比,平均相對誤差小,即構(gòu)建的代表性工況能準(zhǔn)確地反映城市道路上汽車的實(shí)際行駛工況。

[1] 呂林,鎖國濤.不同路況下汽車行駛工況的測試與分析[J].武漢理工大學(xué)學(xué)報(bào),2006,30(5)：781-784.

[2] 張建偉,李孟良,艾國和,等.車輛行駛工況與特征的研究[J].汽車工程,2005,27(2)：220-224.

[3] 杜青,楊延相,鄭偉,等.機(jī)動(dòng)車實(shí)際道路排放特性及若干影響因素的研究[J].內(nèi)燃機(jī)學(xué)報(bào),2002,20(4)：297-302.

[4] Lin J,Niemeier D A.An Exploratory Analysis Comparing a Stochastic Driving Cycle to California's Regulatory Cycle[J].Atmospheric Environment,2002,36(38)：5759-5770.

[5] Eisinger D S,Niemeier D A,Stoechenius T,et al.Collecting Driving Data to Support Mobile Source Emissions Estimation[J].Journal of Transportation Engineering,2006,132(11)：845-854.

[6] 姚志良,賀克斌,王岐東,等.IVE機(jī)動(dòng)車排放模型應(yīng)用研究[J].環(huán)境科學(xué),2006,27(10)：1928-1933.

[7] John K,Mitch C,Megan B,et al.MOVES2004 Software Design Reference Manual Draft[M].Ann Arbor：Environmental Protection Agency,2004.

[8] 陳超,鄒瀅.SPASS 15.0中文版常用功能與應(yīng)用實(shí)例精簡[M].北京：電子工業(yè)出版社,2009.

[9] 張濟(jì)民,王開文,池茂儒.基于AIC信息準(zhǔn)則法的擺式列車傾擺伺服系統(tǒng)建模研究[J].中國鐵道科學(xué),2003,24(6)：6-9.

[10] Zhen Dai,Niemeier D A,Eisinger D.Driving Cycles：a New Cycle-building Method that Better Represents Real-world Emissions.Task Order No.66[R].Berkeley：Department of Civil&Environmental Engineering University of California,2008.

[11] Lee T C,Judge G G,Zellner A.Estimating the Parameters of the Markov Probability Model from Aggregate Time Series Data[M].Amsterdam：North-Holland Publishing Company,1970.