Volume 12, No.4, December 2007, pp. 381-401
利用多尺度匹配自動套疊衛星雷達影像
洪浩倫
1王志添
2陳錕山
3摘要
合成孔徑雷達具有全天候特性，可及時監測提供災害資訊，已經成為環境監測的重要工具。本文使 用不同時期、有相同地區重疊的重複軌道拍攝之雷達影像，利用重複軌道拍攝幾何關係相似特性，像元 攝像入射角度接近，在地表起伏地區，雷達影像幾何變形特徵一致，經過以多尺度匹配法為主的處理過 程，以及共軛點篩選，最後進行濾除斑駁雜訊的步驟，有效地在雷達影像上找到均勻分布且精確可靠的 共軛點，進行影像自動套疊，建置衛星雷達影像自動化套疊與變遷分析，以減少人為誤差及提昇效率。
本研究成果將可提升國土監測敏感度與滿足全天候山區崩塌監測實務需求。
關鍵詞: 合成孔徑雷達、自動套疊、多尺度匹配法
1. 簡介
台灣地處環太平洋地震帶，是典型的板塊碰撞 下產生之大陸邊緣島嶼，造山運動迄今仍在持續的 進行，經常引發台灣密集的地震活動；此外，台灣 西隔台灣海峽連接歐亞板塊，東臨太平洋，位處於 世界季風最發達區域內，季風常隨季節轉變而有顯 著改變，梅雨季與颱風帶來旺盛雨量，極易造成山 區土石流等地質災害。而當災害產生時，必須即時 掌握災害在何處發生以及災害面積，進行災害評 估，尤其在偏遠山區，可能因為環境地形因素，無 法深入山區，而人力資源有限，也無法短時間內探 察災情。而近年來，隨著不同功能與週期衛星發射 升空，而衛星上所載的感測器包含不同光譜解析 力、空間解析力、時間解析力以及輻射特性等，能 快速獲取大範圍資料，可對全球進行環境監測 (Chalermwat, 1999)。
在 遙 測 技 術 上 ， 合 成 孔 徑 雷 達 (Synthetic Aperture RADAR, SAR)是一種主動式的遙測系統，
主要用於測定大範圍地區的地表資料。使用電磁波
譜裡的微波波段，具有良好的大氣穿透力，適用於 日夜的偵測、以及任何天氣條件下獲取資料的能 力，因此，在雷達影像的實務應用上，扮演重要角 色，如變遷偵測、合成孔徑雷達干涉法求定數值高 程模型(Digital Elevation Model, DEM)…等，都需要 對影像進行控制點的量測，即共軛影像上的匹配點 選取(陳良健，1998; 黃金城，2002)。但由於雷達 系統成像性質的影響，導致影像解析力不如航測影 像，且斑駁雜訊干擾嚴重，造成選取共軛點位的困 難(Wang and Chen, 2003)；因此本文主要目的乃發展 自動化雷達影像套疊技術以降低人為誤差並提高 效率。
2. 理論基礎
2.1 影像套疊基本原理
隨著不同功能與週期衛星發射升空，而衛星上 所載的感測器包含不同光譜解析力、空間解析力、
時間解析力以及輻射特性等，對全球進行環境監
收到日期:民國 96 年 09 月 14 日 修改日期:民國 96 年 12 月 20 日 接受日期:民國 96 年 12 月 21 日
1國立中央大學太空科學研究所碩士生
2國立中央大學太空及遙測研究中心副研究工程師
3國立中央大學太空及遙測研究中心及太空科學研究所教授
測。而多重遙測影像之資訊必須藉由數位影像間幾 何位置套疊來建立，即影像套疊，其利用座標轉換 來建立兩張或多張影像間的幾何套疊關係，以定義 參考影像(reference image)或主像(master image)的座 標 系 為 參 考 座 標 系 (Schenk,1999; Hong and Schowengerdt,2005; Barbara and Jan, 2003)。而影像套 疊的方法是根據處理過程中匹配方法不同而分成 以區域為基礎的套疊(area-based registration)，以及
以特徵為基礎的套疊(feature-based registration)( Dare and Dowman, 2001; Liu, B W., 2003)。
影像匹配的目的在於搜尋出影像間相對應特 徵向量的位置，亦即共軛點，共軛點是用來連接或 套疊相鄰的影像，所選取的共軛點必須是在影像重 疊區域中明顯易分辨的點位，如圖 1 (紅色十字位 置)，且應盡可能地均勻分佈於影像周圍。以下將 介 紹 影 像 匹 配 技 術 之 區 域 匹 配 (area-based matching)。
2.2 區域匹配
使用區域匹配於影像中進行匹配尋找共軛點 是航測中最常見的方法之一，其概念是利用地表物 的回波能量值，以成像於影像上的輻射資訊及灰度 值為匹配實體，在兩影像中各取一子視窗，直接對 視 窗 內 灰 度 資 訊 進 行 計 算 作 為 相 似 性 評 估
（Similarity Assessment）的指標，其優點是計算簡 單，可適用於自動化，且配合良好的匹配策略可得 到精確與可靠的成果，又依評估方法不同，分為標 準化互相關法(Normalized Cross Correlation, NCC)與 最小二乘匹配法(Least Squares Matching, LSM)。在 輻射及幾何差異不大的影像上，區域匹配法之精度 較高，但必須事先給定較精確的初始位置(Lin, J. B., 1996; Wang, 2004)。
在區域匹配法中是以影像的灰度值為評估影 像間相似性的計算單元，在本研究中，實驗影像為 重複軌道衛星雷達影像，由於蒐集來自地物相近的 回波資訊，使得其成像的灰度值差異接近，唯影像 斑駁雜訊干擾為隨機，仍需剔除，而幾何變形特徵 (如前坡縮短、疊置)也相似，針對此類型資料，可 適用影像的灰度值進行影像間相似性評估，來找到 影像中的共軛點，作為套合時的初步共軛點對。
3. 研究方法與流程
本研究主要是針對重複觀測軌道的前後期衛 星雷達影像，進行自動化套疊程序，流程如圖 2 所 示。整個作業流程建立在多重尺度匹配處理上，主 要可分為初始套疊以及精細套疊兩部分，而最終目 圖 1. 共軛點連結影像資料(紅色十字位置)
的是希望自動化於全張衛星雷達影像上找到儘可 能均勻分佈、且可靠度高的大量影像共軛點，以供 後續影像套疊或環境變遷偵測之用。
第一部分，初始套疊。在此利用全自動方法。
首先，分別對前後期影像作小波轉換，使用的小波 基底函數為 HAAR 小波(因運算效率高，符合自動 化，故選用 HAAR 小波)，來保留低頻影像，可達 到降取樣，節省運算時間，以及濾除斑駁雜訊效果 (顏世坤，1994)。再分別從前後期低頻影像選擇一 地塊影像，進行標準化互相關分析，可求得低頻影 像間的初始位移量。利用此位移量，在前後期低頻 影像間，進行區域匹配，求得大量的共軛點對，即 共軛點的影像位置表。
進行精細套疊步驟之前，必須先將初始套疊得 到的共軛點進行篩選，保留可靠度高的共軛點，並 逐層逐層地代回原始影像層。共軛點篩選可分為逆 方向區域匹配之檢驗、空間特徵分析、影像金字塔 匹配以及高相關性以及均勻分布，最後可剔除匹配
失敗的共軛點，將其餘共軛點均勻分佈在原始影像 層。
第二部份，精細套疊。透過初始套疊處理之後 已建立起前後影像間的初始對應關係，此時再利用 共軛點所處的地塊影像，作補零處理，找到更精確 的共軛點對之位置，達到次像元精度。將上步驟得 到的共軛點位對，使用二次多項式，以最小二乘法 來求得轉換矩陣，為最佳解。利用轉換矩陣，得到 原影像位置與轉換後的新影像位置，進行粗差偵 測，利用經驗值，大於 3 個像元，迭代剔除。將迭 代剔除後的共軛點位置進行重新取樣，將後期影像 精細套疊在前期影像上，得到精細套疊影像。
後續處理，本研究套疊成果經斑駁雜訊濾 除處理後的套疊影像，進行變遷偵測，紀錄變遷處 的影像位置，並導入數值高程模型(DEM)，可得知 發生變異區的大地二度 TM 座標，即時提供資料供 地質災害、環境變遷時的參考。
圖 2. 本研究之自動化套疊衛星雷達影像作業流程
3.1 影像的初步套疊
初步套疊處理，針對雷達影像的共軛點人工點 選困難問題，利用了自動化套疊方法可降低人為誤 差與增加效率，初步套疊流程如圖 3 所示。
3.1.1 影像金字塔
因本研究之實驗影像為衛星雷達影像，影像範 圍廣大，當空間解析度約 20 公尺時，影像尺寸可 能達一萬像元單位級，且影像中存在斑駁雜訊，若
直接將前後期影像進行匹配，其搜尋範圍大導致計 算量變大，且錯誤匹配的機率也會提高。在此便引 入上述之階層式的影像金字塔策略，利用小波多重 解析力，影像在經過金字塔結構的約化處理後，匹 配搜尋範圍縮小，斑駁雜訊對於影像特徵影響也降 低了，如圖 4 所示，因此將影像約化至最頂層後(粗 解析力)，先以此層影像進行影像平移估算，再以 估算量為輔助，進行區域匹配處理能得到較可靠的 結果，利用這組共軛點作為下一層匹配起始對應位 置，由上到下直到原始影像層為止。
圖 3. 影像初步套疊流程
圖 4. 實驗影像經過小波轉換測試之結果 (a)原始影像 (b)小波轉換後之結果 (左上為低頻部份，其餘為高頻部份)
3.1.2 估算影像間平移量
利用前步驟影像約化後的低頻影像資訊，於前 期影像選擇搜尋視窗，後期影像取出一範圍較小的 目標視窗，進行標準化互相關分析，數學模式如式 (1)，求得標準化互相關性分佈圖，找出峰值的位 置，以估算兩張影像間的平移量，作為下一步區域 匹配的依據。為了測試標準化互相關分析能否正確 估算出影像間的平移量，故選用一張影像當測試影 像，再將此張影像邊緣區域裁切至相反邊緣，也就 是左平移至右，下平移至上的方法，得到平移後之 影像，如圖 5 所示意，
將測試影像與平移後之影像帶入標準化互相 關數學模式，如式(1)，得到標準化互相關分佈圖，
再計算峰值位置與分佈圖中心的位置之差值，即為
兩張影像間的平移量，最後再與原假設平移量作比 較，檢驗平移量是否正確。
測試影像尺寸為 X、Y 方向分別為 492 像元、
502 像元，而 X、Y 方向的平移量皆為 10 個像元寬，
圖 6 分別為測試影像與平移後之影像。
再以測試影像中心為視窗原點，取出搜尋視 窗，同樣的，以平移後的影像中心為視窗原點，取 出一範圍較小的目標視窗，進行標準化互相關分 析，求得標準化互相關性分佈圖，如圖 7 所示，得 到標準化互相關性分佈圖後，找出分佈圖的峰值位 置 (201,201) ， 峰 值 為 1.0 ， 分 佈 圖 中 心 位 置 (191,191)，峰值與中心點位置差值為(10,10)，代表 兩張影像間的平移量為(+10,+10)，與原假設平移量 相同。
圖 5. 影像平移測試流程之示意圖
圖 6. 測試影像(左圖)與平移後之影像(右圖)
圖 7. (a)搜尋視窗 (b)目標視窗 (c)標準化互相關性分佈圖
圖 8. 區域匹配之目標視窗、搜尋視窗與平移量估算值相關示意圖
圖 9. 逆方向區域匹配之目標視窗、搜尋視窗相關示意圖
3.1.3 區域匹配
在本研究中，實驗影像為重複軌道衛星雷達影 像，由於蒐集來自地物相近的回波資訊，使得其成 像的灰度值差異不大，唯影像斑駁雜訊干擾為隨 機，仍需剔除，而幾何變形特徵(如前坡縮短、疊 置)也相似。
故針對此類型資料，可先利用前步驟所得之影 像平移估算值來作為搜尋視窗位置之依據，如示意 圖 8，步驟說明:
1. 根據視窗大小(本研究訂為 16 x 16 像元)，於後 期影像自訂起始點，得到目標視窗中心之影像 位置(x1,y1)，為 tp1。
2. 於前期影像中，由影像位置(x1,y1)，再加上影像 平移量估算值(△x, △y)，得到(x1+△x, y1+△y) 影像位置，以此為視窗中心，取出搜尋視窗(本 研究訂為 32 x 32 像元)。
3. 利用搜尋、目標視窗進行影像相似程度評估，
利用標準化互相關分析法找出互相關性分佈 圖峰值的位置，以求得最佳的匹配位置，即共 軛點，TP1。
4. 將目標視窗由左下至右上依序移動，視窗移動 間格為等距(本研究訂為 16 像元)，重複 1~3 步 驟，依次在前期影像中找到一個共軛點，TP2、 TP3…，最後獲取一組初步共軛點對{TPn , tpn} 並紀錄之，n 為共軛點個數。
5. 獲取{TPn , tpn}新共軛點對。
使用平移估算量為區域匹配之輔助，優點是減 少蒐尋共軛點的範圍，不需要整張影像搜尋，只需 在某一區域搜尋，可大量降低電腦運算時間。
3.1.4 共軛點篩選
將區域匹配得到的初始共軛點進行離群値剔 除，也就是將不可靠的共軛點剔除，以下三種簡單 規則來得到可靠度高且分佈均勻的共軛點對:逆向 匹配、空間特徵分析以及高相關性篩選。
另外使用金字塔匹配法，將較高階(粗解析度) 影像層之共軛點對，代回上一層低階影像層(細解
析度)，重複此方法，一層一層代回原始影像層，
得到原始解析力影像之共軛點對。
3.1.4.1 逆方向區域匹配之檢驗
步驟目的為檢驗共軛點對之區域唯一性，也就 是局部相似程度。示意圖如 9。
步驟說明:
1. 由上步驟所得共軛點對{TPn , tpn}，於前期影像 以第一個共軛點(TP1)位置為中心，取出目標視 窗(本研究訂為 16 x 16 像元)。
2. 再於後期影像中，由第一個共軛點(tp1)位置，
以此為視窗中心，取出搜尋視窗(本研究訂為 32 x 32 像元)。
3. 使用相同方法，利用搜尋、目標視窗進行影像 相似程度評估，利用標準化互相關分析法，如 式(1)，找出互相關性分佈圖峰值的位置，以求 得最佳的匹配位置，共軛點，tp1′。
4. 若 tp1′與 tp1為後期影像中相同位置時，代表 局 部 相 似 程 度 最 高 ， 保 留 第 一 個 共 軛 點 對 (TP1,tp1)，反之則剔除第一個共軛點對。
5. 依順序將每一個共軛點對，重複 1~4 步驟，檢 驗每一個共軛點，最後獲取一組初步共軛點對 並紀錄之，為逆向匹配最後結果。
6. 獲取{TPn , tpn}新共軛點對。
3.1.4.2 空間特徵分析
通過逆匹配檢驗之共軛點雖然皆具有局部最 相似性，但因為找到的共軛點位置可能存在於均質 區域，如海洋、湖泊、平原地區，屬於弱空間特徵 區域，特徵不明顯，相似性局部最高，但可信度低，
所以必須剔除此類共軛點對，保留強空間特徵區 域，如山脊線夾角處。步驟說明如下:
1. 由前節所得共軛點對{TPn , tpn}，針對前期影像 中第一個共軛點 TP1，以此設定一目標視窗(本 研究訂為 16 x 16)，將視窗內灰值代入自相關 數學模式，如圖 11，輸出結果以等值線模式顯 示。
2. 根據等高線模式來分析結果，紀錄數值大於 0.9(因自相關分析，峰值為 1.0)的像元個數佔
自相關分析結果之總像元的百分比，如圖 10 所示之百分比。
3. 依順序將前期影像中共軛點，重複 1~2 步驟，
將所有共軛點根據等高線模式來分析結果，求 得百分比值並紀錄之。
4. 由步驟 3，將所有的百分比值以直方圖分佈來 表示，保留百分比值小於門檻值的共軛點，也 就是保留強空間特徵點的資訊，為最後自相關 分析結果。
5. 獲取{TPn , tpn}新共軛點對。
最後歸納，從等高線模式中發現，當共軛點區 域為強空間特徵時，等高線近峰值區域會趨向非對 稱(asymmetric)，而弱空間特徵時，等高線近峰值區 域則會近乎對稱(symmetric)，故利用此空間特性來 分析，此為定性分析，如圖 10。
而定量分析，可由百分比值可得知當共軛點區 域為強空間特徵時，百分比會偏低，而弱空間特徵 時，百分比則會相對偏高，如圖 11 之百分比值。
將所有的百分比值以直方圖分佈來表示，保留百分 比值小於門檻值(本研究訂為總個數之 95%)的共軛 點，也就是保留明顯的空間特徵點資訊，如圖 12。
圖 10. 共軛點區塊影像與自相關分析結果
圖 11. 明顯空間特徵(左圖)、弱空間特徵(右圖)與其自相關分析結果
3.1.4.3 高相關性以及均勻分布
在區域匹配步驟所匹配出的共軛點，雖然經過 前述共軛點篩選(逆方向區域匹配、空間特徵分 析)，但因為保留之共軛點為局部最相似，為相對 最大值，但由整體共軛點看來，其中可能涵蓋了低 相關性，如低於 0.6 之共軛點對，而又顧及共軛點 需於影像中均勻分布(Hiep et.al., 2004)，本研究以下 述規則來篩選。
本步驟分為兩部分，第一部份:全區域共軛點 的高相關性篩選，篩選門檻由共軛點對，取其最高 相關程度之 90%為篩選標準，保留門檻值以上之共 軛點對；第二部份:由於第一部份篩選後之共軛點
於影像中可能有某些地區是屬於絕對低相關性，也 就是無共軛點存在，此結果會導致套疊影像會發生 扭曲現象，影響最後的分析可信度，所以運用一移 動視窗，從影像中屬於絕對低相關性的區域，保留 一個局部最大值，最後保持共軛點對在整張影像中 均勻分布。流程圖如圖 13 所示，其中綠色星號為 相似性程度大於門檻值，紅色三角形為小於門檻 值。
步驟說明:
1. 圖 13(a)為分布於原始影像之共軛點示意圖。
2. 圖 13(a)為保留大於門檻值之共軛點分布，門檻 值由共軛點對，取其最高相關程度之 90%為篩 選標準。
3. 圖 13(c)、(d)為移動視窗過程，由影像左下至 右上，移動間隔為單一像元。
4. 圖 13(e)，當移動視窗搜尋至某區域，發現其區 域內皆為相關性低於門檻值之共軛點，在步驟 2 已刪除，為了保持共軛點均勻分布，故在此 區域選擇依局部最大相似之共軛點對並保留。
5. 獲取{TPn , tpn}新共軛點對。
圖 12.明顯空間特徵的保留
95%
(a) (b) (c)
(d) (e)
圖 13. 高相關性及分布均勻之篩選流程示意圖
3.2 影像的精細套疊
經過初步套疊處理後，獲得影像初步套疊之共 軛點對後進行精細套疊處理如下:
3.2.1 次像元精度
本研究針對共軛點之區域影像進行補零處理 (zero padding)，也就是將影像由空間域利用傅立葉 轉換轉至頻率域，針對訊號尾端補零，再逆傅立葉 轉換，進而得到網格加密影像，如圖 14 所示。當 然補零處理亦可以小波轉換進行，其間優劣須進一 步檢驗。
在遙測影像的應用上，如變遷偵測、影像上特 徵位置量測等，皆要求結果誤差須小於一個像元，
即次像元精度，將前期、後期原始影像上之共軛點 對，經補零處理後之加密影像，以此進行匹配動 作，找出更精確之共軛點對，此步驟可以改善套疊 的精準度，使得變遷偵測之位置能更準確。
3.2.2 多項式擬合
影像匹配得到共軛點對之後，評估兩影像間變 形的情況後，選擇適當的轉換模式進行套疊參數的
解算。一般常見的轉換模式包括：線性相似轉換(四 參數)、仿射轉換(六參數)、投影轉換(八參數)、或 使用更高階的多項式轉換 [Brown，1992]。而不同 的轉換模式對於影像間的變形具有吸收或凸顯的 作用，因此對於轉換模式的選擇必須考量到影像間 成像幾何的因素。
3.2.2.1 最小二乘法
本研究的實驗影像為衛星雷達影像，因為涵蓋 範圍較大，受地球曲率影響，且影像涵蓋山區，所 以在共軛點位關係的建立，使用多項式來計算其轉 換參數，因此使用二次式，再進行下步驟粗差偵 測。如下式(1):
⎩⎨
⎧
+ + +
+ +
=
+ + +
+ +
=
2 5 2 4 3
2 1 0
2 5 2 4 3
2 1 0
i i i i i i i
i i i i i i i
y b x b y x b y b x b b Y
y a x a y x a y a x a a
X (1)
式中，a0~a5、b0~b5 為幾何轉換參數，Xi、Yi 是後期影像之共軛點座標，xi、yi 是前期影像之共 軛點座標。
3.3 粗差偵測
將這些共軛點進行轉換參數的求解計算後，所 有共軛點會有對應的殘差，殘差是某個方向上原輸 入坐標和轉換後坐標之間的距離
2 2
yi xi
i
R R
R
= + (2) 其中，R
i:共軛點的差合,
R
xi:共軛點的 x 方向殘差,Rxi =Xi −Xˆi,R
yi:共軛點的 y 方向殘差,R
yi =Y
i −Y
ˆi. 每個共軛點有其對應之殘差，當殘差大於某値 時，表示為匹配錯誤的共軛點，改稱為粗差，也就 是說，匹配錯誤的共軛點所對應的粗差值通常會較 大，因此所有共軛點的殘差可透過統計的方法進行 粗差偵測，進行整體一致性的檢查。在測量上常將 粗差值大於三倍的套疊精度標準視為匹配錯誤予 以除，而套疊精度標準為一個像元，合乎應用要 求，故在本研究中以粗差值大於 3 個像元視為錯誤 匹配點位，並予以剔除。最後前後期影像整體匹配 圖 14. 補零處理流程 (a)原始影像 (b)傅立葉轉換後之結果 (c)於訊號尾端補零値 (d) 網格加密影像
(a) (b)
(c) (d)
品質， 是以均 方根誤差(Root Mean-Square Error, RMSE)。前後期影像匹配的均方根誤差可計算如 下：
∑
= 1 2
Ri
RMSE n (3) 每個共軛點有其對應之殘差，當殘差大於某定 値時，代表是匹配錯誤的共軛點，改稱為粗差，也 就是說，匹配錯誤的共軛點所對應的粗差值通常會 較大，因此所有共軛點殘差可透過統計的方法進行 粗差偵測，進行整體一致性的檢查。一般在測量上 常將粗差值大於三倍的套疊精度標準視為匹配錯 誤予以除，而套疊精度標準為一個像元，合乎應用 要求，故在本研究中以粗差值大於 3 個像元視為錯 誤匹配點位，並予以剔除。
3.4 影像重新取樣
在經過幾何轉換參數計算之後，影像座標通常 不會恰好落於原先所設定的整數影像格點上，因此 經過轉換後的影像，在處理上會以一內插模式來進 行影像的再取樣，實務應用上最常見的方法為最近 相鄰法、雙線性內插法和立方迴旋法。最近相鄰法 是取最近點的灰度值作為網格上的灰度值，不會影 響原來的對比，且內插後所產生的新灰度值及影像
內的一些細微變化不致於損失過大，不失真，因此 對分類效果與變遷分析較佳，但是原來線性構造如 道路，可能會產生鋸齒狀的構造，其優點是計算時 間最為經濟。雙線性內插法是取最近的四個點的權 值平均，較最近相鄰法不易造成像元間的損失，整 題來說較平滑，處理時間為最近相鄰法的三倍。而 立方迴旋法是取最靠近的 16 個灰度值進行三次方 多項式最適化以定出網格上的灰度值，雖然可以產 生較為圓滑的結果，其視覺效果最佳，但其計算時 間為最近相鄰法的十倍。因本研究為自動化流程，
且後續為變遷偵測分析，故選擇最近相鄰法。
4. 研究試區
本研究所提出的自動套疊流程，乃是針對重複 軌道(repeat-pass)之衛星雷達影像進行處理，以兩組 不同衛星感測器所拍攝，成像後經輻射校正之多觀 點振幅影像為實驗對象，分別為 Radarsat-1 以及 ALOS 衛星雷達影像，進行分析探討。
實驗區影像分別為 Radarsat-1 衛星以及 ALOS 衛星所拍攝之台灣中部地區前後期雷達影像，如圖 15 (Radarsat-1，前、後期影像)與圖 16(ALOS，前、
後期影像)所示。
圖 15. 台灣中部地區 Radarsat-1 衛星雷達影像: 2006/04/28(左)， 2006/06/15 (右)
圖 16. 台灣中部地區 ALOS 衛星雷達影像:2006/12/31(左), 2007/02/15(右)
圖 17. RADARSAT-1 雷達影像之套疊共軛點及檢核點分布圖(a)前期影像 (b)後期影像
(a) (b)
圖 18. RADARSAT-1 雷達影像之套疊結果(R: 2006/04/28 ; G/B:2006/06/15)
5. 研究結果
5.1 RADARSAT‐1 雷達影像自 動化套疊成果
圖 17 為本研究自動化套疊成果之匹配點分布 情形，紅色三角形代表套疊共軛點位置，黃色三角 形代表檢核點位置。
經過李氏濾波器(濾波視窗為 5x5)濾除斑駁雜 訊，圖 18 為自動化套疊成果之套疊影像。依日期 分別套色，前期影像以紅色顯示，後期影像以綠色 /藍色顯示。
圖 19 為由 266 個檢核點之套疊誤差向量圖，
其向量長度經過 400 倍的放大，觀察整體分佈並無 發現顯著的系統性偏移。圖 20 為取樣(samples)、線 列(lines)方向之誤差直方圖，可得知共軛點於取 樣、線列方向分布均勻。表 1 則是像座標檢核精 度。取樣、線列方向均方根誤差皆小於ㄧ個像元, 達到解析力下的合理要求。
sample line 平均誤差(像元) -0.06 -0.11
均方根誤差 3.81 公尺 4.75 公尺
取樣數 266
圖 19. RADARSAT-1 套疊誤差向量圖
表 1. RADARSAT-1 套疊結果之像座標檢核精度
圖 20. RADARSAT-1 套疊 X 方向(sample)、Y 方向(line)誤差直方圖
圖 21 幾何校正後之變遷偵測(R:大於 6dB；G:小於-6dB)
圖 22. 區域 A(台中縣和平鄉大甲溪流域)之(a)前期影像 (b)後期影像 (c)變遷情形
(a) (b)
(c)
馬鞍水壩
圖 23. 區域 A(台中縣和平鄉大甲溪流域)之 SPOT-5 衛星影像(R: NIR, G: Red, B: Green)
圖 24. 區域 B(台中縣和平鄉知馬漢山)之(a)前期影像 (b)後期影像 (c)變遷情形
圖 25. 區域 B(台中縣和平鄉知馬漢山)之 SPOT-5 衛星影像(R: NIR, G: Red, B:Green)
(a) (b)
(c)
將圖 18 濾除雜訊後的套疊影像，轉變成能量 衰減圖，依據變遷之經驗值[ ]，以-6dB 與+6dB 為 變遷門檻 T，並分別套色為 R:大於 6dB；G:小於 -6dB。最後利用農航所所提供之全台 DEM(解析力 40 公尺，內插為 25 公尺)，以及地面控制點，將變 遷位置經過幾何校正後，轉換至二度 TM 座標 (TWD67)，提供變遷區域之地理座標如圖 21 所示。
針對圖 21 中區域 A、B，分別探討。區域 A 之地理位置為為台中縣和平鄉大甲溪流域中游，區 域 A 中心座標為(240836.3mE, 2675625.7mN)，如圖 22(a), (b)所示前後期影像，圖中紅色箭頭位置為和 平鄉之馬鞍水壩；圖 22(c)為前後期影像之變遷情 形，以能量衰減(dB 值)表示。利用 SPOT-5 衛星影 像作為佐證，SPOT-5 之綠光波段可反應出水的存 在，在流域轉折處可明顯看出其變異，如圖 23。
區域 B 之地理位置為知馬漢山，位於台中縣和平 鄉大甲溪之匹亞桑溪流域，B 區域中心座標為 (257472.7mE, 2687974.6mN)，如圖 24(a)(b)，圖 24(c) 為前後期影像之變遷情形。再利用 SPOT-5 衛星影 像作為驗證， SPOT-5 之近紅外光波段可反應出植
被的存在，前期影像中，近紅外光反應較強烈(黃 色圈選區域)，明顯可見其差異，如圖 25 所示。
5.2 ALOS 雷達影像自動化 套疊成果
圖 26 為 ALOS 影像經自動化套疊之成果(已濾 除斑駁雜訊，濾波視窗為 5x5)。依日期分別套色，
前期影像以紅色顯示，後期影像以綠色/藍色顯示。
圖 27 為匹配成果 245 個檢核點之套疊誤差向 量圖，其向量長度經過 500 倍的放大，觀察整體分 佈發現系統性偏移，且集中在影像中間區域，產生
「交錯」現象，與先前成果大不同，發生原因猜測 可能為影像幾何成像時的錯誤，也可從表 2 之檢核 精度中來思考，取樣方向均方根誤差大於ㄧ個像 元，無法吸收掉此方向誤差，才造成誤差向量圖之
「交錯」現象。圖 29 為取樣(samples)、線列(lines) 方向之誤差直方圖。
針對圖 26 之區域 C、D，分別作前後期影像之 比較，區域 C 為彰化縣王功漁港，如圖 29,其中黃 色圈選 1、2 區域皆有明顯變異處。變遷原因可能 從 2006 年底至今年 2007 年 2 月中，觀光旅遊局正 圖 26. ALOS 雷達影像之套疊結果 (R: 2006/12/31 ; G/B:2007/02/15)
於王功地區進行王功漁港景觀相關改善工程。區域 D 為龍巖農場與馬光農場(雲林縣褒忠鄉)，圖 30 所示。此區域為農場，種植作物，常會隨著季節不 同而有所變化，可由相同極化方向之套疊影像，比 較出不同時期之地物變化，可供紀錄作物輪耕週期 之參考。
表 2. ALOS 雷達影像套疊結果之像座標檢核精度 sample line 平均誤差(像元) 0.10 -0.03
均方根誤差 10.68 公尺 2.25 公尺
取樣數 245
圖 27. ALOS 雷達影像套疊結果之誤差向量圖
圖 28 ALOS 套疊結果之 X 方向(sample)、Y 方向(line)誤差直方圖
1
2
圖 29. 區域 C (彰化縣王功漁港)之前期影像(a)，後期影像(b)，套疊影像(c)，R: 2006/12/31 ; G/B:2007/02/15
圖 30 區域 D (龍巖農場與馬光農場)之前期影像(a)，後期影像(b)，套疊影像(c)，R:2006/12/31; G/B:2007/02/15 (a)
(b)
(c)
(a)
(b)
(c)
6. 結論與建議
6.1 結論
根據實驗過程、成果與分析，提出以下結論與 建議:
(1)建置衛星雷達影像多尺度自動套疊處理機 制，目前用於山區崩塌監測，已有不錯的效果。
(2)本研究主要貢獻在於針對重複軌道之衛星 雷達振幅影像，透過全自動化的套疊流程，於 RADARSAT-1 以及 ALOS 的雷達影像上進行全自 動化影像共軛點的篩選，最後利用此共軛點位進行 影像套疊，以及變遷偵測。希望利用此自動化套疊 流程，能有效取代以往人工進行大量共軛點位的選 取，減少人為誤差及提昇整體效率。
(3)在影像初始位移估算中，利用影像金字塔最 高層影像(粗解析力)之低頻資訊以及標準化互相 關分析，可找出前後期影像間的平移量，並以此作 為區域匹配時的搜尋範圍之依據，不必搜尋整張影 像，省下大量時間。再利用影像金字塔由粗而細的 匹配策略下，由金字塔上一層的匹配成果提供下一 層較佳的起始共軛點位位置，直到原始影像，得到 最終的共軛點對，供套疊影像使用。
(4)在影像金字塔匹配過程，經過逆向匹配、空 間特徵分析、高相關性且均勻分布之共軛點篩選規 則，能有效將共軛點所處位置，具有局部特徵唯一 以及非均値區域(如海面、湖泊)保留下來，以及使 共軛點盡可能均勻分布於整張影像，進而提高共軛 點位的相關性，以及降低影像套疊時的誤差，此為 本研究另一貢獻。
(5)自動化套疊成果之精度評估，在誤差向量 圖，除了 ALOS 之成果有「交錯」現象之外，
RADARSAT-1 皆無系統性偏移，以及均方根誤差在 取樣、線列方向皆小於一個像元，故套疊成果之後 續分析(如變遷偵測)，符合一般應用標準。
(6)自動化效率而言，以兩張前後期影像，影像 尺寸皆約為 3000 x 3000 像元大小，整個自動化影 像套疊流程所花費的電腦運算時間約為 34 分鐘，
而當影像尺寸加大至 10000 x 10000 像元時，所花 費的電腦運算時間約為 2 小時，電腦配備為 CPU：
AMD Athlon 64 Processor 1.81GHz、RAM：1GB，使 用軟體則為 IDL 6.2 與 PCI V9.1。
6.2 建議
(1)本研究所使用影像必須經過輻射校正，取得 影像振幅值的整體一致性，其原因有兩點，第一 點，因本研究匹配方法為區域匹配，利用影像振幅 値為運算單元，來尋找最相似之點位，若無輻射校 正會導致匹配失敗。第二點，由於套疊成果後續應 用，如變遷偵測，本研究利用能量衰減情形(dB 値) 來判斷變遷區域，若無輻射校正，會使得變遷偵測 結果產生錯誤。
(2)本研究方法所適用資料為重複軌道之衛星 雷達振幅影像，故重複軌道拍攝幾何關係相似，像 元攝像入射角度接近，在地表起伏地區，雷達影像 幾何變形特徵一致，存在重複軌道時可能產生的平 移量，與些微尺度上的變化，可利用本研究自動匹 配方法找出共軛點對，再利用二次多項式來擬合，
求得轉換後之影像位置。但是當影像間尺度上變化 過大，或者影像分別為昇模式(Ascending mode)、下 模式(Descending mode)時，因不同視角(look angle)，
自動化匹配方法需要再修正。
近年衛星遙測發展形成環境監測工具，其中以 衛星雷達之全天候特性，且具有重複軌道觀測功 能，可及時監測提供災害或變遷位置。本研究之自 動化影像套疊流程，即可利用公尺級超高解析(約 1~3 公尺)衛星雷達 Radarsat-2 與 TerraSAR-X 即將升 空運作之優勢，進行變遷偵測，來提升國土監測敏 感度，達到小區域崩塌及時監測，或水情監測(水 庫集水區、流域等)。除了可利用超高解析之衛星 雷達影像來進行自動化套疊，偵測出前後期影像變 遷位置之外，尚可利用多偏極雷達影像資訊，來獲 取山坡表層滑動之變遷。
參考文獻
陳良健，1998。 衛星影像間全自動化套合之研 究 ，行政院國家科學委員會專題研究計畫成 果報告。
黃金城，2002。衛載 SAR 地塊影像匹配之參數最 佳化，碩士論文，土木工程研究所，國立中央 大學。
顏世坤，1994。Wavelet 轉換於遙測影像之分析及 其應用，碩士論文，太空科學研究所，國立中 央大學。
Barbara, Z. and Jan, F., 2003. Image registration methods: a survey, Department of Image Processing, Institute of Information Theory and Automation, Academy of Sciences of the Czech Republic.
Chalermwat, P., 1999. High Performance Automatic Image Registration for Remote Sensing, Phd thesis, George Mason University, Virginia, USA.
Dare, P. and Dowman, I., 2001. An Improved Model for Automatic Feature-Based Registration of SAR and SPOT Images, ISPRS Journal of Photogrammetry and Remote Sensing, Vol. 56, No. 1, pp 13-28.
Hiep, Q. L., S. Gautama, and W. Philips, 2004.
Automatic Registration of Synthetic Aperture Radar (SAR) Images, Geoscience and Remote Sensing Symposium, (IGARSS'04), Vol. 6, pp. 3864-3867.
Hong T. D. and Schowengerdt, R. A., 2005. A Robust Technique for Precise Registration of Radar and Optical Satellite Images, Photogrammetry Engineering and Remote Sensing, Vol. 71, No. 5, pp585-593.
Lin, J. B., 1996. Interferometric Synthetic Aperture Radar Image Registration with Sub-pixel Accuracy, MEng thesis, Department of Electrical and Control Engineering, NCTU.
Liu, B. W., 2003. A Feature-Based Approach for Robust Registration of NASA/JPL AIRSAR and Aerial Ortho Imgae, MEng thesis, Department of Geomatics, NCKU.
Schenk, T., 1999, Digital Photogrammetry, Volume I, TerraScience, pp.251~266.
Wang, C. T. and K. S. Chen, 2003. A study of multi-temporal processing of SAR imagery, Journal of Photogrammetry and Remote Sensing, vol.8, no.1, pp.47-56.
Wang, C. P., 2004. Investigation on Image Matching by Variance-Covariance Components , MEng thesis, Department of Civil Engineering, NCU.
Automatic Registration of Satellite SAR Image Based on Multi-Resolution Matching
Hao-Lun Hung 1 Chih-Tien Wang2 Kun-Shan Chen3
ABSTRACT
Due to its all-weather capability and stable repeat observations, synthetic aperture radar is able to provide timely information after a natural disaster. In this study, we developed a multi-scale matching technique to efficiently and effectively integrate multi-temporal for fast and accurate change detection.
Experimental test results from RADARSAT and ALOS satellite radar images confirm that the proposed automatic registration scheme based on a multi-scale matching greatly reduces the artificial error and required computation, and thus is highly feasible for operational purpose.
Key Words: Synthetic Aperture Radar, automatic registration, multi-scale matching
Received Date: Sep. 14, 2007 Revised Date: Dec. 20, 2007 Accepted Date: Dec. 21, 2007 1 Graduate student, Institute of Space Science, National Central University
2 Associate research engineer, Center for Space and Remote Sensing Research, National Central University
3 Professor, Center for Space and Remote Sensing Research and Institute of Space Science, National Central University
