國立高雄大學土木與環境工程學系
碩士論文
數位影像相關係數法應用於影像接合之研究
Application of the digital image correlation method to image
stitching
研究生：周祐諒撰
指導教授：童士恒
謝誌
在研究所的求學過程中，遇到重重的難關，也遭遇到挫折，但有壓力， 才有動力，雖然日子過得辛苦，但這段求學的過程，卻是令人難以忘懷的。 能夠順利完成學業並取得學位，最要感謝的是我的指導教授童士恒教 授，如果沒有他的協助與指導，我無法順利的完成這篇論文。在研究的方 法、研究的方向及解決問題的方向上，童士恒教授給我很多的建議，並循 循善誘，讓我逐步的完成這篇論文，在此要向童士恒教授致上最大的謝意。 其次要感謝研究所課程任課教授們，由於他們努力的教學及對課業的 要求，讓我在土木與環工的領域上，學到更多，增加了更多的專業知識。 再者要感謝口試委員，蕭漢威教授及施明祥教授，對我的論文架構及 內容，不吝惜的指導及提供許多寶貴的意見，讓我的論文得以順利的完成。 於研究的過程中，也適時提供幫助，讓我能解決問題，順利的完成研究。 另外也要感謝研究室的夥伴們以及研究所的同學，有他們的幫忙我才 能更有效率的處理問題，遇到困難時互相扶持，互相勉勵，大家互相的幫 助與打氣，提供了堅持下去的動力來源。 最後要感謝我的家人，在背後默默的支持，容忍我的任性，謝謝您們！I
目錄
目錄 ... I 圖目錄 ... IV 表目錄 ... VIII 中文摘要 ... 1 英文摘要 ... 3 第一章 前言 ... 6 1.1 研究動機 ... 6 1.2 研究目的 ... 7 第二章 文獻回顧 ... 9 2.1 數位影像相關係數(DIC) ... 9 2.2 特徵點提取 ... 10 2.3 影像匹配 ... 13 2.4 影像拼接與融合 ... 14 2.5 文獻總結 ... 15 第三章 特徵點提取與影像匹配相關理論 ... 17II 3.1 特徵點提取法 ... 17 3.1.1 特徵點擷取法文獻比較 ... 18 3.1.2 特徵點提取法抵抗影像旋轉能力評估 ... 21 3.1.3 特徵點提取法運算速度比較 ... 24 3.2 SIFT 演算法 ... 28 3.3 SURF 演算法 ... 34 3.4 數位影像相關係數法 ... 39 第四章 空中三角測量 ... 42 4.1 相機參數 ... 42 4.2 空中三角測量原理 ... 45 4.3 PIX4DMAPPER軟體 ... 47 第五章 實驗驗證與討論 ... 49 5.1 影像匹配演算法精度比較實驗 ... 49 5.1.1 試體 ... 49 5.1.2 實驗步驟 ... 50 5.1.3 實驗結果與討論 ... 51 5.2 三維階梯試體實驗 ... 59
III 5.2.1 試體 ... 59 5.2.2 實驗儀器與配置 ... 60 5.2.3 實驗步驟 ... 63 5.2.4 人工選點與 DIC 定位比較 ... 68 5.2.5 連結點影響之結果與討論 ... 69 5.3 三維模型實驗 ... 73 5.3.1 試體 ... 74 5.3.2 實驗儀器與配置 ... 74 5.3.3 實驗步驟 ... 77 5.3.4 人工選點與 DIC 定位比較 ... 80 5.3.5 連結點影響之結果與討論 ... 81 第六章 結論與建議 ... 85 6.1 結論 ... 85 6.2 建議 ... 86 參考文獻 ... 88
IV
圖目錄
圖 3.1 影像拼接流程 ... 18 圖 3.2 旋轉影像試體，影像大小為 1694×1940PIXEL，分別為 0 度(左上)、30 度(右上)、60 度(左下)與 90 度(右下) ... 21 圖 3.3HARRIS 角點偵測法抵抗旋轉能力評估實驗結果，圖中紅色圈位置為 特徵點點位 ... 22 圖 3.4SIFT 演算法抵抗旋轉能力評估實驗結果，圖中紅色圈位置為特徵點 點位 ... 23 圖 3.5SURF 演算法抵抗旋轉能力評估實驗結果，圖中紅色圈位置為特徵點 點位 ... 23 圖 3.6BRISK 演算法抵抗旋轉能力評估實驗結果，圖中紅色圈位置為特徵 點點位 ... 24 圖 3.7 實驗用幾何圖形影像，影像尺寸為 1694×492 PIXELS ... 25 圖 3.8 隨機拍攝影像，改變大小進行演算速度比較，原影像尺寸為 5184× 3456 PIXELS ... 26 圖 3.9 將影像依尺度倍率以 2 倍為一單位分組，並縮小影像為原來的一半， 每一組別中進行高斯平滑，再將同一組相鄰的經高斯平滑處理後的影像 兩兩相減可得其 DOG 影像(LOWE,2004) ... 29V 圖 3.10 中間層的檢測點和它同尺度的 8 個相鄰點以及上下相鄰尺度的 9 × 2 個點，共 26 個點進行比較，以確保在尺度空間和二維圖像空間都檢測 到極值點(LOWE,2004) ... 30 圖 3.11 特徵點描述示意圖(LOWE,2004) ... 33 圖 3.12 積分影像區塊值（BAY ET AL.,2006） ... 35 圖 3.13 高斯二階微分模板及簡化（BAY ET AL.,2006） ... 36 圖 3.14SURF 與 SIFT 提取法濾波器比較圖（BAY ET AL.,2006） ... 37 圖 3.15 影像中的分割區域（BAY ET AL.,2006） ... 38 圖 3.16 特徵描述表示（BAY ET AL.,2006） ... 39 圖 3.17 物體表面之結構性斑紋與次級影像（網格）示意圖（施等,2006） ... 40 圖 3.18 物體表面上之次級影像變形前後相對位置示意圖（施等,2006） . 41 圖 4.1 偏心扭曲之最大切線扭曲軸及最小切線扭曲軸（WENG,1992） ... 44 圖 4.2 空中三角原理示意圖 ... 46 圖 5.1 原始影像，影像大小為 2918×2189 畫素 ... 50 圖 5.2 五種方法於不同X方向位移下之誤差平均值比較圖 ... 57 圖 5.3 五種方法於不同X方向位移下之誤差標準差比較圖 ... 57 圖 5.4 五種方法於不同Y方向位移下之誤差平均值比較圖 ... 58
VI
圖 5.5 五種方法於不同Y方向位移下之誤差標準差比較圖 ... 58
圖 5.6 三維階梯試體，每階層高程差為 5MM ... 60
圖 5.7CANON EF-S60MM F/2.8MACRO USM 與 CANON EOS650D。 ... 61
圖 5.8 大型 XY 移動平台 ... 61 圖 5.9 小型 XY 移動平台 ... 62 圖 5.10 三維階梯試體航空模擬拍攝實驗配置示意圖 ... 62 圖 5.11 三維階梯試體試體控制點與檢核點三維座標量測 ... 64 圖 5.12 所有控制點之示意圖 ... 66 圖 5.13 將控制點標於控制點影像上 ... 69 圖 5.14 檢核點高程誤差平均值比較圖 ... 72 圖 5.15 檢核點高程誤差標準差比較圖 ... 72 圖 5.16 三維模型試體... 74
圖 5.17 CANON EF-S60MM F/2.8MACRO USM 與 CANON EOS650D ... 75
圖 5.18 高程側微計，下降至懸臂尾端針頭碰觸欲量測位置後，即可讀取高 程數據 ... 76
圖 5.19 大型 XY 移動平台 ... 76
圖 5.20 中型 XY 移動平台 ... 77
VII
圖 5.22 現地模型試體控制點與檢核點三維座標量測 ... 79 圖 5.23 檢核點高程誤差平均值比較圖 ... 84 圖 5.24 檢核點高程誤差標準差比較圖 ... 84
VIII
表目錄
表 3.1 特徵點擷取方法比較表 ... 20 表 3.2HARRIS、SIFT、SURF 及 BRISK 特徵點提取法運算時間 ... 25 表 3.3 不同影像大小，SIFT、SURF 及 BRISK 運算時間 ... 27 表 5.1 五種匹配方法耗時比較表 ... 53 表 5.2 五種方法於不同X方向位移下之誤差平均值 ... 53 表 5.3 五種方法於不同X方向位移下之誤差標準差 ... 54 表 5.4 五種方法於不同Y方向位移下之誤差平均值 ... 55 表 5.5 五種方法於不同Y方向位移下之誤差標準差 ... 56 表 5.6 三維階梯試體三個測站所拍攝之影像 ... 64 表 5.7 人工選點與 DIC 定位，誤差平均值與標準差 ... 69 表 5.8 匯入不同數量之利用 DIC 匹配的特徵點所得之檢核點高程誤差值與 標準差 ... 71 表 5.9 匯入不同數量利用 SURF 匹配之特徵點的檢核點高程誤差值與標準 差 ... 71 表 5.10 人工選點與 DIC 定位，誤差平均值與標準差 ... 81 表 5.11 不同特徵點點數利用 DIC 匹配匯入軟體，檢核點高程誤差值與標準IX
差 ... 83
表 5.12 不同特徵點點數利用 SURF 匹配匯入軟體，檢核點高程誤差值與標
1
數位影像相關係數法應用於影像接合之研究
指導教授：童士恒 博士 學生：周祐諒 國立高雄大學土木與環境工程學系 摘要 無人飛行載具(UAV)可於空中進行航空攝影，近幾年於各領域發展快速，本研究利 用數位影像相關係數法(DIC)來提升航拍影像拼接的精準度與三維數值模型的精度。 DIC 在本研究中的應用有兩個方面，第一在於應用 DIC 定出控制點在不同影像上的位 置，藉由提高控制點的精度來提高影像拼接與三維數值模型的精度；第二在於利用 SURF 特徵點提取法提出影像的特徵點，再利用 DIC 比對出影像間相對應的特徵點影像 座標，將這些比對出點做為連結點匯入 Pix4Dmapper 軟體中增加影像拼接與三維模型製 作的精度。 本研究先進行了影像匹配精度實驗，將影像於 x 和 y 方向移動不同移動量，對移動 前與移動後之影像分別進行了 SIFT、SURF、BRISK 與 DIC 匹配，結果顯示 SIFT、SURF 與 BRISK 在 x 與 y 方向的精度皆在 1 到 0.1 個 pixel，而 DIC 則是 0.02 個 pixel，顯示 DIC 的精度確實高於一些常用的匹配方法。另外，研究中比較了以人工選定控制點與人工選點再輔以 DIC 配對的差異，三維 階梯實驗結果顯示，利用 DIC 分析的控制點影像座標，做出來的模型檢核點高程誤差
2 標準差為 0.08mm，只利用人工選點的標準差為 0.201mm，三維模型實驗結果顯示，利 用人工搭配 DIC 分析的控制點影像座標，做出來的模型檢核點高程誤差標準差為 0.48mm，只利用人工選點的標準差為 0.696mm，階梯模型與三維模型結果皆顯示在單 純將 DIC 應用於控制點定位上已能改善精度。 最後本研究進行了兩個室內實驗來評估不同數量的 DIC 連結點與 SURF 連結點對 影像拼接與三維數值模型精度之影響，從階梯模型與三維模型的檢核點誤差得知，加入 DIC 連結點有助於提升拼接影像與三維數值模型的精度，且改善幅度大致隨點數增加而 提高，且效果比 SURF 比對出的連結點更好。具高精度高程之階梯試體之試驗結果顯示， 未匯入連結點的檢核點高程誤差標準差為 0.08mm，加入 DIC 連結點的最佳結果為 0.049mm，而由三維模型之試驗結果可得未匯入連結點的高程誤差標準差為 0.48mm， 匯入 DIC 連結點的最佳結果為 0.27mm。由此顯示，應用 DIC 確實能有效改善影像拼接 與三維數值模型建立之精度。
3
Application of the digital image correlation method to
image stitching
Advisor: Dr. Shih-Heng Tung Student: Yu-Liang Chou
Department of Civil and Environmental Engineering National University of Kaohsiung
ABSTRACT
Unmanned aerial vehicle (UAV) can be used to take aerial photographs, and it has been rapidly developed in many fields in recent years. This study uses digital image correlation method (DIC) to increase precision of image stitching and three-dimensional terrain model. DIC will be applied in two ways. Firstly, DIC is used to locate control points’ image positions in different photos, so that the precision of image stitching and three-dimensional terrain model will increase by raising control points’ precision. Secondly, DIC is used to match feature points which extracted by SURF algorithm and then these points are imported into Pix4dmapper to increase the precision of image stitching and three-dimensional terrain model.
An experiment is carried out to compare the precision of various image matching methods, such as SIFT, SURF, BRISK as well as DIC. The results of different movements in x and y directions show that the precisions of SIFT、SURF and BRISK range from 0.1 to 1
4
pixel and the precision of DIC is 0.02pixel. It shows that DIC is more accurate than the other methods.
Furthermore, this study compares two different methods adopted to obtain the image coordinates of control points. In one of the method, the control points’ image coordinates are positioned manually. The other method will use DIC method to improve the precision of manually positioned control points. Ladder specimen’s results show that the positioning precision by using DIC is 0.08 mm while the manual positioning precision is 0.201 mm, 3D model specimen’s results show that the positioning precision by using DIC is 0.48 mm while the manual positioning precision is 0.696 mm. It shows that applying DIC to positioning control points can improve the precision.
Eventually, this study carries out two indoor experiments to evaluate the influence of different amounts of DIC connection points and SURF connection points on the precision of image stitching and three-dimensional terrain model. According to the altitude error of check points from ladder specimen and 3D model specimen, we can find that importing DIC connection points can increase the precision of image stitching and three-dimensional terrain model. The magnitude of improvement will be raised as the amount of connection points increases, and its effect is better than importing SURF connection points. Ladder specimen’s result shows that the altitude error is 0.08mm without connection points. The best result is 0.049mm while importing DIC connection points. The 3D model specimen’s result shows that the altitude error is 0.48mm without connection point. The best result is 0.27mm while importing DIC connection points. These results show that applying DIC to image stitching
5
6
第一章 前言
1.1
研究動機
臺灣位於亞熱帶地區且位在板塊活動帶上，常年受到颱風與地震 的威脅，而強震與伴隨大量雨水的颱風都會降低山坡地的穩定性，造 成邊坡崩塌、滑動等與土石流等問題。如 2001 年的大興村土石流掩埋 事件，就是因為桃芝颱風帶來的豪雨，再加上之前 921 大地震後造成 土石鬆動及無數的崩塌地，進而引發大規模之土石流災害，造成約 150 萬立方公尺之土石下移，大興全村 184 戶，近 150 戶遭土石掩埋的慘 重災情。另外，2009 年相當嚴重的高雄縣甲仙鄉小林村滅村事件，亦 是因莫拉克颱風大量降雨造成土石鬆動，使獻肚山發生嚴重山崩，部 分土石直接掩埋小林村靠山地區，另有一部分則形成堰塞湖，當堰塞 湖潰堤時直接將小林村覆滅。 災害發生後，常需取得災區的即時資訊以便評估災害規模及規劃 救災行動，唯災害發生常伴隨道路與通訊的中斷，使得災害資訊的獲 得有相當的困難。為克服此困難，衛星影像與航空影像亦逐漸被用來 取得災區的資訊，衛星影像的取得不受天候影響，但因為衛星的高度 高，影像的解析度不佳，且在天候不良的情況下，地表狀況容易被雲 層所遮蔽。相較於衛星影像，傳統航照因飛行安全考量，對天候的要 求較高，但是可取得解析度較高的影像。 近年隨著各領域科技的進步，無人飛行載具（Unmanned Aerial Vehicle，簡稱 UAV）於各領域的發展應用皆相當成熟，常搭載影像擷7 取設備，來獲取空拍影像，相較於傳統的衛星影像及空拍影像，具有 解析度高、取得成本低的優點；且因無人員在飛機上，在緊急必要時， 仍可在天候條件較差的情況下執行任務取得影像。另於災變發生時， 無人飛行載具也可進入較難到達的區域進行拍攝，獲取該區域的空拍 圖，經過影像處理後，得知相關空間資訊，進而提供即時救災、坡地 監測與崩塌量計算等資料。 無人飛行載具在飛行高度方面相對較低，使單張影像所涵蓋範圍 較少，因此需要大範圍區域航照圖時，須拍攝多張影像，並透過影像 接合技術產出拼接影像，供後續應用。回顧過去的影像拼接方式，在 影像匹配部分，大多利用特徵匹配法，如 SIFT、SURF 與 BRISK 等。 然而這些方法精度約為 pixel 或 sub-pixel 等級，因此要將拼接影像利用 於邊坡監測或是崩塌量計算等，便會造成精度上相當的誤差。而數位 影像相關係數法(DIC)進行影像匹配時，可達約 0.01pixel 的精度，因此， 本研究擬將數位影像相關係數法(DIC)應用於影像匹配上，探討此方法 對於航照拼接影像與三圍數值模型精度的改善效果。
1.2
研究目的
現今航照圖拼接操作多採用套裝軟體，如 Photoscan、Pix4Dmapper 與 Acute3D 等，大多是利用特徵點擷取與特徵匹配的方式來進行影像 拼接，並利用影像控制點與外方位資訊為輔助。然而特徵匹配方式所 找出的匹配點座標精度約為 1pixel，所以當控制點數量較少或外方位資 訊不足時，拼接影像的精度與接合準確度會有下降的趨勢，進而影響
8 後續三維模型與崩塌量計算等精度，而數位影像相關係數法能提高影 像拼接精度，改善此問題。 為了探討數位影像相關係數法對影像拼接精度改善情況，本研究 利用 python 與 openCV2 編寫影像特徵點擷取與特徵匹配程式，再搭配 施明祥教授與童士恒教授開發出的數位影像相關係數法程式，將其與 常用的特徵匹配法如 SIFT、SURF 與 BRISK 進行精度比較，觀察利用 數位影像相關係數法是否提升影像匹配精度。且進一步將演算出的匹 配點匯入套裝軟體 Pix4Dmapper 進行影像拼接與三維數值模型製作， 比較匹配點匯入前後拼接影像接合準確度與三維數值精度的差異，探 討改善的成果。
9
第二章 文獻回顧
2.1
數位影像相關係數(DIC)
數位影像相關係數法（Digital Image Correlation method，DIC）是 以“找尋演算法”為基礎，使用數位影像相關係數，比對兩張影像的 局部相關性，判定變形前後對應關係。利用有限元素法的概念將變形 前後之試體分割成小網格，藉比較變形前後網格內灰階值之相關性， 在變形後影像中尋找變形前網格的相對位置。經由變形前後網格的相 對位置，即可計算出各網格之應變量。以下為相關研究： (1) Peters et al（1982）提出利用計算機分析物件受力前後之數位影像， 求得表面位移量與應變的方法。 (2) Chu et al（1985）提出以結合變形理論及數位影像相關係數法之量 測技術，及以內插理論來擴展此技術之應用範圍，並以實際案例來 驗證此技術之實用性。 (3) Lu et al（2000）對於以影像關係法進行應變分析時，使用高次內插 函數模擬位移場，考慮高次位移項之梯度的影響進行探討，實驗結 果顯示在小變形時，忽略二次位移項之梯度影響不大，但忽略二次 位移項之梯度的誤差會隨著變形的增加而愈來愈大。 (4) Vellinga et al（2000）所發表的研究中提到，由於數位影像相關係 數法之解析度只與像素有關，與實際之長度大小無關，所以結合掃 瞄式電子顯微鏡與數位影像相關係數法，成功的觀測到微小範圍內 的應變，換算成實際空間中的長度，其量測到的位移大約是微米等
10 級，並認為有可能達到次微米的解析度的發展空間。 (5) Lopez et al（2004）提出以兩台高解析度CCD鏡頭試驗性的進行三 維的DIC量測，此法成功判別出複合材料中各材料行為的差異。 (6) Tung et al（2005）利用DIC量測剛體位移，其x方向與y方向的位移 誤差值為0.004個像素。 (7) Lecompte et al（2006）的研究，提出了不同大小的分析網格和不同 尺寸的斑點記號之間的關係對DIC分析精度的影響，並建議最佳化 網格大小與記號大小，以提高DIC分析的精度。 (8) Shih et al（2010）的研究中，利用一台攝影機進行動態的量測，結 果顯示一懸臂樑之模態可利用DIC量測出，並可以此來判斷是否有 缺陷。 (9) Tung et al（2011）的研究，利用單一一台相機進行三維DIC量測， 可避免因兩臺相機不同的機械與光學性質造成的誤差，經扭曲校正 後影像，物距精度可達0.0043%。
2.2
特徵點提取
影像拼接的第一步為尋找影像中具有代表性與獨特性的特徵點， 主要利用數學演算法來進行影像特徵點的提取，數學演算法中具代表 性的有 Harris 角點偵測演算法、SIFT 演算法、SURF 演算法與 BRISK 演算法。此外，亦有不少針對上述方法所提出的改善方法，有些研究 為了提升上述方法的運算效率，只對粗估重疊區域進行特徵點提取來11 減少運算時間，有些為了增加演算法穩健性(robustness)，對演算法進 行修正，以下回顧特徵點提取相關方法: (1) Harris (1988) 提出 Harris 角點偵測演算法，主要用來提取影像中的 角點做為特徵點，角點為輪廓上的連接點，是一種改變視角時能保 持穩定的特徵點，與其任意方向鄰近點差異性大且俱獨特性，為不 錯的特徵點。角點的提取方法主要利用搜索小網格在影像中做小距 離滑動，觀察影像強度變化，如果偵測為角點則網格於任一方向搜 索，強度值皆有相當大的變化。此方法的優點為，偵測效率高且具 有較高的重複率，而缺點為，當影像尺度改變時特徵點位也會隨之 改變的情況與運算時間較長的問題。
(2) Lowe (2004) 改良Harris演算法，提出了SIFT(Scale-Invariant Feature Transform)演算法，改善Harris演算法無法抵抗影像尺度改變的問題。 此方法先對影像做高斯差分( Difference of Gussian, DOG )與影像金 字塔(Image pyramid)分層的方式，使影像特徵點俱備抵抗影像尺度 改變的能力，以及對特徵點做向量化描述來抵抗影像旋轉。此方法 優點為，當影像產生旋轉、尺度與視角改變時的有相當的穩定性。 缺點為，運算量較大相對花費時間較多。
(3) Bay (2006) 提出了SURF（Speeded Up Robust Feature）演算法，改 善SIFT演算法的特徵點提取法，具有抵抗影像旋轉與影像尺度改變 能力且運算快速的方法。此方法快速的主要原因為，利用積分影像 搭配Hessian matrix探測器找出特徵點，且藉由改變filter大小而非影
12 像大小來創造尺度空間，並利用X、Y方向的一階Haar wavelet response來描述特徵點。優點為，運算快速，且特徵點重複性高， 抵抗影像視角改變、尺度變化與模糊能力強。缺點為，延伸的描述 子會喪失。 (4) Yang (2011) 提出了一個改良SURF的演算法，先對影像進行仿射變 化(affine transform)校正，改善當視角有大改變時，SURF演算法會 匹配錯誤的缺點，以及重建尺度空間來改進當尺度大於原影像10 倍以上時SURF演算法特徵點數會減少的情況。優點為，改良SURF 法對於影像的視角變化、旋轉及亮度等皆比SURF有良好的不變 性。 (5) Fathimaa (2013) 提出了結合不變矩和SIFT演算法去做影像處理， 先將圖像分成好幾區塊，再計算各區塊的彎矩值，此彎矩值不會受 平移、旋轉與尺度變化影響，將各區塊與目標影像的區塊進行比較， 找出相似重疊的區域，再對此區域進行SIFT特徵提取法，來提升運 算的速度， 此篇論文發現有70%的點不在重疊區域內，所以用此方 法減少83%的特徵點數量，進而減少了35%的時間，且對影像的旋轉、 尺度變化有其不變性，然而會有13%的可用特徵點將會被踢除。 (6) 魏台軍 (2006) 提出了結合重疊區域及Harris演算法的特徵點提取 法，利用影像間重疊區域的快速定位策略，透過快速的定位法則協 助系統取得重疊區域，以加速後續的特徵擷取步驟，以及區域化的 特徵擷取。有別於過往的全域性特徵擷取，只於局部的區域(重疊
13 區域)進行處理，此方法可有效濾除非有效點及減少特徵點配對上 的計算複雜度進而減少運算時間，但此方法因為利用角點偵測法， 所以不能抵抗尺度變化。 (7) Leutenegger (2011) 提出了BRISK演算法，藉由比較區域內的灰階 值找出特徵點，為較省時的特徵點提取與匹配方法，特徵點運算的 速度較SURF和SIFT快，並具有不錯的尺度和視角改變抵抗能力。 其描述子為二進位的向量，匹配時的運算速度相當快。此方法雖然 減少了運算時間，但對於較大的旋轉抵抗能力較差，而且誤差也較 大。
2.3 影像匹配
提取完影像中具代表性的特徵點後，需將不同影像中的特徵點做 相互匹配，找出匹配的點位，影像匹配會利用特徵點本身的梯度值、 距離差或是灰階強度等做相互比對，提供後續拼接使用。此步驟對拼 接影像的精度影響甚鉅，精度高的影像匹配法，製作出的拼接影像較 準確，三維模型的製作與崩塌量計算等精度較高。以下為相關影像匹 配方法。
(1) Lowe (2004) 提出的次近鄰居法(Second-Closest Neighbor)，可用於 如SIFT、SURF與BRISK等特徵點演算法之特徵描述子比對，將影 像中的特徵點置入另一張影像中尋找匹配時，利用描述子的向量， 計算兩描述子的最近距離與次近距離，利用這兩個值評估特徵點匹 配性。
14 (2) 謝銘倫 (2002) 提出彩色碼比對法，將影像梯度方向做分類並量化 產生彩色碼做為特徵角點比對資訊。影像比對工作便可經由彩色碼 的比較來達成，若兩點的彩色碼越相似，則表示此兩個角點非常有 可能是互相對應的兩點。透過距離限制和雙向對應一致性的要求， 以及與相鄰點關係的比較，並計算出一評估對應關係的分數，最佳 的對應點即可被選擇。優點為，準確性高且效率高。缺點為，特徵 點的選取上必須避免擷取到空間中的歪斜線在影像上所形成的交 點，或是不屬於角點的特徵點，否則對應時發生錯誤的機率會較高； 當場景中有些物體會隨時間改變位置，導致參考已對應組合時，因 物體運動方向不一致，使得對應產生錯誤。 (3) Chang (1995）提出利用仿射變化的四個參數尺度s、旋轉角θ和x與y 方向的平移tx、ty，來找出兩影像間的匹配點。此方法可以在影像轉 換、旋轉及尺度不同情況下正常操作。缺點為，花費的時間多，而 且無進行內插計算，匹配出的點會有精度上的問題 。 (4) Rankov (2005) 提出交互關聯匹配法，為一種基於交互關係的匹配 方法，用來找出兩影像間最相似的點。主要利用搜尋視窗在欲偵測 之影像中進行搜索，並利用相關係數來比對出相互匹配的點，此篇 研究中並無探討非整數點的內插，所以精確度相對較低，且無利用 特徵點做匹配，運算量較大，費時較長。
2.4 影像拼接與融合
影像匹配完後，可利用匹配點進行影像間相互關係的運算，利用15
此關係將影像進行拼接，拼接後於影像結合處會有明顯的接縫產生， 融合影像演算可用來處理接合縫問題，以及處理模糊帶。融合的方法 有權重分配法、多重帶混合法等。以下回顧相關的拼接融合方法: (1) 王君如 (2014) 提出影像拼接幾何校正法，先利用除錯方法，將錯
誤的特徵點做初步的刪除，再計算均方根誤差（Root Mean Square Error；RMSE）值的方式，挑選出誤差小的特徵點作為幾何校正轉 換模型的參數，將特徵點代入幾何校正模型將單張影像進行坐標轉 換，將不同影像校正至同一坐標系來拼接影像，最後再進行影像重 新採樣。從此文獻中的幾何校正模型發現，無法處理徑向扭曲問題， 且文獻中誤差值計算方法較為繁複，花費時間較多。 (2) Rankov (2005) 提出以梯度演算法為基礎的互補影像融合法，將影 像接合處的地方做調整，利用影像間重合部分的梯度權重來做分配 混合，主要用來消除拼接影像尖端的強度變化。利用此方法混合出 的影像，容易會有雙影影像模糊情形產生，且對不同光照下影像需 適當處理光線補償問題。 (3) Burt (1983) 提出多重帶混合法，製作影像拉普拉斯影像金字塔， 並對金字塔中各層影像做相互融合，並依照權重分配，來處理影像 間接合帶問題。此方法能確實消除影像接合處的模糊帶，且不會有 雙影現象產生。
2.5 文獻總結
數位影像量測技術是一種新興之技術，由於近年來相機的快速發16 展，使用此方法進行量測之精度日益提高，而且應用範圍相當的廣泛， 例如工程界中經常需要材料的應力應變關係。本研究希望將數位影像 相關係數法應用於航拍圖拼接，改善影像拼接精度與接合度。 先找出重疊區域再進行特徵點擷取來加快運算速度的方式，對於 航拍影像拼接精度與接合度帶來的效益不大。準確的拼接影像與三維 模型需仰賴精度較高的影像匹配技術，然而利用高精度影像匹配技術 提升拼接影像與三維模型精度的研究相對較少。所以本研究希望利用 數位影像相關係數法高精度的優點，增進拼接影像的精度與接合度， 進而提升三維模型精度以利崩塌量計算等處裡。雖然目前有商用軟體 可以進行航拍圖的拼接等操作，但影像匹配方式大多採用特徵匹配方 式，精度約為1pixel，所以需要外方位資訊與控制點來提高準確度，本 研究會將數位影像相關係數法算出的點匯入套裝軟體中，探討數位影 像相關係數法是否有助提升影像拼接接合度與精度。
17
第三章 特徵點提取與影像匹配相關理論
本研究主要目的為利用數位影像相關係數法進行影像特徵點匹配， 比較數位影像相關係數法與較常用的特徵匹配法精度差異，並將數位 影 像 相 關 係 數 法 與 特 徵 匹 配 法 所 演 算 出 的 點 匯 入 套 裝 軟 體 Pix4Dmapper內執行影像拼接與三維模型建立，觀察不同方法所製作出 的拼接影像接合度，與三維地表模型檢核點高程精度。本章節將依序 介紹，特徵點提取法、數位影像相關係數法。3.1 特徵點提取法
影像拼接的流程如圖3.1，分為三個步驟，第一步為影像特徵點提 取，找出影像中較獨特且具代表性的特徵點，如角點或灰階變化較劇 烈的地方等，加快影像匹配的速度。第二步，對第一步搜尋出的影像 特徵點進行匹配，利用特徵點灰階值或特徵向量等，比對相關性高的 匹配點，利用這些匹配點做接圖處理。最後將影像進行接合與融合， 利用第二步找出的匹配點，演算出影像相對位置關係，將影像套合在 一起，再利用影像融合技術將影像間的接縫進行處理，接合後的圖像 便可用來製作三維模型及計算崩塌量等。影像特徵點提取方法常用的 有Harris角點偵測法、SIFT、SURF與BRISK演算法，本研究先對相關 特徵點提取法進行文獻比較，再針對這四種方法進行抵抗影像旋轉能 力評估與運算速度比較，決定本研究所採用的特徵點提取法。18 影像擷取 特徵點提取 特徵點匹配 影像接合與融合 圖 3.1 影像拼接流程 3.1.1 特徵點擷取法文獻比較 以下對第二章前人所提出的特徵點提取法，比較各方法的相關指 標，指標包括對於旋轉、亮度與視角變化的穩定性、演算所需時間、 抵抗影像扭曲性以及特徵點點數。回顧的演算法如下: 1. Harris角點偵測法，找出影像中的角點做為特徵點，但找出的特 徵點無法抵抗尺度變化及旋轉帶來的影響。 2. 尺度不變的特徵點擷取法(SIFT)，利用影像金字塔和向量描述符 改善Harris角點偵測法的缺點，但運算花費時間較多。
19 3. 速度快且穩定的特徵點擷取法(SURF)，利用積分影像和box filter， 提高演算法的運算速度，並維持強健的穩定性。 4. 快速的 BRISK特徵點擷取與匹配方法，特徵點運算的速度較 SURF和SIFT快，但抵抗旋轉能力較差。 5. 利用粗估的重疊區域和直方圖均化加強影像的強度，對此區域進 行特徵點提取，來增加特徵點提取速度，及使用改良的SURF演 算法，處理SURF演算法無法抵抗大視角改變帶來的影響 。 6. 結合不變矩以及尺度不變的特徵點擷取法，利用不變矩粗估重疊 區域，對此區域進行特徵點提取，改善SIFT演算法花費較多時間 的缺點。 7. 利用重疊區域定位及Harris角點偵測法，先找出重疊的區域，再 利用Harris演算法對區域進行特徵擷取，加快速度與提高準確 度。 文獻上所述以上各方法在抵抗旋轉、尺度和視角變化能力、抵抗 影像扭曲能力、特徵點點數與運算速度之優劣如表 3.1 所示。
20
表 3.1 特徵點擷取方法比較表
演算法 Harris SIFT SURF BRISK
Improved SURF Combined Moment 區域定位 抵抗旋轉、尺度 和視角變化能力 ◎ ◎ ○ ◎ ◎ 抵抗影像扭曲能力 ◎ ◎ ◎ ◎ ◎ ◎ ◎ 特徵點點數 ○ ◎◎ ◎ ○ ◎◎ ◎ ○ 運算速度 ○ ◎ ◎◎ ◎◎ ◎◎ ◎◎ 註：◎◎表示非常良好，◎表示良好，○表示普通，表示不良 從比較的結果可以發現，特徵點演算法改善，主要針對提升運算 速度與增強抵抗視角變化能力為主，但常用的 SURF 演算法演算速度 已相當快速，而 SIFT 雖然花費時間較多，但對一張 700 萬畫素的影像 進行處理也可於幾秒之內完成，且因為拼接航拍圖重疊區域需 70%以 上，重疊區域範圍接近整張影像，所以直接對整張影像提取特徵點即 可，無須費時粗估重疊區域再對此區域進行特徵點提取。基於上述原 因本研究主要對 Harris 角點偵測法、SIFT、SURF 與 BRISK 演算法進 行抵抗旋轉能力與演算速度比較，決定本研究所採用的特徵點提取 法。
21 3.1.2 特徵點提取法抵抗影像旋轉能力評估 當UAV飛航方向改變時，拍攝到的航照圖相較於前一張航照圖， 影像會產生旋轉的現象，為了提高拼接影像的精度，特徵點提取法須 具備抵抗影像旋轉的能力。為了評估Harris、SIFT、SURF及SURF四種 特徵點提取法抵抗影像旋轉的能力，利用python搭配openCV2程式對不 同旋轉角度之影像，進行四種特徵點提取法演算，觀察不同角度下各 演算法所提取的特徵點重複性。 (1) 試驗影像 圖3.2 旋轉影像試體，影像大小為1694×1940pixel，分別為0度(左 上)、30度(右上)、60度(左下)與90度(右下)
22 (2) 實驗步驟
1. 匯入不同角度影像
2. 利用 python 搭配 openCV2 程式進行 Harris、SIFT、SURF 及 BRISK 演算法提取不同角度影像的影像特徵點 3. 將影像分成九宮格，選取每一格內特徵值較大的特徵點 4. 比較不同角度下特徵點位置重複性 (3) 實驗結果與討論 0 度 30 度 60 度 90 度 圖 3.3 Harris 角點偵測法抵抗旋轉能力評估實驗結果，圖中紅色圈位置 為特徵點點位
23 0 度 30 度 60 度 90 度 圖 3.4 SIFT 演算法抵抗旋轉能力評估實驗結果，圖中紅色圈位置為特 徵點點位 0 度 30 度 60 度 90 度 圖 3.5 SURF 演算法抵抗旋轉能力評估實驗結果，圖中紅色圈位置為特 徵點點位
24 0 度 30 度 60 度 90 度 圖 3.6 BRISK 演算法抵抗旋轉能力評估實驗結果，圖中紅色圈位置為 特徵點點位 從圖 3.4、圖 3.5 及圖 3.6 中不同旋轉角度下影像特徵點點位，可 以發現 Harris、SIFT 與 SURF 三種運算方法抵抗旋轉的能力皆相當不 錯，對於 Harris 角點而言，可能角點反應值較大的點抵抗旋轉能力不 錯，反應值越低，抵抗能力越差。從圖 3.6 中可以發現，不同的角度下， BRISK 演算法選取的特徵點不盡相同，證明了 BRISK 演算法抵抗旋轉 能力較弱，也顯示此演算法為了節省時間選用較簡單的運算造成問題 產生。 3.1.3 特徵點提取法運算速度比較 良好的特徵點提取法除了須具備抵抗各種影像變化的能力之外，
25 運算的速度也相當重要，為了評估Harris、SIFT、SURF及BRISK演算 法運算所花費時間，將四種方法對一簡單的幾何圖形影像(圖3.7)進行 特徵點提取運算，比較花費時間。 (1) 試驗影像 圖 3.7 實驗用幾何圖形影像，影像尺寸為 1694×492 pixels (2) 實驗步驟 1. 利用python搭配openCV2程式中Harris、SIFT、SURF及BRISK演 算法對圖3.7進行提取特徵點提取運算 2. 比較四種方法進行特徵點提取演算花費的時間 (3) 實驗結果與討論 從表3.2可以發現，Harris角點偵測法花費的時間相較於SIFT、SURF 與BRISK多，這是因為Harris演算法需對影像中的每一像素進行比較， 需花費較長的時間。 表 3.2 Harris、SIFT、SURF 及 BRISK 特徵點提取法運算時間
特徵點提取法 Harris SIFT SURF BRISK
26 從幾何圖形實驗結果發現，Harris 演算法花費相當多的時間，所以 將其剔除。針對另外三種演算法，選一隨機拍攝的影像(圖 3.8)，將其 縮放不同大小，利用 SIFT、SURF 與 BRISK 進行特徵點提取，觀察三 種演算法所花費的時間，比較它們的演算速度。 (1) 試驗影像 圖 3.8 隨機拍攝影像，改變大小進行演算速度比較，原影像尺寸為 5184× 3456 pixels (2) 實驗步驟 1. 改變影像(圖3.8)大小為100萬、200萬…到700萬畫素，由於影像 超過700萬畫素時，python搭配openCV 2內建之SIFT函數會有錯 誤產生，所以影像最大只取到700萬畫素 2. 利用python搭配openCV2函式庫中SIFT、SURF及BRISK演算法 對不同大小影像進行提取特徵點提取運算 3. 比較三種演算法於不同大小影像所花費的運算時間
27 (3) 實驗結果與討論
從表 3.3 可以發現不論影像大小為何，BRISK 的演算速度較 SIFT 與 SURF 來的快，與表 3.2 結果不同，顯示出雖然處理簡單的幾何圖形 影像是 SIFT 與 SURF 較為快速，但對於一般較複雜的影像，BRISK 演 算法於運算速度上具有相當的優勢，從表 3.3 也可以發現隨著影像變大， 耗時差異也隨著增加，到 700 萬畫素時已有 24%的時間增量，顯示 BRISK 確實是運算速度相當快的演算法。然而為了達到此優點，使用 了較簡略的運算方式，使此方法對於影像旋轉等變化底抗能力較差。 從影像旋轉抵抗能力評估與演算速度比較，發現 SIFT 與 SURF 為較適 合本研究的特徵點提取法，接著詳細介紹兩種特徵點提取法。 表 3.3 不同影像大小，SIFT、SURF 及 BRISK 運算時間
畫素 SIFT(s) SURF(s) BRISK(s)
SURF 運算時間減 BRISK 運算時間 (s) SURF 運算時間對 BRISK 運算時間 增加時間比 7110082 2.36 1.91 0.549 1.361 2.48 6450140 2.15 1.75 0.537 1.213 2.26 5419751 1.78 1.52 0.477 1.043 2.19 4478976 1.47 1.28 0.431 0.849 1.97 3627815 1.21 1.06 0.418 0.642 1.54 2866268 0.95 0.86 0.379 0.481 1.27 1612535 0.58 0.5 0.316 0.184 0.58
28
3.2
SIFT 演算法
SIFT(Scale-Invariant Feature Transform)演算法(Lowe,2004)改良 Harris演算法，改善Harris演算法無法抵抗影像尺度改變的問題，並具 備抵抗影像旋轉的能力。先對影像做不同尺度的縮小，得到一系列縮 小 後 的影 像， 將影 像 尺寸 由大 到小 排 列構 築成 影像 金 字塔 (Image pyramid)，接著對金字塔中每一層的影像做高斯平滑，減少影像圖像雜 訊以及降低細節層次，增強影像在不同比例大小下的影像效果，如邊 緣效果，再將相鄰的高斯影像相減，做高斯差分( Difference of Gaussian, DOG )，利用DOG來近似高斯拉普拉斯(Laplace of Gaussian, LOG)，LOG 可用來檢測影像邊緣，但運算量較大，所以利用DOG來近似，減少運 算時間。 利用影像金字塔搭配高斯差分，使特徵點具備尺度不變的能力， 接著對高斯差分影像做極值檢測，找出區域極值為特徵點，然後過濾 低對比度特徵點，最後進行特徵點向量化描述，描述特徵點的梯度方 位及強度值。其方法步驟如下: (1) 假設L( x, y, σ )表示不同高斯尺度下的圖像，G( x, y, σ )為不同尺度下 的高斯函數，I( x, y )為輸入的影像灰階值，σ為高斯影像尺度改變參 數，Lowe在其研究中進行不同尺度下特徵點的重複性，得出最佳值 為1.6，L( x, y, σ )與G(x, y, σ )表示如下: 2 2 2 ( )/2 2 1 ( , , ) 2 x y G x y e (3.1) ( , , ) ( , , ) ( , ) L x y G x y I x y (3.2)
29 (2) 將金子塔每一層中相鄰的高斯影像做高斯差分( DOG )(圖 3.9)，DOG 可用來找出影像的邊緣與斑點位置，DOG 的算法如式(3.3)，其中 k 為常數，設定 1/ 2 s k ，s 為影像尺度
( , , ) ( ( , , ) (x, y, )) , ( , , ) ( , , ) D x y G x y k G I x y L x y k L x y (3.3) 圖 3.9 將影像依尺度倍率以 2 倍為一單位分組，並縮小影像為原來的一 半，每一組別中進行高斯平滑，再將同一組相鄰的經高斯平滑處理後 的影像兩兩相減可得其 DOG 影像(Lowe,2004) (3) 利用同一組內 DOG 金字塔中相鄰兩層圖像間進行比較可以獲得初 步的特徵點，若某像素的 DOG 值與相鄰的 26 個像素相比為極大值 或極小值時，則將此像素視為區域極值所在的位置，也就是初估的 特徵點，如圖 3.10。為了在每組中能檢測 S 個尺度的極值點，則 DOG30 金字塔每組需 S+2 層圖像，而 DOG 金字塔是由高斯金字塔相鄰兩 層相減得到，所以高斯金字塔每組需 S+3 層圖像。 圖 3.10 中間層的檢測點和它同尺度的 8 個相鄰點以及上下相鄰尺度的 9 × 2 個點，共 26 個點進行比較，以確保在尺度空間和二維圖像空間都 檢測到極值點(Lowe,2004) (4) 為了提高特徵點的精度到子像素，並將低對比度或將不夠顯著的 特徵點位做過濾，利用 DOG 函數於尺度空間中以初估特徵點為原 點的泰勒展開式(曲線擬合)D(x)直到二次項，計算特徵點的子像素 偏移量，其中 x 為到此初估特徵點的偏移量: 2 2 1 ( ) D 2 T T D D D x x x x x x (3.4) 將 D(x)對 x 做一次微分並令它等於零，得到極值ˆx 1 2 2 ˆ D D x x x (3.5)
31 當ˆx在任一維度的絕對值大於 0.5，表示極值已偏移到鄰近點上， 故會將原本極值點移往該點位，並重新計算子像素位置與差值，反 覆迭代直到收斂。 接著計算D x( )位於ˆx的值D x( )ˆ ，D x( )ˆ 過小則特徵點容易受雜訊影 響，所以將其剔除，D x( )ˆ 的式子如下 ˆ 1 ˆ ( ) 2 T x D D x D x (3.6) 文獻中設定D x( )ˆ 門檻值為 0.03，若值小於 0.03 則表示此特徵點的 穩定性不足，將予以刪除。 (5) 因為 DOG 函數對於偵測邊緣上的相當敏感，但邊緣上的點穩定性 較差，為了刪除位於邊緣上的特徵點，可以利用特徵點的主曲率 進行評斷，對於邊上的特徵點，因為穿過邊緣方向的主曲率會遠 大於沿著邊緣方向上的主曲率，因此穿過邊緣方向與沿著邊緣方 向的主曲率比值遠大於位於角落的比值。為了計算主曲率比值， 可以利用 Hessian matrix 的特徵值比例與主曲率比例相同的特性來 去除位於邊緣上的點。Hessian matrix 如下: xx xy xy yy D D H D D (3.7) 矩陣的特徵值α與β代表 x 和 y 方向的梯度
2 xx yy xx yy xy Tr H D D Det H D D D (3.8)32
Tr H 代表矩陣對角元素之和，Det H
代表矩陣行列式值，令 ，
2
2 2 2 2 1 Tr H Det H (3.9) 為了剔除邊緣特徵點，定一門檻值th ，通常為 10，需滿足(3.10) 才保留
2 2 1 th th Tr H Det H (3.10) (6) 為了使特徵點俱備旋轉不變性，於特徵點描述前，先找出特徵點 主梯度方向，利用圖像局部特徵給每一個特徵點分配一個基準方 向做為特徵點向量，以相鄰像素的梯度方位及強度值找出特徵點 的主梯度方向。m x y
, 梯度強度及
x y, 梯度方向分別如式(3.11) 及(3.12)所示。 2 2 ( , ) ( ( 1, ) ( 1, )) ( ( , 1) ( , 1)) m x y L x y L x y L x y L x y (3.11) 1 (x, y 1) L(x, y 1) ( , ) tan ( ) (x 1, y) L(x 1, y) L x y L (3.12) 完成梯度的強度與方向計算後，使用直方圖統計鄰域內像素的梯度 和方向。梯度直方圖將 0~360 度的方向依每柱 10 度，分為 36 個柱， 為了減低離特徵點較遠點位的影響程度，每個相鄰像素依據其高斯 圓形權重(強度值大小與梯度方向)加入特徵點直方圖中，最後直方 圖中峰值的方向即為此特徵點的方向，若直方圖統計中，有其他值 大於或等於峰值的 80%時，代表此點可能屬於多方向的情形，此時33 新增一個與此點同樣位置同樣尺度，但主要方向為不同方向的特徵 點，通常僅有 15％的特徵點被賦予多個方向。 (7) 為特徵點建立一個描述子向量，使特徵點在不同光線與視角下能保 持其不變性，以及使特徵點間可以更容易區別，做為影像間特徵點 匹配依據。以特徵點為中心取 16×16 窗口，再將其分為 4×4 子窗口， 計算每個子窗口內 8 個維度(即 0、45、90、135、180、225、270、 315 度)的梯度方向直方圖，最終以 4×4×8 共 128 維的向量來表示(圖 3.11)。特徵向量形成後，為了去除光照變化的影響，需要對它們進 行正規化處理，假設得到的描述子向量為H
h h1, 2,...,h128
，正規化 後的特徵向量為L
l l1, ,...2 l128
則 128 1 , 1, 2,3,... i i j j h l i h
(3.13) 為了減少非線性亮度的影響，把大於 0.2 的向量值設為 0.2，再進行 一次正規化處理，提高特徵的鑑別性。 圖 3.11 特徵點描述示意圖(Lowe,2004)34
3.3
SURF 演算法
SURF（Speeded Up Robust Feature）演算法（Bay et al.,2006）是
一改善SIFT演算法的特徵點提取法，穩健性高且運算快速的方法，此 方法快速的主要原因為，先製造出積分影像與利用Hessian matrix探測 器找出特徵點，接下來藉由改變filter大小來創造尺度空間快速獲得邊 界 和 小 斑 點 ， 之 後 利 用 3× 3× 3 的 相 鄰 非 極 大 抑 制 (Non-Maximum Suppression)來使影像中的特徵點停留於一處，最後利用x、y方向的一 階Haar wavelet response來描述特徵點。方法步驟如下:
(1) 積分圖像 利用積分圖像概念，將圖像與高斯二階微分模板的濾波轉化對積分 圖像作加減運算，積分圖像為先將影像的每個像素轉為灰階強度 I， 而每一個點 X(x,y)的積分值 IΣ(X)如式(3.14)。 0 0
( , )
j y i x i jI
I i j
(3.14) 再計算任一窗口內的積分值(圖 3.12)35 圖 3.12 積分影像區塊值（Bay et al.,2006） (2) 利用 Hessian matrix 行列式的極大值檢測斑點方法，給定圖像 I 中 任意一個點p x y ,在點 p 處尺度為 σ 的 Hessian matrix H(p,σ)定義
, 為: ( , ) ( , ) ( , ) ( , ) ( , ) xx xy xy yy L p L p H p L p L p (3.15) (3.15)式中 Lxx(p,σ)式為高斯二階微分 2 ( ) 2 g x 在點 p 處與圖像 I 卷積， 其中 g(σ)代表不同尺度下的高斯函數，Lyy(p,σ)和 Lxy(p,σ)有相同涵 義。 由於二階高斯微分模板被離散化和裁剪，導致圖像在旋轉 π/4 的奇數倍時，即轉動到模板的對角線方向時，特徵點檢測的重複性 降低。而在π/2 時，特徵點檢測的重複性最高，但這一小小的不足 不影響我們使用 Hessian matrix 進行特徵點檢測。36 為了將模板與圖像的卷積轉化成盒子濾波(Box Filter)運算，需 對高斯二階微分模板進行簡化，使得簡化後的模板只是由幾個矩形 區域組成，矩形區域內填充同一值，如圖 3.13 所示，在簡化模板 中白色區域的值為 1，黑色區域的值為-1，灰色區域的值為 0。 圖 3.13 高斯二階微分模板及簡化（Bay et al.,2006） 對於σ=1.2 的高斯二階微分濾波，設定模板的尺寸為 9×9 大小， 並用它作為最小尺度空間值對圖像進行濾波和斑點檢測。使用 Dxx、 Dxy與 Dyy，表示模板與圖像進行卷積的結果，將 Hessian matrix 的 行列式簡化如公式(3.16) 2 det(Happrox)D Dxx yy(wDxy) (3.16) 其中相對權重 w 如下: (1.2) (9) 0.912 0.9 (1.2) (9) xy F yy F yy xy F F L D w L D (3.17) (3) 尺度空間表示 尺度空間表示想要獲取不同尺度的斑點，必須建立圖像的尺度 空間金字塔。SURF採用盒子濾波和積分圖像，因此不需像SIFT演
37
算法直接建立圖像金字塔，而是採用不斷增大盒子濾波模板尺寸的 間接方法(圖3.14)。透過不同尺寸盒子濾波模板與積分圖像求取 Hessian matrix行列式的反應圖像(response map)，然後在反應圖像 上採用3D非極大抑制，求取各種不同尺度的斑點。
圖 3.14 SURF 與 SIFT 提取法濾波器比較圖（Bay et al.,2006）
(4) 特徵點方向分配
為了進行影像特徵點匹配，需從特徵點中提取特徵向量，利用 向量相似程度判斷影像間互相匹配的特徵點，因此需要先計算特徵 點的方向，首先以特徵點為中心，以 6s(s 為特徵點的尺度)為半徑 的圓形區域內，對圖像進行 Haar wavelet response。即利用積分圖 像對影像進行梯度運算，提高計算圖像梯度的效率。依照距離遠近 賦於高斯權重係數，使越靠近特徵點的 Haar wavelet response 值貢 獻越大，距離越遠的 Haar wavelet response 貢獻越小，如此特徵點 描述可用一坐標（x,y）表示，其中 x 表示 x 方向的響應強度，y 表 示 y 方向上的響應強度。接著計算主方向，採用一大小為 π/3 的滑
38
動窗口，計算窗口內水平和垂直 Haar wavelet response，相加形成 一局部方向向量，遍歷整個圓形區域後，選擇最長的向量方向為該 特徵點主方向。
(5) 藉由 Haar wavelet response 提取描述子
以特徵點為中心產出一個長為 20s(s 為特徵點的尺度)的區域且以 前一步驟產生的方向為方向，如圖 3.15 圖 3.15 影像中的分割區域（Bay et al.,2006） 再將這些區域分成 4×4 個子區域然後每個子塊利用尺寸 2s 的 Haar wavelet 模板進行響應值計算，然後對響應值進行統計dx、 dx 、dy、 dy 形成特徵向量，如圖 3.16 所示。
39 圖 3.16 特徵描述表示（Bay et al.,2006） 然後以特徵點為中心，對 dy 和 dx 進行高斯加權計算，其 σ=3.3s。 最後，分別對每個子塊的響應值進行統計，得到每個子塊的向量
,
,
,
v
dx dy
dx
dy
(3.18) 共有個 4×4 的子塊，因此特徵描述子由 64 維特徵向量組成。3.4
數位影像相關係數法
數位影像相關係數法(Digital Image Correlation method，DIC)的原 理是利用分析物體表面的灰階分佈特徵，對變形前後的影像進行特徵 的比對來獲得變形前後影像的相對位置，並藉此相對位置計算出影像 中任一位置的位移向量，其它物理量如應變場、剪應變都可以利用此 位移向量求得。 數位影像相關係數法常應用於影像識別技術等相關領域，以物體 表面上的結構性斑紋為影像特徵，結構性斑紋如圖 3.17 中散佈的黑色
40 斑點所示，此特徵則被用來做為分析比對變形前後兩張影像局部相關 性的依據。我們以變形前後兩張影像的局部相關性為目標函數，此相 關性影響變形前後兩張影像的局部對應關係，故以對應函數的參數為 最佳化的變數，以數值方法進行最佳化求解，即當局部相關性最佳時， 其對應函數的參數就是最接近實際值的參數，茲以圖 3.18 為例說明其 對應關係，圖中變形前中心點位置在點 P，變形後點 P 位置改變為 * P ， 兩者之間的函數關係如式(3.19)、(3.20)所示： *
( , )
x
x u x y
(3.19) *( , )
y
y v x y
(3.20) 其中 ( , )x y 為變形前影像座標， * * ( ,x y )為變形後影像座標， ( , )u x y 、 ( , ) v x y 為變形量。 圖 3.17 物體表面之結構性斑紋與次級影像（網格）示意圖（施等,2006）41 圖 3.18 物體表面上之次級影像變形前後相對位置示意圖（施等,2006） 對於未變形的影像，可以利用有限元素法的觀念將其分割成數個 次級影像，如圖 3.17 中一個方形的格子即代表一個次級影像。假設變 形前影像為影像 A，變形後影像為影像 B，二者之間存在著前述的對 應關係，則我們可以用數位影像相關係數來判定二者間的相關程度， 數位影像相關係數的定義（施等,2006）如下式： 2 2 ij ij ij ij g g COF g g
(3.21) 其中g
ij 及g
ij 分別為影像 A 在( , )
i j
座標上及影像 B 在( , )i j 座標上的 灰階值，而( , )i j 座標為影像 A 上( , )
i j
座標依函數計算後所得在影像 B 上之對應座標 以最佳化程序求解出每一次級影像最大相關係數來決定對應函數 的最佳參數，即可獲得次級影像變形前後的相對應座標，再運用這些 相對應座標即可計算出個別的位移量，進而獲得位移場與應變場。42
第四章 空中三角測量
空中三角測量用於解算航照影像的各種資訊，例如相機內外方位、 地面座標等，目前常藉由 UAV 或是空拍機拍攝地面影像，再拼接影像 製作全景圖或是 3D 全景圖等用途，亦可作為還原影像的地面座標的方 法，藉此製作三維地形模型，本研究主要利用 Pix4Dmapper 軟體進行 相關空中三角測量解算，以下介紹空中三角測量原理、相機參數及套 裝軟體 Pix4Dmapper。4.1
相機參數
相機參數包含了內方位與外方位，內方位參數是相機與鏡頭內部 的機組零件所產生，例如焦距、像主點偏移量與扭曲參數，其中扭曲 參數是用來將因鏡頭上的透鏡造成的影像扭曲變形座標轉換成未扭曲 座標的參數；外方位參數則是當相機進行拍攝時所在的空間座標以及 其對於XYZ軸的旋轉量。 1. 內方位參數：指相機內部之參數，包含了像主距 f 、透視中心到像 主點之偏移量X Y0, 0 與鏡頭所產生的扭曲參數。 在扭曲校正的實驗中，通常需要求得扭曲係數以還原影像，而扭曲 係數會因扭曲方程式定義不同而異，在此引用Weng等1992年提出之 扭曲方程式，表示如下： ' ( , ) x x x
x y (4.1) '( , y)
yy
y
x
(4.2)43 其中 ( , )x y 為理想影像座標(無扭曲)， ' ' ( ,x y)為拍攝影像座標(有扭 曲)， ( , x y)為x、y方向之鏡頭總扭曲變形量。 (1) 徑向扭曲（Radial distortion） 鏡頭扭曲絕大部分是因為徑向扭曲的關係，其所造成之扭曲會 使得影像沿著由中心點往外發散的徑向直線變形，令座標點拉 近中心或是遠離中心。 徑向扭曲是徑向距離的函數，在任何角度都是一樣，所以徑向 扭曲也是一個對稱的扭曲，也叫做輻射扭曲。扭曲變形的程度 為座標點到中心距離的函數，其方程式如下： p x x x (4.3) p y y y (4.4) 2 2 2 r x y (4.5)
2 4 6
1 2 3 ... xr k r k r k r x (4.6)
2 4 6
1 2 3 ... yr k r k r k r y (4.7) xr 、yr為逕向扭曲變形量，r為到扭曲中心的距離，x
p、yp 為像主點座標k1、k2、k3為逕向扭曲參數。 (2) 偏心扭曲(Decentering distortion) 由於相機中的透鏡並未完全的共線，使得實際的光學系統受到 不同程度的偏心。造成像點位移的方向垂直於徑向線，造成所 有的通過中心點的徑向線彎曲變形所以此扭曲包含徑向扭曲 與切向扭曲。唯有某一個方向上的徑向線沒有變形的情況，這44 條線即被稱為最小切線向扭曲軸，而與其垂直的則為最大切線 向扭曲軸，如圖 4.1 所示，偏心扭曲之變形量如下：
2 2
1 2 2 2 xp p r x p xy (4.8) 2 2 2( 2 ) 2 1 yp p r y p xy (4.9) xp
、
yp為偏心扭曲量， p1、 p2為偏心扭曲參數。 圖 4.1 偏心扭曲之最大切線扭曲軸及最小切線扭曲軸（Weng,1992）(3) 薄菱鏡扭曲(Thin prism distortion)
因為非完美的鏡片設計或者製造時使得鏡片或感光元件些微 傾斜。薄菱鏡扭曲會造成類似偏心扭曲的影像變形，也包含了 徑向扭曲以及切線向扭曲。 其偏心扭曲之變形量如下：
2 2
1 xp s x y (4.10)45
2 2
2 vp s x y
(4.11) xp 、yp為偏心扭曲量，s
1、s
2為薄菱鏡扭曲參數。 2. 外方位參數: 是指相機在地面座標內之位置與方位，相機相對於 XYZ 方向之旋轉角、、k 與真實座標XL、YL、ZL。可利用 GPS 定位、測量或是共線式等方式去求得。4.2
空中三角測量原理
空中三角測量基本原理為一共線式，即透視中心 o 之空間座標
X Y ZL, L, L
與影像座標
x y0, 0
、任意物點 A 之影像座標
x ya, a
、任意物 點 A 在真實空間座標
X Y ZA, A, A
，理論上位於同一條線上，如圖 4.2 所 示(何維信,1995)，利用此共線原理可得出式(4.12)與(4.13)，其中 f 為 焦距，m為對於 XYZ 軸的旋轉參數如式(4.14)，旋轉角 , , 為相機相 對於 XYZ 軸旋轉。
11 12 13 0 31 32 33 A L A L A L a A L A L A L m X X m Y Y m Z Z x x f m X X m Y Y m Z Z (4.12)
21 22 23 0 31 32 33 A L A L A L a A L A L A L m X X m Y Y m Z Z y y f m X X m Y Y m Z Z (4.13) 11 12 13 21 22 23 31 32 33cos cos cos sin sin sin cos sin sin cos sin cos
cos cos cos cos sin sin sin sin cos cos sin sin
sin sin cos cos cos
m m m m m m m m m m (4.14)
46 圖 4.2 空中三角原理示意圖 空中三角測量常用的數學模式主要分為三種，一種是光束法平差， 其次是自率光束法平差，最後為附加參數光束法平差。光束法即為式 (4.12)與(4.13)而自率光束法是利用非量測型相機取得影像後，相機可以 事先透過率定取得內方位參數，也就是 f 、x0、y0視為定值並帶入共線 式解算。因為當內方位參數不夠正確時，其誤差會轉移到其他參數， 因此通常單純使用光束法平差時，解算的結果較不理想。而自率光束 法將內方位參數為定值，解算時不會一起做誤差的分配，通常可有較 均勻的平差效果；附加參數自率光束法將相機內方位參數為定值與外 方位參數及鏡頭扭曲參數於平差計算時一同解算，可提升測量平差的 精度，加入附加參數後共線式改寫為下式(4.15)與(4.16)，x、y為影 像座標改正的扭曲量。
11 12 13 0 31 32 33 A L A L A L a A L A L A L m X X m Y Y m Z Z x x f x m X X m Y Y m Z Z (4.15)47
21 22 23 0 31 32 33 A L A L A L a A L A L A L m X X m Y Y m Z Z y y f y m X X m Y Y m Z Z (4.16) 利用共線式原理，代入已知條件去以梯形法求出未知參數之結果， 當有控制點之影像座標與控制點之真實座標即可求得內外方位參數： 當內外方位參數與控制點之影像與真實座標皆有時，可求得其他影像 座標之真實座標，此方式可應用至整張影像進而產生三維地形模型。 在進行航空影像拍攝時，空間中必須有平面及高程控制點，藉由這些 控制點來還原影像在空間上之三維座標，進而求得拍攝區之任意點點 位地面座標，且必須確保影像中有足夠的控制點以及一個拍攝航道中 每張影像需重疊足夠的控制點，控制點越多可確保結果越精準。經由 量測得到地面控制點座標與利用數位影像相關係數法分析控制點在影 像中的影像座標後代入共線式，利用最佳化程序可求得未知內外方位 與扭曲係數；或是得到已知內外方位反推影像之真實座標，利用影像 的真實座標製作三維地形模型。相機內外方位可以從空中三角測量過 程中得出，外方位參數的近似值亦可由 GPS 記錄提供，而地面控制點 即為已知的地面真實座標位置，經由測量方式或是 GPS 可量測設置。4.3
Pix4Dmapper 軟體
Pix4Dmapper 為一套專業航測影像處裡軟體，能快速、自動化處理 任何相機拍攝之航空影像、傾斜影像、地面影像、支援各種波段，能 達公分級的三維精度，另外能迅速產製正射影像、數值地形模型、外48 方位參數等成果。軟體的特徵點提取與匹配方法，為 SIFT 演算法(蘇 柏軒，2014)、(莊大賢，2015)。 軟體便利且自由度高，使用者可以依需求匯入不同資料，其中一 項功為可以匯入手動連結點輔助影像拼接操作，這是其他軟體沒有的 功能，利用這項功能可以將數位影像相關係數法與特徵匹配法所演算 出的點匯入軟體中，觀察匯入前後拼接影像與三維模型，比較拼接準 確度與精度做為評估的依據，此特點為本研究選用這套軟體最主要的 原因。
49
第五章 實驗驗證與討論
當災害發生後，為了能及時救災與當邊坡滑動時能準確計算崩塌 量，精確的拼接影像與三維數值模型是必須的，本研究的目的在於利 用DIC來增加影像拼接與三維數值模型的精度，為了探討DIC對影像拼 接與三維數值模型精度的影響，在本章將先就影像匹配精度與其他常 用的影像匹配方法進行比較。接著，利用三維階梯模型與縮尺地形模 型進行三維量測精度的分析，分別就DIC用在控制點定位與匯入DIC連 結點與SURF連結點，對於Pix4Dmapper軟體產製的影像拼接與三維模 型精度影響。5.1
影像匹配演算法精度比較實驗
特徵點的匹配準確度對於影像拼接及三維模型精度有著相當的影 響，為了比較DIC與其他研究及實務方面較常用的影像匹配法（SIFT、 SURF、BRISK）的定位精度差異，本實驗先平移一張空拍影像，分別 往x方向與y方向移動，接著利用上述四種方法進行影像匹配，比較匹 配結果，評估演算法精度。以下詳細介紹本實驗試體、實驗步驟、結 果與討論。 5.1.1 試體 為了比較DIC與其他匹配方法的精度，並確保影像只有單純的位移， 不會有其他如影像扭曲、光照不同等影響，所以選取一戶外UAV空拍 影像(圖5.1)，影像大小為2918× 2189畫素，作為原始影像，再將同一張50 影像往x方向移動123、123.1…每次多移動0.1pixel直到125pixel，共21 張影像，y方向則移動234、234.1…每次多移動0.1pixel，直到236pixel， 共21張影像，x與y方向影像共42張，作為目標影像，供不同影像匹配 方法進行匹配實驗。 圖 5.1 原始影像，影像大小為 2918×2189 畫素 5.1.2 實驗步驟
利用五種影像匹配方法，分別為 SIFT 演算法、SURF 演算法、SIFT 特徵點搭配 DIC、SURF 特徵點搭配 DIC 及 BRISK 對移動前後兩張影 像進行特徵點匹配實驗，再計算移動前與移動後匹配點的 x 方向與 y 方向之相對距離，獲得原始影像與目標影像相對位移，比較各種匹配
