Coseismic Deformation of Chi-Chi Earthquake as Detected by Differential Synthetic Aperture Radar Interferometry and GPS Data

1 Department of Civil Engineering, National Chiao-Tung University, Hsinchu, Taiwan, ROC * Corresponding author address: Dr. Tian-Yuan Shih, Department of Civil Engineering, National ChiaoTung University, Hsinchu, Taiwan, ROC; E-mail: tyshih@mail.nctu.edu.tw A rupture in the Chelungpu fault caused an M w 7.6 earthquake on 21 September 1999 near Chi-Chi in central Taiwan. This earthquake was the most destructive experienced in Taiwan for the past century along this fault. In this study, we examined the earthquake-induced surface deformation pattern using differential synthetic aperture radar interferometry (D-InSAR) combined with global positioning system (GPS) data regarding the footwall of the Chelungpu fault. Six synthetic aperture radar (SAR) scenes, approximately 100 × 100 km each, recorded by the European Remote Sensing Satellite 2 (ERS-2), spanning the rupture area, were selected for study. The data were used to generate a high-resolution, wide-area map of displacements in flat or semi-flat areas. Interferograms show radar line contours indicating line-of-sight (LOS) changes corresponding to surface displacements caused by earthquake ruptures. These results were compared to synthetic interferograms generated from GPS data. Displacements shown by GPS data were interpolated onto wide-area maps and transformed to coincide with the radar LOS direction. The resulting coseismic displacement contour map showed a lobed pattern consistent with the precise GPSbased displacement field. Highly accurate vertical displacement was determined using D-InSAR data using the coordinate transform method, while GPS data was effective in showing the horizontal component. Thus, this study confirmed the effectiveness of the D-InSAR method for determining the coseismic deformation caused by the Chi-Chi earthquake at the footwall of the Chelungpu fault. (


INTRODUCTION
Synthetic Aperture Radar Interferometry combines complex images recorded by antennas at different locations or at different times. The resulting interferograms permit the determination of differences in the range direction. This technique has been used widely to produce highly accurate digital elevation model (DEM) and topographic maps (Zebker and Goldstein 1986;Madsen et al. 1993) and to measure displacement fields or terrain motions (Massonnet et al. 1993;Zebker et al. 1994). This technique also has been applied to ocean current measurements, hazard mapping and glacier motion detection (Massonnet and Feigl 1998;Rosen et al. 2000).
A radar sensor above the Earth detects tiny changes on the ground by accurately measuring changes in the time delay, or phase, of a radar echo. D-InSAR and its spatially dense, accurate deformation measurements have advanced studies of the Earth's crust. Importantly, this technique provides a comprehensive view of detectable motion for the entire area affected. Results from such analyses can supplement ground-based measurements taken at a limited number of locations. Because radar satellites make global observations, we can study deformation processes anywhere at little additional expense once a satellite is in operation. The great advantage of the radar interferometer system is that it allows deformation measurements at very fine spatial spacings, creating a visual image of the deformation distribution (Zebker 2000).
A disastrous M w 7.6 earthquake struck central Taiwan at 01:47 (local time) on 21 September 1999. The Seismology Center of the Central Weather Bureau located the epicenter near the town of Chi-Chi (23.85°N, 120.78°E) in Nantou County. To monitor fault activity, it is important to characterize the crustal deformation caused by this quake. Following the earthquake, measurements were taken at 99 GPS stations. These measurements suggested horizontal displacements of 0.1 to 8.5 m at the sites (Chang 2000). Yang et al. (2000) studied the three-dimensional displacements of 285 geodetic control stations using GPS observations collected before and after the earthquake. Yu et al. (2001) applied annually repeated GPS data from 1992 to 1999 to predict the crustal deformation rate; data within 3 months after the main shock were used to estimate the coseismic displacements.
In addition, Suga et al. (2001) used D-InSAR data to analyze deformation caused by the Chi-Chi earthquake. The land displacement pattern extracted from D-InSAR data was generally consistent with results of a GPS survey conducted by researchers Yang et al. (2000). Pathier et al. (2003) corrected interferograms using GPS data to create a precise map of the D-InSAR component of coseismic displacement in the LOS direction. Chang et al. (2004) and Liu et al. (2004) applied the interferometric method to detect earthquake displacement. The interferometric results were precisely examined using GPS data in the LOS direction.
Major earthquake damage results from the vertical and horizontal components of fault motion. The advantage of the D-InSAR technique is that it can be used to characterize coseismic deformation over a larger area than that covered by GPS surveys. However, the inverse transformation from LOS displacement to vertical and horizontal displacement is not unique. In this study, we attempted to obtain the vertical displacement component using the coordinate transform method and the horizontal component using GPS horizontal data. The results con-firmed that ERS/SAR D-InSAR data indeed effectively determined the coseismic deformation created by the Chi-Chi earthquake at the Chelungpu fault footwall.

GEOLOGICAL BACKGROUND
The Chi-Chi earthquake was caused by a rupture of the Chelungpu fault. This earthquake, the most destructive in Taiwan in a century, struck the west-central island on 21 September 1999 (Kao and Chen 2000). For approximately 90 km along the fault, extensive surface ruptures with vertical thrust and left lateral strikes-slip offsets occurred (Angelier et al. 2001;Lee et al. 2002). The vertical offsets averaged approximately 2 m along the southern half of the fault and about 4 m along the northern segment (Yang et al. 2000;Yu et al. 2001). Horizontal offsets of 5 to 8 m occurred along the northern part of the major fault (Yang et al. 2000;Yu et al. 2001). Tens of thousands of buildings collapsed, resulting in 10000 injuries and over 2300 fatalities.
Central Taiwan has several series of N-S trending and east-dipping thrust faults located in the Western Foothills; from east to west, they are Shuangtung fault, Chelungpu fault and Changhua fault, respectively (Fig. 1). The major rupture of the Chelungpu fault is near N-S trending that follows the boundary between the Western Foothills and the Taichung piggyback basin . The earthquake fault of the Chi-Chi earthquake has been attributed to the Chelungpu fault in the western part; however, fault branches east of the Fengyuan area trending nearly NE-SW were also suggested to form a Cholan segment . For simplification, we consider this segment as part of the Chelungpu fault.
The original boundary of the foreland basin was the Shuangtung fault, but due to continuing northwestward collision from the Philippine Sea plate, the foreland basin moved westward, causing the Chelungpu fault and the Changhua fault to develop in sequence. Consequently, the Deformation Front in the central Taiwan is now located in the Changhua fault, west of the Chelungpu fault.
The Chelungpu fault extends along the western front of the Taiwan Mountain belt in central Taiwan. The footwall of the Chelungpu fault is almost flat or semi-flat, including some urban areas. Thus, these areas are suitable for D-InSAR techniques. In contrast, the hanging wall of the Chelungpu fault lies in a steep mountainous area covered by dense vegetation.

SAR Interferometry
When two radar scans have the same viewing angle but different viewing times, a small change in the target position can create a detectable change in the phase of the reflected signal. The phase is exactly proportional to the measured time delay and effective path length of the signal. The path differences of two signals can be determined to sub-wavelength accuracy by observing the phase differences of echoes. Phase differences may be caused by topographic effects, surface displacements and signal variations due to atmospheric effects (Zebker et al. 1997).
In spaceborne SAR interferometry, two images of the same scene are acquired by a single sensor or by two different antennas passing above the area along two successive orbits. In practice, it is quite difficult to ensure that the radar platform returns to exactly the same position for the second observation. However, this additional displacement is invertible in recovering topography from an interferogram.
The interference pattern of two SAR images of the same region acquired from two passes depends on the topography. If movement occurred between the two passes, the total displacement would consist of two parts: δρ ‚ for the real surface motion between observations and ∆ρ ‚ for the apparent motion due to parallax. The law of cosines permits calculation of the total displacement: where the separation of the satellites, or baseline length, is B, the look angle is θ , and the baseline angle with respect to horizontal at the sensor is α . Neglecting second-order terms, we obtain: An interferogram thus contains phase variations dependent on topography and surface deformation. These effects must be separated if we are to interpret the interferograms unambiguously. The topographic phase signature can be removed from an interferogram during processing if a preexisting DEM is available (sometimes referred to as the two-pass method) or if a second interferogram containing topographic signals only is acquired (the three-or four-pass method). Massonnet et al. (1993) initially developed the two-pass method, in which two SAR images are used to form an interferogram and a DEM is used for topographic reference. If a DEM does not exist for the area of interest or is of poor quality, the three-or four-pass method can be used. Gabriel et al. (1989) first proposed this method, and Zebker et al. (1994) subsequently applied it to measure coseismic deformation. This method use three or four images to form two interferometric pairs: one contains only topographic fringes and is used as a topographic reference, the other contains topographic and displacement fringes.

Vertical Displacements
D-InSAR data show the slant range displacement (SRD) of the satellite rather than actual height changes. However, with accurate GPS horizontal displacement data and the geometric relationship, it is possible to obtain vertical height changes. The geometric relationships of GPS, SRD, and real deformation data demonstrated in Fig. 2.
As shown in Fig. 2, the GPS horizontal displacement vector V equals (a, b, 0), and the included angle between the satellite orbit and the northern direction N equals φ . It is necessary to project the GPS vertical displacement vector (a , b , c ) ′ ′ ′ in the satellite direction using the following equation: (3) Next, the GPS horizontal displacement vector is projected to the satellite direction; for the projection of real deformation, the SRD is subtracted along the satellite direction: where, ∆r = SRD, namely, the differential interference pattern. The vertical displacement h is obtained from Equation (5).

Process Workflow
A SAR processor that employs a range-Doppler algorithm and the Stanford interferometry-processing package was used to process the images (Zebker 2000). In the range-Doppler algorithm, raw radar measurements are first match-filtered to remove signal modulation and generate high-resolution signals in the across-track, or range, dimension. Generating an interferogram involves co-registration and matching the slave image in a master frame to subpixel resolution. In addition, the image data for one of the two radar passes are resampled during processing, so similar ground points in both images coincide in location. Figure 3 shows the flowchart applied in this study. After fine co-registration, every image point corresponds to another at the same position. A complex interferogram is constructed by complex multiplication of corresponding pixels in both image pairs. A multilook process is necessary to decrease the noise in the interferogram. The multilook process can be performed simultaneously with the complex multiplications and is often applied to a range-azimuth ratio that yields approximately square pixels, such as 1 : 5 or 2 : 10.
The phase change of the interference fringe pattern comprises the surface height and changes in distance differences between the point positions and radars. Thus, the effect of distance differences must be removed to calculate the surface height. The distance from different earthsurface points to the two satellite positions may vary slightly. This difference must the removed for the interference pattern to accurately reflect surface information. Given that the image range obtained from satellite images is 100 × 100 km, this relatively long-range length is greatly influenced by Earth flattening, which must also be removed.
The terrain effect is removed to calculate earth-surface changes because interference fringe patterns reflect terrain and surface deformations. Removal can be accomplished in two ways: using existing surface information, such as from a DEM, or using radar data to create a third radar image and any optional image to form a second interference fringe pattern. With the linear combination of two interference fringes in coordination with the base length, it is possible to generate differential phase interference patterns wherein the phase value represents the SRD of the satellite.
Phase based interference patterns are caused by terrain deformation effects. Because this information is given as modulo 2π π π to ( ) − , there is an ambiguity problem in calculating the correct integer of phase cycles, which must be added to each phase measurement to obtain the correct value. The solution to this problem is referred to as phase unwrapping. Although the forward problem of wrapping the absolute phase to the −π π to interval is straightforward, the inverse problem is one of the major difficulties in radar interferometry (Goldstein et al. 1988). Variable phase noises, as well as geometric problems, such as foreshortening and layover, are the main reasons many proposed techniques do not perform as desired (Gens 2003). Numerous methods have been proposed and implemented to address this complex issue (Ghiglia and Pritt 1998). In the process of image acquisition, phase residues may occur due to thermal noise or terrain effects, such as shadows, foreshortening and layover. These residues influence the calculation of actual phase differences for point positions in full-width images. Therefore, the method for dealing with these residues and accurately calculating phase differences between point positions is referred to as phase unwrapping.
We used the SNAPHU program designed by the Stanford Radar Interferometry Group, based on the statistical-cost, network-flow algorithm for phase unwrapping proposed by Chen and Zebker (2000). The heights derived from SAR interferometric data sets were calculated geometrically, determined by the means of orbital parameters. Using interferometric orbital parameters, we derived the physical and geometrical relationships between the two phase observations to obtain estimates of topographic height and surface deformation. The radar coordinate system differs substantially from geographic coordinate systems. The surface data generated by radar images can only be converted into familiar geographic coordinates for data comparison and analysis by performing geocoding calculations. Satellite orbits are given in the conventional terrestrial system. Geocoding refers to a transformation of radar coordinates (range-azimuth-height) to coordinates in a convenient geodetic reference system.
The result is a complex image in which phase components constitute the interferogram. The phases are ambiguous to within an integer of 2π cycles. Precise estimates of the interferometric baseline can be used to separate topographic effects from the deformation phases.
In this study, two interferograms were used to form difference interferograms: one that spanned a period before and after the earthquake (6 May -28 October 1999) and another prior to the earthquake (6 -7 March 1996). The phase of the non-earthquake interferogram was unwrapped using the SNAPHU phase unwrapping algorithm proposed by Chen and Zebker (2000). The algorithm casts phase unwrapping as a maximum a posteriori probability estimation problem in which the objective is to find the most likely unwrapped solution given the wrapped phase as well as other observable quantities, such as image intensity and interferogram coherence.

Data Selection
In Taiwan, the use of SAR interferometry for ground motion measurements is difficult because of coherence problems created by extensive tropical vegetation. The ERS-2 C-band SAR with 35-day repeat cycle produced six images covering the earthquake area close to 21 September 1999; these images were recorded on 6 May, 10 June, 15 July, 19 August, 23 September, and 28 October 1999. Orbit reconstruction provided by precise Delft orbit information enabled us to determine the baselines for the pairs.
Unfortunately, most of the ERS-2 images close to 21 September had poor baselines, suggesting that it was necessary to account for longer time intervals before and after the Chi-Chi earthquake. For the interferograms, we chose pairs with the shortest possible baseline to reduce topographic effects that may have obscured surface deformation. The shortest baseline was 38 m for the 6 May -28 October pair; this pair spanned 175 days. Other pairs had longer baselines.
The four-pass method was selected to detect Chi-Chi earthquake deformation because the pairs before the earthquake had longer baselines. In choosing interferometric pairs for topographic correction, shorter time interval pairs are preferred because these pairs display less surface deformation. The 6 -7 March tandem pair in 1996 over the same area was chosen; this pair was separated by only 1 day and had a baseline of 78 m.

Interferograms
The two interferograms were co-registered before calculating the difference and forming a differential interferogram. To remove the topographical signature from deformation interferograms, the tandem pair was subtracted from the deformation pair for the formation of a differential interferogram. The interferogram was then rectified. Figure 4 shows the geocoded differential interferograms produced by the four-pass method. The interferograms reveal the phase variations associated with the Chi-Chi earthquake surface deformation. The fringe pattern shows 11 clear fringes between the west coast and the fringe center, which corresponds to a deformation of about 2.8 cm on the footwall of the fault line. No fringes are evident on the hanging wall of the image; this hanging wall of the fault had low coherence due to the mostly mountainous terrain covered with dense vegetation. Displacement was also larger on the hanging wall, which contributed to the lower surface coherence.

GPS Data
After the earthquake, the Ministry of the Interior took measurements at GPS control points in affected areas as part of a cadastral survey. We selected 276 control points within our image area and analyzed the data. The accuracy of GPS data in horizontal components was 1 -2 cm, and in vertical component was 6 -7 cm. The horizontal displacement was then calculated: the maximum was 9.23 m, and the average 1.41 m. Large deformation was found between the Chelungpu fault and the Shuangtung fault. In general, movement on the hanging wall of the Chelungpu fault was greater than that on the footwall. Points on the hanging wall showed larger deformation, with the direction of deformation slightly north to northwest as illustrated in Fig. 5. Deformation on the footwall was rather small and oriented towards the southeast.

Fig. 5. Coseismic displacement vectors as measured by GPS data in central
Taiwan. The vector directions are shown in the horizontal direction of motion. Large deformation was found between the Chelungpu fault and the Shuangtung fault. In general, movement on the hanging wall of the Chelungpu fault was greater than that on the footwall. Points in the east showed larger deformation, with the direction of deformation slightly north to northwest. Deformation on the footwall was rather small and oriented towards the southeast.
GPS point data were interpolated onto the grid image and color shaded to illustrate the horizontal magnitude of deformation, so the overall deformation diagram could be displayed as a contour image, as shown in Fig. 6. Maximal deformation is shown in the upper right corner of the image. The magnitude gradually decreases with distance.
To compare effectively the deformation data obtained by GPS and D-InSAR, we changed the coordinates of the GPS data to match the radar coordinates for D-InSAR and then transformed the GPS deformation data into deformation in the range direction and thereby simulate the interferogram. Figure 7 is a flowchart of this transformation processing.  As noted, the radar technique is sensitive to the LOS component of motion. For GPS motion vectors in the direction of the ground projection of the radar sensor boresight, this component is the vector from the sensor to a point on the Earth's surface. For the radar measurements, because the LOS direction is not in the plane defined by the local Earth surface, the equivalent horizontal surface motion to yield the observed SRD is which relates the horizontal displacement ∆y to the SRD and the incidence angle θ inc . Figure 8 shows synthetic interferograms obtained from GPS data for horizontal displacement. The fringes are contour lines of equal ground displacement along the satellite LOS direction. Unlike the interferograms from differential radar processing, which lack some information in low coherence areas, simulation of GPS data produces region-wide information. There are clear fringes on the hanging wall. The Chelungpu fault approximately marks the boundary between the western footwall small deformation and the hanging wall large deformation. The change is greater on the hanging wall than on the footwall, and largest deformation in the region is concentrated on the area between Chelungpu fault and Shuangtung fault. The real fringes, in response to changes in the river area, exhibit discontinuity trend with plane area. The simulation fringe from the GPS data is a continuous curve in the river area. The different patterns of real fringes and simulation fringes can also be found in areas of Nantou and Yuanlin, where the changes are more intensive trends. The two areas are special cases in which the changes may have produced differences due to geological conditions or other causes of deformation. The D-InSAR and GPS methods showed similar deformation trends in the SRD direction and for fringe patterns. However, the D-InSAR data showed 11 fringes, but the GPS data showed only 8 fringes, corresponding to a difference in deformation of about 8.4 cm (2.8 cm per fringe).

Comparison of D-InSAR and GPS Data for Vertical Displacement
To compare the accuracy of vertical displacement calculated using differential interferometry, we selected GPS data from Yu et al. (2001) measured by the Satellite Survey Division and Land Survey Bureau, Ministry of the Interior, the Central Geological Survey, and the Institute of Earth Sciences, Academia Sinica. The accuracy of 28 GPS stations in vertical component are 1 -2 cm, the locations and difference in vertical displacements between GPS and D-InSAR data are shown in Fig. 9. Fig. 9. Difference map in vertical displacements between GPS data and D-InSAR data associated with Chi-Chi earthquake shown in vertical vectors at 28 sites. The vector face upward shows the vertical displacement from D-InSAR larger than vertical displacements from GPS and the vector face downward shows the vertical displacement from D-InSAR smaller than vertical displacement from GPS. Figure 10 shows the results of vertical displacement for the GPS and D-InSAR data at 28 GPS stations. Except for large change differences at some points, the remaining points were largely consistent. The root mean square (rms) difference between GPS and D-InSAR in vertical displacements was 3.3 cm.

CONCLUSIONS
In this study, D-InSAR was used to detect surface deformation caused by the Chi-Chi earthquake. Deformation on the hanging wall of the Chelungpu fault line could not be detected due to low coherence caused by the mountainous terrain covered with dense vegetation and the large displacement in that area. However, deformation on the footwall of the fault was clearly detectable by the D-InSAR technique; the D-InSAR-based results proved useful in interpreting the deformation caused by the earthquake. While GPS observations can show a three-dimensional deformation field, data are only observed at a limited number of locations; D-InSAR provides better spatial resolution. GPS observations are useful for locations with low coherence and thus not applicable for D-InSAR observations. Using a combination of both GPS and D-InSAR techniques, this study showed that the rms value calculated from the differences in vertical displacements between the methods was 3.3 cm. Thus, both techniques can be used to detect that magnitude of displacement. Both D-InSAR and GPS techniques have high precision. In our study, the results from using both techniques showed good correspondence, suggesting that these two approaches can be used as complementary tools.