Joint tomographic inversion of first‐arrival and reflection traveltimes for recovering 2‐D seismic velocity structure with an irregular free surface

Irregular surface flattening, which is based on a boundary conforming grid and the transformation between curvilinear and Cartesian coordinate systems, is a mathematical method that can elegantly handle irregular surfaces, but has been limited to obtaining first arrivals only. By combining a multistage scheme with the fast‐sweeping method (FSM, the method to obtain first‐arrival traveltime in curvilinear coordinates), the reflected waves from a crustal interface can be traced in a topographic model, in which the reflected wave‐front is obtained by reinitializing traveltimes in the interface for upwind branches. A local triangulation is applied to make a connection between velocity and interface nodes. Then a joint inversion of first‐arrival and reflection traveltimes for imaging seismic velocity structures in complex terrains is presented. Numerical examples all perform well with different seismic velocity models. The increasing topographic complexity and even use of a high curvature reflector in these models demonstrate the reliability, accuracy and robustness of the new working scheme; checkerboard testing illustrates the method's high resolution. Noise tolerance testing indicates the method's ability to yield practical traveltime tomography. Further development of the multistage scheme will allow other later arrivals to be traced and used in the traveltime inversion.

Most modern tomographic techniques addressing an irregular surface are based on unstructured grids, which is computationally inefficient in model parameterization and traveltime calculation (Kimmel and Sethian, 1998;Sethian, 1999;Sethian and Vladimirsky, 2000;Rawlinson and Sambridge, 2004a, b;Qian JL et al. 2007a, b;Kao CY et al., 2008;Lelièvre et al., 2011). Unlike these techniques, the structured grid-based schemes usually handle ir-regular surfaces by using model expansion (Vidale, 1988;Reshef, 1991;Hole, 1992;Zhang Z, 2014a, b, 2015) or irregular surface flattening (Haines, 1988;Lan HQ andZhang ZJ, 2011a, b, 2013a, b;Lan HQ et al., 2012). The former scheme usually employs a stair-step approximation of the irregular surface or a flat low-velocity layer covering the surface, which may cause accuracy loss and even distorted images. The latter scheme takes advantage of a boundary conforming grid and a transformation from Cartesian to curvilinear coordinates, allowing the irregular topographic model to be mathematically flattened. In this method, first-arrival traveltimes can be obtained by solving the topography-dependent eikonal equation (hereafter TDEE) (Lan HQ and Zhang ZJ, 2013a).
On the other hand, reflected waves revealed by deep seismic reflection profiling may be numerous and even have large amplitude if the seismic impedance is great; they can be very useful in probing the crust, and are therefore potentially a valuable resource for seismic imaging. Reflected waves also contain information of both the velocity of the elastic medium and the seismic reflector and hence can be used to make joint inversion for velocity structure and interface depths (Li SL and Mooney, 1998;Knapp et al., 2004;Scarascia and Cassinis, 1997;Zhang ZJ and Klemperer, 2010), which is a natural complement of migration imaging (Rawlinson and Goleby, 2012). Having a larger set of seismic first-arrival and reflection data and the right tool to carry out joint inversion of all of them is a task worth addressing. A number of grid-based schemes have been proposed for tracking reflected waves, such as the two-way approach and the multistage scheme, but most of them are based on the assumption of a flat surface and therefore cannot be used with irregular topography (Benamou, 1996;Symes and Qian JL, 2003;Podvin and Lecomte, 1991;Riahi and Juhlin, 1994).
The aim of this paper is to obtain accurate seismic velocity models for irregular free surfaces, by using a mathematical surface flattening strategy and a multistage working scheme and by taking advantage of first-arrival and reflection traveltime data provided by deep seismic experiments. To obtain first-arrival traveltimes we use the fast-sweeping method (FSM) to solve TDEE formulated in curvilinear coordinates. To track reflected waves on a crustal interface we use a multistage scheme via first-arrival wavefront solved by FSM. We then present a joint inversion method of first-arrival and reflection traveltimes to determine velocity models for an irregular surface. We have arranged this paper as follows: first, we briefly describe the irregular surface flattening scheme, including the boundary conforming grid, the topography-dependent eikonal equation, and the method for computation of first-arrival traveltimes and ray paths. Next, we introduce the method for computation of reflection traveltimes and ray paths, and then the joint inversion of first-arrival and reflection traveltimes (JI-FRRT), including the inversion strategy. Finally, we present several numerical examples to illustrate the performance of JI-FRRT, which demonstrate the potential applications of this methodology to reveal the internal structure of the earth through modeling seismic velocity.

Topography-Flattening Scheme and Topography-Dependent Eikonal Equation in 2-D
The topography-flattening is achieved by transforming the Cartesian coordinate system to a curvilinear coordinate system, using grids conforming to the irregular surface to describe the model (Figure 1). Such a grid is called a "boundary-conforming grid" (Thompson et al., 1985;Hvid, 1994), and has been used by many researchers (Fornberg, 1988;Zhang W and Chen XF, 2006;Appelo and Petersson, 2009;Lan HQ andZhang ZJ, 2011a, b, 2013a, b). Under such a transformation, the curvilinear coordinates (q, r) in the mathematical space are mapped onto the Cartesian coordinates (x, z) in the physical space. Hence, a topographic physical model (x, z) is converted to a flattened mathematical model (q, r) (in the mathematical space) (Figure 1). Neighboring elements in the physical space are also adjacent in the mathematical space. Velocity nodes of the models are given at the curvilinear coordinates, but transformed and expressed in the physical space (x, z). Note that the irregular surface should be second-order derivable to insure an efficient computation.
Using the topography-flattening method, the classical eikonal equation expressed in Cartesian coordinates becomes a new generalized equation in the curvilinear coordinate system, which is topography-dependent (Lan HQ and Zhang ZJ, 2013a, b): with .
The "coefficients" A, B and C are topography-dependent. T is traveltime, s is slowness, J is the Jacobian that is written as and denotes , and similar in other cases.

Coordinate
Free surface The Lax-Friedrichs sweeping scheme, a fast-sweeping method (hereafter FSM), is applied to solve equation (1) and to obtain the first-arrival times. We determine the ray paths following the steepest traveltime gradient in the curvilinear system from the receiver to the source (Ma T and Zhang ZJ, 2014). Traveltimes at receivers and shots are calculated by interpolation to trace their positions exactly.

Multistage Scheme for Reflection Tracking
Traditional grid-based eikonal schemes for computing traveltimes are usually limited to obtaining first arrivals only (Rawlinson et al., 2004a), not excluding the topography-dependent eikonal equation. A multistage scheme is an efficient method for tracking later arrivals in layered media, and has been successfully used by some researchers in combination, respectively, with the fastmarching method (Rawlinson and Sambridge, 2014a, b) or the shortest-path method (Bai et al., , 2010Huang GJ et al., 2012). The former has not considered an irregular surface; the latter can handle the irregular problem but is time-consuming. In this paper, we try to combine the multistage scheme with a topography-dependent eikonal equation to trace reflected phase exactly in an irregular model. The two principal difficulties involved in introducing a multistage scheme is (1) how to combine it with first-arrivals' numerical solution processing (FSM); (2) how to represent accurately the reflecting interface, which generally varies with depth and does not conform to the regular velocity grid.
Faithful representation of the reflection interface was solved by introducing a set of nodes that irregularly lie on the interface and are independent from the velocity grid nodes. A connectivity between the irregular interface nodes and the regular grid nodes was constructed by local triangulation in the neighborhood of the interface; traveltimes at the interface nodes were attained by a linear interpolation scheme of the triangle with two points in the velocity grid ( Figure 2b). We adopted the multistage method, in combination with the FSM, for the numerical solution of the TDEE (called multistage FSM), to calculate reflected phases in layered media with an irregular surface. The multistage FSM scheme involves four stages ( Figure 2): the first stage initializes FSM at the source and tracks the incident wave front to all points on the reflecting interface ( Figure 2a); the second stage records traveltimes of first arrivals along the sampled interface, applying the local triangle linear-interpolation to calculate traveltimes on interface nodes from those on velocity nodes ( Figure 2b); the third stage tracks the reflected wave front by reinitializing FSM from the sampled interface ( Figure 2c); the fourth stage tracks ray paths along the negative gradient direction of the obtained traveltime field, from receivers to the reflector and then to sources (Figure 2d).
Since the multistage FSM is established by initializing the wave front from interfaces, later arrivals that propagate in multiple layered media can be correctly computed by following the propagating step. In this paper, we solve the topography case with an irregular surface on the top simply by extending the multistage FSM to 2-D curvilinear coordinates with a single reflector.

Back-Projection Algorithm δs
The inversion is carried out by a back-projection algorithm (Hole, 1992), in which the traveltime residues are uniformly projected along ray paths. Thereafter, the slowness perturbation in each grid cell can be written as: δt k here K is the number of rays passing through any of the neighboring cells, and l k are respectively the traveltime residue and the Interface nodes (black circles) are independent from the velocity grid nodes (gray diamonds), and traveltimes at interface nodes are obtained from those at velocity nodes by a local triangle linear-interpolation; (c) Reflected wavefront is tracked to receivers on the surface; (d) Ray paths are tracked along the negative gradient direction of the obtained traveltime field, from receivers to the reflector and then to sources. Source (S) and receivers (R) are denoted by star and triangles, respectively. total length of the k-th ray. The slowness perturbation is calculated in the ray tracing process, which belongs to the forward process. As the model is parameterized at grid points, the slowness perturbation at each grid point is found by taking the average of the values obtained from the cells surrounding this grid point. The inversion grid is a larger grid resampled from the grid used in the forward calculation; a smoothing factor is introduced to spread the slowness perturbation over a wider area and decrease local velocity anomalies. An appropriate choice of the resample and the smoothing factor ensures that the velocity model will be well updated. The velocity model is updated by iteration until the traveltime residue is sufficiently small and does not change significantly with the next iteration.

Joint Inversion Strategy
The use of first-arrival waves together with reflected waves identified from a deep seismic sounding experiment requires joint inversion strategies. For the sake of simplicity, in this article we assume a specified reflector to take care of exploring the seismic velocity structure of the medium; in other words, in the inversion process we set the reflection interface and try to obtain an accurate velocity image.
In the inversion process, a back-projection algorithm is applied, respectively, to the first-arrival and the reflection phases, and the slowness perturbations are calculated along each phases' paths, respectively. We then sum up the two classes of perturbations in each cell, to update the velocity (slowness) model. The use of the total slowness perturbation to update the velocity model at each iterative step allows us to obtain the final velocity model. Figure 3 shows a flow-chart that describes the JI-FRRT process.
In order to avoid large perturbations and over fitting of noise when applying the inversion routine in the computational domain, we first take the smoothing factor in the horizontal direction equaling to n x (the number of nodes in the horizontal direction), to keep the velocity fixed in this direction and only updated in the vertical direction. After an applicable 1-D model is obtained, we adjust the smoothing operator in both directions to update the velocity and finally obtain a 2-D model.

Numerical Examples
Below we develop several numerical examples in order to verify the performance of JI-FRRT. First, we consider three models having increasing topographic complexity, and an undulated reflector to prove the accuracy of the traveltime joint inversion procedure. Second, we analyze the quality of the results through the classical checkerboard test. Third, we also check the efficiency of the method assuming a model with irregular topography and an undulated reflector, together with an intracrustal high-velocity body. Lastly, we test the robustness of our working scheme by checking its noise tolerance. In all tests, the first-arrival and reflection traveltimes at each receiver, i.e., the input data for further seismic modeling, are calculated by forward process of JI-FRRT, by using theoretical models those to be inverted.

Different Topographic Complexities
We consider three topographic models labeled as Models 1, 2 and 3. Model 1 has a trapezoidal hill located in the middle of the mod- Model 3 is the one that has the most complex topographic surface, since it combines hills and valleys irregularly (Figure 4c). The topographic gradients were designed to increase gradually from Model 1 to Model 3 with the aim of testing the reliability and accuracy of JRFTT under different topographic complexities. We made several tomographic inversions using different grids, and found that the optimal grid spacing was 0.5 km in both vertical and horizontal directions. Figure 4 shows the model mesh 100 × 51 (number of grid nodes in the horizontal and vertical directions) displayed with a density reduction factor of 3, which represents the boundary-conforming grid for the three models, with a grid spacing of about 1.5 km. As before, we fired 10 shots that were recorded by 46 receivers evenly arranged on the irregular topographic surface. In order to include a non-flat reflective surface, thereby increasing the complexity of the initial model and bringing us closer to a real case, we introduced an undulated interface in the velocity field lying at the depth range 21-33 km. Then we carried out several inversions with different re-gridding and moving average factors, and finally found the optimal factors to be, respectively, (4, 2) and (5, 3) for all three tests.
We pick 460 first-arrival and 460 reflection traveltimes to image the model, which is parameterized by grids of size (101 × 51). Figures. 5a, b, and c show the tomography results for Models 1, 2, and 3, respectively. In all plots the continuous line denotes the wavy reflection interface. We can see the theoretical velocity models (upper panels), the initial models (middle panels), and the inverted models (lower panels) obtained by JI-FRRT. The results reconstruct the theoretical models with very little error, even with increasing topography complexities and starting from an initial model that clearly deviates from the target model.
After obtaining the inversion solutions, we are interested in estimating the quality (accuracy) of the results. Thus we evaluate the results based on the following criteria. First, since our working scheme updates the velocity along the ray paths, we need to display the ray geometry and the ray coverage density, i.e. the number of rays that cross each grid cell. This information is given in Figure 6, which reveals the illumination of the explored medium.  Second, it is necessary to estimate the error made in the calculation of the solution. With respect to the regions with denser ray coverage (>5 rays passing through each cell), Figure 7 shows the relative velocity error regarding the final inverted models (lower panels) and that inherent to the initial models (upper panels), both compared to the theoretical models. As can be seen, the former are much smaller than the latter ones, the initial error being larger than 10% and the final error smaller than 3%, especially for the middle region above the reflection interface. Third, RMS traveltime residues versus number of iterations demonstrate the rapid convergence of the inversion for each of the three numerical examples (Figure 8). After about five iterations, the residues gradually stabilize and converge to approximately 0.004 s in all cases. Both the convergence of traveltime residues and the convergence of velocity errors to very small values demonstrate the consistency of JI-FRRT.
In order to evaluate the vertical and horizontal resolution of our inversion scheme, we also perform the checkerboard test for Models 1, 2, and 3. The models are divided artificially into zones of relatively high and low velocity at a scale of 7.5 km × 7.5 km, so that velocity anomalies described by the formula 0.3 × sin(x) × sin(z) are added alternatively to a common reference velocity for all models (upper panels in Figure 9). The synthetic traveltime is obtained from the checkerboard model; the common reference velocity is taken from the initial model. In the inversion process, the rays are unevenly distributed in the model due to the undulated reflector, so the recovered velocity field is closely related to the weighted ray density. Upper panels in Figure 9 show theoretical checkerboards, while lower panels in Figure 9 are reconstructed checkerboards, which reproduce the target velocity values   (Figures 5a, b, and c).  1, 2, and 3 (Figures 5a, b, and c). The upper panels show original checkerboard models, while the lower panels are the corresponding reconstructed checkerboards. Sources (stars) and receivers (triangles) are positioned identically on the irregular surface. quite well, particularly in zones where there are a larger number of rays crossing the medium. This illustrates that our scheme indeed has high resolution in the vertical and horizontal directions.

Comparing with First-Arrival Tomography
As the reflection phase is used to probe a deeper area, compared with the first-arrival tomography (Ma T and Zhang ZJ, 2015), we conduct a first-arrival tomography to make a comparison with the muti-phase tomography. The contrasting results and other relevant outputs are shown in Figures 10a-h. From top to bottom: (a) P-wave velocity distribution depicting the theoretical model, (b) initial model, (c) tomographic image obtained by JI-FRRT, (d) tomographic image obtained by the first-arrival tomography, (e) number of rays that intersect each grid cell for JI-FRRT, (f) number of rays that intersect each grid cell for first-arrival tomography, (g) inverted velocity anomaly for JI-FRRT, (h) inverted velocity anomaly for first-arrival tomography. Since the first-arrival rays can travel only in a shallow area (Figures 10e-f), the velocity difference of the first arrival result is larger than that of JI-FRRT, especially at the deeper area near the reflector (Figures 10g-h).

High-Velocity Anomaly
To check the efficiency of the method we now take a model based on previous Model 3 that irregularly combines hills and valleys ( Figure 4c) and has an intracrustal high-velocity body. We add a high-velocity anomaly (+ 2 km/s) with size 14 km × 14 km just below the central valley of the test structure (Figure 4c). The inversion result and other relevant outputs are shown in Figures 11a-g. From top to bottom: (a) P-wave velocity distribution depicting the theoretical model, (b) initial model, (c) tomographic image obtained by JI-FRRT, (d) ray diagram from the 3rd shot point, (e) number of rays that intersect each grid cell, (f) theoretical velocity anomaly, (g) inverted velocity anomaly. A slight focusing of the rays takes place in the anomaly zone ( Figure 11e). All these results illustrate that the velocity anomaly can be reconstructed with success using JI-FRRT.

Noise Tolerance
The inversion algorithm is known to work well with noise-free data, but this might not be the case with noise-contaminated data  (Zhou B et al., 1992). To determine the response of JI-FRRT to data contaminated by noise, we added random Gaussian noise with standard deviation of σ = 0.2 s to the synthetic data computed from Model 3 (Figure 4c). The inversion result and other relevant outputs are shown in Figures 12a-g. From top to bottom: (a) theoretical model, (b) initial model, (c) tomographic model obtained by JI-FRRT, (d) ray diagram, (e) ray coverage, (f) relative velocity error inherent to the initial model, (g) relative velocity error regarding the final inverted model. Although the velocity error is a little larger than in the noise-free case (lower panel in Figure 7c), the new tomographic image (Figure 12c) fits well to that of the theoretical model (Figure 12a) and is comparable to the image obtained before without noise (lower panel in Figure 5c). This further supports the robustness of JI-FRRT by demonstrating that our inversion scheme works well even with noisy data.

Conclusion
We present a tomographic inversion method to determine the seismic velocity structure of a physical domain whose topograph-ic surface is irregular or non-flat. To deal with this physical feature, we employ the mathematical tool called irregular surface flattening, which is based on the use of a boundary conforming grid and the transformation between curvilinear and Cartesian coordinates. The method has been designed to jointly invert first-arrival and reflection traveltimes (JI-FRRT). The fast-sweeping method (FSM) is applied to a topography-dependent eikonal equation (TDEE) in curvilinear coordinates to obtain first-arrival traveltimes; the multistage scheme is introduced to combine with FSM to propagate reflected waves from an interface. The reflection paths are traced along the negative gradient direction of the traveltime field from receiver to source. We calculate slowness perturbations using first-arrival and reflection paths and the back-projection algorithm, and then we sum these results to update the velocity model.
Numerical examples reveal that JI-FRRT performs well. Tests performed with models of increasing topographic complexity and an undulated reflection interface support the reliability and accuracy of JI-FRRT on the basis of the final tomographic model, ray diagram and ray density coverage, small error in the solution, and rapid convergence of the computational algorithm. Moreover, the classical checkerboard test also indicates that the spatial resolution of JI-FRRT is quite acceptable. The efficiency of JI-FRRT in detecting a velocity anomaly in the interior of the crust is another attractive feature of the method. Finally, JI-FRRT shows a good noise tolerance, which is an important issue when one is working with field data.
In the near future we hope to apply the inversion scheme to other later arrival phases and to include the interface inversion in the scheme, to better apply it to more real travel-time processing.