國立高雄大學資訊工程學系研究所
碩士論文
以單眼視覺空拍影像實現本體移動估測
Establishing Ego-motion Estimation from Monocular
UAV Image Sequences
研究生：黃鼎翔 撰
指導教授：陳佳妍 博士
致謝
時光匆匆歲月如梭，轉眼已來到寫這份致謝的時候了，猶記得兩年前剛踏進 VGLab 的那天，是一個不太寒冷的一月，家泓學長熱心的向我介紹實驗室的一 切，就此開啟我豐富且充滿回憶的研究所生活，其中的精采程度可能不是一篇致 謝能囊括的，但要感謝的人實在太多了，請容許用這極短的篇幅娓娓道來，他們 都是我生命中不可或缺的助力，能完成碩士學位，實在萬分感激。 首先感謝佳姸老師，您不僅是我的指導教授更是我人生的導師，不只給了我 論文的方向，更給了我未來了路。兩年的時間，帶著我們參與了國內外無數的研 討會，以學術為基石寓教於樂，其中還認識了紐西蘭這個美麗且純淨的國度，開 拓了我的視野。感謝祥任學長在程式上不厭其煩地給予幫助，在遇到瓶頸時伸出 援手，是我學業道路上的一盞明燈。感謝家泓學長在電腦視覺基礎概念上的開導， 讓我們能從最根本的觀念循序了解，你們真是 VGLab 珍貴的資產。感謝實驗室 的大家:郁媚、振錡、芝君、星宇，伴我左右，共享革命溫情，造就一幕幕生命 繽紛的分鏡。感謝資工系最美麗系助淑真，無時無刻的給予幫助，並且讓系辦的 笑聲不間斷。 感謝我的碩士學位口試委員，陳銘志教授與殷堂凱主任，百忙之中參加口試 並且對於我的碩士論文進行指導與修正，這些寶貴的意見使我的論文更加豐富與 完善。 最後，將此篇論文獻給我的家人，在我漫長的求學生涯中，一路上給我扶持 與鼓勵，讓我可以專心完成學業，家人永遠是我最堅強的後盾。 我是個很幸運的人，一路上遇到許多貴人相助，所有的幫助我都會銘記在心， 期許未來的每一天，我也能為需要我的人做點什麼，讓世界變得更好。 未來，我會努力朝著這個方向前進。 黃鼎翔 謹致於 國立高雄大學 資訊工程學系 中華民國一零五年十二月以單眼視覺空拍影像實現本體移動估測
指導教授：陳佳妍 博士 國立高雄大學資訊工程學系 學生：黃鼎翔 國立高雄大學資訊工程學系研究所 摘要 過去十年內三維技術發展迅速，但是大多數的系統裝置在獲取環境資訊的同 時，卻因缺乏移動感測技術，導致資訊的取得不完整。這將使系統無法在較大範 圍的區域中擷取資訊，也無法與在其他位置的系統做連結。本論文根據單眼視覺 的本體移動估測，使系統裝置在移動的同時得知本身的路徑移動狀態，並在多系 統的情況下將不同時間與不同地點的資訊做完整地結合。另外，隨著時代的演變 與科技的進步，三維重建技術越來越盛行，其中生成三維模型的方式包含雷射測 距、結構光源、立體視覺等，本研究則使用單眼視覺之方法實現本體移動估測並 能運用於產生三維點雲以獲得物體模型。 本論文研究使用架設於無人機上的一組鏡頭進行視覺移動估測。第一步先使 用固定焦距、光圈及曝光等參數之相機擷取影像，接著對影像進行校正與對正後 並使用加速穩健特徵演算法(Speed Up Robust Features,SURF)對影像偵測特徵，再 進行空間上的特徵匹配。特徵匹配步驟追蹤各時間特徵點的位置變化，同時對影 像特徵點進行三角測量獲得旋轉矩陣R 與位移向量 t，並結合相機參數計算出特 徵點於三維空間之座標，算時間點特徵之旋轉、位移變化，最後推演出系統的路 徑。 本研究利用無人機所獲得的空拍影像資訊產生場景之本體移動估測路徑，並 能運用於產生具有紋理的細緻三維模型。本研究之成果可廣泛地應用於大規模環 場的資料擷取、防災救災應用、古蹟三維模型重建、生態與地質地理研究、都市 計畫與土地應用開發等。 關鍵字：三維重建、視覺測程、單眼視覺、無人機、影像序列整合特徵Establishing Ego-motion Estimation from Monocular
UAV Image Sequences
Advisor: Dr. Chia-Yen Chen
Department of Computer Science and Information Engineering National University of Kaohsiung
Student:Ting-Hsiang Huang
Department of Computer Science and Information Engineering National University of Kaohsiung
ABSTRACT
In the past decades, 3D related technologies have undergone significant development. Most of the 3D acquisition systems, however, lack self-localization abilities and may face difficulties in the reconstruction of a large scale 3D environment. In this work, we investigate the application of monocular vision in ego-motion estimation. In particular, with respect to images acquired by a single camera mounted upon an Unmanned Aerial Vehicle. In the proposed system, we first estimate the camera parameters from the acquired images. The images are then calibrated and registered. Features within images are detected by Speed Up Robust Features, and spatial feature matching is performed. The Semi-Global Matching method is selected for spatial feature matching to maintain the quality and number of matched features. The features are then combined with the camera parameters to calculate the 3D coordinates of the feature points, as well as the temporal features' rotation and translation variations. Finally, the ego-motion of the system is acquired by optimizing the obtained movements. Experiments are performed to demonstrate the applicability of the proposed method, in particular, towards the 3D reconstruction from aerial images.
Keywords: 3D reconstruction, Ego-Motion Estimation, Visual Odometry, Monocular
目錄
致謝... i 摘要... iii ABSTRACT ... iv 目錄... v 圖目錄... vii 表目錄... ix 符號定義... x 第一章 緒論... 1 1.1 前言... 1 1.2 研究動機... 1 1.3 研究目的... 2 1.4 研究方法... 2 1.5 PHANTOM 3 PROFESSIONAL 四軸飛行器 ... 3 1.6 四旋翼飛行器搭載鏡頭 ... 5 1.6.1 雲台系統 ... 5 1.6.2 相機系統 ... 6 1.7 論文架構... 8 第二章 文獻探討... 9 2.1 使用單眼視覺進行本體移動估測 ... 9 2.2 無人飛行器之發展 ... 12 第三章 視覺測程 ... 15 3.1 前言... 15 3.1.1 相機模型(Camera Model) ... 16 3.1.2 相機校正... 21 3.1.3 VO 實作流程... 26 3.2 特徵點偵測 ... 27 3.2.1 積分影像(Integral Image) ... 28 3.2.2 快速海森特徵檢測 ... 30 3.2.3 SURF 敘述子 ... 33 3.3 離群資料剔除與修正 ... 35 3.4 特徵匹配... 37 3.5 移動估測... 38 3.6 最佳化... 41 第四章 實驗結果與討論 ... 44 4.1 實驗資料蒐集 ... 454.2 單眼視覺測程實驗結果 ... 49 4.2.1 室內環境... 49 4.2.2 室外街景環境-直線 ... 51 4.2.3 室外街景環境-環繞 ... 54 4.2.4 室外環境街景-直線-鏡頭朝下 ... 58 4.2.5 校園街道 曲線 ... 61 4.3 實驗分析... 64 第五章 結論與未來研究方向 ... 65 5.1 結論... 65 5.2 未來研究方向 ... 65 參考文獻... 67
圖目錄
圖1-1 PHANTOM 3 PROFESSIONAL 四軸飛行器 ... 3 圖1-2 四軸飛行器之遙控器 ... 4 圖1-3 雲台系統 ... 6 圖1-4 鏡頭 ... 7 圖2-1 本體移動概念圖 ... 10 圖2-2 本體移動估測之概略圖 ... 11 圖3-1 內、外部參數轉換關係圖 ... 16 圖3-2 針孔相機原理圖 ... 17 圖3-3 相機針孔模型 ... 18 圖3-4 徑向形變示意圖[22]... 19 圖3-5 包含徑向形變之相機投影模型 ... 20 圖3-6 影像擷取前後工作流程 ... 21 圖3-7 Zhang 式方法使用之校正板 ... 22 圖3-8 校正平面與其座標系 ... 22 圖3-9 單應平面示意圖 ... 24 圖3-10 六十四個不同角度遠近所拍攝之校正板影像 ... 25 圖3-11 校正外部參數所得之二十個校正版相對於相機之位置 ... 26 圖3-12 視覺測程(Visual Odometry)架構圖 ... 27 圖3-13 積分影像示意圖 ... 29 圖3-14 利用積分影像計算矩形區域像素總和 ... 29 圖 3-15 盒濾波器替代高斯二階微分運算 ... 31 圖3-16 尺度空間構造示意圖[28]... 32 圖3-17 特徵點定位示意圖 ... 33 圖3-18 哈爾(Haar)小波濾波器 ... 33 圖3-19 主方向示意圖 ... 34 圖3-20 敘述子產生示意圖 ... 35 圖3-21 RANSAC 架構圖 ... 36 圖3-22 RANSAC 演算法概念圖 ... 36 圖3-23 N 點透視問題 PnP 示意圖[19] ... 39 圖3-24 重投影誤差示意圖 ... 41 圖 4-1 本論文實驗無人機 ... 44 圖4-2 實驗環境 ... 47 圖4-3 戶外實驗環境 衛星地圖[39]... 48 圖4-4 室內環境特徵點 ... 49 圖4-5 室內環境 293 張 frame 特徵點偵測 ... 50圖4-6 室內視覺測程結果之軌跡 ... 50 圖4-7 室外街道環境建築物特徵點 ... 51 圖4-8 1442 個 frame 特徵點偵測 ... 52 圖4-9 室外直線視覺測程結果之軌跡 ... 52 圖4-10 室外直線 GPS 軌跡 ... 53 圖4-11 室外直線估測路徑與 GPS 軌跡疊圖 ... 53 圖4-12 1594 個 frame 特徵點偵測 ... 54 圖4-13 室外街道環境建築物特徵點 ... 55 圖4-14 室外街道環境建築物特徵點 ... 55 圖4-15 室外街道環境移動估測結果 ... 56 圖4-16 室外環狀 GPS 軌跡 ... 56 圖4-17 室外環狀估測路徑與 GPS 軌跡疊圖 ... 57 圖4-18 666 個 frame 特徵點偵測 ... 58 圖4-19 戶外街道建築物特徵點 ... 59 圖4-20 室外街道環境移動估測結果 ... 59 圖4-21 室外直線 GPS 軌跡 ... 60 圖4-22 估測路徑與 GPS 軌跡疊圖 ... 60 圖4-23 1000 個 frame 特徵點偵測 ... 61 圖4-24 校園街道特徵點 ... 62 圖4-25 校園街道移動估測結果 ... 62 圖4-26 校園街道 GPS 軌跡 ... 63 圖4-2 校園街道估測路徑與 GPS 軌跡疊圖 ... 63 圖5-1 Ricoh Theta S 全景相機 ... 66
表目錄
表1-1 PHANTOM 3 PROFESSIONAL 詳細規格 ... 5 表1-2 雲台詳細規格 ... 6 表1-3 鏡頭詳細規格 ... 7 表2-1 定翼機與旋翼機比較 ... 13 表3-1 SIFT、PCA-SIFT、SURF 演算法比較表[26] ... 28符號定義
𝑡 時間點，單位為一個影格。 𝑝(𝑥, 𝑦, 𝑧) 三維空間中一點座標。 𝑝̇(𝑢̇, 𝑣̇) 𝑅 t 𝐶 二維影像平面上一點座標。 旋轉矩陣參數。 外部位移參數。 相機光學中心。 𝑓 𝑢𝑐、𝑣𝑐 𝜊 𝑓𝑢、𝑓𝑣 𝛾 𝑅𝑐 𝑡𝑐 𝐾 𝜋 相機焦距。 影像平面中心座標。 相機中心點。 相機在影像平面X軸與Y軸的等效焦距。 影像平面中X軸與Y軸夾角偏移程度。 相機外部旋轉參數。 相機外部位移參數。 相機內部參數矩陣。 投影矩陣，包含內部與外部參數。 𝑝̃(𝑢̃, 𝑣̃) 點 𝑝̇(𝑢̇, 𝑣̇)投影在正規化影像平面的位置。 𝑝̆(𝑢̆, 𝑣̆) 𝑟2 𝜅𝑖、𝛾𝑖 𝑝̈(𝑢̈, 𝑣̈) 𝐻̇ 𝐼𝐿、𝐼𝑅 𝑝𝐿、𝑝𝑅 點 𝑝̃(𝑢̃, 𝑣̃) 經過逕向形變修正之後的位置。 逕向距離平方。 形變係數。 點𝑝̇(𝑢̇, 𝑣̇)經過逕向形變修正之後的位置。 單應性矩陣，包含內部與外部參數。 立體視覺中左影像與右影像。 立體視覺中來自𝑝(𝑥, 𝑦, 𝑧)投影在𝐼𝐿, 𝐼𝑅的位置。𝐶𝐿、𝐶𝑅 𝑒𝐿、𝑒𝑅 𝐼𝛴(𝑥, 𝑦) 𝐻 Lxx、Lyy、Lxy 𝜎 𝑔(𝜎) 𝐼(𝑥, 𝑦) Dxx、Dyy、Dxy 𝐻𝑎𝑝𝑝𝑟𝑜𝑥 𝜔 𝑉𝑓𝑖 𝑑 𝐶𝑀𝐼(𝑝,̇ 𝑑) 𝑚𝑖 𝐼𝑏, 𝐼𝑚 𝐼𝑏𝑝̇, 𝐼𝑚𝑞̇ 𝑓𝐷 𝑞̇ 𝑒𝑏𝑚 ℎ𝐼(𝑖) 𝑛𝑐 𝑃𝐼(𝑖) 𝑔(𝑖) 𝐷 左右相機之光學中心 極線幾何中左右影像的極點位置。 積分影像在位置(x, y)的值。 海森矩陣。 點𝑥̇在尺度𝜎之高斯二階偏微分於 X 軸與 Y 軸對應的摺積。 特徵點偵測時的尺度。 尺度𝜎下之高斯函數。 影像I 在(x, y)位置的像素值。 X 軸與 Y 軸方向對應的盒子濾波器遮罩。 以盒濾波器遮罩取代。
Lxx, Lyy, Lxy 與 Dxx, Dyy, Dxy 之轉換對應係數。
特徵點𝑖之特徵向量。 歐幾里得位移量。 平面上一點𝑝̇在𝑑距離的 MI 匹配代價。 兩影像中，兩像素點之MI 值。 兩具有視差關係的平面影像。 兩影像中對應的像素點。 視差圖。 平面上一點座標，與𝑝̇有對應關係。 影像𝐼𝑏, 𝐼𝑚之共軛極線。 像素點𝑖在影像𝐼中之熵值。 像素點彼此對應的數量。 影像𝐼中像素點𝑖出現的機率。 像素點𝑖的高斯函數遮罩。 視差圖。
𝐷𝑝̇ 𝐸(𝐷) 𝐶(𝑝̇, 𝐷𝑝̇) 𝑃𝑘1, 𝑃𝑘2 𝑁𝑝̇ 𝐿𝑟𝑑(𝑝̇, 𝑑) 𝑆𝑟𝑑(𝑝̇, 𝑑) 𝐷̌ 𝐷𝑏𝑝̇, 𝐷𝑚𝑞̇ 𝐷𝑖𝑛𝑣 𝑝ℎ 𝑄𝑚 𝑏𝑙 𝑝𝑟 𝑝𝑟̇ 𝑝𝑟́ 𝑅 𝑡 𝑚𝑖,𝑡(𝑢, 𝑣, 𝑑)𝑇 𝑔̇(𝑖) ℎ(𝑖) 𝑒(𝑚𝑖,𝑡−1, 𝑚𝑖,𝑡) 視差圖中𝑝̇的視差值。 視差圖之能量。 平面上一點𝑝̇在𝑑距離的任一方式計算之匹配代價。 對視差變動的懲罰係數。 𝑝̇之鄰近點。 點𝑝̇在方向方向𝑟𝑑、距離𝑑的修正代價。 半全域的比對代價聚合。 代價聚合後之視差圖。 𝐼𝑏, 𝐼𝑚根據𝑆𝑟𝑑選擇視差後產生的視差圖像素。 無效的視差。 三維齊次座標。 用於視差估算深度之相機參數矩陣。 左右相機距離。 影像平面上特徵點根據相機內外部參數投影至三維空間 座標。 特徵點在影像中的二維座標。 特徵點三維座標根據相機內外部參數投影回影像平面之 座標。 特徵追蹤時之旋轉矩陣。 特徵追蹤時之位移向量。 時間點𝑡，左影像上被追蹤到的特徵點集合。 在歐基里德平面上使用三角測量計算特徵點𝑖三維座標。 三維空間上追蹤到的特徵點集合投影至影像。 三維空間中追蹤到的特徵點集合投影至影像上並與影像 上特徵點集合做相減。
𝑤𝑖,𝑗 𝜀𝐹 𝑚𝑖,𝑡 ̅̅̅̅̅ 𝛼, 𝛽 𝑎𝑖 𝜀𝑖𝑛𝑛𝑜 𝜀𝑇 𝑅𝑜𝑝𝑡 𝑇𝑜𝑝𝑡 𝑛𝑝 ⃑⃑⃑⃑ 𝑛𝑐 ⃑⃑⃑⃑ 𝑛𝑙 ⃑⃑⃑ 𝛿(𝑛⃑ , 𝑝𝑖) ∅(𝑝𝑐,𝑖, 𝑝𝑙,𝑖) 𝑏́ 𝑠(𝛽) 𝜙𝐸(𝛽) 布林係數，觀察時間點𝑖的特徵點𝑗是否存在。 匹配關係之兩點與基礎矩陣運算後之誤差容許值。 到時間𝑡為止之前累積的特徵點之平均值。 重投影誤差權重係數。 觀察到整合特徵的次數。 特徵整合過程之閥值。 剔除離群資料所設立之閥值。 經L-M 演算法最佳化後之旋轉矩陣。 經L-M 演算法最佳化後之位移向量。 校正平面法向量。 相機座標系法向量。 雷射掃描儀座標系法向量。 空間上一點𝑝𝑖投影至法向量為𝑛⃑ 之距離。 相機座標系上之點𝑝𝑐,𝑖經旋轉、位移後與雷射掃描儀座標 系之對應點𝑝𝑙,𝑖的誤差。 位移向量。 最小平方曲線擬合問題，𝛽為最佳化參數。 使用最佳化參數𝛽產生的偏移量之平方和。
第一章 緒論
1.1 前言
科技與時俱進瞬息萬變，各種現代化技術蓬勃發展，其中電腦視覺技術儼然成為重 要之課題，目前已廣泛應用在各科技領域，並且深入人類的生活中，舉凡在電影工業、 遊戲創作、娛樂設施以及汽車、造船等高科技工業，乃至考古、博物館展覽數位典藏等 文化等等，由此可見電腦視覺(Computer Vision) [1]的重要性也因此而提升。如何讓系統 得知自己在環境中的所在位置與移動狀態，即電腦視覺的範疇，也就是本論文所提出的 核心重點：本體移動估測(EGO-Motion Estimation)。另外，在電腦視覺(Computer Vision) 中，生成物體三維資訊主要有兩類方法，一類為使用幾何建模軟體透過人機互動控制下 繪出物體三維模型，另一類是透過電腦視覺、影像處理或其他技術獲取一實體物件之形 體輪廓幾何形狀。前者在現今技術已經非常成熟並有非常多種類的商用軟體，例如： Maya、3D Max、AutoCAD 等…。後者物體三維重建(3D Reconstruction)是指利用二維投 影恢復物體三維資訊（形狀等）的數學過程和電腦技術，包括本論文提出的特徵點偵測、 比對，追蹤等技術。1.2 研究動機
在導航領域，里程計是指通過一定設備獲取執系統運動時的數據，並據此估計系統 隨時間的位置變化，如汽車上使用輪圈的周長乘上旋轉的圈數，即可知道走了多遠。傳 統的方法適用於很多基於輪子或履帶運動的載具，但並不適用於使用非標準移動方式的 移動系統，如步行機器人等。此外，由於輪子可能在地板上滑動，這會產生一個與正常 轉動相比不標準的移動，導致測程法存在誤差問題。當系統在不光滑表面上運動時，這 個誤差會變得更加複雜，測程法的誤差會隨著時間累積和混合，使得測得的結果越來越 不準確。有非常多的方法來計算系統所移動的軌跡，但在本論文研究所專注的即以電腦視覺 為基礎的視覺測程（Visual Odometry,VO），使用視覺測程所獲取的資訊不僅是距離， 而是得到整個系統的移動軌跡。 視覺測程法利用相機所擷取的序列圖像分析其內容，並估計行進的距離。視覺測程 法可以在任何運動方式和任意種類的平面上提高導航準確度。
1.3 研究目的
本論文研究之目的，在於改良傳統使用雙相機系統之移動估測技術，僅獨立使用一 顆相機鏡頭來計算出裝置移動軌跡，發展可靠的本體移動估測。除此之外，並且改良原 本車載系統只能局限於地面環境，所能獲取之資料有限，為了獲取更廣泛的環境資訊， 將相機搭載於四軸無人飛行器上，如此能大範圍錄製校園、街道、城鎮等影像，最後運 用在建立具有高重建精準度以及色彩紋理的三維點雲模型之目的。1.4 研究方法
為了更精簡獲取精確的本體移動估測結果，本研究使用單眼視覺(Monocular -Vision) 方式，並以單相機錄製影片資訊進行分析，再利用全球定位系統(Global Positioning System，GPS)資料比對評估其準確性。 單眼視覺本體移動估測部分，首先將掛載之相機做內部與外部參數之校正，獲得相 機參數與移除鏡頭畸變。第二步驟對所擷取之影像序列進行特徵偵測(Feature- Detection)， 因特徵點演算法具尺度與旋轉不變特性，所以當影格畫面移動造成的遠近、旋轉、角度 的不同時，依然可找到特徵。第三步驟，當相機系統移動時，影格畫面產生變化，此步 驟追蹤各時間特徵點的位置變化，對影像特徵點進行三角測量獲得特徵點的旋轉矩陣 R 與位移向量 T，如此便可算出特徵點實際在空間中的三維座標，估測時可算出相對之移 動距離。第四步驟進行特徵追蹤(Feature Tracking)，追蹤各時間點間特徵點位置變化， 再利用特徵點在影像中之移動關係與相機參數估測出移動路徑，最後以最佳化演算法對估測的路徑進行調整以獲得最佳結果 最後利用本體移動估測之結果，並且結合所錄製的影片材質，可運用在重建出高精 準度且具有色彩紋理之三維環境模型。
1.5 PHANTOM 3 PROFESSIONAL 四軸飛行器
PHANTOM 3 PROFESSIONAL 是大疆創新[2]公司所生產的一款無人機，透過四個 旋翼產生向下氣流能如直升機般停懸於空中，透過控制器，其飛行距離最遠達到 5000 公尺以上；使用15.2 伏特 4800 毫安時鋰聚合物電池，飛行時間為 13 分鐘，並內建 GPS 感應器可精準穩定停旋，圖1-1 為本論文所使用之無人飛行器： 圖1-1 PHANTOM 3 PROFESSIONAL 四軸飛行器 PHANTOM 3 PROFESSIONAL 四軸飛行器外部組成為機體與四個馬達、四個螺旋 槳組成，內部組成包含：飛行控制系統、GPS 及 GLONASS 雙衛星導航系统、ESC 電 子調速器，IMU 慣性測量單元模塊內置六軸陀螺儀及加速度計，實時檢測飛機運動過 程中的微小姿態變化並反饋給主控，主控可由此做出相應的姿態補償和調節，時刻保持穩定飛行。其他詳細相關規格皆列在表1-1 中
圖1-2 四軸飛行器之遙控器
PHANTOM 3 PROFESSIONAL 四軸飛行器使用 2.4GHZ 無線電遙控器(Radio- Control)控制其運動，遙控器與飛行器之間除了姿態的控制，並內建高解析度之圖像傳 輸系統，可即時回傳搭載於無人機上相機所錄製的影像，透過傳輸線連接智慧型手機或 平板電腦，即時觀看實驗錄製資訊，綜觀以上優點，使得四軸無人飛行器成為錄製大範 圍實驗影像之首選設備，本研究論文也選用此飛行器作為錄製實驗影像之器材。
表 1-1 PHANTOM 3 PROFESSIONAL 詳細規格[2] 重量 1280克 對角線尺寸(軸距) 350mm 最大上升速度 5 米/秒 最大下降速度 3 米/秒 旋停精度 垂直: +/- 0.1 米(超聲波工作範圍內); +/- 0.5 米；水平: +/- 1.5 米 最大飛行速度 16米/秒 飛行滯空時間 約13分鐘 電池資訊 15.2V 4480mAh
1.6 四旋翼飛行器搭載鏡頭
單眼視覺相機系統部分使用PHANTOM 3 PROFESSIONAL 四軸飛行本身所搭載的 雲台相機，並搭載感光元件 1/2.3”CMOS 鏡頭，以進行本體移動估測。 1.6.1 雲台系統 一體化相機系統，具備三軸雲台穩定，透過陀螺儀晶片補償因機身晃動產生畫面震 動，且畫面幀數每秒可達 60 FPS(Frame Per Second)，如圖 1-3，其他詳細規格見表 1-2。圖1-3 雲台系統 表 1-2 雲台詳細規格[2] 穩定系統 三軸穩定 可控制範圍 俯仰 -90° 至 +30° 控制模式 可透過遙控器控制俯仰 補償範圍 3-軸 (俯仰, 橫滾, 偏航) 1.6.2 相機系統 此鏡頭擁有一千兩百四十萬高畫素，如圖 1-4，比先前研究所採用車載鏡頭更精細 穩定，且拍攝範圍更廣，符合本研究大範圍錄製的特點，其他詳細規格見表 1-3。
圖1-4 鏡頭 表1-3 鏡頭詳細規格[2] 影像傳感器 1/2.3" CMOS 有效像素 1240 萬 (總 像素 1276 萬) 鏡頭 FOV 94° 20 mm (35 mm 格式等效) 焦距 f/2.8 對焦點無窮遠 電子快門速度 8 秒 至 1/8000 秒 錄影解析度 UHD: 4096x2160p 24/25, 3840x2160p 24/25/30 ISO 範圍 100-3200
1.7 論文架構
本論文研究架構分為六章，如以下敘述。 第一章為緒論，內容分成研究目的、研究方法及使用儀器介紹，最後是本文架構。 第二章為文獻探討，本章將分為兩部分對現有文獻進行討論。第一節為對本體移動 估測之技術演變進行探討；第二節則針對無人機之組成與架構來探究。 第三章為視覺測程，各小節分別提到本體移動流程架構、特徵偵測與特徵匹配演算 法推導，最後是本節實驗討論。 第四章為實驗結果，實作本研究最後成果。實際將相機掛載於無人機上錄製影像， 並利用第三章所提出的估測方法與演算法進行多組實驗並分析。 第五章為結論與未來研究方向。第二章 文獻探討
本論文研究目的在於以單眼視覺(Monocular Vision)空拍影像實現本體移動估測，將 單一相機掛載於無人飛行器錄製影像，並分析計算影像序列，最後得到高精確度移動估 測軌跡之技術。本章文獻將分為兩大重點，第一部分針對使用單眼視覺進行本體移動估 測之相關研究文獻；第二部份針對無人飛行器的組成架構技術相關文獻做探討。2.1
使用單眼視覺進行本體移動估測
本論文研究之目的，在於使用單眼視覺完成一個視覺測距系統。視覺測距系統除了 可以獨自應用外，亦可以輔助各種缺乏移動感測的裝置系統獲得移動路徑。 要如何讓系統知道自己在環境中的移動路徑? 如圖 2-1 所示，系統從 A 地移動至 B 地，感測裝置只知道周圍環境之影像有所變化，卻不清楚實際上本身之移動距離。因此， 本文使用單眼視覺方式，將連續時間點之影像做運算，取得相機在環境中的移動路徑， 而解決這類問題稱之為本體移動(Ego-motion)。圖2-1 本體移動概念圖
解決本體移動問題的技術，類似於電腦視覺中的移動恢復結構法(Structure from Motion, SfM)，稱之為本體移動估測(Ego-motion Estimation)或視覺測程(Visual Odometry,
VO) [3]，是一種使用不同時間點中的影像，進行影像匹配估測移動的技術。「本體移動」
一詞於 1994 年由 M. Irani 等人所定義[4]，利用二維影像計算出相機在三維空間環境中 的移動，以解決本體移動計算的問題。但早在 1980 年，美國太空總署(The National Aeronautics and Space Administration, NASA)已開始研究[5]，在火星探險任務中，火星探 測漫遊者號（Mars Exploration Rover, MER）[6]用來計算無人駕駛探測器在火星上的移 動路徑。另一種類似的估測技術為使用場景與當前影像的匹配來估測移動，希望系統從 未知環境的未知地點出發，在運動過程中通過重複觀測到的地圖特徵（比如，牆角，柱 子等）定位自身位置和姿態，再根據自身位置增量式的構建地圖，從而達到同時定位和 地圖構建的目的。稱之為視覺即時定位與地圖構建(Visual Simultaneous Localization and Mapping, V-SLAM)。大部分的學者們將本體移動估測與 V-SLAM 通稱為即時定位與地 圖構建(Simultaneous Localization and Mapping, SLAM)的問題研究。2006
年，Durrant-Whyte 等人對 SLAM 有詳細的介紹[7][8]；2010 年 B. Williams 等人也對本體移動估測 與 V-SLAM 作介紹並共同實作[9]。 本體移動估測主要重點步驟如圖2-2 所示。首先進行影像擷取和校正，接著同時對 影像進行特徵偵測與匹配計算得出特徵點在三維空間中的座標，並將偵測到的特徵點進 行匹配，再利用特徵點在影像中之移動關係與相機參數估測出移動路徑，最後進行最佳 化處理使路徑結果更為完善。 圖2-2 本體移動估測之概略圖 移動估測為本體移動估測中最重要的步驟，使用相機進行移動估測計算移動矩陣時， 根據特徵點屬性主要可分為三種方式：1.二維影像特徵點的不同時間相互對應關係(2D-2D)；2.三維結構與影像特徵點對應關係(3D-2D)；3.三維資料的對齊(3D-3D)。最早開始 的移動估測，為三維資料的對齊(3D-3D)，主要是使用三維資料進行對齊，其中迭代最 近點(Iterative Closest Point, ICP)演算法進行空間中座標點的對齊便為常使用的方法[10]； 2012 年，Lui 等人使用了微軟公司的體感深度感測機(Kinect)進行深度的擷取[11]，依然 使用到了 ICP 演算法來進行各時間點的空間座標點對齊的工作。而本論文為單眼視覺，
使用第一種方式(2D-2D)僅需使用一台相機在不同時間點的兩張影像進行估測。單眼視 覺測程是二維影像特徵點在不同時間相互對應的關係(2D-2D)，但透過此方法所計算出 的估測路徑缺乏實際距離單位，須以影格序列中的最初前兩張做為相對的參考標準。 而對於三維空間與二維影像對應的估測方式，主要是 2004 年的 D. Nister 等人提出 VO 的概念並提出相關可能進行的研究後[3]，才逐漸的興起。2007 年的 A. Comport 等 人將這種基於二維影像與三維空間點進行估測的方式，視為新的技術且發表了對於移動 估測的研究[12]。此外，越來越多人使用各種儀器結合電腦視覺技術的方式來進行估測， 如 2011 年的 L. Wei 等人使用全球定位系統(Global Positioning System, GPS)結合影像進 行車輛於城市中的移動估測與地圖建構[13]。 2010 年 D. Scaramuzza 等人提出了一份完整的研究報告，統整了近年來對於本體移 動估測的研究與討論[14][15]；2013 年 H. Badino 提出了一種抗雜訊的 VO 特徵點追蹤： MFI[16]演算法，它強調誤差不隨時間增加而累積變大，為當時 KITTI[17]上純視覺 VO 最好的演算法，因此，本研究在第三章將採用其方法來進行本體移動估測。此外 2013 年 P. S. Huang 以限制視窗方式剔除不必要的特徵點，並實作本體移動估測[18]；2014 年 J. H. Zhang 以適應性方法增強特徵點匹配與追蹤正確性，提高本體移動估測之精確 度[19]。
2.2
無人飛行器之發展
無人飛行載具（Unmanned Aerial Vehicles,UAV）[20]或者稱之無人飛機，是近幾年 來航空產業極度熱們的話題。之所以稱為無人飛行器，即是沒有飛行員在載具上控制， 而是使用其他方式執行任務(如遙控或自動飛行)。最早期的無人飛行載具主要是以軍事 為主，如靶機、無人偵查機、無人攻擊機。而因現在科技發展迅速，電子陀螺儀、GPS 、 遙控無線電和飛行控制系統的發展，使無人飛行載具成長迅速。民間也開始流行無人飛 機，用於娛樂、攝影等領域。 無人飛行載具一般可分為「定翼型」(Fixed-wing aeroplane)與「旋翼型」(Autogyro、
Gyroplane)兩大類，定翼型飛行器即大眾熟知的飛機，外型通常是機身兩側有一對固定 的機翼，當飛行於空中時會產生浮力滑翔於空中，飛行方向只能向前進，無法定點於空 中；另一旋翼型飛行器即以直升機為主要外型，機身含有至少一組主旋翼，主旋翼向下 產生推力，且能藉由改變旋翼速度、角度來達到姿態的改變，並且能定點停懸於空中。 近幾年因科技的進步，旋翼型飛行器已從單旋翼直升機發展成多旋翼(Multicopter)， 如四旋翼(Quadcopter)、六旋翼、八旋翼等。因對稱的旋翼可相互抵銷力矩，且加上飛控 陀螺儀的輔助控制，因此飛行控制相對簡單許多，大大降低操作門檻，兩大類飛行載具 之比較如表所示。 表 2-1 定翼機與旋翼機比較 由於多軸飛行器門檻低、結構簡單、攜帶方便穩定等特性，被接受度越來越大，應 用層面也越來越廣。 本論文四軸空拍飛行器的操作單元組成一般包含：機體系統、飛行控制系統、穩定 雲台系統、影像擷取設備、遙控設備、圖像傳輸與監視設備，以上是可執行任務的最低
要求。 各單元設備詳細如下： A. 機體系統 機體組成包掛機身、四顆無刷馬達、電子調速器、電池、兩對正反螺旋槳，因考 量重量問題，多數零件為塑膠製品。 B. 飛行控制系統 多軸飛行器之運動姿態主要依靠馬達的轉速改變來達成，因此需要使用飛控系 統來協助控制，在電路板中整合陀螺儀、加速度計、氣壓計、GPS、 運算核心進而達到 GPS 定點定位、自動導航、自動返航等功能。 C. 穩定雲台系統 空拍機顧名思義就是將相機掛載於飛行器上，但因機體姿態不斷改變與飛行時 馬達的震動，造成畫面晃動模糊不清，將相機裝載於雲台上即可大幅改善此狀 況得到穩定的畫面。本論文所使用的雲台為三軸雲台，具備三個無刷馬台穩定， 透過陀螺儀晶片補償因機身晃動產生畫面震動 D. 影像擷取設備 影像擷取設備主要負責空中攝影，一般常見為掛載單眼相機、類單眼相機、運動 相機等。另外依照專業需求，如勘災、電力檢測等也可掛載紅外線相機，熱影像 相機、多光譜儀、雷射測距儀等。本論文僅使用單一相機進行實驗估測，且最高 畫素可達到 4K(4096x3840)。 E. 遙控設備 遙控設備為操控飛行員與無人機之間的橋樑，包含遙控器與無人機上的接受器， 常見使用的頻率以 2.4Ghz 為主，其控制距離最遠可達 5 公里。 F. 圖像傳輸與監視設備 因執行飛行任務時常飛出視距外，透過圖像傳輸與監視設備(圖傳系統)即可實時 觀看飛行器之狀態，獲取飛行資訊，改變飛行路徑，對實驗所需之路徑取景。
第三章 視覺測程
本章主要在於介紹單眼視覺感測系統之架構及實作流程。此章節之規劃如下，3.1 節 簡 介 單 眼 視 覺 (Monocular Vision) 、 相 機 模 型 (Camera Model) 、 相 機 校 正 (Camera Calibration)及 VO 實作流程(Visual Odometry)，並說明本章之實作流程架構；3.2 節講解 本研究之特徵點偵測(Feature Detection)；3.3 節說明特徵匹配匹配(Features Matching)； 3.4 節為離群資料剔除與修正；3.5 節說明移動估測實作，最後 3.6 節為估測路徑之最佳 化。
3.1
前言
本論文主要研究基於電腦視覺的視覺測程 Visual Odometry (VO)，視覺測程是透過 單相機、雙相機或多個相機擷取影像序列，並透過分析關聯的圖像來測定系統(如機器 人、車輛、儀器)的位置和方向的過程。視覺測程目前被廣泛應用在機器人系統、擴增實 境、汽車導航、火星探測漫遊者號[5]。一般實現視覺測程的方法較常見的為雙眼視覺法； 雙眼視覺意指生物在雙眼視野範圍互相重疊之下，所產生的視覺，這種視覺信號經傳送 至大腦後，即可判斷眼睛到物體之間的精準距離關係。類似的觀念也應用於電腦領域中， 使用兩台平行放置的相機對物件進行拍攝，使電腦與人類一樣有能力判斷眼前物體深遠 的程度，進而產生對三維立體的認知，就是所謂的電腦立體視覺(Computer Stereo Vision)， 根 據 兩 台 相 機 之 相 對 位 置 及 影 像 之 中 重 疊 部 分 進 行 兩 影 像 之 像 素 對 應 分 析 (Correspondence Analysis)，推算影像中物件深度資訊以達到景深的計算[21]。相較於雙 眼視覺方法需使用兩台攝影機來擷取影像，單眼視覺方法僅使用一台攝影機的單眼視覺 法可提升便利性與降低成本，故本論文專注於單眼視覺測程之研究。在單眼視覺法中， 相對運動與 3D 結構都必須從 2D 的序列資料計算而來，且一開始的絕對尺度(Absolute- Scale)是未知的，必須經由最初所輸入的前兩個影格(frame)經由三角測量(Trifocal Tensor) 計算得到相對尺度。單眼視覺測程可分為三個類別：以特徵點為基礎、以外觀為基礎、 或者混合以上兩種方法。以特徵為基礎的方法中，透過追蹤影格中顯著的特徵點並找出
特徵點的描述向量，在空間中建立相對關係；而以外觀為基礎的方法透過檢測每一個像 素(Pixels)，這將會需要龐大的運算，所以本論文採用基於特徵點為基礎的視覺測程。 3.1.1 相機模型(Camera Model) 所謂的相機模型是使用幾何的方式描述真實世界中物體在相機影像平面上的成像 方式。該成像之投影關係可由相機的幾何參數描述，這些參數又分為以下參數： 內部參數(Intrinsic Parameters)與外部參數(ExtrinsicParameters)。相機之內部參數包含： 鏡頭焦距(Focal Length)、影像主軸點(Principle Point)、像素長寬比(Pixel Aspect Ratio)與 稜鏡扭曲(Lens Distortion)；而外部參數則是用來描述世界座標(World coordinate)和相機 座標(Camera coordinate)之間的座標系統轉換關係，通常可透過一組旋轉矩陣(Rotation Matrix)與位移向量(Translation Vector)將之描述，通常表示為 R 與 t，如圖 3-1 所示： 圖3-1 內、外部參數轉換關係圖 其中，C為相機光學中心，𝑝(𝑥, 𝑦, 𝑧)為空間上一點座標，𝑝̇(𝑥, 𝑦)為影像平面上一點座標， (𝑅, 𝑡)為兩座標系統間之旋轉、位移參數。 一般來說，物體直接放置於成像平面前，會因光線漫射的緣故，使得物體無法在影 像平面上產生清晰影像，如圖 3-2 (a)所示，當光線擊中物體上之 A 點後，一部分被物體 吸收外，其餘光線會朝四面八方漫射，使得 A 點無法在成像在平面上形成唯一的對應
光點，導致成像一片模糊。 (a) (b) 圖3-2 針孔相機原理圖 (a)一般情況成像 (b)加上隔板後之成像 解決此問題只需在物體和成像平面之間放一塊隔板，隔板上有個極小的孔洞，如圖 3-2 (b)所示，孔洞可以過濾從 A 點漫射出來的光線，使得 A 點在成像平面上只有唯一的一 個對應光點，此時的影像是最清晰的；反之，孔洞越大，越多光線可通過，影像也越模 糊。 一般數位相機成像就是利用針孔相機(Pinhole Camera)原理，透過針孔成像原理將三 維空間的座標點投影至二維影像平面，如圖 3-3 為針孔相機模型的示意圖，與圖 3-2 不 同的是成像平面移到物體前方，成像平面上的影像也不會上下顛倒。
圖3-3 相機針孔模型 假設 𝑓 為相機之焦距，且在空間中有一點 𝑝 ，𝑝之三維座標為 (𝑥, 𝑦, 𝑧) 點𝑝經過 透視投影(Perspective Projection)至影像平面後，投射在影像空間上的一點 𝑝̇(𝑢̇, 𝑣̇) ，其 中𝑢̇代表影像中 X 軸方向座標、𝑣̇代表 Y 軸方向座標。其投影相關之矩陣可由(3.1)式所 示： 𝑝̇ ( 𝑢̇ 𝑣̇ 1 ) ~ ( 𝑓𝑢 𝛾 𝑢𝑐 0 𝑓𝑣 𝑣𝑐 0 0 1 ) (𝑅𝑐 𝑡𝑐)3×4( 𝑥 𝑦 𝑧 1 ) (3.1)
𝑓𝑢與𝑓𝑣為 X 軸與 Y 軸之等效焦距(Effective Focal Lengths)，𝛾則代表影像平面中 X
軸與 Y 軸之夾角於相對直角的偏移程度，通常忽略為 0。𝑢𝑐與𝑣𝑐為影像平面中心座標， (𝑅𝑐 𝑡𝑐)為相機之外部參數，分別為𝑅𝑐3×3旋轉矩陣與𝑇𝑐3×1位移矩陣。將(3.1)整理後可得 (3.2)，其中的𝐾即為內部參數矩陣，而 𝜋 為包含內部參數及外部參數之 3×4 投影矩陣 (Projection Matrix)，~ 表示左右向量不考慮長度時相等。 𝑝̇ ( 𝑢 𝑣̇ 1 ̇ ) ~ 𝐾(𝑅𝑐 𝑡𝑐)3×4( 𝑥 𝑦 𝑧 1 ) ~ 𝜋 ( 𝑥 𝑦 𝑧 1 ) (3.2)
針孔相機的缺點是經過針孔的光線量少，因此需要較長的曝光時間；現代的鏡頭則 搭配鏡片(Lens)，它具有選擇光線的能力，且其直徑大於針孔，單位時間可讓更多光線 進入，大幅降低曝光時間，但使用鏡片會造成影像或多或少出現徑向形變 (Radial Distortion)的問題。鏡片是一種透鏡，透鏡存在曲率，徑向形變是由於光線穿過鏡片時， 受稜鏡曲率影響，造成光線折射，產生影像變形之現象，越到邊緣，就越嚴重，徑向形 變又分為枕狀形變(Pincushion Distortion)和桶狀形變(Barrel Distortion)，正常的影像如(a)， 如圖 3-4[22]示： (a) (b) (c) 圖3-4 徑向形變示意圖[22] (a) 正常的影像 (b)枕狀形變 (c)桶狀形變 要解決徑向形變問題，首先將點 𝑝̇(𝑢̇, 𝑣̇) 投影在正規化影像平面(Normalized Image Plane)為點 𝑝̃(𝑢̃, 𝑣̃) ，如圖 3-5 所示：
圖3-5 包含徑向形變之相機投影模型 點 𝑝̃(𝑢̃, 𝑣̃) 可由(3.3)式所得： ( 𝑢̃ 𝑣̃ 1 ) ~ ( 1 0 0 0 1 0 0 0 1 0 0 0 ) (𝑅𝑐 𝑡𝑐 0 1) ( 𝑥 𝑦 𝑧 1 ) (3.3) 校正形變後的點 𝑝̆(𝑢̆, 𝑣̆) ，則可由(3.4)式計算： ( 𝑢̆ 𝑣̆ 1 ) ~ ( 𝑢̃ 2𝑢̃𝑣̃ 𝑟2 + 2𝑢̃2 𝑣̃ 𝑟2+ 2𝑣̃2 2𝑢̃𝑣̃ 0 0 0 0 0 1 ) ( 1 + 𝜅1𝑟2+ 𝜅2𝑟4+ 𝜅3𝑟6 𝛾1 𝛾2 1 ) (3.4) 其中， 𝑟2 = 𝑢̃2+ 𝑣̃2 為徑向形變(Radial Distance)之平方，即理想點與影像中心點之距 離，而 (𝜅1, 𝜅2, 𝜅3, … , 𝛾1, 𝛾2) 為形變係數(Distortion Coefficients)。最後點 𝑝̆(𝑢̆, 𝑣̆) 可投 影回影像平面，得到徑向形變修正轉換後之點𝑝̈ (𝑢̈, 𝑣̈) ，如(3.5)式： Center of projection Image Plane 𝑥𝑤 𝑦𝑐 𝑧𝑐 𝑥𝑐 (0, 0) (0, 0)
Normalized Image Plane
𝑧𝑤 𝑦𝑤 World 𝑧 = 1 𝑧 = 𝑓 𝑝(𝑥, 𝑦, 𝑧) 𝑝̇ 𝑝̃ 𝑝Ƽ 𝑝̈
( 𝑢̈ 𝑣̈ 1 ) ~ ( 𝑓𝑢 𝛾 𝑢0 0 𝑓𝑣 𝑣0 0 0 1 ) ( 𝑢̆ 𝑣̆ 1 ) (3.5) 3.1.2 相機校正 本論文將操控無人機實地拍攝街道景色取得影像資料。本研究實驗所分析之影像序 列皆完成校正。 首先進行相機校正與影像擷取，本研究採用 1280*720 解析度之彩色影像，並以每 秒 30FPS 之速度來擷取影像，基本上拍攝前後之流程大致上如圖 3-6 所示： 圖3-6 影像擷取前後工作流程 本研究使用 Zhang 式相機校正法進行校正[23]。使用此方法的原因在於校正方法的 彈性與校正物件容易取得，使用之校正樣板如圖 3-7 所示。 首先，必須先算出相機與校正平面之間的轉換關係，我們定義校正平面法向量為
𝑛𝑝 ⃑⃑⃑⃑ = (0,0, −1)，如圖 3-8 所示。 圖3-7 Zhang 式方法使用之校正板 圖3-8 校正平面與其座標系 校正後可得到相機與校正平面之外部參數𝑅𝑐3×3與𝑇𝑐3×1，接著將外部參數代入相機 投影模型，並推導出相機與校正平面投影模型，如(3.6)式：
( 𝑢 𝑣 1 ) ~ ( 𝑓𝑢 0 𝑢𝑐 0 𝑓𝑣 𝑣𝑐 0 0 1 ) (𝑅𝑐 𝑇𝑐)3×4( 𝑝𝑥 𝑝𝑦 𝑝𝑧 1 ) (3.6) 其中𝑓𝑢和𝑓𝑣分別為水平與垂直方向焦距，(𝑢𝑐, 𝑣𝑐)為影像中心點，(𝑝𝑥, 𝑝𝑦, 𝑝𝑧)為校正 平面上之點，(𝑢, 𝑣)則是(𝑝𝑥, 𝑝𝑦, 𝑝𝑧)投影到相機平面上的點座標。 校正平面法向量𝑛⃑⃑⃑⃑ 經過外部參數𝐸𝑝 𝑐 = (𝑅𝑐, 𝑇𝑐)轉換到相機座標系中之平面法向量𝑛⃑⃑⃑⃑ ，𝑐 再利用齊次座標系(homogeneous coordinate system)將整個問題都向上增加一維度進行推
導，以利用算。因此𝑛⃑⃑⃑⃑ 可改寫成𝑛𝑝 𝑝 = ( 𝑛𝑝,𝑥 ‖𝑛⃑⃑⃑⃑⃑ ‖𝑝 , 𝑛𝑝,𝑦 ‖𝑛⃑⃑⃑⃑⃑ ‖𝑝 , 𝑛𝑝,𝑧 ‖𝑛⃑⃑⃑⃑⃑ ‖𝑝 , −‖𝑛⃑⃑⃑⃑ ‖)，其上有一點P = (𝑝 𝑝𝑥 𝑝𝑦 𝑝𝑧 1 )， 並使用向量內積來表示，可將點𝑝在法向量𝑛𝑝平面上寫成(3.7)式： 𝑛𝑝𝑡𝑝 = 𝑛 𝑝 𝑡𝛪𝑝 = 𝑛 𝑝 𝑡(𝐸 𝑐−1𝐸𝑐)𝑝 = 0 (3.7) 其中外部參數𝐸𝑐 = (𝑅𝑐 𝑇𝑐 0 1)，點 p 轉換到法向量𝑛⃑⃑⃑⃑ 平面上的點為p𝑐 ′，得(3.8 )式： (𝑛𝑝𝑡𝐸𝑐−1)(𝐸𝑐𝑝) = (𝐄𝐜−𝐭𝐧𝑝)t𝑝′= 0 = 𝐧𝐜t𝑝′ (3.8) 根據(3.8)式可得法向量𝑛⃑⃑⃑⃑ 轉換到法向量𝑛𝑝 ⃑⃑⃑⃑ 可表示成(3.9)式： 𝑐 nc = Ec−tnp = ( R−1c 0 −R−1c Tc 1) np (3.9)
(3.9)的結果需再轉回卡式座標(Cartesian coordinate system)，卡式座標系統中之法向量𝑛⃑⃑⃑⃑ 𝑝
與平面法向量𝑛⃑⃑⃑⃑ 轉換關係式表示為(3.10)式： 𝑐 ω(Ec, 𝑛⃑⃑⃑⃑ ) = n𝑝 ⃑⃑⃑ =c Rc𝑛⃑⃑⃑⃑ 𝑝 ‖𝑛⃑⃑⃑⃑ ‖𝑝 [(Rc𝑛⃑⃑⃑⃑ + T𝑝 c) Rc𝑛⃑⃑⃑⃑ 𝑝 ‖𝑛⃑⃑⃑⃑ ‖𝑝 ] (3.10)
ω函數代表原法向量經過剛性轉換後得到的新法向量。經過(3.10)式轉換後之法向量即 為後續剛性轉換中相機的校正資料。
Zhang 方法的校正物體為圖 3-7 之校正板，而成像所在的影像也是一個平面物件， 將空間中某平面上的一點轉移到另一平面上，在電腦視覺領域稱此轉移為平面單應轉移 (Planar Homography Transform)，圖 3-9 為其示意圖：
圖3-9 單應平面示意圖 若校正物體為一平面物件，可將(3.1)式簡化為(3.11)式： 𝑝̇ ( 𝑢̇ 𝑣̇ 1 ) ~ ( 𝑓𝑢 𝛾 𝑢𝑐 0 𝑓𝑣 𝑣𝑐 0 0 1 ) (𝑅𝑐, 𝑡𝑐) ( 𝑥 𝑦 1 ) (3.11) 其中𝐻̇包含相機內部參數、旋轉矩陣和位移向量，又可將(3.11)式改為(3.12)式： 𝑝̇ ( 𝑢̇ 𝑣̇ 1 ) ~𝐻̇ ( 𝑥 𝑦 1 ) (3.12)
Zhang 方法以兩種以上不同角度對校正板進行拍攝，建立與拍攝數量等量的單應性 矩陣𝐻̇(Homography Matrix)，最後再進行最佳化參數的動作，找出最適合目前相機之內 外部參數，圖 3-10 為本實驗校正拍攝的六十四個不同角度校正版，拍攝校正影像過程 中，校正版需盡量均勻遍布整個影像平面，使校正能夠更準確，如此估測結果才能準確。 圖3-10 六十四個不同角度遠近所拍攝之校正板影像 圖 3-11 為校正後之外部參數計算出校正位置，其彩色方格為二十個校正版相對於 相機之位置。
圖3-11 校正外部參數所得之二十個校正版相對於相機之位置
3.1.3 VO 實作流程
完整的視覺測程(Visual Odometry)如圖 3-12 所示。首先進行相機校正，獲得相機內 部與外部參數並移除鏡頭畸變，接著輸入影像序列的最初兩張影格，對最初前兩張影格 做特徵偵測並追蹤相應的特徵點，透過特徵點計算得到基礎矩陣，接著從基礎矩陣分解 出旋轉矩陣(Rotation Matrix)與位移向量(Translation Vector)，並從以上資訊或的尺度， 串聯得到估測路徑。重複以上步驟直到處理完所有的影像序列，最後最佳化之本體移動 資訊。
圖3-12 視覺測程(Visual Odometry)架構圖
3.2
特徵點偵測
本論文研究使用加速穩健特徵演算法(Speeded Up Robust Features, SURF)進行特徵 點偵測之演算。為了改善加速尺度不便特徵轉換演算法 (Scale-invariant feature transform, SIFT)的執行速度與減少特徵點敘述子的向量維度，SURF 演算法由 Herbert Bay[24]等人 在 2006 年發表，其概念及步驟均建立在尺度不變特徵轉換演算法(Scale-invariant feature transform, SIFT)[25]之上，在多幅圖片下更加穩健且有效率，是現今電腦視覺中用來偵 測及描述區域特徵之代表性演算法。
2009 年， Juan 和 Gwon 做出了對 SIFT 演算法與 SURF 演算法的比較[26]，其中包 含 SIFT 的改良：主成份分析－尺度不變特徵轉換演算法(Principal components analysis-Scale-Invariant Feature Transform, PCA-SIFT)，結果如表 3-1 所示：
表3-1 SIFT、PCA-SIFT、SURF 演算法比較表[26]
Method Time Scale Rotation Blur Illumination Affine SIFT common best best common common good
PCA-SIFT good good good best good Good
SURF best common common good best Good
單眼視覺系統下，尺度(Scale)與旋轉(Rotation)之變化較少，所以不考慮使用 SIFT； 而因為 SURF 在速度及其他方面表現較為良好，故本論文研究採用 SURF 作為特徵偵測 演算法。
SURF 是以積分影像(Integral Image)加速運算，再用可改變尺度之海森矩陣(Hessian Detector)快速地偵測特徵，再以哈爾(Haar)濾波器求出主方向，將子區域旋轉至同一方 向對齊後建立 SURF 敘述子(Descriptor)。 3.2.1 積分影像(Integral Image) 積分影像[27]之目的在於加速矩形遮罩的摺積(Convolution)運算，當計算影像中相 近的區域時，不需整個區域重復再進行計算，只需取出四個角落值，就可以得到任何模 板內之影像值。由影像要做積分後才能用於海森矩陣的計算，若沒有積分影像，則需進 行相當耗時的連加運算，造成運算浪費。因此使用同一張影像，只需計算一次的積分， 矩形遮罩之摺積運算則可重複使用區塊像素累加成果，藉此大幅降低運算量。如圖 3-15 所示，其中圖 3-13(a)代表從左上點往右下點所延伸之一組長方形藍色面積，藍色面積內 的灰階總和紀錄在右下紅點，其計算公式如(3.13)式： 𝐼∑(x, y) = ∑ ∑ 𝐼(𝑖, 𝑗) 𝑗≤𝑦 𝑗=0 𝑖≤𝑥 𝑖=0 (3.13)
(a) (b) (c) 圖3-13 積分影像示意圖 (a)影像積分概念 (b)影像灰階值 (c)影像積分值 一旦積分影像建立後，就可以直接利用查表的方式，快速找出任意矩形區域之像素 總和。如圖 3-16 所示，其中𝐼∑1 = ∑𝐵𝑙1、𝐼∑2 = ∑(𝐵𝑙1+ 𝐵𝑙2)、𝐼∑3= ∑(𝐵𝑙1+ 𝐵𝑙3)和𝐼∑4 = ∑(𝐵𝑙1+ 𝐵𝑙2+ 𝐵𝑙3+ 𝐵𝑙4)，若要計算區塊𝐵𝑙4的像素總和，可藉由積分影像中位置的值𝐼∑1、 𝐼∑2、𝐼∑3和𝐼∑4得出：區塊𝐵𝑙4像素總和 = 𝐼∑4 − 𝐼∑3− 𝐼∑2+ 𝐼∑1。 圖3-14 利用積分影像計算矩形區域像素總和
3.2.2 快速海森特徵檢測
在 SURF 中海森矩陣(Hessian Matrix)是核心 ，其用途在於偵測特徵點。海森矩陣是 由二階偏導數組成，我們可以用海森矩陣的行列式值來判斷點(x, y)是否為極值點。假設 有一函數𝑓(𝑥, 𝑦)，圖像中某個像素點的海森矩陣定義為如下(3.14)： 𝐻(𝑓(𝑥, 𝑦)) = | 𝜕2𝑓 𝜕𝑥2 𝜕2𝑓 𝜕𝑦𝜕𝑥 𝜕2𝑓 𝜕𝑦𝜕𝑥 𝜕2𝑓 𝜕𝑦2 | det(𝐻) = 𝜕 2𝑓 𝜕𝑥2 𝜕2𝑓 𝜕𝑦2− ( 𝜕2𝑓 𝜕𝑦𝜕𝑥) 2 (3.14) 行列式中， det(𝐻) > 0且𝜕2𝑓 𝜕𝑥2 > 0，𝑓(𝑥, 𝑦)是局部極小點；det(𝐻) > 0且 𝜕2𝑓 𝜕𝑥2 < 0， 𝑓(𝑥, 𝑦)是局部極大點；det(𝐻) < 0時則𝑓(𝑥, 𝑦)是鞍點。鞍點指函數在該點時，𝑥軸方向與 y軸方向的一階導數為零且兩方向的二階導數正負號相反；相對極值點表示該點一階導 數為零且兩方向的二階導數正負號相同。因此，凡是區域的相對極值點，其海森矩陣行 列式值必大於零。利用此一特性，可在不同尺度空間上計算海森矩陣，在尺度𝜎下，點 𝑝 = (𝑥, 𝑦)處之海森矩陣如下所示： 𝐻(𝑝, 𝜎) = [𝐿𝑋𝑋(𝑝, 𝜎) 𝐿𝑋𝑌(𝑝, 𝜎) 𝐿𝑋𝑌(𝑝, 𝜎) 𝐿𝑌𝑌(𝑝, 𝜎)] (3.15) 將高斯函數的二階偏微分與影像𝐼中之點𝑝(𝑥, 𝑦)處之像素卷積的結果如下： 𝐿𝑋𝑋(𝑝, 𝜎) = 𝜕2 𝜕𝑥2𝑔(𝜎)⨂𝐼(𝑥, 𝑦) , 𝑔(𝜎)為高斯函數 (3.16) 同理𝐿𝑋𝑌、𝐿𝑌𝑌也是如此計算，可算出圖片上所有點的海森行列式值。為加速運算，
藉由積分影像可減少摺積計算時間，SURF 演算法以盒濾波器(Box Filter)將高斯二階微 分做近似處理，以DXX DYY與𝐷𝑋𝑌 取代海森矩陣的元素，如圖3-15 所示，左上𝐿𝑌𝑌是𝑦方 向的二階高斯微分模板，右上𝐷𝑌𝑌是採用盒濾波器的結果，數字表示對應顏色之權重值， 灰色區域為 0，其行列式近似為下式(3.17)： 圖 3-15 盒濾波器替代高斯二階微分運算 det(𝐻𝑎𝑝𝑝𝑟𝑜𝑥) = 𝐷𝑥𝑥𝐷𝑦𝑦− (𝜔𝐷𝑥𝑦)2 (3.17) 其中𝜔為方程式平衡的權重值，由於不同尺度下，𝜔對應的值不同，若以 9×9 的盒 濾波器進行偵測，會與高斯二階導數選用𝜎 = 1.2相似，𝜔將近似於 0.9，如式(3.18)： ω = |𝐿𝑥𝑦(1.2)|𝐹|𝐷𝑦𝑦(9)|𝐹 |𝐿𝑦𝑦(1.2)|𝐹|𝐷𝑥𝑦(9)|𝐹 = 0.912 ≈ 0.9 (3.18) SIFT 建立尺度空間是以原始影像作為最底層，固定高斯遮罩，原始影像尺寸會以 金字塔之方式縮小，如此容易造成影像失真，且每層建立必須等到上層構造完畢後才可 進行處理，彼此依賴性強，相當耗時，如圖3-16(a)所示；SURF 則是採用原始影像尺寸
不變，只改變盒濾波器的大小，並用不同大小之遮罩對原圖進行處理，可上下並行運算，
且每層影像同時處理，較有效率，如圖3-16(b)[28]所示：
(a) (b)
圖3-16 尺度空間構造示意圖[28]
(a) SIFT 尺度空間 (b) SURF 尺度空間
SURF 採用9 × 9的盒濾波器(近似於𝜎 = 1.2的高斯濾波器)對原圖做處理，產生之影 像當成最底層的尺度，越往上層濾波器大小也跟著增加，15 × 15、21 × 21、27× 27… 等，如此可產生多組尺度空間，以獲取更多特徵，過程中使用線性內插法，降低階(Octave) 和階之間的尺度變化，避免尺度取樣因變化大而過於粗糙。 獲得尺度空間後，接著就從這些尺度中搜尋可能的特徵點，也就是搜尋空間中的極 值(Extrema)，每層影像上的每個點對其範圍內之對應點(同層有 8 個，上、下兩層尺度 空間有2 × 9 = 18個，共 26 個點)做海森矩陣行列式值之比較，如圖 3-17，若該點為鄰 近26 個點之最大值或最小值，即判定為特徵點。
圖3-17 特徵點定位示意圖 3.2.3 SURF 敘述子 找到特徵點後，則需要對特徵點進行描述，為保持旋轉不變特性，須找到特徵點之 主方向(Major Orientation)。首先以特徵點為圓心，半徑大小 6𝜎(𝜎為特徵點所在尺度值) 之圓形區域內，對每一點進行哈爾(Haar)小波遮罩運算(小波邊長為 4𝜎)，實際上就是計 算圓域內每像素之水平與垂直的一階導數值，如圖 3-18 所示。 圖3-18 哈爾(Haar)小波濾波器 再以圓心角60°的扇形視窗環繞特徵點掃描一圈，統計扇形視窗內一階導數值(𝑑𝑥 與 𝑑𝑦)的加總，擁有最大加總值的掃瞄角度作為特徵點的主要方向，如圖 3-21 所示，
斑點為哈爾小波遮罩運算一階導數值，灰色為扇形掃描視窗，箭頭長度為經過加總後的 大小，箭頭方向則表示特徵點方向。 圖3-19 主方向示意圖 前幾小節所描述的步驟完成後，特徵點應已被賦予位置、尺度及主方向，其周圍像 素點需以此主方向為基準來建立敘述子。以特徵點為中心，特徵點方位為方向，建立一 個邊長為20 𝜎之正方形區域。再將此方形區域分割成 4×4 個方形子區域(每邊長為 5 𝜎)， 在每個子區域中，以邊長為2 𝜎之哈爾小波遮罩進行運算，也就是說，以左上角第一個 像素點起，每隔𝜎就計算一次，因此在一個子區域共計算了 25 次。最後每個子區域內之 𝑥、𝑦方向的哈爾小波遮罩總和Σ𝑑𝑥、Σ𝑑𝑦與其向量長度總和 Σ|𝑑𝑥|、, Σ|𝑑𝑦|共四個量值， 又共有16 個子區域，所以總計向量維度為 4×4×4 = 64，如圖 3-20]所示，式(3.19)為特 徵點的特徵向量。所有的特徵點建立敘述子後便能使用描述子來達成特徵點的配對。
圖3-20 敘述子產生示意圖
𝑉𝑓𝑖 = (Σ𝑑𝑥, Σ𝑑𝑦, Σ|𝑑𝑥|, Σ|𝑑𝑦|) (3.19)
3.3 離群資料剔除與修正
在增強特徵集合中，我們使用隨機抽樣一致（random consensus sampling, RANSAC） [37]技術進行離群值排除，進行最佳化旋轉與位移矩陣。RANSAC 是一種迭代且非正確 性的演算法，在一組資料中分成內群與離群，推算出一個數學模型，內群經由模型參數 敘述所分佈的資料，離群則代表不適合模型化資料。圖 3-22 中，左邊黑色點為未分群 之資料，右邊為已分為內群與離群，橘色線代表符合內群資料之模型。關於 RANSAC 大致流程如下圖：
圖3-21 RANSAC 架構圖 圖3-22 RANSAC 演算法概念圖 RANSAC 首先以隨機的方式從初步匹配抽取樣本產生符合之模型，計算匹配樣本 與模型之間的差異，繼續重複上述兩步驟，直到迭代結束或得到比預期更小之差異值， 選擇擁有最小差異的模型作為本次的標準，僅留下較好的特徵匹配組合。設定一閥值後 便能分出錯誤或脫離常規之離群值，留下穩定強健的匹配。用於特徵匹配 RANSAC 中 的模型，是一基礎矩陣(Fundamental Matrix)，透過基礎矩陣，匹配關係的兩個特徵點位 置會有以下關係：
𝑝̇𝐹𝑞̇ < 𝜀𝐹 (3.26) 點𝑝̇與𝑞̇為互相匹配的兩特徵點座標，𝐹為基礎矩陣，𝜀𝐹則是一個誤差容許值，理想 狀況下該值應該趨近於零。基礎矩陣𝐹由選定的五個匹配組合基於五點算法[38]進行估 計而來。並使用該矩陣進行評估所有匹配組合之 Sampson 距離，以 Sampson 距離作為 篩選的閥值，通常此閥值設定為0.5。 RANSAC 演算法只考慮非離群之資料以計算適合的模型而不受離群資料的影響， 但無法保證所獲取的內群資料結果最好，故須謹慎定義內群與離群參數。
3.4 特徵匹配
所擷取得到的影像序列經由SURF 特徵點偵測後，可得到特徵點序列，下個步驟要 透過特徵點推估影像的變化，首先進行相鄰序列間的特徵點匹配。 在時間點𝑡，定義影像上被追蹤到的特徵點集合為𝑚𝑖,𝑡 = (𝑢, 𝑣, 𝑑)𝑇，𝑖 = 1,2 … n；在時間 點𝑡 − 1，亦定義影像上被追蹤到的特徵點集合為𝑚𝑖,𝑡−1 = (𝑢′, 𝑣′, 𝑑′)𝑇，𝑖 = 1,2 … n，(u, v) 為影像被追蹤到的特徵點座標，𝑑為其視差值，一次完整視覺測程目標就是要找出各影 像序列中被追蹤之特徵點的旋轉、位移關係： 𝑔̇(𝑚𝑖,𝑡) = 𝑹𝑔̇(𝑚𝑖,𝑡−1) + 𝐭 (3.21) 其中，𝑹為旋轉矩陣，𝐭為位移向量，𝑔̇( )為在歐基里德平面上使用三角測量計算特 徵點三維座標，接著將三維空間上追蹤到的特徵點集合投影至影像上，如下式： 𝑚𝑖,𝑡 = 𝑟(𝑚𝑖,𝑡−1) (3.22)其中，r：ℛ3 → ℛ2，令ℎ() = 𝑔̇−1()，式(3.22)又可表示成： 𝑟(𝑚𝑖,𝑡−1) = ℎ(𝑹𝑔̇(𝑚𝑖,𝑡−1) + 𝐭) (3.23) 接著將三維空間中追蹤到的特徵點集合投影至影像上並與影像上特徵點集合做相 減，可將式(3.22)和(3.23)改為(3.24)： 𝑒(𝑚𝑖,𝑡−1, 𝑚𝑖,𝑡) = 𝑤𝑖,𝑗‖𝑟(𝑚𝑖,𝑡−1) − 𝑚𝑖,𝑡‖ (3.24) 其中權重值𝑤𝑖,𝑗代表在時間點𝑖時特徵點𝑗是否存在，如果存在，𝑤𝑖,𝑗 = 1；若不存在， 𝑤𝑖,𝑗= 0。為使重投影誤差平方最小化，因此可將式(3.24)改為(3.25)： min ∑ 𝑒(𝑚𝑖,𝑡−1, 𝑚𝑖,𝑡)2 𝑛 𝑖=1 (3.25) 利用式(3.25)並使用高斯-牛頓演算法最佳化[29]後可得到最小的誤差結果。
3.5 移動估測
使用相機進行移動估測，根據點的資訊分成三種方式，二維影像特徵點於不同時間 相互對應關係(2D-2D)、三維結構與影像特徵點對應關係(3D-2D)以及三維資料的對齊 (3D-3D)，本研究採用 2D-2D 的方式。 2D-2D 給予𝑡 − 1之三維空間點座標𝑝 = (𝑥𝑖, 𝑦𝑖, 𝑧𝑖)，與𝑡時二維影像中的投影位置𝑞 = (𝑢𝑖, 𝑣𝑖)，在相機參數已知的情況下，該如何求出彼此的轉換方式(𝐑𝒄 t𝑐)與對應關係， 這樣的問題被稱之為 N 點透視問題(Perspective-N-Point Problems, PnP)[19]，如圖 3-23 所 示。圖3-23 N 點透視問題 PnP 示意圖[19] 第一步，在空間中選擇任意4 個能包含所有(𝑥𝑖, 𝑦𝑖, 𝑧𝑖)的參考點{𝑐1, 𝑐2, 𝑐3, 𝑐4, }，三維 空間點可藉由中心坐標𝛼𝑖改寫為： (𝑥𝑖, 𝑦𝑖, 𝑧𝑖)𝑇 = ∑ 𝛼 𝑖𝑗𝑐𝑗 4 𝑗=1 = Γ3×4𝛼𝑖 (3.27) 其中∑4𝑗=1𝛼𝑖𝑗 = 1，Γ = [𝑐1𝑐2𝑐3𝑐4]，且重心座標不會因相機位置而改變。 {𝑐1′, 𝑐2′, 𝑐3′, 𝑐4′}表示為相機在空間中參考點的座標，即： ( 𝑢𝑖 𝑣𝑖 1 ) ~ ( 𝑓𝑢 0 𝑢𝑐 0 𝑓𝑣 𝑣𝑐 0 0 1 ) ( Γ1′ Γ2′ Γ3′ ) 𝛼𝑖 (3.28) 在(3.28)中Γ′ = [𝑐1′ 𝑐 2′ 𝑐3′ 𝑐4′]，Γk′表示第k行之參考點。 將結果轉換為非齊次座標(如：將第一與第二行同除以第三行)： Γ1′(𝑓 𝑢𝛼𝑖) + Γ3′(𝑢𝑐− 𝑢𝑖)𝛼𝑖 = 0 (3.29) Γ2′(𝑓𝑣𝛼𝑖) + Γ3′(𝑣 𝑐− 𝑣𝑖)𝛼𝑖 = 0 (3.30)
觀察𝑁 ≥ 6個點與關係後，如(3.29)與(3.30)般，可建構出齊次線性齊次線性方程： [ 𝑓𝑢𝛼1𝑇 0 (𝑢𝑐− 𝑢0)𝛼1𝑇 0 𝑓𝑣𝛼1𝑇 (𝑣𝑐 − 𝑣0)𝛼1𝑇 ⋮ ⋮ ⋮ 𝑓𝑢𝛼𝑁𝑇 0 (𝑢𝑐 − 𝑢0)𝛼𝑁𝑇 0 𝑓𝑣𝛼𝑁𝑇 (𝑣 − 𝑣0)𝛼𝑁𝑇 ] ( Γ1′𝑇 Γ2′𝑇 Γ3′𝑇 ) = 0 (3.31) 亦可表示為𝐴𝑥 = 0，其中未知數𝑥 = [Γ1′ Γ 2′ Γ3′]𝑇而𝐴為(3.31)左側矩陣。該線性方程
可透過𝐴的奇異值分解(Singular Value Decomposition, SVD)，或者以更有效的方式計算 使𝐴𝑇𝐴𝜉 = 0的向量𝜉~[Γ 1′ Γ2′ Γ3′]𝑇。由於𝐴是一𝑁 × 12之矩陣，當𝑁非常大時，計算成本也 非常高，故改以計算𝐴𝑇𝐴使大小降為12 × 12進行有效率的運算。 要求出𝑀(𝑅𝑐 𝑡𝑐)，須先還原尺度係數。該係數可透過解開𝜉中的[𝑐1′ 𝑐2′ 𝑐3′ 𝑐4′]並與參 考點{𝑐1, 𝑐2, 𝑐3, 𝑐4, }比較彼此距離獲得。𝑀(𝑅𝑐 𝑡𝑐)可由Horn 所提之封閉式[30]以關係𝑐𝑗 → 𝑐𝑗′, ∀𝑗 ∈ {1,2,3,4}求出。最後再以高斯-牛頓最佳化演算法(Gauss–Newton Algorithm)進行 最佳化調整。 α ∑ 𝑒(𝑚𝑖,𝑡−1, 𝑚𝑖,𝑡)2+ 𝑛 𝑖=1 β ∑ 𝑎𝑖𝑒(𝑚̅̅̅̅̅̅̅̅, 𝑚𝑖,𝑡−1 𝑖,𝑡)2 𝑛 𝑖=1 (3.32) (3.32)式用來計算最小重投影誤差，同時也在最佳化旋轉矩陣𝑅與位移向量𝐭，利用(𝑹, 𝐭) 可知道兩時序影像間特徵點的移動關係，也就是計算移動估測。[34][35][36]
3.6 最佳化
由於本體移動估測只針對時間點𝑡 − 1與𝑡的三維及二維對應座標進行運算，故誤差 容易因距離而造成累積，產生路徑估測偏移，此時就需做最佳化調整，使重投影誤差 (Reprojection Error)最小化，如圖 3-24，最常使用的演算法就是光束調整法(Bundle Adjustment, BA)[31]。 圖3-24 重投影誤差示意圖 其中，𝑝𝑟̇ 是特徵點在影像中的二維座標，𝑝𝑟為影像平面上特徵點根據相機內外部參 數投影至三維空間之座標，𝑝𝑟́ 是特徵點三維座標根據相機內外部參數投影回影像平面之 座標，會與𝑝𝑟̇ 形成些許差距，橘色線段即為重投影誤差。 根據投影模型，可定義出用於BA 之重投影誤差𝜙𝐵𝐴𝑅𝑃𝐸如下式：
𝜙𝐵𝐴𝑅𝑃𝐸(𝑝̇𝑠𝑒𝑡; 𝐸1, 𝐸2, … , 𝐸𝑡, 𝑝ℎ𝑠𝑒𝑡) = ∑ ∑ 𝑣𝑖𝑗(𝐾𝐸𝑖𝑝ℎ𝑖𝑗− 𝑝̇𝑖𝑗)2) 𝑁𝑇 𝑗=1 𝑁𝑡 𝑖=1 (3.33) 其中，𝑝̇𝑠𝑒𝑡指追蹤到的特徵點在影像中之位置集合，𝑁𝑡為時間長度，𝐸𝑡指時間𝑡之移 動估測，𝑝ℎ𝑠𝑒𝑡則指追蹤之特徵點在空間中估測之位置集合，𝑣𝑖𝑗檢查時間𝑖的特徵點𝑗是否 存在，若存在𝑣𝑖𝑗 = 1，若否𝑣𝑖𝑗 = 0，而𝐾𝐸𝑖𝑝ℎ𝑖𝑗將𝑝ℎ𝑖𝑗根據時間𝑖之移動估測和相機參數𝐾 將影像投影至影像平面，與𝑝̇𝑖𝑗計算差平方距離。 BA 之作用在於調整移動參數與估測之三維特徵點座標，使重投影誤差𝜙𝐵𝐴𝑅𝑃𝐸最小 化，要達成這個目標，可利用萊文貝格-馬夸特演算法(Levenberg-Marquardt L-M)[32]。 L-M 演算法為一種迭代演算法，在進行光束調整後的重投影誤差𝜙𝐵𝐴𝑅𝑃𝐸的最佳化運 算，能夠解決非線性問題找出數值最佳解， L-M演算法結合梯度下降法(Gradient Descent algorithm)與高斯-牛頓演算法(Gauss-Newton algorithm)優點：梯度下降法收斂速度較緩慢， 很容易陷入區域解而失去找到最佳解的機會；高斯-牛頓演算法則有較佳的收斂速度，但也 因收斂太快，錯過最佳解。因此L-M 演算法在初期離最佳解較遠時使用收斂速度較快的 高斯牛頓法，在接近最佳解時則採用梯度坡降法，降低收斂速度。L-M 演算法包含一係 數λ，λ越大，L-M 演算法的特性越接近梯度坡降法；反之，λ越小，L-M 演算法的特性
越接近高斯牛頓法。通常用於處理最小平方曲線擬合(Least Squares Curve Fitting)問題，如
(3.34)式： 𝑠(𝛽) = ∑[𝑦𝑖− 𝑓(𝑥𝑖, 𝛽)] 𝑘 𝑖=1 (3.34) 其中𝑓代表曲線，有𝑘組資料，𝑥代表輸入資料，𝑦代表𝑥的輸出資料，而 L-M 演算法 可將參數𝛽最佳化，使𝑠得到最佳化後之最小值。 S. W. Huang 在 2012 年提出的雷射掃描儀內部參數最佳[33]，也是使用 L-M 演算 法，將(3.27)式套入(3.34)式可推得最佳化(3.35)式：
𝜙(𝛽) = 𝑚𝑖𝑛 ∑ 𝛿(𝑛⃑⃑⃑⃑ , 𝑝𝛽 𝑖,𝛽)2 𝑘 𝑖=1 (3.35) 𝛽為最佳化參數，用來求得(3.35)式之最小值。本節最佳化外部參數𝑅與𝑇之過程與 S. W. Huang 所討論的方法相同，因此可推得相機之外部參數最佳化目標成為(3.36)式： 𝜙𝐸(𝛽) = 𝑚𝑖𝑛 ∑ 𝛿[𝜔(𝐸𝛽, 𝑛⃑⃑⃑⃑ ), 𝑝𝑐 𝑖,𝛽]2 𝑘 𝑖=1 (3.36) L-M 演算法會逐步調整偏移量，調整一個能使𝜙𝐸得到最小的最佳化參數𝛽，再將𝛽 轉換回最佳化的外部參數𝐸𝑜𝑝𝑡 = (𝑅𝑜𝑝𝑡, 𝑇𝑜𝑝𝑡)。最後將留下來的點進行二次調整後，便可 以得出最佳化之二維座標點雲以及移動估測。
第四章 實驗結果與討論
本章節使用第三章相機視覺系統所提出的視覺測程方法進行實驗與分析。本實驗使 用無人機搭載相機錄製影像，如圖 4-1，；本實驗將操控無人機進行可移動式的資料擷 取，並使用 GPS 資料比對，以達到移動估測之目的。實驗所錄製之影像已使用 3.1.2 節 的方法進行校正。 圖 4-1 本論文實驗無人機4.1 實驗資料蒐集
本實驗分別對室內、室外街道環境與校園街道進行五組實驗資料採集， (1) 實驗室 室內 相機鏡頭朝前 (2) 戶外街道 直線 相機鏡頭朝前 (3) 戶外街道 環繞 相機鏡頭朝前 (4) 戶外街道 直線 相機鏡頭朝下 (5) 校園街道 曲線 相機鏡頭朝前 如圖 4-2。室內環境為本實驗室(a)，空間較狹隘且複雜；街景為楠梓區大學二十六 街、大學二十三街、大學二十二街與大學西路(b)(c)(d)，環境空間大且開放，具有較多 的影像特徵，但易受其他環境因素影響資料。校園街道為高雄大學內環道(e)。衛星影像 如圖 4-3，取自 Google Map[39]。 (a)實驗室 室內環境(b)戶外街道環境 大學二十六街 相機鏡頭朝前
(c)戶外街道環境 大學二十二街
(d)戶外街道環境 大學二十六街 相機鏡頭朝下
(e) 校園街道 高雄大學內環道
(a)戶外街道
(b)校園街道
