小波理論應用於影像套合之研究

全文

(1)小波理論應用於影像套合之研究. 學生：陳維勇. 指導教授：史天元. 國立交通大學土木工程學系摘要本研究藉由小波轉換進行特徵點萃取，並與 TDGO 興趣運算元比較，測試兩者對於旋轉量及雜訊之反應。續以改良自 Hsieh et al. (1997)所發展，應用小波特徵萃取之影像套合流程，與洪偉嘉（1999）所發展，以 TDGO 運算元配合最小二乘匹配法（LSM）之流程，進行數值測試並分析比較。經實驗發現，就特徵萃取而言小波方法對雜訊之忍受度較 TDGO 運算元為高；對影像間旋轉量之表現差異則並不明顯。就套合流程而言，洪偉嘉（1999）之流程雖具備較嚴密之特徵點分佈控制，但仍出現鑲嵌成果精度不均之現象。改良自 Hsieh et al. (1997)之流程控制點配對品質不如洪偉嘉（1999）之流程，但經適當門檻值設定後，可達成完全自動化之成果。. I.

(2) A Study on the Wavelet Theory Applied to Image Registration Student : Wei-Yung Chen. Advisor : Tian-Yuan Shih. Department of Civil Engineering National Chiao Tung University. Abstract. In this study, a wavelet-based edge detector is adopted for extracting feature points from image.. In comparison with the TDGO interest operator, the. performance on noise tolerance and rotation invariance is investigated. Furthermore, the image registration scheme derived from Hsieh et al. (1997) with a wavelet-based edge detector as well as the one developed by Hong (1999) with the TDGO interest operator and the least square matching method are tested with real images. The wavelet-based method is proven to have higher noise tolerance than the TDGO operator.. In regard to rotation invariance, the performances of both. methods are about the same.. The scheme developed by Hong (1999) is. inplemented with a strict control over feature point distribution, however, local inaccuracy still can be observed. The scheme derived from Hsieh et al. (1997) provides fully automation with properly given thresholds.. II.

(3) 誌謝本論文承蒙恩師史天元教授兩年來在各方面的悉心指導，使論文得以順利完成，您對學術的嚴謹態度與研究精神令學生深感欽佩，在此謹向恩師致上最誠摯的謝意與最崇高的敬意。感謝組上陳春盛教授、黃金維教授及李振燾教授在相關課業的啟發與教誨。以及口試委員王蜀嘉教授、陳良健教授、曾義星教授及蔡展榮教授提供寶貴的意見，使本論文更臻完善。感謝宏揚學長陪我走過這最艱苦的半年，不但是我研究上最佳伙伴，也是我連線遊戲的最佳戰友。嘉聲、立信、松安及淼祥學長適時地於課業及研究瓶頸中給予協助，在此一併致謝。感謝交大羽球隊教練廖威彰老師，羽球隊的歷練與梅竹賽的經驗，是我學生生涯中最重要也最難忘的一段歷程。當然，羽球隊嚴苛的體能訓練，也是我夜以繼日、不斷前進的基石。我的好戰友阿進、壕哥、泰哥、小中中、阿衡、阿賓、白猴、大呆、老方、胖瑋、花輪，以及所有可愛的學長姐弟妹們，感謝你們陪我走過這一段難忘的日子。感謝我情同兄弟的好同學傳盛、程傑、輝欽、信璁、奕峰、育智及欣玉，陪我走過每一個苦悶、歡欣、糜爛的日子。兩年研究生活中，也少不了同窗好友宇俊、麟霂、建廷、宗融及振偉在課業、生活上的幫助與鼓勵，謝謝你們。好了，我知道一定漏了一些人，別罵我，不管你們在哪裡，我都謝謝你們。. 最後，謹以本論文獻給我的家人，以及所有關心、支持我的朋友們。. III.

(4) 目錄中文摘要 ........................................................................................................................I 英文摘要 ...................................................................................................................... II 誌. 謝 .........................................................................................................................III. 目. 錄 .........................................................................................................................IV. 圖目錄 .........................................................................................................................VI 表目錄 ......................................................................................................................VIII 第一章. 前言 ........................................................................................................... 1. 1-1 研究動機與目的 ............................................................................................ 1 1-2 研究方法及內容 ............................................................................................ 1 1-3 論文架構 ........................................................................................................ 2 第二章. 特徵萃取 ................................................................................................... 6. 2-1 TDGO 運算元 ................................................................................................ 6 2-2 小波特徵萃取 ................................................................................................ 9 2-2-1 小波轉換 .................................................................................................. 9 2-2-2 Canny 邊緣偵測法 ................................................................................. 10 2-2-3 多尺度（Multiscale）邊緣偵測法........................................................ 10 2-2-4 特徵點選取 ............................................................................................ 12 2-3 特徵萃取成果比較分析 .............................................................................. 14 2-3-1 標準測試 ................................................................................................ 14 2-3-2 旋轉量測試 ............................................................................................ 18 2-3-3 雜訊測試 ................................................................................................ 23 2-3-4 成果歸納與分析 .................................................................................... 27 第三章. 影像匹配 ................................................................................................. 29. 3-1 角度累計直方圖（ANGLE HISTOGRAM）................................................... 29 3-2 標準化互相關法（NORMALIZED CROSS-C ORRELATION, NCC） ............... 30 3-3 最小二乘匹配法(LEAST-SQUARES MATCHING, LSM)................................. 31 3-3-1 LSM2 平差模式 ..................................................................................... 31 3-3-2 LSM2 運算流程 ..................................................................................... 34 IV.

(5) 3-4 相似度量測 .................................................................................................. 36 3-5 配對一致性檢查 .......................................................................................... 36 3-6 TAU 測試...................................................................................................... 38 第四章. 數值測試 ................................................................................................. 40. 4-1 航照立體像對 .............................................................................................. 40 4-2 近景相片 ...................................................................................................... 49 4-3 近景數值影像 .............................................................................................. 54 第五章. 結論與建議 ............................................................................................. 59. 參考文獻 ..................................................................................................................... 61 附錄 A：Spline 小波原理[Mallat, 1992].................................................................. 63 作者簡歷 ..................................................................................................................... 66. V.

(6) 圖目錄圖 1- 1 HSIEH ET AL. (1997)小波影像套合流程 ........................................................3 圖 1- 2 本研究採用之小波影像套合流程...............................................................4 圖 1- 3 洪偉嘉（1999）興趣運算元影像套合流程..............................................5 圖 2- 1 影像之位元型式編碼...................................................................................7 圖 2- 2 定義特徵點之圖型（PATTERN）.................................................................8 圖 2- 3 逐步擴大比較窗之非局部最大壓縮...........................................................9 圖 2- 4 一維訊號邊緣相關計算成果[HSIEH ET AL., 1997]....................................12 圖 2- 5 ISPRS 立體像對 ..........................................................................................15 圖 2- 6 ISPRS 立體像對經 TDGO 運算元萃取之特徵點位圖 .............................15 圖 2- 7 ISPRS 立體像對經小波方法萃取之特徵點位圖 ......................................16 圖 2- 8 翡翠水庫衛星影像特徵萃取成果圖.........................................................17 圖 2- 9 含旋轉量之樣本影像.................................................................................19 圖 2- 10 含旋轉量影像經 TDGO 運算元萃取之特徵點位圖 ..............................20 圖 2- 11 含旋轉量影像經小波方法萃取之特徵點位圖 .......................................21 圖 2- 12 含旋轉量影像之邊緣強度影像...............................................................22 圖 2- 13 含 GAUSSIAN 雜訊之樣本影像 ................................................................24 圖 2- 14 含雜訊影像經 TDGO 運算元萃取之特徵點位圖 ..................................25 圖 2- 15 含雜訊影像經小波方法萃取之特徵點位圖............................................26 圖 2- 16 含雜訊影像之邊緣強度影像..................................................................27 圖 3- 1 雙線性內插重新取樣................................................................................35 圖 4- 1 航照立體像對.............................................................................................40 圖 4- 2 航空立體像對經洪偉嘉（1999）之流程獲得控制點對成果圖.............41 圖 4- 3 航照立體像對經洪偉嘉（1999）之流程殘差分佈圖.............................41 圖 4- 4 航照立體像對掃瞄影像經洪偉嘉（1999）之流程鑲嵌合成圖.............43 圖 4- 5 航照立體像對小波特徵萃取成果圖.........................................................44 圖 4- 6 航照立體像對特徵點初步配對成果圖...................................................44 圖 4- 7 航照立體像對特徵點配對初步除錯成果圖.............................................45 圖 4- 8 航照立體像對特徵點配對 TAU 測試成果圖 ...........................................45 VI.

(7) 圖 4- 9 航照立體像對經小波流程殘差分佈圖.....................................................46 圖 4- 10 航照立體像對掃瞄影像經小波流程鑲嵌合成圖...................................47 圖 4- 11 航照立體像對鑲嵌成果局部放大比較圖 ...............................................48 圖 4- 12 野外地質調查拍攝之近景相片...............................................................49 圖 4- 13 近景相片經洪偉嘉（1999）之流程獲得控制點對成果圖...................49 圖 4- 14 地質調查近景相片經洪偉嘉（1999）之流程殘差分佈圖....................50 圖 4- 15 地質調查掃瞄影像經洪偉嘉（1999）之流程鑲嵌合成圖...................50 圖 4- 16 地質調查近景相片小波特徵萃取成果圖...............................................51 圖 4- 17 地質調查近景相片特徵點初步配對成果圖............................................51 圖 4- 18 地質調查近景相片特徵點配對初步除錯成果圖....................................52 圖 4- 19 航照立體像對特徵點配對 TAU 測試成果圖 ..........................................52 圖 4- 20 地質調查近景相片經小波流程殘差分佈圖...........................................53 圖 4- 21 地質調查掃瞄影像經小波流程鑲嵌合成圖...........................................53 圖 4- 22 近景數值影像...........................................................................................54 圖 4- 23 近景數值影像經洪偉嘉（1999）之流程獲得控制點成果圖...............54 圖 4- 24 近景數值影像經洪偉嘉（1999）之流程殘差分佈圖...........................55 圖 4- 25 近景數值影像經洪偉嘉（1999）之流程鑲嵌合成圖...........................55 圖 4- 26 近景數值影像小波特徵萃取成果圖.......................................................56 圖 4- 27 近景數值影像特徵點初步配對成果圖...................................................56 圖 4- 28 近景數值影像特徵點配對初步除錯成果圖...........................................57 圖 4- 29 近景數值影像特徵點配對 TAU 測試成果圖 .........................................57 圖 4- 30 近景數值影像小波流程殘差分佈圖.......................................................58 圖 4- 31 近景數值影像經小波流程鑲嵌合成圖...................................................58. VII.

(8) 表目錄表 2- 1 ISPRS 立體像對特徵萃取成果記錄 ..........................................................14 表 2- 2 翡翠水庫衛星影像特徵萃取成果記錄.....................................................16 表 2- 3 含旋轉量影像特徵萃取成果記錄..............................................................18 表 2- 4 含雜訊影像特徵萃取成果紀錄.................................................................23 表 4- 1 航照立體像對經洪偉嘉（1999）之流程求解八參數係數表.................42 表 4- 2 航照立體像對經小波流程求解八參數係數表.........................................46 表 4- 3 地質調查近景相片經洪偉嘉（1999）之流程求解八參數係數表.........50 表 4- 4 地質調查近景相片經小波流程求解八參數係數表.................................53 表 4- 5 近景數值影像經洪偉嘉（1999）之流程求解八參數係數表.................55 表 4- 6 近景數值影像經小波流程求解八參數係數表.........................................58. VIII.

(9) 第一章前言 1-1 研究動機與目的影像套合（Image Registration）是指將含有同一景物之二張或多張影像進行幾何糾正（Geometric Correction），以消除各影像間因攝影條件不同而出現之幾何差異。一般而言，影像套合處理可分為兩個重要步驟： (1). 在兩共軛影像間選取適當分佈且數目足夠的控制點（Control Point），或稱之為共軛點（Conjugate Point），並分別量測這些控制點之影像座標。. (2). 利用上述控制點座標解算兩共軛影像間之映射函數（ Mapping Function），利用該映射函數進行座標轉換及影像重新取樣。. 就利用特徵進行影像匹配之程序而言，步驟(1)中所指控制點即為影像中之特徵點，一般可透過興趣運算元（Interest Operator）獲取。本研究擬應用小波邊緣偵測法，萃取位於影像邊緣線上之特徵點，並輔以前述興趣運算元萃取得之特徵點，進行比較分析。續以兩套分別利用小波方法及興趣運算元配合最小二乘匹配法之匹配流程，進行影像套合運算，期瞭解小波特徵萃取於影像套合之適用性及優缺點。. 1-2 研究方法及內容 Schenk (1997)將影像匹配技術大致分為三類：第一類為區域匹配（ Areabased Matching），區域匹配實際上即是利用影像灰度值進行匹配，藉由分析影像間區域性的灰階分佈狀況，進行相似度量測以為匹配依據，常見者如標準化互相關法（Normalized Cross-Correlation, NCC）、最小二乘匹配法（Least Square Matching, LSM）等；第二類為特徵匹配（Feature-based Matching），即藉由自影像中萃取出之特徵進行匹配，於攝影測量領域中多以點特徵為主，而線特徵則為計算機視覺等領域所偏好，特徵匹配較強鈍（Robust），亦不需嚴格之初始條件；第三類為符號式匹配(Symbolic Matching)或稱關連式（Relational）匹配，近似於特徵匹配，但影像間進行比對之基礎，得自於影像特徵或灰度值特性，稱之為符號化敘述（Symbolic Description）。 1.

(10) 本研究探討之小波方法匹配流程，參考自 Hsieh et al. (1997)，並以洪偉嘉（1999）之流程為對照，茲分別說明如下： Hsieh et al. (1997)發展之流程，乃分別針對兩幅影像進行小波轉換後萃取特徵點，以角度累計直方圖估計旋轉量，再以互相關係數（Cross Correlation Coefficient）決定特徵點間配對關係，以一致性檢查剔除錯誤配對後迭代求解，如圖 1- 1。該流程所採用之轉換模式為四參數相似轉換，不適用於傳統相機攝得之影像，故於本研究中改採八參數透視投影轉換，並加入 TAU 測試除錯機制，如圖 1- 2。洪偉嘉（1999）發展之流程，乃利用興趣運算元對影像進行特徵萃取，續以標準化互相關法、最小二乘匹配法進行不同精度層級之計算，經 TAU 測試除錯後求解並進行影像鑲嵌，如圖 1- 3。. 1-3 論文架構本論文共分為五章，各章主題說明如下：第一章：說明本論文之研究動機與目的；研究方法與內容；並說明本論文之組織與架構。第二章：說明以興趣元進行特徵匹配之型態，並針對萃取方法做一陳述，並利用不同型態之樣本影像，測試兩萃取方法標準成果，以及對雜訊及旋轉量之反應。第三章：針對兩方法中所採用的角度累計直方圖法、標準化互相關係數法、最小二乘影像匹配法、相似度量測、一致性檢查及 TAU 測試等理論，做一說明。第四章：利用不同性質影像進行數值試驗，並對成果進行分析。第五章：對本研究作一結論與建議。. 2.

(11) 讀入兩影像並進行小波特徵萃取. 建立角度累計直方圖，以估計感測影像旋轉量. 以旋轉量估值旋轉感測影像，以互相關係數進行特徵點配對. 進行特徵點對一致性檢查，剔除錯誤特徵點配對. 利用奇異值分解（SVD），迭代求取影像相似轉換參數. 感測影像重新取樣. 鑲嵌感測影像與套合影像. 圖 1- 1 Hsieh et al. (1997)小波影像套合流程. 3.

(12) 讀入兩張影像，進行小波特徵萃取. 建立角度累計直方圖，藉以估計初始旋轉量依旋轉量估值將影像旋轉，以互相關係數進行特徵點配對進行特徵點對一致性檢查，剔除錯誤特徵點配對. 劃分子影像，控制特徵點配對分佈. 進行Tau測試，進一步剔除錯誤點對. 初步求解八參數. 迭代求解. 感測影像重新取樣. 鑲嵌感測影像與套合影像圖 1- 2 本研究採用之小波影像套合流程. 4.

(13) 人工選取兩對共軛點（重疊區的左上與右下）. 以共軛點對為邊界於影像內劃分方格，利用共軛點求解正形轉換係數，作為兩張影像起始的對應關係利用興趣運算元萃取特徵點，以 Spread函式約制候選的準控制點分佈位置。利用轉換參數預估其在感測影像上之近似共軛位置. NCC影像匹配，求取共軛點對應整數像元. 最小二乘影像匹配，求取共軛點的對應次像元 N=N+1 當N>3時： 1.變換以八參數轉換，預估特徵點在感測影像上之對應位置。 2.縮小搜尋窗。. 求解八參數轉換之係數，並用TAU測試剔除不適當的控制點，將所得到的殘差求取RMSE。. N>4 n=N與n=N-1時的八參數中的平移係數差小於1個像元。且n=N與n=N-1的RMSE差小於 n=N-1與n=N-2的RMSE差。. 不是. 是感測影像重新取樣. 鑲嵌感測影像與套合影像. 圖 1- 3 洪偉嘉（1999）興趣運算元影像套合流程. 5.

(14) 第二章特徵萃取特徵（feature）乃指影像中灰度值具有明顯變化的位置，利用影像上存在之特徵進行影像匹配，能有效減少影像中所含的過多資訊，增加匹配結果之可靠度。一般而言，點特徵可藉由興趣運算元獲得，例如：Moravec 運算元[Moravec, 1977]、Foerstner 運算元[Foerstner, 1984]、Lue 運算元[Lue, 1988]、LCF 運算元[Lin and Lee, 1995]及 TDGO（Target Defined Ground Operator）運算元[李良輝， 1990] 等，其特色如下[洪偉嘉， 1999]： 1. 點位分佈： (1). Moravec、Lue、TDGO 三種運算元主要找的為線型交錯點，如屋角。. (2). Foerstner 運算元找到的為顯著點、圖形對稱點、線型交錯點。. (3). LCF 運算元主要找區域灰度最大值的特徵點。. 2. 運算速率：. Lue ≈ TDGO > Moravec > Foerstner > LCF 相較之下 TDGO 具有較佳之運算速率及萃取明銳點位之特性，故洪偉嘉（1999）選用 TDGO 作為特徵點萃取之興趣運算元。除上述以興趣運算元萃取特徵點之方法外，本研究另採用小波（Wavelet）方法[Mallat and Zhong, 1992]進行特徵點萃取，將影像進行小波轉換後，計算其小波轉換模數（Wavelet Transform Modulus）之局部最大值（Local Maxima），藉以求得分佈於邊緣線上之特徵點。以下分別對 TDGO 運算元及小波方法作一概述：. 2-1 TDGO 運算元 1. 基礎運算元（Ground Operator）利用相鄰像元之周圍環境條件，進行位相編碼（Topological Encode），藉以萃取特定形狀之點特徵。以 3x3 影像窗為操作單元，經二元化處理後進 6.

(15) 行影像之位元型式（Bit Pattern）編碼。如圖 2- 1所示，對相鄰像元進行下列處理：. 圖 2- 1 影像之位元型式編碼 IF G 0 − G i > thd. THEN. Bi =2i-1. ELSE. Bi =0. (2- 1). i=1,2,… … 8 其中 G0 ：中心像元之灰度值 Gi ：相鄰像元之灰度值. Thd ：門檻值 Bi ：相鄰像元位置代碼故其位元型式值為： 8. BV = ∑ Bi. (2- 2). i =1. 對不同的環繞條件而言，BV 值是唯一的，且共有 256 種組合，建立特定之點特徵之對照表（Look Up Table, LUT），即可有效且快速萃取指定之點特徵。本文使用之點特徵共有 24 種形式，如圖 2- 2：. 7.

(16) 圖 2- 2 定義特徵點之圖型（Pattern） 2. 興趣點（Interest Point）選取由基礎運算元所萃取的影像點，進一步萃取真正的特徵點，由於經過基礎運算元的篩撿，本階段處理所需處理的點數已大幅減少，計算各點之興趣值（Interest Value, IV）如下： 8. IV = ∑ G0 − Gi. (2- 3). i =1. 3. 壓縮局部非最大值（Supression of Local Non-maximum）對每個興趣點比較其與周圍興趣點之興趣值，以消除在高紋理區域內有連續而過多之興趣點。首先定義壓縮窗之大小，以逐步擴大之方式進行比較，若該窗內存在一興趣值大於中心點之興趣值，則將中心點之興趣值設為 0，並停止比較，如圖 2- 3所示：. 8.

(17) 圖 2- 3 逐步擴大比較窗之非局部最大壓縮. 2-2 小波特徵萃取傳統之頻率域分析技術如傅力葉分析，乃利用多頻率之正弦及餘弦函數疊加，將訊號分解以進行分析。然而，由於其平滑之特性，正弦函數面對變化劇烈之訊號時，擬合效果不甚理想，對於不連續及區域性之變化敏感度較低[Graps, 1995]。小波分析以類似原理運作，但利用一有限且不規則之小波母函數（Mother Wavelet），變化其尺度後對訊號進行擬合。小波分析利用不規則之母函數進行分析，故對不連續之偵測能力較傅力葉分析為佳，影像中之邊緣線即具備灰度值之不連續性。 2-2-1 小波轉換小波轉換利用經擴張（dilated）及平移之小波將訊號進行分解，一小波 ψ 為一平均值為 0 之函數，即： +∞. ∫ψ (t )dt = 0. (2- 4). −∞. 經由尺度變數 s 擴張，平移量 u 後，得： ψ u , s (t ) =. 1 t−u ψ( ) s s. (2- 5). 故函數 f 於時刻 u 及尺度 s 之小波轉換為：. 9.

(18) Wf ( u, s ) =. +∞. ∫. f (t). −∞. 1 * t −u ψ ( ) dt s s. (2- 6). 2-2-2 Canny 邊緣偵測法 Canny 邊緣偵測法 [Canny, 1986] 係利用計算影像之梯度向量模數（Modulus），以偵測影像內灰度值出現尖銳變化之點，梯度向量定義如下式：  ∂f    v ∂x ∇ f =  ∂f     ∂y   . (2- 7). v 函數 f 沿向量 n = ( a, b ) 方向之導數為：. v v ∂f ∂f ∂f a+ b (2- 8) v = ∇f ⋅ n = ∂n ∂x ∂y v v 由上式可知當 n 與 f 共線時，該導數之絕對值為最大值，即代表 ∇f 與函數 f ( x , y ) 出現最劇烈變化之方向平行。若一點 ( x0 , y 0 ) 計算得之梯度模數出現局部. 最大值，則 ( x0 , y 0 ) 視為一邊緣點。 2-2-3 多尺度（Multiscale）邊緣偵測法多尺度邊緣偵測法係針對訊號進行不同尺度之平滑化，利用訊號之一階及二階導數進行尖銳變化點之偵測。若函數 θ ( x) 之積分值等於 1，且函數值於 x 趨近無限大時收斂至 0，則稱其為一平滑化函數（Smoothing Function）。假設 θ ( x) 之一、二階導數俱存在且有定義（defined），則令兩小波函數分列如下： ψ a ( x) =. dθ ( x ) d 2θ ( x) ，ψ b ( x ) = dx dx 2. (2- 9). 另定義尺度如下： f s ( x) =. 1 x f( ) s s. 故 f (x ) 於位置 x、尺度 s 狀況下，對於母函數 ψ a ( x) 之小波轉換為： W sa f ( x) = f *ψ sa ( x) = f * (s. dθ s d )( x) = s ( f *θ s )( x) dx dx 10. (2- 10).

(19) 其中，”*”表摺積（Convolution）同理，對母函數ψ b (x ) 之小波轉換為： W sb f ( x ) = f *ψ sb ( x) = f * (s 2. 2 d 2θ s 2 d )( x ) = s ( f *θ s )( x) dx 2 dx 2. (2- 11). 將上述推導引申至二維影像，則兩小波函數對應至 θ ( x, y) 於 x、y 方向之偏導數如下： ψ1 =. ∂θ ( x , y ) ∂θ ( x, y ) ，ψ 2 = ∂x ∂y. (2- 12). 重新定義尺度如下： f s ( x, y ) =. 1 x y f( , ) 2 s s s. 故 f ( x , y ) 於尺度 s 之小波轉換得以下兩分量： W s1 f ( x, y) = f *ψ 1s ( x, y ) ， W s2 f ( x, y ) = f *ψ s2 ( x, y ). 應用於影像時，令尺度 s 為 2 j ，由一維推導可知： ∂  W f ( x, y )   ∂x ( f * θ 2 j )( x, y )  v   = s  = s∇( f * θ j )( x, y) 2 W j f ( x, y )   ∂  2   ( f * θ 2 j )( x, y)   ∂y  1 2j 2. (2- 13). v 由上式可知， f ( x , y ) 經小波轉換後之兩分量與梯度向量 ∇( f *θ 2 j )( x , y) 間存在比例關係，故當尺度 s = 2 j 時，其梯度向量模數與下式成比例： M 2 j f ( x, y ) =. 2. W21j f ( x, y) + W 22j f ( x, y ). 2. (2- 14). 該梯度向量之方向可由下式求得：. A2 j f ( x, y ) = argument (W21j f ( x, y ) + iW 22j f ( x, y)). (2- 15). 如同 Canny 邊緣偵測法， f *θ 2 j ( x, y) 之 M 2 j f ( x, y) 在 A2 j f ( x, y ) 方向出現局部最大值時，即視為尖銳變化點；於影像中，可視其為邊緣線。然而，上述方法實際應用時仍存在雜訊問題，Xu et al. (1994)利用一邊緣相關（Edge Correlation）法，達到有效區分特徵及抑制雜訊之目的。該方法之理論基礎在於：實際邊緣線會在數個不同尺度之轉換成果中，出現較強訊號；而雜訊 11.

(20) 之反應卻會隨著尺度增加而迅速衰減。透過多尺度之小波轉換乘積，可將邊緣之訊號放大抑制雜訊，邊緣相關計算式如下[Hsieh et al., 1997]： n −1. Rn ( j , x, y ) = ∏ M 2 j+i f ( x, y). (2- 16). i =0. 其中 n 為計算中所涵蓋尺度數目 j 為計算之初始尺度值圖 2- 4為邊界相關計算之範例，由範例中可見，真實邊緣線經由乘積凸顯而出，而雜訊部分受到抑制。. (a) 原始訊號 f (x ). (b) W21 f ( x). (c) W22 f ( x). (d) W23 f ( x). (e) R2 (1, x ) = W21 f ( x ) ⋅ W2 2 f ( x) 圖 2- 4 一維訊號邊緣相關計算成果[Hsieh et al., 1997] 2-2-4 特徵點選取影像經多尺度邊緣偵測獲得影像之邊緣資訊後，需透過一機制自邊緣資訊 12.

(21) 中匹配選取所需之特徵點，任一特徵點 P(x,y)需符合以下條件[Hsieh et al., 1997]： (1). M 21 f ( x, y ) > Thd. (2). R 2 (1, x, y) > M 21 f ( x , y ). (3). M 21 f ( x, y) = max {M 21 f ( x′, y′)} ( x ′ , y′ )∈N P. 其中 Thd 為選取邊緣線之門檻值，經本研究測試後建議值為 25.0。為確保 R2 與 M 21 為同一數量級，故將 R2 進行正規化得 R2 ，定義如下：. R2 ( j , x, y ) = R2 ( j , x, y). ∑. x, y. R2 ( j , x , y ). 2. ∑. x, y. M 2i f ( x , y ). 2. N p 代表 P 點之鄰近區域，此處指以 P 點為中心，大小為 n*n 之窗型區域。. 13.

(22) 2-3 特徵萃取成果比較分析本研究使用 ISPRS 提供之航空立體像對 suburb_l.tif、suburb_r.tif 以及 SPOT 衛星紅光波段（XS2）之翡翠水庫影像 Ired.tif 作為實驗樣本，分別對 TDGO 運算元及小波方法兩者進行特徵萃取測試。其後選取 suburb_l.tif，經由適當角度之旋轉及加入隨機雜訊，俾便瞭解其對於特徵匹配過程中常見變數：旋轉量及雜訊之反應。成果表中初選點位欄位，代表符合 TDGO 運算元指定灰階分佈型態之點位數目，而結果點位代表經壓縮局部非最大值處理後之特徵點數目。 2-3-1 標準測試樣本影像 suburb_l.tif、suburb_r.tif 大小為 900*900，如圖 2- 5。分別以 TDGO 運算元及小波方法萃取成果如表 2- 1。表 2- 1 ISPRS 立體像對特徵萃取成果記錄 TDGO 運算元 suburb_l.tif. 小波方法. suburb_r.tif. suburb_l.tif. suburb_r.tif. 搜尋窗邊長. 3. *. Thd. 55. 110. 壓縮窗邊長. 5. 5. 運算時間（秒）. 24. 10. 初選點位. 879. 432. 結果點位. 192. 121. * 198. 146. 前文 2-2-4 中提到小波方法門檻值約為 25.0，但此處為避免特徵點數目過多不易判讀，故將門檻值提高為 110。. 14.

(23) (a)suburb_l.tif. (b)suburb_r.tif 圖 2- 5 ISPRS 立體像對. 圖 2- 5之影像經 TDGO 運算元萃取，選取出 IV 值較大的 100 點，標註於原影像後，即為圖 2- 6。. (a)suburb_l.tif. (b)suburb_r.tif. 圖 2- 6 ISPRS 立體像對經 TDGO 運算元萃取之特徵點位圖標準樣本影像經小波方法進行特徵萃取，選出小波轉換梯度向量模數較大之 100 點並標註於原影像後，可得圖 2- 7。. 15.

(24) (a)suburb_l.tif. (b)suburb_r.tif. 圖 2- 7 ISPRS 立體像對經小波方法萃取之特徵點位圖將試驗樣本改為翡翠水庫影像，大小為 512*512，以兩方法分別進行萃取，成果如表 2- 2。表 2- 2 翡翠水庫衛星影像特徵萃取成果記錄 TDGO 運算元. 小波方法. 搜尋窗邊長. 3. *. Thd. 13. 85. 壓縮窗邊長. 5. 5. 運算時間（秒）. 4. 2. 初選點位. 632. *. 結果點位. 185. 198. 將表 2- 2 TDGO 運算元萃取成果選出 IV 值較高之前 100 點，及小波方法萃取成果中小波轉換梯度向量模數較大之前 100 點，分別標註於原影像上，成果如圖 28。. 16.

(25) (b)TDGO 法萃取成果. (a)原始影像. (c)小波方法萃取成果圖 2- 8 翡翠水庫衛星影像特徵萃取成果圖. 17.

(26) 2-3-2 旋轉量測試以 suburb_l.tif 分別順鐘向旋轉 30 度及 45 度，與原影像一併經適當裁切後，成為大小 608*608 之樣本影像，分別以 TDGO 運算元及小波方法進行特徵萃取後，成果如表 2- 3。表 2- 3含旋轉量影像特徵萃取成果記錄 TDGO 運算元. 小波方法. 旋轉量. 0. 30. 45. 0. 30. 45. 搜尋窗邊長. 3. 3. 3. *. *. *. Thd. 40. 40. 40. 90. 90. 90. 壓縮窗邊長. 5. 5. 5. 5. 5. 5. 運算時間（秒）. 6. 6. 6. 3. 3. 3. 初選點位. 1425. 294. 325. *. *. *. 結果點位. 224. 132. 138. 152. 189. 169. 18.

(27) (a) 原始影像. (b) 旋轉 30 度. (c) 旋轉 45 度圖 2- 9 含旋轉量之樣本影像圖 2- 9三幅影像經 TDGO 運算元進行特徵萃取，選取 IV 值較高之前 100 點標示於原影像，經人工辨識後，將原影像與旋轉後影像萃取成果相差在 1 像元內之點位視為重疊，以紅色十字標表示之，如圖 2- 10。為使比較之區域不受旋轉影響，僅以像幅中心為圓心，半徑 450 像元圓形區域為比較範圍（測試影像均裁切自該圓形區域內）。 19.

(28) (b)旋轉 30 度. (a)原始影像. (c)旋轉 45 度圖 2- 10 含旋轉量影像經 TDGO 運算元萃取之特徵點位圖觀察圖 2- 10，旋轉量為 30 度的狀況下，重疊點數為 14 點，及至旋轉量為 45 度時，重疊點數上升至 20 點。相同影像經小波方法萃取，選取小波轉換梯度向量模數值較大之前 100 點標註於原影像，經人工辨識後，將原影像與旋轉後影像萃取成果相差在 1 像元內之點位視為重疊，以紅色十字標表示之，得成果如圖 2- 11。. 20.

(29) (b)旋轉 30 度. (a)原始影像. (c)旋轉 45 度圖 2- 11 含旋轉量影像經小波方法萃取之特徵點位圖觀察圖 2- 11，旋轉量為 30 度的狀況下，重疊點數為 14 點，及至旋轉量為 45 度時，重疊點數維持在 15 點。就小波方法而言，若將梯度向量模數值大於門檻值 Thd 之點，將其梯度向量模數值對應一 LUT，可得邊緣強度影像如圖 2- 12。. 21.

(30) (b)旋轉 30 度. (a)原始影像. (c)旋轉 45 度圖 2- 12 含旋轉量影像之邊緣強度影像觀察圖 2- 12可知，高強度之邊緣線如建物邊緣等，於三幅影像中仍具一致性，旋轉量影響並不顯著。. 22.

(31) 2-3-3 雜訊測試本測試仍利用 suburb_l.tif 作為實驗樣本，大小為 900*900，於樣本中混入隨機雜訊，藉以測試兩方法雜訊之敏感度。本研究所選用的雜訊型態為 Gaussian 雜訊，平均值為 0，變方分別為 0.01 及 0.05（8 位元灰階值經正規化至 0 與 1 間），雜訊影像萃取測試成果如表 2- 4。表 2- 4 含雜訊影像特徵萃取成果紀錄 TDGO 運算元. 小波方法. 變方. *. 0.01. 0.05. *. 0.01. 0.05. 搜尋窗邊長. 3. 3. 3. *. *. *. Thd. 55. 117. 225. 110. 110. 125. 壓縮窗邊長. 5. 5. 5. 5. 5. 5. 運算時間(秒). 12. 12. 12. 5. 5. 5. 初選點位. 879. 253. 216. *. *. *. 結果點位. 192. 198. 197. 198. 206. 198. 試驗樣本影像如圖 2- 13。. 23.

(32) (b)變方 0.01 雜訊影像. (a)原始影像. (c)變方 0.05 雜訊影像圖 2- 13 含 Gaussian 雜訊之樣本影像圖 2- 13經 TDGO 運算元進行特徵萃取，為便於判讀，故將特徵點標註於不含雜訊之原始影像上，經人工辨識，將原影像與雜訊後影像萃取成果相差在 1 像元內之點位視為重疊，以紅色十字標表示之，如圖 2- 14。. 24.

(33) (b)變方 0.01 雜訊影像. (a)原始影像. (c)變方 0.05 雜訊影像圖 2- 14 含雜訊影像經 TDGO 運算元萃取之特徵點位圖自圖 2- 14可發現，變方為 0.01 的雜訊影像中，萃取成果重疊者有 8 個點，但在變方為 0.05 的狀況下，受到雜訊干擾的狀況較嚴重，重疊點數為 0 點。而相同樣本影像經小波方法萃取後，為便於判讀，故將特徵點標註於不含雜訊之原始影像上，經人工辨識後，將原影像與旋轉後影像萃取成果相差在 1 像元內之點位視為重疊，以紅色十字標表示之，如圖 2- 15。. 25.

(34) (b)變方 0.01 雜訊影像. (a)原始影像. (c)變方 0.05 雜訊影像圖 2- 15含雜訊影像經小波方法萃取之特徵點位圖觀察圖 2- 15發現，在變方為 0.01 之雜訊影像下，點位重疊者有 41 點，至變方提高為 0.05 狀況下，重疊點位下降至 28 點。同前，利用小波轉換所獲得之邊緣強度影像如圖 2- 16。. 26.

(35) (b)變方 0.01 雜訊影像. (a)原始影像. (c)變方 0.05 雜訊影像圖 2- 16 含雜訊影像之邊緣強度影像由圖 2- 16可知，隨著雜訊數量的增加，梯度向量模數計算成果內雜訊干擾亦較顯著，進行區域性觀察後，可知高強度部分於三幅影像中仍保有一致性，若採用梯度模數值為選取特徵點依據，可有效降低雜訊影響。 2-3-4 成果歸納與分析由前數節列表及成果影像可以歸納得以下數點： (1). 就現有選點方式而言，對於影像旋轉後特徵點位的不變性（rotation 27.

(36) invariance），小波方法及 TDGO 運算元並無明顯差異。 (2). 就雜訊影響程度而言，TDGO 運算元所依靠之點特徵型態受到隨機雜訊干擾，故在雜訊存在狀況下，萃取效果較差，在變方為 0.05 狀況下，重疊點數甚至為 0。. (3). 小波方法由於以小波轉換梯度向量模數值作為選點依據之一，而不考慮點特徵之灰階組成型態，參考圖 2- 16邊緣強度影像，即使雜訊存在且嚴重，經由適當的門檻值修正，仍可萃取得強度較高之特徵點。. (4). 參考圖 2- 12之邊緣強度影像，可得類似於前述(3)中之狀況，於旋轉量存在之狀況下，強度較高之邊緣仍存在。小波轉換之成果原本即為一連串之邊緣線，透過 2-2-4 節中所陳述之選點機制選取特徵點，但壓縮窗移動時為依影像邊界逐次移動，故在旋轉後影像內所有特徵均移位的狀況下，可預期點位重複之機會必然下降。. (5). 觀察圖 2- 8之翡翠水庫萃取成果，兩方法均萃取出明暗變化較大的點位，如水陸交界及山脊稜線，但小波方法由於前述(4)之狀況限制，較少出現連續而密集之點位。. 28.

(37) 第三章影像匹配影像匹配依據匹配方式可分為：區域匹配、特徵匹配及符號式匹配等，本研究採用之小波套合流程及對照之洪偉嘉（1999）流程，於影像匹配部分皆採區域匹配。本章即依套合流程中，影像匹配之前置處理（3-1 節）、匹配原理（32 至 3-4 節）以及後續之除錯（3-5 至 3-6 節），說明其原理。. 3-1 角度累計直方圖（Angle Histogram）影像匹配前之初始方位判定，可利用自影像萃取之線特徵，計算線特徵之走向（Orientation），藉比較線特徵走向差值製成一角度累計直方圖，該直方圖即可提供兩影像間概略旋轉角度資訊[Heipke, 1997]。本研究以點特徵作為匹配依據，必須對特徵點進行線性擬合計算，取得各特徵點之走向資訊，以製作角度累計直方圖。利用一大小為 ( 2M + 1) 2 的運算窗，對邊緣點 P(x,y)而言，位在運算窗內可決定 P 點走向之點集合 ∑ pi ( xi , yi ) 需滿足兩項約制條件[Hsieh et al., 1997]： 1. 連接性：若 q k = ( ~xk , ~ y k ) 為特徵點 pi ( xi , yi ) 之鄰近點，若 ~ x k − xi > 2 或 ~ y k − yi > 2 ，則不符合連接性約制，即 q k ∉ ∑ pi ( xi , yi ) 。 2. 方向一致性：同上，若 pi −1 pi 與 p i q k 間夾角大於 p i−1 pi 與 p i pi +1，則不符合方向一致性，即 q k ∉ ∑ pi ( xi , yi ) 。透過上述約制條件，可得近似擬合於直線 yi = mxi + b 之點集合. {pi. = ( xi , yi )}i=1, 2,...., N ，並可透過以下函數評估其擬合程度：.  y − b − mxi   ξ = ∑  i σ i2 i =1   N. 2. (3- 1). σ i 為點 ( xi , yi ) 之權值，當 ξ 出現最小值時為最佳擬合解，故可由 ξ 最小化得(3- 2) 29.

(38) 式： N N ∂ξ ( y − b − mxi ) xi ∂ξ y − b − mxi = −2∑ i =0， = −2∑ i =0 2 ∂m σi ∂b σ i2 i =1 i =1. (3- 2). 令 σ i ＝1，求解上述二式可得 m 如(3- 3)式： N. m=. N X Y − ∑ xi y i i =1 N. (3- 3). NX X −∑x. 2 i. i =1. 其中 X=. 1 N. N. ∑x. i. ，Y =. i =1. 1 N ∑ yi N i=1. 求解 m 即可得過 P 點直線之走向。. {. }. {. }. 若 FPf1 = pi = ( p ix , p iy ) i =1, 2,..., N ，FPf2 = q j = ( qxj , q yj ) f1. j =1 , 2 ,..., N f 2. 各代表自兩幅影像. 中萃取出之特徵點集合，令 A(u)代表特徵點之走向，則自 FPf1 取一點 pi 、 FPf2 中取一點 q j ，則兩點之走向差值為： θ ji = A( q j ) − A( pi ). 將所有點位所計算出之角度累計成為一直方圖，則數量最大之角度即為兩影像間旋轉角的估計值 θ 。. 3-2 標準化互相關法（Normalized Cross-Correlation, NCC）標準化互相關法主要計算灰度向量之夾角，將每個運算窗內灰度之標準偏差視為一向量，故兩運算窗間之關係可利用兩向量之夾角表示，當兩向量最接近時，相關係數最高（接近 1），但此方法對於影像間旋轉與比例尺之變化相當敏感，亦即影像間出現上述兩狀況時，有可能匹配失敗。NCC 計算公式為： n. ρ TM. σ TM = = σ T ⋅σ M. ∑(g i =1. Ti. n. ∑(g i =1. Ti. − u T )( g M i − u M ). − uT ). (3- 4). n. ∑(g i =1. Mi. − uM ). 式中： g Ti , g M i ：分別為目標窗（Target Window）與遮罩窗（Mask Windows）內 30.

(39) 各像元之灰度值 uTi , u Mi ：分別為目標窗與遮罩窗內所有像元之平均值. ρ TM ：互相關係數於本研究中，利用目標窗與對應於遮罩窗，可求得一相關係數，將所有相關係數做一排序並選取最大值，當此值大於門檻值（例如：0.75），則視為匹配成功，記錄該點座標，作為最小二乘匹配之起始值。. 3-3 最小二乘匹配法(Least-Squares Matching, LSM) 利用最小二乘影像匹配技術，可求取影像座標至次像元精度，Rosenholm (1987)將最小二乘匹配分為二種模式，其一為僅考慮幾何平移參數者，即 GT ( X , Y ) = G S ( x + x 0 , y + y0 ) + n( x, y ). (3- 5). 而另一種為同時考慮輻射及幾何變形參數者，即 GT ( X , Y ) = h1 + h2 ⋅G S ( x, y; A) + n( x0 , y 0 ). (3- 6). 其中 GT , GS 分別為目標窗及搜尋窗之灰度值函數. h1 , h2 則為灰度值之平移與尺度參數 x 0 , y 0 為幾何平移參數 A 為仿射轉換設計矩陣本研究所使用者為考慮輻射及幾何變形之最小二乘相關法 (LSM2)。 3-3-1 LSM2 平差模式若以 GT ( X i , Yi ) 代表目標窗像元之灰度函數， GS ( x i , y i ) 代表搜尋窗像元之灰度函數，於(1)搜尋窗與目標窗不大，(2)影像的幾何變形不甚嚴重，二個基本前提下，目標窗與搜尋窗間之轉換模式可用六參數仿射轉換(Affine Transform)來描述，另外由於攝影時光線及陰影的影響，同時將灰度值之平均尺度參數 h1 , h2 列入考慮，故整個方程式可寫成[Ackermann, 1984]：. 31.

(40) T.   a GT ( X i , Yi ) = h1 + h2 ⋅ GS   1  a 4 . a2 a5.  1  a3    x  = h1 + h2 ⋅ GS ( x i , y i ; A) a 6   i    y i  . (3- 7). 其中：. a A=  1 a 4. a2 a5. a3  a6 . 而 GS ( x i , y i ; A) 表示由灰度值 GT ( X i , Yi ) 之 ( X i , Yi ) 座標經 A矩陣轉換到搜尋窗中所在位置之灰度值。LSM2 即解算 A 矩陣中之六元素與 h1 , h2 共八個參數，但(3-6)式為非線性，故必須將其線性化，得(3-8)式後再進行迭代計算： 6.  0 ∂G S  GT ( X i , Yi ) = h + h G S ( xi , y i ; A ) + ∆h1 + G S ( xi , y i ; A ) ⋅ ∆h2 + ∑  h2 ∆a n  ∂a n   0 1. 0 2. 0. 0. n =1. (3- 8) 其中 h10 、 h 20 、 A0 分別為 h1 、h2 、 A 之近似值；∆h1、∆h2 、∆A 則為 h1 、h2 、 A 的變化量，即平差時的未知數； ∆a n 的係數則為搜尋窗灰度函數對仿射轉換中. 的六參數偏微分而來，參考(3-6)式可得： • • • ∂GS ∂GS ∂x i ∂GS ∂y i = + = G Sx ⋅ 1 + G Sy ⋅ 0 = G Sx ∂a1 ∂xi ∂a1 ∂y i ∂a1 • • • ∂GS ∂G S ∂x i ∂GS ∂y i = + = G Sx ⋅ X i + G S y ⋅ 0 = X i G Sx ∂a 2 ∂x i ∂a2 ∂y i ∂a 2 • • • ∂GS ∂GS ∂xi ∂GS ∂y i = + = G Sx ⋅ Y j + G S y ⋅ 0 = Y j ⋅ G Sx ∂a3 ∂xi ∂a 3 ∂y i ∂a3 • • • ∂GS ∂GS ∂xi ∂G S ∂y i = + = G Sx ⋅ 0 + G S y ⋅1 = G Sy ∂a 4 ∂xi ∂a 4 ∂y i ∂a4 • • • ∂GS ∂GS ∂xi ∂GS ∂y i = + = G Sx ⋅ 0 + G S y ⋅ X i = X i ⋅ G Sy ∂a5 ∂xi ∂a 5 ∂yi ∂a 5 • • • ∂GS ∂G S ∂xi ∂GS ∂y i = + = G Sx ⋅ 0 + G S y ⋅ Yi = Yi ⋅ G Sy ∂a6 ∂x i ∂a6 ∂y i ∂a6 •. (3- 9). •. 其中 ( xi , y i ) 為搜尋窗內像元之坐標， G Sx 及 G Sy 分別為灰度函數 G S ( x i , y j ) 在 X 、 Y 方向之一次導數，亦即其於 X 、 Y 方向之梯度。實際上，灰度函數是 32.

(41) 無從得知的，影像灰度值乃是灰度函數加上雜訊而來: g T ( X i , Yi ) = GT ( X i , Yi ) + nT ( X i , Yi ) g S ( x i , y i ) = GS ( x i , y i ) + nS ( xi , yi ). (3- 10). 令各未知參數之初始值為：. h10   0 h2  a10   0 a 2  = a 0   30  a 4  a 0   5 a 60 . 0 1   0   0 1   0 0   1. (3- 11). 將（3-8）、（3-9）、(3-10)式代入(3-6)式，得:. {g. T. ( X i , Y j ) − g S ( xi , y j )} + {nT ( X i , Y j ) + n S ( x i , y j )} =. ∆h1 + G S ( xi , y j ) ⋅ ∆h2 + •. •. •. •. •. •. G Sx ( X i , Y j ) ⋅ ∆a1 + X i G Sx ( X i , Y j ) ⋅ ∆a 2 + Y j G Sx ( X i , Y j ) ⋅ ∆a3 + G Sy ( X i , Y j ) ⋅ ∆a4 + X i G Sy ( X i , Y j ) ⋅ ∆a5 + Y j G Sy ( X i , Y j ) ⋅ ∆a6. (3- 12). 令 l ( X i , Y j ) = g T ( X i ,Y j) + g S ( xi , y j ) 為觀測量 V ( X i , Y j ) = nT ( X i , Y j ) + nS ( x i , y j ) 為改正數. 則(3-11)式可改寫為間接平差之觀測方程式 V = AX − l.  V ( X i , Y j ) = 1 GS . •. •. •. •. •. GSx. GSy. X i GSx. Y j GSx. X i GSy.  ∆h1   ∆h   2  ∆a1    •  ∆a 2  Yi GSx  − l ( X i ,Yj )   ∆a 3    ∆a 4   ∆a   5 ∆a 6  (3- 13). 33.

(42) 假設觀測值是等權不相關的，則權矩陣 P = I ，故法方程式係數矩陣 N = AT A ， U = AT l ，未知數矩陣 X = N −1U ，即可求出各項轉換參數。迭代計. 算使得未知數之增量逐漸收斂至臨界值以下，則此時迴圈結束，求得仿射轉換六參數 a1 ~ a 6 ，如此即可將目標窗像元轉換至搜尋窗中。 3-3-2 LSM2 運算流程進行最小二乘相關匹配時，需注意以下項目： 1. 目標窗與搜尋窗的大小: 一般而言，目標窗大小與影像品質有關，當原始影像反差甚大時，目標窗的尺寸不須太大，但若影像品質不佳，則有必要加大目標窗尺寸以納入更多鄰近像元資訊。又，搜尋窗大小，則可根據目標窗大小加以定義，王聖鐸（1997）建議搜尋窗邊長為目標窗邊長之 1.5 至 2.5 倍。 •. •. 2. 計算灰度梯度 G Sx ( x i , y i ) 與 G Sy ( xi , yi ) ： •. •. 平差時，因設計矩陣 A 中之 G Sx 與 G Sy 無法直接求得，故必須以灰度觀測值 g s 之一次差分代替[陳典熙，1988]： •. •. G Sx ≈ g Sx = =. {g S (x + 1, y ) − g S (x, y )}+ {g S ( x, y) − g S (x − 1, y )} 2. g S ( x + 1, y ) − g 2. S. ( x − 1, y ). 同理 •. •. G Sy ≈ g Sy =. g S ( x, y + 1) − g S ( x , y − 1) 2. (3- 14). 3. 重新取樣求得轉換參數的變化量後，再加上其初始值便可求出轉換參數，利用此轉換參數可將目標窗之中心像元座標轉換至搜尋窗中對應之座標，但因轉換後之座標往往恰非整數，即會無法找到對應之灰度觀測值，如圖 3- 1，因 34.

(43) 此必須經重新取樣以內插得對應之灰度值。內插方法種類甚多，本研究所使用者為雙線性內插（Bilinear Interpolation）： g S ( x i + α , y j + β ) = (1 − α )(1 − β ) g S ( x i +1 , y j +1 ) + (1 − α ) ⋅ β ⋅ g S ( x i +1 , y j ) +. α ⋅ (1 − β ) ⋅ g S ( x i , y j+1 ) + α ⋅ β ⋅ g S ( xi , y j ). (3- 15). 圖 3- 1 雙線性內插重新取樣 4. 迭代結束條件 LSM 迭代過程結束與否，可由以下數項判斷： a. 目標窗內之代表像元位置以轉換參數進行座標轉換後，是否超出搜尋窗範圍？ b. 轉換參數之變化量是否小於臨界值？ c. 迭代之次數？若 a 條件成立，則疊代結果為發散，迴圈結束，若 a 不成立則考慮 b 條件。一般於 b 條件可考慮轉換參數中的平移參數 a1 ， a 4 的變化量，若 a1 與 a 4 的變化量小於所設的臨界值，則迴圈結束，匹配完成。c 條件主要目的在 35.

(44) 當 a, b 條件都不成立時，為避免陷入無窮迴圈，故設迭代次數上限，當迭代超過此次數之後，則迴圈結束。整體而言，唯有 b 條件成立，方代表最小二乘匹配成功，否則皆為失敗。. 3-4 相似度量測一般而言，影像 f 1( x, y ) 及 f 2 ( x , y) 內兩點 ( m, n) 及 (u , v ) 間之匹配關係，可藉由相似度量測獲得。常用之相似度量測機制甚多，計算最大互相關係數之方式由於不受線性輻射變形（Linear Radiometric Distortion）影響，故本研究採用之。定義如下[Zheng and Chellappa, 1993]： C f1 f2 ( m, n; u, v) =. i , j =ω. m 1 ∑ [ f1 (m + i, n + j) − µ1] ⋅ [ f2 (u + i, v + j) − µ 2 ] σ 1σ 2 ( 2ω m + 1) 2 i , j = −ω m (3- 16). 其中 µ i 及 σ i 分別代表兩影像之局部平均值及變方 ( 2ω m + 1) 2 為匹配計算窗之面積. (3-16)式對於旋轉量之影響非常敏感，將旋轉量影響加入(3-16)式後，得下式 [Hsieh et al., 1997]：. [. i , j =ω m 1 C f1 f2 ( m, n; u, v ) = ∑ [ f (m + i, n + j) − µ1 ]⋅ f 2 (u + iˆ, v + ˆj) − µ2 σ 1σ 2 ( 2ω m + 1) 2 iˆ, jˆ= −ω m 1 (3- 17) ˆˆ. ]. 其中 i = iˆcosθ − ˆj sin θ j = iˆsin θ + ˆ j. 若可事先求得(3-17)式中之旋轉量 θ ，可提升點對匹配關係計算之正確性。. 3-5 配對一致性檢查特徵點配對僅由相似度量測方式決定，可能出現錯誤配對之狀況，藉由配對一致性檢查，可有效剔除錯誤配對，Hsieh et al. [1997]提出以下檢查機制： 36.

(45) 今有兩影像 f 1( x, y ) 及 f 2 ( x , y) ，兩影像間概略旋轉角已知，將 f 1( x, y ) 經該旋轉角修正後得 fˆ1( x, y ) 。令 MP = {p i ⇔ q i }i=1, 2,..., N 代表匹配點對集合， p i = ( x ip , y ip ) 為影像 fˆ1( x, y ) 中一特徵點， qi = ( xiq , yiq ) 代表影像 f 2 ( x , y) 中一特徵. 點，若點對集合 MP 中所有配對均正確，則下式應成立： pi = sqi + T. (3- 18). 其中 s 為尺度關係 T 為平移量假設 {pi ⇔ qi }及 {p j ⇔ q j } 各為點對集合 MP 中兩正確之匹配點對，則 fˆ1( x, y ) 與 f 2 ( x , y) 間之尺度關係 s 可由(3-19)式求得：. s = d2 d. (3- 19) 1. 其中 d1 = (xip − x jp ) 2 + (yip − y pj )2 d 2 = (xiq − x qj )2 + (yiq − y qj )2 待求得 s 後， pi 與 qi 間之平移量 Ti 可由上式求得，同理可求得 T j 。若 {pi ⇔ qi }與 {p j ⇔ q j }均為正確點對，則 Ti 與 T j 間差值應甚小。故藉由檢查 Ti 與 T j 之差值大小，吾人可得知 {pi ⇔ qi }與 {p j ⇔ q j }是否一致。假設 S (i ) 代表點對 {pi ⇔ qi }與其他點對間達成一致之次數，則： IF. Ti − T j < Thd. THEN. S (i ) = S (i ) + 1. 俟其與 MP 集合內所有點對進行一致性檢查後，得累計值 S (i ) 。當 S (i ) 大於門檻值時，將點對 {pi ⇔ qi }視為正確配對；反之，當 S (i ) 小於門檻值時，視為錯誤配 37.

(46) 對並剔除。. 3-6 TAU 測試在匹配時為剔除粗差，本研究進一步在轉換參數求取之過程中，配合幾何轉換模式，以 TAU 測試法剔除不適當之觀測量。方式如下： 1. 計算得某一觀測量之殘差 V 及其後驗中誤差 σ v 。 V=AX+L=AX-(-L). (3- 20). Xˆ = ( AT PA) −1 AT P (− L). (3- 21). 將(3-19)代入(3-18)得：. {. }. (3- 22). {. }. (3- 23). V = I − A( AT PA) −1 AT P L 定義 M = I − A( AT PA) −1 AT P 則. {. }. ˆ = P −1 − A( AT PA) −1 AT σˆ2 ∑ v 0. (3- 24). 其中 σˆ02 =. V T PV n−t. (n：觀測數，t：未知數個數， r = n − t. r：多餘觀測數) (3- 25). 假設 P=STS 則 M 可轉換為 M = I − A ( A T A ) −1 A T. (3- 26). 將式(3-23)代入式(3-24)中，可得： σˆV2i = (. mii V T PV )( ) pi r. (3- 27). 但因(3-25)式中之 m 計算時較為複雜，故 Pope [1975]提出一 σˆV2i 近似值計算方式： σˆV2i =. n − u σˆ0 n Pi. (3- 28). 2. 已知信心水平 (1 − α ) ，計算 Tau 值 τ α ：. 38.

(47)  rF α  1− ; r −1   n τα =   r −1 + F α  1− ;1, r −1  n  . 1/ 2. (3- 29). r：多餘觀測數， α ：信心區間，n：觀測數 3. 若. V > τ α ，則該觀測量為一錯誤。 σˆV. 39.

(48) 第四章數值測試本章利用不同型態影像進行套合流程之數值測試，三組測試影像各為航照立體像對、近景相片以及近景數值影像，經洪偉嘉（1999）之流程及本研究採用之小波流程進行運算，詳列其參數值及鑲嵌成果，以進行分析比較。. 4-1 航照立體像對. (a)1831.tif (1582*1873). (b)1841.tif (1582*1873). 圖 4- 1 航照立體像對第一組測試影像為林務局農林航空測量所於民國 86 年 5 月 15 日所拍攝，航線 86p26 編號 183 及 184 之航空照片，攝影航高 6500 呎，焦距 15mm，比例尺約 1：13000，以解析度 200dpi 掃瞄而成，涵蓋地區約在新竹市交通大學附近，如圖 4- 1。. 40.

(49) 指定左上角點對（904,866；910,227）及右下角點對（1417,1507；1435,863），門檻值 70，壓縮窗邊長 5，搜尋窗邊長 51，目標窗邊長 13，圖 4- 1以洪偉嘉（1999）之流程經四次迭代，得 35 組控制點對，如圖 4- 2。. (a) 183.tif. (b) 184.tif. 圖 4- 2 航空立體像對經洪偉嘉（1999）之流程獲得控制點對成果圖其殘差分佈如圖 4- 3。 RMSE=0.56724 2 1.5 1 0.5 Vy. 0 -2. -1.5. -1. -0.5. 0. 0.5. 1. 1.5. -0.5 -1 -1.5 -2 Vx. 圖 4- 3 航照立體像對經洪偉嘉（1999）之流程殘差分佈圖 41. 2.

(50) 求解得之八參數如表 4- 1。表 4- 1 航照立體像對經洪偉嘉（1999）之流程求解八參數係數表係數. a1. b1. c1. a2. 值. 0.9876. 0.0241. -1.1323. -0.0183. b2. 1.0053 -627.0061. 透過求解得參數將影像重新取樣並進行鑲嵌後，得圖 4- 4。. 42. c2. a3. b3. -1E-5. 1e-5.

(51) 圖 4- 4 航照立體像對掃瞄影像經洪偉嘉（1999）之流程鑲嵌合成圖. 43.

(52) 相同像對經本研究採用之小波流程，設定梯度向量模數門檻值 120.0，互相關係數門檻值 0.75，一致性檢查門檻值 10 進行自動化匹配後，各獲得 605 及 761 個特徵點，分佈如圖 4- 5。. (a) 183.tif. (b) 184.tif 圖 4- 5 航照立體像對小波特徵萃取成果圖. 特徵點經互相關配對後，得 124 組特徵點對，分佈如圖 4- 6。. (a)183.tif 圖 4- 6. (b) 184.tif 航照立體像對特徵點初步配對成果圖 44.

(53) 經配對一致性檢查，初步剔除錯誤配對後，得 14 組特徵點配對，如圖 4- 7。. (a)183.tif (b) 184.tif 圖 4- 7 航照立體像對特徵點配對初步除錯成果圖經 TAU 測試後，進一步剔除錯誤配對，得 13 組特徵點配對，如圖 4- 8。. (a)183.tif. (b) 184.tif. 圖 4- 8 航照立體像對特徵點配對 TAU 測試成果圖圖 4- 8中所示 13 組控制點對，其殘差分佈如圖 4- 9。 45.

(54) RMSE=1.063151 3. 2. 1. Vy. 0 -3. -2. -1. 0. 1. 2. 3. -1. -2. -3 Vx. 圖 4- 9 航照立體像對經小波流程殘差分佈圖求解得八參數係數如表 4- 2。表 4- 2 航照立體像對經小波流程求解八參數係數表係數. a1. b1. c1. a2. 值. 0.9987. 0.0143. -2.3743. -0.0120. b2. c2. 0.9938 -621.9988. 經由表 4- 2所示參數，將影像重新取樣並鑲嵌後，可得圖 4- 10。. 46. a3. b3. -3E-6. 2E-6.

(55) 圖 4- 10 航照立體像對掃瞄影像經小波流程鑲嵌合成圖. 47.

(56) 值得注意的是，洪偉嘉（1999）之流程強調控制特徵點分佈，但鑲嵌成果中仍出現精度不均之現象，對照同區域小波流程成果放大圖，如圖 4- 11。推測其成因應為分佈於該區域特徵點對數目不足，故造成精度不均之結果，可由圖 42得證。. (a) 洪偉嘉（1999）之流程成果區塊放大. (b)小波流程成果區塊放大圖 4- 11 航照立體像對鑲嵌成果局部放大比較圖. 48.

(57) 4-2 近景相片. (a)ll.tif (992*647). (b)rr.tif (992*647). 圖 4- 12 野外地質調查拍攝之近景相片第二組測試影像為民國 85 年 2 月 13 日進行野外地質調查，以 135 相機拍攝之近景相片，以 200dpi 解析度掃瞄而成，如圖 4- 12。以洪偉嘉（1999）之流程，指定左上角點對（387,244；35,190）及右下角點對（925,647；547,612），門檻值 28，壓縮窗邊長 5，搜尋窗邊長 71，目標窗邊長 13，經過五次迭代運算後，共得 37 組控制點對，如圖 4- 13。. (a) ll.tif. (b) rr.tif. 圖 4- 13 近景相片經洪偉嘉（1999）之流程獲得控制點對成果圖其殘差分佈如圖 4- 14。. 49.

(58) RMSE= 0.655785 2 1.5 1 0.5 Vy. 0 -2. -1.5. -1. -0.5. 0. 0.5. 1. 1.5. 2. -0.5 -1 -1.5 -2 Vx. 圖 4- 14地質調查近景相片經洪偉嘉（1999）之流程殘差分佈圖解得八參數係數如。表 4- 3。表 4- 3 地質調查近景相片經洪偉嘉（1999）之流程求解八參數係數表係數. a1. 值. 1.2682. b1. c1. a2. -0.0831 -424.7561 0.1619. b2. c2. 1.181314 -137.3658 0.0003. 以。表 4- 3參數將影像重新取樣並進行鑲嵌，得圖 4- 15。. 圖 4- 15 地質調查掃瞄影像經洪偉嘉（1999）之流程鑲嵌合成圖 50. a3. b3 5E-6.

(59) 相同像對經本研究採用之小波流程，設定梯度向量模數門檻值 25.0，互相關係數門檻值 0.8，一致性檢查門檻值 10 進行自動化匹配後，各獲得 981 及 956 個特徵點，分佈如圖 4- 16。. (a) ll.tif. (b) rr.tif 圖 4- 16 地質調查近景相片小波特徵萃取成果圖. 特徵點經互相關配對後，得 79 組特徵點對，分佈如圖 4- 17。. (a) ll.tif. (b) rr.tif 圖 4- 17地質調查近景相片特徵點初步配對成果圖. 51.

(60) 經配對一致性檢查，初步剔除錯誤配對後，得 33 組特徵點配對，如圖 4- 18。. (a) ll.tif. (b) rr.tif. 圖 4- 18地質調查近景相片特徵點配對初步除錯成果圖經 TAU 測試後，進一步剔除錯誤配對，得 28 組特徵點配對，如圖 4- 19。. 圖 4- 19航照立體像對特徵點配對 TAU 測試成果圖圖 4- 19中 28 組控制點對，其殘差分佈如圖 4- 20。. 52.

(61) RMSE=0.473165 2 1.5 1 0.5 Vy. 0 -2. -1.5. -1. -0.5. 0. 0.5. 1. 1.5. 2. -0.5 -1 -1.5 -2 Vx. 圖 4- 20 地質調查近景相片經小波流程殘差分佈圖解得八參數係數如表 4- 4。表 4- 4 地質調查近景相片經小波流程求解八參數係數表係數. a1. b1. 值. 1.2754. -0.067. c1. a2. -432.8725 0.1625. b2. c2. 1.1977 -140.0435 0.0003. 經由表 4- 4所示參數，將影像重新取樣並鑲嵌後，可得圖 4- 21。. 圖 4- 21 地質調查掃瞄影像經小波流程鑲嵌合成圖. 53. a3. b3 3E-5.

(62) 4-3 近景數值影像. (a) scl.tif (1152*864). (b) scr.tif (1152*864) 圖 4- 22 近景數值影像. 第三組測試影像為利用柯達 DC210 數位式相機拍攝獲得，內容為新竹交通大學科學二館，如圖 4- 22。以洪偉嘉（1999）之流程，指定右上角點（441,173； 115,184）及左下角點（1061,593；728,608）門檻值 10，壓縮窗邊長 5，搜尋窗邊長 51，目標窗邊長 13，，經三次迭代得 19 組控制點對，如圖 4- 23。. (a) sci1.tif. (b) sci2.tif. 圖 4- 23 近景數值影像經洪偉嘉（1999）之流程獲得控制點成果圖其殘差分佈如圖 4- 24。. 54.

(63) RMSE= 0.089901 0.5 0.4 0.3 0.2 0.1 Vy. 0 -0.5. -0.4. -0.3. -0.2. -0.1. -0.1. 0. 0.1. 0.2. 0.3. 0.4. 0.5. -0.2 -0.3 -0.4 -0.5 Vx. 圖 4- 24 近景數值影像經洪偉嘉（1999）之流程殘差分佈圖所求得之八參數係數如表 4- 5。表 4- 5 近景數值影像經洪偉嘉（1999）之流程求解八參數係數表係數. a1. 值. 1.1109. b1. c1. a2. -0.0102 -370.0271 0.0487. b2. c2. a3. b3. 1.0684. -13.8092. 0.0001. -6E-6. 利用表 4- 5參數將影像重新取樣並進行鑲嵌成果如圖 4- 25。. 圖 4- 25 近景數值影像經洪偉嘉（1999）之流程鑲嵌合成圖 55.

(64) 相同像對經本研究採用之小波流程，設定梯度向量模數門檻值 25.0，互相關係數門檻值 0.9，一致性檢查門檻值 20 進行自動化匹配後，各獲得 806 及 787 個特徵點，分佈如圖 4- 26。. (a) sci1.tif. (b) sci2.tif 圖 4- 26 近景數值影像小波特徵萃取成果圖. 經互相關配對後，得 298 對特徵點對，如圖 4- 27。. (a) sci1.tif. (b) sci2.tif 圖 4- 27 近景數值影像特徵點初步配對成果圖. 經配對一致性檢查，初步剔除錯誤配對後，得 30 組特徵點配對，如圖 4- 28。. 56.

(65) (a) sci1.tif. (b) sci2.tif. 圖 4- 28 近景數值影像特徵點配對初步除錯成果圖經 TAU 測試後，進一步剔除錯誤配對，得 28 組特徵點配對，如圖 4- 29。. (a) sci1.tif. (b) sci2.tif. 圖 4- 29 近景數值影像特徵點配對 TAU 測試成果圖圖 4- 29中 28 組控制點對，其殘差分佈如圖 4- 30。. 57.

(66) RMSE=0.447161 1.5 1. 0.5 Vy. 0 -1.5. -1. -0.5. 0. 0.5. 1. 1.5. -0.5. -1 -1.5 Vx. 圖 4- 30 近景數值影像小波流程殘差分佈圖所求得之八參數係數值如表 4- 6：表 4- 6 近景數值影像經小波流程求解八參數係數表係數. a1. 值. 1.1063. b1. c1. a2. -0.0115 -367.5270 0.0459. b2. c2. a3. b3. 1.066. -11.205. 0.0001. -8E-6. 利用表 4- 6參數將影像重新取樣並進行鑲嵌後，所得成果如圖 4- 31：. 圖 4- 31 近景數值影像經小波流程鑲嵌合成圖 58.

(67) 第五章結論與建議 1. 就特徵萃取而言，小波方法對於雜訊的抵抗力優於 TDGO 運算元，對於旋轉量影響部分，或可藉由選點機制之變更獲得改善，如採取大壓縮窗依模數值選取特定數目高強度值之點位作為特徵點等。 2. 就匹配程序而言，本研究採用之小波流程利用角度累計直方圖之機制，有效解決初始旋轉量問題，互相關係數、小波轉換梯度向量模數及一致性相符點對數，三者之門檻值經適當設定後，該程序可達到自動化套合之效果。 3. 前述門檻值設定，對於套合成果之影響，以下分項討論之： (1). 互相關係數：互相關係數決定特徵點配對關係，一般而言可採用 0.75，但於 4-3 節中之近景數值影像實例，影像內涵多水平及垂直線條，重複性高，若互相關係數門檻值過低，會造成錯誤配對數目過多，導致配對除錯之難度提升，故將門檻值提升至 0.9。. (2). 小波轉換梯度向量模數：採小波方法之流程，採兩影像同時萃取特徵點後進行配對之方式，若兩幅影像中所萃取之特徵點數目過多，易造成組合爆炸之結果。互相關係數計算次數超過 10 萬次時，建議提高梯度向量模數之門檻值，以減少特徵點數目。4-1 節中航照立體像對之影像尺寸較大，以標準 25.0 之門檻值萃取，特徵點數目過高，故將門檻值提升至 120.0。. (3). 一致性檢查相符點數：採小波方法之流程進行特徵點配對初步除錯時，採用一致性檢查方式，建議之門檻值為 10 點，但於 4-3 節中，由於影像內多垂直及水平線條，故需將一致性檢查門檻提高為 20 點，以免一致性檢查結果保留過多錯誤配對。. 4. 仍就匹配程序而言，採洪偉嘉（1999）之流程相較於本研究採用之小波流程增加了：多次切割子影像控制特徵點分佈及最小二乘匹配兩機制。前者確保影像中特徵點對分佈均勻，使其鑲嵌成果精度均勻程度較小波流程為佳；後 59.

(68) 者提升右像控制點精度至次像元。 5. 本研究所使用之八參數轉換模式，僅適用於地表較為平坦地區，對於高差起伏過於劇烈之區域則並不合適，如實例一中之航空影像，建物重疊處即出現接縫偏移的現象。. 60.

(69) 參考文獻王聖鐸，1997。近景攝影測量中人工標點的自動辨識與定位，國立成功大學碩士論文：P32。李良輝，1990。計算視覺方法於 SPOT 衛星影像之幾何特性分析，國立中央大學博士論文：P29-P35。洪偉嘉，1999。自動化影像套合之研究，國立交通大學碩士論文。陳典熙，1988。PK1 連結 CCD 相機實施數值目標定位之研究，國立成功大學碩士論文。 Ackermann, F., 1984.. Digital Image Correlation： Performance and Potential. Application in Photogrammetry, Photogrammetric Record, 11(64), pp. 429439. Canny, J., 1986. A Computational Approach to Edge Detection, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-8, no.6, pp.679-698. Foerstner, W., 1984. Quality Assessment of Object Location and Point Transfer Using Digital Image Correlation, International Archives of Photogrammetry and Remote Sensing, 25(Ⅲ)ISPRS：197-219. Graps, A., 1995. An Introduction to Wavelets, IEEE Computer Science and Engineering, vol. 2, no. 2. Heipke, C., 1997. Automation of interior, relative, and absolute orientation, ISPRS Journal of Photogrammetry & Remote Sensing, 52(1997), 1-19. Hsieh, J. W., H. Y. Mark Liao, K. C. Fan, M. T. Ko and Y. P. Hung, 1997. Image Registration Using a New Edge-Based Approach, Computer Vision and Image Understanding, vol. 67, no. 2, pp. 112-130. Lin and Lee, 1995. A Simple, Fast, and Accurate Method for Image Registration, Taiwan. Lue, Y., 1988.. Research Report for National Science Council, Taiwan, R.O.C. Interest Operator and Fast Implementation, International Archives of. Photogrammetry. and. Remote. Sensing. ,. International. Archives. of. Photogrammetry and Remote Sensing, 27(Ⅱ)：491-500. Mallat, S. and S. Zhong, 1992. Characterization of Signals from Multiscale Edges, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-14, no. 7, pp. 710-732. 61.

(70) Moravec, H.P., 1977.. Towards Automatic Visual Obstacle Avoidance, International. Joint Conference of Artificial Intelligence, Vo1. 1, Tokyo, pp 584. Pope, A. J., 1975. The Statistics of Residuals and the Detection of Outliers. XVI General Assembly of the International Union of Geodesy and Geophysics, International Association of Geodesy. Rosenholm, D., 1987. Multi-Point Matching Using the Least Squares Technique for Evaluation of Three Dimensional Models, Photogrammetric Engineering and Remote Sensing, 53(6)：621-626 Schenk, T., 1997.. Towards automatic aerial triangulation, ISPRS Journal of. Photogrammetry & Remote Sensing, 52：110-121. Xu, Y. S., J. B. Weaver, D. M. Healy and J. Lu, 1994. Wavelet Transform Domain Filters：A Spatially Selective Noise Filtration Techinique, IEEE Transactions on Image Processing, vol. 3, no. 6, pp. 747-758. Zheng, Q. and R. Chellappa, 1993. A Computational Vision Approach to Image Registration, IEEE Transactions on Image Processing, vol. 2, no. 3, pp. 311326.. 62.

(71) 附錄 A：Spline 小波原理[Mallat, 1992] 本研究所使用之小波母函數為一二次雲曲線（Quadratic Spline），如附圖一(b)，細部推導如下。定義一重建小波函數 χ (x ) ，使得此函數的傅力葉轉換滿足： +∞. ∑ψˆ(2 ω )χˆ(2. ω) = 1. (1). f ∗ χ 2 j ( x). (2). j. j. −∞. 則函數 f (x ) 可經由式(2)重建：. f (x ) =. +∞. ∑W. 2j. j =−∞. 令一平滑函數 φˆ(ω ) 使得： +∞. 2 φˆ(ω ) = ∑ψˆ( 2 j ω ) χˆ( 2 j ω ). (3). j =1. 對於任何尺度大於 1 的情況，將(3)式改寫為(4)式： J. 2 2 φˆ(ω ) − φˆ(2 J ω ) = ∑ψˆ(2 j ω ) χˆ( 2 j ω ). (4). j =1. 將此一平滑函數之傅力葉轉換寫成一無限乘積型式如式(5)： +∞. φˆ(ω ) = e −iωω ∏ H ( 2 − pω ). (5). p =1. 其中 H (ω ) 為週期 2π 的可微函數使得 H (ω ) + H (ω + π ) ≤ 1 且 H (0) = 1 2. 2. (6). 上式(6)為在 L2 (R) 中使 φ ( x) 為一平滑函數之充分條件[Mallat, 1989]，式中的 ω 代表取樣位移期使 φ ( x) 以 0 為對稱，由式(1)可推廣： φˆ( 2ω ) = e − iωω H (ω )φˆ(ω ). (7). 定義一小波函數ψ (x ) ，其傅力葉轉換ψˆ( x ) 可由： ψˆ( 2ω ) = e −iωω G (ω )φˆ(ω ). (8). 其中 G (ω ) 為週期 2π 的函數使得： 2. 2. φˆ(ω ) − φˆ( 2ω ) = ψˆ( 2ω ) χˆ( 2ω ) 令 63. (9).

(72) χˆ(2ω ) = e iωω K (ω )φˆ(ω ). (10). 其中 K (ω ) 為週期 2π 的函數. 將式(7)、(8)、(10)代入(9)中可得： 2. H (ω ) + G (ω ) K (ω ) = 1. (11). 則. K (ω ) =. 1 − H (ω ). 2. (12). G(ω ). 若得知 H (ω ) 及 G (ω ) ，便可依式(12)、(10)推求出 K (ω ) 乃至 χˆ(ω ) ，以便進行函數之重建。函數的選取，欲使小波函數ψ (x ) 為平滑函數 θ ( x) 的一階導數，且盡可能的具備規律性（regular）及反對稱性（antisymmetrical），可以一週期為 2π 之函數家族滿足上述限制： H (ω ) = e iω / 2 (cos(ω 2)) 2n+1. (13). G(ω ) = 4ie iω / 2 sin( ω / 2). (14). K (ω ) =. 1 − H (ω ). 2. (15). G(ω ). 由式(5)、(8)可推得  sin( ω / 2)  φˆ(ω ) =    ω /2 . 2 n +1. (16).  sin( ω / 4)  ψˆ(ω ) = i ω    ω/4 . 2 n+ 2. (17). 由上述的推導可得原函數 θˆ(ω )  sin( ω / 4)  θˆ(ω ) =    ω /4 . 2n+ 2. (18). 為使：1.小波函數不以 0 為中心成對稱，2.尺度函數為以 0 為中心成對稱， 1 於本例中將平移常數設定為 ω = ，由此可以證明出小波函數ψ (x ) 為一二階雲曲 2. 線，而函數 θ ( x) 為一三階雲曲線，且其積分值為 1。同時 H (ω ) 、 G (ω ) 、 K (ω ) 可 64.

(73) 視為三轉換函數，代表具備有限脈波響應（finite impulse reponse）特性的離散式濾波矩陣，相對應的脈波響應量如附表一所示。附圖一. (a) 三次雲曲線. (b) 二次雲曲線. 附表一 n. H. G. K. -3. 0.0078125. -2. 0.054685. -1. 0.125. 0.171875. 0. 0.375. -2.0. -0.171875. 1. 0.375. 2.0. -0.054685. 2. 0.125. -0.0078125. 3. 65.

(74) 作者簡歷姓名：陳維勇籍貫：江蘇省海門縣出生日期：民國 65 年 4 月 17 日學歷：台灣省立花蓮高級中學畢業國立交通大學土木工程學系畢業國立交通大學土木工程學系碩士班肄業. 66.

(75)