Volume19, No3, January 2015, pp. 171-187
1國立成功大學測量及空間資訊學系 碩士 收到日期:民國 103 年 03 月 01 日
2國立成功大學測量及空間資訊學系 副教授 修改日期:民國 103 年 05 月 07 日
3國立成功大學測量及空間資訊學系 碩士 接受日期:民國 103 年 05 月 15 日
＊通訊作者, 電話: 06-2370876 ext.838, E-mail: tsayjr@mail.ncku.edu.tw
改良型 PS-InSAR 測量法求定臺灣中部地區的 地層下陷量
楊佳祥
1蔡展榮
2*蘇柏宗
3摘 要
本文提出改良型永久散射體雷達干涉測量法(PS-InSAR)，使用臺灣中部地區的 ALOS PALSAR 雷達 影像來進行 PS-InSAR 測量，用來求定此區的地層下陷量，此區涵蓋大範圍的植被區和山區，找出此區可 用的 200 個精密水準點資料，與 PS-InSAR 計算成果進行比較，兩者求定的高程變動速度差值之絕對值，
平均值為 1.4cm/year，最小值為 0.0cm/year，最大值為 4.5cm/year，均方根值 RMSD 為 1.4cm/year，考量 其顯著性下，兩者測量結果可視為相當吻合。實驗結果驗證了 PS-InSAR 配合 ALOS PALSAR 雷達影像在 測量臺灣中部地區的地層下陷量之應用潛力，ALOS PALSAR 目前已經除役，未來採用品質更好的雷達 影像如 TerraSAR-X 或 TanDEM-X，將具有提高 PS-InSAR 成果品質的良好潛力。
關鍵詞：合成孔徑雷達、永久散射體雷達干涉測量、地層下陷
1. 前言
根 據 臺 灣 地 層 下 陷 資 料 庫 (http://www.subsidence.org.tw/index2.aspx) 歷 年 來 的監測資料顯示，臺灣中部一直有地層下陷的現象，
尤其以彰化和雲林地區較為嚴重。地下水位和土層 含水量降低，導致整體的土層結構變得較為鬆散，
因土層重量往下壓而呈現大範圍或局部的地表下 陷。經濟部水利署(2011)的報告顯示，造成地層下 陷的原因除了超抽地下水之外，還包含地層結構及 豐枯水期的影響。臺灣的豐水期發生在每年的 5 月到 10 月，這期間包含梅雨與颱風季節，整年度 的降雨量幾乎集中在這 6 個月中；另外六個月則為 枯水期，降雨量只占年度總降雨量的一小部分。理 論上，豐水期間雨量豐沛，土層含水量上升，會減 緩地層下陷的速度；反之，枯水期的地層下陷速度 會增快。
如 地 層 下 陷 資 料 庫
(http://www.subsidence.org.tw/theorem.htm) 的 內 容
所述，到目前為止，全區域整體地層下陷現況之調 查是以精密水準測量為最普遍且精度較優的方法。
比較歷年的檢測結果，亦可觀察出地層下陷範圍之 變化及其與產業分佈或土地利用之關係；水準測量 時間較冗長、耗費人力、不易自動化，且測量之時 間間隔較大，不易獲得連續之地層下陷資料。目前 水利單位於沿海地層下陷較嚴重之地區施行傳統 水準測量，考量所需的人力經費和時間之限制，故 平均 2~3 年方得輪流施測一次，因此一般精密水準 測量較難個別評估豐枯水期的地層下陷，且水準點 的數量和密度並不足以反映整個大範圍監測區的 地表下陷。
採用衛載雷達資料配合永久散射體雷達干涉 測量(Permanent Scatterers Interferometric Synthetic Aperture Radar, PS-InSAR)可克服上述問題並得到 大範圍的地表下陷速度。相較地面測量作業，
PS-InSAR 的優勢有：1.主動式的雷達衛星可以不 分日夜進行地面資料收集，2.所發出和接收的訊號 為波長較長的微波波段，因此訊號可以穿透雲霧直
接對地表進行觀測，此外，3.衛載雷達資料可以在 較短的時間間隔內取得大範圍的地表面下陷資 料。
本 文 採 用 日 本 的 ALOS (Advanced Land Observation Satellite) PALSAR (Phased Array type L-band Synthetic Aperture Radar)衛載雷達影像(本 文簡稱為 ALOS 影像)，並以 PS-InSAR 計算影像 取得期間的地表下陷速度。第一個計算例(第 4 節) 是探討使用 ALOS 影像進行大範圍監測地表下陷 的可行性，測試區面積約 33 km×54 km，涵蓋地表 下陷嚴重的彰化和雲林地區，計算結果顯示大範圍 密集分布 PS 點的下陷速度和地層下陷資料庫中的 地真資料具有一致性，顯示未來可使用雷達資料來 監測地表下陷。第二個計算例(第 5 節)挑選地表下 陷較嚴重的虎尾地區來進行 PS-InSAR 測量，面積 約 12 km×12 km，分別計算其豐枯水期的地表下陷，
結果顯示豐水期的地表下陷速度確實較枯水期來 得趨緩，顯示 PS-InSAR 測量臺灣地區地表下陷的 應用潛力，可以突破地面測量的限制，可以測量大 面積地表下陷速度場且大幅縮短測量時間間隔。
2. PS-InSAR 計算方法
本 文 提 出 的 改 良 型 PS-InSAR 測 量 法 是 在 Mora et al.(2003)的方法中加入一個有利於挑選永 久散射體候選點(Permanent Scatterer Candidate，本 文簡稱為 PSC)的機制，同時結合垂直基線門檻值 和同調性門檻值來挑選 PSC，透過 Delaunay 三角 網將 PSCs(Permanent Scatterer Candidates)連結起 來 (Delaunay,1934; Costantini and Rosen,1999) ， PSCs 的位置分布取決於地形效應、地表變動效應、
大氣效應和訊號雜訊等因素，利用一個線性觀測數 學模型來擬合 PSCs 的相位結構，接著透過最小二 乘平差法來解算此線性觀測數學模型，最後萃取出 永久散射體 PSs(Permanent Scatterers)並得到其變 動速度。本文使用敝團隊楊佳祥先生改寫國立台灣 大學 胡植慶教授領導的研究團隊開發的 Matlab 程式來進行改良型 PS-InSAR 計算。
2.1 挑選 PSC
當兩張 SAR 影像套合完成後，一對共軛點的 訊號可分別表示成(Lee and Pottier,2009)
1 1
1
x j y
P
(1)2 2
2
x j y
P
(2) 其中j
代表 1；x
1、x
2、y1、y2為實數。以(1) 為 例 ， 訊 號 的 強 度 可 表 示 為 強 度 值 (Intensity) (Lillesand et al.,2004)2 1 2
1
y
x
(3) 或振幅值(Amplitude) (Ulaby et al.,1988)2 1 2
1
y
x
(4)，相位值可表示為(Santitamnont,1998)
1 1 1
tan x
y
(5)根據 Rosen et al.(2000)，(1)式和(2)式兩共軛點的同 調性值可表示為
2 2*
* 1 1
* 2 1
P P E P P E
P P E C
(6)其中𝑃1∗ 和𝑃2∗分別為 P1 和 P2的共軛複數；𝐸[ ]為 期望值。同調性值的值域從 0 到 1，可反應出目標 區塊的散射穩定度，同調性越高的目標區塊代表其 散射特性較穩定，並可在干涉計算中提供較可靠的 相位資訊。
長時間維持高同調性的 PSC 有較大的機率被 選為 PS 並能在 PS-InSAR 的計算過程中被萃取出 來。Mora et al.(2003)的方法結合同調性影像和同調 性門檻值，以簡單、迅速和較少的 SAR 影像對來 挑選出足夠可用的 PSCs，其經驗顯示只需 7 對 SAR 影像對即可得到不錯的計算結果，然而當 SAR 影像對少於 5 對時，則成果品質不佳；Ferretti et al.(2001) 採 用 振 幅 散 佈 指 標 (Amplitude
Dispersion Index)挑選出影像強度較強的目標區塊 為 PSC，此方法需要使用至少大於 30 張的 SAR 影 像才具有統計意義(Mora et al., 2003)；Zhang et al.(2011)利用雷達散射波在空間上的偏移特性來 挑 選 出 散 射 特 性 較 穩 定 的 目 標 區 塊 ， 稱 為 Temporarily Coherent Point (TCP)，模擬資料的計算 結果顯示，只要兩張 SAR 影像即可獲得數量足夠 且可靠的 TCPs 來進行 PS-InSAR 的計算。
SAR 影像對在垂直基線較短的條件下，較能 避免受到干涉幾何條件的影響而降低本身的同調 性(Ferretti et al.,2000; Mora et al.,2003)。改良型 PS-InSAR 延伸自 Mora et al.(2003)所提出的方法，
適合應用在計算條件不佳的區域，其本身只採用垂 直基線低於某門檻值的同調性影像來挑選 PSC，當 某目標區塊在所有同調性影像上的同調性值皆大 於或等於同調性門檻值時，則被視為 PSC，接著不 考慮垂直基線門檻條件，將所有 PSCs 在差分干涉 圖上的相位資料全部拿來做 PS-InSAR 的計算。和 Mora et al.(2003) 所 提 出 的 方 法 不 同 ， 改 良 型 PS-InSAR 並不是採納同調性平均值，而是當所有 同調性值皆需大於特定門檻值時才能被挑選為 PSC，此方法較符合 PS 須能長時間維持散射穩定 的特性，避免因某次特定極高的同調性值而拉高整 體同調性平均值，導致挑選出不可靠的 PSC。
2.2 萃取 PSs 並計算其變動速度
本文主要參考 Mora et al.(2003)所提出的方法，
以下公式推導包含(7)式到(18)式主要參考並修改 自(Mora et al.,2003)，使用最小二乘平差法來進行 PS-InSAR 的解算，且只計算地表線性變動，暫不 考慮地表突發性、快速性和不規則的變動。以 Delaunay 三角網法將所有 PSCs 連結起來，某兩連 結 PSCs 干涉相位值的差值可表示為
noise atm
mov topo
flat
int
(7)其中等號右邊為𝛿𝜙int的五個主要組成因子，𝛿𝜙𝑓𝑙𝑎𝑡
代表目標區塊到 SAR 載具之間地面距離的影響；
𝛿𝜙𝑡𝑜𝑝𝑜代表地形效應的影響；𝛿𝜙𝑚𝑜𝑣代表斜距方向 上的變動效應影響；𝛿𝜙𝑎𝑡𝑚代表大氣效應的影響；
𝛿𝜙𝑛𝑜𝑖𝑠𝑒為訊號雜訊所引起的相位因子。𝛿𝜙𝑓𝑙𝑎𝑡、
𝛿𝜙𝑡𝑜𝑝𝑜和𝛿𝜙𝑚𝑜𝑣可分別表示為
tan
4 B r
flat
r
(8)
sin
4 B h
topo
r
(9)nonlinear nonlinear
linear mov
T
v
4
(10)
其中
為雷達訊號的波長；r代表斜距；B
為空 間基線長，代表兩張 SAR 影像取像位置之間的空 間距離； r
為兩連結 PSCs 之間r的增量值；
為 觀測視角，代表垂線方向和斜距方向的夾角； h
和
v
分別代表兩連結 PSCs 之間高程的增量和斜 距變動速度的增量；𝛿𝜙𝑙𝑖𝑛𝑒𝑎𝑟和𝛿𝜙𝑛𝑜𝑛𝑙𝑖𝑛𝑒𝑎𝑟分別為 𝛿𝜙𝑚𝑜𝑣的線性和非線性分量；T 代表時間基線長。當𝛿𝜙𝑓𝑙𝑎𝑡和𝛿𝜙𝑡𝑜𝑝𝑜在差分干涉計算中被移除後 (Zebker et al., 1994; Madsen and Zebker, 1998)，兩 連結 PSCs 的差分干涉相位值的差值可由(7)式轉換 成
noise atm
mov topoerror
dif
(11)(11)式中的𝛿𝜙𝑡𝑜𝑝𝑜𝑒𝑟𝑟𝑜𝑟代表差分干涉計算中移除地 形效應時所殘留的地形誤差相位因子，可將其表示 為
sin
4
B
topoerror
r
(12) 其中
代表兩連結 PSCs 的地形誤差增量值。大 氣條件會隨著 SAR 影像取像的時間不同而改變，𝛿𝜙𝑎𝑡𝑚可視為隨機變量(Hanssen, 1998)，且𝛿𝜙𝑎𝑡𝑚和
𝛿𝜙𝑛𝑜𝑖𝑠𝑒彼此獨立不相關。
(11)式可進一步詳細表示成
n xx yy ii n xx yy ii
i y x i y x
y x y
r x B
y x v y x v T i
y x y x
n n dif m m dif
n n dif m m dif
n n dif m m dif
n n dif m m dif i i
i
n n dif m m dif i n
n m m dif
, , ,
,
, , ,
,
, , ,
,
, sin ,
4
, 4 ,
, , , ,
(13)
其中(xm,ym)和(xn,yn)分別代表兩連結 PSCs 的局部
平面坐標；
i
代表差分干涉圖的編號；vdif代表線性 的斜距變動速度；𝜀𝑑𝑖𝑓為差分干涉計算所產生的地 形誤差量；假設𝑣𝑑𝑖𝑓和𝜀𝑑𝑖𝑓在每個 SAR 影像對皆為 固定值；𝛽𝑑𝑖𝑓、𝛼𝑑𝑖𝑓和𝑛𝑑𝑖𝑓分別代表非線性的斜距 變動量、大氣效應和訊號雜訊所引起的相位因子。根據 Hanssen(1998)所提出的理論，水平距離 1km 以內的兩連結 PSCs 擁有相似的大氣條件，可 將兩者的大氣差異視為空間上的低頻訊號。將水平 距離超過 1km 的兩連結 PSCs 之間的相位關係剔除，
接著兩連結 PSCs 之間彼此的大氣條件視為相似且 在(13)式中可表示成
x
my
mi
dif x
ny
ni
dif
, , , ,
(14) 將(14)式代入(13)式、移除非線性斜距變動速度的 影響並將相位雜訊因子視為最小二乘平差計算中 的殘差，最後可得到線性觀測數學模型，即(15) 式。其中λ、𝑇𝑖、𝑟𝑖、𝐵𝑖和𝜃𝑖視為已知值；𝛿𝜙𝑀𝑜𝑑𝑒𝑙為 觀測值；𝑣𝑀𝑜𝑑𝑒𝑙和𝜀𝑀𝑜𝑑𝑒𝑙則為待求的未知數，可採 用最小二乘平差法解算得到兩者在兩連結 PSCs 之 間增 量值 的最 佳估 值。 模 型同 調性 函數 (Model Coherence Function，本文簡稱為 MCF)可表示為(16) 式(Ferretti et al.,2000)。
其中𝑁𝑑為差分干涉圖的數量，上述最小二乘平差
法的解算條件就是要最大化 MCF 值(值域從 0 到 1)，每對連結的 PSCs 都可得到一個代表平差結果 品質的 MCF 值，MCF 值等於 1 代表最完美的平差 結果；而 0 則代表平差計算得到完全錯誤的結果。
待平差計算完成後，先剔除兩 PSCs 之間 MCF 值小於 0.7 的連結，接著再剔除無任何其它連結的 獨立 PSC，萃取剩下的 PSCs 成為 PSs。MCF 門檻 值 0.7 參考自 Mora et al.(2003)所使用的數值，可根 據不同的計算條件做調整，提高 MCF 門檻值雖然 可以萃取出品質較高的 PSs，但 PSs 的數目會因此 減少而較難正確評估整體計算區的變動；相反地，
降低 MCF 門檻值則可以提升 PSs 的數目，不過卻 也降低了計算結果的可靠性，建議在計算條件良好 的區域如建物密集的都市區，可提高 MCF 門檻值 來獲得品質較高的 PSs，卻仍然可保有密度足夠的 PSs。
兩連結 PSs 之間𝑣𝑀𝑜𝑑𝑒𝑙的增量值可表示為(17) 式，其中𝛾𝑚𝑎𝑥代表 MCF 值已達最大化。利用 Xu and Cumming(1999)所提出的一種相位復原區域成長 法可計算得到每個 PS 的𝑣𝑀𝑜𝑑𝑒𝑙，先以 MCF 值較高 的 PSs 為起始資料來源的種子點，接著依序計算出 每個 PSs 的𝑣𝑀𝑜𝑑𝑒𝑙，其計算式可表示成(18)式。
Model m m Model n n
i i
i
n n Model m
m Model i n
n m m Model
y x y
r x B
y x v y x v T i
y x y x
, sin ,
4
, 4 ,
, , , ,
(15)
Ndi Model m m n n
n n m m dif d
n n m m Model
i y x y x
i y x y x y N
x y x
1
, , , ,
, , , 1 ,
1 exp ,
,
,
(16)
,
max, ,
,
,
m n n Model m m Model n n m
estimated
x y x y v x y v x y
v
(17)
p
p p Model p
p estimated p
p estimated
p
p p Model estimated
y x y x y
x y x v
y x v
y x y x y
x v
, , , ,
, , ,
, , , , 1
( 1 8 )
其中
x, y
代表目前所處理的 PS 的局部平面坐標；p
代表與目前處理的 PS 有所連結的其它 PSs；𝛾𝑀𝑜𝑑𝑒𝑙(𝑥, 𝑦, 𝑥𝑝,𝑦𝑝)可 視 為 用 來 計 算 𝑣𝑒𝑠𝑡𝑖𝑚𝑎𝑒𝑑(𝑥, 𝑦) 的權，前者越大則對於後者的貢獻也越多。
3. 測試資料
表 1 列出本文計算所採用的 20 張 ALOS 影像 的資訊，資料取得時間為 2006 年 12 月 31 日至 2011 年 2 月 26 日止，原始影像的像幅約為 60 km×70 km(方位方向×距離方向) ，主要涵蓋地表下陷較嚴 重的彰化和雲林地區，方位方向（Azimuth Direction）
表 示 雷 達 天 線 移 動 的 方 向 ； 距 離 方 向 （ Range Direction）則與方位方向垂直。影像資料模式包含 Fine Beam Single Polarization (FBS)和 Fine Beam Dual Polarization (FBD)，FBS 提供單極化的 HH 或 VV 影像，FBD 則提供雙極化的 HH+HV 或 VV+VH 影像，本文測試資料的極化方式統一採用 HH。在 極化方式的表示中，H 和 V 分別表示水平和垂直 的極化方向，以 HV 為例，前者 H 代表感測器發 射水平極化的訊號脈衝波；後者 V 則代表感測器 接收垂直極化的訊號脈衝波。本文根據不同計算案 例的需求，個別挑選所需要的影像資料。
表 1 測試所用的 ALOS 影像資訊 影像編號 衛星 感測器 取樣時間
(年/月/日) 模式 極化 格式等級 軌道編號 像幅編號 1 ALOS PALSAR 2006/12/31 FBS HH 1.1 446 460 2 ALOS PALSAR 2007/02/15 FBS HH 1.1 446 460 3 ALOS PALSAR 2007/08/18 FBD HH 1.1 446 460 4 ALOS PALSAR 2007/10/03 FBD HH 1.1 446 460 5 ALOS PALSAR 2008/01/03 FBS HH 1.1 446 460 6 ALOS PALSAR 2008/05/20 FBD HH 1.1 446 460 7 ALOS PALSAR 2008/07/05 FBD HH 1.1 446 460 8 ALOS PALSAR 2008/08/20 FBD HH 1.1 446 460 9 ALOS PALSAR 2008/10/05 FBD HH 1.1 446 460 10 ALOS PALSAR 2009/01/05 FBS HH 1.1 446 460 11 ALOS PALSAR 2009/07/08 FBD HH 1.1 446 460 12 ALOS PALSAR 2009/08/23 FBD HH 1.1 446 460 13 ALOS PALSAR 2009/10/08 FBD HH 1.1 446 460 14 ALOS PALSAR 2010/01/08 FBS HH 1.1 446 460 15 ALOS PALSAR 2010/02/23 FBS HH 1.1 446 460 16 ALOS PALSAR 2010/07/11 FBD HH 1.1 446 460 17 ALOS PALSAR 2010/08/26 FBD HH 1.1 446 460 18 ALOS PALSAR 2010/10/11 FBD HH 1.1 446 460 19 ALOS PALSAR 2010/11/26 FBD HH 1.1 446 460 20 ALOS PALSAR 2011/02/26 FBS HH 1.1 446 460
圖 1 顯示本文測試案例的計算區域，並將 ALOS 影像套疊到 Google Earth 上，所有 ALOS 影 像皆透過多觀點處理，將像元的地面尺寸統一轉換 為 31m×31m，因此單一 PS 點的變動速度代表對應 地面範圍 31m×31m 在斜距(slant range, 本文簡稱 SR)方向上的平均變動速度 VSR。為了簡化計算，
假設每一個 PS 點上的雷達入射角等於視角(look angle)，則垂直變動速度等於 VSR x cos。圖 1 藍 色方框區塊 (33 km×54 km) 為大範圍 PS-InSAR 計算的區域，包含彰化和雲林地區，並以「水利雲 88」水準點做為 PS-InSAR 計算中穩固不動的基準 點。紅色方框區塊(12 km×12 km)位於虎尾地區，
用來做為分析豐枯水期對地層下陷影響的計算區 域，基準點為「大東堤防 0+00」水準點。採用基
準點是因為 PS-InSAR 的計算結果為相對垂直變動 速度，因此必須依靠基準點將其轉換成絕對垂直變 動速度，實際做法是在 PS-InSAR 的計算過程中，
將最接近這兩個水準點的 PS 點設為變動速度 0 的 基準點。經濟部水利署(2011)的報告顯示豐水期與 枯水期會影響地層下陷的速度，因此若基準點位於 地層下陷區域，會造成無法正確分析豐枯水期對地 層下陷的影響，因此所挑選的基準點應該為較穩固 不動的固定點。圖 2 和圖 3 顯示這兩個水準點「水 利雲 88」與「大東堤防 0+00」在 2007 到 2010 年 期間的高程變動量分別為 1mm 和 2.4cm，高程變 動相對較小，為較穩定的固定點，用來做為本文 PS-InSAR 計算的基準點。
圖 1 計算區域
距離方向 方
位 方 向
圖 2 「水利雲 88」水準點在 2007~2010 年期間的高程變動情形
圖 3 「大東堤防 0+00」水準點在 2007~2010 年期間的高程變動情形
4. 使用 PS-InSAR 測量大 範圍地層下陷與成果評
估
此計算例採用的 13 組 ALOS 影像對資料如表 2 所示，以編號 13 的 ALOS 影像為主影像，並與 其它副影像組成影像對資料。垂直基線為兩張影像 所對應的兩天線基線長的垂直分量，隨著垂直基線 越長，影像對品質也會隨著空間降相關的影響而降 低，造成最後 PS-InSAR 計算成果的品質下降。時 間基線為影像對中兩影像資料取得的時間間隔，負 號和正號分別代表副影像的取樣時間在主影像之 前和之後。
根據實際的計算經驗顯示，當採用 ALOS 影 像進行小區域的 PS-InSAR 計算時，使用垂直基線 長小於 1750 m 之影像對組合即可得到合理的計算 結果；然而針對大範圍地區的計算，由於計算的不 確定因子增加且空間降相關的影響，則需要使用品
質更高的影像對資料。簡單來說，使用更短的垂直 基線門檻值能有更大的機會得到更好的計算結果，
本計算例即採用垂直基線長小於 1300 m 之影像對 來進行計算。
表 2 ALOS 影像對資料 影像
編號
取樣日期 (年/月/日)
垂直基線 (m)
時間基線 (days) 2 2007/02/15 330.3 -966 3 2007/08/18 657.9 -782 4 2007/10/03 333.8 -736 5 2008//01/03 580.6 -644 7 2008/07/05 721.8 -460 10 2009/01/05 1223.9 -276 12 2009/08/23 240.6 -46
13 2009/10/08 / /
14 2010/01/08 90.1 92 15 2010/02/23 641.1 138 16 2010/07/11 964.5 276 17 2010/08/26 1262.4 322 18 2010/10/11 1254.6 368 19 2010/11/26 1206.1 414 2007年08月 2008年05月 2009年06月 2010年5月
水利雲88 79.208 79.209 79.209 79.209 76
77 78 79 80 81 82 83
高程 (m )
2007年08 月
2008年05 月
2009年06
月 2010年5月 大東堤防0+00 31.549 31.538 31.525 31.529
28 29 30 31 32 33 34 35
高程 (m )
將 13 組影像對進行干涉計算得到 13 張同調性 影像和 13 張主要包含地形效應和變動效應的干涉 圖。同調性影像如圖 4 所示，同調性值的值域為 0 到 1，單一像元的同調性越高代表該地面目標區塊 在兩張影像上的相關性越高，能夠提供更可靠的相 位資料進行 PS-InSAR 計算。本測試例先將在所有 同調性影像中其同調性值大於 0.26 的點挑選為 PS 候選點 PSC，配合特定過濾 PSC 的演算法，進一 步挑選出 10386 個 PSC 點。上述同調性門檻值 0.26 的決定來自該地區 PS-InSAR 的計算經驗，能在 PSC 點數和電腦計算負擔能力間取得平衡。
將 干 涉 圖 中 的 相 位 資 料 扣 除 掉 SRTM-3
Version 4 Digital Elevation Model (DEM)資料計算 得到的地形效應相位分量，得到主要包含變動效應 的差分干涉圖，如圖 5 所示。據經濟部水利署(2011) 的報告指出，目前彰化地區地層下陷中心有二林鎮、
溪湖鎮與溪州鄉；而雲林地區的主要下陷區域為虎 尾鎮、土庫鎮、元長鄉與褒忠鄉，並往南投方向逐 漸趨緩。從部分時間基線較長之差分干涉圖可明顯 看出其干涉條紋顯示的變動趨勢與上述的情形相 同，可知雷達干涉測量具有偵測出地表下陷變動趨 勢的能力。最後將所有 PSC 點在所有差分干涉圖 中的相位值連結後，進行 PS-InSAR 計算，總共萃 取出 10360 個 PS 點並轉換得到其垂直變動速度。
2007/02/15 2007/08/18 2007/10/03 2008/01/03 2008/07/05
2009/01/05 2009/08/23 2010/01/08 2010/02/23 2010/07/11
2010/08/26 2010/10/11 2010/11/26
圖 4 同調性影像(日期為對應副影像的取樣時間)
2007/02/15 2007/08/18 2007/10/03 2008/01/03 2008/07/05
2009/01/05 2009/08/23 2010/01/08 2010/02/23 2010/07/11
2010/08/26 2010/10/11 2010/11/26
圖 5 差分干涉圖 (日期為對應副影像的取樣時間) 圖 6 顯示大範圍 PS-InSAR 計算成果的垂直變
動速度圖，圖 7 使用不同的顏色來表示各水準點的 下陷速度等級。每一個水準點每一年會有一筆高程 資料，但各水準點測量的月份未必都相同，其中，
有的水準點比較新，有的後來遺失，有的資料不完 整，所以我們只挑出 95 年 10 月至 99 年 5 月每年 皆有資料的 200 個水準點來進行下陷速度之計算 與分析比對，最後再直接使用 95 年 10 月與 99 年 5 月的高程值來計算每個水準點在這段期間的平 均年下陷量。比較圖 6 與圖 7 並檢視相關的統計數 據可發現，在彰化和雲林北部地區，PS-InSAR 數 據顯示約每年下陷 2~4cm 而精密水準點顯示約每 年下陷 2~6cm，證實 PS-InSAR 測量地表下陷的良
好應用潛力。進一步分析兩者的成果，可看出圖 6 顯示的地層下陷情形較為趨緩，推論原因為 ALOS 影像所採用的 SAR 訊號波長為 23.6 cm，一個干涉 條紋所涵蓋的 SR 方向上的相對變動量為 11.8cm(=
波長/2)，波長越長，干涉條紋越不易顯示微小的 變動量，對地表變動的偵測就越不靈敏，因此 ALOS 影像所採用的波長為造成地層下陷情形被 低估的原因之一。未來建議採用 TerraSAR-X 資料 進行地層下陷的監測作業，其訊號脈衝波採用波長 更短的 3.1cm 的 X 波段，將可偵測到更細微的地 表變動量。
圖 6 垂直變動速度之位置分布圖
圖 7 各水準點下陷速度的位置分布圖
將圖 6 之成果利用克利金法（Deutsch and Journel, 1997）內插得到圖 8 的成果。比較圖 8 與 圖 7 可明顯看出，PS-InSAR 與水準測量得到的地 表年下陷量分別約 2cm~4cm 與 2cm~6cm，全區的 下陷趨勢相當一致，兩者都在濁水溪南北側各有下 陷較嚴重的局部區域。另外附上圖 9 為水準點下陷 量內插圖，比較圖 8 與圖 9 可更清楚發現兩者求定 的全區整體的下陷趨勢大致相同，但是 PS-InSAR 得 到 的 平 均 年 下 陷 量 (2cm~4cm) 較 水 準 測 量 者 (2cm~6cm)略小(1.4cm)。在 95 到 99 年每年都有高 程資料的 200 個水準點上，以每一個水準點為中心，
分別往東、西、南、北四個方向延伸 200 公尺，成 為一個 400m x 400m(約 13 x 13 個像元)的正方形區 塊，計算這區塊內所有 PS 點的平均年下陷量，並 與水準點計算出來的年平均下陷量做一比較。兩者 求定的高程變動速度差值之絕對值，平均值為 1.4cm/year ， 最 小 值 為 0.0cm/year ， 最 大 值 為 4.5cm/year ， 均 方 根 值 RMSD(root mean square
difference)為 1.4cm/year。表 2 顯示 10360 個 PS 點 上的年平均下陷量是使用 96 年 2 月 15 日到 99 年 11 月 26 日期間共 13 個不同的取樣日期取得的 SAR 影像經 PS-InSAR 計算而得，而 200 個水準點 上的平均年下陷量則是使用 95 年 10 月與 99 年 5 月的兩個日期的高程值計算得到，考量兩者的測量 時間點、取樣個數與測量精度之不同，在檢視數據 的顯著性下，兩者測量得到的年平均下陷量大抵可 視為無顯著差異。
總結來說，採用 SAR 影像資料配合 PS-InSAR 可用來監測大範圍的地層下陷，未來可採用品質更 好且尚在服役中的衛載雷達資料如 TerraSAR-X 來 進行實際的監測作業，將可計算得到更精細地面目 標且精度更高的地表下陷速度，並以地面測量結果 如精密水準測量資料和 GPS 觀測資料做為輔助資 料，建立基準點並協助驗證 PS-InSAR 的計算成 果。
圖 8 內插得到的垂直變動速度場分布圖
圖 9 95~99 年水準點下陷量內插圖
5. 豐枯水期的地層下陷之 分析比較
本測試例分別 挑選豐枯 水 期間取樣得到 的 ALOS 影像，並以編號 13 的 ALOS 影像為主影像，
其取樣時間在豐枯水期交界時間的 10 月，避免混 合豐水期或枯水期的資料而影響計算結果。將主影 像與其它多張副影像分別組成兩組影像對資料，簡 稱為豐水組和枯水組，如表 3 所示。將豐水組和枯 水組分別以同調性門檻值 0.28 與 0.27 挑選出 10046 個和 10426 個 PSC 點，最後萃取出 9878 個 和 10172 個 PS 點並得到 PS 點上的垂直變動速度。
上述兩個同調性門檻值 0.28 與 0.27 為根據計算經 驗所決定的值，能在 PSC 點數和電腦計算負擔能 力間取得平衡。
圖 10 和圖 11 分別為豐水組和枯水組的垂直變 動速度圖，兩者垂直變動速度採用相同的顏色尺度 俾利分析，從整體的趨勢分析來看，枯水組的地表
下陷情形確實比豐水組來得嚴重。根據統計資料，
豐水組的 PS 點中，最大上升速度為+5.9cm/year，
最大下陷速度為-4.5cm/year，平均垂直變動速度 Vavg為+0.2cm/year，垂直變動速度和 Vavg的差值之 均方根值 RMSD 為 1.2cm/year；在枯水組的 PS 點 中，最大上升速度為+4.7cm/year，最大下陷速度為 -7.2cm/year，平均垂直變動速度 Vavg為-1.1cm/year，
垂直變動速度和 Vavg的差值之均方根值 RMSD 為 1.8cm/year。
圖 12 顯示豐水組和枯水組中垂直變動速度的 統計直方圖，由圖中可明顯看出枯水組相較豐水組 來說，其 PS 點的整體下陷速度較大。前述的測試 成果與分析顯示，PS 點密度比水準點密度高，高 精度的水準資料可做為 PS-InSAR 計算的基準點和 檢核點，PS-InSAR 加上數量足夠且時間分布合宜 的衛載雷達資料可用來分別計算豐枯水期間的大 區域地層下陷量，具有良好的應用潛力。
表 3 豐水期和枯水期的影像對資料
豐水期 枯水期
主影像 編號
副影像 編號
垂直基線 (m)
時間基線 (days)
主影像 編號
副影像 編號
垂直基線 (m)
時間基線 (days) 13 3 659.099 -782 13 1 1961.321 -1012 13 7 721.307 -460 13 2 326.149 -966 13 8 2567.308 -414 13 5 581.225 -644 13 11 1388.46 -92 13 10 1224.048 -276 13 12 237.23 -46 13 14 90.067 92 13 16 969.216 276 13 15 641.108 138 13 17 1264.509 322 13 19 1204.837 414
/ / 13 20 2142.353 506
圖 10 豐水組垂直變動速度位置分布圖
圖 11 枯水組垂直變動速度位置分布圖
圖 12 豐水組和枯水組中垂直變動速度的統計直方圖
6. 結論與建議
第 一 個 測 試 案 例 使 用 ALOS 影 像 並 配 合 PS-InSAR 來計算大範圍彰化和雲林地區的地層下 陷，成功得到 33 km×54 km 範圍內共 10386 個 PS 點的垂直變動速度，其 PS 點密度約為 5.8PS/km2， 相較地面測量方法如精密水準測量或 GPS 測量，
可以更低的成本得到分布在大範圍且密度更高的 監測點的垂直變動速度，提供更多且均勻分布的監 測點進行可靠的大範圍區域地層下陷監測之實際 作業。
一般考量到成本與時間，地面測量方法在實務 上較難個別針對豐枯水期來做地層下陷的監測工 作。第二個測試案例分別挑選豐枯水期的 ALOS 影像對資料，並利用 PS-InSAR 成功計算得到豐水 組和枯水組的地表下陷，其結果顯示枯水組的整體 下陷速度較豐水組來得嚴重，可見 PS-InSAR 使用 衛載雷達資料分別計算豐枯水期間的地層下陷之 應用潛力。
目前 ALOS 已經除役，建議未來可採用尚在 服役中且更高解析度、性能更好的衛載雷達資料 (例如 TerraSAR-X)來測量地表下陷，俾得到更精細 地面目標的地表變動速度。
致謝
本研究承蒙兩位匿名審查委員提供珍貴的修 正意見，國立中央大學 吳究教授與國立台灣大學 胡植慶教授以及達雲科技有限公司對本研究的大 力支持，特此一併致以由衷謝忱。
參考文獻
經濟部水利署，2011。100 年度多重感應器應用於 臺北、彰化及雲林地區地層下陷監測與機 制探討。
Costantini, M. and Rosen, P. A., 1999. A Generalized Phase Unwrapping Approach for Sparse Data.
Proceedings of IEEE, 1999 International Geoscience and Remote Sensing Symposium, Hamburg, Germany, Vol. 1, pp. 267-269.
Delaunay, B., 1934. Sur La Sphere Vide. Bulletin of Academy of Sciences of Union of Soviet Socialist Republics, pp. 793-800.
Deutsch, C. V. and Journel, A. G., 1997. Gslib:
Geostatistical Software Library and User’s Guide, Oxford University Press, New York.
Ferretti, A., Prati, C. and Rocca, F., 2000. Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry. IEEE Transactions on Geoscience and Remote Sensing, Vol. 38, No. 5, pp. 2202-2212.
Ferretti, A., Prati, C. and Rocca, F., 2001. Permanent Scatterers in SAR Interferometry. IEEE Transactions on Geoscience and Remote Sensing, Vol. 39, No. 1, pp. 8-20.
Hanssen, R., 1998. Atmospheric Heterogeneities in ERS Tandem SAR Interferometry. DEOS Report No.98.1, Delft University press, Delft, Netherlands.
Lee, J. S. and Pottier, E., 2009. Polarimetric Radar Imaging: From Basics to Applications. CRC Press, Boca Raton, USA.
Lillesand, T. M., Kiefer, R. W. and Chipman, J. W., 2004. Remote Sensing and Image Interpretation.WILEY Press, New York, USA.
Madsen, S. N. and Zebker, H. A., 1998. Synthetic Aperture Radar Interferometry: Principles and Applications. John Wiley and Sons, New York, USA.
Mora, O., Mallorqui, J. J. and Broquetas, A., 2003.
Linear and Nonlinear Terrain Deformation Maps from a Reduced Set of Interferometric SAR Images. IEEE Transactions on Geoscience and Remote Sensing, Vol. 41, No. 10, pp.
2243-2253.
Rosen, P. A., Hensley, S., Joughin, I. R., Li, F. K., Madsen, S. N., Rodriguez, E. and Goldstein, R.
M., 2000. Synthetic Aperture Radar Interferometry - Invited Paper. Proceedings of IEEE, Vol. 88, Issue 3, pp.333-382.
Santitamnont, P., 1998. Interferometric SAR Processing for Topographic Mapping. Doctor’s Thesis, Fachbereich Bauingenieur- und Vermessungswesen, Leibniz Universität Hannover , Deutschland.
Ulaby, F. T., Haddock, T. F. and Austion, R, T., 1988.
Fluctuation Statistics of Millimeter-Wave Scattering from Distributed Targets. IEEE Transaction on Geoscience and Remote Sensing, Vol. 26, No. 3, pp. 268-281.
Xu W. and Cumming I., 1999. A Region-Growing Algorithm for InSAR Phase Unwrapping. IEEE Transactions on Geoscience and Remote Sensing, Vol. 37, Issue 1, pp. 124-134.
Zebker, H. A., Rosen, P. A., Goldstien, R. M., Gabriel, A. and Werner, C. L., 1994. On the Derivation of Coseismic Displacement Fields Using Differential Radar Interferometry: The Landers Earthquake. Journal of Geophysical Research, Vol. 99, No. B10, pp. 19617-19634.
Zhang, L., Ding, X. and Lu, Z., 2011. Ground Settlement Monitoring Based on Temporarily Coherent Points between Two SAR Acquisitions.
ISPRS Journal of Photogrammetry and Remote Sensing, Vol. 66, Issue 1, pp. 146-152.
1 Master, Department of Geomatics, National Cheng Kung University Received Date: Mar. 01, 2014
2 Associate Professor, Department of Geomatics, National Cheng Kung University Revised Date: May. 07, 2014
3 Master, Department of Geomatics, National Cheng Kung University Accepted Date: May. 15, 2014
*.Corresponding Author, Phone: 886- 6-2370876 ext.838, E-mail: tsayjr@mail.ncku.edu.tw
Improved PS-InSAR Algorithm for Determining Subsidence in Central Taiwan
Chia-Hsiang Yang 1 Jaan-Rong Tsay 2* Po-Tsung Su 3
ABSTRACT
This paper presents an improved PS-InSAR algorithm and preliminary test results on subsidence determination in central Taiwan. PS-InSAR with ALOS PALSAR images is used to evaluate the land surface subsidence in central Taiwan mostly covered by vegetation and mountains. The PS-InSAR results are analyzed and then evaluated by comparing with precise leveling data. Both results determined by PS-InSAR and precise leveling show similar subsidence pattern. The statistic figures derived from the absolute differences in vertical displacement velocity vectors on 200 PSs and 200 benchmarks show that the absolute differences have the mean 1.4 /year, the minimum 0.0 /year, the maximum 4.5 /year and the root mean square difference (RMSD) 1.4 /year. The conclusion infers that monitoring large-scale subsidence in central Taiwan by means of PS-InSAR is feasible and applicable, and it has a good application potential for a national-level regular subsidence monitoring project particularly with advanced radar data.
