3-D multi-parameter type traveltime tomography in a spherical coordinate frame: comparison of double and triple class simultaneous inversions

1.   Introduction

Earth structure, and hence the seismic velocity field generally varies both vertically and laterally. Traveltime tomography for double-parameter class inversion, such as the velocity model and the hypocenter locations in local and regional earthquake tomography (i.e., Pavlis and Booker, 1980; Thurber, 1983; Protti et al., 1996; Asad et al., 1999), or for the velocity model and the reflector geometry in the reflection tomography (i.e., Williamson, 1990; Zelt and Smith, 1992; McCaughey and Singh, 1997; Rawlinson et al., 2001; Zelt et al., 2003) has been fully developed in recent years, especially at the local scale in Cartesian coordinates. However, triple-parameter class inversion has seldom been done, especially at the regional or global scale in spherical coordinates (i.e., Sambridge, 1990; Preston, 2003).

The main obstacles for simultaneous inversion for multiple classes of parameters at the regional or global scale are that one needs a fast and accurate multi-phase ray tracing algorithm working in a real 3-D spherical coordinate frame, and there is a strong bias towards emphasizing one class of model parameters over the others during the inversion process (for example, traveltime derivatives with respect to velocity variations, depth changes of the reflectors and deviations of source locations), because their units are quite different and the sensitivity magnitudes may vary substantially.

To deal with the first issue, there have been developed different approaches to multi-phase arrival two-point 3-D ray tracing algorithms such as the perturbation ray theory method (Bijwaard and Spakman, 1999), the dynamic ray tracing method (Dahlen et al., 2000), the revised pseudo-bending method (Zhao, 2001) and the generalized ray theory method (Ballard et al., 2009). These two-point ray tracing methods require that the ray tracing process is repeated for every source-receiver pair. As a result, it consumes a significant amount of CPU time when on the order of a million rays for different phases are needed in iteratively updated nonlinear tomographic reconstruction. Furthermore, for the ray bending or pseudo-bending algorithms, there exists a local minimum problem, whereas for the dynamic ray tracing method or the ray perturbation theory method, the difficulty of finding source-receiver raypaths increases with the complexity of the velocity model.

To tackle these problems, there are four kinds of simultaneous inversion algorithms for balancing the tradeoff between the double-parameter class (velocity and hypocenter location or velocity and reflector depth) updating process: (1) separate the different parameter types in the model space (i.e., SVD method, Pavlis and Booker, 1980); (2) introduce basis vectors in the different directions of the model space, corresponding to the gradient variation of the misfit functional with respect to each type of model parameter (i.e., subspace method, Kennett et al., 1988); (3) define different weighting factors to balance the two types of model parameters in the inversion (i.e., weighting factor method, Zelt and Smith, 1992) and (4) normalize the different sub-Jacobian matrices corresponding to each parameter class (i.e., traveltime derivatives with respect to velocity variations, depth changes of the reflectors and deviations of source hypocenter locations) before the inversion (i.e., normalizing factor method, McCaughey and Satish, 1997; Hobro et al., 2003). There also exists the possibility of some modified or combined version of the above methods, which may be incorporated with the regularization scheme (i.e., to seek a smallest model or to find a smoothest model, Greenhalgh et al., 2006).

To date only two research papers (Sambridge, 1990; Preston, 2003) have tackled the triple-parameter class traveltime inverse problem. Sambridge (1990) worked with a combined passive and active seismic data set and for the forward modeling used a 3-D boundary value ray tracing algorithm (Sambridge and Kennett, 1989). The latter was incorporated with a subspace inversion solver (Kennett et al., 1988) to yield the velocity model, the reflector configuration and the earthquake hypocentres. Ten years later, Preston (2003) carried out simultaneous inversion for the same three classes of model parameter by combining the first arrival times with wide-angle reflection traveltime data. The forward modeling involved two parts: a finite-difference eikonal equation solver (Hole, 1992) was used to find the first arrivals and a two-way approach (Hole and Zelt, 1995) was applied to predict the primary reflected arrivals. For the inversion, the objective function was regularized and the three different classes of Jacobian elements were weighted differently and a conjugate gradient least squares algorithm (Paige and Saunders, 1982) was subsequently applied.

With the recent increase in computational resources, grid/cell-based wavefront propagation methods (referred to as multiple-point ray tracing to distinguish from the above two-point ray tracing approach) have been rapidly developed for regional or global ray tracing. This includes the multistage fast marching method (referred as FMM, de Kool et al., 2006) and a multistage irregular shortest-path method (referred as ISPM, Huang et al., 2013). These two grid/cell based wavefront propagation methods are advantageous over the traditional two-point ray tracing schemes for some applications (such as traveltime tomography) for the following reasons: (1) they are capable of finding first arrivals within continuous media, whereas for traditional two-point raytracing algorithms there is usually no guarantee that the predicted arrival is a first or later arrival; (2) they guarantee that a global solution will be found regardless of the velocity complexity; (3) they are highly efficient for common source point gathers when compared with two-point raytracing algorithms and (4) they are able to compute the raypaths and the corresponding traveltimes for later arrivals at all sampled nodal points within the complex velocity model by using an outward expanding wavefront scheme.

As is well known, the efficiency of tomographic inversion is heavily dependent on the efficiency of the ray tracing scheme (more than 90% of run time is consumed in the forward modeling part). For this reason, it is worthwhile to examine the performance of the newly developed ray tracing algorithm when incorporated with a suitable inversion subroutine to enable simultaneous traveltime inversion. To do so, we developed a new simultaneous inversion algorithm for triple-parameter class updating in Cartesian coordinates (Bai et al., 2015) and a new simultaneous inversion algorithm for double-parameter class updating in spherical coordinates (Huang et al., 2014). In this paper we go further and present an approach for three-class parameter inversion in spherical coordinates and then compare it with the two-parameter class inversion approach. The results and comparisons of synthetic tests between the two and three model parameter class inversions are used to establish the feasibility, efficiency and possible limitations of regional and global scale tomography.

2.   Forward modeling method

Previously we extended the functional of the multistage irregular shortest-path method in Cartesian coordinates (Bai et al., 2010) to a spherical coordinate system (Huang et al., 2013), and were able to predict traveltimes and corresponding raypaths for more than 49 global seismic phases with an absolute time error less than 0.1s when compared with the AK135 global traveltime tables (Kennett et al., 1995). Here we briefly summarize the algorithm.

2.1   Earth model parameterization

Basically we can utilize irregular cells near subsurface interfaces or velocity discontinuities and regular cells elsewhere in the model. For the parameterization, we first divide the subsurface model into several layers according to the major interfaces and parameterize each layer with regular or irregular cells. Secondary nodes are then uniformly distributed along the cell surfaces so that each (both primary and secondary) node belongs to a specific (regular or irregular) cell (see Figure 1a). In 3-D spherical coordinates we use a trapezoidal prism (Figure 1b) to divide the 3-D spherical earth model, except for the global and other irregular subsurface interfaces, where a trapezoidal cone or a hexahedral cell is used (Figures 1c and 1d). By similar considerations, the secondary nodes are inserted along the sides of the trapezoidal prism or trapezoidal cone or hexahedral cell. Inside each cell there are no nodes at all but sources and receivers may be located there.

Figure  1.  Model parameterization in 3-D spherical coordinates (a), for clear representation, the secondary nodes are not depicted, and three kinds of cells used to divide the spherical model (b–d), in the figure the larger circles are primary nodes and the smaller circles are the secondary nodes, respectively

DownLoad: Full-Size Img PowerPoint

2.2   Velocity interpolation function

In general, in 3-D spherical coordinates the velocity field is sampled at the primary nodes of the cell and a specified velocity function is defined across the cell, which links the primary and secondary nodes (including the source and receiver positions in a cell). For a 3-D trapezoidal prism cell, a trilinear Lagrangian interpolation function is used.

where $x_k^e$ and $V(x_k^e)$ $(k = 1, 2, \ldots , 8)$ are the vector co-ordinates and the sampled velocity values at the primary nodes of the cell, respectively. For other special cells (Figures 1c–1d), the velocity values at the secondary node positions are determined by a volume interpolation function. For example, in trapezoidal cone cells the following velocity interpolation function is used to determine the velocity values at the secondary node positions.

where ${v_j}$ is the velocity value of the primary node at that trapezoidal cone cell which includes the ith node and ${u_j}$ is the volume coordinates for that trapezoidal cone cell (for details, refer to Huang et al., 2013).

2.3   Minimum traveltime computation

In calculating the minimum travel times and locating the associated ray paths for all gridded nodes, we can gradually expand outwards from the source to the volume of the computed nodes by continually adding the undetermined neighbouring nodes of the cells to the computed nodes. The wavefront expansion is realized cell by cell. The connection between the nodes is restricted within each cell. In this process one should start with the node that has a minimum traveltime in the subset ${N_j}$ (where ${N_j}$ is the total number of computed nodes in the current wavefront) so as to keep track of the first arrival times for the undetermined nodes. An interval sorting method (Klimeš and Kvasnička, 1994) is used so that the larger traveltimes are deleted, only the minimum traveltime and the related raypath are retained. The minimum traveltime from a source node $i$ to an undetermined node $j$ in a cell can be expressed as

where ${D( {{x_i},\;{x_j}} )}$ is the distance between the current source node $i$ and the undetermined node $j$ , and ${V( {{x_i}} )}$ and ${V( {{x_j}} )}$ are the velocity values at the ith and the jth node positions. Furthermore, in this process the order number of the incident point (node) i* giving the minimum traveltime to the node j is recorded for the coordinates of the related raypath.

2.4   Multistage scheme of computation

According to the layered spherical earth model (Figure 1a), we can divide the earth model into several spherical layers (or shells) corresponding to the major velocity discontinuities. For illustrative purposes consider as an example the Moho discontinuity which separates the crust and mantle. A simulated downwind wavefront (P₁ phase) is propagated through the crust (or upper computational domain) until it impinges on all sampled nodes of the Moho subsurface interface (in our notation convention, P or S represents a compressional or shear wave, respectively, and the subscript or superscript indicates a downwind or upwind seismic wave propagating in different regions, respectively). At this stage the independent computational domain is halted at the upper computational domain and we are left with a narrow band of traveltime values defined along the sampled (Moho) subsurface interface. From here, a downwind propagation of a transmitted wave (P₁P₂) or transmitted and mode-converted branch (P₁S₂ phase) can be simulated by reinitializing it, starting at the sampled node position of the (Moho) interface with the minimum travel time (i.e., according to Huygens Principle, the node is treated as one new source point in the wavefront). Meanwhile, an upwind-propagating wavefront consisting of a reflected branch (P₁P¹ or PmP phase) or a reflected and converted branch (P₁S¹ or PmS phase) can now be obtained by reinitializing the wavefront and starting at the sampled node position with the minimum travel time, from the narrow band of the Moho interface into the incident layer (or crust). Different velocity models (i.e., P or S) are used if wave-mode conversion occurs at the subsurface interfaces. In summary, the multiple arrivals are the different combinations or conjugations, via velocity discontinuities (i.e. subsurface interfaces), of the incident, reflected and converted phases which obey Snell’s Law, Fermat’s Principle and Huygens Principle.

Meanwhile, the computational efficiency is comparable to the currently used 3-D multiple arrival ray tracing algorithm in 3-D spherical coordinate, for example, the multistage FMM scheme (de Kool et al., 2006). Figure 2 shows an example for multiple raypaths (direct P, reflected P₄₁₀P and P₆₆₀P) computed by the above multistage ISPM ray tracing algorithm in 3-D spherical coordinates. More complicated multi-phase arrivals can be predicted in this way (for details, please see Huang et al., 2013).

Figure  2.  Multiple ray paths in a 3-D regional velocity model. In the figure the source is located at a depth of 100 km below the model surface center and the receivers are uniformly arranged along the boundaries of model surface, except for the invisible (back) boundary. The direct P arrivals (blue lines), the primary reflections (P₄₁₀P) from the 410 km subsurface interface (grey lines) and the primary reflections (P₆₆₀P) from the 660 km subsurface interface (black lines) are shown. The 3-D velocity model used is shown in Figure 5

DownLoad: Full-Size Img PowerPoint

3.   Inversion method

3.1   Subspace inversion solver

Previously we compared the subspace inversion solver (Kennett et al., 1988) with the DMNLS-CG inversion solver (the damped minimum norm, least-squares and constrained problem with a conjugate gradient approach, Zhou et al., 1992) for simultaneous inversion of two classes of model parameter in local reflection tomography (Huang et al., 2012). The same two inversion solvers can be used in triple-parameter class inversion (Bai et al., 2015), in which case we found the subspace inversion method to perform better than the DMNLS-CG approach. For this reason, we prefer to use the subspace inversion method (e.g., Kennett et al., 1988; Sambridge, 1990; Williamson, 1990; Kennett and Sambridge, 1998; Rawlinson et al., 2001), because it minimizes the misfit or objective function between the calculated (or theoretical) and observed travel times, by simultaneously adjusting the three types of model parameters (velocity field, reflector depths, hypocenter coordinates), subject to regularization constraints. Our description of the procedure follows closely that given by Rawlinson et al. (2001).

Assume we have a model vector m of length M and a data vector d of length of N, then the discrepancy between d=g (m) and d_obs ≈g(m_true), where g is the forward model operator governing the physics of the process, specifies the accuracy of the inverted model m. Subject to Gaussian errors in the observed data d_obs, the objective function S(m) to be minimized is given by

where C_d and C_m are covariance matrices of the data and model, respectively, and m₀ is an initial model. The damping factor $\varepsilon$ determines the trade-off between the data fit and closeness to the initial model. $\eta$ is a smoothing factor and D is a second-order differential operator. Normally, the last term in the right hand of equation (4) is a smoothing term (structural constraint), which seeks a smoothed model (Sambridge, 1990). Because our model parameterization (e.g., interface sampling spacing) is naturally smooth and the nodal spacing is relatively uniform, we prefer to apply such a smoothing function after each iterative inversion. So in the following discussion, we neglect the smoothing term in the objective function as follows

With the usual quadratic approximation about the current model, we can write for the perturbed objective function (i.e., Rawlinson et al., 2001),

where $\delta {{m}}$ is the perturbation to the current model, ${\overset{\frown} {{\gamma}}} =$ $\partial S/\partial {{m}}$ and ${\overset{\frown} {{H}}} = {\partial ^2}S/\partial {{{m}}^2}$ are the gradient vector and the Hessian matrix.

The main difference between the subspace method and other gradient-based methods is that the subspace method works by constraining the minimization of the quadratic approximation (or objective function) to an n-dimensional subspace of the entire model space,

where ${{A}} = [{{{a}}^j}]$ is the M $\times$ n projection matrix and ${{{a}}^j}$ is a basis vector in the j direction. The quantity ${\mu _j}$ determines the length of the corresponding vector ${{{a}}^j}$ that minimizes the quadratic form of S(m) in the space spanned by ${{{a}}^j}$ . The minimum of the objective function in the n-dimensional subspace is obtained when

Knowing ${{\gamma }}$ , H and the projection matrix A determines the model update ( $\delta {{m}}$ ) as the solution of a small n $\times$ n system of linear equations. The corresponding expressions for the gradient vector ${{\gamma }}$ and Hessian matrix H are

where ${{G}} = \partial {{g}}/\partial {{m}}$ is the Fréchet derivative matrix (or Jacobian), which includes all three kinds of travel time partial derivatives with respect to the changes of velocity, reflector depth and source coordinates. If g(m)–d_obs is not too large and locally linear, then one can neglect the second term in the H matrix and the new or pseudo (approximate) H matrix is then simply written as follows

For weakly non-linear problems, as may arise in global seismology, we can use an iterative approach to update the model,

where m_i and m_i+1 are the current model and the next updated model, respectively.

If the errors are uncorrelated, and then we can define the covariance matrices ${{{C}}_d} = \{ {\delta _{ij}}{(\sigma _d^j)^2}\}$ and ${{{C}}_m}=\{ {\delta _{ij}}$ ${(\sigma _m^j)^2}\}$ . The square root of each non-zero element in C_d and C_m thus indicates the estimated uncertainty in the corresponding traveltime and initial model parameter, respectively. In real data inversion, the $\{ \sigma _d^j\}$ can be estimated from the picking error of each traveltime. The $\{ \sigma _m^j\}$ estimates are based on a priori information on the error associated with the initial model.

Through a suitable choice of basis vectors, the subspace method can avoid some problems encountered in the gradient method, such as a strong dependence on the model parameter scaling and poor convergence when multiple classes of model parameters have to be simultaneously inverted for. If the basis vector { ${{{a}}^j}$ } is constructed in such a way that we choose the steepest ascent vector in the model space ${\overset{\frown} {{\gamma}}} = {{{C}}_m}{{\gamma }}$ at each iteration, then the three search directions (or basis vectors) a¹, a² and a³ in the model parameter space of velocity, reflector depth and hypocenter location can be formulated by partitioning this vector on the basis of model parameter class,

More basis vectors (or search directions) can be obtained by determining the rate of change of the ascent vectors. This has the effect of increasing the dimensions of the subspace, and thus the rate of convergence of the solution. For example, a further nine basis vectors are obtained by pre-multiplying a¹, a² and a³ with the model space Hessian ${\overset{\frown} {{H}}} = {{{C}}_m}{{H}}$ . Similarly, additional basis vectors can be produced by repeating the process of pre-multiplication of the latest set of vectors by the model space Hessian. With a suitable number of basis vectors thus determined, they are then orthonormalized, e.g., with a Gram Schmidt algorithm as follows

where ${{a}}_{\rm orth}^i$ is the ith orthogonal normalized basis vector. Hence, the number of basis vectors is a compromise between computational effort and rate of convergence. The iterative process for the subspace method can be summarized as follows

(1) determine the number (p) of basis vector, $\left\{ {{{{a}}^j}} \right\}$ (j=1, 2, …, p), by using equation (9) to construct ${{{a}}^1} = {{{C}}_m}{\overset{\frown} {{\gamma}}}$ , and if p $\ge$ 1, then using equation (11) to obtain

and then to use some kinds of methods (for example, SVD) to orthogonally normalize the above basis vector ${{a}}_{\rm orth}^i$ ;

(2) to use LU to decompose inverse matrix ${\left[ {{{{A}}^{\rm T}}{{\overset{\frown} {{H}}} { A}} } \right]^{ - 1}}$ , and then obtain the current perturbed model $\delta {{{m}}_i}$ by using equation (8);

(3) to iteratively obtain the current updated model ${{{m}}_i}$ by using equation (12) until the observed traveltimes agree with the model predictions to within the expected noise level or the gradient of the objective function reaches a specified small value. Finally we find a solution ${{{m}}_i} = {{{m}}_{i - 1}} +$ $\delta {{{m}}_{i - 1}}$ .

The entire subspace inversion process is iterative, starting from a suitable initial model. At each iteration ray tracing is used to determine the new ray paths and the predicted traveltimes, and the Fréchet matrix is computed explicitly. Subspace inversion is then used to update to the new model by an amount of δm.

3.2   Jacobian (or Fréchet Derivative) matrix

The travel time from source i to the receiver j via a subsurface interface along the ${R_{i \cdot j}}$ ray path in a spherical coordinate system ( $\theta, \phi, z$ ) is given by the ray integration equation

or in discrete summation form for a parameterized cell model

where dl is the segment of the ray trajectory in a specific cell, ${V_c}(\theta, \phi, z)$ is the velocity distribution of the cell, $k$ is the index of the ray segment in the cell, and n is the total number of segments of the ray.

When the inversion involves three classes of the model parameters, the Jacobian matrix comprises three parts (or sub-Jacobians): one is the traveltime derivative with respect to velocity variations; the second is the traveltime derivative with respect to the depth changes of the reflectors and the third is the traveltime derivative with respect to the source hypocenter location changes.

where ${v_k}$ and ${z_k}$ are the velocity distribution in the kth cell and the depth of sampled kth reflector, respectively. $\Delta {v_k}$ and $\Delta {z_k}$ are the relative perturbations of ${v_k}$ and ${z_k}$ . The quantity $\Delta {\theta _k}$ (k=1, 2 and 3) is the relative perturbations of source hypocenter locations in the latitude ( $\theta$ ), longitude ( $\phi$ ) and depth (z) directions, respectively. The N is the number of velocity nodes along the i-j raypath.

The first derivative with respect to the velocity variation within a cell can be approximated by replacing the integral (equations 15 and 16) with a summation over a series of the averaged values (traveltime derivative with respect to velocity variation) between the two end values on a segment of the ray path across the cell (for simplicity, here we omit the source index i),

where $k$ is the number of segments of the ray path across the cell, and m_i and m_j are the two end points on the mth segment in the cell.

When a ray reaches a subsurface interface, the second term of the derivative with respect to the depth change of reflector can be approximated by the following equation (Sambridge, 1990)

where ${z_n}$ is the depth coordinate of the interface node, ${h_{\rm{int} }}$ is the displacement normal to the interface at the point of intersection by the ray and ${z_{\rm{int} }}$ is the depth coordinate of the intersection point. This is shown in Figure 3, along with the unit vectors ${{{w}}_n}$ and ${{{w}}_j}$ . The derivative of equation (19) can be approximated by

where ${v_j}$ and ${v_{j + 1}}$ are the velocities on the two sides of the interface. If only reflections are considered here, then ${{{w}}_{j + 1}} {{{w}}_j} = - {{{w}}_j} {{{w}}_n}$ and ${v_j} = {v_{j + 1}}$ . Equation (20) then takes on the simpler form.

The third type of Fréchet derivative specifies the arrival time derivative with respect to the source parameters (location and origin time) changes. The derivative of the arrival time t ( $t = {t_{\rm arriv}} - {t_{ 0}}$ ) with respect to origin time is simply given by

If the perturbation of the source location is $\Delta {z_h}$ (see Figure 4 for details), the first-order accuracy expression for the partial derivative for the three spatial directions are as follows (i.e., Thurber, 1983; Eberhart-Phillips, 1986; Sambridge, 1990)

where θ is indicated in Figure 4, ${\omega _1}$ and ${\omega _2}$ are the angles the horizontal projection of the ray subtends with the latitude and longitude directions, respectively.

Figure  3.  Plane wave incident on a perturbed planar interface for first-order approximation of interface Fréchet derivatives with respect to depth change of reflector

DownLoad: Full-Size Img PowerPoint

Figure  4.  Schematic diagram showing how a change in depth Δz_h of the hypocenter results in a change in the path length Δl of the ray, which can be used to find the derivative of the traveltime derivative with respect to the source location change

DownLoad: Full-Size Img PowerPoint

In addition, the discretization of the model is achieved through the primary and secondary nodes. All the nodes are involved in the forward modeling, but only the velocity values at the primary nodes are updated during the inversion. In this way, it is possible not only to maintain the accuracy of the ray tracing, but also to reduce the number of inversion parameters. In similar fashion to the inversion for the velocity field, the reflector depths are updated only through the primary nodes. If the reflector is not located at a primary node position (i.e., at a secondary node position), then the updated value at this secondary node is distributed to the primary nodes of that cell through a bilinear Lagrange interpolation function. Note the source location updating can be at any position, but restricted within the model volume.

It is well known that the units of the three different types of partial derivatives forming the Jacobian matrix are different and their magnitudes may vary substantially, because one class of model parameter is the reflector depths and the other classes are the velocity values and the source locations. This would result in a strong bias towards changing one type of model parameter over the others during the inversion process unless normalization is first applied. Instead of normalizing the entire model parameters, here each type of model parameter is normalized by the factor (Hobro et al., 2003)

where (i=1, 2, 3) indicates the three types of model parameters, and the summation is over all Fréchet derivatives of that parameter type. The new Fréchet derivatives are obtained by multiplying with the reciprocal of n_i. This balances the gradient of the misfit function evenly between the three types of model parameters. The three norms each contain a single normalization constant which is designed to bring all three terms close to unity. This prevents the number of data and the model geometry from having a strong effect on the different model parameters, and is required to achieve a stable inversion (e.g., Huang et al., 2012).

4.   Synthetic inversion tests

In order to numerically analyze the 3-D simultaneous inversion for the triple class of parameters and compare it with the double-parameter class inversion, we selected a regional model at a scale of 20°×20° with depth of 700 km. The two-dimensional cross section of the heterogeneous velocity field is shown in Figure 5a, it is uniformly extended in the latitude direction from 0° to 20°. Two irregular subsurface interfaces were introduced (see Figure 5b), one is located at a depth of the 410 km and the other at 660 km, to simulate the two major discontinuities in the mantle. The maximum amplitude of the vertical fluctuation of both subsurface interfaces is ±20 km. Note that these two velocity discontinuities in the upper mantle are quite complicated, sometimes they maybe appear as gradient velocity zones. Moreover, the reflected phases P₄₁₀P and P₆₆₀P are rarely identified and used in real applications. Currently, the P₄₁₀S and P₆₆₀S reflected and converted phases can be used in the teleseismic receiver-function method (i.e., Bostock, 1995; Bostock and Cassidy, 1996; Shen et al., 1996; Vinnik, 1977) by using a complicated phase-stacking processe (i.e., Rawlinson and Kennett, 2004). For numerical simulated data inversion, we still use these two velocity discontinuities as the reflecting interfaces, but other velocity discontinuities are also possible, such as the Moho. A 1°×1° with depth of 50 km cell is used to parameterize the 3-D regional velocity model, and a secondary node spacing of 10 km (corresponding 9 secondary nodes per cell surface) is applied to improve the ray angular coverage. The grid spacing on the subsurface interfaces is the same as the secondary node spacing. The 30 sources are located within the depth range 20–100 km (seeFigure 6a for source locations projected onto the horizontal plane and Figure 6b for source locations projected onto the longitude-depth cross section). The 441 receivers are uniformly arranged over the top “recording surface” of the model. A combination of direct P arrivals and primary reflected P arrivals from both the 410 km (P₄₁₀P) and 660 km (P₆₆₀P) discontinuities are used to conduct the simultaneous inversion.

Figure  5.  Velocity model showing one cross-section (a) and two undulating subsurface interfaces (b). Note that the 3-D regional velocity is obtained by extending the cross-section velocity field over the latitude range 0°–20°

DownLoad: Full-Size Img PowerPoint

Figure  6.  Thirty randomly located sources within the depth range 20–100 km used in the simultaneous inversion (a) sources projected onto the horizontal plane and (b) sources projected onto the longitude-depth cross section

DownLoad: Full-Size Img PowerPoint

For a fair comparison, three different schemes are trialed for the two cases of noise free and noise added. The first is to simultaneously update both the velocity field and the source locations (referred to as the VS method) with the two subsurface interfaces fixed at their true positions in the inversion process. The second is to simultaneously invert for both the velocity field and the reflector positions (referred to as the VR method) with the source locations fixed at their true positions in the inversion process. The third scheme is our newly proposed method to update all three model parameter classes (referred as the VSR method) simultaneously, with perturbed source locations being used as the initial or starting model values. The initial velocity model is the background flatlayer (1-D) velocity distribution. The initial configuration of the two subsurface interfaces is set as flat at depths of 410 km and 660 km, respectively. To simulate the source location uncertainty, we randomly perturb the real source locations by maximum amounts of ±2°, ±2° and ±10 km in latitude, longitude and z directions, respectively (these are large deviations for such a model, especially in the epicentral coordinates).

The numerical tests with random noise added to the arrival times have maximum errors ±0.2 s, ±0.3 s and ±0.3 s for the direct P arrivals, and the primary reflections P₄₁₀P and P₆₆₀P, respectively. In both noise-free and noisy data inversions, the search dimensions (number of basis vectors) for the subspace solver is 6 (we tried values from 2 to 8, but 6 was the best, see Huang, 2014) and the damping factor is 0.01 (soft one) for all three inversion approaches. For the noise-free case the inversion iterations are terminated when the travel time RMS reaches 0.05 s, whereas for the noise added case the inversion iterations are terminated when the travel time RMS reaches 0.15 s. In the noise-free case the travel time RMS starts at 1.83 s, 9.15 s and 9.65 s for the VR, VS and VRS methods, respectively, whereas in the noise added case, the corresponding values are 1.96 s, 9.68 s and 10.15 s, respectively. In general, the results of the noise free case are better than those of the noise added case. So in the following discussion, we only analyze the results of the noise added case.

Figure 7 shows the inverted velocity fields in cross-section view (latitude of 10°) using the three different simultaneous inversion algorithms. It can be observed that the reconstructed velocity fields have roughly the ideal shape, similar to the real (true) velocity model (a), regardless of whether the VRS (b), VR (c) or VS (d) method is used, but is only partially recovered in the deeper parts of the model due to the relatively low ray density there. To analyze the recovered velocity amplitude, we show the same cross-section view as above in Figure 8, but with the percentage difference between the real and the inverted velocity field. Theoretically, if the simultaneous inversion schemes works perfectly, the above percentage difference values would be zero everywhere (see Figures 8b, 8c and 8d). But in practice the reconstructed velocity field is dependent on the initial model being close to the background velocity of the real model. We see in Figure 8 that only small number of anomalies is left with the same shaped anomalous structures. The less the residual anomaly is present, the better the simultaneous inversion algorithm performs. In this sense, the image of the double parameter type (VR or VS) inversion is slightly better than that of the triple parameter type (VRS) inversion. It is not just the residual anomalies left but also the anomaly shape which should be evaluated for the tomographic procedure. From Figure 8 we see that the result of the VR inversion (Figure 8b) yields superior anomaly shape when compared to the result of the VS inversion (Figure 8c). This can be mainly attributed to the initial source positions with relatively large deviations (maximum ±2°, ±2° in latitude and longitude directions) from the true source positions used for the VS simultaneous inversion. In such situations, the tomographic effort should be much more on adjusting the source positions, which in turn may bring some artifacts in the reconstructed velocity fields, due to the trade-off between the two classes of parameter. Interestingly, such a problem is partially compensated in the result of the VRS simultaneous inversion (Figure 8d), which may have resulted from introduction of the third model parameter type in process. Then the trade-off is between the velocity field, the source locations and the reflector geometries.

Figure  7.  The reconstructed velocity fields in one cross-sections (latitude of 10°) obtained by three different simultaneous inversion algorithms in the noise added case (a) real velocity model; (b)–(d):the inverted velocity model with the VR, VS and VRS method, respectively

DownLoad: Full-Size Img PowerPoint

Figure  8.  The percentage differences between the real and the initial velocity fields (a), and between the inverted and the real velocity fields obtained by the three different simultaneous inversion methods (b: VR method; c: VS method and d: VRS method) in one cross-section (latitude of 10°) in the noise added case

DownLoad: Full-Size Img PowerPoint

The inverted reflector geometries for the two subsurface interfaces are shown in Figure 9. Similar to the velocity field reconstruction, the inverted configurations of the two interfaces are similar to the true model, regardless of whether the VS or the VRS simultaneous inversion method is used. From Figure 9 it is hard to tell which inversion method performs better–the double parameter class inversion (Figure 9a) or the triple parameter class inversion (Figure 9b). To further analyze the performance of the above two simultaneous inversion schemes, we present, in Table 1, the averaged distances between the initial and true subsurface interfaces, and between the inverted (VR and VRS) and true subsurface interfaces. In general, the updated interface depths of the VR inversion are much closer to the known depths than those of the VRS inversion (see Table 1). This is to be expected because the triple parameter type inversion has more degrees of freedom.

Figure  9.  The inverted depths for the two subsurface interfaces by the double model parameter type inversion (a) and the triple model parameter type inversion (b)

DownLoad: Full-Size Img PowerPoint

Table  1.  Statistical results of deviations between the initial and true subsurface interfaces and between the updated and real subsurface interfaces for both VR and VRS simultaneous inversion.

Subsurface interface in depth (km)	Initial distance from true ones (km)	Distance after the VR inversion (km)	Distance after the VRS inversion (km)
410	12.76	1.42	2.39
660	12.76	1.55	2.63

 | Show Table

DownLoad: CSV

In the simultaneous VS and VRS inversions of the noise-free and noise added cases, the initial (starting) source hypocenter locations are randomly perturbed by large amounts along the latitude and longitude directions from the true positions. In order to analyze the performance of two different inversion schemes for earthquake location, we give histograms in Figure 10, which shows the number of sources (or frequency) versus the perturbed values of actual (hypocentral distance), horizontal (epicentral distance) and vertical distance (focal depth) from the true source coordinates in Figures 10a, 10b and 10c, respectively. The comparable plots for the inverted locations, and their deviations from the true values, are depicted in Figures 10d, 10e and 10f for the VS inversion method and in Figures 10g, 10h and 10i for the VRS inversion method. From Figure 10, it can be seen that the initial source positions have large deviations from the true source positions (the maximum deviation is 320 km, see Figure 10a), especially in horizontal or epicentral distance (the maximum deviation is 240 km, see Figure 10b). After simultaneous inversion the updated source positions converge to the real source positions within 10 km deviation for the VS inversion (Figure 10d, most within 5km) and within 15 km deviation from true source positions for the VRS inversion (Figure 10g, most within 5 km). The parameter updating of source positions for the two different inversions mostly occurs in the horizontal direction (see Figures 10e and 10h), and not the vertical or focal depth direction (see Figures 10f and 10i). In general, the double parameter (VS) inversion performs slightly better than that of the triple parameter (VRS) inversion scheme.

Figure  10.  Perturbed distribution of number of initial sources against the real source locations along the hypocentral distance (a), epicentral distance (b), focal depth distance (c), and the corresponding distribution of number of inverted sources converged to the real source locations after the two simultaneous inversion methods (d, e and f: for the VS method; g, h and i: for the VRS method)

DownLoad: Full-Size Img PowerPoint

From the above numerical tests and comparisons, it is evident that the newly developed simultaneous triple parameter (VRS) inversion method is capable of recovering all three sets of model parameters (velocity field, reflector geometry and source positions) from using multiple phase traveltime information.

It should be noted that in real earthquake applications, not only the subsurface interfaces are poorly known, but also the source locations have uncertainties due to the uncertainties in the underlying velocity field. Furthermore, the origin time of the earthquakes has also a certain degree of uncertainty. We can of course use arrival time (or travel time) differences to remove the incorrect origin time influence. For the imprecise source locations, one may use two kinds of simultaneous inversion, for example, to update both the velocity and the source hypocenters for several iterations and then invert both the velocity and reflector geometry simultaneously for several more iterations. This can be referred to as a sequential or cascaded inversion approach. Alternatively, one can directly invert all three model parameters types simultaneously as we do in this paper within a spherical coordinate frame for regional/global tomography or in a companion paper (Bai et al., 2015) within a Cartesian coordinate frame for local earthquake tomography.

5 Conclusions

In this paper we have combined the spherical coordinate multistage shortest-path ray tracing algorithm (Huang et al., 2013) with a popular subspace inversion solver (Rawlinson et al., 2001) to formulate a simultaneous inversion scheme to update the velocity field, the reflector geometry and the hypocenter locations, using travel times from multiple classes of arrivals. Numerical tests are performed for the triple model parameter class inversion and compare it with the double model parameter class inversion scheme for two cases (noise free and noise added). It is clear that the VRS simultaneous inversion algorithm has nearly the same capability to concurrently update both the velocity field and the hypocenter locations when compare with the VS inversion algorithm, and roughly equal performance to simultaneously invert both the velocity field and the reflector geometry when compare with the VR inversion approach for the two cases. The important finding is that with large random deviation of the initial source positions from the true values the simultaneous triple parameter type inversion scheme is still able to obtain comparable image results to the double model parameter type updating procedure. All three schemes (the VR, VS and VRS inversions) are quite tolerant to modest errors in the travel time picks, which make them potentially robust for practical applications. The results show that the VRS simultaneous inversion method is an efficient and feasible way to recover velocity information together with reflector geometry and hypocenter locations for regional tomographic problems provided that arrival time information is available for multiple phases (e.g., direct and reflected). Furthermore, the output model of the triple class parameter inversion could be a good initial model for later full waveform inversion.

The benefits of the proposed three class model parameter inversion algorithm are that (1) it can be used to further tune the model parameters to approximate the real values by including as many kinds of different arrivals as possible subsequent to two class (velocity and reflector geometry or velocity and source locations) model parameter inversion. For instance, reflection or refraction tomography, is primarily conducted to seek a first guess model; (2) the output model of the triple model parameter inversion can be used as an initial model for later full waveform inversion, where a close initial model to the true model is needed due to the high nonlinearity of the inverse problem.

Ackowledgement

This study was partially supported by the Doctoral Programming Research Fund of Higher Education, Chinese Ministry of Education (No. 20110205110010).