智慧型停車場系統-車輛追蹤子系統

全文

(1)第一章簡介 1.1 研究動機 1.1.1 現有的停車問題隨著時代的演進，機動車輛已逐漸成為一般民眾的必需品之一，而民眾所擁有的機動車輛也與日俱增。民國 97 年交通部統計處提供的統計數據，每年機動車輛登記數不斷的向上攀升（表 1.1），此一現象也造成許多交通上的問題，例如：交通壅塞、交通事故、停車糾紛…。表 1.1 全台機動車輛登記數統計年代. 該年車輛登記數(輛該年車輛登記數輛). 2008. 21,068,657. 2007. 20,711,754. 2006. 20,307,197. 2005. 19,862,807. 2004. 19,183,136. 而在都會區中，停車問題尤其嚴重，以現有停車場所提供的停車格數量早已供不應求，以台北市為例，由台北市交通局提供的數據顯示，目前台北市現有停車位共有 1,001,826 個，其中包括機車停車位 463,223 個以及汽車停車位 538,603 個（表 1.2），而台北市現有機動車輛總數達 1,801,245 輛，其中包含機車 1,080,091 輛以及其他各式車種 721,154 輛（表 1.3），從數據中可看出台北市汽車約有 25%. 1.

(2) 沒有停車位，而機車更是缺少了 57%的停車位。因此要如何解決停車問題，以避免車輛因尋找停車位而浪費時間，甚至造成交通阻塞，或是不必要的交通事故及交通問題，已是刻不容緩。表 1.2 停車位數量統計時間. 汽車. 機車. 總計. 97/08. 531,666. 451,262. 982,928. 97/09. 534,975. 457,113. 992,088. 97/10. 538,603. 463,223. 1,001,826. 表 1.3 機動車輛數量統計時間. 其他各車種. 機車. 總計. 97/08. 723,101. 1,076,342. 1,799,443. 97/09. 722,224. 1,079,044. 1,801,268. 97/10. 721,154. 1,080,091. 1,801,245. 為了解決停車問題，各地政府也在道路兩旁提供了為數不少的路邊停車格以供停車。但在都會區由於人口密集，且土地有限，以台北市為例，都市人口密集度達 9,666 人/每平方公里，而高雄市人口密度更是高達 9,827 人/每平方公里，因此無法提供足夠的路邊停車格。所以在都會區要解決停車問題唯有在大型公有停. 2.

(3) 車場上著手，使每位車主皆能充分而有效的使用這些停車場，以徹底解決停車問題。因此政府單位也利用網路提供了全台公有停車場位置的線上查詢，亦可提供目前該停車場現有停車空位數目，以使每位車主皆能充分掌握各停車場的位置及現況，以避免因盲目尋找停車場所衍生出來的交通阻塞、及交通事故等問題。. 1.1.2 研究目的目前現有的停車場所提供的停車場管理系統（Intelligent Parking System, IPS），最常見的有空位統計系統與監視系統（Surveillance System）（圖 1.1）。其中監視系統主要目的是提供監視人員能即時的看到停車場內各種異常狀況的資訊（例如：人為蓄意縱火、車輛發生擦撞…等），並隨時加以錄影以便日後有任何事故或糾紛可以做為佐證。而空位統計系統則是用來顯示各停車場內尚存的停車位（圖 1.2（a）），或是停車場內各停車區域尚存的停車空位（圖 1.2（b）），使車主能很快的找到可停車的停車場，避免車主進入停車場內盲目的尋找車位。. 圖 1.1 監視系統 3.

(4) （a）. （b）圖 1.2 空位統計系統. 註 1：出自 Google 圖片搜尋。監視系統與停車空位顯示，可說是目前停車場的基本配備，但此兩系統卻是各自運作，無法做一整合，以致所能提供給車主的服務是非常有限的。在監視系統方面，它僅能提供管理者一個監視的畫面，管理者必須緊盯著每一個螢幕才能確實的掌握目前停車場的即時狀態，這對管理人員來是個嚴酷地考驗，且無法給予車主任何即時的幫助。而停車空位顯示，雖然可以讓車主知道那一個區域有空的停車位，但前提是車主必須把車開到該停車區域附近才看得到該顯示，若該停車區域已經停滿車輛，則車主必須再開到另一個區域的入口，才知道該區域是否可以停車，因此停車空位的顯示對車主的幫助實在是有限的。在停車的尖峰時刻，車主難免還是要盲目的開著車在停車場內尋找停車位。本研究的目的就是為了解決停車空位顯示及監視系統的缺失，利用影像處理 4.

(5) 的技術將這兩個系統做一個有效的整合，使得這兩個系統能發揮最大的效益，並提供車主最大的服務。首先在監視系統方面，將加入車輛追蹤的功能，對停車場內的所有車輛做追蹤，以確實掌握停車場內每一輛車的動態，並將所有車輛的動態整合在一個停車場的平面圖中(如圖 1.3)，使得管理人員能夠完全的掌握停車場內所有的車輛現況。而車輛追蹤的目的除了是給管理者一個整合的停車場車輛動態的監視畫面外，另外可結合車位的統計，根據各區域車位剩餘的狀況，來啟動車輛前方的方向指示器(如圖 1.4)，即時的給予每一輛進入停車場的車輛一個停車方向的指示，以指示該車輛就近往有各車位的停車區域前進，來停放車輛，如此可以避免車輛盲目的尋找停車位。. 圖 1.3 實驗停車場平面圖. 5.

(6) 圖 1.4 方向指示器. 1.1.3 監視攝影機的選擇目前對於追蹤常用的固定式攝影機有 PTZ 攝影機（Pan, Tilt, Zoom Camera）（圖 1.5）以及環場攝影機（Omni-direction Camera）（圖 1.6）。兩種攝影機在功能上各有其優缺點，對 PTZ 攝影機而言，可以根據使用者需求不同改變攝影機焦距，拍攝角度，而得到較為清晰的影像，但是相對之下，拍攝的場景也較為狹隘受到侷限，拍攝時可能會產生拍攝死角；相對的，環場攝影機所拍攝到的場景，可以涵蓋將近 360 度的範圍（圖 1.7），同時可以避免像 PTZ 攝影機所產生的死角問題，但由於該攝影機焦距為固定，無法因使用者需求不同而有所改變，拍攝出的影像也相對的較為模糊。. 圖 1.5 PTZ 攝影機 6.

(7) 圖 1.6 環場攝影機. 圖 1.7 環場影像在本研究中採用的是環場攝影機，原因有以下三點： 1. 本研究主要只針對停車場中的車輛移動去做追蹤的動作，並不會對車輛的. 7.

(8) 細部特徵去做分析，即便拍攝到更清晰的畫面，對於實際研究上並不會提供更多的助益； 2. 考量到安裝攝影機的成本，由於 PTZ 攝影機拍攝視角較為狹隘，所以在安裝上所耗費的數量會遠比環場攝影機多上許多，所造成的花費是很可觀的； 3. 本研究最終目的是希望可對多輛的車子同時進行追蹤，而 PTZ 攝影機主要是應用在對單一物體進行移動追蹤，然而在停車場中不可能永遠只對單一的車輛進行追蹤，而是要針對全部停車場中所有移動的車輛同時進行追蹤，因此在此先天條件之下，要使用 PTZ 攝影機對多台車輛進行追蹤就會出現問題；假若為了對單一特定事件進行追蹤（例如：要鎖定某個可疑人士進行追蹤），則可使用環場攝影機再搭配少數 PTZ 攝影機，即可以最少的成本達到最佳的追蹤成果。. 1.2 文獻探討為了得到停車場中各車輛的行車資訊，有不少研究是在停車場的各角落分別架設了 PTZ 監視攝影機，從拍攝到的停車場畫面中擷取出移動車輛的部份，再對偵測出的車輛進行追蹤。在這樣的情況下進行追蹤的結果，會容易出現覆蓋（Occlusion）的情況，為了解決這個問題，在[Bla 02]中，利用在停車場架設多台攝影機，從不同的角度對同一位置的車輛進行追蹤，藉由攝影機之間的相對位置，來得到出現在畫面中車輛的正確位置，如此以來可以避免從某個角度拍攝時產生的覆蓋問題。此方法雖然可以解決覆蓋的問題，卻會造成在架設時的成本提. 8.

(9) 高許多，且在追蹤時的計算量會過於龐大，無法在實際的停車場中進行更多複雜的即時追蹤。在[Wan 06]中採用了在停車場上方架設環場攝影機來解決這項困擾，利用這種俯視（Top-view）的監測方式，可以大幅的降低覆蓋問題的出現。但是環場攝影機拍攝出的畫面通常畫質並不高，在辨識度上會出現問題，且拍攝的角度是由上而下，若要應用在安全監控上，只由環場攝影機所提供的資訊仍顯不足，在[王 05]中提出同時利用 PTZ 攝影機及環場攝影機來進行安全監控的動作。由於環場攝影機在硬體結構設計上的影響，會造成環場攝影機拍攝出的畫面與真實的場景產生扭曲的現象，此一現象會對日後停車場中車輛的追蹤造成極大的誤差，在 [Wan 06]中，對於環場攝影機拍攝出的扭曲畫面有提出解決的方法。本論文根據 [Wan 06]的系統架構，更進一步修正追蹤的準確度，並提出更多的追蹤方法提高追蹤的準確度。在停車場追蹤系統的架構主要有兩層，分別是車輛偵測以及車輛追蹤兩個部份。在車輛偵測上最重要的動作就是前景物的擷取，現在常見的前景物擷取方法有背景相減法（Background Subtraction）[Als 04] [Kim 08] [Unn 07]、連續影像相減法（Temporal Difference Method）[Lip 98] [Kuh 06] [Wu 06]。背景相減法要先從拍攝影像畫面中建立出不含任何前景物資訊的影像作為背景影像，接著將一連串輸入的影像與背景影像相比較，即可得到前景物的資訊。連續影像相減法，使用的方法是利用前後時間差的連續影像（Frame）互相比較，得到移動物體的資訊，. 9.

(10) 此一方法對於環境改變的適應性較佳。由於一般的連續影像相減法所得到的前景物資訊會有殘缺不全的現象，[羅 08]中提出使用水平集方法（Level Set Method）的技術來解決這個問題。近年來有不少技術被應用在追蹤與預測上，例如：卡曼濾波器（Kalman Filter） [Mas 01][姚 08]、Mean Shift 以及粒子濾波器（Particle Filter）。簡單來說，卡曼濾波器是一個”最優化回歸數據處理演算法”（Optimal Recursive Data Processing Algorithm），對於解決很多大部分簡單的問題，效率都是各方法中最高的，但是在[林 02]中提到卡曼濾波器是基於高斯分佈來預測穩定移動的物體，對於像車輛這種可能突然有加速、減速或改變方向的情況，並不適用。Mean shift 在追蹤上雖然執行速度上沒有卡曼濾波器快，但是在穩定度上卻比卡曼濾波器來的穩定，而粒子濾波器對於移動較快或是運動變化劇烈的物體，可以進行複雜的追蹤，與前兩者相比，得到欲追蹤物的機率較大。 Mean Shift[Com 99][Che 08][張 05]的追蹤方法是透過不斷的計算 Mean Shift Vector，使 Vector 逐漸往 Bhattacharyya 係數區域的最大值靠近，由於在計算的區域範圍內所得到的最大值只有一個，因此利用 Mean Shift 進行追蹤所得到的結果也比較準確、穩定；但因為 Mean Shift 是用上個時間點所得到的目標物位置作為起始位置，來預測現在時間點目標物的位置，再利用 Mean Shift Vector 會逐漸往區域最大值靠近的特性，求出的結果會落在距離起始位置不遠的區域最大值處；且因為物體的移動會具有時間及空間連續性，所以在下一個時間點所得到的目標. 10.

(11) 物位置，必定會落在目前時間點目標物的位置附近。然而，若移動物體的移動速度太快，在下一個時間點時移動已超過了區域範圍，那麼在計算區域最大值時，便會出現求得的結果並非實際目標物的情形，亦即宣告追蹤的失敗。粒子濾波器[Boe 07][Ekm 07][張 05]也可以稱作是”序列蒙地卡羅法” （Sequential Monte Carlo Methods, SMC），是利用在影像上隨機選取粒子的方式來預測下一個時間的位置，用前一個時間點目標物的位置作為起始位置，以起始位置為中心向外隨機的分散粒子，透過不同的粒子具有不同的權重，不斷的去預估粒子權重的特性，來達到最終追蹤的目的。雖然起始位置附近所散落的點數較多，但由於是隨機分散的粒子，也有可能會散落在距離起始位置較遠處，因此對於移動速度較快的物體而言，與 Mean Shift 相比較之下會有較高的機率可以追蹤的到，但因為是採用隨機取樣的方式，也造成了在預測的過程中必須帶有運氣的成分，追蹤結果也相對的較不穩定，雖然最後可以利用散落更多粒子的方法來提高追蹤到物體的準確性，相對受到影響的就是執行的速度也會因此變慢；與 Mean Shift 相比，不論在準確性、穩定性比較，還是執行的速度上都還是有些許的落差。根據不同的研究，粒子濾波器也發展出許多不同的變形，例如：擴展卡曼粒子濾波器（Extended Kalman Particle Filter, EKPF）、Unscented 粒子濾波器（Unscented Particle Filter, UPF）、混合卡曼粒子濾波器（Mixed Kalman Particle Filter, MKPF）、 Rao-Blackwellised 粒子濾波器（Rao-Blackwellised Particle Filters, RBPF）。綜合各種技術的分析，在偵測前景物方面本研究採用連續影像相減法，目的. 11.

(12) 是考量到停車場中隨時會有停放的車輛，若採用背景相減法會造成不必要的偵測結果，在擷取到前景物之後，在使用 level set method 來填補擷取到物體的殘缺區域，接著用連接組件（Connected Component）得到各個前景物的資訊。到了追蹤階段，參考[Sha 04][Xu 08]分別對手部追蹤以及臉部追蹤提出將兩種技術結合使用的想法，同時擷取出兩種技術的優點，利用 Mean Shift 對物體進行追蹤時具有較高的準確度，粒子濾波器對物體追蹤時具有較高的成功率，彼此彌補對方技術上的缺點，透過這項技術達到本研究最終追蹤的目的。. 1.3 文章架構本論文在第二章中，對本研究的系統架構和執行環境做更詳細的解釋，以及追蹤系統的運作流程圖。第三章接著說明環場攝影機拍攝影像，空間與成像平面轉換的原理，以及該原理應用在本研究上的校正動作（Calibration）。在第四章則會先介紹什麼是連續影像相減法以及其優點，並且介紹如何對停車場中移動的車輛擷取出資訊。第五章將針對本研究中所採用的追蹤方法，一般粒子濾波器（Traditional Particle Filter）以及聯合粒子濾波器（Joint Particle Filter）兩種方法分別做詳細的介紹及比較。第六章為實驗數據及追蹤得到的結果。最後第七章會提出本研究過程中，遇到的困難以及未來研究要繼續努力的方向。. 12.

(13) 第二章系統架構 2.1 架設環境一般來說，監控裝置的執行環境可分為兩種，室內場景及室外場景。室外場景受到許多因素造成的影響也遠比室內場還要嚴重複雜，這些因素包含： ‧ 影像明亮度（例如：晴天、陰天、雨天…） ‧ 光線色彩變化（例如：清晨、中午、黃昏、夜晚） ‧ 光線因空氣懸浮粒子產生的折射現象（例如：起霧、下雪、煙…） ‧ 場景的背景因物體移動或不安定因素造成的晃動（例如：物體產生的陰影、樹葉受到風的吹動產生的晃動…） ‧ 因光源變換造成的現象（例如：固定架設在某處的攝影機，因太陽東起西落的自然現象影響，產生拍攝時有面光或背光的問題出現）與室外場景這些不安定的因素相比，室內場景幾乎不存在這些因素，以理論來看似乎是處於一個拍攝穩定狀態，但實際上卻也具有不少的潛在問題： ‧燈光出現閃爍情形（例如：照明光源本身會閃爍、照明光源發生故障造成閃爍…） ‧ 場景中同時有多個照明光源彼此互相干擾 ‧ 室外造成的干擾（室外的光線透過窗戶造成的變化） ‧ 週遭環境變化產生的影響（例如：人的走動、開關門…） ‧ 不同照明光源對同一物體會產生不同的影響（例如：白幟燈、鹵素燈、. 13.

(14) 日光燈、燈泡…等） ‧ 背景環境的影響（背景表面光滑程度造成的折射，例如：金屬表面造成的反光現象遠比毛布料表面造成的反光現象還要來的嚴重）本研究使用的實驗環境為室內實驗停車場（圖 2.1），停車場上方裝置有共八台環場攝影機，本研究主要是在將環場攝影機架設於同一水平面下的情況下來做研究討論，即便在實際停車場場景會因建築設計的問題可能造成無法架設於同一水平面，本研究的理論架構依舊是可順利執行的；在硬體設備上提供兩台主機做為監視平台（圖 2.2），每台主機內皆裝置了監視系統及追蹤系統，各負責四支環場攝影機的監控，接著將監控的畫面經過主機處理後所得到的資訊，統一彙整至使用者操作平台，將使用者所需資訊提供給使用者，整個系統示意圖，如圖 2.3 所示。. 圖 2.1 室內實驗停車場. 14.

(15) 圖 2.2 監視平台. 圖 2.3 系統示意圖由圖 2.1 可觀察到本研究是將環場攝影機架設於停車場上方，會如此架設環場攝影機的原因，主要是為了將環場攝影機的優點徹底的運用，同時考量到要有. 15.

(16) 效利用環場攝影機的可拍攝範圍。環場攝影機的特性之一就是可以拍攝 360 度的場景，假若裝置環場攝影機時，鏡頭與地面呈現一個非九十度的夾角，那麼拍攝出的影像畫面會有過多不必要的拍攝資訊（例如：天花板也會被拍攝到）。再者，停車場進行監控的對象，主要不外乎是車子與行人，這些欲偵測物體產生移動行為的空間也僅限於地平面，即便拍攝到再多的非地面的部份，亦無法提供任何有效的資訊給系統使用，為了降低這種情況造成環場攝影機無法有效的發揮優點，才決定將環場攝影機以垂直地面的方式裝置。為了提高環場攝影機的有效使用，同時也提高追蹤的正確性，必須先針對停車場中裝置的環場攝影機建立各種相關性。首先必要條件必須先得到裝置於停車場中各環場攝影機的空間座標位置（X, Y, Z）、各攝影機的代表編號（num）以及各攝影機可拍攝的有效拍攝半徑（R），其中 R 是實際拍攝空間的拍攝半徑，而非影像畫面上的拍攝半徑。利用前兩個已知條件，可以得到各環場攝影機之間的方位關係矩陣（directional matrix）（表 2.1）表 2.1 相鄰關係矩陣 num. 0. 1. 2. 3. 4. 5. 6. 7. 0. 0. 1. 1. 1. 0. 0. 0. 0. 1. 1. 0. 1. 1. 0. 0. 0. 0. 2. 1. 1. 0. 1. 1. 1. 0. 0. 3. 1. 1. 1. 0. 1. 1. 0. 0. 4. 0. 0. 1. 1. 0. 1. 1. 1. 5. 0. 0. 1. 1. 1. 0. 1. 1. 6. 0. 0. 0. 0. 1. 1. 0. 1. 7. 0. 0. 0. 0. 1. 1. 1. 0. 16.

(17) 利用 0~7 分別代表八個不同的方位，依序是東、東北、北、西北、西、西南、南、東南（表 2.2）。表 2.2 方位代碼代號. 方位角度. 0. -22.5° ~ 22.5°. 1. 22.5° ~ 67.5°. 2. 67.5° ~ 112.5°. 3. 112.5° ~ 157.5°. 4. 157.5° ~ 180.0° -180.0° ~ -157.5°. 5. -157.5° ~-112.5. 6. -112.5° ~-67.5. 7. -67.5° ~-22.5. 而攝影機的有效拍攝半徑，須經由的反扭曲（Dewarpping）計算的公式，以得到區域座標系統與全域座標系統之間的對應關係，此對應關係的計算另在 3.2 節中會有詳細的解釋。在得到將影像座標（區域座標系統）對應到場景座標（全域座標系統）的結果後，再利用（2.1）即可求得有效拍攝半徑。. R = Qx 2 + Qy 2. （2.1）. Qx 和 Qy 分別代表全域座標系統的 x、y 座摽值。接著由第一個已知條件可以求出空間中各個環場攝影機彼此之間的距離（ Dij ），透過（2.2）便可以得到各環場攝影機之間的相鄰關係矩陣（Adjacent matrix），此為一對稱矩陣（Symmetric Matrix）（表 2.3） 17.

(18) 2 Dij − R > ε  2 Dij − R < ε. rel = 1 rel = 0. （2.2）. ε 為一門檻值（Threshold），rel 表示各攝影機之間的相關性，1 表示彼此是具有方位相關性，0 則表示彼此不相關，最後將所得到的兩個矩陣合併，可得到關係矩陣（表 2.4）。關係矩陣中，有相關的環場攝影機彼此之間會出現拍攝範圍疊覆（Overlap）情形（圖 2.4），這種疊覆情形是有意義的（Significant），因為當追蹤的物體超出環場攝影機的有效拍攝範圍時，可透過追蹤時得到的資訊（例如：位置、速度…）及關係矩陣來判斷該由何台環場攝影機接手繼續進行監控。有效拍攝範圍內發生疊覆的部份所得到的資訊，將這些資訊傳遞給使用者操作平台時，可能會同時產生兩個或多個不同的結果，造成使用者在操作時的困擾，因此本研究在最後這個步驟必須確保結果的唯一性。. （a）. （b）圖 2.4 具相鄰關係的影像。. （a）（b）兩圖是左右具有相鄰關係的攝影機所拍攝出的影像，（a）圖右半部與（b）圖左半部的拍攝範圍會出現疊覆的現象. 18.

(19) 表 2.3 方位關係矩陣 num 0 1 2 3 4 5 6 7 6 0 7 0 0 0 0. 0. 1 0 0 0 0 0. 1. 2. 2. 4 5. 3. 3 4 2. 4. 4 4 4 5. 5. 4 4 3 4 2. 6. 4 4 4 4 4 5. 7. 4 4 4 4 3 4 2. 6 0 7 0 0 1 0 0 0 6 0 7 1 0 6. 表 2.4 關係矩陣 num 0 1 2 3 4 5 6 7 6 0 7. 0. 1 0. 1. 2. 2. 4 5. 3. 3 4 2. 6 0 7 1 0. 4. 4 5. 6 0 7. 5. 3 4 2. 1 0. 6. 4 5. 7. 3 4 2. 6. 2.2 系統流程圖本研究主要的系統架構，著重於在監視平台中所做的處理過程，至於圖 2.3 中的第一層以及最後一層，只是單純的影像輸入及資料輸出，在研究上並沒有做任何的運算或處理，因此在本論文中不加以討論。. 19.

(20) Image Sequences. Vehicles Detection. Vehicles Tracking. Coordinate Transformation. Floor-plane Mapping. Output. 圖2.5 系統流程圖圖 2.5 為系統流程圖，首先將停車場中所裝置的環場攝影機拍攝到的畫面輸入，由於攝影機可能會出現偏差，因此必須先對輸入的影像進行修正，這個部份. 20.

(21) 會在 3.1 節中提到。接著對輸入的影像進行情景物的擷取，在這邊因考量到停車場的變化性很大，因此使用連續影像相減法來取得前景物，前景物擷取的方法會在下一節進行討論。假若沒有取得前景物資訊，就不會進入下一個步驟；相反的，若有找到前景物，則對已經取得前景物的影像進行反扭曲（Dewarpping）的動作，這個動作在 3.2 節中會有詳細的介紹。研究過程中，在這個步驟會有出現兩種爭議情況，一種情況如上所述，另一種情況則是執行過程對調，也就是先對影像進行反扭曲的動作，再對反扭曲後的影像進行前景物的擷取。選擇採用前者的原因，主要是考量到系統執行的即時性以及對前景物擷取的意義性，由於進行反扭曲的動作所需耗費的運算量極大，因此在多台攝影機同時要進行這個動作時，往往會造成嚴重的影像流失及主機負擔，無法有效做出即時的追蹤；再者，透過以上兩種不同順序分別取得的前景物資訊（圖 2.6），很明顯的可以觀察到，後者所得到的資訊遠比前者所得到的資訊要少，對接下來的追蹤也會產生較大的影響。. （a）. （b）圖2.6 前景物資訊。 21.

(22) （a）先做前景物偵測後，再做反扭曲（b）先做反扭曲之後，再做前景物偵測接下來的步驟就是要進行目標物的追蹤，本研究提出兩種追蹤方法來進行討論，分別是一般粒子濾波器以及聯合粒子濾波器。本研究中有嘗試使用Mean Shift 進行追蹤，與粒子濾波器相比容易出現追蹤失敗的情形，研判出現錯誤的條件有二： ‧當目標物與背景出現顏色相近的時候，會發生追蹤失敗。由於Mean Shift 是根據這一個時間偵測到的目標物RGB Histogram值為基準，到下一個時間的影像中尋找相近的RGB Histogram值做為預測下一個時間點位置的方法，當出現背景與目標物顏色相近時，就無法判定出實際的目標物到底在影像中的何處，造成在追蹤時可能找到非目標物的結果。 ‧物體移動速度過快的時候，前一個段落中有提到，當Mean Shift在追蹤時遇到欲偵測物移動速度較快的情況，可能會在下一個時間點超出附近區域最大值範圍，此時便會出現追蹤失敗的訊息。上述執行過程只能針對單一環場攝影機拍攝的影像進行處理，當移動物體超出了環場攝影機的拍攝範圍時，就必須利用前一節中所提到的關係矩陣來做環場攝影機的切換，以達到時間與空間連續性的追蹤，最後再將各支環場攝影機每次經過追蹤過後的結果統一顯示在使用者的監視平台上，即完成整套流程。. 22.

(23) 第三章環境地圖建立本研究的最終目的是希望將攝影機所拍攝到的影像透過電腦進行追蹤，將得到的車輛位置座標可以顯示在使用者介面上。由於本研究是使用環場攝影機進行影像的拍攝，因此所得到的車輛座標位置是在環場影像上，必須透過透視投影的動作將其座標轉換到透視影像上，再將區域座標位置轉換為全域座標位置，作後再對座標位置進行修正，轉換過程如下： Omni-directional Image. Perfective Image. Environmental Image. Environmental Map. 圖 3.1 環境地圖建立流程圖. 3.1 攝影機成像環場攝影機是由一台 CCD（Charge Coupled Device）攝影機及一個半雙曲面的凸型鏡面所構成的（圖 1.6）。根據環場攝影機成像結果分析，空間中實體影像透過環場攝影機的拍攝會產生扭曲變形的現象（Warpping）（圖 3.2（a））。由於環. 23.

(24) 場攝影機拍攝影像扭曲的相當嚴重，無法對本研究在監控平台上的處理直接提供有效的資訊（例如：移動物體的正確位置），這種情況會對最終實驗的結果產生很大的干擾，透過反扭曲（Dewarpping）的技術可以對此一現象做大幅改善，將扭曲的影像恢復成較接近真實場景的透視影像（Perfective View）（圖 3.2（b））。另一方面，由於本研究的室內實驗停車場裝置了八支環場攝影機，透過不同的環場攝影機所拍攝出各自不同的影像，可以方便日後繪製停車場全景圖，在這之前必須先將這些環場攝影機所拍攝出的影像進行拼湊（Mosiacing）的動作（圖 3.3）。為了達到以上這些目的，首先要將影像平面的像素座標從區域座標系統透過雙曲鏡面的反射，重新對應到新的全域座標系統（圖 3.4）。. （a）. （b）圖 3.2 （a）環場影像（b）透視影像. 24.

(25) 圖 3.3 環境影像首先定義雙曲線公式（3.11），a 和 b 分別代表鏡面雙曲線的兩軸（Semiminor and Semimajor），值分別設定為 67.082（cm）及 150.000（cm）。. Z2 X 2 +Y 2 −1 = b2 a2. （3.1）. ，圖 3.6 中，OM 與 OC 分別代表雙曲面的內焦（Inner Focus）及外焦（Outer Focus）座標點分別設為（0,0,c）及（0,0,-c），c 表示原點到內焦（或原點到外焦）的距離，即 OM O 。由（3.6）可得，c = 165.316768（cm）。. c = a 2 + b2. （3.6）. 接著定義區域座標系統中的攝影機成像平面對應到全域座標系統的 Z 軸值為. f − c ，即影像平面上所有的點皆可以（x,y,f-c）表示，其中 x、y 為攝影機成像平面的區域座標值，f 表示該環場攝影機的焦距，值為 1592（cm）。通過 OC 與攝影機成像平面上的任意點（x,y,f-c）畫出一條直線 L1 ，與鏡雙曲面交於一點 Qn ，在影像平面的座標系統通過一點座標（x,y）。 L1 直線方程式定義為 25.

(26) X = t  Y = ty / x  Z = (tf / x) − c . （3.7）. +. ，其中 t 是一個參數， t ∈ R 。. L2. L1. 圖 3.4 環場影像成像原理. 26.

(27) 我們將（3.6）與（3.7）代入（3.5）運算，可以得到另一方程式：. (. (. tf t2 y2 − c)2 t2 + 2 x x −1 = 2 2 b a t2 f 2 tf t2 y2 2 2 − 2 c + c ) t + x2 x x2 − 1 = b2 a2. t2 f 2 tf t 2 y2 2 2 2 2 2 a ( 2 − 2 c + c ) − a b = b (t + 2 ) x x x 2. a 2 f 2t 2 b 2 y 2t 2 a 2 ft 2 2 −b t − −2 c − a 2b 2 + a 2 c 2 = 0 2 2 x x x t2(. a2 f 2 b2 y 2 a2 f 2 − b − ) − t (2 c ) − a 2b 2 + a 2 c 2 = 0 2 2 x x x. f2 y2 f 2 t [ a 2 − b (1 + 2 )] − t (2a 2 c ) + a 4 = 0 x x x 2. 2. 整理可得. 2a 2 c t=. f f f2 y2 ± (2a 2 c ) 2 − 4[a 2 2 − b 2 (1 + 2 )]a 4 x x x x 2 2 f y 2[a 2 2 − b 2 (1 + 2 )] x x. 2 2 f y2 4 2 f 6 f 4 2 2a c ± 4a c 2 − [4a 2 − 4a b (1 + 2 )] x x x x t= 2 2 f y 2[a 2 2 − b 2 (1 + 2 )] x x 2. 2 2 f y2 2 2 f 4 f 2 2 a c ± a a c 2 − [a 2 − a b (1 + 2 )] x x x x t= 2 2 f y a 2 2 − b 2 (1 + 2 ) x x 2. 27.

(28) 2 f y2 2 2 2 f 2 2 a c ± a a (c − a ) 2 + a b (1 + 2 ) x x x t= 2 2 f y a 2 2 − b 2 (1 + 2 ) x x 2. 2 f y2 2 2 f 2 2 a c ± a a b 2 + a b (1 + 2 ) x x x t= 2 2 f y a 2 2 − b 2 (1 + 2 ) x x 2. a 2c t=. f f 2 + y2 ± a a 2b 2 (1 + ) x x2 f2 y2 a 2 2 − b 2 (1 + 2 ) x x. （3.8）. +. ，可得到該直求出兩根解，我們只取正值（ t ∈ R ），以 t1 表示。將 t1 代回（3.7）線與鏡雙曲面的交點 Qn （ t1 , t1. f y , t1 − c ）。通過 Qn 與 OM 得到另一條直線 x x. L2 ，將此直線方程式定義為. X = t  Y = ty / x  Z = tf / x − 2tc t + c 1 . （3.9）. +. 假設影像平面對應到全域座標系統，Z 軸值為 k（ k ∈ R ），L2 交該平面於 Q 點，假設該點座標為（ Qx , Qy ,k）令 tf. t(. / x − 2tc t1 + c = k. f 2c − ) = k −c x t1. t=. t (k − c) x k −c = 1 f 2c t1 f − 2 xc − x t1 28.

(29) t1 (k − c) x  Q = x  t1 f − 2 xc   Q = t1 ( k − c) y  y t1 f − 2 xc. （3.10）. 將（3.8）所求得的 t1 代回（3.10）.  f (b 2 − c 2 ) x Qx = 2 2 (b + c )k − 2bc x 2 + y 2 + k 2   f (b 2 − c 2 ) y Q =  y (b 2 + c 2 ) k − 2bc x 2 + y 2 + k 2 . （3.11）. 透過最後的公式（3.11），便可將原本拍攝到嚴重扭曲變形的影像，修復回與實際場景影像接近的影像。由於修復過的影像周圍會產生許多模糊不清的放射狀情況（圖 3.2（b）），這些資訊對於本研究並沒有太大的幫助，為了提供之後追蹤系統有效的資訊，對於這些影像只擷取中間較為清晰的畫面（圖 3.5）。. 3.2 影像校正上一節中利用（3.11）將扭曲影像（Omni-directional Image）恢復成透視影像（Perspective Image），並且透過使用者所定義的有效範圍將影像可拼湊如圖 3.2 所示，該影像稱為環境影像（Environmental Image），再將此環境影像以人工手繪的方式可以得到圖 3.6 的結果，稱作環境草圖（Environmental Manuscript）。. 29.

(30) 圖 3.5 有效區域. 圖 3.6 環境初稿. 圖 3.7 環境地圖 30.

(31) 環境初稿與環境地圖（Environmental Map）（圖 3.7）比較之下還是會有出現誤差的情況，為了降低這種現象對本研究產生的影響，在這一節提出電腦視覺（Computer Vision, CV）中的幾何攝影機校正（Geometric Camera Calibration）及 Linear Least Square Method 來解決這個問題。進行攝影機的校正時，必須考慮兩種攝影機參數的影響，分別為內部參數（Intrinsic Parameters）以及外部參數（Extrinsic Parameters），在本研究中僅探討外部參數所帶來的影響，針對全景影像 C（圖 1.3）與反扭曲拼湊影像 W （圖 3.5）找出兩者之間的對應關係。首先定義攝影機校正公式（3.12）， C. ，其中. C. P = M WP. P = [x '. （3.12） T. y ' 1] 表示全景影像上的像素座標位置， W P = [ x. 表示在反扭曲拼湊影像上的像素座標位置， M. 的校正矩陣， R. T. y 1]. R O =  ， M 為一個 3 by 3 0 1. t t  =  1 2  為旋轉矩陣（Rotation Matrix）， O = [t5  t3 t 4 . T. t6 ] 為位. 移向量（Translation Vector）， {ti }i =1,...,6 ∈ R 。接著從透視影像中挑選出的數個特徵點座標 ( x, y ) ，以及與這些特徵點對應至全景影像的座標 ( x′, y ′) ，將這些點座標一一代入（3.12），可以得到一個六元一次聯立方程組（3.13）。. 31.

(32)  x1 ' = x1t1 + y1t2 + 0t3 + 0t4 + t5 + 0t6   y1 ' = 0t1 + 0t2 + x1t3 + y1t4 + 0t5 + t6  x2 ' = x2t1 + y2t2 + 0t3 + 0t4 + t5 + 0t6   y2 ' = 0t1 + 0t2 + x2t3 + y2t4 + 0t5 + t6  M   xn ' = xn t1 + yn t2 + 0t3 + 0t4 + t5 + 0t6   yn ' = 0t1 + 0t2 + xnt3 + yn t4 + 0t5 + t6. （3.13）. 將（3.13）以矩陣可表示為，. UT = V  x1 0   x2  0 U = M   xn 0   M. （3.14）. T = [t1 t2 V =  x1'. 0 x1 0 x2 M 0 xn M. y1 0 y2 0 M yn 0 M. t3 y1'. t4 x2 '. 0 y1 0 y2 M 0 yn M. t5. 0 1  0  1 M  0 1  M . 1 0 1 0 M 1 0 M. T. t6 ]. y2 ' ... x3'. y3' .... T. U 和 V 會隨著輸入的點數不同而有所改變，U 是一個 m by 6 的矩陣，T 為所要求的解， V 是 m by 1 的矩陣，其中 m 是輸入的點數。由於（3.12）中包含有 t1 到 t6 六個參數值，因此需要六條以上的方程式才可求出這些參數，即至少需要三對以上的點座標，在本研究中為了得到更佳精確的實驗結果，將所有透視影像中取得的特徵點（共 116 個）皆代入計算。表 3.1 為利用 Linear Least Square Method（3.15）所求得的數據。 32.

(33) T = (U T U ) −1U T V ，其中 U. +. （3.15）. = (U T U ) −1U T 為 U 的 pseudo inverse。表 3.1 校正矩陣.  2.0823 -0.0495 -142.2398   0.1111 2.0022 -162.0952     0  0 1. 1.8135 0.1168 -143.6494   0.1333 1.8666 81.6688     0  0 1. 攝影機一號. 攝影機二號. 1.9884 -0.0966 147.7995   0.1231 2.0189 -158.0205     0  0 1.  1.9862 -0.0142 137.8453  -0.0392 1.9988 67.3252     0  0 1. 攝影機三號. 攝影機四號.  2.0256 -0.0768 398.8315   0.1407 1.9866 -172.3693    0  0 1.  2.0720 -0.0483 388.1822   0.0389 2.1418 1.6456    0  0 1. 攝影機五號. 攝影機六號.  2.1753 -0.1131 671.6380   0.0317 2.0981 -165.1985    0  0 1.  2.2743 -0.0977 634.5074   0.0406 2.1550 24.0031     0  0 1. 攝影機七號. 攝影機八號. 33.

(34) 第四章車輛偵測在這一章會針對本論文是如何偵測出停車場中移動車輛的過程詳加介紹。在 4.1 節會先介紹使用的背景相減法，接著在 4.2 節會提出注意力區域的名詞並且加以定義解釋以及注意力區域的尋找方法，緊接著 4.3 節會介紹如何找出現在時間的背景邊緣影像，最後在 4.4 節使用水平集方法找出前景物，得到結果，流程圖如下：. Image Sequences. ROA Generation. Edge Detection. Object Detection. Vehicles 圖 4.1 車輛偵測流程圖. 4.1 連續影像相減法前景物擷取的技術，現在常見的方法有兩種，分別是背景相減法（Background Subtraction Method）以及連續影像相減法（Temporal Difference Method）。背景相. 34.

(35) 減法在技術上是較為簡單的，但在實作上受到的限制也比較大，因為所建立出的背景影像為固定的，當下面這些狀況發生時，就有可能會出現偵測出非欲偵測物的錯誤資訊： ‧光線出現閃爍（例如：日光燈）或光源有變換（例如：白光變成黃光），此時影像與建立的背景影像差異過大，會造成嚴重的錯誤資訊。 ‧攝影機拍攝畫面中，出現與初始建立的背景影像不一樣的物體時，在影像相減之後會被視為移動的物體，但實際上該物體應為背景的一部分。 ‧在建立好背景影像後，攝影機因各種的因素造成拍攝的角度或是架設位置出現偏移，拍攝得到的影像與起初架設時拍攝出的畫面有差異，此時若使用之前建立的背景影像來進行背景相減的動作，會發生嚴重的錯誤。為了解決背景相減法出現的這些問題，本研究使用的技術是連續影像相減法，目的是為了隨時應付停車場中場景的變化，例如：停車場中原本已經停放的車輛，會被認定做背景，假如使用背景相減法，會出現很大的問題；若使用適應性背景相減法，雖然可以解決上述的問題，但因為過程中需要一段時間做修正，對於本研究需要達到即時追蹤的目的，沒辦法立刻提供有效資訊可能造成追蹤失敗。. 4.2 注意力區域為了提高偵測到物件的正確率，本研究使用連續影像相減法決定注意力區域. 35.

(36) （Region of Attention, ROA）。所謂的注意力區域，就是指影像中車輛可能出現的區域，對於這個區域提高更多的注意力去尋找移動的車輛，如此不但可以減少影像上搜索的範圍，也可以提高尋找到車輛的準確度。首先對連續輸入的影像 It 和 It −1 （圖 4.2）做比較，企圖找出移動中的車輛，並且利用（4.1）得到比較後的影像資訊 Dt ，.  ( x, y ) = 1, Dt =  ( x, y ) = 0,. if I t ( x, y) − I t −1 ( x, y) ≥ π if It ( x, y) − I t −1 ( x, y) < π. （4.1）. 其中 π 為根據當時拍攝環境所訂下的固定門檻值， Dt 顯示出這兩張連續影像的差異部分（圖 4.3）。藉由 Dt 可以判斷出移動的物體在影像中的移動範圍，接著將變動區域 Dt 與前一張影像 It −1 中的偵測到的物件（圖 4.4）作結合（4.2）。從圖 4.3 可以觀察到影像中出現許多不必要的雜點（noise），因此再將圖 4.3 和圖 4.4 結合之後要再做由 Erosion 及 Dilation 形成的 Opening 和 Closing 來去除雜點，最後可以得到一個全新的注意力區域 At （圖 4.5）。. At = {Dt ( x, y ) ∪ Tt −1 ( x, y)}. （4.2）. 圖 4.2 連續輸入影像 36.

(37) 圖 4.3 相減影像. 圖 4.4 物件邊緣影像. 圖 4.5 注意力區域. 4.3 邊緣偵測（Edge Detection）） 4.3.1 背景邊緣建立因為車輛與背景會出現強烈的邊緣對比，因此選擇使用邊緣偵測的方法來偵測車輛，邊緣偵測的方法這邊採用 Canny Edge Detection。由於對車輛進行追蹤前，必須要先找到影像中移動的車輛，為了要把這些移動的車輛抓出來，就必須要先知道影像中哪些邊緣屬於車輛，哪些邊緣屬於背景。首先對前一張影像 It −1 作 Canny 邊緣偵測得到邊緣影像 Et −1 （圖 4.6），利用 37.

(38) 圖 4.6 邊緣影像. 圖 4.7 背景邊緣影像. ，可以利用條件式邊緣影像 Et −1 與上一節 4.2 最後得到的注意力區域 At （圖 4.5）（4.3）建立出現在影像 I t 的背景影像 Bt （圖 4.7），假設 B0 為 t = 0 時候的背景影像。.  Bt −1 ( x, y ) Bt ( x, y) =   Et −1 ( x, y). if ( x, y ) ∈ At otherwise. （4.3）. 4.3.2 物件邊緣建立在 4.2 節中有提到，要得到全新的注意力區域必須要有前一個時間偵測到的物件資訊，為了得到這個物件的邊緣資訊，必須將邊緣影像 Et 中的背景邊緣部份移除。關於影像背景必須注意到背景資訊並非永遠固定不變的，因為環場攝影機在拍攝過程中輸入影像是連續的，前後時間觀察到的背景影像看似相同，但實際上彼此之間仍會產生些微的偏差，為了解決這種背景影像出現偏差的現象，可用（4.4）進行修正，. d ( X , Y ) = V [( X − Y )]. （4.4）. X 、 Y 表示像素鄰近的邊緣方向分佈， V 為此分佈的變化程度。將具有邊緣影像 38.

(39) Et 的像素周圍 5*5（pixel）區塊範圍內的背景邊緣 Bt 與場景邊緣的垂直及水平邊緣比較差異程度。 d 是一個距離函數用來顯示差異的程度，假如距離小於一個門檻值，即判斷此邊緣像素為背景邊緣，將此像素從 Et 中移除；反之則保留此資訊。在執行車輛偵測的過程中，每次都必須要對輸入的整張影像進行運算，如此會浪費許多不必要的時間，為了提高系統實作的效率，可以利用注意力區域 At 輔助，來提高執行的效率。因為移動的車輛資訊只會存在於注意力區域中，也就是不會在非注意力區域中出現，根據這個特性可以將邊緣影像 Et 中不屬於注意力區域的像素皆從 Et 中移除，得到只剩車輛的資訊（圖 4.8）。. Et -1. Et. 圖 4.8 物件邊緣影像 39.

(40) 從圖 4.8 可以發現得到的邊緣可能會出現缺口，造成邊緣不連續的情形，為了解決這個問題，必須對此不連續的邊緣進行邊緣收縮的動作，來得到移動車輛的資訊，這邊提出使用水平集方法（Level Set Method）來達到邊緣收縮的效果。. 4.4 水平集方法水平集方法於 1980 年代由美國數學家 Stanley Osher 和 James Sethian 提出，可應用在影像處理、計算機圖學、計算幾何、最佳化等方面，在這邊主要探討它在影像處理上的應用。所謂的水平集方法是一種追蹤物件介面和形狀的接口技術，這麼方法的優點在於可以輕易的解決物件形狀出現拓墣變化的問題（例如：物件分裂成兩個、出現缺口…等）（圖 4.9）。. 圖 4.9 水平集方法註 2：出自維基百科。. 40.

(41) 在圖 4.9 中的第一列出現了物件分裂的情形，要以參數來描述這樣的形變是十分困難的，從下一列可以觀察到，水平集函數僅僅做了上下的移動就可以達到所要的結果。因此水平集方法也可以解釋成是將低維曲線轉換成高維曲面的方法。首先定義平面上的一條封閉曲線（Closed Curve）Γ ，在三維空間上可視為一個封閉曲面（Closed Surface），即為水平集；而三維座標系統內的圓錐體稱作水平集函數。將曲線當作邊界可分為三個部份，曲線本身、曲線外部及曲線內部。接著在平面上定義一個符號距離函數 ϕ （Signed Distance Function, SDF）（4.5）. 0  ϕ ( x, y ) =  + d  −d . if ( x, y ) ∈ Γ if ( x, y ) ∈ Γin. （4.5）. if ( x, y) ∈Γ out. 其中 Γin 表示封閉曲線內的區域，Γ out 表示封閉曲線外的區域，d 表示像素點 ( x, y) 到封閉曲線的最短距離，距離函數符號取決於該點是位於封閉曲線內部或外部。一般來說，定義曲線內部點的符號為負，反之取正。當 ϕ ( x, y ) = 0 時，表示點 ( x, y) 在曲線上，此時水平集稱作零水平集（Zero Level Set），可以（4.6）表示。. Γ = {( x, y ) | ϕ ( x, y) = 0}. （4.6）. 儘管這種轉換會使得問題在形式上變得更加複雜，但是在問題的求解上卻帶來了很多優點，其中最大的優點在於曲線的拓墣變化能夠得到最自然的處理，並且可以獲得唯一滿足熵（Entropy）條件的解。水平集方法的演算法則就是通過不斷的更新封閉曲線，使得曲線可以達到零. 41.

(42) 水平集的狀態，此時得到的收縮邊界即為物件邊緣。為了要說明曲線的更新（運動）過程，將水平集函數 ϕ 以 Chain Rule 的方式來推導 ϕ 的偏微分方程（Partial Differential Equation）. ∂ϕ (Γ) ∂ϕ (Γ) ∂Γ ∂ϕ = + = ∇ϕΓς + ϕς = 0 ∂ς ∂Γ ∂ς ∂ς. （4.7）. ς 表示更新的間隔時間， Γς 可視為封閉曲線的傳遞速度。定義曲線的傳遞方向為指向封閉曲線內部的方向。令 Γς = F ⋅ −. ∇ϕ ∇ϕ ，其中 − 表示指向曲線內部的法 ∇ϕ ∇ϕ. 向量， F 為速度方程式，代入（4.7）可得到. ∇ϕ ( F ⋅ −. ∇ϕ )ς + ϕς = 0 ∇ϕ. ϕς = F ⋅ ∇ϕ. （4.8）. （4.8）表示零水平集以速度 F 在正常的方向上移動（法向量， ∇ϕ ），此方程稱為 Hamilton-Jacobi Equation，亦為偏微分方程。速度是有正負之分的，如果速度 F 為負，表示曲線的運動方向是朝外部的，如為正，則表示曲線運動方向是朝內部進行的，曲線的運動方向具有單一性，不可能同時往兩個方向運動。速度方程式受到作用力的影響，可定義為三個部份，分別是外部力 Fext 、內部力 Fin 以及圖像力 Fimg （4.9）. F = Fext + Fin + Fimg. （4.9）. （1）外部力（External Force） Fext ：定義用來將封閉曲線 Γ 向內壓縮的常數。外部力通常為一個極小常數，在物件的邊界上可將其設為 0。. 42.

(43) （2）內部力（Internal Force） Fin ：又稱作法向曲率力，將內部力以曲線的曲度表示，（4.10）. Fin = α Γ′′. （4.10）. 其中 Γ′′ 表示 Γ 的二階導數，表示封閉曲線的曲度，（4.11）. Γ′′ = nin − n ⋅ (2n + 1). （4.11）. nin 表示在像素 i 鄰近範圍 (2n + 1) ⋅ (2n + 1) 內屬於 Γin 的像素點個數。α 為權重，當曲線向外凸出時， α > 0 （圖 4.10 中藍色箭頭處）；當曲線向內凹陷時， α < 0 （圖 4.10 中紅色箭頭處）；可根據 n 來調整 α 。（3）圖像力（Image Force） Fimg ：目的是為讓曲線可以被物件的邊界所吸引，定義 Fimg （4.12）. Fimg = β Eobj. （4.12）. β 為權重，當 Fext + Fin > 0 時， β > 0 ；當 Fext + Fin < 0 時， β < 0 。 Eobj 為物件的邊緣強度。為了避免曲線收縮時受到損毀，可固定 β 值為負。. 圖 4.10 邊緣受力示意圖 43.

(44) 實際的執行過程，將影像的像素分為兩個集合 Γin 和 Γ out ，將 Γin 的像素（ ni ∈ Γin ）連接到 Γ out 的像素（ n j ∈Γ ext ）計算速度方程式 F ，若 F > 0 ，則表示. ni 會被移到曲線外部 Γ out ，即 Γin 內的像素數量會變少，使得曲線向內收縮。檢查每次曲線更新的結果，假如曲線的改變已趨於穩定，即表示曲線更新完成。此時，封閉曲線內的像素點集合即為前景物，也就是車輛偵測的結果（圖 4.11）（圖 4.12）。. 圖 4.11 前景物黑白影像. 圖 4.12 前景物彩色影像. 44.

(45) 第五章車輛追蹤 5.1 粒子濾波器簡介粒子濾波器（Particle Filter, PF），亦稱作循序蒙地卡羅法（Sequential Monte Carlo Method, SMC），是利用粒子來模擬中子在空間中運動的方式，進而對物體的運動來進行模擬。由於中子在運動時會相互碰撞，產生不規則隨機的運動軌跡，當粒子在模擬中子進行運動時，主要在模擬中子的運動模式。所謂的蒙地卡羅法，是近代興起的一種科學計算方法，常被用來計算出不容易計算出的數學方程式以及各種自然現象。這個方法是在二次世界大戰末期時，由當時著名的物理學家和數學家：范紐曼（Vonneumann）、鄂拉姆（Ulam）、費米及梅卓普立斯（Metropolis）所發展出來的，因為這方法和賭博的機率、亂數的產生有關，因此才以最著名的賭場蒙地卡羅來命名。原理很簡單，假設我們算一個函數範圍大於. f 的面積，可對 f 找一個類似、更簡單且. f 的函數 g ，可以挑選 g 的面積較容易求得的形狀，例如：方形。然. 後在 g 的範圍內丟入足夠多的隨機亂數（粒子），將落在. f 內的數量（ N f ）除. 以全部亂數的數量（ N ），最後再把這比值乘上 g 的面積（ Ag ），就會得的大約值了。. Af =. Nf N. 如圖5.1所示：. 45. ⋅ Ag. f 面積.

(46) g. 圖5.1 蒙地卡羅法註3：出自維基百科。也就是說，當所灑入的粒子數越多，落在得的比值就會更接近. f 內的粒子數量也會逐漸提高，最後求. f 與 g 之間的比例，即越接近正確的結果。. 粒子濾波器是卡曼濾波器（Kalman Filter）的演進，兩者之間的不同在於基本卡曼濾波器是假設物體的運動是線性（Linear）、呈現高斯分佈的，而粒子濾波器的物體運動是呈現非線性（Nonlinear）、呈現非高斯分佈的，利用不同的函數所得到的期望值作為結果。粒子濾波器主要是利用粒子的密度來表示機率（Probability Density），並利用這些以某種特性分佈（例如：高斯分佈、常態分佈）隨機產生的粒子推出事後機率（Posterior），再將最後求得的期望值來當作預估的結果。由於一般日常生活中物體的運動模式多數為非線性的，為不失一般性，因此使用粒子濾波器與卡曼濾波器比較，所得到的追蹤效果也相對的比較好。使用粒子濾波器之所以會得到比較高的精確結果的原因，是因為粒子濾波器在對物體進行追蹤時，會對追蹤物的狀態一次提出多個隨機預測狀態（Estimate State），並且 46.

(47) 會保留多個可能性較高的預測狀態，而非只保留一個預測狀態。根據這個特性，即便當物體的運動狀態突然出現變化時，造成前一個時間點的預測出現錯誤，亦可以藉由其他保留可能性較高的預測狀態對錯誤進行容忍。. 5.2 傳統粒子濾波器 5.2.1 基本觀念粒子濾波器是根據貝式定理（Bayesian Theorem）探討的一個追蹤問題，將其狀態定義為（5.1）. sk +1 = f ( sk ) + g k. （5.1）. sk 表示在時間 k 的狀態， g k 表示一個非高斯雜訊，狀態轉換函數（State Transition Function）. f 用來將前一個時間得到的狀態轉換成現在時間的狀態。另外再定義. 一個非線性的觀測函數（Measurement Function）（5.2）來表示時間 k 時的觀測值. zk = h( sk ) + ek. （5.2）. zk 表示在時間 k 時得到的觀測值， ek 也是一個非高斯雜訊，函數 h 表示將現在時間得到的狀態透過觀察所得到觀測值 zk 。為了將這兩個函數應用在實作上，將（5.1）及（5.2）以機率的型式來表示，分別如下：. p ( sk | sk −1 ). （5.3）. p ( zk | sk ). （5.4）. 公式轉換過程以（5.1）為例，假設物件狀態包含有位置 ( x, y , z ) 、速度 (vx , v y , vz ). 47.

(48) 六個特徵，以向量 sk 表示. sk =  x. y. z vx. vz . vy. T. 等加速度直線運動公式以及速度公式如下. 1 D( 位移 ) = v0t + at 2 2 v = v0 + at 將（5.1）以上述兩個公式表示. sk +1 = P ⋅ sk + Q ⋅ g k. （5.5）. 、 Q （6 by 3）分別定義兩個矩陣 P （6 by 6）. τ 2  ⋅I Q=2     I ⋅τ .  I I ⋅τ  P=  0 I . 其中 I 為 3 by 3 的單位矩陣， τ 表示單位時間。.  I I ⋅τ  P ⋅ sk =    x I 0   =  x + vx ⋅τ. y. z vx. y + v y ⋅τ. vy. vz . z + vz ⋅τ. T. vx. vz . vy. T. τ 2  ⋅I Q ⋅ gk =  2  gk    I ⋅τ  1 =  axτ 2 2 其中 g k. = [ ax. ay. 1 a yτ 2 2. 1 azτ 2 2. axτ. a yτ.  azτ  . T. az ]T ，表示物件加速度所形成的向量。透過上面的公式可. 以找出物件新的狀態中的位置 ( x′, y ′, z ′) 和速度 (vx′ , v y′ , vz′ ) 。 48.

(49) 1 1 1 x′ = x + vx ⋅τ + axτ 2 ， y′ = y + v y ⋅τ + a yτ 2 ， z′ = z + vz ⋅τ + azτ 2 2 2 2 vx′ = vx + axτ ， v y′ = v y + a yτ ， vz′ = vz + azτ （5.2）推導到（5.4）的過程亦相同。利用 Chapman–Kolmogorov Equation（馬可夫鏈，Chain Rule）可將事前機率. p ( sk | z1:k −1 ) 表示如下， p ( sk | z1:k −1 ) = ∫ p ( sk | sk −1 ) p ( sk −1 | z1:k −1 )dsk −1 在（5.6）式中，. （5.6）. p ( sk | sk −1 ) 為轉移模型（Translation Model），表示從前一. 個時間狀態轉換為現在狀態的機率密度，. p ( sk −1 | z1:k −1 ) 表示在前一個時間所需. 要的機率密度函數（Probability Density Function），z1:k −1 表示在連續時間 1 到 k − 1 累積的所有觀測值。當時間 k. = 1 ， p ( sk | z1:k −1 ) = p ( s1 ) ； k = 0 ，. p ( sk −1 | z1:k −1 ) = p ( s0 ) 。而 p ( sk −1 | z1:k −1 ) 須為已知，（5.6）式才可計算。接著將（5.2）所得到在時間 k 時的觀測值 zk 定為已知條件，將求得的事前機率. p ( sk | z1:k −1 ) 及觀測模型（Observation Model） p ( zk | sk ) 利用貝氏定理. （Bayes Rule）及聯合機率（Joint Probability）兩個觀念，可結合得到事後機率（5.7）. p ( sk | z1:k ) =. 其中. p ( zk | sk ) p ( sk | z1:k −1 ) p ( zk | z1:k −1 ). （5.7）. 1 為正規化常數（Normalizing Constant），定義此常數 p ( zk | z1:k −1 ). p( zk | z1:k −1 ) = ∫ p ( zk | sk ) p ( sk | zk −1 )dsk （5.7）證明推導過程如下： 49. （5.8）.

(50) p ( sk | z1:k ) =. p ( z1:k | sk ) p ( sk ) p ( z1:k ). =. p ( zk , z1:k −1 | sk ) p ( sk ) p ( zk , z1:k −1 ). =. p ( zk | z1:k −1 , sk ) p ( z1:k −1 | sk ) p ( sk ) p ( zk | z1:k −1 ) p ( z1:k −1 ). =. p ( zk | z1:k −1 , sk ) p ( sk | z1:k −1 ) p ( z1:k −1 ) p ( sk ) p ( zk | z1:k −1 ) p ( z1:k −1 ) p ( sk ). =. p ( zk | sk ) p ( sk | z1:k −1 ) p ( zk | z1:k −1 ). 從（5.6）及（5.7）交互觀察可以發現，（5.7）表示根據前一個時間點的資訊得到現在的事前機率，同時更新現在的事後機率的狀態，（5.6）再利用更新過後的事後狀態來更新下一個時間的事前機率，如此重複交互運算，以達到狀態不斷更新的目的，以圖 5.2 表示. observed. zk − 2. zk. zk −1. p( zk | sk ). unobserved. sk − 2. sk −1. p(sk | sk −1 ). sk. 圖 5.2 物件狀態更新示意圖在本研究中的事後機率是由粒子分佈的密度所推導出，可以將連續時間 1 到. k 得到的所有觀測值以及目前狀態 sk 所形成的事後機率密度函數以離散的方式 50.

(51) 表示如下 l. p ( sk | z1:k ) = ∑ wki δ ( sk − ski ). （5.9）. i =1. p ( sk | z1:k ) 為一個連續機率分佈（Continuous Probability Distribution）， l 表示粒子個數，當選取的粒子數量越多（ l. → ∞ ）時，所得到的結果也會越接近理想。. wki 表示每個粒子 i 的權重。 δ 是一個 Dirac delta Function，表示在 sk 空間上如果 i. 出現與粒子 sk 吻合（hit）的現象時， δ ( sk. − ski ) = 1 。根據此離散的式子，可以. 表示出事後機率的分佈，亦可表示為權重分佈函數。 i. i. 每個粒子的權重 wk 並非一個定值，它可以根據前一個時間的權重 wk −1 、轉移模型. p ( ski | ski −1 ) 和觀測模型 p ( zk | ski ) 不斷的更新粒子的權重， i k. i k −1. w ∝w. p ( zk | ski ) p ( ski | ski −1 ) q ( ski | ski −1 , z1:k ). （5.10）. 函數 q 表示在連續時間 1 到 k 得到的觀測值及前一個時間的狀態下，產生現在狀態的重要性密度函數（Important Density Function），根據此重要性密度函數的分佈特性隨機選取粒子，並計算各個粒子的權重重要性。在部分的研究中會假設事前機率即為重要性密度函數，. q ( ski | ski −1 , z1:k ) = p ( ski | ski −1 ) 因此可將（5.10）簡化為. wki ∝ wki −1 p ( zk | ski ) 其中觀測模型. （5.11）. p ( zk | ski ) 可以用粒子與追蹤物之間的相似程度（Likelihood）來定 51.

(52) 義，在下一小節會針對本研究中所提出的三種相似程度定義進行說明。透過上面這些式子，可對每一個時間的狀態計算其機率密度函數，並利用期望值預估現在可能的狀態。演算法如下：步驟一：定義空間上粒子的數量 l ，並給予事前機率初始值機率分布），由. p ( s0 ) （起始的. p ( s0 ) 定義出每個粒子的初始狀態 s0i 以及初始權重 w0i 。. 步驟二：粒子根據轉移模型. p ( ski | ski −1 ) 開始移動，找出現在每個粒子的狀態. ski 。步驟三：重新計算每個粒子的權重（5.12），在本研究中不考慮現在時間粒子權重與前一個時間粒子權重之間的影響，可將（5.11）再簡化為. wki = p ( zk | ski ). （5.12）. 步驟四：將權重進行正規化（5.13）. i k. w =. wki （5.13）. l. ∑ wki i =1. 步驟五：求出事後機率. p ( sk | z1:k ) ，找出期望值狀態作為輸出結果（5.14）. l. p ( sk | z1:k ) = ∑ wki δ ( sk − ski ) i =1 l. E ( f ( sk )) = ∑ wki f ( ski ). （5.14）. i =1. 步驟六：重新取樣（Resample），利用一個門檻值（Threshold）剔除權重值. 52.

(53) 較小的粒子，保留剩餘權重值較大的粒子，再根據剩餘粒子權重大小的比例，重 i. 新分配成 l 個粒子，使得 sk 可以符合事後機率步驟七： k. p ( sk | z1:k ) 分佈。. = k + 1 ，回到步驟二。. 以示意圖表示，如圖 5.3，. Prediction Move particle Resample. Weight Sample. 圖 5.3 粒子濾波器流程圖圖 5.3 中可利用（5.10）求得其權重分佈函數，根據轉移模型子移動的動作。. 53. p ( ski | ski −1 ) 進行粒.

(54) 5.2.2 色彩直方圖在粒子濾波器中觀測模型. p ( zk | sk ) ，在本研究中定義為粒子所形成的樣本. 與追蹤物之間的相似程度（Likelihood）。為了判斷由一個粒子代表的矩形區域與追蹤物是否相似，在這邊使用三種特徵作判斷，分別是（1）粒子矩形面積大小與車輛面積大小之間的相似程度（面積相似程度）、（2）粒子矩形中心位置與車輛預測位置之間距離的相似程度（距離相似程度）以及（3）粒子矩形形成的色彩直方圖（Color Histogram）與車輛形成的色彩直方圖的相似程度（色彩相似程度）三個部份。本節主要探討該如何進行色彩直方圖的相似程度比較。假設在現在的時間，影格上有抓到一個物件（Target），首先定義特徵值 z 表示該物件的顏色（Color）及模型（Model）特徵，物件中心座標位置為 y ，偵測到的物件其特徵分佈函數以. pz (y) 表示，樣本物件的特徵分佈函數為 q z ，要得到相似程度比較的結果，也. 就相當於要判斷出這兩個特徵分佈函數的差異度，這邊採用 Bhattacharyya coefficient 的方法，表示如下：. ρ (y) ≡ ρ [ p(y), q ] = ∫ pz (y)qz dz. （5.15）. 當兩個物件模型越類似的時候，所得到 Bhattacharyya coefficient 的值就會越大；反之，兩個物件模型越不相似的時候，所得到 Bhattacharyya coefficient 的值就會越小。以色彩直方圖表示該式，則令 z. = 1,..., v，v 為色彩直方圖的切割數（v-bin 色. 54.

(55) ˆ，彩直方圖），接著定義物件在座標位置 y 時的預估函數 Pˆ (y) 和樣本預估函數 Q v. ρ (y) ≡ ρ[ Pˆ (y), Qˆ ] = ∑ pˆ z (y) qˆ z. （5.16）. z =1. v. Pˆ (y) = { pˆ z (y)}z =1,...,v （ ∑ pˆ z (y) = 1 ） z =1 v. Qˆ = {qˆ z }z =1,...,v （ ∑ qˆ z = 1 ）。 z =1. 假設 {x ′j } j =1,...,n′ 表示樣本物件中每一個像素的位置， n′. qˆ z = C ′∑ k ( x ′j )δ [ h(x ′j ) − z ] 2. （5.17）. j =1. 其中函數 δ 為 Kronecker Delta Function， C ′ 為正規化常數（Normalizing v. Constant），目的是為了使. ∑ qˆ. z. = 1 ，函數 k 為一個核函數（Kernel Function）。. z =1. 核函數 k 包含有三種形式，分別是 Epanechnikov Kernel、Uniform Kernel、Normal Kernel（圖 5.4），定義如下：. 圖 5.4 核函數註 4：出自維基百科。. Epanechnikov Kernel. 2 n ≤1 cd (1 − n ) k E (n) =  otherwise  0. 55.

(56) Uniform Kernel.  c n ≤1 kU (n) =  d otherwise 0. Normal Kernel. k N (n) = cd ⋅ exp( −. 1 2 n ) 2. 其中 cd 為一常數。假設在現在時間得到的物件中心位置為 y ，物件每個像素在影像上位置為. {x i }i =1,...,n ，接著再定義一個函數 h : R 2 → {1...v} ，表示透過函數 h 可將物件在影像上位置所代表的顏色值轉成以色彩直方圖表示，透過函數 δ 及 h 可以得到（5.18） 2. n. y - xi pˆ z (y) = C ∑ k ( )δ [h(x i ) − z ] r i =1. （5.18）. v. C 亦為正規化常數，使得 ∑ pˆ z (y) = 1 。根據 Bhattacharyya coefficient，的方法 z =1. ˆ 為在前一個時間的時候對現在時間所做的預將這兩個函數作差異度的比較，令 y 估位置結果，將（5.16）以泰勒展開式（Taylor Series）可以得到. ρ [ Pˆ (y), Qˆ ] ≈. 1 v 1 v qˆ z ˆ ˆ ˆ p (y) q + pˆ z (y) ∑ ∑ z z ˆ 2 z =1 2 z =1 pˆ z (y). （5.19）. 將（5.18）代入（5.19） 2. y − xi 1 v C n ˆ ˆ ˆ qˆ z + ∑ wi k ( ρ [ P(y), Q ] ≈ ∑ pˆ z (y) ) 2 z =1 2 i =1 r v. 其中 wi = ∑ δ [ h(x i ) − u ] z =1. （5.20）. qˆ z w i ˆ ，定義 i 為現在時間像素的權重。 pˆ z (y). 56.

(57) 5.2.3 重新取樣（）重新取樣（Resample）在傳統粒子濾波算法中存在的粒子退化（Particle Degeneracy）現象和重新取樣過程中所產生的粒子匱乏（Particle Impoverishment）問題。所謂的粒子退化，是指在重新取樣過後，會有少數的粒子被複製下來，而大多數的粒子沒有後代，使得粒子喪失多樣性的現象。隨著反覆（迭代）運算次數的增加，在所有粒子中僅有少部分的粒子會保有其權重值，其餘大多數的粒子的權重值皆會趨近於零，亦即花費了大量不必要的運算時間在計算那些對於估計結果幾乎毫無意義的粒子權重更新上。由於重新取樣的過程中，會偏好那些適合的粒子（權重值較大的粒子），因此那些權重值逼近於零的粒子就會被捨棄掉，從而引發出粒子退化的現象。粒子匱乏的問題，即為採樣枯竭的問題，指的是在經過多次的反覆重新採樣之後，所有的粒子都會逐漸收斂到一個樣本的附近。這對於整個系統在動態時的預測結果會造成較大的偏差，甚至出現不收斂的現象。. 5.3 聯合粒子濾波器聯合粒子濾波器（Joint Particle Filter, JPF）是粒子濾波器的另一種變形，主要的目的是為了解決傳統粒子濾波器在多輛車追蹤時所出現的不穩定現象。傳統粒子濾波器是針對在前一張影格中所偵測到的各個追蹤物（車輛）一一進行計算，進而可以得到各自獨立的預測結果，也就是說傳統粒子濾波器對於物. 57.

(58) 件是以採用各自獨立（Independent）的處理方式。在計算上也顯得比較簡單，只需要考慮該物件過去所累積下來的資訊，但相對的在結果預測上也容易出現錯誤，例如：所預測到的位置可能會與其他物件預測的位置重疊。聯合粒子濾波器一般用於多物件追蹤（Multi Target Tracking）上，使用時必須對前一個影格中所有偵測到或已預測出的追蹤物同時進行計算，並且考慮各個物件在計算過後位置是否會出現重疊的問題。與前者相比較之下，進行運算所耗費的時間會較多，但是追蹤的結果會比前者來得穩定許多。聯合粒子濾波器必須去預測所有欲追蹤物的位置，並且同時考慮到所有物件的狀態關聯。首先對每個物件定義其狀態函數. ski +1 = f ( ski ) + g ki 其中 i. （5.21）. = 1,...,τ ， τ 為物件的個數。另外定義觀測函數 zk = h( ski ) + ek. 將物件在時間 k 的狀態以集合 S k. （5.22）. = {s1k , sk2 ,..., sτk } 表示，以集合 z k = {zkj }Mj =k1 表 k. 示在時間 k 的觀測結果，其中 M k 為觀測值個數。再以一個集合 Z 表示經過了 k 時間所累積記錄下來的觀測結果， Z. k. = {z 0 , z1 ,..., z k } 。最後定義所有物件的. 狀態與所有觀測結果之間的相連關係事件（Joint Association Event） φ ，每個事件都是由數對 ( j , i ) 所構成，其中 ( j , i ) ∈ {0,..., M k } × {1,...,τ } ，例如：圖 5.5 可以用 (1,3) 、 (3,1) 、 (2, 2) 、 (4, 4) 表示其中的一種相連關係事件。將所有的相連關係事件 φ 以集合 Φ j ,i 表示。在此對應關係上 M k 與 τ 的數量不一定要相. 58.

(59) 同，例如：圖 5.6 表示第三個物件並沒有得到相對應的觀測值，可能為發生偵測錯誤的現象。在傳統粒子濾波器中考慮粒子權重的部份，主要分為（1）面積相似程度、（2）位置相似程度以及（3）色彩相似程度三個特徵來考慮；在聯合粒子濾波器除了. 1 1 3 3. 2 4. 4. 2. 物件狀態. 觀測結果圖 5.5 相聯關係事件. 1 1 3 3. 2 4. 2. 物件狀態. 觀測結果. 圖 5.6 相聯關係事件錯誤情形考慮上述的三種特徵之外，還必須多考慮到（4）在現在的時間所得到的觀測結果中，各種相連關係事件 φ 可能出現的機率來當作權重之一（5.23）. l j .i =. ∑ φ. P (φ | Z k ). （5.23）. ∈Φ j ,i. 59.

(60) 考慮單一發生相連關係事件 φ 的機率. P (φ | Z k ) = P (φ | z k , Z k −1 ). （5.24-1）. = P (φ | z k , S k ) =. （5.24-2）. 1 P( z k | φ , S k ) P (φ | S k ) c. （5.24-3）. 。其中 c 為一正規化常數（Normalizing Constant）（5.24-1）可利用馬可夫性質（Markov Property）可以得到（5.24-2）。關於馬可夫性質的說明如下：假設有兩個骰子 U 和 V ，骰子 U 有五面是 V ，一面是. U ；骰子 V 則有五面是 U ，一面是 V ，遊戲規則是一次只能丟一顆骰子，丟出那面所顯示的字母為下次要擲得骰子。假設每次丟骰子的動作都是獨立的，若第一次丟擲骰子 U ，則有 1. 6 的機率為 U ，有 5 6 的機率為 V ，因此在丟擲第二. 次骰子時，得到 U 和 V 的機率為. 1 5 P (T2 = U ) = P (T1 = U ) ⋅ + P (T1 = V ) ⋅ 6 6 5 1 P (T2 = V ) = P (T1 = U ) ⋅ + P (T1 = V ) ⋅ 6 6 其中 T1 和 T2 分別代表第一次和第二次丟擲的結果。在丟擲 n 次之後，會發現到第 n 次得到的結果只會與第 n − 1 次得到的結果有關，以條件機率表示為. P (Tn = U | Tn −1 = U ) =. 1 6. P (Tn = U | Tn −1 = V ) =. 5 6. 60.

(61) P (Tn = V | Tn −1 = V ) =. 1 6. P (Tn = V | Tn −1 = U ) =. 5 6. 一旦 Tn −1 為已知，就可以決定出 Tn 的機率分佈。要注意的是 Tn 的機率分佈與 Tn −1 有關，而 Tn −1 的機率分佈又跟 Tn − 2 有關，因此 Tn 與 Tn − 2 並不算是互相獨立的事件，但他們之間又沒有直接的相依性，而有以下的性質. P (Tn | Tn −1 ) = P (Tn | Tn −1 , Tn − 2 ) = P (Tn | Tn −1 , Tn − 2 , Tn −3 ,...) 這個特性稱之為馬可夫性質。. Tn. Tn − 2. Tn −1 圖5.7 馬可夫性質. （5.24-2）可以再利用貝氏定理與聯合機率（Joint Probability）可得到（5.24-3），推導過程如下：. P (φ ) P ( S k | φ ) P ( z k | φ , S k ) P (φ | z , S k ) = P( Sk ) P( z k | Sk ) k. P (φ ) =. P (φ | S k ) P ( S k ) P( z k | φ , S k ) P (φ ) P( Sk ) P( z k | Sk ). P (φ | S k ) P ( z k | φ , S k ) = P( z k | Sk ) 分別探討（5.24）中 P ( z. k. | φ , S k ) 以及 P (φ | S k ) 兩個部份。在時間 k 的物 61.

(62) 件狀態 S k 之下可以得到的相連關係事件 φ 機率. P (φ | S k ) = PDτ − z0 (1 − PD ) z0. （5.25）. 其中 PD 表示物件的偵測率（Probability of Detection）， z0 表示物件沒有對應到觀測值的物件個數。計算 P (φ | S k ) 及計算物件1可能是物件a的機率，物件2可能是物件b的機率…依此類推。圖5.8（a）（b）中， z0 1. 2. = 0 ；（c）（d）， z0 = 2 。 a. b. c 3. d. 4. e （a）. 1. （b）. 2. a. 5. 3. b. c. 4 （c）. （d）圖 5.8 物件偵測率示意圖. 圖5.8（a）（c）表示在前一個時間的對現在影像預測的結果，（b）（d）表示現在時間影像中偵測到的物件。 k. 在出現物件狀態 S k 並且發生事件 φ 的條件之下，得到觀測結果 z 的機率表. 62.

(63) 示如下， M k − (τ − z0 ) P ( z k | φ , S k ) = PFA. ∏φ P( z. j k. | ski ). （5.26）. ( j ,i )∈. 其中 PFA 表示發出錯誤訊息（Fault Alarm）的機率，也就是說理論上觀測值數量（ M k ）與實際上物件有對應到觀測值的數量（ τ. − z0 ）會出現落差。以圖5.9為. 例，假設現在有偵測到四個物件，在理論上應該會出現四個對應的觀測結果，但實際上卻只有出現 (1,3) 、 (3,1) 、 (2, 2) 三種對應的關係，觀測值4並沒有被對應到，會發出錯誤訊息，. 1 1 3 3. 2 4. 4. 2. 物件狀態. 觀測結果. 圖 5.9 相聯關係事件發出錯誤訊息. P ( zkj | ski ) 表示物件 i 有對應到觀測值 j 的機率， j k. i k. P( z | s ) =. 1 (2π ) M 2 Σ j ,i. 12. exp(. −ι 2j ,i 2. ). （5.27）. 其中 Σ j ,i 為innovation covariance matrix，ι j ,i 為normalized innovation，可以（5.28）表示. ι 2j ,i = ε Tj ,i Σ −j 1,iε j ,i. （5.28） 63.

(64) j. 其中 ε j ,i 表示觀測結果 zk 與預測結果 zˆk ,i 之間的誤差，. zˆk ,i = h( sˆki ) + ek. ε j ,i = zkj − zˆk ,i = zkj − h( sˆki ). （5.29）. 將上述推導後的結果合併至（5.24-3）可以得到. 1 P ( z k | φ , S k ) P (φ | S k ) c. P (φ | Z k ) =. =. 1 P ( z k | φ , S k ) PDτ − z0 (1 − PD ) z0 c. =. 1 M k −(τ − z0 ) PFA P ( zkj | ski )PDτ − z0 (1 − PD ) z0 ∏ c ( j ,i )∈φ. =. 1 τ − z0 M k − (τ − z0 ) PD (1 − PD ) z0 PFA P ( zkj | ski ) ∏ c ( j ,i )∈φ （5.30） j. 在本研究中針對 P ( zk. | ski ) 的部份，改以計算相似程度的比較，計算的條件與一. 般粒子濾波器相同，使用面積相似程度、位置相似程度以及色彩相似程度三種特徵。整個公式的推導流程可以以流程圖5.10表示：. 64.

(65) l j .i =. ∑ φ. P (φ | Z k ). ∈Φ j ,i. P (φ | Z k ) =. M k − (τ − z0 ) P ( z k | φ , S k ) = PFA. 1 P( z k | φ , S k ) P(φ | S k ) c. ∏φ P( z. j k. | ski ). P (φ | S k ) = PDτ − z0 (1 − PD ) z0. ( j ,i )∈. j k. i k. 1. P( z | s ) = (2π ). M 2. Σ j ,i. 12. exp(. −ι 2j ,i 2. ). 圖5.10 公式推導流程圖在聯合粒子濾波器中，必須要將在過去觀測過程中，物件狀態有關聯的建立其相關性，即為連續物件狀態向量（Joint State Vector）。假設現有物件個數為 τ ，. n （ n = 1,..., N ）表示第 n 個樣本，將物件狀態 S k 定義成整張影像上所有偵測到的物件所形成的連續物件狀態向量， S k (n) N n =1. 連續物件狀態向量 {S k. }. = {s1k , sk2 ,..., skτ } ，定義樣本所形成的. = {sk( n ),1 , sk( n ),2 ,..., sk( n ),τ }nN=1 ，其中 N 表示樣本. （Sample）個數。 m. 圖5.11為一個多維狀態空間（State Space） D 的示意圖，維度的大小可依當時給予狀態特徵的多寡來決定，例如：位置（Location）、面積（Area）、速度 65.

(66) （Velocity）、加速度（Acceleration）、顏色（Color）、邊緣資訊（Edge）、角度（Degree）…等。影像中偵測到的各個物件，可以找出該物件在樣本中對應的物件狀態，依照其求得的物件狀態一一標示在此狀態空間上，空間上每一個點都代 ( n ),i. 表一個物件的狀態 sk. ，點的個數即為物件數量，共有 τ 個（ i. = 1,...,τ ），圖. 5.11中的每個藍點表示在第 n 種樣本的物件狀態。為了增加各個物件之間的關聯性，此特性即為聯合粒子濾波器比傳統粒子濾 (n). 波器功能來的強大的原因，將這些粒子樣本向量收集起來，以一個集合 S k 來表 n. 示，並將集合中的各點分別標示於另一個多維狀態空間 D （Joint State Space）上，即為圖5.12中的紅點。圖5.12中黃點表示在時間 k 時候的物件狀態向量，即為 Sk. = {s1k , sk2 ,..., skτ } ，紅點表示在此空間上對該物件狀態向量進行取樣所得到 (n) N n =1. 的樣本向量， {S k. }. = {sk( n ),1 , sk( n ),2 ,..., sk( n ),τ }nN=1 。. 圖5.11 狀態空間. 66.

(67) 圖5.12 聯合狀態空間 n. 將蒙地卡羅法的觀念應用在聯合粒子濾波器上，在多維空間 D 上任意或依 (n). 某種分佈特性（高斯分佈）散落 N 個樣本（即為紅點），每個樣本以 S k 表示， (n). 其中 n = 1,..., N ， S k 為一狀態向量，分別代表在下個時間點各種不同的預測結果。根據上述觀念可將（5.27）及（5.30）改寫如下，便可得到在不同樣本的狀況下，每個樣本存在的可能性（權重）。 j k. ( n ),i k. P( z | s. )=. P (φ | Z k )( n ) =. 1 (2π ) M 2 Σ j ,i. 12. exp(. −ι (j n,i),2 2. ). （5.31）. 1 τ − z0 M k − (τ − z0 ) PD (1 − PD ) z0 PFA P ( zkj | sk( n ),i ) ∏ c ( j ,i )∈φ （5.32）. l (jn,i) =. ∑. P (φ | Z k ). (n). （5.33）. φ ∈Φ (jn,i). 演算法如下： 67.