孫 劍，馮 羽

（同濟大學 道路與交通工程教育部重點實驗室，上海201804）

車輛出行OD（origin－destination）是交通系統(tǒng)規(guī)劃、設(shè)計及運行管理的核心基礎(chǔ)信息.傳統(tǒng)的OD 矩陣獲取方法是進行大規(guī)模的車輛出行抽樣調(diào)查，為了避開出行調(diào)查技術(shù)困難、經(jīng)費巨大及耗時長的問題，基于交通檢測信息的OD 估計理論一直是研究的熱點.自1978 年Van Zuylen 和Willumsen首次運用極大熵理論通過檢測路段交通量來“反推”O(jiān)D矩陣以來，OD 估計研究已經(jīng)有30余年的歷史.估計方法包括最小二乘法、狀態(tài)空間模型以及信息論方法等［1－3］.

隨著信息檢測技術(shù)的迅速發(fā)展，交通信息采集技術(shù)已經(jīng)迅速從以環(huán)形線圈為代表的“斷面型”檢測技術(shù)向“廣域”的車輛自動識別技術(shù)（automatic vehicle identification，AVI）轉(zhuǎn)變.AVI技術(shù)包括視頻牌照自動識別技術(shù)、基于高精度GPS的車輛定位技術(shù)、路側(cè)信標與車輛通信識別技術(shù)等，其核心是可檢測車輛ID（車牌號）、通過時間以及車輛位置等信息.伴隨著AVI技術(shù)的應用，國內(nèi)外學者對AVI環(huán)境下的OD 估計進行了探索，國外Dixon［4］等通過擴展卡爾曼濾波模型，同時將交通量和AVI之間的行程時間作為觀測變量對車輛OD 進行了估計；Zhou［5］在利用非線性最小二乘估計中考慮了行程時間影響；Barcelo［6］等通過車載藍牙信息預測行程時間，并將之引入卡爾曼濾波進行OD 估計.國內(nèi)Deng［7］等通過整合3S信息資源進行動態(tài)OD 矩陣的研究；陳大山［8］等將神經(jīng)網(wǎng)絡嵌入到卡爾曼濾波中進行動態(tài)OD 矩陣研究.然而上述方法仍然存在如下問題和挑戰(zhàn)：

（1）對于AVI環(huán)境下的動態(tài)OD 估計，目前的方法仍然是通過對經(jīng)典OD 估計模型進行改進，加入AVI檢測時間信息來提高OD 估計精度.事實上，AVI檢測的路徑信息對OD 估計的精度至關(guān)重要.

（2）OD 估計不僅與網(wǎng)絡拓撲、路段流量相關(guān)，還與先驗信息的精度有密切的關(guān)系，以往的研究通常都假設(shè)已獲得可靠的先驗信息來推算OD，未對不同質(zhì)量的先驗信息環(huán)境下估計算法的實用性和可靠性進行精細化分析.

（3）在實際路網(wǎng)中，受投資預算、安裝條件、AVI識別率等多方面的影響，AVI檢測精度和覆蓋率都有一定限制.如何在有限精度／覆蓋的檢測條件下獲得較高精度的動態(tài)OD 矩陣是一個巨大的挑戰(zhàn).

針對上述問題，本文在文獻［9］靜態(tài)OD 估計的基礎(chǔ)上，引入AVI檢測的部分路徑信息、動態(tài)行程時間信息以及檢測器可測性（measurability）判據(jù)，通過改進的粒子濾波算法，從微觀的視角對路網(wǎng)中任意車輛的起迄點以及進入路網(wǎng)的時間進行動態(tài)估計，然后將路段AVI檢測流量作為輔助參數(shù)校正OD 估計結(jié)果，進而獲得較高精度的動態(tài)OD 矩陣.

1 動態(tài)OD估計流程

文獻［9］靜態(tài)OD 估計的方法是以AVI檢測的流量和路徑信息為基礎(chǔ)，基于貝葉斯估計理論修正更新路段－路徑流量關(guān)系，確定可能路徑以及OD.若要進行動態(tài)的OD 估計，還需要解決車輛進入路網(wǎng)的時間判斷、不同時段的動態(tài)行程時間以及車輛路徑估計問題.本文的方法是以車輛路徑信息、車輛到達時間信息為基礎(chǔ)，基于檢測器可測性判斷依據(jù)，通過迭代更新貝葉斯算法來縮小車輛缺失路徑搜索范圍，修正可選路徑的后驗概率；同時在行程時間快速預測方程基礎(chǔ)上確定任意車輛的路網(wǎng)進入時間以及所有缺失路徑的行駛時間；最后通過蒙特卡洛隨機仿真模擬任意車輛的路徑選擇來實現(xiàn)動態(tài)的OD估計.其算法的數(shù)學基礎(chǔ)是粒子濾波思想，即通過尋找一組在狀態(tài)空間傳播的隨機樣本對概率密度函數(shù)進行近似，以樣本均值代替積分運算，從而獲得狀態(tài)最小方差分布的過程［10］.

動態(tài)OD 估計算法描述為如圖1 所示流程，共包含8步驟.

圖1 動態(tài)OD 估計流程圖Fig.1 Flow chart of dynamic OD estimation

Step1根據(jù)AVI采集得到的路徑完整程度，將車輛路徑信息數(shù)據(jù)分成以下3類：第一類是路徑全知型；第二類是路徑半知型，包括已知完整路徑中的兩個或兩個以上斷面，以及兩個斷面之間的完整路徑；第三類是路徑單知型，只知道其出行路徑中的某一個斷面.在考慮AVI檢測精度的基礎(chǔ)上，根據(jù)估計時段內(nèi)AVI檢測的實際捕獲的車輛數(shù)，將先驗信息放樣至與實際捕獲車輛數(shù)相同.

Step2將第h時段中獲得的第一類信息（路徑全知型信息）與AVI先驗信息ψ進行聯(lián)動遞減，即將先驗信息中的OD 對流量減去已知完全路徑的OD信息，獲得新的先驗信息ψ′.當先驗信息中的OD 對流量小于已知的實際路徑流量時，先驗信息中的OD對流量保持為最小流量1.

Step3將檢測獲得的AVI信息中的到達時間依據(jù)AVI檢測器布設(shè)進行分類劃分，獲得任意兩個相鄰的AVI檢測器之間的行程時間，并依據(jù)AVI間行程時間構(gòu)建路段－路徑行程時間快速預測方程.

Step4提取第h時段中AVI信息中的第二類路徑信息（路徑半知型信息），根據(jù)檢測器的可測性判斷依據(jù)和車輛的行程時間估算，用貝葉斯估計對任意車輛的可能缺失路徑進行重判斷，重新確定缺失路徑選擇范圍并修正新選擇路徑后驗選擇概率；最后通過蒙特卡羅隨機仿真對單個車輛進行缺失路徑的估計，進而獲得完整路徑估計數(shù)據(jù).

Step5將第二類信息的路徑估計結(jié)果進行保存并與AVI先驗信息ψ′進行聯(lián)動遞減獲得新的AVI先驗信息ψ″，并基于第三類信息（路徑單知型信息）的捕獲量以及識別誤差，再次進行先驗擴樣信息，更新得到先驗信息ψ?.

Step6第三類路徑信息的處理方法與第二類路徑信息的處理相似，并將此路徑信息進行保存.

Step7 將三類信息進行整理累加獲得第h時段估計的初始修復車輛OD 矩陣.

Step8將初始修復車輛OD 矩陣與AVI檢測器獲得的路段流量進行反向校正驗算，在路段流量的允許誤差范圍內(nèi)修正初始修復OD 對流量.獲得最終第h時段路徑OD 矩陣，并以此作為h＋1 時段OD 先驗矩陣.

1.1 貝葉斯估計算法

貝葉斯估計模型可以在基于不完全情報條件下，對部分未知的狀態(tài)利用主觀概率估計，然后利用貝葉斯公式對發(fā)生概率進行修正，最后利用期望值和修正概率做出最優(yōu)決策.其數(shù)學模型可以描述成如下形式：

式中：Dh表示在第h時段所獲得的實測數(shù)據(jù)；Hi表示第i條路徑；P（Hi／Dh）表示在第h時 段，基于AVI數(shù)據(jù)D與先驗概率后第i條路徑的后驗選擇概率表示在第h時段，選擇第i條路徑的先驗概率；表示在第h時段，基于AVI數(shù)據(jù)可測性判據(jù)對第i條路徑的修正概率；?表示在貝葉斯估計中存在的噪聲.

其中模型的分母為常量，因此模型可以簡化為

1.1.1 行程時間快速判斷方程

動態(tài)的交通出行OD 矩陣估計中，車輛行程時間用于描述車輛進入路網(wǎng)和離開路網(wǎng)的時刻.因此行程時間估計精度將對提高動態(tài)OD 矩陣估計至關(guān)重要.由于AVI檢測的是AVI點對的行程時間，需要將其轉(zhuǎn)化為路段（路徑中的兩個節(jié)點）行程時間.AVI點對與路段行程時間轉(zhuǎn)化數(shù)學方程可描述為

式中：tn1－n2表示任意兩個AVI之間的行程時間；表示第e輛車經(jīng)過任意相鄰AVI檢測器的時間；LAVI1－AVI2表示任意兩個相鄰AVI檢測器之間的距離；Ln1－n2表示兩個AVI檢測器之間所在的兩個節(jié)點之間的距離；β表示行程時間計算中的隨機誤差.

利用所有車輛的平均行程時間預測任一車輛的路段行程時間并不合理，原因在于每個車輛的駕駛行為差異很大，造成車輛之間的行程時間有較大差別.因此為了更加合理地描述車輛的行程時間，本文將估計時段內(nèi)，所有通過任意兩個AVI檢測器的車輛按行程時間區(qū)間進行分類，并分別計算各類行程時間的概率，以P（1），P（2），…，P（n）表示；再通過隨機模擬的方法，確定任一車輛的行程時間，即可滿足車輛個性化的行程時間計算.其數(shù)學模型描述如下：

如圖2所示，假設(shè)編號為A1，A2，A3，…，Aj的車輛分別經(jīng)過檢測器AVI1和檢測器AVI2，并且已知這些車輛到達兩點的時間 集合根據(jù)公式分別計算任意已知車輛經(jīng)過兩個AVI的行程時間，并依據(jù)車輛總數(shù)獲得車輛平均行程時間Tave；將所有車輛的行程時間與平均行程時間進行比較并進行區(qū)間分類，得 到 區(qū) 間 分 類 集 合），其中表示第r組區(qū)間中有輛車，并通過公式（4）計算各組區(qū)間的選擇概率當存在編號為A 的車輛經(jīng)過檢測器AVI1時，可以通過隨機仿真判斷其到達AVI2的行程時間，而該數(shù)據(jù)是檢測器可測性判據(jù)的重要依據(jù).同時通過公式（3）和公式（5）可以計算車輛A 在兩個AVI或兩個AVI間任意兩個節(jié)點的行程時間，進而判斷車輛是否是在本時段內(nèi)進入路網(wǎng)，這部分信息是動態(tài)OD 估計中的一個重要參數(shù).

圖2 AVI點對以及路段節(jié)點對平均行程時間計算示意圖Fig.2 Schematic of calculate average travel time between AVIs and nodes

1.1.2 檢測器可測性判據(jù)

在現(xiàn)實路網(wǎng)中，AVI檢測存在如下兩個事實：①由于AVI檢測器存在車輛捕獲誤差和識別誤差，因此其檢測精度有限；②由于資金和安裝條件等約束，AVI 無法全覆蓋.據(jù)此提出AVI 可測性（measurability）概念，定義為：基于AVI檢測信息數(shù)據(jù)缺失的特點和已檢測車輛的動態(tài)行程時間信息，分析判斷車輛經(jīng)過某一路徑的可能性.

如圖2所示，假設(shè)估計時段內(nèi)，已知編號為A 的車輛通過AVI系統(tǒng)獲得其部分路徑信息，即車輛A通過檢測器AVI1和AVI2，并且已知經(jīng)過兩檢測點的到達時刻分別為T1和T2，通過式（3）計算得到AVI2至AVI3的行程時間為T2－3，并已知車輛A 的缺失 路 徑 概 率 集 合 為ξ（P（O1D3），P（O1D4），P（O1D5），P（O1D6）），根據(jù)AVI可測性定義，當車輛A 有充足時間到達下一個路段檢測器AVI3，但車輛A 的路徑信息卻不包含檢測器AVI3時，存在兩種情況：①車輛A 通過AVI3，但由于檢測誤差而沒有被記錄；②車輛沒有經(jīng)過布設(shè)AVI3的路段，直接到達D3.基于以上兩種情況，分別提出判別模型如下：

如果（T2＋T2－3＜h＋Δt）則

式中：Δt表示估計時段長度；P′（OmDn）表示修正后的缺失路徑的車輛選擇概率；P（OmDn）表示初始的缺失路徑的車輛選擇概率；ε表示AVI設(shè)施的檢測誤差，包含捕獲誤差與識別誤差.

AVI檢測器可測性判別依據(jù)對于車輛路徑的上游路段是同樣適用的，將檢測器可測性判據(jù)根據(jù)行駛路徑的下游缺失路徑和上游缺失路徑分別描述為如下兩種形式：

式（7），（8）中：P（u）表示下游行駛方向中第u組缺失路徑的先驗信息；θ（u）表示上游行駛方向中第u組缺失路徑的先驗信息；εju表示在第u組缺失路徑下游行駛方向，基于j組滿足可測條件的AVI檢測器的可測性修正概率；εku表示在第u組缺失路徑上游行駛方向，基于k組滿足可測條件的AVI檢測器的可測性修正概率；Ti表示在車輛i最后一次被檢測到的時刻；T′i表示在車輛i第一次被檢測到的時刻；h表示估計時段的開始時刻；P′（u）表示下游路徑中第u組路徑的后驗概率；θ′（u）表示上游路徑中第u組路徑的先驗概率.

1.1.3 貝葉斯路徑概率模型

將先驗信息矩陣，以及1.1.1節(jié)的行程時間快速預測方程，1.1.2節(jié)的檢測器可測性判據(jù)作為先驗路徑選擇概率，貝葉斯車輛路徑概率算法模型更新如下：

式中：Pu（Hu／Dh）′表示在第h時段第u條缺失路徑修正后的路徑選擇概率；Pu（Hu）h表示在第h時段第u條缺失路徑的先驗選擇概率；αu表示路徑選擇概率的隨機誤差.

1.2 蒙特卡羅仿真路徑選擇過程

蒙特卡羅方法由于能夠真實地模擬實際物理過程，逼近真實的物理結(jié)果，因此通過蒙特卡羅方法對AVI檢測的樣本粒子模擬基于概率的空間狀態(tài)選擇.

對任意被AVI檢測器所檢測的車輛，其所有的缺失路徑修正概率集合為ξ（P1，P2，…，Pu），并且根據(jù)行程時間快速預測方程計算的各缺失路徑起點的出發(fā)時間集合為ζ（T1，T2，…，Tu），設(shè)定

式中：Wi表示路徑后驗概率的累積頻率；Pi表示第i組路徑后驗概率.

仍以圖2為例，通過不斷產(chǎn)生隨機數(shù)的方法，模擬任一車輛可能的選擇路徑.若隨機數(shù)S小于Wi大于Wi－1，則說明車輛的路徑為第i條缺失路徑.完成路徑選擇后，相應工作是車輛進入路網(wǎng)時間的確定，根據(jù)駕駛員的最后路徑選擇，判斷：

式中：Vh表示第h時段內(nèi)實際進入路網(wǎng)的車輛集合.

1.3 路段流量修正

AVI檢測器同樣可以獲得路段流量數(shù)據(jù)，路段流量與OD 對之間存在如下關(guān)系：

式中：Fi表示路段i上的路段流量；fmn表示從起點m到終點n的OD 量；pimn表示從起點m到終點n的OD 量中選擇路徑i的概率.

通過貝葉斯估計和蒙特卡洛隨機仿真，可以計算所有被檢測車輛的后驗路徑選擇概率，并使車輛作為檢測樣本粒子擬合后驗路徑選擇概率密度獲得后驗OD 矩陣.本文采用動態(tài)的隨機規(guī)劃方法，當路段－路徑流量相對誤差小于5%，可視為收斂.其數(shù)學過程描述如下：

式中：α1表示檢測誤差下限；α2表示檢測誤差上限；fimn表示通過第i條路徑的后驗矩陣中的起點為m，終點為n的OD 對流量.

2 案例分析

2.1 研究區(qū)域及實驗平臺

2008年始，上海市開始在快速路系統(tǒng)規(guī)劃布置視頻牌照自動檢測設(shè)備，本案例以上海市南北高架快速路作為研究對象.在VISSIM 仿真模型中設(shè)置“虛擬檢測器”模擬實際視頻的檢測參數(shù)及檢測精度，利用虛擬的檢測數(shù)據(jù)進行OD 估計并和VISSIM中的“真實”O(jiān)D 做比較，由此來解決OD 估計精度的驗證問題［11］.

2.2 動態(tài)OD 估計結(jié)果分析

2.2.1 OD 估計精度分析

基于虛擬的“AVI”檢測數(shù)據(jù)，得到不同估計時段OD 估計相對誤差值，如圖3，4所示.當AVI覆蓋率為60%，先驗OD 精度為40%時，其相對誤差僅為21.8%；而當AVI覆蓋率為80%，先驗精度為80%時，相對誤差精度提高為17.3%.另外當覆蓋率為50%的情況下，動態(tài)OD 估計的精度出現(xiàn)較大的波動，但其平均的動態(tài)OD 估計相對誤差為30.11%；當AVI覆蓋率為50%，先驗精度60%時，相對平均誤差為28.87%.文獻［8］通過對簡化快速路的研究，在僅考慮較高AVI覆蓋率和高精度的先驗信息的前提下，其動態(tài)OD 估計相對誤差為21.04%，與本文實際路網(wǎng)動態(tài)OD 估計精度相當［8］.但當AVI覆蓋率為45%，OD 估計精度迅速衰減，其平均相對誤差降為41.9%，已經(jīng)較難滿足實際需求.主要原因是當覆蓋率較低時，主要的OD 對流量已經(jīng)不能被檢測到.

2.2.2 對先驗信息的敏感性分析

圖3 80%覆蓋率下相對誤差Fig.3 Relative error in 80%coverage

圖4 60%覆蓋率下相對誤差Fig.4 Relative error in 60%coverage

以往的研究對先驗信息有較大的依賴，而本文引入的檢測器可測性方法則對先驗信息有較強的魯棒性.以較低的AVI覆蓋率（60%）為例，在不同的先驗信息精度下，OD 估計精度變化如圖5所示.從圖5可以看出，在任意先驗信息精度下，本文方法均能較好地估算動態(tài)OD.

圖5 60%覆蓋率下相對誤差變化Fig.5 Relative error changes in 60%coverage

3 結(jié) 論

本文針對當前視頻牌照自動識別系統(tǒng)在我國逐步應用的現(xiàn)實需求，提出了基于AVI信息的動態(tài)OD 估計新方法，并以上海市南北高架快速路網(wǎng)模型對方法進行檢驗.主要結(jié)論如下：

（1）本方法突破了傳統(tǒng)方法中僅依賴路段流量和行程時間作為OD 估計的主要信息的局限，基于粒子濾波的思想提出了基于路徑信息、檢測器可測性、路徑行程時間的OD 估計方法.上海市南北高架的案例分析表明，當AVI覆蓋率為50%，先驗精度為60%時，平均相對誤差為28.87%，仍能基本滿足動態(tài)交通管理的需求.

（2）針對AVI檢測誤差和覆蓋度所提出的檢測器可測性判據(jù)，較好地緩解了對先驗信息的依賴，可以在隨機的先驗信息下推算出較高精度的OD.

（3）此方法對大規(guī)模開放式地面道路網(wǎng)絡同樣具備適用性.通過AVI信息的采集以及交通分配算法得到路段－路徑流量數(shù)據(jù)庫，應用本文所提出的方法，能迅速縮小任意車輛路徑選擇范圍，提高車輛真實路徑的選擇概率，并可通過AVI檢測的流量數(shù)據(jù)對OD 矩陣再次進行修正，實現(xiàn)快速準確的對大規(guī)模開放式路網(wǎng)OD 進行估計.

（4）本文提出的方法提供了一種從微觀的車輛路徑修復角度進行動態(tài)OD 估計的方法.在進一步的研究中，將重點考慮如何結(jié)合本方法優(yōu)化AVI設(shè)施的布局方案，在規(guī)劃層面，指導更科學地布設(shè)，從源頭上為運營管理中獲取高精度動態(tài)OD 提供支撐.

［1］ 常云濤.考慮交通流行駛時間的高速公路動態(tài)OD 矩陣估計模型［J］.同濟大學學報：自然科學版，2009，37（9）：1184.CHANG Yuntao.Dynamic OD matrix estimation model of freeway with consideration of travel time［J］.Journal of Tongji University：Natural Science，2009，37（9）：1184.

［2］ Tsekeris T， Stathopoulos A. Real－time dynamic origindestination matrix adjustment with simulated and actual link flows in urban networks［J］.Transportation Research Record 1857，2003：117.

［3］ Bierlaire M，Crittin F.An efficient algorithm for real－time estimation and prediction of dynamic OD tables［J］.Operation Research.2004，52（1）：116.

［4］ Dixon P M，Rilett R L.Population origin destination estimation using automatic vehicle identification and volume data ［J］.Journal of Transportation Engineering，2005，131（2）：75.

［5］ Zhou X S， Mahmassani H S. Dynamic origin－destination demand estimation using automatic vehicle identification data［J］.IEEE Transactions on Intelligent Transportation Systems，2006，17（1）：105.

［6］ BarcelóJ，Montero L，Marqués L，et al.Carmona，travel time forecasting and dynamic OD estimation in freeways based on bluetooth traffic monitoring［J］.Transportation Research Record 2175，2010：19.

［7］ DENG Bin，HU Jianming，DAI Zhiguang，et al.Using 3S technology to estimate OD matrix［C］／／Proceedings of the 12th International IEEE Conferenceon Intelligent Transportation Systems.St.Louis：IEEE，2009：836－842.

［8］ CHEN Dashan，SUN Jian，LI Keping.Dynamic OD estimation under automated vehicle identification environment［J］.Key Engineering Materials，2011，46（7）：835.

［9］ 孫劍，馮羽.車輛自動識別環(huán)境下的OD 估計新方法研究［J］.同濟大學學報：自然科學版，2011，39（12）：1800.SUN Jian，F(xiàn)ENG Yu.A new method of OD estimation based on automatic vehicle identification Data［J］.Journal of Tongji University：Natural Science，2011，39（12）：1800.

［10］ 胡士強，敬忠良.粒子濾波算法綜述［J］.控制與決策，2005，20（4）：361.HU Shiqiang，JING Zhongliang.Overview of particle filter algorithm［J］.Control and Decision，2005，20（4）：361.

［11］ 孫劍，李克平，魏靜，等.基于視頻牌照識別的動態(tài)交通OD估計仿真優(yōu)化［J］.公路交通科技，2009，26（8）：130.SUN Jian，LI Keping， WEI Jing，et al.Dynamic OD estimation simulation optimization based on video license plate recognition［J］.Journal of Highway and Transportation Research and Development.2009，26（8）：130.