利用永久散射體雷達干涉技術分析地表變形：以臺灣西南部平原為例

(1)

國立臺灣大學理學院地質科學系碩士論文

Department of Geosciences College of Science

National Taiwan University Master Thesis

利用永久散射體雷達干涉技術分析地表變形 : 以臺灣西南部平原為例

Analysis of surface deformation based on PS-InSAR technique: Case studies in coastal plain, SW Taiwan

童忻 Hsin Tung

指導教授：胡植慶博士 Advisor: Jyr-Ching Hu, Ph.D.

中華民國 97 年 7 月 July, 2008

(2)

國立臺灣大學碩士學位論文

口試委員會審定書

利用永久散射體雷達干涉技術分析地表變形 : 以臺灣西南部平原為例

Analysis of surface deformation based on PS-InSAR technique: Case studies in coastal plain, SW Taiwan

本論文係童忻君（R95224206）在國立臺灣大學地質科學所完成之碩士學位論文，於民國九十七年七月九日承下列考試委員審查通過及口試及格，特此證明

口試委員：

（簽名）

（指導教授）

(3)

Abstract

Differential Interferometric Synthetic Aperture Radar (DInSAR) technique has increasingly been applied in geosciences as a powerful tool to monitor land surface deformation. In addition, the large archives of SAR data enable us to monitor ground deformation continually. However, temporal and spatial decorrelations of radar signal have prevented this technique from more frequent utilization. An advanced InSAR processing technique that tracks the signals of discreted point-wise targets called Persistent Scatterers Interferometry Technique (PSInSAR) was herein applied to our research area in order to better understand the deformation patterns. Thus the PSInSAR method allows us to know the deformation of these discrete points that were minimally affected by decorrelation of radar signals through time.

The active orogeny of Taiwan generated large amount of surface deformation that were partitioned throughout the island, which provides a great opportunity for analyzing tectonic displacements through interferometric methods. Here we present a PS-InSAR result deduced from 1996-1999 time spans in and around Tainan Tableland of southwestern Taiwan close to the deformation front. About 1,650 persistent scatterers are deduced from the processing of 33 SAR images (Tract: 232, frame, 3141) in an area of 144 km². The maximum slant range displacement (SRD) rate along the radar line of sight (LOS) toward to satellite is about 8 mm which is close to the borehole data after transferring SRD rate to vertical deformation rate. Besides, the nonlinear deformation evolutions show a respectable correlation with seasonal rainfall on the Tainan Tableland near the Houchiali fault. One of the reasons might be the behaviors of fault creeping during the rainy seasons. These results indicate that the Tainan Tableland is an aseismic area.

(4)

In addition, anthropogenic ground subsidence due to massive pumping of groundwater is one of problems in Taiwan. The Yunlin County located in the southwestern coastal region of Taiwan is one of the most counties with serious land subsidence because of the agricultural needs. Here, we also represented a both stacking DInSAR and PSInSAR results deduced from 1996-1999 time spans for monitoring of land subsidence in this area. Both DInSAR and PSInSAR results show the Baojhou, Tuku Township and northern Mailiao reveal a maximum subsidence rate of about 6 cm/yr along LOS and the Taishi Township located on the coastal area reveals a subsidence rate of 1.58 cm/yr, which is quite coincident with precise leveling result.

These two results have proven that the effective reduction of labor and cost could be achieved by using this technique on monitoring land subsidence in Yunlin County.

Key words: InSAR, DInSAR, PSInSAR, Tainan Tableland, Houchiali Fault, Tawan Lowland, Yunlin County, Surface Deformation, Land Subsidence.

(5)

摘要摘要摘要摘要

近年來，差分合成孔徑雷達干涉技術已經成為地球科學研究上一項利於觀測地表變形的工具。透過所獲取的大量影像，有助於我們更能夠連續性的觀測地表變形。但是，雷達訊號隨著空間基線的增大和時間的增長，降低雷達訊號的相關性，限制了此項技術無法更廣泛的應用。永久散射體雷達干涉技術，利用萃取分散的穩定點目標，能有效克服時間以及空間上低相關的限制；因此，本研究將此項技術應用在所選取的研究區域內，以期能更清楚的了解地表變化的時空演變。

台灣位於歐亞大陸板塊與菲律賓海板塊之間的碰撞縫合帶，構造活動相當劇烈，提供了差分合成孔徑雷達干涉技術一個良好的材料來分析地表變形。本研究利用 1996 年至 2000 年間所接收的 33 張雷達影像，以永久散射體雷達干涉技術分析位於變形前緣的台南台地之地表變形；在 144 平方公里區域內，挑選約有 1,650 個永久散射點做為分析的點目標。分析結果顯示，台南台地上最大視衛星方向上的抬升速率為每年 8 公厘，將其座標轉換至垂直方向上的抬升速率，與層序沉積定年的結果相當接近，暗示此區長期抬升速率與短期抬升速率相近，屬於無震區。

除此之外，靠近後甲里斷層的永久散射體的變形行為皆呈現出與季節性降雨有高度的相關性，其原因可能是斷層在雨季發生潛移活動的行為，呼應了此區為無震區的推論，但此推論還需要更多的測量結果才能確定。

除了構造活動外，人為活動也是造成地表位移的原因之ㄧ。台灣西南沿海的雲林縣是以農業為主的區域，因灌溉而長期抽取地下水造成此地區地層下陷的情況相當嚴重；因此，本研究再利用累加差分合成孔徑雷達干涉技術與永久散射體干涉技術，監測雲林地區 1996 年至 1999 年 6 月間的地層下陷時空演變。兩項結果都顯示出，位於內陸的褒忠鄉、土庫鄉以及麥寮鄉北部區域擁有最大的下陷量，

約為每年 6 公分，沿海地區的台西鄉則呈現出每年 1.58 公分的下陷量，與精密水準測量結果相當吻合。此結果證明了利用此項技術，能夠有效的減少人力以及財力的花費，得到雲林地區地層下陷點狀以及面狀的資訊。

(6)

關鍵字關鍵字

關鍵字關鍵字：：：：差分合成孔徑雷達干涉技術、永久散射體雷達干涉技術、台南台地、後

甲里斷層、大灣低地、雲林縣、地表變形、地層下陷。

(7)

致謝致謝致謝致謝

兩年的研究所生涯，終於告一段落。這一篇論文能夠順利完成，首先要感謝我的指導老師－胡植慶博士；感謝老師在研究以及課業方面的協助，給我極高的自由度，讓我能夠選擇有興趣的研究目標與方向。在遙測的世界裡，我很開心。

再來要感謝我的口試委員，陳文山博士、余騰鐸博士、謝嘉聲博士以及顏君毅博士，於學生口試期間撥冗擔任我的口試委員，對於本論文以及研究上不足之處給予非常多的指導和幫助，使本論文接近完整。陳文山老師在地質解釋方面多次的教導，提供了我遙測以外的資料與觀念。余騰鐸老師是讓我對遙測產生興趣的最初原因，慶幸自己能在大學時遇見了老師，真的非常感謝。謝嘉聲老師總是提供我另外一層思考，也給與我很大的肯定，讓我有信心走下去。顏君毅老師時常一大早就被我用 MSN 攻擊，還要一直忍受我的糾纏與笨，如果沒有你，我應該早就放棄了吧。除此之外，盧佳遇博士與中央大學張中白博士教導了我許多報告上以及研究上的要領，從老師身上，我得到不只是研究上面的精進，還有更多處事的態度。真的感謝各位老師，希望我還有機會能向各位老師學習。

除此之外，要特別感謝香港理工大學丁曉利博士、大陸西南交通大學劉國祥博士以及陳強博士；有了你們的協助，讓我度過了研究途中最大的瓶頸，成功的在今年達成目標，順利畢業。另外提供論文重要資料的吳逸民老師、賈儀平老師，

饒瑞鈞博士以及成功大學的志憲，感謝有你們，我的論文才能完整。

再來，我要感謝黃孟涵學長，幫我一點一點的建立遙測的觀念，就像是我的指導老師ㄧ般；一起工作的週末，台灣中部、南部都有我們工作的紀錄，真的謝謝你對我的百般協助。而四人小組中善解人意的雅琳、聰明認真的黃蓉、氣質出眾的婷兒，加上採購部主任宜芬，真的是我這段期間最美好也最捨不得的回憶。

感謝妳們在我傷心、開心的時候在我身邊，聆聽我的煩惱，解決我的問題，在這段期間備受妳們的照顧。有妳們在身邊，是我每天開心的原動力，超愛妳們的，

希望有機會再一起工作唷。

(8)

另外，實驗室一起奮鬥的顏色、詠恬，以及昭榮學長、育雅老師、宜蓉、佳漢、黃鐘、蚯蚓、冠全、彭彭、Maple、偉聖，謝謝你們這段日子的包容。中央的貓王、小黑、嘉貞、育芳、萬慧、勤怡、小熊、偉恆，認識你們，是我很大的收穫。乙君、踢哀恩、液晶、Ray、鴉哥、阿和、珀儂、信樺、阿天、林桑，因為有你們，我的研究生活才如此的充實與快樂。還有進上、張立、中正的學長姐們以及其他在這裡認識的好朋友們，兩年的時間，承蒙了大家的照顧，也祝福大家都能夠一切順心。

最後我想感謝的是我最重要的家人，爸爸媽媽和姊姊。謝謝爸爸媽媽每天辛苦的工作，讓我的學生生活無後顧之憂，還有忍受我回台北以來的壞脾氣；謝謝姊姊每天載我上學，讓我可以早早到學校工作。還要謝謝楹丰一直以來的陪伴。

謝謝你們在生活上的照顧和支持，謹以此論文獻給我最愛的你們。

另外還有很多曾經幫助過我的朋友，因為有大家的幫助，我才能有今天的成果。謝謝大家。

(9)

List of Figures

Figure 1.1. Tectonic framework of Taiwan in 3-D ... 3

Figure 1.2. Average subsidence rate in Taiwan Island from 2002 to 2006 ... 5

Figure 1.3. SAR images of track 232 and frames 3141 and 3123... 8

Figure 2.1. Imaging geometry of a RAR acquisition ... 12

Figure 2.2. SAR pulse transmitting and receiving mode... 13

Figure 2.3. Amplitude image of ERS SAR in Northern Taiwan ... 14

Figure 2.4. Radar interferogram of a region near San Francisco, California ... 15

Figure 2.5. Imaging geometry for repeat- and single-pass radar interferometry... 16

Figure 2.6. Imaging geometry for satellite radar interferometry... 17

Figure 2.7. Interferometric configuration ... 19

Figure 2.8. Imaging geometry for differential interferometry... 20

Figure 2.9. Three-pass differential interferometry configuration... 22

Figure 2.10. Four-pass differential interferometry configuration... 23

Figure 2.11. Phase simulated for a scatterer pixel ... 25

Figure 2.12. Numerical simulation results for D and phase stability ... 28 _A Figure 2.13. Example of Delauney triangulation ... 31

Figure 2.14. Layout of linear deformation and DEM error extraction ... 33

Figure 2.15. Layout of nonlinear deformation and atmospheric artifacts extraction .. 35

Figure 3.1. Geodynamic framework and structural map of Taiwan ... 39

Figure 3.2. Active deformation in Tainan Tableland ... 40 Figure 3.3. Geological structural setting and location of GPS and leveling stations,

elevation measurement and uplift rate measurement of precise leveling41

(13)

Figure 3.4. Three schematic structural models of the Tainan Tableland ... 43

Figure 3.5. Typical subsurface geological column and stratigraphy in Tainan area.. 45

Figure 3.6. Totally 32 interferograms with phase wrapped ... 48

Figure 3.7. Selected points and the Delaunay Triangulations ... 49

Figure 3.8. The SRD rate over PS points ... 50

Figure 3.9. Linear SRD rate and Topographic profile ... 51

Figure 3.10. Time series of deformation ... 52

Figure 3.11. Deformation evolutions of points... 53

Figure 3.12. Locations of profiles, cores, trough ditch, and outcrops ... 54

Figure 3.13. Long-term deformation rate and topographic profile... 55

Figure 3.14. The geometrical relationship of the SRD and vertical deformation ... 56

Figure 3.15. GPS components at station CK01 during 1996 to 2006. ... 57

Figure 3.16. Location map of precipitation stations and groundwater well ... 58

Figure 3.17. Monthly rainfall data at Tainan and Yongkang ... 59

Figure 3.18. Map of seismicity in Tainan city and surrounding areas... 60

Figure 3.19. Numbers of earthquakes happened in study area... 61

Figure 3.20. Nonlinear deformation of five PS points ... 62

Figure 3.21. Comparison of nonlinear deformation evolution and seasonal rainfall .. 63

Figure 4.1. Average land subsidence area during 2002 to 2006... 68

Figure 4.2. Cumulative subsidence in Yunlin County from 1992 to 2007 ... 69

Figure 4.3. Geographical location and environment of Yunlin coastal area... 70

Figure 4.4. Conceptual hydrogeologic profile of the Chuoshui River alluvial fan ... 71

Figure 4.5. Stratigraphic profile and subsidence measuring sensorsin... 72

(14)

Figure 4.6. Average land subsidence of Yunlin area during 1992 to 1994 ... 73

Figure 4.9. 17 interferograms produced with phase wrapped. ... 79

Figure 4.10. Unwrapped SRD rate and land subsidence contour... 80

Figure 4.11. 32 interferograms produced with phase wrapped. ... 81

Figure 4.12. Selected points and Delaunay Triangulations ... 82

Figure 4.13. The rate of SRD from PS points of the study area ... 83

Figure 4.14. Location of precise leveling station ... 85

Figure 4.15. Deformation evolutions of three selected PS points ... 86

(15)

Lists of Tables

Table 3.1. Descending orbit data processed for the Tainan area ... 47

Table 4.1. Descending orbit data processed for the Yunlin area... 77

Table 4.2. Differential Interferograms for Yunlin area ... 78

Table 4.3. Average SRD rate of the downtown areas in study area. ... 84

Table 4.4. Precise leveling measurements for four stations in Yunlin area ... 86

(16)

Lists of Abbreviations

ALOS advanced land observing satellite CNES Centre National D'Etudes Spatiales CSA Canadian Space Agency

CWB Central Weather Bureau DEM digital elevation model

DEOS Delft Institute of Earth oriented space research

DIAPASON differential interferometric automated process applied to survey of nature

DInSAR differential InSAR

EMD empirical mode decomposition

ENVISAT European Space Agency environment satellite

ERS European Space Agency Earth remote sensing satellite ESA European Space Agency

EU Eurasian plate

GPS global positioning system

InSAR interferometric synthetic aperture radar ITRI Industrial Technology Research Institute JERS-1 Japanese Earth resources satellite-1 JPL Jet Propulsion Laboratory

LOS line of sight

NSDAJ National Space Development Agency of Japan NASA National Aeronautics and Space Administration PALSAR phased array type L-band synthetic aperture radar PS persistent scatterers

PRF pulse repetition frequency

(17)

PSP Philippine Sea plate Radar radio detection and range RAR real aperture radar

SAR synthetic aperture radar SBAS small baseline subset SLAR side-looking RAR

SNAPHU statistical-cost, network-flow algorithm for phase unwrapping SNR signal to noise ratio

SRD slant range displacement

SRTM shuttle radar topography mission

StaMPS Stanford method for persistent scatterers THSR Taiwan High Speed Rail

WRA Water Resources Agency

(18)

(19)

Chapter 1 Introduction

(20)

1.1 Motivation

Recently, InSAR (Interferometric Synthetic Aperture Radar) has been used widely in many fields such as glacial processes, landslide, terrain classifications, ocean current and vegetation parameters [Massonnet and Feigl, 1998; Rosen et al., 2000] due to its “all weather” capability and the day/night access of the sensors. This remote-sensing tool distinguishes itself from other geodetic measurements in that requires no presence in the filed and is available worldwide. Basically, InSAR could map surface deformation at dense pixel of 20-100 m over almost unlimited spatial distances and is particularly sensitive to vertical displacements comparable to the precision of leveling measurements. As a result, InSAR is also applied to monitor crust deformation, volcanism and land subsidence [Bürgmann et al., 2000].

Although spaceborne Differential InSAR (DInSAR) has already proven its potential for mapping ground deformation, the atmospheric perturbation and phase decorrelation have prevented conventional DInSAR from achieving full operational capability. On the other hand, the Persistent Scatterers (PS) technique is a powerful and fully operational tool for monitoring ground deformation on a high spatial density of point-wise targets by an extensive time series of SAR data [e.g. Ferretti et al., 2000, 2001; Colesanti et al., 2003]. Moreover, the analysis of pixel-by-pixel character of the phase stable targets makes it feasible to detect the ground deformation in dense vegetation or mountainous area as long as the men-made structures or exposed rocks are available.

Taiwan Island is located along a segment of the oblique convergent boundary between the Philippine Sea Plate (PSP) and Eurasian Plate (EP), where collision and subduction processes occur (Figure 1.1). The convergent rate between these two plates is about 82 mm/yr [Yu et al., 1999; Zang et al., 2002], due to this oblique convergence,

(21)

the young Taiwan orogenic belt is an ideal natural laboratory to study the tectonic processes based on geodetic measurements. Consequently, the island-wide GPS (global positioning system) network was established from 1991 to meet acquirement of the assessment of seismic hazards [e.g., Yu et al., 1997]. In addition, the conventional DInSAR has been applied to detect the coseismic deformation of the 1999 Chi-Chi earthquake [Pathier et al., 2003; Chang et al., 2004] and the active deformation in Tainan area [Fruneau et al., 2001; Huang et al., 2006]. However, the limitations of conventional DInSAR offer lots of uncertain of deformation measurements. Thus, one of the goals of this study is to apply PSInSAR technique and discuss the feasibility in Tainan area.

Figure 1.1. Tectonic framework of Taiwan in 3-D [Angelier et al., 1999].

(22)

Ground subsidence induced by heavy withdrawal of groundwater has resulted in environmental hazard and potential risk in Taiwan (Figure 1.2), particularly in the Chou-Shui River alluvial fan where the Yunlin section of the Taiwan High Speed Rail (THSR) had been constructed. Seasonal effects of land subsidence occurring in the study area had been estimated using a regression analysis of a series of weekly GPS height solutions. The average rate of ground subsidence in this area over the period of 1995-2001 was 3 cm/yr. Based on data collected at the piezometer, the variation of land subsidence rate appears to be associated with an groundwater level, which drop gradually during winter and either remains constant or rises during summer time.

Consequently, land subsidence rates vary considerably from 1.5 cm/yr for the summer time to 9.0 cm/yr for the winter time [Chang and Wang, 2006].

Therefore, the second aim of this study is to monitor land subsidence using SAR interferometry. However, due to the dense vegetation in Yunlin area, it is hard to get the time series distribution of land deformation by using conventional DInSAR technique. I used the stacking DInSAR and PSInSAR technique in order to obtain the temporal and spatial variations of land subsidence. For the PSInSAR technique, I adopt the idea from Ferretti et al. [2000, 2001] and Mora et al. [2003] to develop the programs for this study.

(23)

Figure 1.2. Average subsidence rate in Taiwan Island from 2002 to 2006 [Land Subsidence Database].

(24)

1.2 Spaceborne Radar History

Seasat is the first SAR mission in space lunched in July 1978 primarily for ocean studies by NASA/JPL (National Aeronautics and Space Administration/Jet Propulsion Laboratory). It was a L-band (23 cm wavelength) synthetic aperture radar system with horizontal (HH) polarization. Unfortunately, it stopped operating after 99 days of lunching. However, the imagery was found to have a great deal of potential for scientific purposes.

Base on the success of Seasat, two shuttle missions of one week duration called SIR-A and SIR-B were operated in November 1981 and October 1984 respectively.

Both SIR-A and SIR-B were L-band synthetic aperture radars with HH polarization.

Unlike Seasat and SIR-A were recorded analogously on tapes, SIR-B was already equipped with full digital recording [Hanssen, 2001].

The European Space Agency (ESA) started to launch its first satellite equipped with a SAR in July 1991, the ERS-1. It was a C-band (5.66 cm wavelength) radar operating in VV polarization. In April 1995, an identical mission satellite named ERS-2 (European Space Agency Earth Remote Sensing Satellite 2) was launched and put in the same orbit as ERS-1. From July 1995 to May 1996, six so-called 'tandem pair' were collected with an exact 24 hours interval between these two satellite.

The JERS-1 (Japanese Earth Resources Satellite-1) was launched in early 1992 by National Space Development Agency of Japan (NSDAJ). The L-band gave choices of analyzing due to its capability of passing through the land surface.

Started in 1995, the Canadian Space Agency (CSA) operated Radarsat, designed to perform ice coverage of Arctic. It was a C-band SAR system with HH polarization.

It has several imaging modes with different incidence angle, resolution, and swath

(25)

width providing various data for observations.

The ENVISAT (European Space Agency Environment satellite) satellite was launched put in the same orbit as ERS-1/ERS-2 in February 2002 by ESA. It also has various imaging modes with different incidence angle, resolution, and swath width. In addition, the satellite has multi-polarization modes.

The Phased Array type L-band Synthetic Aperture Radar (PALSAR) mounted on the ALOS (Advanced Land Observing Satellite) launched in January 2006 is the Japanese second-generation spaceborne SAR following the JERS-1. It is the first L-band SAR system with full-polarization.

1.3 Geodetic Data and Methodology

The SAR images used in this study were took by ERS-1 and ERS-2 launched by ESA. For the Tainan area in Chapter 3, 33 SAR images have been chosen during January 1996 to August 1999 from track 232, frame 3141. For the Yunlin area in chapter 4, 33 SAR images have been chosen during January 1996 to June 1999 from track 232, frame 3123 (Figure 1.3). In addition, the 'Precise Orbit' data are downloaded from the Delft Institute of Earth Oriented Space Research (DEOS).

The procedure of differential InSAR is processed by the software DIAPASON V4.0 (Differential Interferometric Automated Process Applied to Survey of Nature) developed by Centre National D'Etudes Spatiales (CNES), France, and now maintained by Altamira Information, Spain. The phase unwrapping procedure for Chapter 4 is processed by the software SNAPHU (Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping) from Stanford University, USA. Furthermore, the

(26)

PSInSAR processes are compiled by MATLAB 7.1 language.

Figure 1.3. SAR images of track 232 and frames 3141 and 3123 are used in this study.

(27)

1.4 Thesis Roadmap

The aim of this study is to apply persistent scatterers InSAR (PSInSAR) technique on Tainan Tableland and Yunlin County caused by tectonic activity and land subsidence, respectively. First I introduce the principles of PSInSAR technique that have been published in scientific journals, and adopt these ideas to develop the procedurals of PSInSAR processes of our own.

Chapter 2 provides a brief overview of InSAR and then introduces the persistent scatterer InSAR methods used in this study. Chapter 3 is primarily concerned with the exactitude of procedures and concepts for isolating and extracting the deformation signal from persistent scatterer pixels. I apply this method on datasets from Tainan Tableland as a test site, and validate the technique by comparing results with sequence stratigraphy dating, precipitation, groundwater level data and seismicities. Chapter 4 demonstrates both conventional DInSAR and PSInSAR techniques applied in a land subsidence area, the Yunlin County, where dense vegetation causes significant perturbation of measuring by conventional DInSAR in wet seasons. As a consequence, I apply PSInSAR in order to obtain the deformation evolution close to the Yunlin section of THSR. Chapter 5 provides a summary of the thesis and suggestions for future work.

(28)

(29)

Chapter 2 Synthetic Aperture Radar Interferometry

2.1 Synthetic Aperture Radar

Radio Detection and Range (Radar) is a technique to determine the distance between transmitter and target by emitting electromagnetic waves of microwave frequency and receiving the reflected signals. The usage of imaging the Earth by radar arose in late 1950s, but the first time came into scientific use by Seasat satellite built by JPL/NASA in 1987, which was the first task to analyse Earth by radar images [Massonnet and Feigl, 1998; ].

Real Aperture Radar (RAR) is the earliest radar image system, which was focused on measuring Earth surface. The resolution in the flight direction was obtained by using a physically long antenna. Figure 2.1 shows a Side-Looking RAR (SLAR) mounted on an aircraft or satellite to illuminate the Earth surface; hence it could

(30)

receive the information on ground range by delay time. As the spacecraft traveled with a constant velocity along the flight direction, images on nearby area can be received continually. Nevertheless, the resolution in the flight direction (Ra) with constant wavelength λ was decided by the length of antenna (La).

cos hsat

Ra La La r

λ λ

= ⋅ θ =

⋅ (2- 1)

where h is the height of satellite. Taking ERS1 for example, 25 meters azimuth _sat resolution needs an antenna for 440 meters. It is unfeasible to carry an antenna with such size when traveling in the space.

Figure 2.1. Imaging geometry of a RAR acquisition. The pulses illuminate a swath parallel to the satellite track. The resolution of footprint depends on antenna bandwidth and ground range. θ is the looking angle, r represents the range distance (modified from [Hanssen, 2001]).

(31)

In order to improve the azimuth resolution, fixed antenna is no longer used in an aircraft or a satellite. Thus, the Synthetic Aperture Radar (SAR) was developed. SAR system mounted on a moving spacecraft from position 1 to 11 around the Earth with a series of microwave pulses to illuminate the Earth surface (Figure 2.2). As the moving antenna passes, all return signals with different Doppler frequency shift from multiple locations were used to synthesize a large antenna. SAR system significantly improves the resolution of point target in both the cross-track (range) and along-track (azimuth) direction by focusing the raw radar echoes [Elachi, 1988; Curlander and McDonough, 1991].

Figure 2.2. SAR pulse transmitting and receiving mode. Here, a pulse emitted at position 1 is received between position 10 and 11. The time interval between the pulse transmission is the inverse pulse repetition frequency (PRF) [Hanssen, 2001].

(32)

2.2 Interferometry Synthetic Aperture Radar

A single SAR image contains the information of phase and amplitude (Figure 2.3).

By comparing the phase of two different images with different position (two antennas mounted on the same spacecraft received at the same time) or different time (one antenna received at different time), which can acquire the Z-direction representing the height of surface target. After coregistering and calculating the different phase of received echoes, the interferogram could be reproduced (Figure 2.4).

Figure 2.3. Amplitude image of ERS SAR in Northern Taiwan represents by gray color. The areas with brighter signals mean the stronger reflective signals.

This technique is called Interferometric Synthetic Aperture Radar (InSAR).

InSAR technique was developed after the Second World War and was first applied in earth-based observations of Venus to separate ambiguous echoes from the northern and

(33)

southern hemispheres [Rogers and Ingalls, 1969]. Subsequently, the elevation data of the Moon [Zisk, 1972] and the Venus [Rumsey et al., 1974] were acquired by InSAR.

Not until the late 1980s, scientists started using InSAR to measure the elevation of Earth [Zebker and Goldstein, 1986].

Figure 2.4. Radar interferogram of a region near San Francisco, California, representing measured amplitude (brightness of each point) and phase (color) [Zebker and Goldstein, 1986]

Shuttle Radar Topography Mission (SRTM) was performed between 11 and 23 February 2000 by using a single-pass configuration with a fixed 60 m boom (baseline) to carry two radar antennas, which produced a precise global Digital Elevation Model (DEM) by this SAR interferometry (Figure 2.5). Furthermore, with ERS-1 and ERS-2

(34)

launched into identical orbits, but with the ground tracks separated by 24 hours during the tandem operation, the repeat-pass interferometry became feasible. Moreover, these two missions maintained in the 35 day repeat orbit keep providing abundant data to continuously monitor the crust deformation.

Figure 2.5. Imaging geometry for repeat- and single-pass radar interferometry [Knöpfle et al., 1998].

In an interferogram, the phase was determined by the difference between the phase values of the signal in the two SAR images (master and slave). Assuming that no changes in surface properties, surface elevation, or atmospheric path delay occurred during the short time interval between the two satellite passes. The interferogram phase φ is directly related to the difference in path ∆ according to the following formula r (Figure 2.6):

4 4

( )

s m π rs rm π r

φ φ φ

λ λ

= − = − − = − ⋅△ (2- 2)

(35)

Where m is the position of “master” and s is the position of “slave”; r is the range distance.

Figure 2.6. Imaging geometry for satellite radar interferometry. B is spatial baseline between two SAR acquisitions, and B_⊥ refers to the perpendicular baseline.

α is the angle between the baseline vector and the horizontal, θ is the looking angle, and θ_i is the incidence angle [Hooper, 2006].

α

(36)

2.3 Differential Interferometry Synthetic Aperture Radar

Supposing a target with height h_tar above the reference surface moved during two passes, the phase in the interferogram contains not only topographic contributions (Figure 2.7) but also surface displacements (Figure 2.8). For the topographic contributions, the interferometric phase can be written as:

( ) ( )

4 4 4

cos cos

sin htar

r B B

r

π π π

δφ θ α δθ θ α

λ λ λ θ

= − ∂△ = − − ≈ − − (2- 3)

In fact, the sensitivity of topography is controlled by the position of satellites. The baseline between two flights decides the value of ambiguity high [Hanssen, 2001], which means the elevation of a complete fringe (i.e. –π to π) when δφ is equal to 2π:

2

sin

a 2 h h r

π B

λ θ

⊥

= = (2- 4)

where

( )

cos

B_⊥ =B θ ω− (2- 5)

From equation (2-4), the bigger baseline between two flights causes the more sensitive of topography, and increases the phase error by target elevation. It is hard to create clear fringes in interferogram when two flights with big baseline value.

Therefore, baseline value is the most important parameter when picking up two SAR image for InSAR.

(37)

Figure 2.7. Interferometric configuration. Two SAR sensors are separated by a baseline B. The changes in looking angle δθ can be determined from the interferometric phase. Using the range change and the height of satellite it is possible to determine the height of point P (modified from [Hanssen, 2001]).

To conclude this evaluation, taking the influences of topography h_tar and surface displacement h_dis into account. The phase value ∂ to satellite can be written as: φ

1

4

dis sin tar

h B h

r φ π

λ θ

 ⊥ 

∂ = −  + 

  ^{(2- 6)}

For the sake of simplicity, the phase delay due to differential atmosphere conditions and random noises were ignored. The interferogram phase can be defined as the difference between the measured phase φ and the expected phase ϑ from the geometry, so

φ φ ϑ

∂ = − (2- 7)

(38)

where

4π Bsin( )

ϑ θ α

= λ − (2- 8)

According to equations (2-6), (2-7) and (2-8), the measured phase at target would be:

1

4 sin( )

sin( )

dis tar

B h B h

R

φ π θ α

λ θ α

 ⊥ 

=  − − − −  ^{(2- 9)}

Therefore, the displacement along the Line of Sight (LOS) direction by removing the affection of topography by using two interferograms with different time interval could be got. This technique called Differential InSAR (DInSAR) has been proved its great usages with centimeter precision in applying on measuring coseismic deformation [e.g., Massonnet et al., 1993], glacier and ice motion [e.g., Goldstein et al., 1993], and volcanic eruption [e.g., Amelung et al., 2000].

Figure 2.8. Imaging geometry for differential interferometry. Using interferograms to differentiate the DEM or that from tandem pair interferograms, the target P with a displacement h_dis during two passes can be determined.

(39)

DInSAR is based on using two interferograms with different time. The first one is called topographic pair (e.g., only with topography contribution), and the other one is called deformation pair (e.g., with topography and deformation contribution). Thus, as long as subtracting two interferograms, the interferogram only with deformation can be got. There are three ways to construct a differential interferogram:

Two-pass Differential Interferometry

The two-pass method [Massonnet et al., 1993] uses an existed DEM that is converted to radar coordinates for the topographic pair, and chooses an interferogram interested with deformation for the deformation pair. After subtracting the simulated DEM from the interferogram, deformation information is propagated.

Three-pass Differential Interferometry

The three-pass method [Zebker et al., 1994] uses three SAR acquisitions. First, chooses two suitable SAR images S1 and S2 to create the so-called topographic pair.

This pair is assumed to have no deformation. The deformation occurs between A1 and S2' to create the deformation pair (Figure 2.9). After subtracting the topographic pair from the deformation pair, deformation information is propagated. The advantage of this method is the differential interferogram can be created even without an external DEM. Nevertheless, this method can only be applied when both the topographic pair and the deformation pair have a joint image. For some baseline limitation, it is hard to find the suitable pair.

(40)

Figure 2.9. Three-pass differential interferometry configuration and radar imaging geometry. The solid lines show that radar signal paths for the first interferogram pair formed by antennas at S1 and S2. Dashed lines show signal path for second interferogram acquired over the same site with antennas located at S1 and S2' [Zebker et al., 1994].

Four-pass Differential Interferometry

According to the limitation of the three-pass method, four-pass method uses four SAR acquisitions. In this case, the topographic pair and the deformation pair are independent, i.e. they share no common SAR acquisition. Here, there is no deformation between S1 and S2 to create the topographic pair, and the deformation occurs between S3 and S4 to create the deformation pair. After subtracting the topographic pair from the deformation pair, deformation information is propagated. In

(41)

this method, the deformation occurs between S2 and S3 may cause big error of measuring the real deformation.

Figure 2.10. Four-pass differential interferometry configuration and radar imaging geometry.

(42)

2.4 Persistent Scatterers InSAR

2.4.1 Introduction

The Persistent Scatterers (PS) InSAR is an advanced technique in comparison with conventional InSAR technique, which addresses to overcome the problems of decorrelation for generating a time series of phase changes without atmospheric and DEM residual effects. This technique has been developed in the late 1990s by A.

Ferretti, F. Rocca, and C. Prati of the Technical University of Milan (POLIMI). The first algorithm to find out the PS pixels was brought up by Ferretti et al. [2000; 2001], and trademarked it as the ''Permanent Scatterers technique^TM''. After that, similar processing algorithms have since been developed by Colesanti et al. [2003], and Kampes [2005]. Besides, the SBAS (Small Baseline Subset) technique developed by Berardino et al. [2002] and StaMPS (Stanford Method for Persistent Scatterers) developed by Hopper et al. [2004] are also the same idea of the PSInSAR technique with different names.

In a radar image, a signal of a resolution element was decided by the average signal of all individual scatterers within an area. Therefore, there will be either construction or destruction contribution to phase and amplitude in a resolution element developed a lot of noises in a pixel (Figure 2.11a). That is why that good conventional InSAR result can not be revealed in an area with dense vegetation. Persistent scatterers technique uses the largest contributor signals (i.e. bridges, buildings. Figure 2.11b) as the signals of a resolution element which were called persistent scatterers. In the time series, these PS points are stable enough could refer the information in the whole area.

According to the conditions between neighborhood within an area are similar, using the PS point to represent the information of a resolution element can reduce the

(43)

uncertain signal of back scatters and continuously get the phase changes of any PS point easily in the time series.

Figure 2.11. Phase simulated for (a) a distributed scatterer pixel and (b) a persistent scatterer pixel. The cartoons above represent the scatterers contributing to the phase, and the plots below show the phase values [Hooper, 2006].

The main steps of the PSInSAR technique can be summarized as following steps [Ferretti et al., 2000, 2001; Mora et al., 2003].

1) Formation of Differential Interferograms: Giving N+1 SAR images (same area) to form N interferograms with the same master image (detail in 2.4.2).

2) Selection of Stable Point Target: Calculating the stability of each pixel to select PS candidates using analysis steps in 2.4.3.

3) Extraction of Linear Deformation and DEM Residual: Using the information of temporal baseline and perpendicular spatial baseline to generate a deformation

(44)

model of the study area. Then, the linear deformation and DEM residual over each PS pixel is extracted in comparison with generated models and the real interferogram patterns (detail in 2.4.4).

4) Extraction of Nonlinear Deformation and Atmospheric Artifacts: Due to different phase characters in atmosphere and nonlinear motion, the spatial and temporal filters are used to separate different phase contributions (detail in 2.4.5).

5) Deformation of Temporal Evolution Construction: Finally, integrating the linear and nonlinear deformation; the deformation evolution can be obtained.

2.4.2 Choice of “Master” Image

The first step is to choose the one SAR image that minimizes the sum decorrelation, i.e., maximizes the sum correlation, of all the interferograms in the data set. The correlation depends on four terms: temporal baseline ( )T , perpendicular spatial baseline

( )

^B⊥ , Doppler centroid frequency baseline (F_DC) and thermal noise [Zebker and Villasenor, 1992]:

1 ( ) 1 ( ) 1 ( )

total temporal spatial doppler thermal

DC

thermal

c c c

DC

F T B

f f f

T B F

ρ ρ ρ ρ ρ

ρ

⊥

 

 

     

=

≈ − − −

i i i

(2- 10)

where ρ denotes correlation and superscript c denotes the critical value which means the limit of producing a useable interferogram. For the ERS data only, the T = 5 yeras, B_⊥^c = 1100 meters, F_DC^c = 1380 Hz, and ρ_thermal is a constant value [Hooper, 2006].

(45)

2.4.3 Persistent Scatterers Identification

In the processes of PSInSAR, the algorithm is established on PS points selected.

Thus, the selections of PS points are quite important. There are two common methods of selecting PS points, one based on pixel amplitude stability and the other on spatial coherence.

The Amplitude Dispersion Index

( )

DA is defined by Ferretti et al. [2001] as:

A A

A

D σ

≡ µ (2- 11)

where σ_A and µ_A are respectively the standard deviation and mean of a series of amplitude values in each pixel. First, providing a complex reflectivity g and a complex circular Gaussian noise n with real

( )

nR and imaginary components

( )

n . I

Because the distribution of the amplitude value A is given by the Rice distribution:

( )

2 0 2 ⁽^a² ^g²⁾^{/ 2} ⁿ² 0

A

n n

a ag

f a I e ^σ a

σ σ

− +

 

=   >

  ^{(2- 12)}

where I is the modified Bessel function. For the high 0 Signal to Noise Ratio (SNR, / _n 4

g σ > ), fA approaches a Guass distribution. Provided that σ_n _≪ g , then:

A nR nI

σ ≃σ =σ (2- 13)

since the modulus is primarily affected by the noise component parallel tog n

( )

R . The phase distribution

( )

σv can then be estimated starting from the amplitude dispersion:

nI A

v A

A

g D

σ σ

σ ^≃ ^≃ µ ⁼ ^{(2- 14)}

The D is then a measure of phase stability, at least for high SNR values (low_A D ). _A Typically, D < 0.25 is usually chosen as the threshold [Ferretti et al., 2001]. _A

(46)

Figure 2.12. Numerical simulation results for D and phase stability. Small values _A of D are good estimates of the phase dispersion [Ferretti et al., 2001]. _A

For spatial coherence method, following the idea of PSInSAR, the correlation threshold would be the easiest approach. First, two SAR images I x y1

(

,

)

and

( )

2 ,

I x y can be written as:

( ) ( )

⁽ ⁾

( ) ( )

⁽ ⁾

1

2

,

1 1

,

2 2

, ,

j x y

I x y A x y e I x y A x y e

φ

=

= (2- 15)

( )

1 ,

A x y , A x y are the amplitude values, and 2

(

,

)

φ1

(

x y,

)

, φ2

(

x y,

)

are the phase values, then the coherence can be expressed as:

( ) ( ) ( )

⁽ ⁾ ⁽ ⁾

( ) ( )

1 , 2 ,

1 2

2 2

1 2

, ,

,

, ,

j x y j x y

A x y A x y e x y

A x y A x y

φ φ

γ

−

=

∑

∑ ∑

^{(2- 16)}

where ∑ means average within a window. As long as the pixel exhibits

(47)

coherence always greater than the threshold in all data set, provide that the phase stable enough. All these pixels with this character will be selected as PS candidates. Due to the principle of coherence, the value of coherence is based on the estimation window.

The larger the window dimension, the higher the estimator accuracy, but the lower the resolution. Besides, the DEM residual cause significant contribution to phase.

Supposing a SAR pair with big spatial baseline, the coherence value will be low even the pixel with strong stability. Thus, the suitable values of window dimension and coherence threshold are needed to be finding out.

After picking up the PS candidates, only these pixels will be analyzed in the following steps. It is important to point out that the PS candidates which selected by the above-mentioned algorithms were the preliminary result, they are not the final PS points. In the following steps, the points with bad correlation to the real deformation will be rejected to get the real PS points.

2.4.4 Extraction of Linear Deformation

When producing a DInSAR result by two SAR images, the phase can be written as the sum of four terms:

mov topo error atm noise

φ

diff =

φ

+

φ

+

φ

+

φ

(2- 17)

where

φmov: Phase change due to movement of the pixel in the range distance;

topo error

φ : Residual topographic phase due to misfit with the DEM;

φatm: Phase due to different atmospheric condition between passes;

noise

φ : Noise term due to temporal and spatial decorrelation and thermal.

(48)

Here, the movement term

(

φmov

)

has two contributions: linear movement and nonlinear movement only correlate with velocity, and the DEM error is proportional to the perpendicular spatial baseline:

4

mov linear nonlinear π v T nonlinear

φ φ φ φ

= + = λ ⋅ ⋅ + (2- 18)

4

topo error sin

B r

π ε

φ λ θ

⊥⋅

= ⋅

⋅ (2- 19)

where

λ: System wavelength (5.6 cm for ERS);

v: Constant velocity of the linear displacement;

T : Temporal baseline between two SAR acquisitions;

r: Range increment between two SAR acquisitions;

ε : DEM error;

θ: Incidence angle.

Because the phase is wrapped, the phase of individual pixel has an infinite solution. In order to calculate the best solution, the usage of Delauney triangulation is joining to connect all the discrete points for controlling each two points with each other.

This kind of triangulations relates all the neighboring pixels of irregularly girded data generating non-overlapped triangles. In Figure 2.13, all the discrete points are connected to form a network to calculate the best solution.

The phase difference between two points ( ,x y_a _a) and ( ,x y_b _b) in Figure 2.13 can be expressed as:

mov atm noise

diff ε

δφ

=

δφ

+

δφ

+

δφ

+

δφ

(2- 20)

(49)

For the whole data set, it becomes as:

[ ]

( , , , , ) 4 ( , ) ( , )

( )

4 ( , ) ( , )

( ) sin( )

( , , ) ( , , ) ( , , ) ( , , ) ( , , ) ( , , )

diff a a b b i i a a b b

i

a a b b

i i

a a i b b i

x y x y T T v x y v x y

B T x y x y

r T

x y T x y T

a x y T a x y T n x y T n x y T δφ π

λ

π ε ε

λ θ

β β

= ⋅ ⋅ −

+ ⋅ ⋅ −

⋅

+ −

(2- 21)

where

x, y : Position coordinates of the pixel;

Ti: Temporal baseline of the i th interferogram;

β: Nonlinear component of velocity;

a: Atmospheric artifacts;

n: Noise.

Figure 2.13. Example of Delauney triangulation. ( , )x y are the position coordinates of PS point a and b [Mora et al., 2003].

(50)

Here, a limitation that atmospheric perturbations are the same in the neighboring pixels in a small area was added. The maximum distance allowed to connect two separate pixels is limited to approximately 1 kilometer. For this reason, first reject the connection with distance longer than the limitation distance.

Using the information of temporal and spatial baseline that already known and assuming a velocity and DEM error range in (2-18) and (2-19), a model can be assumed as:

[ ]

( , , , , ) 4 ( , ) ( , )

( )

4 ( , ) ( , )

( ) sin( )

model a a b b i i model a a model b b

i

model a a model b b

i i

x y x y T T v x y v x y

B T x y x y

r T δφ π

λ

π ε ε

λ θ

= ⋅ ⋅ −

+ ⋅ ⋅ −

⋅ (2- 22)

Then, Using the following model coherence function to perform maximization

model

γ :

0

1 exp ( ( ) ( ))

N

model diff i model i

i

j T T

γ N δφ δφ

=

 

= ∑  ⋅ −  ^{(2- 23)}

where N is the number of interferograms. This function is equal to 1 when the adjustment to data is perfect, and zero in case of total decorrelation. Once the maximization process has been done for each connection, the result is the following set of velocity and DEM error:

[ ]

( , , , ) ( , ) ( , )

estimated a a b b model a a model b b maximize

v x y x y v x y v x y

x y x y x y x y

γ

ε ε ε γ

∆ = −

∆ = − ^{(2- 24)}

At this point, a new quality test is performed and all the connections with coherence below a threshold are rejected. Good results have been obtained with coherence thresholds larger than 0.7 [Mora et al., 2003]. As a consequence of the quality test, some pixels left without connection will be eliminated. Besides, the

(51)

velocity and DEM error in each connection can be written as:

( ) ( ) ( ) ( )

, ,

estimated a a b b

v v x y v x y

x y x y

ε ε ε

∆ = −

∆ = − ^{(2- 25)}

Then, A least mean square process is necessary to obtain absolute values for each point.

Figure 2.14 shows a detailed layout of linear deformation and DEM error extraction.

Figure 2.14. Layout of linear deformation and DEM error extraction.

Triangulation Pixels Selection Differential

Interferograms

Coherence Images

Linear Model Adjustment

Deformation and DEM Error Maps

Over PS Areas LMS

(52)

2.4.5 Extraction of Nonlinear Deformation and Atmospheric Artifacts

After finishing the linear velocity and DEM error, it is possible to obtain the nonlinear component to complete the whole deformation in interested area. At first, the residual phases obtained by subtracting the previously estimated linear deformation and DEM error from the original differential interferometric phases. Then, the phase residual can be expressed as:

( , , ) ( , , )

residual i nonlinear i

atmos i noise i

x y T x y T

δφ δφ

=

+ + (2- 26)

where two important terms (the nonlinear contribution δφ_nonlinear and the atmospheric artifacts contribution δφ_atmos) can be isolated taking advantage of their different frequency characteristics in space and time [Ferretti et al., 2000]. The mean value of the residual phases (δφ_residual) is an estimation of the atmospheric phase contribution of the master image in each PS point (common to all the differential interferograms, said correlate in temporal). Beside, the nonlinear motion contribution was also considered to reveal a high correlation in temporal. Thus, the estimation of the nonlinear motion contribution should be considered as the low pass component of this function:

( ) _residual _

nonlinear residual Ti LP Time

δφ =^δφ −δφ ^ ^{(2- 27)}

Furthermore, atmospheric retardation is spatial correlated over distance of a specified length scale, said 1 kilometer. Atmospheric perturbations are considered as a low spatial frequency signal in each interferogram, and a random phase in time for each pixel (for slave images). Thus, the atmospheric estimation is then carried out by smoothing spatially the high pass filtered residual phase time series:

[ ]

_{_}

_ _

( ) _residual

atmos residual i HP Time LP Space LP Space

δφ =^ δφ T ^ + ^δφ ^ (2- 28)

(53)

where the first term is an estimation of the atmospheric phase contribution superimposed on the single SAR image acquired at time T_i. The filtered atmospheric component relative to the master acquisition is then added back to get the atmospheric artifacts of each differential interferogram. Figure 2.15 shows a detailed layout of nonlinear deformation and atmospheric artifacts extraction.

Figure 2.15. Layout of nonlinear deformation and atmospheric artifacts extraction.

Differential Interferograms

Phase Model Generation Residual

Phase

Deformation and DEM Error Maps

Over PS area

Temporal High Pass Filtering

Mean

Atmospheric Artifacts

+ _

Spatial Low Pass Filtering

Temporal Low Pass Filtering

Nonlinear Deformation

_ +

(54)

(55)

Chapter 3 Active Deformation in Southwestern Taiwan

In this chapter, I present an application of PSInSAR technique based on SAR images deduced from 1996-1999 time periods in and around Tainan Tableland of southwestern Taiwan, which is an active fault-related folding near the deformation front. First, I want to characterize the deformation pattern due to the evolution of this active structure. Second, the relationship between nonlinear behavior demonstrated by PSInSAR results and precipitation, groundwater level, and microseismicities will be investigated.

Here, I adopt the idea from Ferretti et al. [2000; 2001] and Mora et al. [2003] to develop the procedures of PSInSAR processing by Matlab for this study. In order to ensure the reliability of our programs, I choose the Tainan area as a test site where it had been well studied by using conventional DInSAR technique [Fruneau et al., 2001;

Huang et al., 2006] and the geodetic measurements [Rau et al., 2004]. Thus, the Tainan area is expected to be an ideal urbanized site to validate the PSInSAR results.

(56)

3.1 Introduction

Taiwan is located in an active collision boundary where the EP has collided with the PSP since late Miocene (Figure 3.1). Because of the oblique convergence between the two plates, strong crustal deformation is demonstrated by frequent seismicities and active structures [e.g. Hu et al., 2001]. The Tainan Plain of southwestern Taiwan is incorporated into active deformation by a regional fold-and-thrust belt associated with the convergent plate margin of the Luzon Arc and the Eurasian continent.

The Tainan city surrounded Tainan Tableland, which is considered to be located close to the deformation front in southwestern Taiwan [Sun, 1964]. Previous DInSAR research from Fruneau et al. [2000] revealed a ground motion of 28 mm along the radar line of sight towards the satellite during the period 1996 to 1998. Also, Ho [2006]

showed a decreasing SRD from the southern Tainan Tableland (12 mm/yr) to the northern Tainan Tableland (10 mm/yr). Furthermore, Huang [2006] predicted an average SRD rate of 12.5 mm/yr during the period 1996 to 2001 on the Tainan Tableland which is quite consistent with the precise leveling data measured by Rau et al. [2004] (Figure 3.2). The precise leveling data reveals maximum uplift rates on the Tainan Tableland of about 12 mm/yr relative to TN01 (Figure 3.3). However, a borehole data near TN01 shows a long term subsidence rate is about 6.6 mm/yr [Zhou, 2007]. Thus, the uplift rate predicted by precise leveling by using TN01 as a reference station should be overestimated. In addition, Chen and Lin. [2000] demonstrated that the long term uplift rate on the Tainan Tableland is about 5 mm/yr and revealed an eastward increasing in uplift rates based on radiocarbon ages and relative sea-level curves of late Holocene. These entire results figure out the Tainan Tableland is a significant growing structure.

(57)

Figure 3.1. Geodynamic framework and structural map of Taiwan [after Ching et al., 2007]. Red arrow indicates the relative velocity across Taiwan Island with 82 mm/yr along N309 ° E from the PSP [Yu et al., 1997]. Pink rectangle indicates the study area, the Tainan tableland.

(58)

According to these abundant geologic and geodetic observations in Tainan area, the Tainan Tableland is the most suitable study area to validate our PSInSAR technique.

The goal of this case study is focus on PSInSAR technique to obtain the deformationpattern in the Tainan area in comparsion with geodetic measurements to better understand the correlation between the deformation patterns and geological structure in Tainan area.

Figure 3.2. Active deformation in Tainan Tableland revealed by D-InSAR and dislocation model. (A) SRD of interferogram after phase unwrapping on shaded topography, which can indicate the deformation pattern of a pup-up structure. (B) Inferred structural model as a decollément-related pop-up structure [Huang et al., 2006].

(59)

Figure 3.3. (a) Southwestern Taiwan geological structural setting and location of GPS (triangles and circles) and leveling stations (rectangle and circles). (b) Elevation measurement of precise leveling (relative to TN01) project to AB profile. (c) Uplift rate measurement of precise leveling relative to TN01 [Rau et al., 2004].

利用永久散射體雷達干涉技術分析地表變形：以臺灣西南部平原為例

Department of Geosciences College of Science

National Taiwan University Master Thesis

國立臺灣大學碩士學位論文

口試委員會審定書

本論文係童忻君（R95224206）在國立臺灣大學地質科學 所完成之碩士學位論文，於民國九十七年七月九日承下列考試 委員審查通過及口試及格，特此證明

口試委員：

（簽名）

（指導教授）

Abstract

摘要 摘要 摘要 摘要

致謝 致謝 致謝 致謝

Contents

List of Figures

Lists of Tables

Lists of Abbreviations

Chapter 1

Introduction

1.1 Motivation

1.2 Spaceborne Radar History

1.3 Geodetic Data and Methodology

1.4 Thesis Roadmap

Chapter 2

Synthetic Aperture Radar Interferometry

2.1 Synthetic Aperture Radar

2.2 Interferometry Synthetic Aperture Radar

2.3 Differential Interferometry Synthetic Aperture Radar

( ) ( )

( )

2.4 Persistent Scatterers InSAR

2.4.1 Introduction

2.4.2 Choice of “Master” Image

( )

ρ ρ ρ ρ ρ

ρ

2.4.3 Persistent Scatterers Identification

( )

( )

( )

( )

( )

( )

(

)

( )

( ) ( )

( ) ( )

( )

(

)

(

)

(

)

( ) ( ) ( )

( ) ( )

∑

∑ ∑

2.4.4 Extraction of Linear Deformation

φ

φ

φ

φ

φ

(

)

δφ

δφ

δφ

δφ

δφ

[ ]

[ ]

[ ]

[ ]

[ ]

[ ]

[ ]

[ ]

[ ]

本論文係童忻君（R95224206）在國立臺灣大學地質科學所完成之碩士學位論文，於民國九十七年七月九日承下列考試委員審查通過及口試及格，特此證明

摘要摘要摘要摘要

致謝致謝致謝致謝