Recovery of the parameters of a focal zone by the data of earthquakes

We present an approach for solving the inverse kinematic problem of seismic with internal sources, based on the method of multidimensional data approximation on irregular grids. The times of arrival of elastic waves to the seismic stations are considered as known. The hodographs from earthquake to the stations are approximated for further determining the velocities of longitudinal and transverse waves using the eikonal equation. The ratio of these velocities determines the Poisson’s ratio, and the other elastic parameters of the medium can be found in units of the density. The results of implementation of the approach, based on the real data, are presented.


Introduction
Practical modeling in geophysics is based on very scarce information about the structure of media on the day surface and at the depths available for drilling. The structure of media at great depths, as well as in areas where drilling is not available or is expensive, can be discussed by indirect data using information about earthquakes as sources of elastic disturbances that can be measured at the day surface.
Investigations of the inverse kinematic problem of seismic (IKPS) has started in 1907 by G. Gerglotz and E. Wichert [1]. The posed problem is important both in the theory of inverse problems of mathematical physics [2], [3], and in practical geophysics [4].
Equations of motion in an anisotropic elastic medium [5] have the form where ρ is the density of the medium, u = (u i ) is the vector displacement field; c ijkl are the components of elasticity modules tensor with symmetries c ijkl = c jikl , c ijkl = c ijlk .
Here δ jk is the Kronecker symbol. We restrict ourselves to the following equations of motion in the isotropic elastic medium Both of these waves propagate in an elastic medium. The compression wave (v p ) leads to a relative change of volume and is equal to divu p . The displacement wave (v s ) is characterized by the relation divu s = 0. So the vector field u s of displacements in the transverse direction is solenoidal.
Expressions for such parameters of the elastic medium, as the Young module (elasticity coefficient) E and the coefficient of Poisson (deformation) σ [5], have the following form: ., Here We use the thesis that monochromatic wave in the medium is described by the Helmholtz equation where x ∈ R 3 , k 0 = ω/c is the wave number, ω is the frequency, c is the value of the wave propagation velocity in the medium, n(x) is the refraction index. We are looking for the solution to (3) in the form of a plane wave, assuming that the wave front is locally almost plane and changes weakly [6], Here the amplitude a(x) and wave vector k(x) = k 0 ∇ψ(x), k = (k 1 , k 2 , k 3 ), also changes weakly compared to wave length. Representing the field u(x) as a series and substituting (5) into (3), we obtain infinite system of equations ⟨∇ψ, ∇ψ⟩ = n 2 , 2⟨∇a 0 , ∇ψ⟩ + a 0 ∆ψ = 0, Function ψ is called eikonal, the first equation of the system (6) is eikonal equation. We use eikonal equation as equality connecting approximate values of longitudinal and transverse velocities. At first we state the problem, then we discuss approximation algorithm to realize the approach for numeric solution, and at last we demonstrate certain results of solving the problem based on real earthquake data from a number of seismic stations in a region of high seismicity.
2. An inverse kinematic problem of seismic with internal sources Now formulate the inverse kinematic problem of seismic [4], [7]. Let D be a half-space x 3 > 0 in R 3 with the boundary ∂D equal to the coordinate plane OXY , B ⊂ D be the bounded domain (the focal zone). The half-space D consists of elastic medium, where the waves with a velocity v(x) = 1/n(x), x = (x 1 , x 2 , x 3 ) are propagating. We consider the velocities v p of longitudinal and v s of transverse waves. Let τ (x 0 , y) be the travel time for the wave, propagating from the source at a point x 0 ∈ B to the detector at a point y ∈ ∂D. The function τ (x 0 , y) (= τ (y, x 0 )) is defined as follows, Here γ(x 0 , y) is the geodesic connecting the points x 0 , y, x(s) is the current point of the geodesic. The function τ (x 0 , y) satisfies to the eikonal equation with the condition on infinity The problem is to find v(x) for x ∈ B ⊂ D, if we know τ (x, x 0 ), which satisfy (8) for x ∈ ∂D, x 0 ∈ B. It is assumed usually in the inverse kinematic problem of seismic that the subset B belongs to the boundary ∂D, i.e. B ⊂ ∂D. In our case B is a 3D domain, B ⊂ D, so this problem is the problem of the inverse kinematic problem of seismic with internal sources. The problem statement is based on the following assumptions. The coordinates of the sources and the times of the earthquakes beginning are known.
Besides, the conditions for using the eikonal equation (8) are satisfied. At the first stage for the available set of values of the detected arrival times we construct the continuous travel time function (7) for the points at the focal zone B. Then on the base of the eikonal equation we determine the velocities v p and v s .

Numeric solution of IKPS with internal sources
From here and below we use traditional designations P, Q, . . . for the points of R 3 , and (x, y, z) for their coordinates. Moreover we use the designations conventional for the approximation theory.
At first consider the problem with only one detector. Let a detector, located at the point P 0 (x 0 , y 0 , 0) on the day surface, registers the incoming signals from certain points P i (x i , y i , z i ), i = 1, . . . , n, called sources of disturbances. The sources fill the spatial region inside the Earth (focal zone) in some way, generally speaking, chaotically. The times τ (P i , P 0 ) of wave arrival from the source P i to the detector P 0 are assumed to be known. Since the location of the detector P 0 is fixed we know the values t i of the function T (P ) at the points The eikonal equation determines the dependence of the gradient of the function T (P ) and the function v(P ) of the signal propagation velocity at the point P . Our problem is to determine the velocity distribution in the focal zone. If the coordinates of the sources had a regular location, for example, at the nodes of a three-dimensional rectangular grid, it would be easy enough to calculate the approximate values of the partial derivatives of the function T at the grid nodes, and, consequently, the values of the gradient.
However, in a real problem, the sources of disturbances are located in the hypocenters of earthquakes and, of course, there can be no question of any regularity. In addition, the data contain significant errors. In this case, difference methods do not give acceptable results. It is necessary to construct a function that not only approximates well the available values of t i = T (P i ), but also has a good approximation of derivatives. It is clear that local methods, for example, based on triangulation, are hardly acceptable here. Non-local spline approximation methods based on radial basis functions [8,9] seem to be "good" methods for such a problem.
Splines of this type can be obtained as a result of minimization of some functional J(f ) on the set of functions f from a certain space (condition of minimality of the norm) under the condition of interpolation of a given set of values or their smoothing [10,11]. As a result, the spline is represented in the form Here P i , i = 1, . . . , n, are basis spline nodes (points of multidimensional space), R 3 in our case; R(P, P ′ ) = ϕ(|P − P ′ |) is a radial basis function, ϕ is a univariate function, and |P − P ′ | is the Euclidean distance between points P and P ′ of the space; λ i ∈ R are the spline coefficients; π ∈ P, where P are given finite-dimensional space of functions in the considered multidimensional space, called the spline trend. Usually P is the space of algebraic polynomials of some degree.
Values of the coefficients λ i , i = 1, . . . , n, of interpolation spline s(P ) are determined from the interpolation conditions and conditions on the spline trend Choosing some basis u j ∈ P, j = 1, . . . , k, k = dim P, and representing the function π of the trend P in the form we obtain a linear system of equations Here B is symmetric (n × n)-matrix with the entries b ij = R(P i , P j ), C is (n × k)-matrix with the entries c ij = u j (P i ),λ = (λ 1 , . . . , λ n ) T andμ = (µ 1 , . . . , µ k ) T are vectors of unknown spline coefficients, andt = (t 1 , . . . , t n ) T are given function values. For uniqueness of finding spline coefficients we demand n ≥ k and rank C = k. This conditions are also sufficient, if the function R(P, P ′ ) conditionally positive definite ( [12]), so it plays a role of a reproducing kernel of some functional space. The condition for the minimum norm for a spline also allows us to formulate the problem of constructing a smoothing spline s α (P ), as a result of minimizing the functional where α > 0 is a smoothing parameter, which we can vary. Assuming that the system of equations (11) has the unique solution (n ≥ k and rank C = k), we obtain that the smoothing spline s α (P ) is also determined uniquely and its coefficients satisfies the system of equations Here I is the identity matrix. It is evident That for α → 0 the smoothing spline s α tends to interpolation one, and for α → ∞ as a limit we have a trend function from P, minimizing the  5 We consider the radial basis functions [12] of the form ϕ(r) = (−1) m/2+1 r m ln r, m even, that reduces to D m -splines and thin-plate splines [9,10]). Trend here is a set of three-variable polynomials of degree m − 1. The number m is called spline degree. As follows from [12], trend P can be extended by adding a set of linearly independent basis functions. The conditional positive definiteness of the chosen radial function is preserved, and it will also be the reproducing kernel of the space with an extended trend. Thus, restrictions on the dimension of the polynomial space for trend are removed.
A trend space determines the behavior of the spline outside the domains with data. Therefore, when we use, for example, the classical cubic thin plate splines the trend is a set of polynomials of the first degree. This means that outside the data domain, the spline will behave close to the plane. However, this behavior may not suit us if we know for some reason that the function being approximated, for example, is bell-shaped, i.e. we would like the trend to contain polynomials of degree 2 or higher. As we can see, in our case it is possible. In a number of problems, we have already tested the approximation by such splines with an extended trend in two-dimensional and three-dimensional cases [13][14][15], which showed quite good results.
Practical construction of an interpolation spline is reduced to solving the system of equations (11). Questions of the existence and uniqueness of splines of general form were studied in [16]. Note that restrictions on the initial data that ensure the existence and uniqueness of the splines under consideration are not burdensome and, in essence, reduce to the question of the existence in the set of data points of a set of k points for which the interpolation polynomial π(P ) exists and is unique. For example, for a trend consisting of polynomials of the first degree, for the existence and uniqueness of a three-variable spline, it is sufficient that among the n initial data points there are three points that do not lie on one straight line.
If the initial data are specified with an error, then usually interpolation splines do not give an acceptable result. In such cases, it is necessary to construct smoothing splines, i.e. it should be solved the system (13) instead of the system of equations (11). Changes in the matrix are only in the elements of the main diagonal of the first block.
Construction of both interpolation and smoothing splines is reduced to solving a system of linear equations with an almost dense symmetric matrix, the size of which is mainly determined by the number of interpolation points P i . For sufficiently large values of n, such a system may be ill-conditioned. Despite the fact that the matrix of the system (11) or (13) is symmetric, the Cholesky method cannot be used to solve it, since these matrices are not positive definite. We use Aasen method [17], in which the decomposition replaces the diagonal matrix with the tridiagonal one.
Note that although the complexities and difficulties of constructing interpolation and smoothing splines with radial basis functions are similar, the nonzero smoothing parameter α improves the conditioning of the system of equations. The smoothing parameter adjusts the level of data smoothing. The value α = 0 corresponds to the interpolation spline, i.e. s(P ) = s 0 (P ). With an increase in α, the behavior of the spline becomes smoother, oscillations decrease or even disappear, but at the same time, generally speaking, the deviation from the initial data increases, as already noted, in the limit as α → ∞ the spline tends to some element of the trend. Let us add that, generally speaking, the degree of smoothing can be adjusted at each point. For this, in the minimized functional (12), we need to add appropriate weights to each term, what does not lead to any complications in the calculations.

Numerical modeling of an elastic medium based on real data
Solving a problem of reconstruction of medium elastic characteristics we use our algorithm and a software package SEISMIC, perfectly performed on synthetic data.
The data for the proposed algorithm are the coordinates of earthquake hypocenters P i (x i , y i , z i ), i ∈ I, time values of happened earthquakes and registration times of the longitudinal and shear waves recorded by various situated at the surface of the Earth seismic detectors with known coordinates Q j (x j , y j , z j ), j ∈ J.
The differences τ (P i , Q j ) between the moments of earthquakes in the hypocenters P i and the moments of arrival of the waves to the detectors Q j are considered as discrete hodograph. So hodograph for every Q j as a function of P is given by its values on a finite chaotically located set of points P i in the focal zone B. Note, that the sets of such points P i , where we have the values of τ (P i , Q j ) are, generally speaking, different for different receivers Q j .
Earthquakes almost never occur simultaneously. In order to collect data, sufficient for executing our algorithm, we divide observation period into time intervals of weak seismic activity in the region before some strong shock and after it. So the time of the shock plays a role of common time boundary of this intervals. Such time intervals provide results for precise enough estimation of average state of the medium in the focal zone before the strong shock and after it. The study of seismic activities in dynamics can be carried out by analyzing of several such time intervals. These results can provide a basis for the current assessment and prediction of the behavior of the internal structures in the Earth.
Operating real data from one well known seismically active region of Kamchatka we outline corresponding focal zone and performed calculations using above algorithm. Our experiment concerns a strong earthquake with great ecological consequences occurred in the region of Karymsky Volcano and Akademia Nauk Volcano on January 1, 1996, further we call it Event. In particular, the former freshwater Karymsky Lake was filled with a high concentration salt solution. In order to analyze the effect of strong shocks on the state of medium in the depth we took the moment of this event to appoint the boundary of two time intervals (1.01.95-31.12.95) and (1.01.96-31.12.96) for data sampling. We choose the Cartesian coordinate system with the origin O, located at a point on the day surface (near the Karymsky Volcano, further we call it simply Volcano). The Oz axis is directed towards the center of the Earth. The axes Ox and Oy lie in the plane, tangent to the day surface at the point O. One tangent axis (Oy) passes through the epicenter of Event, located about 10 km to the south of the origin. The Ox axis is directed to the east. In this coordinate system, the line of interest (passing through two named volcanos near the Lake) defines the vertical section x = 0. To simplify demonstration of 3D results we show graphs just in this plane section, which have a form of rectangular.
The calculations were carried out on the data obtained from three stations, located outside the focal zone at a distance of about 100 km from the epicenter of the Event. So, applying the algorithm SEISMIC, based on the described above mathematical methods, we get 3D-model of the velocity characteristics from collected in the focal zone data, that give us decreet hodographs of P -waves and S-waves, registered by every station. That makes it possible to determine the functions of the velocity values v p (x, y, z) and v s (x, y, z) already at any point in the focal zone. Hence, it becomes possible to determine at each point main elasticity characteristics of medium such as the ratio Lame parameters λ and µ, (1), (2). The density value in medium structure is usually determined by dependence on the depth and some empirical information.
At great enough depths we can assume, without loss of generality, ρ = 3 (T /m 3 ) and calculate the elasticity coefficients in units of this value of ρ. Since the sources of recorded earthquakes are located at depths of up to 60 km, we define the focal zone in depth by the boundaries of 0 ≤ z ≤ 60 (km). So we are to compare the state of the medium before and after the Event, in the direction along the line, passing through one of the Volcano and the Lake. In the vertical section x = 0 of the cartesian coordinate system we show the graphs in the rectangular with the horizontal side −40 ≤ y ≤ 10 (km) and the vertical side 0 ≤ z ≤ 60 (km) ( figure 1 -figure 5).
The Volcano is located near the origin of the coordinate system and the Lake is located near the mark of −10 (km). Marks of depth (in km) are given on the vertical axis Oz. Now one can see and analyze the great difference between two parts of following figures, containing information about the medium before the Event (left pictures), and after the Event Before the Event, the number of recorded earthquakes for which the arrival times of longitudinal and transverse waves were recorded is about 100. After the Event, the number of the same relatively weak earthquakes is about 400. It remains unknown for us why this numbers greatly differ. But we can propose that evident great changes illustrated on the right graph may be associated just with the Event and became simply a consequence of the strong shock that caused further multiple increased number of seismic shocks.
The distribution of the dimensionless Poisson's ratio before and after the Event is shown on the pictures a) and b) in the figure 2.
Artifacts in the depth of 35 (km) may be an evidence of abnormal changes of the media state, and these changes distributes up to the Lake in the region of the mark from 5 to 10 (km) in a depth and from −10 to 10 (km) in lateral.
The Young module (in GPa, under the assumption that the density is uniform and equal to 3T /m 3 ) is shown on the pictures a) and b) in the figure 3.
We illustrate the distribution of Lame parameters λ and µ (in GPa and the same assumptions on the density) in figure 4 and figure 5 correspondingly. We calculate v p , v s and v p /v s independently for data from every of three taken seismic stations. All these functions are changed within reasonable limits. Therefore averaging over three stations of functions v p /v s , σ, E, λ and µ also have reasonable for analysis nature-consistent values. In the parts of focal zone, where hypocenters are distributed extremely rarely or outside the focal zone we get an extrapolation, which also can be reasonably analyzed. The approximation problem is reduced to solving a system of linear algebraic equations with a dense matrix. To solve the system, the reliable and fast methods with control of the condition number are applied. The eikonal equation gives an approximations for the velocities of longitudinal and transverse waves in the focal zone, since the approximation of the travel times is a differentiable function.

Conclusion
The results of the numerical implementation of the approach, carried out on the basis of the data of earthquakes between 1995 and 1997 are presented.
The presented mathematical model of waves propagation show the sensitivity to a strong earthquake on 01.01.1996 in the area of Kronotsky Bay. Velocity changes in the medium structures, both near the epicenter and at different distances from it, are significant both in depth and laterally, but do not fall outside the reasonable limits of the values of the kinematic and elasticity parameters. This sensitivity gives reason to believe that the considered model can be used for prognostic purposes. It seems the prospects for further research of the problems, arising from the obtained results, look promising.
In this work, the times of the onset of earthquakes and the coordinates of the hypocenters were considered as known, but it would be interesting to carry out the data refinement procedure based on the approaches described in Section 2 and in [14].
It seems the study of the functionals, proposed there, will help to solve the problem of extrapolation of the focal zone, including up to the day surface, in combination with a seismic tomography tools and a-priori data on wave velocities in the near-surface zone.
It looks very tempting to include a time variable in the model for the purpose of real-time monitoring of the state of the deep structures of the Earth, including for the purpose of predicting the possibility of strong earthquakes.