Grid-Based Moment Tensor Inversion Technique by Using 3-D Green ’ s Functions Database : A Demonstration of the 23 October 2004 Taipei Earthquake

Moment tensor inversion is a routine procedure to obtain information on an earthquake source for moment magnitude and focal mechanism. However, the inversion quality is usually controlled by factors such as knowledge of an earthquake location and the suitability of a 1-D velocity model used. Here we present an improved method to invert the moment tensor solution for local earthquakes. The proposed method differs from routine centroid-moment-tensor inversion of the Broadband Array in Taiwan for Seismology in three aspects. First, the inversion is repeated in the neighborhood of an earthquake’s hypocenter on a grid basis. Second, it utilizes Green’s functions based on a true three-dimensional velocity model. And third, it incorporates most of the input waveforms from strong-motion records. The proposed grid-based moment tensor inversion is applied to a local earthquake that occurred near the Taipei basin on 23 October 2004 to demonstrate its effectiveness and superiority over methods used in previous studies. By using the grid-based moment tensor inversion technique and 3-D Green’s functions, the earthquake source parameters, including earthquake location, moment magnitude and focal mechanism, are accurately found that are sufficiently consistent with regional ground motion observations up to a frequency of 1.0 Hz. This approach can obtain more precise source parameters for other earthquakes in or near a well-modeled basin and crustal structure.


InTrODUcTIOn
The moment tensor inversion procedure has long been used to determine earthquake source parameters, including moment magnitude and focal mechanism.But the inversion quality is usually controlled by the prior knowledge of an earthquake location and the suitability of velocity model used for computing Green's functions.These problems become considerable when earthquakes occur in a region with a complex velocity structure such as the Taipei metropolitan area located on a low velocity sedimentary basin and surrounded by mountains and a volcano.
On 23 October 2004 (22:04:27.7 local time), an earthquake (M 3.8) occurred near the city of Taipei (the "1023 Taipei earthquake," Fig. 1).The hypocenter reported by the Central Weather Bureau of Taiwan (CWB) was 121.58°E, 25. 02°N (8.8 km deep).Although the earthquake was small, the hypocenter was very close to the city center (about 10 km).People in the city felt very strong shaking lasting about 5 seconds.The seismic intensities in most parts of the city were about 10 -30 gal on average.This earthquake has been the largest event to occur in the Taipei metropolitan region since 1988.The CWB report posted a first-motion focal mechanism as determined by the Taiwan Strong-Motion Instrumentation Program (TSMIP) (Liu et al. 1999) in the Taipei basin.The result shows a normal fault strike of 275°, dip of 33°, and raking at -105°.Another focal mechanism was reported by the Broadband Array in Taiwan for Seismology centroid-moment-tensor (BATS CMT) (Kao et al. 1998) that indicated a dip-slip fault strike of 3.9°, dip of 84.07°, and raking at 106.5°.Because only two broadband stations were able to well record the event and used in the full waveform inversion for this moderate earthquake, the resolution of BATS CMT was limited to poor station azimuthal coverage.Furthermore, the determinations of these two sets of focal mechanism were both based on the earthquake location reported by CWB.This hypocenter was determined from a 1-D Taiwan average crust model (Chen and Shin 1998) which cannot provide sufficient resolution for the study of local earthquakes.
Recent studies of this earthquake, including Lin (2005) which indicated this earthquake was possibly triggered by the construction of the world's then tallest building, the Taipei 101 building.Chen et al. (2010) relocated this earthquake by using the JHD method (Douglas 1967) and found that the hypocenter moved toward the south and the earthquake was related to a blind normal fault beneath the basin.
Here we perform a forward simulation by using the CWB first-motion solution based upon a 3-D spectral-element method (SEM) (Komatitsch and Vilotte 1998;Komatitsch et al. 2004) simulation to examine the suitability of source parameters.As shown in Fig. 2, most of the synthetics cannot fit the observations in both timing and polarization.Except the uncertainty from 3-D velocity model, the shift in time might due to an error in source location and the mismatch of phase polarization is likely to have come from the focal mechanism.
In order to evaluate the safety of Taipei city, it is necessary to investigate the source parameters of the 1023 Taipei earthquake carefully, including the hypocentral location, focal mechanism and its magnitude.In this study, we perform a grid-based moment tensor inversion technique using the SEM for 3-D Green's functions.The scheme is comparable to Grid MT (Kawakatsu 1998;Fukuyama and Dreger 2000;Ito et al. 2006;Tsuruoka et al. 2009) in which real-time monitoring determines the moment tensor in parts of the Japan area.The drawbacks of the existing Grid MT methods are the use of longer period waveform and a simplified velocity model used to calculation Green's functions.When considering near field strong motion data, these constraints will significantly influence the moment tensor inversion result.Here we apply a different inversion algorithm and consider 3-D Green's functions at a high frequency (1.1 Hz) to accurately locate the earthquake source and obtain focal mechanisms at one time.The 1023 Taipei earthquake which occurred near the city of Taipei is undoubtedly a good candidate to examine this new moment tensor inversion procedure.

DATA AnD METhOD
Records from the TSMIP are used in this study to provide full waveform information in the inversion.The stations tracking the 1023 Taipei earthquake are shown in Fig. 1.The total number of triggered TSMIP stations was 59, most of them located north of the epicenter, in the middle part of the basin and along the basin edge.Records from stations located in the western part of the basin have a poorer quality with high background noise signals.Here, we focus on stations with a better record quality which are mostly located in the eastern part of the basin and close to the epicenter (as indicated by the study area in Fig. 1).Another data set, the downhole array operated by the Institute of Earth Sciences, Academia Sinica (IES), is also considered in the inversion analysis (Fig. 1).Four sites compiled 14 records of downhole array data triggered by the Taipei earthquake.This data set provides good data quality and absolute timing which can help to control the moment tensor inversion result.All the used data are integral to velocity and band-pass filtered between 0.1 and 1.1 Hz.

Grid-Based Moment Tensor Inversion
Here, a grid-based moment tensor inversion technique is utilized to accurately locate an earthquake hypocenter and obtain source focal mechanisms simultaneously.The technique combines a full-waveform CMT inversion and grid search procedure which determines the best-fit solution found within the Green's functions database.The inversion flowchart is presented in Fig. 3. Using the 1023 Taipei earthquake as an example (Fig. 1), we define a 3-D grid around the original hypocenter first, and then calculate the 3-D Green's functions for six moment tensor dipoles with a unit moment of 1.0 × 10 25 dyne-cm at each grid point.The small open circles show the location of Green's functions grids.This 3-D grid covers about 6 × 6 × 10 km 3 , each side has 10 grid points and thus a total of 1000 grids are set.The inversion is then based on this grid distributed 3-D Green's functions database, searching for the best-fit with a minimum misfit between synthetic and observed waveforms at all grid points.The result provides a location with 600 m horizontal and 1000 m vertical resolution.The characteristics of 3-D Green's functions and the inversion method are discussed in the following sections.

3-D Green's Functions
The effects of lateral velocity variation can influence the determination of a source location especially in a region with a complex tectonic background like the Taipei basin.The source location will further affect the determination of the focal mechanism.In order to account for the possible 3-D effect in the Taipei basin (Lee et al. 2008a), we consider Green's functions in full 3-D scale.The Green's functions database is calculated using the spectral-element method (SEM), and a single-couple point source with a half-duration of 0.39 sec is considered on each grid.Figure 4 shows 1-D layered and 3-D SEM mesh models as well as examples of their synthetic waveforms.By the characteristics of SEM, the surface topography (40 m DEM), largescale 3-D velocity (Kim et al. 2005), basin geometry and low velocity materials in the basin (Wang et al. 2004) are carefully incorporated in the 3-D mesh model.The slowest compressional-and shear-wave speeds in the Taipei basin model are 1.5 and 0.2 km s -1 , respectively.More details about the construction of Taipei basin SEM mesh can be found in Lee et al. (2008b).For a comparison of synthetic waveforms, 3-D synthetics have more complex waveforms and larger amplitudes compared to 1-D results, especially for the S-wave and later phases.The travel-time difference can also be found among two sets of synthetic waveforms, even if the hypocentral distance is less than 20 km.All of these discrepancies between 1-D and 3-D Green's functions will affect the inversion result which can help to clarify the issue of how source parameters are affected by the velocity model used.

Inversion Method
We use a non-negative least square (NNLS) inversion technique to invert the densities of six moment tensor dipoles.The NNLS inversion is usually used in finite fault source inversion studies which has been used successfully to investigate the slip distribution of many large earthquakes worldwide (such as Hartzell 1983;Wald and Heaton 1994;Lee and Ma 2000;Ma et al. 2001).Here we reformulate  this inversion algorithm and apply it to the moment tensor inversion.The inversion can be formed as a linear problem and written as Ax = b, where A is the matrix of Green's functions, x is the densities of six moment tensor dipoles that we want to invert and b is the observed data vector (Fig. 5a).This system of equations can be solved by the standard least-squares method by adding the constraints for moment tensor, including isotropic and minimum solution.However, due to the non-negative solution in NNLS, the negative moment tensor density can not be solved for the x vector.To deal with this problem, we expand matrix A from 6 to 6 × 2 rows, the extended elements are mapped from the original A matrix and multiplying by (-1) as negative Green's functions.At the same time, the x vector also expands to 6 × 2 elements.In this approach, both the positive (x 1 to x 6 ) and negative (x 7 to x 12 ) moment tensor density can be considered in a NNLS inversion.We use a misfit function, defined as Ax b b 2 2 -^h , to evaluate the quality of a solution.
Since matrix A is symmetrical from positive and negative scales of Green's functions, we can apply a parallel NNLS inversion technique to improve inversion performance (Lee et al. 2006(Lee et al. , 2008c)).We put all positive Green's functions in one computing node and all negative Green's functions in another one, then use the Message Passing Interface (MPI) (Gropp et al. 1996) as the communicant between the two MPI processes during the inversion (Fig. 5b).By applying the parallel NNLS technique, the computing time of each inversion can be greatly reduced.It is worth noting that the relatively small amount of time required on each inversion will shorten the overall grid search timing.The trend of multiple processors can be straightforwardly applied to the parallel NNLS on a single multiple-core workstation or desktop computers.

Benchmark Tests
Two benchmark tests are performed to examine the grid-based moment tensor inversion procedure.The first benchmark assumes a vertical strike-slip source (M 12 ) occurred at grid point (in this case, 121.575°E, 25.043°N, depth 9.75 km).Three components in Green's functions at station TAP020 (with a moment of 1.0 × 10 20 dyne-cm) are multiplied by 10 5 and then set as input records.Three different random noise levels, 0%, 20%, and 40% of the maximum waveform amplitudes, are added on the input data to test the stability of this inversion procedure.The inversion result shows that not only the inverted synthetic waveforms can have a good fit with the input records, but also the moment, hypocenter location and focal mechanism are well determined (Fig. 6a).The test shows that the grid-based moment tensor inversion is a useful procedure which can help to determine the source parameters utilizing limited data.
The second benchmark tests the sensitivity between the timing of synthetic waveform and source parameters in both location and focal mechanism.The north component velocity seismogram recorded by TAP026 during the 1023 Taipei earthquake is considered as the input waveform.We shift it in time by +1.0, 0.0, -1.0, and -2.0 sec with respect to the original record.The inversion results are shown in Fig. 6b.Most of the synthetic waveforms can fit well with the shifted records.When the record has a negative shift in time, the hypocenter location will move toward station TAP026.This conclusion is reasonable in that with less travel-time, a shorter hypocentral distance is needed.The focal mechanisms are also varied with the shifted waveforms.This test indicates that the inverted source parameters are sensitive to the input waveform, especially in the timing.In other words, we can use the characteristics of local records in both timing and amplitude of waveforms to locate the earthquake and obtain its focal mechanism.

ExAMplE OF 23 OcTOBEr 2004 TAIpEI EAr-ThqUAkE M 3.8
Once the stability and accuracy of grid-based moment tensor inversion has been thoroughly validated, we aim now to investigate source parameters of a real earthquake, i.e., the 23 October 2004 Taipei earthquake (M 3.8).For the waveforms used in the inversion, we adopt most of the downhole data to control the inversion quality, several TSMIP records are also considered to provide better station azimuthal coverage.To ensure good data quality in the inversion, the north and east component velocity waveforms are used, but vertical component which usually has more complex waveform is temporarily held in abeyance.We consider an 11.5-second time window in the inversion which starts from the original time of the 1023 earthquake.Both Green's functions and observed data are band-pass filtered between 0.1 and 1.1 Hz.

source parameters from Inversion result
After a large number of inversions during the grid search, the best-fit solutions for source location, seismic moment, and its focal mechanism are found.The best fault plane solution and better solutions among different depths at every position are shown in Fig. 7.Note that the focal mechanisms become more similar when their locations approach to the best-fit solution.The inversion result indicates that the location of the Taipei earthquake should be moved to the left, about 3 km from that given by the CWB report.The depth becomes shallower which moves from 8.8 to 7.75 km.The new epicenter is located at the basin edge which is very close to the Taipei 101 building.A normal fault strike of 71.8°, dip of 59.6° and rake of -34.9° is indi- cated from the inversion result.The moment magnitude is M 3.7 which is among the reports from CWB and BATS.The synthetic waveforms determined from the best-fit solution can fit well with the observed data in both timing and amplitude of the major phases.

Forward Examinations
To examine the inverted source parameters, we perform a full 3-D wave propagation simulation based on SEM.Comparisons between observed velocity records from the TSMIP and the corresponding synthetic seismograms are shown in Fig. 8.All the waveforms are band-pass filtered between 0.1 and 1.0 Hz.The synthetics for the horizontal (North and East) components explain both the major phases and the later arrivals quite well.On the smaller-amplitude vertical component, however, there are larger discrepancies.In general, the synthetics can explain most of the ground motion characteristics recorded in the seismograms.
We further analyze the maximum particle acceleration to obtain the information on peak ground acceleration (PGA) distribution (Fig. 9a).The observed PGA values are obtained from the norm of all three components (low-pass filtered with a corner frequency of 1.0 Hz) to compare with the corresponding simulated results.The simulated PGA pattern can explain the characteristics of observed PGA distribution well.Significant amplification is observed in the basin with relatively large PGA values compared to its surroundings.Because the epicenter is located close to the basin, source radiation dominates the PGA distribution in the near-field records.Combined source radiation and basin amplification effects produce large PGA values in the southeastern part of the basin.Outside the basin the PGA values exhibit a symmetric pattern centered on the epicenter with a strike determined by the focal mechanism.The simulated PGA distribution determined from the CWB fault plane solution is also shown in Fig. 9b for comparison.The results from the CWB report cannot explain the characteristics of observed PGA values, such as northern basin exhibits relatively small observed PGA and south eastern basin should have a stronger PGA compared to other places in the basin.The forward examinations indicate that source parameters determined from our study can better explain the observed data in both waveforms and PGA distribution.

DIscUssIOn
We placed a great deal of emphasis in obtaining precise source parameters in this study using 3-D Green's functions.From the grid-based moment tensor inversion result of the 1023 Taipei earthquake, the source location is a bit changed compared to the CWB report.The earthquake source location reported by the CWB is determined from short-period records using 1-D velocity model to calculate the travel-time.In general, as the seismic velocity has been underestimated in the model, wherein, the determined earthquake source location will move toward low velocity area; in the opposite direction, the source has been moved to a high velocity area when velocity was overestimated in the model.The error in source location will lead to a further misinterpretation of its focal mechanism.This problem is considerable because the earthquake occurred in a region with a complex velocity structure.

Influence of Velocity Model
To investigate the influence of a 3-D velocity structure on the determination of source parameters, we consider both 3-D and 1-D velocity models as shown in Fig. 4  -^h , for the 1-D model is 0.717 which is larger than that obtained from the 3-D model (misfit 0.658).The discrepancies between 1-D and 3-D inversion results, especially in the hypocenter location, are mainly dependant on whether Green's functions can reproduce the correct travel-time information.In this case, the 1-D model can not represent the low velocity sedimentary in the basin, and thus this causes the source to move to the mountain area with a higher seismic velocity because of the Tertiary outcrop.

comparisons of Focal Mechanisms
Figure 11 shows the locations and focal mechanisms determined by different agencies and studies.The detailed source parameters are presented in Table 1.For the result of focal mechanisms, the JHD (Chen et al. 2010) and gird-based moment tensor inversion (both in 1-D and 3-D) are similar.For the hypocenter location, the discrepancies among most results are not more than ±5 km, except for the BATS which adopts the epicenter location from CWB and then finds the best focal mechanism at a depth of 17 km.The location error within ±5 km is acceptable in an ordinary earthquake report.But sometimes this error is crucial for the determination of the focal mechanism especially when using near field records.By the use of a grid-based moment tensor inversion and 3-D Green's functions, we improve source parameter accuracy in both location and focal mechanism.This result points out that adopting a more realistic   velocity model to calculate Green's functions is important which plays an important role in the determination of source parameters.

cOnclUsIOn
We have developed a grid-based moment tensor inversion technique which is based on a parallel NNLS inversion, grid-search procedure and 3-D SEM Green's functions da-tabase.By using this new technique, the source parameters, including hypocentral location, seismic moment, as well as focal mechanism can be precisely obtained all at one time.A serial of benchmark tests have shown that this technique is useful for application to moderate, regional earthquakes.The inversion scheme is applied to the 23 October 2004 Taipei earthquake where a reasonably set of strong-motion records are obtained.The inverted source location and focal mechanism are in good agreement with high frequency (1.0 Hz of velocity waveform) strong-motion records, including TSMIP and IES downhole array data.Inversion result shows that the hypocenter has moved by a few km from the CWB, BATS and the JHD locations to the eastern boundary of the Taipei basin where the Taipei 101 building is located.This hypocenter move is also associated with appreciated changes in moment tensor solution.

Fig. 1 .
Fig. 1.A grid-based moment tensor inversion as an example of the 23 October 2004 Taipei earthquake (M 3.8).The large rectangle shows the study area and small rectangle indicates the area of Green's functions grid.The depth distribution (in meters) of the Taipei basin basement is shown by the dashed lines.Open circles present the location of Green's functions grid.The Taipei earthquake is indicated by a star and its focal mechanism is shown by a beach ball.The stations triggered by this earthquake including the TSMIP stations (solid square) and IES downhole stations (open triangle) are also shown.

Fig. 2 .
Fig. 2. Forward simulation results using source parameters from the CWB report.Black lines are observations recorded by the TSMIP and gray lines are synthetic waveforms calculated by SEM.The beach ball shows the location and focal mechanism of the 1023 Taipei earthquake.All the waveforms are band-pass filtered between 0.1 and 1.1 Hz.

Fig. 4 .
Fig. 4. Comparison between a 1-D and 3-D Green's functions.The left diagram shows the 1-D layered and 3-D velocity model used for the calculation.The surface topography (40 m DEM), large-scale 3-D velocity (Kim et al. 2005), basin geometry and low velocity material in the basin (Wang et al. 2004) are carefully considered in the 3-D SEM mesh model.The right diagram shows the waveforms determined by 1-D layered (black lines) and 3-D (red lines) velocity models.

Fig. 5 .
Fig. 5.The Ax = b linear equation system.(a) Matrix system for moment tensor inversion where A is the matrix of Green's functions, b is the observed data vector, and is the solution vector of the density of moment tensor dipoles.(b) Application of Parallel NNLS on moment tensor inversion.The Message Passing Interface is applied as the communicant between two computing nodes in the parallel computing process.

Fig. 6 .
Fig. 6.Benchmark tests of grid-based moment tensor inversion.(a) A stability test using the Green's functions at grid no.587 for station TAP020 as the input data.Three different random noise levels, 0%, 20%, and 40% of the maximum waveform amplitudes, are added on the input data to test the stability of inversion.(b) A sensitivity test shifting the time in the observed waveform at station TAP026.The input waveforms are presented by black lines and synthetic waveforms are shown by gray dotted lines.All the waveforms are band-pass filtered between 0.1 and 1.1 Hz.Inverted focal mechanisms and their locations are shown by beach balls.Tables in (a) and (b) show the inversion results of source parameters and misfits.

Fig. 7 .
Fig. 7.A grid-based moment tensor inversion result of the 1023 Taipei earthquake.Beach balls at different locations show the better CMT solutions among different depths.The normalized misfit, as defined in the right-hand legend, is used to show the inversion result.The larger beach ball indicates the smaller misfit between synthetic and observation.The hollow star is the original earthquake location reported by the CWB.The best-fit solution is marked by an open square.Black lines represent observed velocity data used in the inversion and colored lines are synthetic waveforms determined by best-fit solution.All the waveforms are band-pass filtered between 0.1 and 1.1 Hz.

Fig. 9 .
Fig. 9. Comparison between simulated and observed peak ground acceleration (PGA) using source parameters from (a) grid-based moment tensor inversion results and (b) the CMB report.The observed PGA values are obtained from the norm of all three components (low-pass filtered with a corner frequency of 1.0 Hz) to compare with the corresponding simulated results.The observed PGA is shown by the open circles, and the simulated results are drawn using the rainbow scale shown at the right.The epicenter and source mechanism of the Taipei earthquake are represented by the "beach ball." to calculate Green's functions and then invert the source parameters of the Taipei earthquake.The grid-based moment tensor inversion results are shown in Fig. 10.The focal mechanisms determined from the two models are similar and both show normal faulting fault plane solutions.But their locations are disparate, the result from 1-D moves to the south about 3 km distant from the basin.Chen et al. (2010) used the JHD method to relocate the earthquake that also obtained a similar result.The waveform misfit, defined as Ax b b 2 2

Fig. 10 .
Fig. 10.Comparison of the grid-based moment tensor inversion result using 1-D and 3-D Green's functions.(a) Results determined by 3-D Green's functions; (b) results obtained from 1-D Green's functions.The caption of grid-based moment tensor inversion result is the same as that in Fig. 7.

Fig. 11 .
Fig. 11.The locations and focal mechanisms determined by different agencies and studies.Their source parameters are listed in Table1.