Walkaway vertical seismic profiling first-arrival traveltime tomography with velocity structure constraints

1 Introduction

An ideal velocity model is essential in seismic exploration [1]. Currently, walkaway vertical seismic profiling (WVSP) technology has the advantages of one-way wave attenuation, high signal frequency, geophone depth positioning, and high data accuracy and is therefore an effective exploration tool in complex areas [2]. Traveltime tomography using WVSP first-arrival traveltime is a valuable tool for obtaining a deep velocity model. Luming et al. [3,4] utilized the resolution advantage of the first-arrival tomography inversion method in the longitudinal and transversal directions to realize the inversion of a complex surface model capable of adapting to arbitrary changes in velocity. Later, the method was used in three-dimensional surface models and achieved good results. Shijun and Jianzhong [5] used a nonlinear conjugate gradient optimization method to realize the first-to-traveltime laminar inversion, which led to a stable convergence of the inversion process. Wei [6] combined first-to-traveltime velocity inversion and reflection wave traveltime velocity inversion, which achieved a better inversion effect. Aiyuan et al. [7] proposed a double-weighted tomography inversion algorithm for VSP traveltime based on confidence and ray length in order to overcome the difficulty of reflecting wave velocity modeling in complex mountainous regions.

In the velocity modeling process of WVSP first-arrival traveltime tomography, the tomography sensitivity matrix is built based on ray distribution, and the amount of modification of the model parameters is related to the coverage density of the rays [8]. In addition, the validity of the field, the finiteness of the observations, and the errors in the observations and calculations are factors that make the first arrival inversion an ill-posed problem. In order to alleviate the ill-posed problem, it is necessary to introduce additional information and incorporate reasonable constraints (smoothing and regularization). The smoothing constraint is primarily utilized to filter out the impact of high-frequency noise either during the inversion process or in the final results. This enhances the realism of laminar outcomes. In contrast, the regularization constraint enhances the function’s nature throughout the laminar process, acting as an intrinsic smoothing that surpasses typical laminar profile smoothing. [9]. Here, we used regularization as the constructive information constraint.

Since Tikhonov proposed regularization to address the ill-posed problem, regularization constraints have been an important topic in theoretical research in the field of inversion. Zhou et al. [10] first introduced the regularization strategy into the study of cross-hole radar tomography; Fomel used regularization to achieve an overall smoothing effect on the model [11]; Liu et al. [8] introduced a priori information into the first-arrival traveltime tomography through regularization, replacing the traditional external constraint model. Clapp et al. used regularization to incorporate information on the dip angle of the formation to improve the accuracy and convergence of the inversion [12,13]; Qin et al. incorporated acoustic and time difference logging data into regularization to constrain the inversion process, resulting in more accurate inversion results [14]; Hong et al. combined the position and dip angle information of geological formations to construct a Gaussian smooth matrix for the Gaussian bundle lamination inversion with regularization constraints, which improves the stability of the inversion process and inversion speed, and renders the inversion results more geologically reasonable [15].

These studies show the necessity and effectiveness of constraining the original inversion process by including prior information in the process through regularization. The advantage of WVSP velocity modeling is that it can be further calculated based on the ground seismic processing results. Therefore, in this study, we extracted seismic information and performed tectonic block division of the model based on the ground seismic layer and fault information as a priori data to constrain VSP velocity modeling. Considering the effectiveness of the regularization constraint, we transformed the tectonic block division information into regularization factors, added the structure information to the WVSP first-arrival traveltime inversion process through regularization, established a regularized tomography inversion equation to improve the complexity of the WVSP first-arrival traveltime tomography inversion equation, and updated the velocity model using multiple iterations. Finally, we used theoretical simulation experiments to prove the effectiveness of the constructive information regularization constraint for the WVSP arrival traveltime tomography inversion.

2 Methods

2.1 WVSP first-arrival traveltime tomography

Under the high-frequency approximation theory, seismic waves can be considered to propagate along a ray path within a slowness field, where slowness is the reciprocal of the velocity. The first-arrival traveltime can then be considered a radon-positive transformation along the ray path:

(1) t ( s , r ) = ∫ L ( s , r ) m ( x , y ) d l .

In the above equation, t ( s , r ) represents the traveltime from the excitation point s to the receiver point r and L ( s , r ) is the ray path from point s to point r , where m ( x , y ) is the slowness of the model at the ( x , y ) position. The slowness model needs to be discretized for calculation in the tomography process, and the general discretization form is grid discretization. Therefore, the above equation has the discrete form

(2) t ( s , r ) = ∑ j N M j Δ l j ,

where Δ l j is the length of the ray in the jth grid and m j is the slowness of the jth grid. As shown in Figure 1, for the corresponding transmitter and receiver, the sum of the product of the slowness of the ray between the excitation receiver and the length of the ray passing through the grid is the traveltime between the two points. Thus, for a single excitation point for each geophone, there exists a ray passing through a number of grids, and if the total number of rays is I and the number of discretized grids is J , then a system of ray traveltime equations can be obtained as follows:

(3) a 1 , 1 a 1 , 2 … a 1 , J a 2 , 1 a 2 , 2 … a 2 , J ⋮ ⋮ ⋱ ⋮ a I , 1 a I , 2 … a I , J m 1 m 2 ⋮ m J = t 1 t 2 ⋮ t J ,

where a i , j denotes the length of the ith ray passing through the jth grid, m i denotes the slowness of the jth grid, and t i denotes the traveltime for the ith ray. The above equation set can be written in a linear form as follows:

(4) A m = t ,

where A is the ray grid matrix, m is the slowness matrix, and t is the traveltime matrix. To facilitate the inversion calculation, given the initial inverse model, the above equations can be transformed into the following incremental equations:

(5) A Δ m = Δ t ,

where Δ t represents the difference between the real detection ray traveltime and that calculated using the initial model.

(6) Δ t = t obs − t cal = A true m true − Am ,

where t cal denotes the ray-traced orthorectified traveltime, t obs denotes the actual detected traveltime, A true denotes the ray propagation path length matrix in real geology, and A denotes the ray-traced orthorectified ray path length matrix in the given model. Under the assumption of “small perturbation of the model parameters with constant ray path” (linear approximation), A can be approximated as A true ; thus, the above equation can be written as follows:

(7) Δ t = A ( m true − m ) = A Δ m .

Figure 1

Schematic diagram of the ray path.

Given the initial model, the ray-tracing forward is computed to obtain the A-matrix and t cal as t obs can be obtained in the detection; thus, Δ t can also be obtained by equation (6). Δ m can be iteratively corrected to minimize Δ t , which will be infinitely close to the real velocity model, and this process is the tomography velocity inversion process. This process can be considered as an optimization problem, and the objective is to find the velocity model that minimizes Δ t . According to the definition of the objective function of the nonlinear inversion problem by Tarantola and Valette [16], the objective function of this optimization can be written as follows:

(8) minimize Φ (Δ m ) = ( A Δ m − Δ t ) T ( A Δ m − Δ t ) .

According to ∂ Φ ∂ Δ m = 0 , it is obtained that

(9) Δ m = ( A T A ) − 1 A T Δ t .

This equation is the least-squares solution of equation (5) [11,16–25].

3 Regularization constraint

To reduce the uncertainty caused by the inversion mixing problem, it is necessary to introduce additional information and incorporate reasonable constraints. The advantage of WVSP over ground seismic traveltime tomography is that it can perform further calculations based on ground seismic processing results and provide richer a priori information. In addition, detailed tectonic interpretation information can be obtained during ground seismic processing by picking up layer and fault information. The tectonic interpretation information is crucial owing to existing tectonic knowledge, experience, and understanding of the work area, which will aid smooth incorporation of information accumulated in the process into the interpretation information. In view of this, we can extract the tectonic information from the ground seismic interpretation and establish a tectonic regularization term to constrain the velocity within the same structure to reduce the influence of “zero space” on the inversion solution and reduce the inversion uncertainty.

The structure interpretation information contains the layer configuration and the location of faults; according to this information, the target modeling area can be divided into structure blocks. As the velocity variation within the same structure block is usually small, the velocities can be considered similar within the same block; therefore, a separate smoothing constraint can be applied within each structure according to the structure block division. The structural blocks can be divided by geosurface fitting, and for the structure data shown in Figure 2, the block division was performed as follows:

Figure 2

Schematic diagram of structural delineation. (a) Tectonic interpretation information, red for layer data and green for fault data; (b) structural delineation; (c) discretization of the grid; (d) marking of the grid within the structural block; (e) classification of the grid markers on the boundary line.

(1) Fitting the interpreted information and ensuring the closure of all structure blocks.

(2) Discretizing the target region, initializing the grid, and labeling the fault and laminar line grids.

(3) Dividing and labeling the grid according to the breadth-first traversal algorithm and classifying the fault and laminar grids as adjacent arbitrary block grids.

After the division is completed, the entire velocity model is considered a collection of different structural blocks:

where a single structure block consists of N divided meshes:

where index i , j represents the jth grid in the ith constructive block, after getting the grid number under the division of all constructive blocks, we can develop the regularization matrix for each constructive block, and for the ith constructive block P i :

The matrix R i is a sparse matrix of size ( N i − 1 ) × M , where N i is the number of meshes of the ith structure block and N is the total number of meshes of the model. The physical meaning of this sparse matrix is to make the solutions of two adjacent positions inside the same constructive block equal and perform regional spatial smoothing of the model in terms of the overall model in terms of the constructive body. Thus, for two adjacent meshes on the model, their constraint positions in the matrix are also adjacent; for example, index i , j is 1 and index i , j + 1 is –1. For meshes on the model that are at the upper and lower levels, their constraints in the matrix are given according to the mesh number; for example, index i , j is 1 and index i , k is –1, and k is the number of the next mesh in the model for that construct. Most of the matrix constraints are adjacent to each other.

The overall structure constraint regularization matrix is obtained by setting the regularization matrix for each structure block as follows:

where R i is the regularization matrix of a single constructive block, R is the complete constructive constraint matrix of size ( N − k ) × N – where N is the number of overall model-division meshes and k is the number of constructive body divisions. The above regularization matrix adds a flat constraint within the structural block to the solution in order for the inversion process to contain structural information from the interpretation data. However, even for the same medium in the same tectonic layer, the propagation velocity of seismic waves varies, which is not as drastic as that of the upper and lower layers of the reflection interface. However, the velocity of the top and bottom of the same layer also exhibits some variation in the large-scale structural block. Therefore, further sublayers are required based on the structural block division to match the velocity variation.

As depicted in Figure 3, the tectonic block can be categorized into two scenarios. The first involves a quadrilateral structure, wherein the upper and lower interfaces serve as the vertical boundaries, while the fault or working zone boundary constitutes the lateral boundaries. According to the characteristics of the geologically stratified cover, this structure is generally divided uniformly into flat layers (Figures 2–4a). The other is a triangular structure, where the top–bottom interface intersects at one end. This structure can be divided isoparametrically (Figures 2–4b) or as flat (Figures 2–4c) based on information related to the underlying layers. In this study, we performed small-layer construction by isometric division. We obtained the regularization matrix acting on the small-layer flat using the same method after small-layer division. Compared to the flat within the construction block, the small-layer division can portray the velocity field on a smaller scale and reduce the effect of the regularization term on the convergence of the solution.

Figure 3

Schematic diagram of small layer division. (a) Quadrilateral configuration with left and right boundaries; (b) triangular configuration with isometric division; (c) triangular configuration with flat layer division.

After obtaining the regularized constraint matrix for complete subdivision, the objective function of equation (8) can be optimized as follows:

where R*delta(m) is added to the objective function as a regularization term and ε is the weight coefficient of this regularization term, which can be adjusted during experimentation, whose least squares solution is shown in the following equation:

This equation is the inversion equation of traveltime tomography with the constraint of constructive information regularization. The initial velocity model is modified by solving the equation, and the modified model continues to be modified to minimize the difference between the real traveltime and the traveltime calculated by the model. The final model is the inversion result.

4 Results

To verify the effectiveness of the regularization constraint on tectonic information, we designed shallow and deep undulation models in this study. We then compared the inversion results with the regularization constraint of velocity structural information with the inversion results of the overall smoothing constraint. The two models were of the same size, with 200 grid points in both the horizontal and vertical directions and a grid spacing of 20 m. The velocity ranges of Models I and II were 1,500–4,500 and 1,500–4,000 m/s, respectively. The observation system was a WVSP, and the shot point spacing was 10 grid points, i.e., 200 m, with a total of 20 shot points. The receiver was located 2,000 m laterally, with a spacing of eight grid points, i.e., 160 m.

Model I was a combination of a shallow undulating structure and a deep flat structure to verify the effectiveness of the regularization constraint of the structural information on the undulating geological structure. Model I and its inversion results are shown in Figure 4.

Figure 4

Model I and inversion results. (a) True model; (b) result of the overall smoothing constraint; (c) result of the velocity structure constraint.

Model I had obvious undulations at a depth of approximately 1,000 m, with a uniform velocity structure below and second-layer and third-layer velocities of 3,000 and 4,500 m/s, respectively. As shown in the inversion results, the overall smooth-constrained inversion exhibited many sharp shapes in the undulating region, and its velocity values also had a large gap from the real model, such as at lateral 800, 1,500–2,500, and 3,500–4,000 m. Overall, the constructive constraint inversion results are able to reflect the real model structure, while the overall smoothing inversion results have large differences from the real model results. We selected lateral velocity slices at a depth of 1,000 m for comparison, and the comparison plots are shown in Figure 5. Figure 6 shows the speed error percentage compared with that of the true model.

Figure 5

Comparison of lateral velocity in undulating areas.

Figure 6

Speed error percentage comparison with the true model.

As shown in Figures 5 and 6, except for the large difference in the velocity values of the inversion results of the overall smoothing constraint in the above area, the velocity values of the inversion results of the structural smoothing constraint in other undulating areas were also closer to those of the real model, indicating the higher accuracy of the structural information constraint for the inversion of the undulating terrain. In addition, we selected the lateral velocity at 3,200 m for comparison to verify the effectiveness of the tectonic constraint for the deep region, and its comparison figure is shown in Figure 7. Figure 8 shows the speed error percentage compared to the true model at 3,200 m.

Figure 7

Comparison of deep lateral velocity.

Figure 8

Speed error percentage comparison with the true model.

As shown in Figures 7 and 8, the overall inversion velocity of the structural constraint is closer to that of the real model in the deep region, whereas the overall smoothing constraint inversion is closer to the structural constraint effect in the middle region where the rays are denser. However, in the regions on both sides where the rays are sparse, the inversion results are different from those of the real model and mainly depend on the initial velocity model. Evidently, the structural information constraint also has a good effect on the inversion of the deep region and can effectively improve the dependence of the ray-sparse region on the initial model.

Model II was a deep undulation model used to verify the effectiveness of structural constraints in deep undulating complex areas. The model and inversion results are shown in Figure 9.

Figure 9

Model II and inversion results. (a) True model; (b) result of the overall smoothing constraint; (c) result of the velocity structure constraint.

Model II was a four-layer model with undulations at depth, and the upper layer of the undulations is a tilted layer, which increases the difficulty of geological inversion of the undulations. Overall, the structural constraint inversion results were closer to those of the real velocity model, and each layer was smoother with the constraint of the velocity structural information. The overall smoothing constrained inversion, on the contrary, was mainly influenced by ray density. There is an inverted triangular shape in the overall WVSP initial-to-ray distribution, and the velocity varies drastically, indicating that the overall smoothing constraint is not effective. Figure 10 shows the lateral velocity comparison at a depth of 1,200 m, and Figure 11 shows the speed error percentage comparison with the true model. The geology of this layer is flat above and inclined below. The accuracy of the structural constraint on the shallow velocity inversion was verified by comparing the inversion results of this layer.

Figure 10

Comparison of lateral velocity in the shallow part.

Figure 11

Speed error percentage comparison with the true model.

As shown in Figures 10 and 11, the overall smoothing constraint inversion results in this layer were poor, the velocity variation was sharp and complicated, and the difference from the real model was significant. The results of the structural constraint inversion are closer to the real velocity value in the whole lateral range, and the velocity variation is smooth, especially in the area on both sides where the ray coverage is sparse; the velocity structural constraint inversion also achieves better results.

Figure 12 shows a comparison of the lateral velocity slices at 3,200 m depth, and Figure 13 shows the speed error percentage comparison with the true model.

Figure 12

Comparison of lateral velocity in deep undulating areas.

Figure 13

Speed error percentage comparison with the true model.

As shown in Figures 12 and 13, the constructively constrained inversion results were closer to those of the true velocity model in the deep undulating region. In particular, in the sparse ray distribution area on both sides, the overall smoothed inversion results depend on the initial model, whereas the constructively constrained inversion results enhance the accuracy of the real velocity model results. In the middle region with a concentrated ray distribution, the overall smoothed inversion results were also closer to that of the true velocity model, whereas constructive constraint inversion yielded better results. Evidently, the velocity structural information regularization constraint inversion has a higher accuracy than that of the traditional overall smoothing constraint inversion in both the overall model and the undulation detail region, and it is effective to use the velocity structural information regularization term as a constraint.

5 Discussion

Our view is that the accuracy of velocity inversion can be improved by utilizing seismic data as the constraint information for WVSP velocity inversion. The effectiveness of the velocity structure as a regularization constraint is verified by numerical experiments. Compared with the overall smoothing regularization [11], which is also derived from the Tickhonov regularization method [26], the velocity structure regularization can contain more information. This allows the velocity structure regularization to produce more reliable results than the overall smoothing regularization. In velocity structure regularization, the choice of the regularization weight parameter ε is also an issue to be considered. ε is too large to lead to bias in the solution and too small to have enough constraining power. The choice of the parameter is explored by Engl and Ramlau [27].

6 Conclusion

In this study, we obtained the velocity structure of complex structures by constructing geotectonic models and modeling the layer sequence based on depositional patterns, using ground seismic tectonic interpretation data to extract stratigraphic and fault tectonic information. We used the velocity structure to divide the model and then constructed the block smoothing regularization operator and transformed it into a regularization term incorporated into the inversion process to constrain the WVSP first-arrival traveltime tomography. The processing results of the theoretical model verify that, compared with the traditional overall smoothing constraint, the accurate structure information regularization constraint can effectively improve the inversion accuracy and render the model structure clearer and smoother with a better constraint effect. The accuracy of the velocity structure information should be considered during the actual processing.