Probabilistic local waveform inversion for moment tensor and hypocentral location

SUMMARY 
 
We have developed a probabilistic method to determine simultaneously earthquake-source mechanism, hypocentral location and source time function (STF) from the inversion of short-period waveform data of weak local events (epicentral distance less than 100 km). The procedure takes into account the effects of the random noise contained in the seismograms, the uncertainty of the hypocentre determined from arrival times and the inaccurate knowledge of the velocity structure, while estimating the error affecting the derived focal parameters. 
 
 
 
The forward modelling of the seismic ground motion is computed for frequencies up to 5 Hz in horizontally layered elastic media. Assuming that all uncertainties can be described by Gaussian probability density functions (PDFs), a priori information, measurement errors and theoretical errors are estimated. Then, the proposed probabilistic approach maps the posterior PDFs for the hypocentral coordinates, the moment tensor, and the STF. The final estimates for the focal parameters are given by the maximum likelihood points of the PDFs, while solution uncertainties are presented by scatter density plots. The estimated uncertainties in the moment tensor components are plotted on the focal sphere in such a way, that the significance of the double couple, the compensated linear vector dipole and the volumetric parts of the source can be assessed. 
 
 
 
Tests of the method were first performed on synthetic data generated for a configuration of a local seismic network in Hungary. Simple point sources as well as complex mechanisms were successfully reconstructed. In the next step, real data from the region were inverted to obtain the mechanisms of three local earthquakes. The resulting moment tensor solutions have negligible volumetric part, implying the tectonic nature of the events. The retrieved source mechanisms are in agreement with the available clear readings of first-arrival P-wave polarities and with the main stress pattern published for the epicentral region.


I N T RO D U C T I O N
Since the work of Gilbert & Dziewonski (1975), who inverted free oscillation data to retrieve the source moment tensor of large earthquakes, seismic waveform inversion has routinely been applied to teleseismic observations of earthquakes in order to determine their focal mechanisms and source parameters. Provided that a structural model of the region under study is available, synthetic seismograms generated by various types of source can be constructed and compared with observed records in a suitable inversion scheme. The main advantage of such a procedure, in comparison with the classical body-wave polarity studies, is that complete seismograms are used in the analysis. This allows the retrieval not only of the faultplane solution and hypocentral location, but also of the source time history.
Several methods have been developed to simultaneously determine the point source mechanism and source time function (STF) of teleseismic events. Techniques for inverting normal-mode data (Gilbert & Dziewonski 1975;Buland & Gilbert 1976), surface waves (McCowan 1976;Aki & Patton 1978;Kanamori & Given 1981), and body waves (Dziewonski et al. 1981;Langston 1981;Sipkin 1982;Oldenburg 1982) have been presented. For regional earthquakes the inversion of very broad-band seismograms has been proposed (Fan & Wallace 1991;Zhao & Helmberger 1994;Ritsema & Lay 1995). Some of these approaches allow determination of the source depth as well. Sophisticated source retrieval methods also exist for the inversion of near-source strong motion recordings (Ji et al. 2002;Sekiguchi et al. 2002). These approaches make use of near-source data recorded by dense accelerometer arrays.
When studying local seismicity, we deal with seismic events that occur in local epicentral distances. In areas of low seismicity the weak local events provide the only key to determine fault parameters and small-scale tectonic structure. These earthquakes are usually too small to provide enough information (above noise level) at long periods but are well recorded in the relatively high-frequency band. For this type of earthquakes, however, the seismic network is usually not dense enough and the knowledge of the medium is not detailed enough to allow the reconstruction of a finite-source model as in near-source studies. In comparison with teleseismic waveforms, seismograms of local events contain much higher frequencies making the application of teleseismic approaches dubious (Koch 1991a,b). The use of higher frequencies requires the use of more detailed models of inhomogeneous media for which the synthetic seismograms, specifically the Green's functions, should be computed. Insufficient knowledge of the medium is detrimental to the processing of high-frequency signals, because short wavelengths are more affected by small-scale heterogeneities of the medium, which are only rarely known well.
Attempts have been made to retrieve source parameters from weak, local, and noisy waveforms. Mao et al. (1994) have proposed a linearized algorithm. Sileny et al. (1992) used the method of modal summation to compute Green's functions in a vertically inhomogeneous medium with very fine structure. To improve the results, they used dynamic relocalization of the depths of the events, that is, they used amplitudes to relocate event depths determined from the kinematics. The reliability of their solution can be assessed by the procedure of Sileny et al. (1996) or Sileny (1998). The method has been successfully applied for the inversion of weak volcanic earthquake data (Panza & Sarao 2000;Sarao et al. 2001).
In this study we introduce a probabilistic non-linear waveform inversion method to retrieve the source parameters of weak local earthquakes. Assuming that all uncertainties can be described by Gaussian probability density functions (PDFs), we estimate the posterior PDFs for the hypocentral coordinates, the moment tensor components and the STF. The final estimates for the focal parameters are given by the maximum likelihood points of the PDFs, while solution uncertainties are presented by scatter density plots. Scatter density plots are obtained by drawing samples from the PDF with the number of samples proportional to the probability. Then, the estimated uncertainties in the moment tensor components are plotted on the focal sphere in such a way, that the significance of the double couple (DC), the compensated linear vector dipole (CLVD), and the volumetric parts of the source can be assessed.
The method is thoroughly tested by the inversion of synthetic data sets generated for simple point sources as well as for complex mechanisms. Then, it is applied to three local events that occurred in the central part of the Pannonian basin.

T H E I N V E R S I O N P RO C E D U R E
Using the moment tensor representation for body waves, the jth component of the displacement field for a point source approximation may be written as where M kl (t) are the components of the time dependent moment tensor, G jk,l (r, h, t) the elastodynamic Green's functions, r is the position of the receiver, and h the hypocentre (Aki & Richards 1980, eq. 3.22). The symbol denotes temporal convolution. If only farfield terms are used for the description of the seismic displacement, then only those terms in the Green's functions are left, which are coupled toṀ kl (t), the time derivative of the moment tensor. Owing to the symmetry ofṀ kl (t), the far-field displacement can be written as (Sipkin 1982) where m(t) is the vector containing the six independent moment tensor rate functions (MTRFs), that is, and the g jk (r, h, t) are the corresponding Green's functions. Most inversion schemes depend on this linear relationship. For weak events, when the focal mechanism is considered as constant in time during the rupture process, all of the m k (t) components have the same time dependence s(t). Then the far-field displacement seismogram is where m is the vector containing the six independent, time-invariant elements of the moment tensor and s(t) is the STF. According to eq. (2), if the Green's functions (i.e. the velocity structure and hypocentral location) are known, the earthquake waveforms can be inverted for the unknown MTRFs by a linear procedure. Using a suitable forward modelling method, the synthetic Green's functions are computed for the given earth model and hypocentral coordinates. Once the Green's functions have been determined, the MTRFs are estimated by a least-squares linear inversion of eq. (2). In this study we apply the truncated singular value decomposition (TSVD) algorithm for this purpose (see, e.g. van der Sluis & van der Vorst 1987; Xu 1998). Alternatively, one can also use the method of Sipkin (1982) extended to local waveforms by Koch (1991a).
In case of short epicentral distances, routinely determined hypocentres are usually not accurate enough to be used in focal mechanism inversion. Therefore, in the first step of the proposed inversion procedure, we perform a Bayesian inverse calculation to estimate the hypocentral location from waveform data. As has been developed by Tarantola (1987), the Bayesian solution of this nonlinear inverse problem is the a posteriori PDF of the model parameters (hypocentral coordinates), which is denoted by σ (h). This PDF is the product of two terms. The a priori PDF ρ(h) incorporates information about the hypocentral position that is independent of the observed data. The second term, the likelihood function L(h), measures the fit between the observed data and data predicted from the model parameters. If the Gaussian model is an acceptable approximation for all uncertainties of the problem, the a priori distribution is given by while the likelihood function is given by for a given hypocentre. Forward modelling consists of two steps. Firstly, the MTRFs are estimated for the given hypocentre by linear inversion. To be consistent with the Bayesian approach, eq. (2) is inverted using a weighting matrix W, for which W T W = C −1 D . In the second step, the estimated MTRFs are used to calculate the predicted waveform data.
In this study, the software package NonLinLoc (Lomax et al. 2000) is applied to estimate the a priori PDF ρ(h) using arrival times. Then the a posteriori PDF σ (h) is mapped by the oct-tree importance sampling algorithm (Lomax & Curtis 2001) using waveform data. The oct-tree algorithm gives accurate, efficient and complete mapping of the posterior PDF. It uses recursive subdivision and sampling of cells in three dimensions to generate a cascade of sampled cells, where the number of sampled cells follows the values of σ (h) at the cell centre, thus leading to higher density of cells in areas of higher σ (h). The final hypocentral location is given by the maximum likelihood point of the a posteriori PDF, while location uncertainties are presented by scatter density plots. Scatter density plots are obtained by drawing samples from σ (h) with the number of samples proportional to the probability, that is, high density of samples denotes high probability. Given the samples of the scatter plot, several statistical properties of the hypocentral coordinates can be determined, such as standard deviation, mean or median.
For all points in the hypocentre density plot, the MTRFs are calculated as well. Their distribution represents the uncertainty of the MTRF due to that of the hypocentre. However, measurement errors and modelling errors lead to MTRF uncertainty even for a fixed source position. In the second step of the proposed non-linear inversion procedure, this type of uncertainty is estimated by the bootstrap resampling technique (Efron & Tibshirani 1986). The key concept of bootstrapping is that the original data set is resampled to form a large number of new data sets. After inverting the newly generated bootstrapped data sets, the resulting multiple estimates of the model give information on model variance. This method is completely insensitive to the statistical properties of the data and can equally be used for both linear and non-linear problems. The short discussion of the mathematical basis of bootstrapping and its application to a geophysical inverse problem can be found in Tichelaar & Ruff (1989).
Thus, a large number of bootstrapped data sets are generated and linearly inverted using hypocentres randomly chosen from the posterior hypocentral PDF σ (h). The MTRF solutions obtained in this way can be considered as samples from the posterior PDF of the MTRF, so statistical properties, such as the maximum likelihood point, the mean, or the median, can be determined. However, when the Green's functions are deconvolved from the seismograms, the features of the records not modelled in the Green's functions appear in the independent MTRFs as not correlated signals across the set of the six functions. As a result, the TSVD inversion of eq. (2) yields, in general, focal mechanisms that are variable in time and, in particular, may contain backward slip motion.
As a consequence of the assumption that the focal mechanism of weak events is constant in time, the MTRFs are modelled by the product of their correlated part and six multipliers (Sileny et al. 1992;Sileny 1998). The correlated part represents the STF, while the multipliers describe the six moment tensor components of the average source mechanism (eq. 4). This decomposition of the MTRFs is a non-linear problem and solved iteratively by imposing constraints such as positivity of the STF and the requirement of a mechanism consistent with clear readings of first-arrival polarities. The positivity constraint guarantees a physical sense of the STF, because it allows only forward slips. As Kravanja et al. (1999) point out, considering the MTRFs as independent functions during the TSVD inversion leads to an overparametrization of the problem, which is advantageous for reducing the effects of poor modelling of the structure, while the subsequent decomposition of the MTRFs reduces the bias due to non-exact Green's functions since it keeps only their coherent part.
According to the above discussion, in the last step of the proposed inversion procedure, each MTRF solution is decomposed into the product of the time-invariant moment tensor and the corresponding STF. The distribution of the resulting moment tensors and STFs approximates well the uncertainty of the focal mechanism. Once the moment tensors are retrieved, their principal axes are deduced. Then each moment tensor is decomposed into an isotropic part, representing an explosive or implosive component, and into a deviatoric part, containing both the DC and the CLVD components. Finally, the resulting mechanisms are plotted on the focal sphere in such a way, that the statistical significance of the DC, CLVD and isotropic parts of the source can be assessed.
The different stages of the proposed inversion scheme are shown in Fig. 1.

P R A C T I C A L C O N S I D E R AT I O N S
Efficient treatment of the forward problem is crucial for any inversion problem. As it is performed many times, it is essential that the calculation be as fast as possible. Since, in the present work, we use horizontally layered earth model to describe the real medium, elementary Green's functions (ELGFs) are constructed by a propagator matrix-wavenumber integration method (Wang & Herrmann 1980;Herrmann & Wang 1985), which allows for the calculation of the full waveform at high frequencies and short epicentral distances. The ELGFs, which are independent of the event-station azimuth, are determined for a discrete set of source depth and epicentral distance and stored on disc. In the course of the inversion, for intermediate values of depths and epicentral distances, new ELGFs are computed using linear interpolation. Finally, the Green's functions needed for the inversion of eq. (2) are determined by a linear combination of the interpolated ELGFs taking into account the event-station azimuth as well. As the ELGFs are calculated only once for a given earth model, the time required to treat the multiple forward problem does not depend on the complexity of the medium. If the strategy of precalculating and storing the azimuth independent ELGFs were not used, the synthetic Green's functions would to be calculated each time when needed and this would significantly prolong the total time of the inversion. In this study, when calculating the ELGFs, we use discretization intervals of 1 and 5 km for source depth and epicentral distance, respectively.
In the course of the first step of the proposed inversion procedure, the data covariance matrix C D in eq. (6) has to be estimated. This matrix represents the measurement errors of the observed waveforms and errors in the calculation of the theoretical seismograms. Measurement errors can be estimated from the observed waveforms themselves, while the most appropriate estimate of the modelling errors due to mismodelling of the medium can be achieved through synthetic experiments. Our simplistic approach consists of generating synthetic seismograms for several reasonable earth models and calculating their variance. This variance gives the rough estimate of the theoretical error part of C D .
Model parameter uncertainties of probabilistic inversion results are often shown as confidence contours, which can be estimated from the posterior PDF. This is only possible, however, if the complete PDF is available, as obtained by, for example, a grid-search algorithm. In global sampling methods such as the bootstrapping technique, the complete PDF is not available and the uncertainties of the solution are shown by scatter density plots. The proper illustration of the density plots is very important from the interpretational point of view. It is especially true for the visualization of the moment tensor uncertainties.
In this study, the method of Riedesel & Jordan (1989) is employed to display the scatter plot of the moment tensor solution. The principal vectors of a moment tensor define the tension (T), neutral (N) and compression (P) axes, while the principal values give their magnitudes. In the principal axis system, various unit vectors can be constructed using various linear combinations of the principal vectors. The vector that describes a general source mechanism is MT, a DC source mechanism has the vector representation DC, the vector corresponding to a purely isotropic source is the vector ISO, and two possible CLVD vectors can also be defined (Jost & Herrmann 1989). The scatter plot for the MT vector, together with the DC, ISO and CLVD vectors corresponding to the best moment tensor solution are then plotted on the surface of the focal sphere. The great circle that connects the DC and CLVD vectors on the unit sphere defines the subspace on which MT must lie for a deviatoric source. The distribution of the MT scatter plot with respect to the DC, ISO, and CLVD vectors informs us on the significance of the DC, ISO and CLVD parts of the solution: if the vector DC lies within the MT scatter plot, the mechanism is a DC; for a reliable CLVD solution, the MT scatter plot lies on top of one of the CLVD vectors; and when the MT scatter plot lies off the deviatoric great circle, the isotropic portion is reliable.
The horizontally layered earth model used in this study is only an approximation of the true medium. The effects of inadequacies of the structural models have been investigated by Sileny et al. (1992) and Kravanja et al. (1999). They have proved, by synthetic tests, that a poorly known velocity structure: (1) introduces apparent non-DC components in the moment tensor solutions, contaminating particularly the CLVD component; (2) influences the orientation of the DC; (3) leads to spurious peaks in the STF and (4) may introduce 20-30 per cent error in the scalar moment.
However, with the error analysis that we apply, it is possible to discriminate spurious non-DC components from physically meaningful solutions. Focal depth may also vary with velocity structure within about 2-3 km (see, e.g. Fan & Wallace 1991), but this variation still provides a much better focal depth constraint when compared with that derived from the inversion of arrival times read from sparsely distributed seismograph stations, which is the case in Hungary.
If P, SV and SH amplitudes and arrivals are known, two threecomponent stations can, in principle, determine the system of equations to be solved in the inversion. However, as is well known, the horizontal components can be sensitive to site effects not adequately modelled by a horizontally layered earth structure. Unless the velocity structure is adequately known, the simultaneous waveform match on all three components is almost impossible (Saikia & Herrmann 1985). Moreover, Koch (1991a) points out that the inversion for the MTRFs becomes unstable when using more SH seismograms than absolutely necessary, because they make the inversion problem ill-posed. For these reasons, whenever the available data make it possible, we use only vertical-component seismograms in the present work. If only vertical-component waveforms are used, a minimum of three stations are required to obtain the complete MTRFs (P and SV data for each station). Nevertheless, even if the number of stations is sufficient, the inversion problem can be illposed, so as many seismograms should be used as possible with good signal-to-noise ratio.

T E S T S O N S Y N T H E T I C D ATA
In this section, we apply the proposed inversion method to synthetic data generated for the station configuration shown in Fig. 2. The stations belong to the Hungarian short-period seismic network operating near a nuclear power plant. The hypothetical hypocentre coincides with that of a real earthquake that occurred near Szabadszállás, Hungary (Fig. 13). This event will be thoroughly discussed in the next section.
In the course of the tests, synthetic displacement seismograms are calculated for a set of model parameters (hypocentral coordinates, source mechanism and STF) and for an assumed earth model. These seismograms are considered as 'observed' signals, and the model parameters used in the computation are the 'correct' model. Then the waveforms are inverted for the model parameters as described in the previous sections.
For the generation of the 'observed' signals, we use the horizontally layered earth model specifically constructed for the inversion of local waveforms recorded in Hungary (Fig. 14). The waveforms are calculated by a propagator matrix-wavenumber integration method (Wang & Herrmann 1980;Herrmann & Wang 1985) assuming triangular STF with a total duration of nine samples (about 0.2 s). White noise with variance reaching 10 per cent of the peak amplitude of the signal is superimposed on the synthetic seismograms, which are then filtered by a causal bandpass filter from 1 to 5 Hz. In order to estimate the prior hypocentral coordinates, the P arrivals are picked up and kinematic localization is performed by NonLinLoc (Lomax et al. 2000). During the waveform inversion procedure, the Green's functions are calculated using the same velocity model as that used for generating the observed seismograms.

Simple sources
The first source used for generating the observed records is a vertical strike-slip with strike 90 • , dip 90 • and rake 0 • , situating at a depth of 13 km. Figs 3-5 illustrate all information the proposed inversion method can yield.
The MTRFs in Fig. 3 show some fluctuations due to the noise contaminating the signals and to modelling errors originating from the fact that the Green's functions are calculated by interpolation. At the 95 per cent confidence level, however, only the m 2 (=Ṁ 12 ) component differs significantly from zero. The synthetic waveforms computed using the best (maximum likelihood) source parameters are very close to the observed seismograms (Fig. 4). The retrieved mechanism is also in good agreement with both the firstarrival P-wave polarities and the strike-slip mechanism to be reconstructed (Table 1). The histograms and scatter density plots of the source parameters are shown in Fig. 5. The MT scatter plot is concentrated around the DC vector and the best MT vector practically coincides with the pure DC solution. The principal axes are also well defined. The only peak of the STF closely resembles the filtered triangular STF used to generate the observed seismograms. Summing up, the inversion of the vertical-component waveforms has been able to successfully reconstruct the correct source parameters.
The other simple source used to test the proposed waveform inversion method is a thrust faulting mechanism with strike 315 • , dip 45 • , rake 90 • , and focal depth 13 km (Fig. 6a).
When inverting the vertical seismograms generated by this thrust faulting source, the fit of the maximum likelihood mechanism with the true one is worse than in the previous strike-slip case (Fig. 6b). The result does not become better even if horizontal-component waveforms are included in the inversion (Fig. 6c). The poorer solution may be the consequence of the less pronounced azimuthal dependence of amplitudes compared to the strike-slip mechanism. Note that for the strike-slip event we have data from all four quadrants of the P-wave radiation pattern, while for the thrust faulting event most of the data correspond to only two dilatational zones However, if the inversion is performed with the constraint that the solution be purely deviatoric, the resulting mechanism becomes very close to the true one, even if only vertical components are used (Fig. 6d, Table 1). The uncertainty of the solution is larger than that for the strike-slip mechanism (Fig. 7). Both the DC vector and one of the CLVD vectors lie within the MT scatter plot, implying that the solution may have considerable amount of CLVD component. Including horizontal-component seismograms in the inversion does not change this result. The best MT vector, however, practically coincides with the pure DC solution.
Closer inspection of the solution reveals that the mechanisms with considerable CLVD component in the scatter plot correspond to focal depths less than 12 km. Since the earth model has a velocity contrast at that depth (Fig. 14) and the correct source is situated below this layer boundary, the modelling error for hypocentres above 12 km becomes significant leading to significant amount of spurious CLVD component. Compared to the previously examined strike-slip source, for the given acquisition geometry and source orientations, the investigated thrust faulting event seems to be more sensitive to this type of error. This example demonstrates that imposing deviatoric constraint during the inversion may be necessary even if the number of the recording stations seems to be sufficient to perform unconstrained inversion.

Complex sources
Any source must be considered as spatially extended, unless the dominant wavelength observed is much larger than the fault extension. On short-period local seismograms, however, extended sources produce complex waveforms. In this subsection, we wish to discover whether the point-source approximation used by the pro-posed inversion method can yield useful results for such complex sources.
To model the effects of extended sources, two subevents are considered with the same seismic moment at hypocentral depth of 13 km. The first subevent has a vertical strike-slip mechanism with strike 90 • , dip 90 • , and rake 0 • (Fig. 8a). The second subevent is situated in the slip direction of the first one. The spatial distance  between the two subevents is 0.25 km, while the time offset is 0.1 s, resulting in a rupture velocity of 2.5 km s −1 . The ratio between the rupture and S-wave velocities is then 0.7 (v S = 3.58 km s −1 ). The mechanism of the second subevent is defined by strike 70 • , dip 70 • and rake 20 • (Fig. 8b).
Since the STF has a time duration of 0.2 s, the bandpass filtered seismograms cannot resolve the complex nature of the above described event. In this case the result of a waveform inversion should be a combination of the two subevents. This composite mechanism can be approximated by a DC component with strike 80 • , dip 77 • , and rake 9 • in addition to an about 13 per cent CLVD part (Fig. 8c, Table 1).
The solution of the proposed inversion technique approximates well the composite mechanism of the two subevents (Fig. 8d,  Table 1). Because the first subevent has just the same strike-slip mechanism as the one investigated in the previous subsection, it is worth comparing the two results (Figs 5 and 9). The uncertainties of the hypocentral coordinates and scalar moment do not differ considerably for the two events. The inversion method is certainly not able to resolve the different hypocentres of the subevents. The uncertainties of the MT vector and the principal axes for the complex event, however, are larger than those for the pure strike-slip mechanism. According to our expectations, the STF has been broadened, but the two subevents cannot be resolved.
When we enlarge the spatial distance between the subevents from 0.25 to 0.5 km, the second subevent is delayed by 0.2 s with respect to the first one. This time delay is as long as the STF, leading to more complicated waveforms than in the previous case. The inversion results are illustrated in Figs 10-12.
Due to the more significant modelling errors (point-source approximation of an extended source), the inverted MTRFs are far more complex than those calculated for the pure strike-slip event (cf. Figs 3 and 10): their oscillatory character is more pronounced and the m 1 (=Ṁ 11 ) and m 3 (=Ṁ 22 ) components have considerable amplitude. Nevertheless, the observed seismograms are fitted very well and the maximum likelihood mechanism is in agreement with the composite mechanism of the two subevents (Fig. 11, Table 1). Comparing the waveforms in Figs 4 and 11 reveals the waveform complexities an extended source can produce.
The uncertainty of the mechanism is even larger than in the previous complex source example (Fig. 12). The well-defined principal axes, however, allow only a small variation of the orientation of the mechanism. The maximum likelihood solution has some spurious volumetric component, but it has no statistical significance. The STF has two peaks significant at the 95 per cent confidence level, the second peak certainly corresponding to the second subevent. This result suggests that multiple peaks in the STF may be the sign of a complex source if its subevents are separated well enough in time.

A P P L I C AT I O N T O O B S E RV E D WAV E F O R M S
In this section we illustrate the proposed method through the waveform inversion of three local earthquakes that occurred in the Hungarian part of the Pannonian basin (Fig. 13).
The waveform data used in this study were recorded by the threecomponent seismological stations shown in Fig. 13. The transfer functions of the velocity instruments are flat above 1 Hz. To remove  low-and high-frequency numerical effects from the records, a causal bandpass filter from 1 to 5 Hz was applied to the time series after transforming them to displacement. The same filter was applied to the displacement Green's functions. For the generation of the synthetics, a propagator matrixwavenumber integration method (Wang & Herrmann 1980;Herrmann & Wang 1985) was used, which allows the calculation of the entire wavefield for horizontally layered earth structures. The construction of a suitably detailed earth model required for the inversion of local waveform data is not an easy task. For this study a rather simple velocity model was derived from measured traveltimes of local earthquakes and controlled seismic sources (Fig. 14).
When the number of the recording stations made it possible, we used only vertical-component records, since the horizontal components are usually more sensitive to site effects not adequately modelled by a simple horizontally layered earth structure. The time window for the seismogram to be processed starts at the arrival of the P phase. Its duration is approximately limited to the first pulse followed by some converted phases. The length of the time window was shortened for some stations when it became evident that the latter parts of the seismograms had not been recovered satisfactorily. Only seismic stations with first P g arrival and good signal-to-noise ratio were used in the inversion process.
In the course of the routine observatory data processing, the hypocentral locations of the selected events were determined by the standard HYPO71 algorithm (Lee & Lahr 1975) using a simplistic two-layer crustal model (Tóth et al. , 1997(Tóth et al. , 1999, while for estimating the fault-plane solutions the computer program FPFIT (Reasenberg & Oppenheimer 1985) was used. Both the routinely calculated source parameters and those determined by the proposed waveform inversion method are summarized in Table 2.
The waveform inversion results for the 1996 March 28, Szabadszállás event (M L = 3.0, Event #1 in Table 2) are presented in Figs 15 and 16. Fig. 15 compares the observed seismograms and synthetic waveforms computed using the best (maximum likelihood) source parameters and the earth model illustrated in Fig. 14. Due to the good signal-to-noise ratio and the azimuthally well-distributed recording stations, the waveforms are fitted very well. The retrieved mechanism is in agreement with the available clear readings of firstarrival P-wave polarities.
The histograms and scatter density plots of the source parameters obtained by the proposed waveform inversion procedure are shown in Fig. 16. Both the hypocentral location and the scalar moment are obtained with low uncertainty. The STF has two peaks significant at the 95 per cent confidence level. The main energy release lasts about 0.15 s, while the second one is concentrated around 0.4 s with about 0.1 s duration. The reliability of the second peak is much less than that of the main one. The MT scatter plot is fairly small in dimension, that is the source parameters are determined with high reliability. It contains the DC vector, thus, a pure DC may be the solution of the inversion. The tightly confined zones of the principal axes allow only a small variation of the orientation of the mechanism. The best DC part of the maximum likelihood moment tensor agrees well with the fault-plane solution obtained by the FPFIT algorithm (Table 2). It should be emphasized here that while the classical FPFIT method seeks the best pure DC mechanism satisfying the observed P-wave polarities, the waveform inversion method presented in this paper does not constrain the solution at all: it may contain CLVD and isotropic components as well. However, the MT scatter plot in Fig. 16 suggests that the non-DC components have no statistical significance and the only reliable DC component has been proved to be almost identical to the FPFIT solution.
The inversion results for the local event that occurred on 1995 June 9, near Szabadszállás (M L = 1.6, Event #2 in Table 2)  presented in Figs 17 and 18. Due to the smaller signal-to-noise ratio and the poorer azimuthal station coverage, the uncertainty of the solution is larger than in the previous case. However, the observed seismograms are fitted fairly well and the maximum likelihood mechanism is in agreement with the available first-arrival P-wave polarities (Fig. 17). The DC vector lies within the MT scatter plot, thus, a pure DC mechanism may also be acceptable as the solution of the inverse problem (Fig. 18). In fact, the non-DC components are very small and statistically insignificant. The azimuth of the P axis is well constrained, it spreads about 20 • . The STF has practically only one peak significant at the 95 per cent confidence level, lasting about 0.1 s.
Because of the small number of polarity data, the FPFIT algorithm was not capable of producing reliable fault-plane solution for Event #2. However, the focal mechanism estimated by waveform inversion is almost identical to that of Event #1 ( Table 2). Considering that the two earthquakes occurred very close to each other, the similarity of the two mechanisms suggests that the two events were generated on the same fault system by the same stress pattern. Actually, the P axis strikes NE-SW for both events, coinciding with the general trend of the stress field in the epicentral region Gerner et al. 1999). The 1998 July 3, Szalkszentmárton event (M L = 2.1, Event #3 in Table 2) was recorded by only four stations. The four verticalcomponent seismograms were not enough to define a well-posed problem, so two horizontal-component waveforms were also included in the inversion. To further improve the stability of the procedure, the inverse calculation was performed with the constraint that the solution be purely deviatoric. The results are illustrated in Figs 19 and 20.
The deviation between the recorded and calculated seismograms is largest for the horizontal components (Fig. 19). This is not surprising since, compared to the vertical components, the horizontal waveforms are usually more sensitive to velocity inhomogeneities not modelled by a simple 1-D earth model.
The oscillatory STF has two peaks significant at the 95 per cent confidence level, but the reliability of the second one is fairly low (Fig. 20). The considerable CLVD component in the maximum likelihood mechanism may be attributed to large modelling errors and/or the complex nature of the source. The CLVD vector, however, just lies off the MT scatter plot (Fig. 20), that is its statistical significance is very low. Moreover, the MT scatter plot touches the DC vector, so a DC may also represent the solution of the problem. The scatter plots of the N and T axes are stretched in such a way that a rotation of the mechanism around the P axis is possible. Nevertheless, the azimuth of the P axis is well constrained, so it can be used in research studies concerning the local stress field. Due to the small number of polarity data, the FPFIT algorithm was not capable of producing reliable fault-plane solution.

C O N C L U S I O N S
The method illustrated in this paper is able to retrieve simultaneously the hypocentral location, the STF, and the full seismic moment tensor of weak local earthquakes from waveform data even when few noisy records are available. The error analysis, carried out by the Bayesian approach and the bootstrap resampling technique, allows us to estimate and display the uncertainties of the source parameters. On the basis of the resulting moment tensor uncertainties, the statistical significance of the DC, ISO and CLVD parts of the solution can be readily assessed. The method is fully non-linear in that it does not require any linearization of the problem. Because the azimuth independent elementary Green's functions are pre-calculated and stored on disc, and the samples from the posterior PDF of the MTRFs are determined by a standard linear inversion technique, the proposed procedure is fairly efficient.  Tests on synthetic data indicate that the presented algorithm gives good results for both simple and complex sources. It should be pointed out, however, that when inverting waveforms for the optimum point source parameters, the resulting solution will be an average of the true source, in both space and time. Thus the solution does not necessarily correspond to the fault-plane solution determined by fitting P-wave polarity data, which represents the source mechanism at the initiation of rupture. Similarly, the hypocentral location retrieved by waveform inversion corresponds to the point where the main energy is radiated, while the hypocentre calculated from arrival times determines the point of rupture initiation. Moreover, for complex events, the average source mechanism may even disagree with some first-motion polarity data.
The method has also been applied to three local earthquakes that occurred in the central part of the Pannonian basin. The inversion of the short-period records was successful. The isotropic component of the moment tensor solutions for the selected events is insignificant, implying the tectonic nature of the events. The statistically significant part of the STFs lasts about 0.1-0.2 s. The retrieved mechanisms are in agreement with the available clear readings of first-arrival P-wave polarities. The principal axes of the source mechanisms also agree well with the main stress pattern published for the epicentral region.

A C K N O W L E D G M E N T S
The reported investigation was financially supported by the Hungarian Scientific Research Fund (No. T042572). The assistance of the Bolyai Research Scholarship is also acknowledged. The author is also grateful to Georisk Ltd. for providing the waveform data in the study. I am also thankful to Anthony Lomax and one anonymous reviewer for making various suggestions, which have helped me in improving this paper. Figures were prepared using the GMT software (Wessel & Smith 1998).