Estimation of Source Parameters by the Inversion of Near Source Strong Motion Wave Forms

A method has been developed allowing the direct inversion of near source · strong motion records to obtain an estimate of the seismic source from . a moment tensor. The inversion is implemented in the time domain, and the least-square method is used to estimate the source mechanism: Significant verification using synthetic seismograms including artificial noise was done. The inversion of synthesized data suggests that with this method it may be possible to retrieve the source information of local earthquakes using only a few high-quality, three-component seismograms. This method has been ap­ plied to the study of the source parameters of the ML= 5.4 event of May 28, 1992 in southeastern Taiwan using only the strong motion data recorded at Chengkung. The retrieved-source mechanism is in agreement with the distri­ bution of the few observed first P motion of short-period network data and is similar to the CMT solution. The advantage of this method is its capability for fast determination of the source mechanisms for locally felt events using strong motion data from the currently distributed seismic stations.


INTRODUCTION
The earthquake source mechanism has been traditionally determined by using P-wave fi rst-motion data and teleseismic wave forms. A fault plane solution was obtained graphically from P-wave first-motion data by trial and error. Such a graphical method, however, is apt to lack objectivity. The basic difficulty of this method is that the distribution of stations from which the P-wave first-motion data are obtained is often inadequate to constrain two nodal planes. For offshore earthquakes near eastern Taiwan, the coverage of the P-wave first motion data from in land stations is not sufficient to determine the source mechanism of the events, and the wave forms which are excited by most of the events with magnitude less than 5.5 are too weak to clearly be recorded by seismic stations on teleseismic ranges. On the other hand, the strong motion instruments near the epicenter are always triggered. These data include not only arrival times but also wave forms which · are useful in determining source parameters.
The wave-form inversion for source mechanism was initiated by Gilbert (1970;1973), who introduced moment tensors for calculating the displacement at the free surface; this can be expressed as the sum of moment-tensor elements multiplied by the corresponding Green's function. Thus, a linear relation exists between the moment tensors and Green's function elements. Gilbert (1973) suggested that this property should be utilized in formulating the inversion problem for source mechanism. The concept of seismic-moment tensors was further extended by Backus and Mulcahy (1976), and Backus (1977). The reason that moment tensors are important is that they completely describe, in a first-order approximation, the equivalent forces of general seismic point sources. The equivalent forces can be correlated to physical source models such as shear dislocation and explosion (Aki and Richards, 1980). Moment tensors can be determined from the free oscillation of the earth (Gilbert and Dziewonski;1975), long-period surface waves (Mccowan, 1976;Kanamori and Given, 1981), long-period body waves (Stump and Johnson, 1977) or short-period body waves (Jost and Herrmann, 1989). Dziewonski et al. (1981) used the Centroid-Moment Te nsor (CMT) inversion method to simultaneously obtain the moment tensor and the hypercentral parameters of the best-point sources, called the centroid (Backus and Mulcahy, 1976), by using teleseismic long-period seismic records. Presently, earthquakes with M0 greater than 5x1023 dyne-cm have been routinely processed by the Harvard group. These solutions are published in the Preliminary Determination of Epicenters (PDE) listings. However, for earthquakes with Mo less than the low limits, the resolution of teleseismic records are limited to determine earthquake source mechanisms. Some wave form inversion methods were proposed to determine the source mechanisms from regional seismic records (Dreger and Helmberger, 1991 a, b;Thio and Kanamori, 1991). The aim of this study is to develop a moment inversion methOd to determine the source mechanisms of local earthquakes from strong motion data. The method proposed by Jost and Hemnann (1989) to determine the moment tensors based on the fitness of the short-period SH-wave peak amplitudes was applied to match the three-component wave forms of high quality near-source records. The major task of this study was to describe the formulation and testing of the developed programs to obtain the source mechanism from near source strong motion data recorded by a single or small number of well-calibrated stations. Results of numerical experiments and real-data inversion are presented.

Representation for the Moment Tensor
By using the representation theorem for seismic sources (Aki and Richards, 1980), the observed displacements d( x, t) at arbitrary position x at time t due to a distribution of equivalent body force, fk, in a source region is where Gnk are the components of Green's function containing the propagation effects, and V is the source volume where f k are non-zero. The subscript n indicates the component of dis � lacement.
The comma between indices in (2) describes partial derivatives with respect to the coordinates after the comma. The moment of the equivale1_1t body forces is defined as the force moment tensor: If the conservation of linear momentum for the body forces is assumed, then the term M k does not exist in (3), and only the term, m=l, in equation (3) chosen (Jost and Herrmann, 1989) for equation (1), then becomes : ( 4 ) where * denotes the temporal convolution. Based on the assumption that all the compo nents of the time dependent seismic moment tensor in (4) have the same time function s(t)' (synchronous source, Silver and Jordan, 1982) equation (4) becomes : where Mkj is now a set of constants and represents the components of the second order seismic moment tensor M. The important point of the representation of equation (5) is that displacement d n can be expressed as a linear function of the moment tensor elements and Green's functions.

Linear Moment Te nsor Inversion
Equation (5) can be expressed as a matrix form as in : where the vector d consists of n sampled values of the observed ground displacement at various arrival times, stations, azimuths and components. G is an nx6 matrix containing Green's functions calculated using an appropriate algorithm, and an earth model, for example the half-space Green's functions of Johnson (1974), while m is a vector containing the 6 moment tensor elements; thus mxx• myy• mzz,, mxy• ffixz• and myz• are to be determined (Stump and Johnson, 1977). Usually, the observed data points are appended as a vector from the first to last digital points of a single seismogram, then the other two components are added by the same procedure, then the data of other stations are added sequentially. The Green's functions, in this study, convolved the source time function first and reduced it to six components; Gxx, Gyy• Gzz,, Gxy• Gxz• and Gyz· Each component has n digital points corresPQnding to the element of d. The researcher defines a generalized prediction error as E = eT We e, where the vector e represents the prediction error of each observation and the matrix We stands for the relative contribution of each individual error to the total prediction error. The least-square solution of (6) is obtained by solving the normal equations : The estimated model parameters can be expressed as : The inversion procedure is very close to that of Jost and Herrmann (1989), the difference being in the introduction of the 3-component wave-form fi tting. Jost and Herrmann (1989) used the peak amplitudes of the SH waves only on their inversion code. To determine the elements of moment tensor, observations at more than 5 distinct azimuths are necessary with the Jost and Herrmann method. In this study, for the three-component full wave form information which included the amplitudes, travel times and polarization directions used, the number of necessary stations was reduced. In fact, the moment tensor can be determined by only using data from one-station, 3-component seismograms.

Numerical Implementation
In this study, a numerical program based on equation (5) was implemented. The half space Green's function developed by Johnson (1974) was employed in this inversion. This Green's function is based on the Laplace transform technique and leads to the Cagniard-de Hoop solution. To generate synthetic seismograms, the following source-time function is used (Herrmann, 1979) : The definition of the moment tensor elements followed equations (A5.4), (A5.5) and (A5.6) in Jost and Herrmann (1989). The fault-plane parameters defined in this study are the same as those of Herrmann and Wang (1985) shown in Figure 1. To check the definition of source parameters and the accuracy of moment-tensor Green's functions, the seismograms in Figure  9 of Herrmann and Wang (1985) were reproduced as shown in Figure 2.
The numerical implementation of linear inversion was done using the code by Press et al. (1987). The singular value decomposition (SYD) method was chosen to solve the linear least-square problem of equation (5). The FORTRAN computer code used to form the inversion algorithm is the program SVDAT (Press et al., 1987) to fit the observed data which uses the program SYDCMP to compute the SYD and the program SYDKSB to evaluate the chi-square test. Fig. 1. Definition of the Cartesian coordinates (x,y,z). The origin is at the epicen ter. Strike (t/>8) is measured clockwise from north, dip (6) from horizontal down, and slip (,\) counterclockwise from horizontal. u and v, are the slip vector and fault normal, respectively. ¢ is the azimuth measured clockwise from north (from Aki and Richards, 1980).   (Johnson, 1974), which were used in this study.
In this study, all parameters which were used to compute the synthetic seismograms are following Herrmann and Wang (1985).

Test on Synthetic Data
To prove the ability of the moment tensor inversion method proposed in this study to retrieve a fault plane solution by matching the observed 3-component nearby strong motion wave form, a series of synthetic tests was performed.
The first test was proposed to examine the resolution of inversion based on ground displacement and velocity. The source used for generating the 'observed' records is given by dip=45°, slip=:-67.5°, and strike=0° and situated at a depth of 10 km. The half-space earth model with V p=6.15 km/sec, V ,,=3.55 km/sec, and density=2.8 gm/cm3, which was used by Herrmann and Wang (1985), is used, and remains the same for all tests in this study. The source time function is set as the form of equation (9) with r= 0.5 sec, and is the same for all tests in this study. Both 3-component single-station displacement and velocity of 'observed' seismograms at the epicentral distance of 55 km, with azimuth of 0° were inverted, as shown in Figure 3. The fit of synthetic seismograms is quite impressive for both sets of 'observed data'. However, the source parameters inverted based on ground displacement are more accurate than those based on velocity, as shown in Figure 3.
To investigate the resolution of the observed data from different azimuths to inverse source mechanism , a pure strike-slip fault with slip=0°, and strike=0° (its auxiliary plane with strike=90°, and slip=180°) has been considered. Two sets of 'observed' seismograms from azimuths of 230° and 180° were generated. The source mechanism was retrieved from the wave form inversion of single station 3-component displacements. The fit of synthetic seismograms and the inverted source mechanisms are shown in Figure 4. The source mech anism which was retrieved from the 'observed' data of azimuth 230° is different from than that of azimuth 180°. The data from azimuth 180° exactly reproduced the strike and dip angle of the fault plane; however, the resolution of the auxiliary plane of the observation point located at the node of P and SV wave is poor (Figure 4 b), with a 23° deviation in the dip angle of the auxiliary plane. Although some non double-couple term exists in the inverted source mechanism of the 'observed' data of azimuth 230°, the source mechanism retrieved from the best-fit double couple solution is very consistent with the original fault.
The next inversion was performed as a test of the retrieved accuracy of source mech anism from multiple station data. The pure strike slip fault proposed above was followed to generate 'observed' data. Of the two sets of seismograms, one included 3 observations from azimuth 40°, 130° and 230°, whereas the other one included 6 observations from 0°, 85°, 185°, 230°, 275°, and 320°. The retrieved source mechanisms from both data sets are shown in Figure 5. The source mechanism inverted from multiple stations ( Figure 5) is better than that from a single station. However, in Figure 5, the source mechanism retrieved from the 3 stations is more accurate than that from the 6 stations. To learn the reason for that, the observations were checked and 4 observed points of Figure 5 (a) were found close to the nodes of P and SV waves . As shown in Figure 4, since the resolution for the source mechanism retrieved from those data are poor, the azimuth distribution of the observed data affects the resolution of source inversion.
To study the influence of random noise, 'observed' data contaminated with noise were examined. In this test, the pure strike-slip fault as used in Figure 4 was followed. The seismogram of an epicentral distance of 75 km and on an azimuth of 230° was computed as noise-free observed data. The data were added with random noise of amplitude varying from 10 to 30 per cent of the peak-signal amplitude. The inverted source mechanisms from different noise level data are shown in Figure 6. The source mechanisms retrieved from all      noise data were still consistent with the true solution. Figure 7 shows the fit of synthetic seismograms to noise data with amplitude at 30 per cent of the peak signal level; herein although the noise wave forms were.heavily distorted, the inverted source mechanism was still consistent with the true solution. Hence, the inversion code is very stable for the purpose of source mechanism retrieve.

Inversion of Observed Wave Forms
The method presented in the above section has been applied to the inversion of strong motion seismograms recorded at Chengkung (CHK), a permanent station of the Central Weather Bureau (CWB) of eastern Taiwan, during the earthquake of 1992 May 28 (Origin time: 5/28/1992 23:19:35.6 GMT; ML=5.4; 23.15N, 121.35E; depth=13.7). These records are interesting because the station was only about 5 km from the epicenter, not only the far-field but also the near-fi eld displacements were clearly recorded. Figure 8 shows the acceleration records by the Geotech A800 Force Balance Accelerometer at CHK station. The maximum acceleration at this station was 215 gals. To apply our inversion technique to the observed waveforms, the N-S and E-W components were rotated into radial and trans verse components. The displacement trace was then obtained from the acceleration record by time-domain integration with high pass filtering at 0.4 Hz to remove the baseline drift, as shown in Figure 9. The far-fi eld pulses and the near-field displacement (displacement w w TAO, Vol.5, No.1, March 1994  Random noise is added with an amplitude of 10%, 20%, and 30% of the maximum signal amplitude of (a).
between P and S) are clearly recorded. Since the distance to the hypocenter is very short (�= 4.76 km, azimuth=198.4°, backazimuth=18.4°), the arriving phases can be assumed to be directly propagated from source and the near-field displacements should be considered. However, the half-space full wave form Green's function (Johnson, 1974) is suitable for wave-form synthesis. For simplicity, the earth model used to compute Green's functions of this event is assumed to be a homogeneous half-space (P velocity=6.15 km/sec, S veloc ity=3.3 km/sec, density=2.8 g/cm3). For the digital wave-form fit required by the inversion scheme, the observed data were resampled to the same time interval of Green's functions (10 samples/second, in this study). To reduce the high-frequency scattering which is not included in synthetic wave forms, the data experienced band-pass filtering between the 0.4 and 2.0 Hz frequency bands, and the same filtering was done on the theoretical Green's functions for consistent frequency content on wave-form fit. Figure 10 (a) shows the retrieved source mechanism obtained by the wave-form inversion using only the 3-component displacements of the CHK station. The best-fit double-couple mechanism had the following parameters: for plane 1, strike = 281° dip = 72° and slip = 170°; for plane 2, strike = 14 °, dip = 81 ° and Synthetic observed data which included random noise of 30% of peak sigrial amplitudes (upper) and synthetic seismograms (lower) reproduced from source parameters after moment tensor inversion (Figure 6(d)).

21
slip=l8°. The limited first P motion data identified from the Central Weather Bureau short period network is also shown. Figure 11 compares the observations and synthetic seismo grams computed by the source mechanism from the wave form inversion. Although there are some later phases after the S wave with small amplitudes on the observed data which cannot be fit using the half-space Green's functions, the overall observed and synthetic seismograms agree for each component. Figure 10 (b) represents the source mechanism of this event obtained from the PDE monthly listing, which is determined by the long�period body-wave inversion method (Sipkin, 1982) or the CMT inversion method (Dziewonski et al., 1981) using data from global networks. Both methods got similar results. It is well known that the source mechanism retrieved from regional network data (high frequencies) is need not to have exactly the same solution as from the global network data (long-period) for the different frequency bands of wave forms used in both methods. Therefore, the difference between the solution of this study and the CMT solution is reasonable.

DISCUSSIONS
The determination of the source parameters of moderate sized (ML <5.5) earthquakes of Taiwan and its nearby area is an important seismological problem. Earthquakes of moder ate size are widespread in Ta iwan and in the nearby area, and these provide major clues to the   , Vol.5, No.1, March 1994 active tectonics of Taiwan; however, such events are too small to be well recorded teleseis mically. Traditionally, these source mechanisms were studied from the fi rst P motion data · of a local short-period network. Unfortunately, for earthquakes occurring offshore of eastern Taiwan, the first P motion on the focal sphere is quite ill constrained. The focal mechanisms of these events are difficult to obtain using only the sparsely distributed P motion data. The recommendation of the inversion procedure developed in this study for source parameters in the Taiwan area is founded on several reason�. For the seismograms recorded from local to regional distance range, the structure of the crust has a strong effect on the body wave forms. To estimate the source parameters accurately, the complete Green's functions from the detailed structure are necessary. However, on the complex, active tectonic region of Tai wan (Ho, 1982), the crust structure for regional wave-form modeling is still not available. In other words, the wave form inversion methods (Wallace and Helmberger, 1982;Patton, 1988) which were suggested for regional body wave and surface wave are too underdeveloped to apply to the Taiwan area. The other reason is that many high dynamic range strong motion instruments have been deployed on the major seismically active areas of Taiwan, and there are many new instruments which will be installed in the near future by the CWB. It is easy to find several strong motion records for most moderate sized earthquakes with epicentral distances of less than 10 km; thus, the seismic waves incident upon the seismic station are near vertical. In this case, the wave forms are not sensitive to the crust structure but sensitive to near-field effects; hence, the half-space full wave form Green's function (Johnson, 1974), which were provided by rapid computation, are indeed suitable for wave form inversion as used in this study.
For the ability as shown in the numerical test ( Figure 5), therefore, the accurate source parameters can be determined from the sparse stations, and this will greatly improve the information that can be recovered from portable instrumentation, such as a PDAS recorder (with Force Balance Accelerometer). The latter might be used for monitoring aftershocks and for studying the mechanism of earthquake sequence. For several portable stations deployed in the aftershock area, small earthquakes with a local magnitude of less than 2 can be clearly recorded, and the source parameters which are usually impossible to retrieve from a local' short-period network can be determined by near-field wave-form inversion.
Another potential application of the near-field wave form inversion is the rapid determi nation of source parameters for broadcasting by the CWB. For example, the seismograms of the May 28, 1992 event were retrieved from the low gain channel of CWB network. These data can be dialed up by a telephone system and can provide a rapid estimation of source mechanism only a few minutes after an earthquake.
For the half-space model considered in this study, only the near-source station was selected for the moment tensor inversion. However, when detailed crust velocity structure is available, using the Green's functions developed for the layer model, either the generalized ray method (Helmberger and Harkrider, 1978) or the wavenumber-frequency integration method (Wang and Herrmann, 1980), or the lateral inhomogeneous model; either the finite difference method (Vidale et al., 1985) or the finite element method (Huang, 1992), may be applied to invert the source mechanism from the seismograms of the regional earthquake or the very shallow events without any modifi cations. Otherwise, with the least-squares inversion used in this study, the effects of noise on inversion is commonly addressed by examining the covariance of model space. The resolution kemal of Green's function is the key for accurate analysis of results (Aki and Richards, 1980). To completely define accuracy of results, more detailed analysis is nt'.Cessary in the near future.

CONCLUSIONS
The waveform inversion scheme developed in this study is a fast way to retrieve the moment tensor of a point source from the near-field strong motion data. The method has demonstrated its utility in retrieving source parameters from the three-component wave forms of a few sparsely-distributed or single-station strong motion records. Numerical tests have shown that this method can be used to retrieve source mechanism from noise data which included random noise of less than 30% of peak signal amplitudes. The inversion of the single station (CHK) record of the May 28, 1992 event revealed that it can be used to rapidly determine the source mechanism of a local earthquake.