-
Notifications
You must be signed in to change notification settings - Fork 16
Method
ISOLA is based on multiple point-source representation and iterative deconvolution method, similar to Kikuchi and Kanamori (1991) for teleseismic records, but her 8000 e the full wavefield is considered, and Green's functions are calculated by the discrete wavenumber method of Bouchon (1981). Thus the method is applicable for regional and local events. Instrumentally corrected band-pass filtered velocity records are used. The code transforms velocity into displacement, inverts the displacement, and provides synthetic displacement. Seismogram s(t) is approximated by a combination of elementary seismograms ei(t), i=1,2,...6, i.e. s(t)= sum(ai.ei(t)). The elementary seismograms correspond to basic focal mechanisms as follows:
i strike dip rake 1 0 90 0 2 270 90 -90 3 0 90 90 4 90 45 90 5 0 45 90 6 isotropic
Seismograms (3-component), from a set of stations, represent the data, d, while the coefficients ai are the parameters to be found, m. The linear inverse problem d= G m is solved by the least-square method: m = (GT G)-1 GT d, where GT is G transposed and ()-1 is the inverse of the so-called system matrix (GT G). Eigenvalues of the system matrix (GT G) are calculated to understand how well-posed (or ill-posed) the inverse problem is. Ratio of min/max eigenvalue is reported. The lower the value, the worse the results are (although, formally, the quality of the fit may be good, e.g. variance reduction can be large). Question how the individual stations contribute to the well-posed inversion has no simple answer. It is to be carefully analyzed case by case.
Optionally, the weighted inversion is performed, with weights prescribed by user, different weights for different stations and components. As a rule, up-weighting records with small amplitudes, and down-weighting large amplitudes, is recommended. Another strategy is to up-weight stations whose seismograms appear to be the most sensitive with respect to small modifications of the problem.
The coefficients ai are related to moment tensor Mpq (where p,q refer to geographic coordinates x (>0, N), y (>0, E), z (>0, up)):
Mxx=-a4+a6 Myy=-a5+a6 Mzz=a4+a5+a6 Mxy=Myx=a1 Mxz=Mzx=a2 Myz=Mzy=-a3
Moment tensor components in spherical coordinates r, t(=theta), p(=phi), used in Harvard catalogue, or in GMT, are given by:
Mtt=Mxx Mpp=Myy Mrr=Mzz Mtp=-Mxy Mrt=Mxz Mrp=-Myz
Eigenvectors of the moment tensor provide strike, dip and rake. Eigenvalues provide scalar moment Mo and decomposition of the moment tensor into three parts: double-couple (DC), compensated linear vector dipole (CLVD), and voluminial (VOL).
Possible moment tensor (MT) inversion modes are as follows:
-full MT inversion: all six basic focal mechanisms, 1,2,...6; DC+CLVD+VOL -deviatoric MT inversion: basic focal mechanisms 1,2,...5; DC+CLVD; VOL%=0; -DC-constrained MT inversion: DC only; VOL%=CLVD%=0 -known and fixed DC moment tensor (only position and time is searched)
Opposed to strike, dip, rake and Mo, the DC% is strongly unstable. It closely relates with the fact that a large deviation from pure DC has, as a rule, a very small effect upon the seismograms. Indeed, it can be demonstrated by many examples in which synthetics for, say DC%=50 and DC%=100, cannot be distinguished by eye from each other, and quantitative measure of their fit with real data is nearly identical (e.g. variance reduction, introduced below, differs by 1%). Therefore, we can say that the DC% value provided by the routine moment tensor inversion is rarely of any physical meaning. The same applies for the volume percentage VOL%.
Optionally, the inversion can be constrained through the Lagrange multipliers to give a nearly 100% DC solution. The (non-linear) condition that the moment tensor must have zero-valued determinant is solved iteratively. It slows down the calculation and may create problems. There is no special advantage in this option, because the strike, dip and rake of the DC-constrained solution is usually very near to that of the (unconstrained) deviatoric MT inversion.
Thus the most recommended inversion mode is deviatoric inversion (VOL%=0) in which, however, the resulting DC% is considered with great caution.
It is assumed in the above explanation that the source position and time are known. This code allows their optimization through a grid search over set of trial source positions and time shifts, pre-defined by the user. The search seeks the trial source position and time for which the residual error of the least square method is minimized. It is equivalent to maximize correlation between the observed and synthetic seismograms. The concept of optimizing position and time is applicable to each point element (the so-called subevent) of the generally more complex source model. If the whole earthquake is represented (in the low-frequency range) by a single source, the optimum solution provided by the grid search corresponds to the centroid (the centroid position, centroid time and the centroid MT).
User can check the correlation values as a 2D function of the trial source position and time. The code has an interactive part in which user can use automatic search of the optimum correlation, or can make his preferred choice. The preferred choice may be the source position or time with a somewhat lower correlation, close to formal maximum, in which however some other conditions (constraints) are better satisfied. In this way user can apply, for example, his knowledge about the first-motion polarities, geologically deduced fault dip, various assumptions about rupture propagation, etc. The 2D correlation may also help in identifying poor resolution of some parameters, for example the down-dip source position.
Narrow band-pass inversions around period T can often provide two correlation maxima, in time t and t+T/2, characterized by rake and rake+180° respectively. This is just the case where even a single guaranteed first motion polarity may help to decide between proper time, t or t+T/2.
The code works iteratively. A point source contribution is called subevent. Using data code searches for subevent 1. Once subevent 1 is found, its contribution is subtracted from the data to get the so-called residual data. In the next step the residual data are processed in the same way as the original data during the first step. Subevent 2 is found, etc. Original data minus the last residual seismogram provide the resulting synthetic seismogram.
Optionally, the search of spatial and temporal distribution of subevents can be performed with fixed (100% DC) moment tensor. The user prescribes the known (or assumed) strike, dip and rake, and these values are kept fixed for all subevents. It can be recommended when the unconstrained moment tensors vary chaotically from one subevent to the other.
When the solution is represented by several subevents of different focal mechanism, the resulting scalar moment cannot be obtained as the sum of the subevent scalar moments. That is why we also provide an output file where, subsequently, the subevent contributions are summed up as tensors and the scalar moment of the resulting (cumulative) tensor is shown.
Denoting data as d, and synthetics as s, the fit between them is measured by the so-called variance reduction varred= 1 - |d-s|2/|d|2, where |.| denotes the L2 norm, e.g., |d|2=suma(di2). The norm is calculated as a single number summing up over all time series, all components, and all stations. The optimum value of varred is 1. Wrong synthetics may give even varred < 0 (no physical meaning).
The variance reduction is determined repeatedly, “when building up” the synthetics, i.e. after adding every subevent. In other words, when s1, s2,… are synthetics corresponding to individual subevents, then first we calculate varred from d and s=s1, then from d and s=s1+s2, etc. As such, varred is a measure of the convergence of the iterative deconvolution process. Slow increase of varred with adding subevents indicates that additional subevents do not further improve the fit.
At this point it is useful to clarify how the variance reduction relates with the correlation. The relation is simple when all stations and all components are used, and all weights equal 1, and when we calculate only the first subevent. Take the observed seismogram as data d, and the first subevent as synthetics, s=s1. It is an internal property of the least-square minimization of |d-s|, during the corresponding matrix operations, G, GT G, etc., that we automatically get the value of c2, where c is the correlation between d and s, without numerically calculating the correlation integral. (Speaking more precisely, we get sqrt(c2)=abs(c), i.e. not the correlation c, but only its absolute value. Anyway, hereafter, we simply say only “correlation”). At the same time, it can be shown that in c2= 1 - |d-s|2/|d|2 = varred.
When we calculate already the second subevent, s2, the least-square method minimizes |r-s2|, where r is the residual seismogram defined as difference between the original data and the first subevent, d-s1. Therefore, from this least square run, we get automatically, and we report on the code output, (the absolute value of) the correlation between r and s2. At the same time, after adding the second subevent, the code calculates varred using d and s=s1+s2, i.e. varred= 1- |d-(s1+s2)|2/|d|2. Obviously, the varred value reported after the second subevent is not the same as the (squared) correlation between r and s2. Moreover, this varred is also not the same as the (squared) correlation between d and s=s1+s2, because we did not find the synthetic s=s1+s2 by a single least-square minimization of |d-(s1+s2)|.
When non-unit weights are used, and/or some stations are deselected from a particular inversion run, the relation between correlation and variance reduction gets even more complicated.