雷達干涉技術 - 應用多重感應器監測雲林地區三維變形

過去地層下陷的地表檢測方法是水準測量與GPS 測量綜合並行，以降低工作成本與時間，但上述兩種測量方法皆為點狀測量，而地層下陷主要為面狀下陷，因此以往進行地層下陷監測，必須佈設大面積的監測網與監測樁，進行地面測量後，再進行統計內插，以獲取面狀下陷資訊，此作法必須耗費相當大的人事成本與時間，因此取樣的頻率與密度常會受限於監測的經費，進而影響整體監測的準確度。

近年來隨著航太技術的快速發展，人們開始大量利用飛機及人造衛星等各式航空器酬載的儀器，達成許多從前不容易辦到的大範圍地表變形監測。本章節將介紹合成孔徑雷達差分干涉技術

（Differential Synthetic Aperture Radar Interferometry, 以下簡稱 DInSAR）原理、誤差來源、分析方法及應用於雲林地區地表變形監測之結果，藉以輔助及驗證GPS及其他大地測量方法。

5-1 雷達干涉技術原理 5-1-1 InSAR原理

InSAR(Synthetic Aperture Radar Interferometry, 合成孔徑雷達干涉)技術係利用不同時間或不同位置的雷達天線所獲得之兩幅或多幅複數影像做干涉計算，藉由相位差異獲取地表三維資訊。為能了解其原理，將InSAR 技術之簡單成像幾何表示如圖 5-1-1。其中， S1、S2為獲取地表資訊時天線的空間位置、B為兩天線之基線長、H為載體航高、 θ為視角、 R為第一個天線至目標物之距離、R+δR為第二個天線至目標物之距離、 α為基線與水平線之夾角、h為目標物高程。

圖5-1-1 InSAR成像幾何示意圖(修改自張中白等，2003)。

2 2

圖5-1-2 DInSAR成像幾何示意圖(修改自張中白等，2003)。

如圖5-1-2所示，由不同時間雷達天線 S1與S2所拍攝之兩張影像做干涉計算，即可獲得第一組干涉影像對。假設於兩張影像相距的時間內，任意點由P移動至 P'(P, P', S1 三點共線 )，其中， S1

至P之距離為 R1=R2+dR、 S2至P'之距離為 R2+δR、 δR為沿目標物與天線連線(S P'2 )方向上，地面位移所產生之位移向量。兩張影像經過套合，則可得到相對應之共軛點。根據干涉雷達量測理論，S1與S2兩影像間P與 P'之相位差 φ 為：

2 1 2 1

4 4

= - = π (R + R-R )= π (-dR+ R)

φ φ φ δ δ

λ λ (5-10)

式(5-10)中，φ 為相位回復後之相位差，由圖 5-1-2可知 (Zebker et al., 1994)：

dR Bsin ( - )≈ θ α (5-11)

dR B≈ ll (5-12)

由式(5-10) 可知，相位差 φ 是由地形效應產生之相位差 (-⁴π dR

λ )，與地表變形產生之相位差 ( ⁴π δR

λ ^⋅ )之總和，由式 (5-10) 可知：

R= ( +4 dR) 4

λ π

δ φ

π λ (5-13)

為消除由地形效應產生之相位差(-⁴π dR

dh=R sin

二軌跡差分干涉處理中，影響干涉處理品質之好壞包

括：利用干涉法產生之地形對必須進行相位回復，處理過

雷達訊號在發射與接收過程中都必須經過大氣層，大

果產生附加之高程誤差。此外，二軌跡差分干涉法中，若

若以低角度入射角接觸目標物表面，其回波強

3. 散射體 (Diffuse Reflector)

中，因幾何條件不同會產生成像結果的差異。

(iv) 相位變化值符合要求

干涉的條件中，鄰近像元間之相位變化不能超過一條干涉條紋(2π)，若超過此臨界點，像元間就會產生不相關。以Envisat衛星為例，進行偵測變形時，一個條紋代表變形2.8公分，若變形事件於像元解析度範圍內之雷達視距方向變化超過2.8公分，則像元間產生不相關且變形量無法於干涉圖中顯示。

圖5-2-1 Tuscon, Arizona地區大氣效應造成DInSAR相位變化範例

(Hoffmann, 2003)

(a) (b)

圖5-2-2 大氣壓力變化與溼度變化對不同波段之相位計算影響 由圖(a)及圖(b)得知，溼度變化對相位遲滯的影響量遠大於大氣壓力變化造成之影響量，且波長越短之波段，其訊號傳遞時間遲滯現象越顯著，

SIR-C/X-SAR 三種波段依序為：L(λ= 24 cm)、C(λ= 5.67 cm)、X(λ= 3.1 cm) (Zebker et al., 1997)

圖5-2-3 L波段相對濕度與大氣遲滯造成DInSAR誤差關係圖 L波段於相對溼度20%時，大氣遲滯對二軌跡法與三軌跡法造成約10公分與 14公分之雷達斜距變形誤差(Zebker et al., 1997)

圖5-2-4 垂直基線長與同調性參數關係圖 (Hooper, 2006)

圖5-2-5 雷達成像幾何變形示意圖

（a）側視雷達與地表目標物觀測幾何（b）前坡縮短（c）疊置（d）陰影 (陳卉瑄，2001)

0 Coherence 0.7

Short perpendicular baseline (176 m)

5 km

Longer perpendicular baseline (917 m)

圖5-2-6 地表植被分布與DInSAR結果關係圖

地表植被分布區域容易使雷達波散射，造成影像同調性降低且無法產生干涉條紋 (Hooper, 2006)

圖5-2-7 時間基線長與同調性參數關係圖 (Hooper, 2006)

Amelung, Jonsson et al. [2000] NASA, Visible Earth website

4 months 8 years

5 km

0 Coherence 0.4

5-3 雷達干涉技術處理流程 5-3-1 資料來源與使用軟體

本研究利用DInSAR技術，對台灣雲林地區地層下陷情況進行評估。其中，衛星資料取自歐洲太空總署所發射之Envisat衛星，並擷取軌道編號track 232, frame 3123之圖幅 (圖 5-3-1)，另外搭配40公尺空間解析度的數值高程模型 (Digital Elevation Model, DEM)進行差分干涉計算，且配合法國太空總署 (CNES)所研發之 Diapason V.4.1 版軟體進行處理。

圖5-3-1 Envisat衛星經過台灣地區之軌道路徑圖 紅色方框為本研究區 (track/frame, 232/3123)

Envisat衛星於西元 2002年 3月發射，主要目的為延續 ERS-1/2 之地球觀測任務。因此，同於ERS-1/2為一太陽同步資源衛星，飛行高度約為799.8公里，軌道傾斜角度約為 98.55度，重複觀測週期為35天。Envisat衛星上，共酬載 10個載具，其中包括一合成孔徑雷達系統，名為 ASAR (Advanced Synthetic Aperture Radar)。 ASAR為一多偏極雷達，可提供 HH & VV, HH & HV, VV

& VH資料，其衛星籌載波段為 C波段，波長為 5.6公分。

5-3-2 處理流程

本研究針對雲林地區，以Envisat衛星所擷取之雷達影像為主要資料來源，篩選出衛星軌道基線差較小的影像配對，進行 DInSAR計算。處理流程主要包含以下三大項 (圖 5-3-2)：

(1) 干涉影像對規劃。

研究區域決定後，影像對規劃與選取標準即成為初期工作最重要的關鍵。選取影像對的標準如下：

(a) 軌道基線標準：在影像對軌道基線(baseline)選取規劃上，垂直基線長小於400公尺，使其空間基線不相關造成之誤差較小。

(b) 影像對時間標準：根據經驗，將超過二年的影像對日期加以配對後，得到的影像干涉同調性結果較差，同時影像模糊且較難準確套疊。因此在規劃時，兩組(或多組)衛星資料日期最好能選擇在二年以下，使其時間基線不相關造成之誤差較小、以得到較佳的干涉成果。

(2) DInSAR計算。

本研究利用法國太空總署 (CNES) 所研發之 Diapason 差分干涉作業軟體，藉以調整及設定雷達影像參數值、雷達影像套合及影像差分干涉處理。

Diapason 軟體之影像處理流程大致可分為 7 個步驟 (圖 5-3-3)：

(a) 將主影像及副影像之原始回波資料經由SAR 處理子程式，分別處理成主影像及副影像的單觀點複數影像。

(b) 將ESA 所提供之 Envisat 衛星精密軌道資料 (Precise Orbit Data)加入程式中，以消除軌道所造成之誤差。 (c) 將主影像之單觀點影像及副影像之單觀點影像取樣

套合後，經干涉處理產生干涉圖。

(d) 將涵蓋雷達影像範圍之數值高程模型資料輸入程式中，產生模擬全相位圖像。

(e) 套合步驟(c)之干涉圖與步驟 (d)之模擬全相位圖像。 (f) 再次干涉處理，去除數值高程模型模擬而來的全相位

圖像，得到去除地形效應之差分干涉圖。 (g) 干涉相位濾波處理。

利用 Goldstein and Werner (1998) 所提出之 Goldstein radar interferogram filter處理雜訊。其主要概念為：將干涉圖像中相位變化較為平滑之區域進行較強的濾波處理以濾除雜訊；將相位變化亂度較高之區域進行較弱之濾波處理。這個方法目前公認對於干涉影像雜訊濾除十分合理且有效 (Goldstein and Werner, 1998)。

(3) 差分干涉圖像之相位回復。

利用美國 Stanford Radar Interferometry Research Group 所研發之 SNAPHU 解相位軟體進行相位回復 (Phase Unwrapping)(Chen and Zebker, 2002)。

圖5-3-2 DInSAR處理流程圖

圖5-3-3 Diapason軟體進行DInSAR處理流程圖

(Diapason, 2006)

5-4 應用DINSAR技術監測雲林地區下陷研究

本研究利用不同時間Envisat 衛星所拍攝雷達複數影像計算，得到之相位值差異以獲取雷達斜距方向之地表變形量

（line-of-sight displacements）。

在差分干涉的計算中，有許多因素會影響相位因而導致地表變形計算之誤判。其中，影響最顯著的為衛星軌道誤差及非重覆性的誤差(Massonnet and Feigl, 1998)。衛星軌道誤差會在差分干涉圖中產生明顯平行排列之干涉條紋，嚴重影響資料計算。因此，本研究採用ESA所提供之 Envisat精密軌道參數以消除此項誤差。而非重覆性的誤差(non-stationarity error)包含兩大類：第一類誤差為隨機的雜訊（random noise），沒有特定分佈，例如數值高程模型(Digital Elevation Model, DEM)、 SAR處理過程、影像本身及地表特徵物所造成的誤差；第二類非重覆性的誤差主要來自大氣的遲滯效應(atmospheric conditions)。雷達波雖然可以穿透大氣中的雲、霧及水氣，但卻可能使雷達波傳遞遲滯，因而造成計算時的誤差。Hoffmann (2003)指出，大氣遲滯效應所造成之干涉誤差一般難於差分干涉計算中消除，因而造成應用DInSAR監測地表變形之限制。

本研究選取2006年 9月 21日至 2007年 6月 28日之 8幅 Envisat影像，為降低前述差分干涉計算時可能造成之誤差，配對條件設定為垂直基線長小於400公尺，共計產生 13組影像對，影像對之規劃表及規劃圖如表5-4-1與圖 5-4-1。選取之 Envisat影像格式為： ASAR image single look complex中之 VV偏極雷達，且選取降模式及IS2 之攝像模式，其影像幅寬 105公里、空間解析度約 30 公尺、雷達入射角度為19.2~26.7度。

表5-4-1 雲林地區 Envisat影像干涉規劃表

Pair

No.

Master image (date)

Slave image (date)

︱B⊥︱

(m)

Day difference (days)

1 2007/03/15 2007/01/04 130 70 2 2007/04/19 2007/02/08 141 70 3 2007/04/19 2007/03/15 190 35 4 2007/02/08 2006/10/26 306 105 5 2006/10/26 2006/09/21 362 35 6 2007/03/15 2007/02/08 331 35 7 2007/04/19 2007/01/04 320 105 8 2007/05/24 2007/02/08 88 105 9 2007/05/24 2007/03/15 243 70 10 2007/05/24 2007/04/19 53 35 11 2007/06/28 2007/02/08 58 140 12 2007/06/28 2007/03/15 273 105 13 2007/06/28 2007/04/19 83 70

圖5-4-1 雲林地區Envisat影像(track/ frame, 232/3123)干涉規劃圖

資料處理部分，使用Diapason軟體以二軌跡法解算差分干涉

值而非實測值，因此DInSAR可以比水準測量測得更精確的下陷中心，同時展現更細微的變化，進一步凸顯DInSAR於監測地表變形之高空間解析的特性。

為進一步求取DInSAR及水準測量之差異量，因此在圖 5-4-2 與圖5-4-3中，以原雲林水準網上之水準點（雲林地區之水準網參考圖3-1-5），連結成 A-B（東西向）及 C-D（南北向）兩條剖面線，並以兩者量得之垂直變化量進行比對，經比對之成果顯示，AB剖面線之均方根差異量（ RMSE）為 1.2cm/yr， CD剖面線之均方根差異量（RMSE）為 1.1cm/yr。

由於DInSAR 及水準測量的測量時間有些微的差異，同時 DInSAR會受到大氣效應與植被的影響，造成 DInSAR計算之空間及時間不相關誤差，推測其為造成DInSAR及水準測量的成果產生差異的可能原因。但由本研究資料顯示，應用DInSAR技術探

在文檔中應用多重感應器監測雲林地區三維變形 (頁 119-152)