行政院國家科學委員會專題研究計畫成果期末報告
自然災害災區㆔維空間資料之快速蒐集建立應用暨災變前後變遷分析（I）-子
計畫㆓：資料快速獲取、處理、計算、與呈現之技術
Fast Acquisition of 3D-Data and –Information after Cataclysm in Disaster Area and its Application
on Change Detection and Analysis
期末報告(1/3)
衛星雷達測量臺灣九㆓㆒大㆞震災區災變前後高程變動之測試結果
Test Results of Satellite Radargrammetry for Determining Vertical Displacements after Chi-Chi
Earthquake
計畫編號：NSC 89-2218-E-006-126
執行期限：89 年 8 月 1 日至 90 年 10 月 31 日
主 持 ㆟：蔡展榮 研究助理：呂建興
執行機構及單位名稱：國立成功大學測量工程學系
㆒、㆗文摘要
本文使用涵蓋臺灣㆞區的㆓號歐洲遙測衛 星資料和衛星雷達測量法的差分干涉技術來求 定 臺 灣 九 ㆓ ㆒ 大 ㆞ 震 災 區 災 變 前 後 高 程 變 動 量。首先分析相關的理論內容，數學公式證明了 ㆒項特性：「如果合成孔徑雷達干涉法 INSAR 能求定出數公尺等級的數值高程模型精度，則其 差分法 D-INSAR 可以求定出數公釐等級的高程 移位精度。」其次，將實驗成果和高精度衛星定 位測量法求定的高程移位量做㆒比較，顯示兩者 平均較差約6cm，如果衛星雷達資料合宜而且後 續資料處理恰當，則使用歐洲遙測衛星雷達資料 和 D-INSAR 技術來量測臺灣都市㆞區的高程變 動極具潛力。關鍵詞
：雷達量測、數值移位量模型、合成
孔徑雷達差分干涉法、㆔像法
Abstract
In this paper, some second European remote sensing satellite (ERS-2) data in Taiwan area and the radargrammetry technology of differential interferometric synthetic aperture radar (D-INSAR) are utilized to determine vertical displacements after Chi-chi earthquake. Some basic theoretic aspects are firstly analyzed. It verifies that D-INSAR can provide a digital displacement model (DDM) with an accuracy of a few millimeters, if
INSAR can creates a DEM with an accuracy of a few meters. Comparing with the precise displacements determined using GPS shows that the mean difference of the displacement values estimated by D-INSAR and the ones determined by GPS is about 6cm. It indicates that three-pass method and ERS data processed properly can detect height surface deformation in urban area in Taiwan.
Keywords
:Radargrammetry, Digital Displacement Model, D-INSAR, Three-pass method㆓、緣由與目的
「主動式雷達量測系統」直接對探測區放射 出特定波長的電磁波，並且量測和記錄從探測區 回來的回波。在任何㆝氣㆘，只要載臺可以升空 飛行，無論是白㆝或晚㆖，任意時刻皆可以使用 雷 達 來 量 測 。 合 成 孔 徑 雷 達 SAR(synthetic aperture radar)是㆒種的雷達量測技術，透過合成 技術虛擬出㆒枝比實際的雷達㆝線更長的長㆝ 線，因此，即使選用長波長的雷達波，仍然可以 獲得方位方向（即載臺飛行方向）㆖的高解析力 雷達資料(Henderson & Lewis, 1998)。SAR 不僅獲 取 探 測 區 ㆞ 表 覆 蓋 物 的 雷 達 散 射 係 數(radar scattering coefficient)，也記錄探測物的相位(phase) 資 料 。 合 成 孔 徑 雷 達 差 分 干 涉 法 D-INSAR(differential interferometric SAR)是㆒種 的 SAR 技術，它使用相同探測區在不同時間㆘ 的差分干涉圖像(differential interferogram)來計算 探測區的物表面變形量，它可以求定（例如由㆞震或土石滑動所造成的）瞬間快速的高程變動、 與（例如㆞層㆘陷或緩慢的㆞殼移動所產生的） 慢速度的高程變動。而㆔像法(three-pass method) 是㆒種的 D-INSAR 技術，本文分析此法的計算 方法，並從理論㆖扼要㆞評估其成果之品質，再 提出臺灣都會㆞區的測試成果。
㆔、㆔像法
本文提出的測試成果皆使用美國 Vexcel 公 司的3DSAR SAR 處理系統的㆓個子系統來進行 資料處理與計算，這㆓個子系統分別名為FOCUS 和 PHASE，以㆘簡要㆞來闡述本文使用相關演 算法和參數值，詳細內容請參閱(Vexcel, 2000)。FOCUS 是美國 Vexcel 公司的 3DSAR 系統 裡的SAR 資料處理模組，它包含了㆔個子系統：
SAR 參 數 之 估 計 、 都 卜 勒 演 算 法(range Doppler algorithm) 、 CEOS (=Committee on Earth Observation Satellites)資料轉換。有名的都 卜勒演算法是用來處理常見的 CEOS L0 格式的 原 始 輸 入 資 料 。 SAR 參 數 估 計 運 用 MLBF(multi-look beat frequency) 和 MLCC(multi-look cross correlation)這兩種演算法 來 求 定 都 卜 勒 奇 異 數 (Doppler ambiguity number)，其㆗，MLBF 和 MLCC ㆓法分別最適 用於高反差區和低反差區，兩者分別使用在影像 的不同區域，最後挑選出最㆒致的奇異數值。都 卜勒估值器(field Doppler Estimator)使用了「相關 都 卜 勒 估 值 器 CDE(Correlation Doppler Estimator)」或「符號都卜勒估值器 SDE(Sign Doppler Estimator)」來求算都卜勒矩心(Doppler centroid) 的 奇 異 量 ， 而 以 距 離 相 關 法 (range correlation algorithm)來求定都卜勒速率(Doppler rate)。基本㆖，要求算的 SAR 參數是「都卜勒矩 心(Doppler centroid)」和「都卜勒速率(Doppler
rate)」，兩者用來決定 azimuth chirp function，此
㆒函數用於方位壓縮(azimuth compression)㆖。前 述結果用來執行兩個㆒維快速傅立葉變換 1D FFT(Fast Fourier Transform)，首先，由衛星㆖的 感 測 器 暨 波 束 性 質 來 決 定 range chirp coefficients，這些係數表達㆒個 range chirp 資 料，將這些資料和原始輸入資料的距離方向的線 資 料 相 乘 ， 依 此 方 式 來 進 行 距 離 壓 縮(range compression)。緊接著運用先前取得的 azimuth chirp coefficients 來 進 行 方 位 壓 縮 (azimuth compression)。 另 外 ， 也 可 進 行 距 離 曲 率 改 正 (range curvature correction)、輻射率定(radiometric calibration)、多視處理(multi-looking)、和㆞距㆖ 的再取樣(resampling to ground range)。最後，
FOCUS 產生各式標準的 CEOS L1 格式的輸出資 料，例如，本文使用的歐洲遙測衛星 ERS1/2 所 要的CEOS SLC(single-look complex image)單觀 點複數影像資料。 ㆔像法(three-pass method)使用㆔張 SLC 影 像（又稱為pass）做為主要的輸入資料來產生㆒ 個 數 值 位 移 量 模 型 DDM(digital displacement model)，因此，此法稱為「㆔像法」。它是 PHASE 系統裡的㆒個子系統，而 PHASE 是㆒個完整的 SAR 干涉處理系統。㆔像法假設位移是㆒種短時 瞬間快速的移位、或是㆒種隨著時間產生等速率 均勻變化的移位。在短時瞬間快速移位的情況 ㆘，㆔張SLC 影像的選擇暨編號順序、以及它們 和短時瞬間快速移位事件之時間先後關係應予 以適當定義，如圖1 所示，㆞震是㆒種典型的短 時瞬間快速移位事件之㆒例，它應發生於所選定 的參考影像(reference SLC image)和第㆔張影像 取得時刻之期間，而不是發生於參考影像和第㆓ 張影像取得時刻之期間。另外，將㆒張SLC 影像 的各像元資料乘以第㆓張 SLC 影像相應的同名 像元資料之共軛複數，所得各像元之新數據產生 所謂的「干涉圖像(interferogram)」。㆒張 SLC 影像裡的每㆒個像元均有其自己相應的衛星位 置，即取得該像元資料之時刻的衛星位置，因 此，干涉圖像裡的每㆒個像元皆有其相應的基 線，此㆒基線是指產生該干涉像元的兩個衛星空 間位置的連線，即兩個相應的SLC 影像像元對應 的兩衛星位置之連線。 ㆒般常用多視平均(multi-looks average)處理 來減少 SLC 影像的班駁雜訊(speckle noise)及相 位量測值雜訊量(Zebker et al. 1994b)，卻也同時 降低空間解析力及失去若干相位資訊。為了克服 前述之缺點，若干研究尚待進行，例如（王志添 等，2001）提出雷達影像多時處理(multi-temporal processing)之初步測試成果，證明多時平均處理 應可有效降低班駁雜訊，改善原有單張影像之品 質，同時也保有原影像之空間解析力。 圖 1. ㆞震事件、㆔張影像取得㆕者可能之時間 先後順序（修改自(Vexcel, 2000)） 在前述的㆔張影像時間關係㆘，由參考SLC 影像和第㆓張SLC 影像產生的干涉圖像稱為「參 考干涉圖像(reference interferogram)」，它反映了
探測區的㆞形資訊，另㆒方面，由參考SLC 影像 和第㆔張SLC 影像產生的干涉圖像稱為「第㆓張 干涉圖像(second interferogram)」，它反映了探測 區的原㆞形資訊與因（㆞震）事件引起的位移資 訊之和，因此，兩者相減可得純由（㆞震）事件 引起的位移所產生的相位數據，進㆒步計算處理 可得探測區位移量之數值模型DDM。 ㆔像法使用㆔張 SLC 影像產生㆒張探測區 的位移圖，其資料計算處理步驟依序如㆘： 1.對第㆓張 SLC 影像做再取樣(resample)，以取得 它在參考影像裡每㆒個相應像元位置㆖的新 影像資料。 2.使用參考影像和經再取樣處理產生的第㆓張影 像 來 計 算 取 得 參 考 干 涉 圖 像 (reference interferogram)。 3.過濾參考干涉圖像。 4. 對 過 濾 後 的 參 考 干 涉 圖 像 進 行 相 位 復 原 (unwrapping)處理。
5.使用㆞面控制點 GCP(ground control points)對 參考干涉圖像進行幾何改正。 6.對第㆔張 SLC 影像做再取樣(resample)，以取得 它在參考影像裡每㆒個相應像元位置㆖的新 影像資料。 7.使用參考 SLC 影像和經再取樣處理產生的第㆔ 張 SLC 影 像 來 計 算 取 得 第 ㆓ 張 干 涉 圖 像 (second interferogram)。 8.將產生的第㆓張干涉圖像減去參考干涉圖像的 復原相位，以得到移去㆞形影響分量後的「移 平干涉圖像(flattened interferogram)」。 9.過濾此時的移平干涉圖像。 10.對過濾後的移平干涉圖像進行相位復原之計 算。 11.修正第㆓張干涉圖像的基線長。 12.產生位移圖(displacement map)。 前述資料處理步驟1-7 和 9-10 請參閱（蔡展榮、 陳鴻緒，2001）。在步驟 8 裡，使用步驟 5 由㆞ 面控制點求定的精確基線長來改正步驟7 產生的 第㆓張干涉圖像，然後，再將它減去參考干涉圖 像的復原相位，以得到移去㆞形影響分量後的 「移平干涉圖像(flattened interferogram)」。這㆒ 張移平干涉圖像在步驟9 過濾後應無太大的改正 （更動），也應可以使步驟 10 的相位復原更易 處理成功。如果第㆓張干涉圖像之相位值不是採 用由㆞面控制點求定的精確基線長來改正，而是 採用基線長之粗略估值來改正，則可能會有「躍 進式」的相位粗差(residual phase ramp)存在於干 涉圖像裡，並且可能有㆒些殘留的㆞形訊號存在 於干涉圖像裡。 在步驟 11 裡，使用精確基線長來修正前述 步驟產生的干涉圖像之相位值，以得到修正後的 干涉圖像，這㆒張修正後的干涉圖像之計算只用 來 研 判 精 密 基 線 長 是 否 成 功 ㆞ 移 除 ㆞ 形 訊 號 (topographic signal)，後續的資料處理並不使用這 ㆒張修正後的干涉圖像。如果沒有完全移除㆞形 訊號，㆒般建議先剔除不良的㆞面控制點，然後 再重新計算精密基線長。㆒般建議不使用落在位 移大的區域裡的㆞面控制點。 在步驟12 裡，可供選用的位移模式有兩種： 短時瞬間快速移位(abrupt change)模式、和 隨 著時間產生等速率均勻變化的移位(continuous change)模式，前者 假設位移事件發生於參考 SLC 影像和第㆔張 SLC 影像取得時刻之期間， 而在參考 SLC 影像和第㆓張 SLC 影像取得時刻 期間不發生任何位移事件，而後者 假設在整個 ㆔張影像取得時刻期間裡，發生等速率位移。 簡言之，兩張干涉圖像是使用㆒張共同的參 考SLC 影像計算產生的，其㆗的參考干涉圖像記 錄著由㆞形訊號引起的相位值，第㆓張干涉圖像 記 錄 著 由 ㆞ 形 訊 號 和 ㆞ 面 移 動 引 起 的 相 位 之 和，因此，由㆞形產生的相位可由參考干涉圖像 計算而得，將第㆓張干涉圖像減去此㆒㆞形分量 就可以估算出由㆞面移動引起的相位值，也可以 進㆒步計算出數值移位量模型DDM。
㆕、誤差來源
到目前為止，雖然已有許多成功的應用案 例，但是，SAR 技術還不能完全實現㆞震學和㆞ 球科學領域所要求的高精度㆞表面移位訊息之 量測需求，其主要原因是肇因於㆘列各項輻射、 幾何、和計時誤差來源(Curlander & McDonough, 1991)： (1) 雷達系統： 計時控制、和頻率控制、 類 比式無線電波頻率RF(Radio Frequency)電子 設備、 數位電子與資料傳輸、 雷達㆝線。 (2) （衛載、或空載）載臺和資料㆘接系統： 波 段 誤 差(channel errors) 、 資 料 ㆘ 接 (downlink)速率降低、 資料壓縮、 浮點量 化處理(block floating point quantization)。 (3) ㆞面資料處理系統： 關聯器演算法和結構(correlator algorithms and architectures)、 輻 射改正、 幾何改正、 影像資料瀏覽。
(4) 模式誤差：各種理論假設和真實狀況之不㆒ 致性，例如，和其他的㆝文及太空大㆞技術
㆒樣，合成孔徑雷達干涉法INSAR 也受到電 磁波通過大氣層所產生的延遲效應之影響， 此㆒大氣延遲現象在不同的空間位置和不同 的時間都會改變，各種的大氣條件造成對流 層延遲(tropospheric delay)作用，進而導致雷 達波整體行進路徑之延遲，也增加了約公分 級的誤差量。另外，電離層擾動(ionospheric disturbance)又加入了額外的雷達波傳播延遲 (Zebker et al., 1997)。另㆒方面，干涉法有㆒ 個重要的假設，即：雷達回波(radar echo)的 相 位 值 必 須 具 有 可 重 覆 取 得 之 特 性 (reproducibility)，簡言之，如果在完全相同的 同㆒個位置再次量測雷達波相位值，則這些 量測值應該彼此相等。然而，事實㆖，這些 重新量測取得的相位值將會彼此不相等，所 以也會引起所有的INSAR 成果的誤差。這些 重覆量測的相位值之差異性稱為相位值的 「 降 相 關(decorrelation)」 ， 並 以 「 同 調 性 (coherence)」、或「相關性(correlation)」指標 來表示之。會造成相位降相關的主要因子是 時間降相關(temporal decorrelation)、 基線 降相關(baseline decorrelation)、和 斜視降相 關(squint decorrelation)。和雷達脈衝頻率值 (pulse repetition frequency)相較，當兩個不同 的 相 位 量 測 所 用 的 都 卜 勒 頻 率(Doppler frequencies)之差值不是很小的時候，此時的 斜視降相關將會相當顯著，也導致相位資料 問題。如果降相關性太顯著，則會導致相位 資料無法復原，這樣的相位資料就變成無用 的資料。
五、品質分析
假設在圖2 所示的㆒個簡易合成孔徑雷達量 測幾何情況㆘， 、 、 分別表示㆔張 SLC 影像，在取得 和 期間，㆒個㆞面點 P 移到 P′，再假設 、P、P′㆔點在同㆒條直線㆖，此時， 可得㆘列㆓式： 1 S S2 S3 1 S S2 1 S dh R B d ≈ − − ⋅θ
λ
α
θ
π
φ
sin ) cos( 4 1 (1)(
R d d δ λ π φ = 4 ⋅)
(2) 在本文使用的測試資料裡，雷達波長λ≈0.057m， 基線長B≈220m，斜距 ≈785000m，視角θ≈23°。 在使用臺灣㆞區的歐洲遙測衛星 ERS1/2 資料 ㆘，P 點㆖的相位誤差 dφ和高程誤差dh 之關係 曲線如圖3 左所示，而相位誤差 dφ和移位量求定 值誤差 d(δR)之關係曲線如圖 3 右所示。當α=0° 時，1m 的高程誤差將引起約-9°的相位誤差，而 等量的移位量求定值誤差 d(δR)=1m 將引起約 12631°的相位誤差，因此，如果合成孔徑雷達干 涉法INSAR 能產生㆒個數公尺(m)精度等級的數 值高程模型DEM 時，則 D-INSAR 求定出的相對 應之數值移位量模型DDM 的精度將可達數公釐 (mm)等級(Zebker et al. 1994a)！1 R 圖 2. ㆒個簡化的合成孔徑雷達干涉法之幾何， 其㆗ 【修改自(Zhou et al., 1998)】 R d R dR R R1= 2+ = 3+ ′ 圖3. 在使用臺灣㆞區的 ERS1/2 資料㆘，相位誤 差 dφ和高程誤差dh 之關係曲線（左：實線表示 α=0°之關係曲線）、相位誤差 dφ和移位量求定值 誤差d(δR)之關係曲線（右）
六、測試結果
本文的全部測試使用㆕張原始CEOS 格式的 ERS-2 資料，它們隸屬於 232 號軌道、frame 編 號3123，資料涵蓋臺灣㆗部㆞區，即民國 88 年 九月21 日凌晨 1:47 時發生的集集大㆞震主要災 區，資料取得時間分別是同㆒年（民國 88 年） ㆒月 21 日(1/21)、五月 6 日(5/6)、九月 23 日 (9/23)、十月 28 日(10/28)，測試區在台灣㆗部都 會區，此區離九㆓㆒集集大㆞震的震央不遠。參 考干涉圖像是採用㆞震前的 5/6、和 1/21 兩張 SLC 影像計算產生，第㆓張干涉圖像是分別採用 5/6、和 9/23 兩張 SLC 影像以及 5/6、和 10/28 兩 張SLC 影像計算產生的，㆞震發生於產生第㆓張 干涉圖像的兩張SLC 影像取得時刻之間，藉著這 些資料來計算位移量值。首先，依序執行”估算 粗略基線長”、”第㆓張 SLC 影像之再取樣”、” 計算產生干涉圖像”， 表 1 顯示㆔組 SLC 影像對 (SLC pairs)的㆒些基本計算成果，包括它們的垂 直基線長估值、水平基線長估值、以及套合結果 （連結點之點數、訊雜比 SNR 門檻值、方位與 距離方向㆖的套合誤差），圖4 顯示它們在各個 像元位置㆖的同調性值，計算成果顯示了東部山 區的同調性值低於 0.2，都會區的同調性值較 高，其值約 0.5~0.7，而西部海面的同調性值為 零。圖5 顯示了它們的干涉圖像，由於第㆔組影 像對（5/6、10/28）的同調性較第㆓組者（5/6、 9/23）高（請參閱圖 4），所以第㆔張干涉圖像 的干涉條紋比第㆓張者更清晰。因為前述兩組的 SLC 影像對至少皆包含著由集集大㆞震所引起 的降相關，所以第㆒組SLC 影像對（5/6、1/21） 所產生的參考干涉圖像具有最清楚的干涉條紋。 表1. ㆔組 SLC 影像對(SLC pairs)和它們的垂直 基線長估值、水平基線長估值、套合結果（連結 點之點數、訊雜比 SNR 門檻值、方位與距離方 向㆖的套合誤差） 套合結果 參考 SLC 影像 第㆓ 張 SLC 影像 垂直 基線 長估 值(m) 水平 基線 長估 值(m) 連結點 之點數 訊雜比 SNR 門 檻值 方位方 向套合 誤差(像 元) 距離方 向套合 誤差(像 元) 5 月 6 日 1 月 21 日 89 -89 570 11 0.127 0.101 5 月 6 日 9 月 23 日 -206 49 167 12 0.290 0.150 5 月 6 日 10 月 28 日 16 -41 107 12 0.212 0.071 圖 4. ㆔張同調性影像：5/6、1/21（左，垂直基 線長 ≅97m），5/6、9/23（㆗， ≅-220m）， 5/6、10/28（右， ≅5m） ⊥ B B⊥ ⊥ B 圖 5. ㆔張干涉圖像：5/6、1/21（左，垂直基線 長 ≅97m），5/6、9/23（㆗，B ≅-220m），5/6、 10/28（右， ≅5m） ⊥ B ⊥ ⊥ B 緊接著進行”過濾參考干涉圖像”、”復原參考 干涉圖像”、和”參考干涉圖像之幾何改正”，參考 干 涉 圖 像 幾 何 改 正 所 採 用 的 ㆞ 面 控 制 點 之 位 置，如圖6 所示。而圖 7 分別顯示了過濾後、相 位復原後、幾何改正後的參考干涉圖像。 然後，執行”第㆔張 SLC 影像之再取樣”、” 計算產生第㆓張干涉圖像”，兩張第㆓張干涉圖 像如圖5（㆗、右）所示。接著進行”第㆓張干涉 圖像之移平”處理，換言之，就是將第㆓張干涉 圖像扣除㆞貌引起的相位分量，此㆒部份的相位 值是從相位復原後的參考干涉圖像計算而得，由 此產生的兩張移平干涉圖像，如圖8 所示。此時， 再 繼 續 施 行” 移 平 干 涉 圖 像 之 過 濾 和 相 位 復 原”，由此產生的兩張復原後的干涉圖像，如圖 9 所示。接著對此時的第㆓張干涉圖像做基線長的 改正，圖 10 顯示這樣的幾何改正後的第㆓張干 涉圖像。 圖6. ㆞面控制點 GCP 之位置分佈，每㆒個星號 表示㆒個GCP
圖 7. 過濾後（左）、相位復原後（㆗）、幾何 改正後（右）的參考干涉圖像
圖 8. 扣 除 ㆞ 形 相 位 分 量 後 的 移 平 干 涉 圖 像 (flattened secondary inteferograms) ： 5/6 、 9/23 （左），5/6、10/28（右） 圖 9. 相位復原後的移平干涉圖像(unwrapped secondary inteferograms)：5/6、9/23（左），5/6、 10/28（右） 最 後 ， 計 算 產 生 位 移 量 圖(displacement map)。㆒個數值位移量模型(digital displacement model, DDM)儲存著㆞表面各處的位移值，這些 位移值表示沿著衛星的視方向（雷達㆝線相位㆗ 心到㆞面探測點之方向）㆖的位移向量之大小， 而且當此㆒位移向量指向衛星時，其值為正。隨 後採用和(Takeuchi et al., 2000)相同的方法，將這 些值（位移向量）轉換到探測區的局部高程系統 ㆖，再將這些轉換後的高程位移量來和快速靜態 GPS 觀測法(GPS in fast static survey mode)求定 的高程位移值做㆒比較。圖 11 分別顯示第㆒組 5/6、1/21、9/23（左㆒、左㆓）以及第㆓組 5/6、 1/21、10/28（右㆓、右㆒）以彩色影像表示的數 值位移量模型 DDM 以及和它相應的正射影像 (orthorectified power image)，其㆗，兩張彩色 DDM 影像裡的每㆒個色彩條紋循環間隔表示㆒ 個單位的3cm 位移量。 圖 10. 幾何改正後的干涉圖像(geometry-refined secondary inteferograms)：5/6、9/23（左），5/6、 10/28（右） 圖 11. 以彩色影像表示的數值位移量模型 DDM 以 及 和它相應的正射影像(orthorectified power image)：第㆒組 5/6、1/21、9/23（㆖左、㆖右）， 第㆓組5/6、1/21、10/28（㆘左、㆘右）。其㆗， 兩張彩色DDM 影像裡的每㆒個色彩條紋循環間 隔表示㆒個單位的3cm 位移量 將㆖列 D-INSAR 求定出的高程位移值和落 在求定出的DDM 範圍裡的 7 個㆞面 GPS 點施測 求定出的高程位移值做㆒比較，得出兩者平均差 值約6.2cm（第㆒組）、6.5cm（第㆓組），最大 差值均小於 18cm，這些測試成果顯示了：「如 果有合宜的ERS 資料、而且資料計算處理恰當， 則㆔像法和 ERS 資料可以應用在臺灣都會㆞區 來求定高程移位訊息。」所謂「合宜的資料」至 少指的是「兩張SLC 影像的同調性(coherence)不 能太低，降相關性(decorrelation)不能太大」，因 此，此㆒條件可以列為選購衛載雷達影像之依 據。
七、結論
1. 本文初步的測試實驗結果顯示，如果有合宜 的資料而且後續各項的資料處理恰當，則使用㆔像法和歐洲遙測衛星雷達資料求定出的 臺灣都會㆞區的數值位移量模型之精度具有 達到數公分等級之潛力。 2. 位移較大的區域裡的㆞面控制點不宜使用。 3. 為了提高 D-INSAR 求定的數值位移量模型 DDM 之品質，空載合成孔徑雷(airborne SAR, AIRSAR)達似乎可行。在臺灣，譬如山區和 密植被區㆒直是航空攝影測量作業困難之㆞ 區，AIRSAR 能在這些㆞區發揮多大的互補 功能？AIRSAR 能在這些㆞區求定精確 DEM 和 DDM 的區域大小如何？其可以達到的精 度又是多少？在㆒般㆞區，航測受限於「可 見」的條件，不能在能見度差的環境（例如： 黑夜、多雲霧、多懸浮粒之混濁空氣）㆘執 行空照，此時，AIRSAR 的成果品質能達到 那㆒種精度等級？其執行快速資料取得之效 率如何？…等等，這些細節內容和問題尚待 後續研究和測試驗證之。
參考文獻
1. 王志添、陳錕山、陳家堂、郭進民，2001， 雷達影像多時處理研究，第㆓十屆測量學術 及應用研討會論文集，第 239-244 頁，國立 ㆗央大學太空及遙測研究㆗心。 2. 蔡展榮、陳鴻緒，2001，使用合成孔徑雷達 技術和歐洲遙測衛星資料來產生台灣㆞區的 數值高程模型之實驗成果，八十九學年度期 末研究成果研討會論文集。3. Curlander, J.C., and McDonough, R.N., 1991. Synthetic Aperture Radar – Systems and Signal Processing. John Wiley & Sons, New York, pp. 249-499.
4. Henderson, F.M., and Lewis, A.J., 1998. Principles and Applications of Imaging Radar. In: Manual of Remote Sensing, 3rd edition. John
Wiley & Sons, New York, pp. 2-4.
5. Takeuchi, S.Y., Suga, Y., Oguro, Y., Chen, A.J., and Yonezawa, C., 2000. Verification of INSAR Capability for Disaster Monitoring – a Case Study on Chi-Chi earthquake in Taiwan. ACRS 2000, Vol. 2, pp. 738-743.
6. Vexcel, 2000. 3dSAR SAR Processing System, User’s Manual.
7. Zebker, H.A., Rosen, P.A., Goldstein, R.M., Werner, C.L., and Gabriel, A., 1994a. On the Derivation of Coseismic Displacement Fields
Using Differential Radar Interfermetry: the Landers Earthquake. Journal of Geophysical Research, Vol. 99, No. B10, pp. 19617-19634. 8. Zebker, H.A., Werner, C.L., Rosen, P., and
Hensley, S., 1994b. Accuracy of Topographic Maps Derived from ERS-1 Radar Interferometry. IEEE Transactions on Geoscience & Remote Sensing, Vol. 32, No. 4, pp. 823-836.
9. Zebker, H.A., Rosen, P.A., and Hensley, S., 1997. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps. Journal of Geophysical Research, Vol. 102, No. B4, pp. 7547-7563. 10. Zhou, Y.Q., Kless, R., Van Gendren, J.L., and
Li, D.R., 1998. Differential SAR Interferometry: Priinciples and applications. Spatial Information Science Technology and its Applications, pp. 227-236.
行政院國家科學委員會專題研究計畫成果期末報告
自然災害災區㆔維空間資料之快速蒐集建立應用暨災變前後變遷分析（I）-子
計畫㆓：資料快速獲取、處理、計算、與呈現之技術
Fast Acquisition of 3D-Data and –Information after Cataclysm in Disaster Area and its Application
on Change Detection and Analysis
期末報告(2/3)
使用合成孔徑雷達技術和歐洲遙測衛星資料來產生台灣㆞區的數值高程模型之
實驗成果
Test Results on DEM Generation in Taiwan by Using INSAR and ERS Data
計畫編號：NSC 89-2218-E-006-126
執行期限：89 年 8 月 1 日至 90 年 10 月 31 日
主 持 ㆟：蔡展榮 研究助理：陳鴻緒
執行機構及單位名稱：國立成功大學測量工程學系
㆒、㆗文摘要
本文分析合成孔徑雷達干涉法 INSAR 產生 數值高程模型 DEM 資料之演算法和其運用歐洲 遙測衛星資料在臺灣㆞區求定 DEM 之成果，初 步實驗成果顯示：「如果可用的（高同調性的） 雷達回波處理恰當，則在臺灣都市㆞區可以求定 ㆒個約23m 精度等級的 DEM 資料，而在較少植 被、或完全裸露㆞處可以求定約5m 精度等級的 DEM 資料。」關鍵詞
：合成孔徑雷達干涉法、數值高程模型、 前後縱排㆓衛星觀測模式、單觀點複數Abstract
This paper analyzes the algorithm of SAR interferometry for generating a DEM in Taiwan. Our preliminary test results show that the DEM determined using INSAR and ERS data has an accuracy of about 23m in average in urban area, and about 5m in bare areas without vegetation and buildings, if all data in areas with usable radar echoes (higher coherences) is processed properly.
Keywords
:INSAR, DEM, Tandem mode, SLC㆓、緣由與目的
如果和航空攝影測量法做㆒比較，雷達影像 測量法(imaging radar)可以視為㆒個新的遙感探 測技術，㆒般常將它使用的探測波之波段包含了 波長λ=1mm~1m 之範圍，而㆒般的水氣雲霧只對 波長λ<2cm 的雷達波具有顯著的影響，雨水對波 長λ>4cm 的雷達波的影響甚小。主動式雷達量測 系統裡的感測器放射出雷達波，傳向探測物表 面，並接收記錄與量測從探測物而來的雷達回波 響應值，它的測量作業在白㆝和夜間皆可實施， 較不受㆝氣的限制（除了惡劣㆝氣會限制飛機安 全飛行條件外）、和㆒㆝ 24 小時都能施測，此 ㆓特性是雷達的兩個著名的、也是最常被㆟提及 的特性(Henderson and Lewis, 1998)。本文探討在 臺灣㆞區使用合成孔徑雷達干涉法 INSAR 和㆒ 號暨㆓號歐洲遙測衛星 ERS1/2 資料來產生數值 高程模型DEM 的㆒些基礎內容和特性。㆔、
INSAR 產生 DEM 之演算法
本文的各項實驗計算均採用本系購自美國 Vexcel 公 司 的 ㆔ 維 合 成 孔 徑 雷 達 處 理 系 統 (3DSAR SAR processing system)的兩個子系統 FOCUS 與 PHASE 之部份計算模組來進行，其計 算步驟和採用的相關參數值簡述如后，詳細內容 請參閱(Vexcel, 2000)。 資料處理的第㆒步驟是使用 FOCUS 計算軟 體將輸入的原始資料做㆒格式轉換，原始的雷達 資料（相位歷史資料）是以 CEOS L0 格式來記 錄儲存的，將它轉換成標準的 CEOS L1 格式的 輸出資料。這㆒個合成孔徑雷達資料處理模組 FOCUS 的詳細分析請參閱（蔡展榮、呂建興， 2001）、(Vexcel, 2000)。在完成前述的資料格式變換處理之後，接著 將每㆒對（兩張）SLC(single-look complex)影像 依照㆘列七個步驟來進行各項計算處理以產生 數值高程模型DEM(digital elevation model)。 1. 估算基線長起始值：計算基線長之起始估值， 用 以 評 估 是 否 值 得 繼 續 進 行 後 續 的 計 算 處 理，只有基線長落在合理值範圍內，才會繼續 執行後續的各項計算。而此時的基線長估值僅 僅只用來評估選取的雷達資料之適用性，後續 的各項計算皆不採用此㆒粗略的基線長起始 估值。 2. 兩張 SLC 影像之套合與再取樣：計算第㆓張 SLC 影像套合(register)到參考 SLC 影像㆖的投 影(mapping)函數，找出參考影像裡的每㆒個像 元 在 第 ㆓ 張 影 像 ㆖ 的 對 應 位 置 並 做 再 取 樣 (resample) ， 以 取 得 取 樣 後 的 第 ㆓ 張 影 像 (resampled secondary SLC)。在套合計算裡，使 用大量的連結點(tie points)，這些連結點在兩 張影像㆖的位置是已知的，並採用仿射函數做 為兩影像之間的投影函數之近似，運用所有的 品 質 良 好 之 連 結 點 來 求 定 六 個 仿 射 函 數 參 數 ， 這 些 品 質 良 好 的 連 結 點 之 訊 雜 比 SNR(Signal-to-Noise Ratio)皆大於㆒個給定的 門檻值，同時，這些品質良好的連結點在水平 和 垂 直 兩 方 向 ㆖ 的 套 合 誤 差 均 方 根 值 RMS_x、RMS_y 滿足㆘列條件者才被視為是 ㆒個良好的套合： n x RMS _ < 0.03 pixels 且 n y RMS _ < 0.03 pixels 其㆗，n 是前述品質良好的連結點之點數。在 再取樣過程裡，對第㆓張 SLC 影像的傅立葉 大小函數(Fourier spectrum)做㆒過濾，只保留 該函數和參考 SLC 影像者兩者共有的分量函 數部分，在長基線或斜視角大的情況㆘，這樣 的頻譜過濾處理(spectral filtering)可以顯著改 善干涉圖像的品質。 3.計算干涉圖像：此時也對參考 SLC 影像做㆒過 濾，只保留它與再取樣後的第㆓張影像兩者共 有的頻譜分量部份，然後，將兩影像對應的兩 像元值 、 做 之乘法運算，其㆗， 是第㆓張影像像元值 的共軛複數，全部 的兩影像像元對(pixel pairs)完成前述的乘法 運算，其結果產生㆒張干涉圖像。另㆒方面， 此時可以用精密衛星星曆資料、或是直接由套 合係數來估算基線長，其㆗，衛星星曆包含了 某㆒時刻的衛星位置和速度。此㆒資料處理階 段也分別輸出參考影像和取樣後的第㆓張影 像 的 多 觀 點 能 量 影 像(multi-looked power images)、同調性影像(coherence image)、和圓 ㆞球相位(spherical earth phase)，圓㆞球相位必 須加㆖干涉圖像裡的相位值以得到完整的干 涉 相 位(full interferometric phase) 。 同 調 性 (coherence)是㆒個表示干涉相位品質的指標， 其值介於0 與 1 之間，同調性等於 1 表示其干 涉相位最精確，同調性等於0 表示干涉相位計 算完全失敗。 r
A
A
s * s rA
A
⋅
A
* sA
s 4. 過 濾 干 涉 圖 像 ： 本 文 的 全 部 測 試 皆 採 用 Goldstein 濾波器(Goldstein and Werner, 1997) 來過濾干涉圖像資料，它對干涉圖像裡的各個 小區塊分別進行㆓維傅立葉變換(2D Fourier transforms)，將每㆒個傅立葉大小 s (spectrum) 依㆘式產生另㆒個傅立葉大小新值s′： s′ = s|s| (1) α 其㆗，α是Goldstein 濾波器參數，其有用的設 定值範圍是α = 0.0-2.0。Goldstein 濾波器處理 快 速 相 位 變 化 的 探 測 區 是 最 快 速 有 效 的 (Vexcel, 2000)。 5.復原干涉圖像：將㆔角函數計算記錄的值域[-π, π] 的 干 涉 相 位 進 行 相 位 復 原 (phase unwrapping)，轉換成復原後的連續相位函數。 PHASE 軟體提供兩種可用的相位復原方法：Goldstein 法(Goldstein et al. 1988; Ghiglia & Pritt, 1998)、 QGBC (= Quality-Guided-Branch -Cuts)法(Flynn, 1995; Gliglia and Pritt, 1998)， 兩 者 都 是 所 謂 的 「 路 徑 追 蹤 復 原 演 算 法 (path-following unwrapping algorithms)」，在干 涉圖像裡找尋奇點(residue points)並建立阻隔 牆(branch cut)。Goldstein 法連接正負奇點來產 生阻隔牆，以阻止復原路徑通過阻隔牆，確保 全幅之相位復原不矛盾，並使得阻隔牆的長度 越短越好。相較於 Goldstein 法建立越短的阻 隔牆之特點，QGBC 法以阻隔牆連接正負奇 點，使得阻隔牆落於干涉相位最不精確可靠 處 。 ㆒ 般 而 言 ，Goldstein 法處理速度遠比 QGBC 法快，但是，在干涉相位不確定處， Goldstein 法可能產生復原錯誤之結果。在低同 調性區域裡，經常會出現相位復原錯誤，為降 低此㆒影響，㆒般建議先使用 Goldstein 法進 行相位復原，如果在興趣區裡出現顯著的復原 錯誤時，可再使用QGBC 法，以期獲取最大的 成功相位復原區域。在前述的兩種路徑追蹤復
原演算法裡，可以選擇高同調性的像元做為起 始像元來啟動相位復原處理，同時，同調性值 低於門檻值的全部像元均不做相位復原之處 理，前述的同調性門檻值常定義在 0.0~0.2 的 範圍裡。 6.使用㆞面控制點進行相關的幾何改正：使用細 心 挑 選 出 的 ㆒ 群 ㆞ 面 控 制 點 GCP (ground control points)資料來計算參考影像相應的精 確衛星軌道（星曆）參數，並改正參考影像和 第㆓張影像兩影像之間的基線長。此時的精密 衛星星曆之計算是供輸出成果的精確正射糾 正(ortho-rectification)之用，而精確基線長之求 定是供後續求定精確高程值之計算用。使用㆒ 些輔助資料（例如㆞形圖、參考影像），經適 當方法可以將影像位置坐標(x,y)轉換到㆞圖 位置坐標(X,Y,Z)，反之亦然。求定的數值高程 模型DEM 之品質取決於選用的㆞面控制點的 精度和控制點之個數，㆒般而言，最少應使用 10 個㆞面控制點，使用越多個㆞面控制點，則 干涉式合成孔徑雷達 INSAR 技術可求定越精 確的數值高程模型DEM。 7.產生高程圖(height map)：此㆒資料處理階段將 輸出㆒個數值高程模型DEM 和㆒張正射糾正 後的合成孔徑雷達 SAR 影像，其㆗，可以選 擇指定的高程基準和水平基準以產生所要的 ㆒張高程圖(height map)。另外，可以設定 X 和 Y 兩方向㆖的高程點間距。本文研究過程 ㆗，也自行撰寫程式來將PHASE 產生的高程 圖檔轉換成 ASCII 格式的數值高程模型資料 檔(X,Y,Z)，俾供後續的 DEM 品質評估分析之 用。
㆕、實驗成果
如 果 與 歐 洲 遙 測 衛 星 ERS1/2 和 加 拿 大 RADARSAR 衛星提供的 C 波段(λ≈5.7、5.6cm) 雷達資料做㆒比較，日本JERS 衛星提供的 L 波 段(λ≈23cm)雷達資料之雷達頻率較低，波長較 長，所以探測物表面的時間變化產生的「時間降 相關(temporal decorrelation)」對 JERS 資料的影 響較小，換言之，㆔者之㆗，JERS 資料最不易 受時間降相關之影響，它是提供衛載雷達干涉量 測的理想資料；但是，由於本文研究期間尚無 JERS 資料可用，所以本文全部的實驗測試皆採 用ERS 資料，而且是「前後縱排㆓衛星觀測模式 (tandem mode)」資料，它們是 1995 年 12 月 8 日原始CEOS 格式的 ERS1 資料（track 461、orbit 22994）、和 1995 年 12 月 9 日原始 CEOS 格式的ERS2 資料（track 461、orbit 3321），其垂直 和水平基線長分別大約是-205m、78m。使用 ERS 資料欲求定好成果，垂直基線長應該介於25m 和 300m 之間(Vexcel, 2000)。如果基線太短，則對 ㆞貌起伏變化之感測將不靈敏。如果基線太長， 則基線降相關(baseline decorrelation)之現象將變 得顯著，相位資料將因此變得低同調性終致成為 無用的資料。圖 1 顯示了參考 SLC 影像(ERS2, 1995/12/09)的能量影像(power image)、和實驗區 的多觀點平均影像(multi-looked image)。 影像套合結果顯示在方位方向和距離方向 ㆖的套合誤差分別是0.14 和 0.09 像元，共有 550 個連結點用來求定將第㆓張 SLC 影像(ERS1, 1995/12/08) 套 合 到 參 考 SLC 影 像 (ERS2, 1995/12/09)的六個仿射轉換套合參數㆖，為了能 夠產生㆒張良好的干涉圖像，這㆒項影像套合誤 差必須要小於1 個像元，而且連結點之點數必須 要多於100 點。如果套合誤差大、或者連結點數 量少，則訊雜比門檻值(SNR threshold)應該要調 整；降低訊雜比門檻值將可取得更多的連結點， 增 大 訊 雜 比 門 檻 值 將 可 取 得 品 質 更 好 的 連 結 點，但是此時的連結點數量將減少。 當PHASE 軟體對第㆓張 SLC 影像做再取樣 時，它同時過濾水平（距離）方向和垂直（沿衛 星飛行）方向㆖的傅立葉頻譜，使得兩張 SLC 影像具有最高同調性，再用這樣高同調性的兩張 SLC 影像來產生干涉圖像。本文使用的 SLC 影 像對經 PHASE 軟體過濾處理後，在距離和沿衛 星飛行方向㆖分別僅保有頻寬裡 58%和 60%的 共有頻譜，如果在距離或沿衛星飛行方向的任㆒ 方向㆖保留的共有頻譜比率低於50%，則表示此 ㆒對的 SLC 影像將無法產生㆒張良好的干涉圖 像。沿衛星飛行方向㆖保留的共有頻譜比率是取 決 於 兩 影 像 取 得 時 刻 的 都 卜 勒 頻 率(Doppler frequencies)之差異大小，而距離方向㆖保留的共 有頻譜比率是取決於基線長。如果取得兩影像時 刻的都卜勒頻率不同，則沿衛星飛行方向㆖保留 的共有頻譜比率低。如果基線太長，則距離方向 ㆖保留的共有頻譜比率低。在本文的實驗計算 裡，使用 2 個距離方向㆖的像元和 10 個方位方 向㆖共 20 個像元取平均，再根據這樣解析力的 兩影像產生干涉圖像。圖 2 顯示了本文實驗影像 的同調性影像和干涉圖像。在同調性影像裡，低 同調性的像元用藍色表示之，較高的同調性之像 元以紫紅色、黃色表示之，而最高同調性之像元 以綠色表示之。在同調性影像裡，同調性值是表 示3 x 3 個像元窗形區平滑後的同調性值。當產 生干涉圖像後，使用套合成果提供的訊息來重新
估算基線長，此時求定出的基線長是比先前估算 基線所得的粗略基線長更加精確。 圖1. 參考 SLC 影像(power image)（㆖），矩形 區是實驗區，此區的多觀點平均影像顯示如㆘圖 圖 2. 同調性影像（左：”藍⇒紫紅⇒黃⇒綠”表示 其同調性”低⇒高”）和干涉圖像（右） 圖2 清楚顯示都會㆞區的同調性值高，山區 的同調性低，而海面區域的同調性值為零；另 外，將同調性影像和干涉圖像做㆒比較可知，同 調性越高，則該處的干涉條紋越清晰。 在進行相位復原之前，先對干涉圖像做平滑 處理以濾除㆒些相位雜訊，俾讓相位復原變得容 易些。在此過濾處理㆗，使用 Goldstein 濾波器 並且採用其參數內定值α=1，此㆒處理將很快完 成。經過這樣的過濾處理產生的干涉圖像如圖 3 所示，將它和未過濾前的干涉圖像（圖2 右）做 ㆒比較，顯示過濾後的干涉相位較原本的干涉相 位平滑多了。 緊接著採用 Goldstein 方法來復原相位，其 ㆗的同調性門檻值設定為 0.1，讓相位復原處理 不進入同調性低於 0.1 的區域裡（例如海面區域 和部份山區）。選擇㆒個㆕周都有良好同調性的 種子像元(seed point)，從這㆒個選定的種子像元 開始進行相位復原之計算。計算結果顯示復原後 的相位和先前過濾後的干涉圖像幾乎相同，只有 ㆒些洞形區因其外圍屬低同調區致使相位復原 無法跨越此區，所以洞形區內部較高同調性處仍 未做相位復原之處理。實驗區裡的復原相位值以 影像表示之，如圖4 左所示，由於選定的種子像 元位於實驗區（臺北都會區）內、而且此區外圍 均是低同調區，所以復原相位的範圍只有涵蓋實 驗區而已，也就是說它只涵蓋圖2 和圖 3 的整張 干涉圖像的北方㆒部份而已。 圖3. 使用 Goldstein 濾波器(α=1)過濾後的干涉圖 像 圖 4. 實驗區裡的復原相位（左）、復原相位影 像和能量影像之套合（右） 在計算數值高程模型DEM 之前，先使用 84 個㆞面控制點來估算㆒個較先前的”計算干涉圖 像”步驟由套合資訊所估算出的基線長還要更可 靠精確的新基線長，並且運用這些控制點來求定 ㆒個加常數，用來對復原相位做進㆒步的改正。 ㆒般而言，12~20 個㆞面控制點是足夠的，每㆒ 個控制點的平面坐標(X, Y)是用來修正衛星的星 曆資料，以供後續各項輸出成果的正射糾正之 用；這些控制點坐標應該要有數十公尺或更佳的 精度，否則衛星星曆資料精度無法提高。 本文採用比例尺1/10000 的實驗區正射影像 (ortho image)和同區域的 5m x 5m 間隔的數值高 程模型 DEM，這兩項資料是由平均攝影比例尺
1/5000 的航空攝影像片資料經航測方法產生 的，平均航高離㆞面約1600m，全部像片的攝影 時間從1997 年 10 月到 2000 年 1 月，也就是本 文實驗用的tandem mode 資料取得日期 1995 年 12 月 8-9 日之後的約 2-4 年期間。由於沒有和實 驗資料取得時間接近的較佳資料可供使用，所以 ㆞面控制點是由這兩項資料取得，並假設在那 2-4 年期間沒有發生顯著的㆞表變形。全部 84 個 ㆞面控制點的位置分佈如圖5 所示，很明顯㆞， 東部區域無任何的㆞面控制點。 在產生數值高程模型 DEM 的過程㆗，投影 方法設定是世界橫麥卡脫投影 UTM，水平面和 高程基準設定為WGS84 系統。在使用全部的控 制點資料之前，先將它們從原來的 1997 年臺灣 基準(Taiwan Datum 1997)TWD97 坐標系統轉換 成WGS84 系統，然後，使用這 84 個控制點計算 DEM，將計算出的 DEM 以高程圖(height map) 表示之，如圖5 右所示。 圖5. 全部 84 個㆞面控制點之位置圖（㆖）、求 定出之高程圖（㆘） 從高程圖可以清楚得知，計算得出的東部山 區DEM 資料是錯的，此區的高程不應低於 0m。 探究其因，此區多山潮濕而且覆蓋緻密植被，是 ㆒ 個 衛 載 合 成 孔 徑 雷 達 干 涉 法(spaceborne INSAR)求定 DEM 作業的困難㆞區，前述的航測 5m x 5m 數值高程模型 DEM 只涵蓋了如圖 7 所 示之範圍，在此區域裡，INSAR 求定的 DEM 之 某㆒個高程格點若落在航測DEM 範圍內，則在 以該點為㆗心的㆒個50 x 50m 窗形區內，採用所 有的航測 DEM 高程格點資料，以距離平方倒數 做為權值，進行加權平均，計算出該點處對應的 航測高程值，由此取得了航測DEM 和 INSAR 求 定之 DEM 兩者共有 93850 個共同點，航測和 INSAR 求定之 DEM 的高程值分別以 、H 表 示之，兩者的高程差dH = - 之平均值 是–3.70m，而其均方根值 RMS(dH)約 23m，圖 6 顯示高程差 dH 之直方圖和其位置分佈圖，實驗 區內部藍色矩形區是航測 DEM 無資料區，所以 無法和衛載INSAR 成果做㆒比較。 AP H AP H SAR SAR H 在這些共同的高程點之㆗，|dH|≤20m 者佔 86%，而且這些點幾乎都位在都市區域裡，圖 7 顯示它們的分佈圖和航測正射影像，將兩者套合 比較可以清楚得知幾乎全部的|dH|≤5m 點皆落於 無緻密植被和無建物之近似裸露㆞面區㆖，其他 的5m < |dH| ≤ 20m 點顯然落於緻密高建築物區 內，將dH 資料以㆔維面表示之，呈現㆒個粗糙 的㆔維城市表面模型。 另㆒方面，將高程圖（圖 5）和航測正射影 像（圖 7）做㆒比較，顯示在高樓圍繞的河川淺 灘㆞面或河面㆖，衛載INSAR 也計算出高程值， 此㆒現象可由此區的雷達多重散/反射效應造 成，當然是求定的DEM 之㆒項誤差來源。 圖6. 航測和衛載 INSAR 求定之高程的差值 dH = - 之直方圖（㆖）和位置分佈圖（㆘） SAR H HAP
圖7. 實驗區裡小於 20m 的高程差|dH|≤20m 分佈 圖（㆖）、和同區的航測正射影像（㆘）
五、結論
臺 灣 是 ㆒ 個 位 於 ㆞ 球 的 亞 熱 帶 ㆞ 區 的 島 嶼，潮濕多山多植被或是開發密集的都會市鎮區 域，本文初步實驗結果顯示，如果能取得合宜（高 同調性）的雷達資料、並且資料處理恰當，則使 用衛載合成孔徑雷達干涉 INSAR 技術和歐洲遙 測衛星 ERS 資料求定位於北臺灣的實驗㆞區的 數值高程模型DEM，其精度約為 23m（都會區）、 或5m（無植被和建築物之近似裸露㆞面）。 臺灣山區覆蓋緻密植被，針對此區使用衛載 INSAR 技術和 tandem mode ERS 資料來求定 DEM 之作業而言，還是㆒個困難區域，對實務 應用需求而言，相關的詳細研究仍是必要的。在 未來幾年裡，我們將直接研究空載合成孔徑雷達 干涉技術，以尋求在航測作業困難㆞區求定精確 DEM 的可行方法，以擴大台灣㆞區高精度基本 資料涵蓋全島面積之比率。參考文獻
1. 蔡展榮、呂建興，2001，衛星雷達測量臺灣 九㆓㆒大㆞震災區災變前後高程變動之測試 結果，八十九學年度期末研究成果研討會論 文集。2. Flynn, T.J., 1995. Consistent 2-D phase unwrapping guided by a quality map, Proceedings of the 1996 International Geoscience and Remote Sensing Symposium, Lincoln, NE, May27-31, IEEE, Piscataway, NJ, p2057.
3. Ghiglia, D.C., and Pritt, M.C., 1998. Two Dimensional Phase Unwrapping, John Wiley & Sons, New York, p103, p136.
4. Goldstein, R., and Werner, C., 1997. Radar Ice Motion Interferometry. 3rd ERS Symposium, Florence Italy.
5. Goldstein, R.M., Zebker, H.A., and Werner, C.L., 1988. Satellite radar interferometry: two-dimensional phase unwrapping, Radio Science, Vol. 23, No. 4, p713.
6. Henderson, F.M., and Lewis, A.J., 1998. Principles and Applications of Imaging Radar. John Wiley & Sons, New York, pp. 2-4.
7. Vexcel, 2000. 3dSAR SAR Processing System, User’s Manual.
行政院國家科學委員會專題研究計畫成果期末報告
自然災害災區㆔維空間資料之快速蒐集建立應用暨災變前後變遷分析（I）-子
計畫㆓：資料快速獲取、處理、計算、與呈現之技術
Fast Acquisition of 3D-Data and –Information after Cataclysm in Disaster Area and its Application
on Change Detection and Analysis
期末報告(3/3)
以
CSG 模型半自動化萃取建物㆔維資訊之研究
A Study on Semi-Automated Building Extraction Based on CSG Modeling
計畫編號：NSC 89-2218-E-006-126
執行期限：89 年 8 月 1 日至 90 年 10 月 31 日
主 持 ㆟：蔡展榮 研究助理：王聖鐸
執行機構及單位名稱：國立成功大學測量工程學系
㆒、㆗文摘要
以預先建立的模型萃取影像㆖建物的㆔維 空間資訊，稱為模型式建物萃取，是㆒種精確而 可靠的萃取建物㆔維空間資訊的方法。本文首先 依據建構實體幾何(Constructive Solid Geometry, CSG)的理論，參考台南㆞區常見的建物類型，設 計出數種參數化的基本元件(Primitive)。萃取建 物㆔維資訊時，則是分別以適當的基本元件套合 建物的某㆒部分，再將所有基本元件依布林運算 以建構實體幾何樹(CSG Tree)的概念加以組合， 即可表達㆒棟完整的建物模型。本研究採取半自 動化量測之策略，將高階的判斷工作由操作員執 行，如：選擇元件、粗略定位、賦予屬性等，而 將低階的計算工作交由電腦運算，如：模型與影 像之精密套合，以求得建物模型之形狀與姿態參 數。此㆒方法可大幅減輕操作員逐點進行立體量 測之負擔，且所得成果以物件實體為單位，易於 導 入 ㆔ 維 空 間 資 訊 系 統 或 其 他 後 續 之 應 用 處 理。經由實驗證明，本法將可有效改善目前航測 製圖的方式，快速而完整㆞獲取並展現建物的㆔ 維資訊。關鍵詞
：建構實體幾何、㆔維城市模型、建物萃 取、模型-影像套合Abstract
Building extraction based on pre-established models has been considered a promising approach to acquiring precise, reliable and complete 3D data of buildings from aerial images. This paper
proposes a semi-automated approach to extracting buildings based on CSG (Constructive Solid Geometry) modeling. In this approach, complex buildings are modeled as a combination of volumetric primitives. The concept of the approach is to let high-level tasks (building detection, model selection, and attribution) be carried out interactively by the operator and optimal model-image fitting be performed automatically by a computer algorithm. The shape and pose parameters associated to a primitive provide a link between perception (images) and prior knowledge (primitive) of a building part, so that the fitting method proceeds to determine the shape and pose parameters so as to fit a primitive with the corresponding images. Having all of the building parts been uniquely represented by parametric primitives, a building can be reconstructed by using CSG Boolean set operators to combine the building parts. This approach is developed with the prospects of releasing the operator from tedious point measurement and efficiently delivering precise and reliable results. In contrast with the traditional point-by-point digitization mapping process, this approach promotes an object-by-object data acquisition procedure.
Keywords
: Constructive Solid Geometry (CSG), 3D City Model, Building Extraction, Model-Image Fitting㆓、緣由與目的
㆔維城市空間資訊系統可以廣泛應用在都 市計畫、市政管理、觀光導覽等，整合、分析並
提供各方面所需的空間資訊，因此如何快速獲取 建物的㆔維空間資訊已成為測量界的重點研究 方向。應用各種全自動或半自動的方法，從航空 相片㆗快速萃取出房屋的㆔維空間資訊，則同時 受 到 航 測 領 域 與 計 算 機 視 覺 領 域 學 者 的 重 視 [Braun, et al., 1995, Grün, 2000, Lang and Förstner, 1996, Mohan and Nevatia, 1989, Shufelt, 1999, van den Heuvel, 2000, Vosselman, 1999]。儘管建物萃取的操作方式會 隨著影像資料類型與比例尺、物件複雜度、所需 詳細程度以及產品類型而異，㆒般的處理程序可 以概分為㆘列㆔個步驟：偵測(Detection)、重 建 (Reconstruction) 以 及 賦 予 屬 性 (Attribution)。目前已經有許多萃取建物的方 法被提出，其㆗㆒些特製的全自動化方法雖然可 以 解 決 某 些 特 殊 的 情 形 [Brenner, 2000, Förstner, 1999, Henricsson, et al., 1996, Hsiao and Wong, 1999]，然而由於自動化影像 判釋仍處於相當粗糙的階段，因此要發展㆒套能 應付各種狀況的全自動化建物萃取方法仍有困 難。基於㆔維空間資訊的精確性、可靠性及完整 性的考量，科學家們於是轉而尋求開發半自動化 的方法[Chio and Wang, 1999, Gülch, et al., 1998]。
模型式建物萃取[Ameri, 2000, Braun, et al., 1995, Brenner and Haala, 1999, Fischer, et al., 1999, Lowe, 1991, Suveg and Vosselman, 2000a, van den Heuvel, 2000, Vosselman, 1998] 的假設是以建物模型來描述㆔維空間㆗的目標 建物，設法在模型與現有影像資料之間建立關 係。大部分的模型式建物萃取多採取半自動化的 程序，其㆗高階的工作像是：建物偵測、模型選 擇、以及賦予屬性等都由操作員以㆟機互動方式 完成，而電腦系統則是用來執行低階的計算工作 如：模型-影像的最佳套合。當模型-影像套合完 成後，建物的空間資訊也就隨之確定，這種㆒個 物件㆒個物件萃取的方式將可以大幅減輕操作 員過去逐點量測的沈重負擔，有效㆞提供精確而 可靠的成果。 模型式建物萃取包括兩大議題：建物模型設 計與模型-影像套合。建物模型設計的議題在於 如何完整建立足以代表各式建物的模型，基於建 物外型的多樣性與複雜性，要為所有建物建立其 專屬模型幾乎是不可能的，較為可行的方法是採 用 CSG 的概念來設計建物模型。CSG 並非針對整 棟建物設計模型，而是依照建物所需詳細程度， 組合數個稱為基本元件的簡單形狀實體（Solid） 構 成 ㆒ 棟 建 物 的 複 合 模 型 [van den Heuvel,
2000]。模型-影像套合的重點則在設計㆒套電腦 程式演算法，自動調整模型的形狀與姿態參數， 使得模型框架投影在影像㆖時可以理想㆞重合 於所對應的邊緣線像元。在這裡我們假設影像的 內、外方位資訊均為已知，且模型的形狀與姿態 參數近似值在最佳套合之前已透過㆟機互動方 式決定。Sester 與 Förstner[1989]曾發展㆒套 叢集演算法(Clustering Algorithm)以強鈍㆞估 計形狀與姿態參數，然而該方法僅適用於單張影 像而非多張影像。Lowe[1991]提出㆒套最小㆓乘 演算法以進行模型-影像的套合，然而其主要目 標是求取投影㆗心參數，而非求取模型參數，亦 即影像的方位參數是未知而待定的。經由修改 Lowe 的最小㆓乘法，本文將詳細㆞描述㆒套針對 模型式建物萃取所發展的模型-影像套合演算 法。
㆔、相關研究
在過去十年間，建物萃取已成為㆒項廣泛研 究的議題，因此已經有許多有用的參考文獻收集 與回顧[Förstner and plümer, 1997, Grün and Dan, 1997]。在參考了各種不同的方法之後，可 以發現建物萃取必須以各式各樣的建物模型為 之，因為當㆔維空間的物件投影到㆓維空間時， 將會失去㆒些建物萃取所需的資訊。尤其是當影 像㆗有用的資訊被不相關的資訊所干擾時，如： 植披、車輛、建物的細部等，使用建物模型便顯 得十分重要，因此模型式萃取會因為所採用的建 物模型以及所採取策略的不同而有所不同。依照 所 採 用 的 模 型 ， 可 以 概 略 區 分 為 ： 多 面 體 (Polyhedral Models) 、 稜 柱 體 (Prismatic Models) 、 參 數 化 多 面 體 (Parameterized Polyhedral Models) 、 以 及 CSG 模 型 (CSG Models)。然而建物重建的策略就難以作很清楚 的分類，以全自動或半自動或許可以將建物重建 的程序概略分為兩類，但是獲取建物的外型並將 其匹配到影像㆖的方法則是相當多樣性。 3.1 建物模型 3.1.1 多面體模型 多面體模型是以邊界線來表達㆒個物件，因 此在模擬建物的外型㆖最具有彈性[Brenner, 1999, Grün, 2000, Henricsson, et al., 1996]，然而由於缺乏建物外型的語意類別，因 此 無 法 以 ㆒ 個 通 用 的 多 面 體 模 型 來 作 建 物 萃 取。㆒方面多面體模型的彈性可能造成重建程序 的不穩定，另㆒方面使用多邊體模型也很難預估 建物被遮蔽的部分。加入幾何或物件本身的約制，如：平行、共面、重直等，或許可以解決部 分㆖述問題，但是這又會引出如何設定約制條件 及加權的新問題。 3.1.2 稜柱體模型 稜柱體是由傳統㆓維製圖概念衍伸而來，允 許建物在㆞面的任㆒處表面㆖重建，但是建物的 牆 壁 必 須 垂 直 且 建 物 的 屋 頂 必 須 為 平 面 [Hendrickx, et al., 1997, Lang, et al., 1995, Mohan and Nevatia, 1989]。稜柱體可以視為多 面體的㆒種特例，儘管建物的外型受到進㆒步限 制，卻仍然沒有進㆒步分類，也沒有特定的參數 足以描述，因此其缺點和㆒般多面體模型㆒樣。 3.1.3 參數化多面體模型 參數化多面體模型是以㆒組參數預先定義 某㆒特定形狀的多面體，藉由指定該組參數的值 來描述物件的外型[Vosselman and Veldhuis, 1999]。這些預先定義的多面體可以是㆒些常見 的建物外型，例如：長方體、屋脊型、L 型、T 型等。每㆒個預先定義的多面體都具有㆒些形狀 參數以調整其長、寬、高，同樣㆞也具有㆒些姿 態參數以描述其位置與方位。參數化多面體隱含 了物件本身的約制條件，因此被遮蔽的建物部分 仍然可以被完整重建出來，而重建好的建物模型 也等於經過分類。而其缺點則是對於多樣的建物 外型缺乏彈性，無法針對各式各樣的建物建立㆒ 套完整的模型資料庫。 3.1.4 建構實體幾何模型 CSG 模型不像參數化多面體模型㆒次建立整 棟建物，而是依布林運算組合㆒組基本元件來表 達 ㆒ 個 複 合 式 的 物 件 [Braun, et al., 1995, Gülch, Mueller, Läbe and Ragia, 1998, Lang and Förstner, 1996, Veldhuis, 1998]。它可 以依照所需的建物詳細程度，以少數幾個基本元 件來表達㆒棟複合式的建物。基本元件是㆒種預 先定義的簡單形狀實體，帶有㆒些轉換參數，可 以調整其比例尺、旋轉及平移，以表達建物的空 間幾何性質，而基本元件之間的結合則是以布林 運算的結合、交集、以及差集來運作。因此㆒棟 複合式的建物可以 CSG 樹來表達，其㆗每㆒個節 點 (Node) 具 有 兩 段 分 支 用 以 代 表 所 結 合 的 元 件，而其末稍的樹葉則是所使用的基本元件。由 於 CSG 模型不僅可以彈性㆞表達建物的外型，同 時也隱含了物件內的約制條件與類別，因此許多 研究均採取 CSG 模型作為模型式建物萃取的基 礎。 3.2 全自動與半自動方法之比較 自動化㆒直是電腦科技發展的目標之㆒，但 是如果無法提供可靠的成果，那自動化就不具實 用的價值。全自動化的建物萃取自然是數位航測 ㆗㆒項極致的目標，在過去十年㆗也有許多方法 被提出[Brenner, 2000, Förstner and Gülch, 1999]，但是都只適用於特定的案例，要開發㆒ 套通用的全自動化萃取程序仍然有許多方面的 困難，其困難源自於目前目前自動化的影像判釋 仍停留在相當初級的階段。由於建物萃取的工作 會因為影像資料類型、比例尺、物件複雜度、所 需詳細程度與成品類型而不同，因此特製的的全 自 動 化 萃 取 程 序 在 不 同 的 狀 況 ㆘ 往 往 無 法 成 功。也因此半自動化的方法越來越受到注意，以 求取更精確、可靠而完整的㆔維空間資訊[Chio and Wang, 1999, Gülch, et al., 1998, Lang, Löcherbach and Schickler, 1995, Rottenstener, 2000]。 影像判釋需要含有高階資訊的知識庫與低 階資料的影像才能進行，經由處理影像來求取高 階資訊，稱之為：由㆘而㆖(Bottom-up)或資料 取向(Data-driven)的方法；反之，經由高階資 訊萃取特徵物並驗證萃取出特徵物與資料之間 的㆒致性，則稱之為：由㆖而㆘(Top-down)的方 法。通常影像判釋會同時包含由㆖而㆘以及由㆘ 而㆖兩種方法，而高階資訊與低階資料會在某㆒ 處理階段會合，以驗證其是否㆒致。㆒般而言， 全自動化建物萃取相較於半自動的方法，多半在 比較高階的階段驗證資訊與資料的㆒致性，以假 說測試或知識工程來決定最適合景象的建物模 型，然而迄今尚未有㆒套通用於建物萃取的假說 測試程序，因此全自動化的方法常因環境的不同 而失敗。因此半自動化的方法是由㆟來決定使用 何種模型，讓電腦不需處理高階的判釋工作，而 自動化㆞完成模型-影像套合的計算。由於㆟類 執行高階判釋工作的成果遠較電腦計算成果要 快且可靠，而低階的最佳套合計算則是由電腦執 行較有效率，因此半自動化的策略將比全自動化 的方法更實用。
㆕、半自動化建物萃取策略
本研究採用「由㆖而㆘」的作法，參考台南 ㆞區常見的建物類型，預先建立基本元件的模型 資料庫，例如：長方體、直角㆔角柱體、圓柱體、 圓錐體等等。每㆒個基本元件都具有特定的形狀參數(Shape Parameters)以描述其外觀，如圖 ㆒ ； 同 時 也 具 有 ㆒ 組 共 通 的 姿 態 參 數 (Pose Parameters)以描述基本元件在物空間㆗的位置 與姿態。 圖㆒、模型資料庫㆗的基本元件：由左而右依序 是(a)長方體、(b)直角㆔角柱體、(c)圓柱體、 (d)圓錐體 操作員在量測時只需要挑出適合的元件，拖 拉到影像㆖建物的對應部分，手動調整元件使其 近似對齊後，交由電腦以最小㆓乘模型-影像套 合演算法自動進行最佳套合即可。若套合成功， 則此時元件的形狀參數即代表某部分建物之外 型，而姿態參數則描述了該部分建物在物空間的 姿態與位置。若建物的外型相當複雜而無法以單 ㆒房屋元件表達時，則重複㆖述的步驟將㆘㆒個 元件拉到視窗套合，直到該房屋的所有元件都重 建完畢，再根據 CSG 的概念將所有元件依布林運 算(Boolean Set Operation)組合為㆒個房屋模 型。整個半自動化房屋萃取流程可以概分為：模 型選取、近似套合、最佳套合以及元件組合等㆕ 個步驟，如圖㆓所示： 操作員 選擇模型 ㆟機互動式 近似套合 電腦計算 模型-影像 最佳套合 布林運算與 元件修改 模型庫 圖㆓、半自動化房屋㆔維資訊萃取之流程 由圖㆓可看出，「半自動化」的分工原則在 於將高階的判釋工作，如：選擇元件、給予元件 初始大小、位置等，交由㆟工判斷；而低階的計 算工作，如：模型-影像套合，則是交給電腦自 動計算。此㆒分工的好處是㆟工判釋提高了元件 類型的正確率，同時給予套合計算良好的初始 值，相較於全自動化萃取法更為可靠；而電腦處 理低階計算工作可以減輕㆟力的負荷，提升房屋 萃取的效率。在組合元件的過程㆗，由於無法避 免的隨機誤差，元件與元件之間很可能出現不符 合現實狀況的縫隙，因此在布林運算時，也同時 提供㆟工介入修改元件的互動式介面以確保建 物模型的正確性。 利用本方法所萃取的建物資訊不再是逐點 數化後連線，而是逐棟將模型與房屋影像套合得 來，因此不再有兩棟房屋間共界線無法切割的問 題，所得房屋資訊更適合㆔維展示與㆔維城市空 間資訊系統之建置。
