Joint Inversion of Geodetic Observations and Relative Weighting—The 1999 M w 7.6 Chi-Chi Earthquake Revisited

: The M w 7.6 Chi-Chi earthquake struck central western Taiwan in 1999. The rupture was complex with several dislocations along the 100-km long Chelungpu thrust fault. Revisiting this earthquake is a challenge, as the precision and coverage of the data sets available are quite poor. Furthermore, the topographic and vegetation coverage complexity of the area prevents coherent radar images. In this paper, radar and optical images, and terrestrial geodetic measurements, were utilised to study the fault. The Particle Swarm Optimization and Okada Inversion (PSOKINV) geodetic inversion package was used with the generalized Akaike’s Bayesian Information Criterion (gABIC) to precisely determine the slip distribution and relative weighting of datasets. Differences in results using the data sets jointly or separately (e.g., under-estimation due to InSAR, inconsistencies in SPOT offsets, smoother slip distribution with gABIC weighting) are observable. Most of the energy was released in the northern part of the fault, where the strike veers toward the east, and mainly at depths less than 4 km. The PSOKINV-gABIC approach is viable for the study of complicated cases such as the Chi-Chi earthquake and can signiﬁcantly beneﬁt the weight determination and physical realism of the fault geometry.


Introduction
Taiwan is the emerged part of an orogen formed by the collision between the Eurasian and Philippine Sea plates [1] entering in convergence at a rate of 82 mm/year with an orientation of 310 • . Two subductions are consuming the island-at the south, the Eurasian plate goes beneath the Philippine sea plate in the east direction at the Manila trench, and the Philippine sea plate subducts towards the north near the Ryukyu trench. It is a high-risk seismic area (Figure 1) with the major and most destructive inland earthquake of the 20th century being the Chi-Chi event in central western Taiwan, occurring on 21st September 1999 with a magnitude of M w 7.6 and at a depth fault surface without introducing overlaps and gaps. It shows similar slip pattern to GPS based models, but the magnitude of the fault slip is smaller and the slip pattern smoother.
The results of Reference [10] based on SPOT imagery pointed out a shallow 20-35 • east dipping ramp with a low dipping decollement at a depth of 6 to 8 km. However, a more recent paper, Reference [16], concluded a decollement at a depth of about 12 km, more consistent with the hypocentral depth of the mainshock and the distribution of aftershocks. According to them, the maximum slip is between 10 and 15 km on the shallow northern part of the fault with a moment of 3.6 × 10 20 N·m but most of the moment was released at shallow depths (between 0 and 5 km). In addition, the area ruptured by the main shock was fully locked before the earthquake. The slip is mainly aseismic for depths above 10 km [16]. No earthquake similar in magnitude to Chi-Chi's is known to have occurred in recent centuries in central western Taiwan. In this study, we use five types (Figure 1) of geodetic observations (GPS, levelling, benchmarks of urban planning, InSAR and SPOT offsets) in a joint inversion in order to estimate the slip distribution with and without precise estimation of the relative weights. In previous models, non-negligible residuals are present due to lack of information (density of the GPS network and no coherence on the hanging-wall for InSAR). This study focuses on showing the importance of combining and using complementary data sets rather than improving the model geometry. The fault geometry was determined by previous studies and slightly modified as a result of our data. The method presented in this paper shows the impact of the different data sets on the slip distribution as well as the necessity to properly determine the relative weights.

Geodetic Data Sets and Processing
The use of geodetic observation data is efficient for the study of earthquakes because of their high precision and their complementarity. In this study, several types have been selected-GPS, ERS-SAR and SPOT-optical images, levelling and benchmarks. GPS data have a high temporal resolution and a precision less than 1 mm/year useful to study slow deformation such as the post-seismic phase but it can also be used to study fast deformation such as the co-seismic phase; however their spatial distribution is sparse. Levelling and benchmarks have also an high spatial precision (few cm or higher) but rely on field measurement campaigns. InSAR data have a high spatial density and sub-cm precision but a low temporal sampling rate depending on the satellite (for instance, 6 days for Sentinel satellites but 35 days for ERS). Finally, SPOT images have a similar spatial density to SAR images, but can often provide useful measurements despite vegetation where SAR images lose coherence and at higher density than a pixel offset map estimated from the SAR amplitude would give. Even though the study area is partially covered by dense vegetation making InSAR output non-coherent, the joint study of InSAR with correlation of optical images and terrestrial measurements is an efficient and precise way to understand earthquakes. In central western Taiwan, from InSAR outputs, information about the footwall of the fault is obtained and SPOT offsets cover the full study area more densely than GPS. Levelling and benchmarks give precision about vertical movements and about the curved part at the north of the fault as benchmarks were measured only in this area ( Figure 1B,C).

GPS Measurements
The GPS network in Taiwan was installed in 1989 and has been surveyed since 1990 covering the whole island. Figure 2 is showing benchmarks and an antenna of the Taiwan network. Specifically for the coseismic study ( Figure 1A), the displacements were extracted from Reference [9]. 41 permanent stations are available (operated by the Central Weather Bureau (CWB); the Satellite Survey Division and Land Survey Bureau, Ministry of Interior (MOI) and the Institute of Earth Sciences, Academia Sinica (IESAS)) and 101 campaign-surveyed stations (conducted by three institutions (the Central Geological Survey, Ministry of Economic Affairs (MOEA, CGS); MOI and IESAS)) covering 1992 to 1999. They separated pre-and post-seismic deformation to estimate the coseismic displacement due to the Chi-Chi earthquake. They were processed using Bernese software v. 4.0 and corrected with precise ephemerides, and residual tropospheric zenith delays were estimated by least squares. The reference station selected is Kinmen, located on the east coast of mainland China (∼257 km from the epicenter) which did not show any coseismic movement. Furthermore, the pre-seismic secular motions were corrected by interpolating the velocity of 24 stations observed from 1992 to 1999. Finally, post-seismic displacements were corrected at the stations surveyed during the three months following the earthquake by interpolation from the nearby stations [17]. GPS, benchmarks and levelling are well correlated with each other. In addition to GPS data, benchmarks of urban planning [18] situated at the north part of the fault, in Fengyuan city, were used (locations shown by the purple dashed box in Figure 1A and the displacements on the footwall and hanging wall in Figure 1B,C respectively). The pre-earthquake benchmarks were measured by total-station and the survey after the earthquake is a mixture of total station trilateration/triangulation and GPS campaigns. The horizontal positions were estimated by total station trilateration/triangulation. The data were provided with respect to the origin of TWD67 System (Transverse Mercator Two-degree Grid system) that we converted to WGS84 (World Geodetic System 1984) to be consistent with the rest of the data sets. The densely spaced distribution of the city-planning benchmarks provides an accurate location of constructions. Before the earthquake, a city map was created giving accurate position of human constructions. As many constructions were damaged during the earthquake, a second survey was run to remap the city after the event, the new map contains all the old city-planning benchmarks as well as Chi-Chi earthquake surface ruptures. The coseismic horizontal displacement data set was extracted from the comparison of these two maps. The horizontal displacements present an uncertainty of 1-2 cm against a mean displacement of 5 m in this area for the coseismic displacements.

Levelling Data
In addition to GPS data, levelling measurements (locations shown by the dark-blue dashed box in Figure 1A and the displacements in Figure 1D) were procured by CGS (MOEA) from 1998 to 2002 and converted from TWD67 to WGS84. They give more information about the vertical displacement, and are consistent with the GPS. Two of the six levelling lines are presented in Figure 3.   Figure 1D for profile locations.

InSAR: ERS-2 Interferometric Processing Chain
ERS images were interferometrically processed using the ESA (European Space Agency) open-software SNAP (Sentinel Application Platform). Their parameters are presented in Table 1 and their coverage shown in Figure 1A. Only the footwall of the fault could be analyzed on the interferograms. The hanging-wall, which experienced the main deformation during this event, is densely vegetated resulting in low coherence. Coseismic interferograms show about 10-11 fringes in the footwall which is equivalent to a surface displacement of up to 30 cm ( Figure 4). Three interferograms were processed using SNAP and SNAPHU (Statistical-cost, Network-flow Algorithm for Phase Unwrapping [19][20][21]) for unwrapping the phase (the interferograms are consistent to each other with a regression coefficient of 0.7-0.8). After the coregistration of the primary and secondary images using Delft orbits and a 30 m SRTM (Shuttle Radar Topography Mission) DEM (Digital Elevation Model), the interferograms were formed and multilooked by 4 × 20 pixels to reduce the speckle effect, which gives a ground resolution of about 80 m. Finally, the interferograms were corrected for the effects of topography using the 1 arc-second (30 m) resolution SRTM DEM and filtered using Goldstein filtering (an adaptative power spectrum filter [22]) before unwrapping the phase and geocoding. Only pixels with a coherence higher than 0.2 were kept to be used in the joint inversion. The low coherence was the main challenge of the InSAR processing especially for the unwrapping of the phase. Furthermore, the hanging-wall of the fault, which is covered by dense vegetation, presents the main deformation as shown by GPS data. Thus, the main area of the deformation is not visible on the interferograms. Due to this reason, we decided to use optical images in addition to the radar interferometry. They have a better spatial resolution (10 m) than pixel offsets from InSAR (about 20 m) and will give us information about the deformation on both sides of the fault.

Sub-Pixel Correlation of SPOT Images
In order to obtain observations of the hanging-wall, the COSI-Corr (Co-registration of Optically Sensed Images and Correlation) software was used to correlate pre-and post-earthquake optical images. SPOT satellite panchromatic images were processed to determine the horizontal coseismic displacement using the subpixel correlation technique ( [23,24]). In this study, we used SPOT 1-2 images with a 10 m resolution and with lower cloud coverage (see Table 2). The locations of the SPOT images are showed by the orange dashed box in Figure 1A. The pre-earthquake image is from 29 January 1999 and 2 months following the event for the second image (23 November 1999). In order to minimize uncertainties in the correction of topographic effects, we selected images covering the area with near-vertical incidence angle (2.9 • for the pre-earthquake image and 2.5 • for the post one). We used a Tandem-X DEM with a spatial resolution of 12.5 m for the processing as the resolution of the images is 10 m (InSAR resolution is 80 m). Even though the images selected are covering 8 months before the event and two months after, no significant pre-seismic deformation was present on the GPS observations and only a few cm of horizontal post-seismic deformation, so we considered that the offsets are essentially representing the co-seismic displacements [10]. In order to obtain the ground deformation, SPOT pairs were processed using COSI-Corr. It is a module developed in IDL (Interactive Data Language) and integrated in ENVI (Environment for Visualizing Images) by the Caltech Tectonics Observatory (USA). The methodology allows an automatic and precise orthorectification and coregistration of satellite images [25]. The main steps of the procedure are an ortho-rectification of the primary image with Tandem-X DEM due to precise selection of ground control points (GCPs), then ortho-rectification of the secondary image using GPS data values in the far field to correct for long wavelength bias. An optimization is applied to improve the ortho-rectification-reduction of the mis-registration unavoidable from the manual tie points selection by optimizing the GCPs of the secondary image (iteration until ground offset between the primary and secondary images is negligible). Finally, the subpixel correlation of the pre and post-event ortho images leads to three bands (East-West ( Figure 5A) and North-South offset fields ( Figure 5B)), and the Signal-to-Noise ratio from which the horizontal ground displacement is retrieved. For the joint inversion, we only kept pixels with a Signal-To-Noise ratio (SNR) above 0.1 and removed the area covered by clouds. An RMS (Root Mean Square) of about 2 m is present between GPS and SPOT offsets ( Figure 6) mainly due to the different number of points and the issues in the far field. Long-wavelength artifacts are observed on the offsets (see Figure 5), mainly for the North-South image, due to satellite attitude changes.  Table 2) acquired before and after the earthquake (scale in m). Grey areas are where the correlation was too low.

Earthquake Modelling
Okada's dislocation theory [26] has been generally used to study and invert seismic fault displacements. However, it only works under the assumption that the Earth crust is a homogeneous elastic half-space and for near-field observations, as the Earth's curvature and radial heterogeneity are not taken into account in Okada's equations. This standard dislocation model calculates an analytical solution for ground deformation due to shear and tensile faults. The fault is modelled by a rectangular geometry (defined by length, width, depth, strike, dip and 3-component dislocation amplitude (rake, slip and opening)) which is decomposed into small dislocation patches. This theory will be used to determine the slip distribution in this study.
From previous models available for the Chi-Chi event in the literature as presented in the introduction of this paper, inconsistencies are observed between rupture processes of the earthquake using one or multiple types of observations. For instance, inconsistencies due to the quality of the data sets: low density of GPS points and low coherence of InSAR, non-negligible residuals between observations and model or the type of dislocation model used introducing overlaps and gaps.
The proposed method in this paper combines five data types in a joint inverse model and thus requires a precise estimation of relative weighting of each data set. In that objective, the generalized Akaike's Bayesian Information Criterion (gABIC) is used.
Furthermore, before estimating the relative weights and inverting the data, in order to make the computation feasible and efficient, InSAR observations are often down-sampled to limited numbers of points (from several thousand points to a few hundred). The Quadtree method [27] was used and displacement gradients are used to determine the sub-sampling window size. Before this, a linear trend is removed from the interferograms (to remove the orbital ramp).

Generalized Akaike's Bayesian Information Criterion (gABIC)
A precise weighting of each data set is fundamental to realize the joint inversion. These weightings strongly affect the result of the inversion (as explained in Reference [28]). For instance, an overweighting of one data set could induce an under-weighting of the other data sets and disturb the estimation of each parameter of the fault geometry. In that goal, the weight determination method developed by Yi ( [29,30]) is used: gABIC based on the entropy maximization principle [31] to determine the regularization parameter (a spatial constraint used to regulate the rupture process), and the Bayesian theorem to estimate the probability distribution of the slip for given observations (GPS, benchmarks, levelling, InSAR, SPOT offsets) and prior constraints (Laplacian constraint in our case). In this generalized version, the type of observations and prior constraints do not need to be specified; only diagonal covariance matrices have to be estimated before running gABIC script.
As we have 8 sets of (pseudo-) observations (3 interferograms, GPS, benchmarks, levelling and SPOT offsets in addition to the Laplacian constraint), 7 hyperparameters need to be determined using the search-based gABIC method. To obtain optimal weight ratios that minimize the ABIC value, large enough ranges of these hyperparameters have to be considered. In practice, searching the seven data sets with a small enough interval is highly time consuming. In order to reduce the processing time, a first estimation of each relative weight was done by calculating it for each data set in respect with the GPS one, which covers the whole area with a low uncertainty. For this first step, the weighting range was determined using a power law method (10 x ) and then refined using a linear method (x) for the interval of inputs (x representing the values to test). Once the interval of hyperparameter was reduced, gABIC could be used using several data sets at once.
The inversion equations [30] are the following (to make the equations more manageable, only one interferogram is mentioned): with: Furthermore, e ∼Q(σ 2 i ) whereQ(σ 2 i ) = σ 2 GPS with Q (i refers to the five data sets) covariance matrices of observations and a similar equation for the constraint). The covariance matrix Q is expressed as follows (blkdiag means block diagonal matrix): The ABIC function was first introduced by Reference [32], using only static GPS observations and Laplacian spatial constraint, as: gABIC is the generalized version, it searches for the hyperparameters (σ 2 i ) that minimize the ABIC value as presented in Reference [29]: where: is the smoothness parameter; • i is a hyperparameter representing the weight of the constraint;

Modelling Strategy: Fault Geometry and Slip Distribution
Regarding the fault geometry (Table 3), this study was based on the previous model designed by References [14,17] and adjusted in accordance with the more recent literature as described in the introduction. Reference [14] is in agreement with the geologically constrained fault model of Reference [33]. Mainly the decollement subfault was adjusted: especially the depth and dip and the rakes of each sub-fault according to the inversion of our data sets. The decollement was fixed arbitrarily to be wide enough so that insignificant slip occurred at one edge of it. Furthermore, for segments 2 and 3, we tried to avoid as many gaps and overlaps as possible. The purpose of the study was to test gABIC on a complicated but known case including several data sets, not to determine the most precise fault geometry. However, for more precision and physical reality a triangular dislocation model would be necessary. Poly3d will be used in a future work [34].
Under the assumption of an elastic homogeneous half-space, the fault geometry and the slip distribution are determined following a 2-step approach using the Particle Swarm Optimization and OKada INVersion package (PSOKINV). It is a geodetic inversion package [35] developed in Matlab which can solve linear and nonlinear problems. Multiple geodetic data sets can be handled by the program.
First, a nonlinear inversion is performed to determine the fault parameters (location, strike, dip, rake, length, width, depth and slip) minimizing the square misfit under the hypothesis of a uniform slip on a rectangular fault using the multi-peak particle swarm optimization (M-PSO) [36], a non linear algorithm approach based on the elastic Okada dislocation models [26]. The slip distribution is then determined as a linear problem, and optimally-smoothed parameters are obtained. This is an iterative technique for estimating optimal dip angle and smoothing factors based on a conjugate gradient least square method. When all the parameters are optimally determined, the fault plane was discretized into rectangular sub-patches of 2 km × 2 km (see Figure 7).  Figure 7. Fault geometry used in this study (Table 3).

Benefit of Multiple Data Sets
The Chi-Chi earthquake is a particularly interesting event to work on a joint inversion methodology as it has been well studied already, its fault geometry is quite stabilised and the island is well monitored so a large variety of data sets is available from satellite images (ERS and SPOT) to field measurements (levelling, benchmarks, GPS in our case).
As explained previously, geodetic observations were selected with their advantages and disadvantages but presenting a good temporal and spatial complementarity. Figure 8 presents the slip distribution obtained for each data set alone to see its impact to the overall slip distribution. For the GPS based slip distribution ( Figure 8A) and adding benchmark and levelling data ( Figure 8B), similar patterns as the one presented in Reference [14] are observed. The maximum of the slip occurs near the surface at the northwestern edge of the fault and negligible slip is present below 10 km depth. The contribution of benchmark and levelling data sets in addition to GPS is useful to get more information in the curved part (where the slip has the highest magnitude and where the geometry of the fault may not be optimal). A shallower slip distribution is visible in comparison to GPS alone (especially subfaults 2 and 3) and we obtain more precision about the vertical movement. From InSAR observations (the result from the interferogram 28 October 1999-15 July 1999 is presented in Figure 8C), the slip is under-estimated in comparison with GPS data due to its asymmetrical distribution, only covering the footwall. This is also why we do not observe any signal on the decollement at depth. Finally, Figure 8D, from only SPOT offsets, is not over-or under-estimating the slip but showing more deformation at depth mainly due to the large uncertainties in the far-field.

gABIC Impact
It is essential in the joint inversion to determine the relative weights between the different data sets and the constraint to properly balance the contributions of each of them. The difference with previous studies is the use of the generalized formula of the ABIC method ( [31]) from Reference [29] (gABIC) to be able to take into account any type of data sets without need of specification and different type of constraints. The point of using gABIC is to obtain a precise value for the weight of each data set that is used for the estimation of the slip distribution. The main issue with this method is that it is computationally intensive. In order to reduce the processing time, the following strategy was followed. First each data set was run separately with the GPS as reference to estimate the relative weight of this data set in relation to the GPS. Then, when a first estimation was made, we could run all data sets together by restricting the search intervals.
The first row of Figure 9 presents the slip distribution obtained from GPS and the three interferograms where first all the weights were set to 1 and secondly the weights were estimated by gABIC. The second row is presenting the slip distribution obtained from all the data sets used jointly following the same steps: first the weights were set to 1, then weights were estimated by gABIC.  (Table 1) with all the weights set to 1; (B) same data sets but with the relative weights estimated by gABIC (GPS 1, 0.09, 0.10, 0.10 for the three interferograms). All the data sets were included in (C,D). After the estimation of the relative weights by gABIC, we can see that the slip distributions are smoother, with a maximum of slip at the curve of the fault of about 12 m. No unexpected slip is present on the edge of the fault. The weights for the benchmarks and levelling is low because the spatial extent of the data sets is small in comparison with the other data sets even though their precision is better.

Seismic Hazard
The seismic moment was estimated based on the inverted slip distribution, with a maximum at 3.9 × 10 20 N·m: µ: Shear modulus taken as 33 GPa, S: Area of the fault, δu: Mean slip from strike slip and dip slip.
The seismic moment release was summed for each 2 × 2 km patch in horizontal rows to obtain the variation in moment release with depth, and summed in vertical columns to get the moment release along strike. The slip distribution and the study of the integrated moment by depth and strike show that the main deformation happened in the north part of the fault, where the fault is curving towards the east with a slip reaching 12.6 m ( Figure 9D) and most of the moment was released at depths shallower than 4 km.
To evaluate the effect of the rupture on the surrounding area, the Coulomb Stress Change (CSC) was calculated using Coulomb 3.3 (Graphic-Rich Deformation and Stress-Change Software) USGS software [37,38]. CSC represents the transfer of stress. It is possible to estimate, from the CSC, the triggering effect of large earthquakes such as the Chi-Chi event. CSC was estimated on a specified fault (strike, dip, and rake specified) at a depth of 11 km with a friction coefficient of 0.4. Figure 10 shows that the increase of Coulomb stress is concentrated at the north of the Chelungpu fault and towards the north-east which indicates that most of the stress was transfered to faults at the north of the Chelungpu fault. As presented by Reference [16], aftershocks (see Figure 11) mainly occured downdip from the rupture area which is consistent with the increase of Coulomb stress that we observe in Figure 10A. We observe in Figure 10B an increase of stress for depth above 10 km which is consistent with the aftershock distribution.

Conclusions
The Chi-Chi earthquake is a complicated event that was caused by the reactivation of the Chelungpu fault whose geometry is quite complex and which occurs in a region of steep topography and dense vegetation.
Five types of geodetic observations (GPS, levelling, benchmarks, InSAR and SPOT offsets) were combined to determine the fault slip distribution of the 1999 Chi-Chi earthquake. A methodology was established in order to precisely estimate the relative weights between the data sets, which is a critical step for the modelling. A wrong estimation will create an under or over-estimation of the slip distribution. For this purpose, the generalized ABIC formula was applied to estimate the relative weights between the observations and the constraint. The slip distribution is mainly determined by the GPS as presented in previous papers, however the incorporation of other geodetic observations constrains the model especially on the curved part of the fault.
The preferred slip distribution including all the data sets and using the relative weights estimated by gABIC shows that the maximum slip was 12.6 m on the north curved part of the fault with a seismic moment of 3.9 × 10 20 N·m, corresponding to a magnitude of M w 7.6. The study of the moment and Coulomb Stress Changes show that most of the slip happened on the north part of the fault at shallow depths and that most the energy released by the event was transferred to the faults at the north of the area, which is now the most probable location of a future event.
In the last decade, satellite images have dramatically improved especially with the Sentinel constellation and the evolution of atmospheric methods, as well as the densification of GPS networks. It is essential to take advantage of all these data, to take into account their quality and coverage by estimating precise relative weights and thus to establish methodology to precisely monitor the earth deformation.  Acknowledgments: Some figures were made using the Generic Mapping Tools by Reference [39]. ERS data were provided by ESA and SPOT images by Airbus Defence and Space.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: