Computation and Modelling of Seismic Refraction Tomography for Anisotropic Medium Based on Pseudo Bending Method

Seismic refraction tomography is one of the imaging techniques in geophysical methods used to remodel the near-surface velocity layer structure of the Earth. In this study, we carried out a new computational approach and modelling of seismic refraction tomography using the pseudo-bending method. The true model of the near-surface is designed to be anisotropic medium which is having a low velocity anomaly distribution. This anomaly is constructed in such a way as to be similar to a model of liquid waste away spreading, which exhibits seismic velocities ranging from 1600 m/s to 1800 m/s. Based on our computations and numerical modelling results, it was found that the ray tracing path using pseudo-bending method displays an asymmetrical trajectory when the positions of the source and geophone are exchanged. Altering the shooting configuration from direct shoot (DS) to reversed shoot (RS) also reveals a significant difference in travel time values. The results of delay time tomography inversion, which represents the difference between travel times in the true model and the initial model using the SIRT method, indicate the presence of a low velocity anomaly that can be interpreted as the distribution of liquid waste.


Introduction
Seismic tomography is part of the geophysical method that plays an essential role in obtaining subsurface images and has become a reliable concept for estimating and imaging of geophysical parameters.In recent years, seismic tomography has been widely applied in the field of geophysics, especially for the hydrocarbon exploration, oil, and natural gas.In the civil and environmental engineering fields, particularly seismic refraction tomography, has been extensively utilized to image near-surface rock layer structures.This is primarily done to design the foundations of buildings, docks, airports, railways, highways, and more.Imaging subsurface layer structures can also be useful for imaging the distribution of liquid waste in the near-surface layers.Seismically, liquid waste is a type of fluid that can deflect the direction of seismic wave propagation and is characterized by low velocity anomalies.
Up to the present time, a number of researchers have been developed seismic refraction tomography methods.This method has developed for shallow subsurface investigation and more effective for site characterization compared to conventional seismic refraction [1].Imaging of subsurface can be conducted using seismic refraction method based on "width-band Fresnel" tomography simulation and utilized homogeneous velocity as an input model [2].Application of seismic refraction exploration for engineering and environmental geophysical can be implemented by employing the travel time tomography [3].In addition, the combination of seismic refraction travel time tomography with seismic reflection methods should be applicable for planning in urban areas and infrastructure development in densely populated regions [4].In general, seismic tomography has been successful applied for imaging the velocity structure both at global and local site investigations.For entire Japanese Islands including the Sea of Japan and the Pacific Ocean, seismic tomography conducted to image the shallow zone along the Sea of Japan by using air gun data [5].Travel time tomography can be also implemented to image the seismic velocity structure beneath Nabro and to analyse post-eruptive state by mapping subsurface fluid distributions [6].
In this study, we implement the MATLAB programming language to image the distribution of liquid waste in the near-surface soil layers.Over the past few decades, the author has prioritized the implementation of MATLAB programming in conducting geophysical computations.We have been focused and committed to using the MATLAB programming language for modelling vertical electrical sounding (VES) [7] and calculating the multimode Love wave dispersion curves [8].In this paper, our efforts focus to image the near-surface velocity structure based on the seismic refraction delay time tomography.This is accomplished by utilizing travel time data from a source to a geophone and employing the pseudo-bending ray tracing algorithm.The input velocity model data used in this research is assumed to be anisotropic, where the direction of seismic wave propagation differs for various rock layers.This approach closely approximates natural conditions where the Earth is composed of multiple rock layers with varying velocities that exhibit seismic anisotropy.
The fundamental issues addressed in the study include: • How to construct the actual algorithm for calculating the seismic refraction travel time when the medium is anisotropic.In this context, the calculation of travel time is based on simulating the propagation model of seismic refraction (ray tracing) from the source to the receiver (geophone).The approach used in ray tracing is the pseudo-bending method [9] and for anisotropic assumption, we then employ weakly anisotropic P-wave velocity equation [10].• The obtained travel time data will serve as input parameters in the inversion process.The subsequent issue to investigate is how to employ optimization using the Kaczmarz inversion method with a modified version of SIRT (Simultaneous Iterative Reconstruction Technique).• The inversion results obtained provide a representation (tomogram) of the conditions or structures of subsurface rock layers.Further, the identification of liquid waste distribution in the near-surface soil layers is carried out.The approach used involves analysing the inversion results (tomogram) and certainly confirming that the liquid waste will expose anomalous layer structures with low seismic velocity values.

Seismic Tomography Stages
The imaging of the seismic tomography is primarily based on the spatial distribution of seismic wave velocities using the inversion technique of delay time data, which is usually obtained from a set of seismograph stations distributed across the Earth's surface.The first step in seismic tomography is to formulate an iterative linear inversion by first defining the initial model and parameters, and performing seismic ray tracing through the model using specific source-receiver pairs.Seismic tomography studies can generally be classified into three main stages [11] namely: 1) defining an initial model and determining its parameters, 2) shooting a set of rays through a model to generate a collection of theoretical travel time data (ray tracing).This stage also involves determining the properties of ray propagation through the medium, from the source (s) to the geophone or receiver (r), known as the forward modelling step, 3) implementing an optimization process with aims to refine the model, minimizing the discrepancies or mismatches between observed travel time and theoretical travel time.This stage constitutes the core of tomography (backward modelling or inversion step).

Seismic Tomography Equation
Seismic tomography equations are based on a system of linear equations that can be expressed seismically as [12]: Here, ΔT represents the vector of travel time residuals, A is the Jacobian matrix which is the matrix of partial derivatives of the changes in travel time for a block with respect to the magnitude of velocity perturbation in the corresponding block (Aij=Ti  Vj), and ΔV is the vector of velocity model perturbation values.From this equation, it is evident that to determine ΔV, it is necessary to first calculate the travel time residuals and also construct the matrix of partial derivatives of travel time with respect to velocity perturbations.If there are M observations and N velocity blocks in the model, then the matrix A will have dimensions of MN.

Ray Tracing Employing the Pseudo-bending Method
Based on the Fermat's principle, the travel time of a ray path with respect to the change in path length is stationary (∂T/∂L = 0).In order for this stationary state to be achieved, the travel time along the ray path is minimized.This minimization equation is used to calculate perturbations by exploiting the characteristics of ray path curvature.The solution to the ray equation is carried out through the pseudo bending method [13].The travel time (T) along the ray path is expressed as a line integral between two end points, that is: where r and L is the path length parameter, and V is seismic velocity parameter.
The computation of travel time is carried out by using the numerical summation along the segments of the ray's trajectory.The travel time equation can be reformulated employing the trapezoid rule in the following scientific manner: where n is the number of points along the trajectory, Xk represents the position vector at the k-th point, and Vk denotes the magnitude of velocity at the k-th point.If travel time is directly minimized by simultaneously perturbing each point along the ray's path, this minimization process will involve solving a large number of nonlinear equations.As a result, an alternative approach is undertaken by employing a three-point perturbation scheme, as illustrated in figure 1.In addition, the ray equation can be written as follow:

Figure 1 Three-point perturbation scheme
where r is the position vector along the ray [14].The second term on the right-hand side of equation ( 4) represents the component parallel to the gradient of velocity with respect to the ray's path.This equation states that the normal component of the velocity gradient with respect to the ray vector consistently opposes the curvature of the ray's path.Since the local direction of the ray at X'k is defined by the line connecting the two endpoints of the path segment, Xk-1 and Xk+1, the normal component of the velocity gradient with respect to this direction indicates the direction of curvature.This direction corresponds to the factual location of X'k that conforms to equation (2).The formulation of this direction is provided by The second term in equation ( 5) constitutes the parallel component of the velocity gradient in relation to the ray's direction, or equivalently, with the normalized unit vector n = n' / |n'|.
The estimation of velocity at the new point X'k is necessary due to the unknown in the previous perturbation values.If a Taylor expansion is employed for the velocity at the mid-point (Vmid), the velocity at the new point V'k can be approximated through the expression: where (V)mid denotes the velocity gradient at middle point Xmid.The actual number of perturbation Rc along the direction n is derived through the minimization of the equation ( 6) along the segments connecting the three points Xk-1, X'k, and Xk+1.Consequently, the resulted accumulation perturbation can be expressed as: 1/2 (7) where

Seismic Anisotropic Equation
Seismic ray tracing in anisotropic media is crucial for exploring the detailed structures of the Earth's crust [15].The seismic anisotropy signifies a variation in seismic velocity, where the velocity is contingent upon the inherent elastic characteristics of the medium, relative to the direction in which the velocity is measured.Seismic anisotropy is explicitly defined as the dependency of velocity on angle [10].
In general, the discourse surrounding anisotropy primarily focuses on scenarios involving the propagation of waves, particularly in cases where the anisotropy present within the rock is weak.The seismic velocity of P-wave for weak anisotropy is formulated as [16] V P (θ)=α 0 (1+δ sin 2 θ cos 2 θ +ε sin 4 θ ) (8) where δ, ε, and γ are the anisotropy parameters (Thomsen parameters),  is the phase angle (0 0 is measured from the vertical direction).The 0 indicates the P-wave velocity at =0.

Tomography Inversion Techniques
The Kaczmarz algorithm is one of the iterative inversion methods that can be employed for solving the large-scale systems of linear equations [12].This algorithm starts with an initial solution x(0), and subsequently moves towards a solution x(1) by projecting the initial solution onto the hyperplane defined by the first equation in the system of equations.Further, the solution x(1) is projected onto the hyperplane defined by the second equation.The process continues until the solution is iteratively projected onto hyperplanes defined by all equations, converging towards the solution of the system.To implement this algorithm, a formulation is required to compute the projection of a vector onto the hyperplane defined by equation i.Let Ai be the i-th row of the matrix A. The projection of x (i) onto the hyperplane Ai+1 x = bi+1 is given by: A modification to the Kaczmarz algorithm leads to the Simultaneous Iterative Reconstruction Technique (SIRT), which practically provides improved images.In the SIRT algorithm, model updates are computed for cell j along each ray i that passes through cell j.The update computation for cell j is is carried out for all m rays, and then the average of all updates is computed before applying the update to cell j.The formulation of the SIRT update can be expressed as follows [12]: not pased ray i+1 (10)

Method
This article presents a numerical procedure based on the quantitative method by implementing a MATLAB programming language.In general, the computational and numerical modelling approach can be absolutely carried out if the algorithm of the actual calculation from an empirical system equation is definitely established.Based on the theoretical study of tomography and lists of equation presented in the literature review section, the actual computation and modelling of seismic refraction tomography can be briefly summarized as follow: • Defining an initial model and setting up its parameters, including the number of blocks in the horizontal axis (nx) and vertical axis (ny), dimension of the blocks in the horizontal axis (blx) and vertical axis (bly), the number of sources and receivers (ns and nr), the positions of source (xs,ys) and receiver (xr,yr), the value of bounding (r_bound), number of segments (nsg), the value of polar anisotropy parameters namely near-vertical P-wave (delta) and horizontal Pwave (epsilon).• Shooting a set of rays through a model (ray tracing) to generate a collection of theoretical travel using pseudo-bending method.At this stage, equations 5-7 will be employed to illustrate the simulation of seismic wave propagation or ray tracing from the source to the receiver.Each time the ray moves between blocks, the computation of the segment's length within the traversed block will be executed.The lengths of these segments within each block are recorded in a kernel matrix.The multiplication of the kernel matrix by the velocity vector will yield the travel time values from a source-i to receiver-j.For an anisotropic medium, the weak anisotropic velocity of P-wave [16] will be employed (equation 8).• The travel time data acquired in stage 2 is employed as a parameter input in running the inversion process.The subsequent issue to investigate is how to employ optimization using the Kaczmarz inversion method with a modified version of SIRT (Simultaneous Iterative Reconstruction Technique).• The inversion results obtained from SIRT provide a representation (tomogram) of the conditions or structures of subsurface rock layers.Further, the identification of liquid waste distribution in the near-surface soil layers is carried out.The approach used involves analyzing the inversion results (tomogram) and certainly confirming that the liquid waste will expose anomalous layer structures with low seismic velocity values (low velocity).
• The entirety of the computational process (numerical modelling) to be carried out in this study is implemented by using the MATLAB programming language.

Travel time computation of anisotropic seismic refraction using pseudo-bending method
In order to ascertain ray tracing and travel time within a seismically anisotropic medium, the seismic wave velocity in each block intersected by the ray is determined through the utilization of the P-wave velocity equation ( 8).The values of the anisotropy parameters are obtained from various literature sources.For the case of weak anisotropy, the parameters are δ=0.05 and ε=0.20 [10,16].The magnitude of the angle θ is measured based on the direction of ray curvature on the respective block.Consequently, if a ray passing through a block from a different direction, the seismic wave velocity will yield different values depending on the incoming ray direction on that block.The depiction of the correlation between P-wave velocity and both the angle θ, as well as the anisotropy parameters (ε and δ), is illustrated in figure 2. Based on the illustration in figure 2 and equation 8, in the context of isotropic media (with δ=0.0 and ε=0.0), the P-wave velocity (VP) in an isotropic medium is independent of the measurement direction.Conversely, in an anisotropic medium, the wave propagation velocity is highly dependent on the angle theta (θ) or the measurement direction.In this context, waves arriving from different directions will yield different travel time calculations.In the perspective of anisotropic media, for the same offset distance, travel time calculations with a direct shoot (DS) configuration (receiver or geophone located to the right of the source) will produce different values compared to a reversed shoot (RS) configuration (receiver or geophone located to the left of the source).This emphasizes the crucial consideration needed in achieving more precise calculations of seismic wave travel times, particularly when employing an approach of the anisotropic nature of the medium.

Seismic Refraction Anisotropic Tomography: Numerical Test
In its implementation, seismic refraction is a geophysical method used to investigate the subsurface characteristics of the Earth.In this study, seismic refraction anisotropic tomography is used to identify the subsurface distribution of liquid waste near the surface.Liquid waste, in its physical state, behaves as a fluid.Seismically, it is characterized by having a relatively low seismic velocity, ranging from about 1500 to 1800 meters per second (m/s).Figure 3 provides a visual representation of a twodimensional (2D) cross-section of the Earth, where in the near-surface region, there is a distribution of low anomaly velocity that can be interpreted geologically as the spread of liquid waste.A legend displayed in the form of a vertical color bar on the right side of the image indicates the seismic velocity intervals in m/s.Bright color indicates lower seismic velocity, while dark color indicates higher seismic velocity.Moving forward in the tomographic process, the next step involves the application of a technique called ray tracing pseudo bending to the true model to obtain travel time data (observed travel time) of seismic refraction waves from the source to the receiver or geophone.This step is crucial for obtaining data related to the travel time of seismic waves, often referred to as observed travel time.Figure 4 shows the ray tracing on the true model with a reversed shoot (RS) configuration.This configuration is chosen based on its suitability for this particular study.Thereafter, the forward modelling process is carried out by creating a 2D cross-sectional initial model of the Earth, as illustrated in figure 5.This model serves as the starting point for the tomographic analysis.It is a representation of the subsurface structure before any anomalies or features are accounted for.Ray tracing using pseudo-bending method is applied to the initial model to obtain travel time data (calculated time) of seismic refraction waves from the source to the geophone.The pattern of anisotropic seismic refraction ray tracing on the initial model is depicted in figure 6.This pattern reflects how seismic waves propagate and interact with the subsurface structures in the model.The input data for the source positions, geophone positions, and shooting configuration remain consistent with the true model.The distance and depth are also kept the same.However, in the initial model, it is assumed that there are no velocity anomalies indicating the presence of liquid waste spread.The Earth cross-section is assumed based on theoretical data, where the Earth is composed of several layers of rock.The deeper the layer, the greater the velocity of the seismic waves.The curves of travel time versus offset which illustrate the seismic refraction travel time data from the source to the geophone for both the true and initial models can be presented in figure 7.These curves represent the relationship between travel time and offset (distance between the source and the geophone).The curve of TX [S1] provides valuable information about the travel time characteristics of the first wave source.Similarly, the curves of TX [S2] and TX [S3] respectively represent the travel time versus offset curves for the second and third wave sources.Each of these curves offers insights into the travel time behavior of their corresponding wave sources.The TX curves encapsulate crucial information for the interpretation of seismic refraction data.They offer a detailed record of how travel times are influenced by spatial separation, which is instrumental in deducing subsurface properties and detecting anomalies, such as the presence of liquid waste.Analyzing these curves aids in the refinement of the subsurface model and enhances our understanding of the geological structures at play.Furthermore, travel time data obtained on the previously stage are then processed by applying an inversion method (inverse modelling process).This process constitutes the primary stage in tomographic modelling.The inversion technique employed is based on the modified version of the Kaczmarz's algorithm known as simultaneous iterative reconstruction technique (SIRT).This inversion technique will be utilized to determine the number of unknown variables, which is equal to the number of blocks (nxny = 200), given the available number of equations, which is equal to the number of rays (nsnr = 30).Specific cases like this are known as under-determined cases.Mathematically, the solution obtained will not have an exact value.
In solving seismic tomography equations, under-determined cases can be addressed by incorporating additional geophysical or geological data.In this study, additional data used in the inversion process are theoretical data based on the assumption that within the same layer, all blocks along the horizontal direction in that layer have the same seismic velocity.This supplementary data is gathered in a matrix known as the horizontal smoothing matrix.This method can be utilized as a priori information that can provide better resolution in the tomographic image.If the available data is fewer than the model parameters being sought, the presence of horizontal smoothing in the inversion process allows for the transformation of under-determined cases into over-determined ones.
The execution results of the inversion program using the SIRT method are illustrated by the tomograms (contours of 2D Earth cross-section) in figures 8, 9, and 10.The figure 8 displays the SIRT inversion outcome for the 50th iteration, while figures 9 and 10 represent the SIRT inversion results for the 150th and 400th iterations.Qualitatively, it is evident that the inversions in figures 8, 9, and 10     The contour plots of the 2D Earth cross-section from the inversion results exhibit layered structures that resemble the true model (refer to Figure 3).Theoretically, it is understood that the obtained inversion results will not perfectly match the true model or the velocity values obtained from the inversion results differ from the actual velocity values in the true model.This implies that there is a relative error percentage or deviation () in the inversion results compared to the true model.Numerically, the deviation () in the inversion results can be calculated using the equation: where N is number of block (nxny).Deviation image of the inversion results for the 50th, 150th, and 400th iterations is performed in figure 11.The calculation results using the equation ( 12) produce respectively deviations of 3.802%, 3.418%, and 3.564% for the 50th, 150th, and 400th iterations.

Conclusion
The simulation of seismic refraction ray tracing using MATLAB programming language with the pseudo-bending method can be utilized for calculating travel times.The computation of travel time from the source (s) to the receiver (r) with an anisotropic medium approach becomes an input parameter in the main stage of seismic tomography.In the perspective of anisotropic media, for the same offset distance, travel time calculations with a direct shoot (DS) configuration produce significantly the different values compared to a reversed shoot (RS) configuration.The implementation of the SIRT inversion technique in anisotropic seismic refraction tomography is capable of revealing images or contours of subsurface layer velocities (near surface).The 2D Earth cross-sectional images obtained from the inversion results show the presence of a distribution of low velocity anomalies at depths ranging from 20 m to 100 m.These anomalies are interpreted as indicative of the presence and distribution of liquid waste.The velocities associated with these anomalies fall within the range of 1600 m/s to 1800 m/s.Furthermore, the assessment of the inversion process includes an evaluation of relative errors.Specifically, for the 50th, 150th, and 400th iterations, the deviations of the inversion results from the true model are quantified.These evaluations demonstrate a remarkable level of accuracy, with relative error percentages of 3.802%, 3.418%, and 3.564%, respectively.This suggests a high degree of fidelity in the reconstructed subsurface velocities obtained through the seismic tomographic process.

Figure 2 .
Figure 2. The correlation P-wave velocity with the angle θ.Graphs of type A, B, and C represent the relationship between P-wave velocity values and the angle theta (θ).These values are computed through the utilization of equation (8), employing distinct combinations of Thomsen parameters: A(δ=0.05 and ε=0.20),B(δ=0.15 and ε=0.15), and C(δ=0.50 and ε=0.10).The initial P-wave velocity, VP(θ=0), is established at 2100 m/s.Based on the illustration in figure2and equation 8, in the context of isotropic media (with δ=0.0 and ε=0.0), the P-wave velocity (VP) in an isotropic medium is independent of the measurement direction.Conversely, in an anisotropic medium, the wave propagation velocity is highly dependent on the angle theta (θ) or the measurement direction.In this context, waves arriving from different directions will yield different travel time calculations.In the perspective of anisotropic media, for the same offset distance, travel time calculations with a direct shoot (DS) configuration (receiver or geophone located to the right of the source) will produce different values compared to a reversed shoot (RS) configuration (receiver or geophone located to the left of the source).This emphasizes the crucial consideration needed in achieving more precise calculations of seismic wave travel times, particularly when employing an approach of the anisotropic nature of the medium.

Figure 4 .
Figure 4. Ray tracing on the true model with a reversed shoot (RS) configuration.

Figure 6 .
Figure 6.Ray tracing on the initial model with a reversed shoot (RS) configuration.

10th
Asian Physics Symposium (APS 2023) Journal of Physics: Conference Series 2734 (2024) 012049 IOP Publishing doi:10.1088/1742-6596/2734/1/01204910 clearly delineate the dispersion of liquid waste at depths ranging from 20m to 100m.This is indicated by the presence of a low velocity anomaly with velocity intervals between 1600 m/s and 1800 m/s.

Figure 8 .
Figure 8. Image of the inversion result for 50 th iteration.

Figure 9 .
Figure 9. Image of the inversion result for 150 th iteration.

Figure 10 .
Figure 10.Image of the inversion result for 400 th iteration.