Seismic Inversion for the Calculation of Velocities Using the Generalized Inverse Linear Matrix, the Wave Equation, and a Non-Reflective-Boundaries Condition

1. Direct Model

The wave equation (1) is classified as a hyperbolic differential equation (Burden et al., 2016) in three dimensions:

where V is the velocity, which depends on spatial coordinates, and ϕ is the wave function. Considering some boundaries and initial conditions, such as:

where E represents the spatial coordinates, and ∂E is the E boundary; also f(E),g(E) are the known functions in the entire domain. In Cartesian coordinates (x, y, z) and considering integers (i, j, k, n), so that

with the soace and time discretized. Using Taylor series expansion with a fourth order approximation we obtain:

for the coordinates y and z we obtain analogous expressions. By replacing the equation system (7) in (1), solving 𝜙(𝑥_i, 𝑦_j, 𝑧_k, 𝑡_n+1) , and following the nomenclature proposed by Yee (1966), we obtain:

Equation (8) solves the wave equation if the values in the boundary and the initial conditions are known. The elements of the wave function in the time step n + 1 depend on the four nearest neighbors in a previous time step n and on the same component in a time step n - 1. The direct algorithm estimates the arrival time of a seismic wave (T _dir ), with propagation from the source j to the station k. The theoretical wave time (P) is obtained aggregating T _dir with the time of the seismic event (T_f), as can be seen in (9):

The theoretical calculation of wave time S assumes that the wave velocities S and P are proportional in the entire medium and behave like a Poison solid (10), with proportionality constant with . Thus, T_s ^teo = αT _dir - T _f and the difference of arrival times of waves S and P is:

As the station locations are fixed, T _dir depends on the velocity distribution and the hypocenter of the seismic event. Therefore T _p ^teo is the only variable that depends on the origin time.

2. Generalized linear inversion

In the inverse algorithm, we compared data measured on the surface with the results obtained from the simulation. Thus, we quantify the difference between the simulated values and the measured observations in the field. Then, it is possible to modify the initial values of the P wave velocity parameters in the field and the coordinates of the hypocenter. With these new initial values, we simulate the arrival times of the P wave for each one of the seismic events and compare the new results with the measured data on the surface, obtaining a new error difference, which will allow us to recalculate the new initial parameters for the next simulation.

This new methodology iterates until the difference between the simulated and the observed data on the surface is lower than the set threshold, in order to invert the velocity, the hypocenter and the origin time simultaneously and the origin time (Crosson, 1976; Ramadan, 2016). In the same way, the hypocenters inversion method was generalized in a simultaneous inversion of velocities and hypo-centers (Franco et al., 2006).

In the model, we consider N seismic sources with hypocenters given by x _j = longitude, y _j = latitude, z _j . = depth, and arrival times t _j , where j = {1,2,3...,N}. Additionally, we split the in L blocks space each block with velocity v. , where l = {1,2,3...,L}. These variables compose the model m⁰:

where 4N + L represents the number of elements.

To determine m ⁰ , we consider M stations (k = 1,2,...,M) located at (x _k ,y _k , z _k , t _k ) and the arrival time observed in the k-th position for the j-th source. A functional relationship T(m ⁰ ) allows us to reproduce the observed times according to the model m ⁰ . If this relation is linear, the sensitivity matrix defined as:

where A _j . is a MX4 matrix given by:

and C is a NMxL matrix, the i-th line of C is:

with i(j,k) = k + (j-l)M and δ _jki is.

To find Δm (16), we need to invert G, but G is not necessarily a square matrix. Therefore, we use the least square method to determine m, which consists of multiplying the matrix by its transpose and then we invert the result. According to this, (16) is transformed into:

where ΔT is the variation between the time observed (measured on the surface) and the time simulated:

The error between both times is calculated as follows:

and Δm is the change of the model, thus:

The matrix G ^T G may have determinant equal zero, in which case the matrix is not invertible; this is means that the system shown in (20) does not have a unique solution, and there is at least one null eigenvalue. As proved by Schwerdtfeger (Schwerdtfeger, 1960), there is a method to obtain the subspace solution to the highest likely dimension, known as Lanczos breakdown. This technique consists of decomposing G through matrices built with eigenvalues different from zero. The generalized inverse found with this method fits better the data on the model. However, the computational method cost makes it difficult when implementing real-world applications where matrices are bigger. Furthermore, there is difficulty to formally organize the eigenvalues, due to the data uncertainty because some of them may have values too close to zero to be significant. Therefore, we use a different approach to overcome this obstacle: a weight matrix Є ² W _m is added to the G ^T G as shown in (21):

where Є is the dampening factor. If e increases, the error curve presents fewer fluctuations because it depends on the iteration. This factor is found empirically, so we obtained a reasonable fit between the theoretical and observed data, and we also found that the variations of the model during each iteration are small. The W _m matrix in the most straightforward cases is the identity matrix. The initial model m ⁰ is proposed randomly but within a range of values that make sense in the geological model. With this model, the arrival time is simulated with the direct model and is obtained the difference with the reference time ΔT ⁰ .

It is possible to define the smoothing matrix (Stein & Wysession, s. f.) to modify the change rate of the velocity model. For a 3D geometry composed by a rectangle with several blocks I = L _x x l _y x L _z , the number of columns for D is L and has the following form:

where the sub-matrices D _x , D _y and D _z represent the first derivative in the directions x,y and z with the number of rows (L _x -1)L _y L _z ,L _x {l _y - 1)L _z and L _x L _y (L _z - 1) , respectively. The smoothing is applied to the velocity distribution, so it does not affect the hypocenter and origin time model, hence:

where I is the identity matrix with dimensions 4N x 4N

D _x , D _y , and D comprise the sub-matrices B _x , B _y , and B _z , defined as:

where,

Usually, B _x , B _y , and B_z, are equal, and the only difference is the separation of the columns between -1 and 1. For B _x there is no separation; for B _y the separation is L _x ; and for B _z separation is L _x L _y .

3.1 Disconnected wave equation

Hastings et. al. (Hastings et al., 1996) proposed a method based on (Berenger, 1994a) that consists of defining potential and decompose the wave equation according to the formulation of (Yee, 1966), is used for the electromagnetic fields.

The wave equation (1) can be split in connected equations, assuming that V is constant. Potential is defined as H = (Hx, Hy, Bz) so that:

also,

where 𝜙 = 𝜙 _x + 𝜙 _y + 𝜙 _z .

Following the scheme used by (Yee, 1966), the following nomenclature for the indeces is defined as:

and,

as suggested by Gonzalez (2012), the operators of centered difference are defined as:

According to the equations (30) to (32), the discretization of (27) in Taylor series with second order truncation error is:

and the discretization for the system shown in (28) is:

Equations (32) and (33) show the evolution rule, as well as the relation between the potential and the wave function with their neighbors.

Each component of the potential contributed to the full wave function in the corresponding direction (Figure 1 ).

Figure 1: Spatial discretization diagram for the wave equation solution: Circles represent potential, Hx (red), Hy (blue) and Hz (green), and the cross shows the wave function. (a) Relation of Hx with its neighbors (the components of the wave function). (b) Association of 𝜙 x , 𝜙 y, and 𝜙 z with its respective neighbors (the potential of each case). When the space split, and the disconnected wave function is in the semi-entire points and around the corresponding potential, which is going to help define its value.

Source: Authors.

3.2. Attenuation factors

Assuming that the space formed by a square and covered by a shell of determined width, we define an inside region called PML and a contact region called interface (Figure 2). So that one wave will expand freely in the interface and will be absorbed in the PML region.

So that the wave attenuation occurs, (Berenger, 1994b) proposes attenuation factors Q _s1 and Q _s2 in the (33) and (34), so the potential is redefined as:

and for the disconnected wave function:

The attenuation factors

where s = {x,y, z} and β = {i, j, k}, respectively, and

where P is the length of the PML region and r(β) is the depth from the interface to the outside edge as shown in Figure 2(a) and 2(b).

3.3. Complete solution

There is a slight change in the interface between the two regions, and the wave is slowly absorbed depending on the length P. Given the calculated limit in (38), when r → 0 then q → 0. So, from the (36) and (37) it is obtained that Q _s1 → 1 and Q _s2 → A _s thus the same (32) and (33) are reproduced. Then, we use the standard solution in the inside region by using finite differences (8) and redefining the indices following the diagram in (Kane Yee, 1966):

Figure 2: Wave space expansion. (a) The cube is cut in the z axis and forms a plane (known as a slice). A dark color represents the PML region. (b) Upper right edge of the cube.

Source: Authors.

In the PML region, we use the solution of the disconnected wave equation (35), as well as the following equation in the interface:

where I _d is the value that the index i takes in the right-side interface (Figure 2). Similarly, the equation of potential Hx is obtained in the left interface and for the rest of the potentials Hy and Hz. In other words, we have both the potential and the disconnected wave function in the PML region. In the interface we have only the potential and only the wave function in the inside (Figure 3). In most outer boundaries, we can use the Dirichlet or Neumann conditions (2).

Figure 3: Slices of the cube for the complete solution of the wave expansion with a non-reflective boundary: Circles represent H x (red) and H y (blue), while crosses represent 𝜙 x in (a), 𝜙 y in (b) and 0 (bold cross). Also, the arrows symbolize which neighbor feeds each term and r(β) is the depth of the PML region from the interface to the edge. (a) Discretization for 𝜙x and H x . (b) Discretization for 𝜙y and H y .

Source: Authors.

3.4. Scope and limitations

There are three possible limitations in the proposed model. First, the velocity blocks, hypocenters, and stations must be located; then it is possible to define the space and its limits based on those locations. Additionally, the space should be small for a short computational time; conversely, ample space contains the location of elements and stations of the model and allows achieving an equilibrium. Second, it is better to have many blocks for the tomography in order to be more detailed; however, there is a limit to the information available. Finally, since the stations are located on the surface and the hypocenter is distributed in the sub-surface, there is significant uncertainty of the velocity solution the larger the distance between the stations and the hypocenter.

4.1. Number of elements in the model

There is a limit for the number of elements considered in the model depending on the information available. According to the number of seismic events (N = 54), 403 observed times are determined (number of equations), including the arrival time of the P wave (TP) and the difference between the arrival times of the P and S waves (TSP). To have a better adjustment of the convergence of the generalized inversion corresponding to the observed information, the number of equations cannot exceed the number of unknown factors; otherwise, we obtain unexpected solutions. The unknown factors are the elements of the model; for every seismic event there are 4 (location x; y; z of the hypocenter and origin time), that is 4 x N = 216. They reduce the number of blocks to consider:

So, the maximum number of blocks is 187. If h_x, h _y , and h _z split up in L_x, L_y, and L _z intervals, then a matrix forms a L = L _x xL _y x L _z total blocks, so that L _x xL _y x L < 187.

4.2. Uncertainty of the velocity with depth

The last limitation is the uncertainty of velocity with depth. The trajectory between the hypocenter and the station is named ray according to the Fermat principle: the way taken by a wave expanding from one point to the other is such that the time used to cover it is minimum. The more significant the depth, the littler the ray density; so, the deepest cubes have a lower possibility of being crossed by a ray; furtherore, the uncertainty of velocity is more significant.