基于短行程特征聚類的移動源排放清單構建方法

李兵兵, 康 宇, 曹 洋, 李亞民, 許鎮義

(1.安徽省生態環境監測中心,安徽 合肥 230071; 2.合肥綜合性國家科學中心人工智能研究院,安徽 合肥 230088; 3.中國科學技術大學 自動化系,安徽 合肥 230026; 4.中國科學技術大學 先進技術研究院,安徽 合肥 230088)

0 引 言

近年來,隨著我國城市化進程和社會經濟的快速發展,機動車保有量快速增長,移動源排放會造成細顆粒物和光化學煙霧污染,已成為我國城市大氣污染重要來源。排放清單作為移動源污染定量核算的重要手段,可以用于大氣污染防治措施制定與溯源分析。

車輛行駛工況是針對某一類車輛,在特定區域和時間段下,描述該類車輛行駛規律的速度-時間曲線。隨著車速增加,燃油車細顆粒物排放物質量濃度會明顯升高,在車輛加速情況下,車輛排放的顆粒物數量、質量濃度也明顯增高[1]。通過研究道路移動源行駛工況來估計車輛排放,對移動源排放清單構建具有重要參考價值。

美國和歐洲開展汽車排放標準研究較早,已建立符合其交通特點與行駛狀態的行駛工況,如新標歐洲循環測試(New European Driving Cycle,NEDC)工況[2]、美國聯邦試驗規程75(Federal Test Procedure 75,FTP-75)工況[3]。NEDC工況測試時間較短,速度分布比較規則,行駛里程較短且加速度比較穩定,其循環時長為1 180 s,行駛里程為11.007 km,平均速度為33.6 km/h,主要在歐洲、中國、澳大利亞等國家和地區使用。FTP-75工況是美國環境保護署用于排放標準制定的測試工況,其由車輛冷啟動階段、車輛熱穩定階段和車輛熱啟動階段3部分組成,工況總時長為1 874 s,行駛里程為17.77 km,平均速度為34.1 km/h,主要使用的國家和地區包括美國、加拿大、南美洲等。

隨著城市路網規模不斷擴大,車輛高速行駛時間不斷增加,歐美和日本排放實驗領域專家共同制定了全球統一的輕型車測試工況(Worldwide harmonized Light-duty vehicle Test Cycle,WLTC)[4]。WLTC工況主要包括低速階段、中速階段、高速階段和額外高速階段,工況運行時長為1 800 s,運行距離為23.3 km,平均速度為46.3 km/h,其充分考慮車輛性能差異,對不同功率質量比的車輛采用不同測試循環,能夠客觀真實地反映車輛的實際排放狀況。

當前行駛工況合成方法主要包括基于馬爾科夫決策的工況構建方法和基于聚類算法的短行程構建方法。

文獻[5]基于馬爾科夫決策方法,使用特定的選取原則進行片段選取,形成聚類片段庫,從而從物理意義上反映車輛的行駛規律;該方法適用范圍廣泛,不限道路類型,但在起始和終止片段的選取方面缺乏統一的選取原則。而起始和終止片段的選取,特別是起始片段的選取,將直接影響所構建工況模型的精確性。

文獻[6]首次提出基于聚類算法的短行程構建方法,該方法可以有效地與主成分分析(principal component analysis,PCA)[7]和K-means[8]等聚類方法相結合,自動選取行程片段,因此受到更多的關注。

在國內,文獻[9]針對不同測試循環工況下的輕型車污染物排放特性進行研究。

由于我國汽車工業起步較晚,目前在工況測試排放標準研究方面大多沿用歐洲的NEDC工況,但是該工況與我國目前的人員駕駛習慣及車輛標準有一定出入。因此,根據我國實際道路駕駛數據,構建符合當地實際路況的行駛工況具有現實意義,對于車輛排放物及燃料消耗等測試與評價有重要參考作用。

同時,現有的移動源排放清單研究大多基于文獻[10]中的排放因子本地化修正方法,采用歐洲環境署開發的宏觀排放模型COPERT(Computer Program to calculate Emissions from Road Transport model)[11]、美國環境保護署應用于汽車排放源計算的IVE模型(International Vehicle Emissions Model)[12]和美國環境保護署開發的移動車輛排放模型(MOtor Vehicle Emission Simulator,MOVES)[13]等排放因子模型方法,初步探討了全國[14]、區域[15-17]、省(直轄市)[18-21]、省會城市[22-23]等多個尺度的機動車污染物排放結構組成、時空分布特征及相應減排情景對策等。

本文采用基于PCA和K-means聚類的短行程法,構建福州市輕型車的行駛工況。首先對原始數據進行預處理、劃分運動學片段并計算其特征參數,運用PCA對劃分后片段進行參數降維,再對其進行K-means聚類分析形成典型片段庫;然后分別根據最小參數偏差法和最大相關系數法挑選片段組成行駛工況,對比后輸出最終工況,并對最終行駛工況進行有效性驗證;最后利用VT-Micro(Virginia Tech microscopic)模型并結合合成工況,計算單車排放因子,構建2020年福州市機動車排放清單。

1 基于短行程特征工況的清單構建

本文提出一種基于短行程特征聚類的移動源行駛工況排放清單構建方法,具體步驟如下:

1) 對獲取的車輛GPS軌跡數據進行預處理。

2) 對軌跡的速度數據進行片段劃分得到運動學片段,并根據運動學片段特征計算車輛運動學特征參數和總數據的總體分布參數。

3) 利用PCA對運動學片段進行運動學特征參數降維。

4) 使用K-means聚類對輸入片段進行分類,得到多種類型的片段庫。

5) 分別按照最大相關系數法和最小參數偏差法從各類片段庫中挑選輸入片段,構建合成行駛工況,并對比2種方法所構建工況的綜合指標。

6) 引入VT-Micro模型,計算合成工況排放清單。

1.1 數據采集

本文采用未規定路線的自主行駛法采集福州市出租車行駛數據。不限定出車時間及地點,將GPS定位的終端設備安裝于實驗車輛上,進行數據遠程傳輸和存儲。為使得采集的輕型車行駛數據可以代表福州市乘用車行駛特點,行駛路線涵蓋福州市主要道路。

福州市區分為5個行政區,下轄5個縣及2個縣級市。5個行政區分別為鼓樓區、臺江區、晉安區、倉山區及馬尾區,前3個為福州市主城區。5個縣包括閩侯縣、閩清縣、永泰縣、羅源縣及連江縣,2個縣級市為福清市、長樂市。福州市區道路網呈“棋盤式+環路”分布。

考慮到大部分車輛行駛軌跡都集中于二環以內,且大部分的醫院、商場、學校等也設于二環內,本文選取的主要實驗區域為鼓樓區、臺江區、晉安區及倉山區;考慮到福州大學、閩江學院及馬尾區也會占有城市部分交通量,選擇二環快速、三環快速、福州繞城高速及主城區內的主干道作為實驗主要線路。

利用車載自動診斷系統(On-Board Diagnostics,OBD)接口采集發動機運行狀態信息,車輛運行參數包括軌跡點采樣時間、GPS車速、X軸加速度、Y軸加速度、Z軸加速度、軌跡采樣點經緯度位置坐標、發動機轉速、扭矩百分比、瞬時油耗、油門踏板開度、空燃比、發動機負荷百分比和進氣流量,累計采集331 550條GPS軌跡數據,采樣時間間隔為1 s,車輛運行參數的4種統計數據見表1所列。

表1 車輛運行參數4種統計數據

1.2 數據預處理

需要對采集的原始GPS軌跡和車輛信息數據進行預處理,以去除無效片段及各類異常值,保證用于構建工況的行駛數據有效可靠。

考慮到軌跡數據采樣時間間隔為1 s,車輛運行前后時刻速度不會突變過大,因此可以根據前后時刻的速度來填充當前時刻速度,若前后時刻速度全部缺失,直接刪除該數據。速度缺失值處理策略見表2所列,插值處理結果如圖1所示。

圖1 GPS車速缺失記錄插值處理結果

表2 速度缺失值處理策略

經車載設備直接采集的原始數據會包含部分異常數據。由于高層建筑或隧道遮擋,部分時刻GPS信號丟失,使得采集數據中的時間不連續,對于此類數據直接進行丟棄處理。若車輛長時間低速且間斷行駛(速度小于10 km/h),則可視為怠速狀態,默認最長怠速時間為180 s,多于180 s則按異常值進行丟棄。

1.3 運動學片段劃分與特征參數計算

由于行駛過程中車輛速度隨時間的變化特點不同,可以根據加速度和速度將不同類型的速度-時間片段定義為不同的運動狀態,便于后續分析和聚類。本文定義以下4種運動狀態:

1) 怠速狀態。車速值為0但發動機轉速值不為0的運動狀態。

2) 加速狀態。車輛加速度大于0.15 m/s2的運動狀態。

3) 減速狀態。車輛加速度小于-0.15 m/s2的運動狀態。

4) 勻速狀態。車輛加速度在-0.15～0.15 m/s2之間且車速值不為0的運動狀態。

將運動學片段定義為前一個停車時刻(速度值降為0)到下一個停車時刻之間的速度-時間曲線,因此運動學片段會由1段速度值為0的怠速片段和1段速度值不為0的行駛片段組成。

由于從上一停車時刻到下一停車時刻,速度必須經過上升和下降過程,車輛必然會出現加速、減速、勻速等狀態,運動學片段一般會出現怠速、加速、減速、勻速4種運動狀態。典型運動學片段組成如圖2所示。不同的運動學片段是車輛在不同地段不同時間段的行駛數據,其速度-時間曲線的時間長度及形狀特點存在差異,而每個運動學片段都由若干個速度-時間點構成,通過定義各類特征值來區分不同曲線,也能定量刻畫不同曲線之間的差異性,為后續的聚類分析提供依據。

圖2 車輛運動學片段組成

同時,在構建最終行駛工況后,還需要一些指標來定量分析和評價其有效性和準確性,以及與車輛真實行駛情況的貼近性。

因此,選取的特征參數必須盡可能全面、準確地描述運動學片段的運動特點。根據參數用途,將定義的特征參數分為2個大類:① 車輛運動特征參數,主要用來描述運動學片段的運動特征,可作為后續聚類分析的分類依據;② 數據總體分布特征參數,主要用于最終構建工況的檢驗與有效性驗證。選取的15個車輛運動特征參數和10個數據總體分布特征參數分別見表3、表4所列。

表3 車輛運動特征參數及其含義

表4 數據總體分布特征參數及其含義

首先對加速度特征參數進行計算,在采樣數據中,陀螺儀加速度是無法直接使用的,需要根據速度數據單獨計算車輛加速度,即

(1)

其中:Δt為采樣時間間隔,Δt=1 s;vt+1、vt分別為t+1、t時刻速度;at+1為t時刻到t+1時刻的加速度。

其次計算各類行程時長特征,計算公式為:

T=N

(2)

ta=Na

(3)

td=Nd

(4)

tc=Nc

(5)

te=N-Na-Nd-Nc

(6)

其中:N為運動學片段的總點數,考慮數據采樣時間間隔為1 s,采樣點數等于運動時長;Na為加速度大于0.15 m/s2的總點數;Nd為加速度小于-0.15 m/s2的總點數;Nc為速度為0的總點數。

然后依次計算S、vmax、vm、vmr、vsd、amax、amin、aa、ad及asd,計算公式分別為:

(7)

vmax=max{v1,v2,…,vN}

(8)

(9)

(10)

amax=max{a1,a2,…,aN}

(11)

amin=min{a1,a2,…,aN}

(12)

aa=

(13)

ad=

(14)

(15)

將所有采樣樣本按照式(1)～式(15)進行參數計算,得到2 818個運動學片段的車輛運動特征參數和樣本的總體分布特征參數。

4種運動狀態時間占比分別為:

加速,25.02%;減速,20.35%;

勻速,30.43%;怠速,24.19%。

速度分布占比分別為:

p[0,20),43.35%;p[20,40),26.43%;

p[40,60),18.99%;p[60,80),6.49%;

p[80,100),4.10%;p[100,120),0.60%。

由運動狀態分布結果可知,福州市輕型車的加速、減速、勻速和怠速4類狀態時間分布較為均勻,其中勻速時間段比例最大(占比30.4%),由此可見福州市輕型車的行駛較平穩,勻速駕駛狀態占比較大。而根據速度分布結果可知,大部分時間車輛速度處于[0,40) km/h的中低速區間,高速([60,120) km/h)占比小(占比11.2%),具有較強的區域特征。因此,構建特定區域或時間段行駛工況對于車輛行駛特征研究是必要的。

在聚類后,會生成K個(K為聚類類別數)不同的典型片段庫,在進行片段挑選構成最終工況時,需要計算能代表典型片段庫的綜合參數,以反映該片段庫的車輛行駛特點。基于一定挑選原則,根據每個片段車輛運動學參數和典型片段庫整體參數的接近程度,來挑選構成最終工況的運動學片段,采用相關系數法挑選片段時需要片段綜合參數。本文選取pa、pd、pc、pe、vmax、vm、vmr、amax、amin、p[0,20)、p[20,40)、p[40,60)、p[60,80)共13個指標作為片段庫的綜合參數。考慮到速度段[80,120) km/h的數據占比不足8%,納入綜合指標會導致較大偏差,因此,保證[0,80) km/h區間的速度分布貼合,即可保證構建工況的準確性。將典型片段庫中每種狀態時間總和除以典型片段庫的總運行時間得到各類狀態占比。

1.4 基于PCA和K-means的典型工況構建

短行程法將不同運動學片段根據特征參數值相似度分為幾類,形成能反映不同交通特點和行駛狀況的典型運動學工況片段庫,再根據一定的原則從庫中挑選出片段組合,形成最終的工況曲線。本文先通過PCA將冗余參數去除,只留下幾個主要特征組合,再通過運動學片段聚類構建典型工況片段。

利用PCA對運動學片段進行運動學特征參數降維,樣本數據包含m個個體、n個評價指標,則相應指標矩陣為:

(16)

重新組合成的新綜合指標為:

(17)

在新的綜合指標中,Fi、Fj為不相關特征向量。方差達到最大并要求其依次遞減,滿足的限制條件為:

(18)

在進行PCA之前,需對各類指標進行標準化處理,消除單位影響。標準化計算公式為:

(19)

(20)

(21)

則標準化后的矩陣為:

X=[X1X2…Xp]

(22)

相關系數矩陣為:

(23)

其中,rij為Xi與Xj的相關系數。rij計算公式為:

(24)

(25)

計算出相關系數矩陣的特征值λ及特征向量E,則有:

|R-λE|=0

(26)

主成分貢獻占比fi、主成分累計貢獻率ak計算公式分別為:

(27)

當ak>85%時,主成分可以代表絕大部分信息,主成分值可計算為:

Z=[Xu1Xu2…Xuk]=[Z1Z2…Zk]

(28)

其中,uk為對應的主成分特征向量。

1.5 基于合成工況的單車排放因子估算

(29)

其中:v為車輛的平均速度vm;Rv,j為車輛在vm下的第j種污染物排放因子;Ri,j為第i個比功率區間的第j種污染物排放率;Pi,v為vm在第i個比功率區間的分布值。

由于無法直接獲取機動車污染物排放率數據,需要利用其瞬時速度和加速度數據進行估算。VT-Micro 模型能夠基于車輛的瞬時速度和加速度計算車輛瞬時排放,其結構簡單,被廣泛用于交通研究中[24]。利用 VT-Micro模型進行交通排放測算,計算公式為:

(30)

其中:v為車輛的瞬時速度;a為車輛的瞬時加速度;Rv,a為車輛在v、a下的排放率;αi,j為vi和aj的回歸系數,根據美國橡樹嶺國家實驗室收集的車輛排放數據[24]得到的回歸系數取值見表5所列。

表5 回歸系數αi,j取值

(31)

2 實驗和討論

對劃分的運動學片段采用PCA進行分析,得到前4個主成分M1～M4,其累計貢獻率為88.64%,大于85%,則前4個主成分可以代表樣本的大部分信息。主成分方差及貢獻率見表6所列,前4個主要成分的載荷矩陣見表7所列。

表6 車輛運動特征參數主成分方差、貢獻率計算結果

表7 前4個主成分特征載荷矩陣

對劃分的運動學片段進行K-means聚類,得到典型片段庫,并根據戴維森-堡丁指數(Davies-Bouldin index,DBI)選擇最優聚類數。DBI通過計算距離來度量相同類之間的相似性、不同類之間的差異性,計算公式為:

(32)

其中:k′為簇別數;avgCi為第i類樣本點到簇中心ui的歐式距離平均值;dcen(ui,uj)為第i類簇中心ui與第j類簇中心uj之間的歐式距離。DBI越小,表示相同類運動學片段之間的相似性越高、不同類運動學片段之間的差異性越大,聚類效果越好。選取DBI值第1次出現局部最小值時的k′值,作為最優聚類數k0。根據最小參數偏差法和最大相關系數法,分別對各類片段進行排序挑選,計算各類片段庫占總體片段庫的時間比,確定各類片段庫中選取的片段時長,根據排序優先級挑選片段形成行駛工況。最終選擇聚類數目k′=3,聚類效果通過t分布隨機鄰近嵌入(t-distributed stochastic neighbor embedding,t-SNE)方法可視化,如圖3所示。

圖3 運動學片段聚類t-SNE降維可視化

進一步分析各類庫所代表的車輛行駛特征,根據參數整體分布情況,可將運動學片段分為低速(類別2)、中速(類別0)、高速(類別1)3類,計算典型運動學片段庫綜合參數,結果見表8所列。

表8 3類典型片段庫綜合參數

類別1運動學片段勻速比例大,速度最大值高達117.100 km/h,且中高速([40,80) km/h)及高速([80,120) km/h)占比最大,可定義為高速行駛片段,且加減速時間均勻、勻速比例大,路況暢通,屬于高速公路或快速道行駛。

類別2運動學片段低速占比(82.74%)非常大,怠速比例大,平均速度僅8.370 km/h,可定義為低速行駛片段,屬于市內擁堵路段或低速次干道行駛。

類別0運動學片段以中低速([0,40) km/h)為主,4種運動狀態時間均勻,平均行駛速度適中,最大速度介于類別1與類別2之間,為中速行駛類型,屬于主干道行駛。低速行駛時長(84 516 s)占比最大,中速時長(63 424 s)占比次之,高速時長(35 673 s)占比最少,可見福州市乘用車整體以低速擁堵行駛情況為主,主干道行駛情況較少,高速公路及快速道行駛情況極少。

所有運動學片段被分為高速、中速、低速3類典型片段庫。為客觀反映車輛真實的駕駛特性,最終的工況應該涵蓋這3類典型運動學片段,并且在最終工況中的各類片段時間比例與整體數據中比例一致。本文采用總體最大相關系數法(方法1)與最小參數偏差法(方法2)2種方法進行片段挑選,3類典型片段2種方法選擇結果如圖4～圖6所示。

圖4 高速片段2種方法選擇結果

圖5 低速片段2種方法選擇結果

圖6 中速片段2種方法片段選擇結果

對于2種片段選擇方法構建出的工況,用平均相對誤差來檢驗工況的有效性。平均相對誤差即所有參數相對誤差的平均值。

2種片段選擇方法所構建工況的相對誤差見表9所列。

表9 2種方法片段合成相對誤差

由表9可知:最大相關系數法(方法1)和最小參數偏差法(方法2)的平均相對誤差分別為7.49%、10.88%,最大相關系數法的平均相對誤差更小,其整體工況更貼近于真實駕駛情況;在運動狀態占比和速度占比方面,最大相關系數法明顯更加貼近總體數據。因此,選擇最大相關系數法輸出工況為最終工況。

最終構建的工況參數在整體上和真實行駛情況很貼近,尤其是在速度段[0,40) km/h比例與4種運動狀態比例方面吻合度極高,在速度段[60,80) km/h和速度標準差方面吻合度欠佳。

根據式(31)計算出單車合成工況排放因子,見表10所列。

表10 單車合成工況排放因子計算結果

根據文獻[25],2020年福州市汽車保有量為1 498 695輛,公路通車里程合計11 616.816 km,得到2020年福州市機動車污染物排放清單,見表11所列。

表11 2020年福州市機動車年均排放清單 單位:t

3 結 論

本文提出一種基于短行程特征聚類的移動源行駛工況排放清單構建方法,采用基于PCA和K-means聚類的短行程法,構建福州市輕型車的行駛工況。首先對原始數據進行預處理、劃分運動學片段并計算其特征參數,運用PCA對劃分后片段進行參數降維,再對其進行K-means聚類分析,形成典型片段庫;然后分別根據最小參數偏差法和最大相關系數法挑選片段組成行駛工況,對比后輸出最終工況,并對最終行駛工況進行有效性驗證;最后利用VT-Micro 模型并結合合成工況,計算單車排放因子,構建2020年福州市機動車排放清單。主要結論如下:

1) 相比于傳統的隨機選擇法,本文采用最大相關系數法和最小參數偏差法,分別以參數向量整體相似度最高和參數數值總體偏差最小為標準進行最終工況的片段挑選。

2) 對于選定參數作出較大調整,主要是以整體參數為準,去掉p[100,120)([100,120) km/h速度段占比不足1%,會引入較大誤差)、vmax、amin等參數,目的是追求整體吻合度更高,大部分數據貼合真實情況,而不是追求極值或者少部分數據擬合效果。

3) 在實驗線路及區域選取時,根據車流量分布特點,按一定比例選取主城區與郊區實驗車輛;并采用“不選定路線的自主行駛法”采集行駛數據,駕駛員可不受限制地自主駕駛,采集數據更具真實性。

本文雖然考慮行駛工況信息進行排放清單構建,但采用的VT-Micro排放因子模型仍然是一種基于統計回歸的數學估算模型,最終的排放清單核算與實際情況仍然存在較大偏差。在未來的研究中,可考慮利用便攜式排放測試系統(portable emissions measurement system,PEMS)獲取的實際道路排放監測數據,對現有排放因子估算模型進行修正,以獲得較準確的排放清單核算。