CN114779335B - Elastic wave direct envelope inversion method based on anisotropic total variation constraint - Google Patents

Abstract

The invention relates to a direct elastic wave envelope inversion method based on anisotropic total variational constraints, which is to obtain high-precision longitudinal and Shear wave velocity structure. Firstly, the elastic wave field is decomposed into the wave field mode to obtain the longitudinal and transverse wave fields respectively, and the envelope field is calculated to obtain the forward longitudinal and transverse wave envelope fields respectively; then, according to the gradient expression, the elastic wave direct envelope inversion can be obtained The P-wave and S-wave velocity gradients can be used to update the P- and S-wave velocity models; after that, an anisotropic total variational constraint is applied to the current P- and S-wave velocity models to obtain the constrained P- and S-wave velocity update models; finally, the anisotropy The result of direct envelope inversion of fully variational constrained elastic waves is used as the initial model, and the high-precision longitudinal and shear wave velocity models of elastic strongly scattering media can be obtained by performing anisotropic full variational constrained elastic wave full waveform inversion.

Description

Direct Envelope Inversion Method for Elastic Waves Based on Anisotropic Total Variation Constraint

Technical Field

The invention relates to a method for applying anisotropic total variation constraints in elastic wave direct envelope inversion to improve the inversion effect of longitudinal and transverse wave velocity structures of strong scattering media.

Background Art

Strong scattering media can usually form a good sealing medium for oil and gas resources. Therefore, high-precision parameter modeling methods for strong scattering media have attracted great attention in the field of exploration geophysics. The full waveform inversion method is currently the most accurate method for parameter modeling in the field of exploration seismic, which can make full use of the kinematic and dynamic information of the seismic wave field. In recent years, people hope to use the full waveform inversion method to perform high-precision parameter modeling on strong scattering media (such as salt domes) encountered in actual exploration, so as to improve the imaging quality of strong scatterers and their shielding areas. However, the conventional full waveform inversion method is based on the Born approximation, a weak scattering approximation, and the actual seismic data often lacks effective low-frequency information, which makes it difficult to directly use the conventional full waveform inversion method to perform parameter modeling of strong scattering media. In recent years, domestic and foreign scholars have improved the theoretical framework of the full waveform inversion method from different angles to adapt to the parameter inversion of strong scattering media. Most of these studies are based on acoustic wave media, mainly including the Laplace-Fourier domain waveform inversion method, the standard setting method, the total variation constraint method, the deep learning method, and the direct envelope inversion method. The research on multi-parameter modeling of strong elastic scattering media is still in its initial exploration stage, and the existing methods are not very ideal for modeling the velocity inside and below strong elastic scatterers. In short, there is currently a lack of a method that can perform high-precision longitudinal and transverse wave velocity modeling for strong scattering media in the absence of low-frequency data and without prior information.

The direct envelope inversion method can model large-scale parameters of strongly scattering media when low-frequency seismic data is missing and there is no prior information, and it has relatively high computational efficiency. The biggest difference between the direct envelope inversion method and the full waveform inversion method is that it defines a direct envelope sensitive kernel function (different from the waveform sensitive kernel function of conventional envelope inversion), which can directly map low-frequency envelope data perturbations to large-scale parameter perturbations of strongly scattering media. However, the current elastic wave direct envelope inversion method is difficult to obtain ideal internal velocity information of strong scatterers when modeling the longitudinal and shear wave velocities of elastic strong scattering media, which affects the velocity modeling effect of the lower boundary of the strong scatterer and below it.

Summary of the invention

The purpose of the present invention is to address the deficiencies of the above-mentioned prior art and provide a new multi-parameter inversion method for elastic waves in strongly scattering media, thereby solving the problem of high-precision modeling of P- and S-wave velocities of strongly scattering reservoirs in seismic exploration of oil and gas resources, and providing high-precision P- and S-wave velocity models for imaging of complex strongly scattering reservoirs.

The idea of the present invention is to make full use of the advantage of anisotropic total variation constraint in sharpening the model boundary and enhancing the uniformity in the layer, and introduce it into the elastic wave direct envelope inversion process. In each iteration of the elastic wave direct envelope inversion, anisotropic total variation constraints are applied to the iteratively updated longitudinal and transverse wave velocity models, respectively, so that the interior of the strong scatterer tends to be uniform and the boundary is more obvious. Finally, the high-precision longitudinal and transverse wave velocity model of the strong scattering medium is obtained by combining the anisotropic total variation constraint elastic wave full waveform inversion technology, overcoming the shortcomings of the prior art.

The objective of the present invention is achieved through the following technical solutions:

First, prepare the preprocessed elastic wave multi-component observed seismic data; use the observed seismic data to estimate the source wavelet and give the initial model of P- and S-wave velocities (without prior information); perform forward simulation on the initial velocity model to obtain simulated seismic data and forward seismic wavefield; perform wavefield mode decomposition on the forward seismic wavefield to obtain the P- and S-wave forward wavefields, respectively, obtain the corresponding envelope fields, and obtain the P- and S-wave forward envelope fields, respectively; make a difference between the observed envelope data and the simulated envelope data to obtain the envelope companion source; perform companion envelope field calculation and mode decomposition on the initial P- and S-wave velocity models to obtain the P- and S-wave forward envelope fields. , and shear wave accompanying envelope fields; use the forward envelope field and the accompanying envelope field to calculate the P-wave and S-wave velocity gradients respectively; obtain the step size and update the P-wave and S-wave velocity models; impose anisotropic total variation constraints on the updated P-wave and S-wave velocity models respectively to obtain the constrained P-wave and S-wave velocity models; perform iterative updates until the stopping conditions are met to obtain the large-scale P-wave and S-wave velocity structure of the strong scattering medium; use the large-scale P-wave and S-wave velocity model of the strong scattering medium as the initial model, perform anisotropic total variation constrained elastic wave full waveform inversion, and obtain high-precision P-wave and S-wave velocity modeling results of the strong scattering medium.

The elastic wave direct envelope inversion method based on anisotropic total variation constraint described in the present invention is implemented through the MATLAB platform;

The elastic wave direct envelope inversion method based on anisotropic total variation constraint of the present invention comprises the following steps:

a. Install MATLAB software platform;

b. Perform static correction and denoising preprocessing on the data to obtain high-quality elastic wave multi-component observation seismic data;

c. Perform wavelet estimation on seismic data and extract the source wavelet of each shot data;

d. By analyzing the background velocity, the approximate range of the background P-wave and S-wave velocities is obtained, and a background velocity model is generated, which does not contain any prior information of strong scatterers, and is used as the initial P-wave velocity model v _p0 and the initial S-wave velocity model v _s0 for inversion respectively;

上角标i表示i方向的分量，对于二维情况，指x(水平)和z(垂直)方向；对观测地震数据取包络得到观测包络数据

通过公式(1)计算弹性波直接包络反演的目标函数σ_EDEI：e. Calculate elastic wave multi-component simulated seismic data based on the initial P- and S-wave velocity models, and obtain simulated envelope data by taking the envelope of the simulated data

The superscript i indicates the component in the i direction. For the two-dimensional case, it refers to the x (horizontal) and z (vertical) directions. The observed seismic data is enveloping to obtain the observed envelope data.

The objective function σ _EDEI of elastic wave direct envelope inversion is calculated by formula (1):

In the formula, the subscript sr of the summation symbol indicates the integration of all sources and receivers, t indicates time, and T indicates the total recording time length;

和横波正传包络场

f. Calculate the simulated seismic wave field on the initial model, decompose the simulated seismic wave field to obtain the positive longitudinal wave field and the positive transverse wave field, and take the envelope to obtain the longitudinal wave positive envelope field

and the shear wave positive envelope field

和横波伴随包络场

g. Calculate the difference between the simulated envelope data and the observed envelope data to obtain the companion source, and then transmit the companion source back to obtain the companion envelope field; perform wave field mode decomposition on the companion envelope field to obtain the longitudinal wave companion envelope field

and the shear wave envelope field

h. The P-wave velocity gradient of direct envelope inversion of elastic waves is obtained by performing zero-delay cross-correlation between the P-wave forward envelope field and the P-wave accompanying envelope field, as shown in formula (2):

Where _vp represents the longitudinal wave velocity and ρ represents the density. The shear wave velocity gradient of the elastic wave direct envelope inversion can be calculated by formula (3):

In the formula, _vs represents the shear wave velocity, · represents the dot product, and μ represents the shear modulus;

和

表示，对二者施加各向异性全变分约束的过程等价于求解公式(4)和(5)所示的最优化问题：i. Select an appropriate step size and use the steepest descent method to update the P-wave and S-wave velocity models. Assume that the current iteration number is m, and the P-wave velocity and S-wave velocity models obtained by the current iteration are respectively expressed as

and

It means that the process of applying anisotropic total variation constraints to the two is equivalent to solving the optimization problem shown in formulas (4) and (5):

和

分别为第m次迭代的纵、横波速度的更新量，λ₁和λ₂分别为对纵、横波速度施加各向异性全变分约束的权系数，||·||表示二范数，||·||_ATV表示各向异性全变分范数。对速度模型v计算各向异性全变分范数的具体表达式如公式(6)所示：Where J ₁ and J ₂ represent the objective functions for applying anisotropic total variation constraints to the P-wave and S-wave velocity models, α ₁ and α ₂ represent the update steps of the P-wave and S-wave velocities, respectively.

and

are the updated values of the P-wave and S-wave velocities at the mth iteration, λ ₁ and λ ₂ are the weight coefficients for applying anisotropic total variation constraints to the P-wave and S-wave velocities, ||·|| represents the bi-norm, and ||·|| _ATV represents the anisotropic total variation norm. The specific expression for calculating the anisotropic total variation norm for the velocity model v is shown in formula (6):

Where nz and nx represent the number of grid points in the vertical and horizontal directions of the model, respectively. Solve the optimization problems shown in formulas (4) and (5) to obtain the P- and S-wave velocity models after anisotropic total variation constraints as the velocity model after the current iteration update;

j. On the updated model, the iteration stop condition is judged; if the stop condition is not met, the updated P- and S-wave velocity models are used as the initial model, and the iterative calculation is continued in step e; if the stop condition is met, the output result is the large-scale P- and S-wave velocity structure v _pt and v _st of the strongly scattering medium after the anisotropic total variation constraint;

k. Using v _pt and v _st as the initial models, anisotropic total variation constrained elastic wave full waveform inversion is performed to obtain the final inversion result, namely the high-precision longitudinal and shear wave velocity structure of the strong scattering medium.

Compared with the prior art, the beneficial effect of the present invention lies in that the present invention introduces anisotropic total variation constraints into the direct envelope inversion process of elastic waves, which can make the internal velocity of the strong scatterer uniform and the boundary information prominent, thereby improving the boundary characterization of the strong scatterer and the velocity modeling effect of the shielded area below, and then obtaining high-precision longitudinal and transverse wave velocity inversion results of the elastic strong scattering medium.

The advantages are as follows: 1. The present invention introduces anisotropic total variation constraints, making it easier to obtain high-quality strong scatterer velocity information during the elastic wave direct envelope inversion process. 2. The anisotropic total variation constrained elastic wave direct envelope inversion involved in the present invention can better characterize the boundary information of strong scatterers, especially the velocity information of the lower boundary of the strong scatterer and the shielded area below it. 3. The elastic wave direct envelope inversion based on anisotropic total variation constraints and the elastic wave full waveform inversion method are connected in series, which can take into account the inversion of large-scale macroscopic structures and small-scale detailed structures of the medium, and finally obtain high-precision longitudinal and transverse wave velocity structures of strong scattering media. 4. The present invention can perform high-precision multi-parameter modeling of elastic strong scattering media in the absence of low-frequency information in seismic data and without model prior information.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a flow chart of the elastic wave direct envelope inversion method based on anisotropic total variation constraint;

Figure 2: True velocity model diagram;

(a) Model diagram of true longitudinal wave velocity, (b) model diagram of true shear wave velocity;

Fig. 3 Source wavelet and its spectrum;

(a) Source wavelet diagram (b) Source wavelet spectrum diagram;

Fig. 4 Initial velocity model diagram;

(a) Initial P-wave velocity model diagram, (b) Initial S-wave velocity model diagram;

Fig. 5. Conventional elastic wave full waveform inversion results;

(a) P-wave velocity result from conventional elastic wave full waveform inversion; (b) S-wave velocity result from conventional elastic wave full waveform inversion;

Fig. 6 Results of direct envelope inversion of unconstrained elastic waves;

(a) P-wave velocity results from direct envelope inversion of unconstrained elastic waves; (b) S-wave velocity results from direct envelope inversion of unconstrained elastic waves; (c) Final inversion results of P-wave velocity; (d) Final inversion results of S-wave velocity;

Fig. 7 Results of direct envelope inversion of elastic waves based on anisotropic total variation constraints;

(a) The results of P-wave velocity inversion using direct envelope inversion of elastic waves based on anisotropic total variation constraints; (b) The results of S-wave velocity inversion using direct envelope inversion of elastic waves based on anisotropic total variation constraints; (c) The final inversion results of P-wave velocity; (d) The final inversion results of S-wave velocity.

DETAILED DESCRIPTION

The present invention is further described in detail below with reference to the accompanying drawings and examples.

The elastic wave direct envelope inversion based on anisotropic total variation constraint of the present invention comprises the following steps:

a. Install the MATLAB software platform under Windows 7 or Linux system. MATLAB R2016a or above is required, and the Parallel Computing Toolbox must be installed.

b. Perform data preprocessing, static correction processing on the data, correct the influence of the undulating surface on the reflection phase axis; denoise the data to remove microseismic, low-frequency and high-frequency background noise and other random noise; remove interference waves, including sound waves, surface waves, industrial electrical interference, ghost reflections, multiple reflections, side waves, bottom waves, mixed reverberation and ringing, etc. Finally, high-quality elastic wave multi-component observation seismic data is obtained.

c. Perform wavelet estimation on seismic data. The estimation method can adopt direct wave estimation method and autocorrelation method to extract the source wavelet of each shot data.

d. Through background velocity analysis, the approximate background P-wave and S-wave velocity ranges are obtained, and a background velocity model is generated, which does not contain any prior information of strong scatterers, and serves as the initial P-wave velocity model v _p0 and initial S-wave velocity model v _s0 for inversion respectively.

上角标i表示i方向的分量，对于二维情况，指x(水平)和z(垂直)方向；对观测地震数据取包络得到观测包络数据

通过公式(7)计算弹性波直接包络反演的目标函数σ_EDEI：e. Calculate elastic wave multi-component simulated seismic data based on the initial P- and S-wave velocity models, and obtain simulated envelope data by taking the envelope of the simulated data

The superscript i indicates the component in the i direction. For the two-dimensional case, it refers to the x (horizontal) and z (vertical) directions. The observed seismic data is enveloping to obtain the observed envelope data.

The objective function σ _EDEI of elastic wave direct envelope inversion is calculated by formula (7):

Wherein, the subscript sr of the summation symbol indicates integration over all sources and receivers, t indicates time, and T indicates the total length of recording time.

和横波正传包络场

其中，波场模式分解采用公式(8)和(9)的方法：f. Calculate the simulated seismic wave field on the initial model, decompose the simulated seismic wave field to obtain the positive longitudinal wave field and the positive transverse wave field, and take the envelope to obtain the longitudinal wave positive envelope field

and the shear wave positive envelope field

The wave field mode decomposition adopts the method of formula (8) and (9):

Among them, u represents the elastic wave displacement vector, _up represents the longitudinal wave field, _us represents the shear wave field, ρ represents the density, λ and μ represent the Lame constants, ▽ represents the gradient operation, ▽· represents the divergence operation, and ▽× represents the curl operation.

g. Calculate the difference between the simulated envelope data and the observed envelope data to obtain the companion source f _s , that is,

和横波伴随包络场

波场模式分解方法见公式(8)和(9)。The companion source is transmitted back to obtain the companion envelope field; the companion envelope field is decomposed into wave field modes to obtain the longitudinal wave companion envelope field

and the shear wave envelope field

The wave field mode decomposition method is shown in formulas (8) and (9).

h. The P-wave velocity gradient of direct envelope inversion of elastic waves is obtained by performing zero-delay cross-correlation between the P-wave forward envelope field and the P-wave accompanying envelope field, as shown in formula (11):

Where _vp represents the longitudinal wave velocity and ρ represents the density. The shear wave velocity gradient of the elastic wave direct envelope inversion can be calculated by formula (12):

Where _vs represents the shear wave velocity, · represents the dot product, and μ represents the shear modulus.

和

表示，对二者施加各向异性全变分约束的过程等价于求解公式(13)和(14)所示的最优化问题：i. Select an appropriate step size and use the steepest descent method to update the P-wave velocity model and S-wave velocity model. Assume that the current iteration number is m, and the P-wave velocity model and S-wave velocity model obtained by the current iteration are respectively

and

It means that the process of applying anisotropic total variation constraints to the two is equivalent to solving the optimization problem shown in formulas (13) and (14):

和

分别为第m次迭代的纵、横波速度的更新量，λ₁和λ₂分别为对纵、横波速度施加各向异性全变分约束的权系数，||·||表示二范数，||·||_ATV表示各向异性全变分范数。对速度模型v计算各向异性全变分范数的具体表达式如公式(15)所示：Where J ₁ and J ₂ represent the objective functions for applying anisotropic total variation constraints to the P-wave and S-wave velocity models, α ₁ and α ₂ represent the update steps of the P-wave and S-wave velocities, respectively.

and

are the updated values of the P-wave and S-wave velocities at the mth iteration, λ ₁ and λ ₂ are the weight coefficients for applying anisotropic total variation constraints to the P-wave and S-wave velocities, ||·|| represents the second norm, and ||·|| _ATV represents the anisotropic total variation norm. The specific expression for calculating the anisotropic total variation norm for the velocity model v is shown in formula (15):

Where nz and nx represent the number of grid points in the vertical and horizontal directions of the model, respectively. Solve the optimization problems shown in formulas (13) and (14) to obtain the P- and S-wave velocity models after anisotropic total variation constraints as the velocity model after the current iteration update.

j. On the updated model, the iteration stop condition is judged; if the stop condition is not met, the updated longitudinal and transverse wave velocity models are used as the initial model, and the iterative calculation is continued in step e; if the stop condition is met, the output result is the large-scale longitudinal and transverse wave velocity structure v _pt and v _st of the strongly scattering medium after the anisotropic total variation constraint.

k. Using v _pt and v _st as the initial models, anisotropic total variation constrained elastic wave full waveform inversion is performed to obtain the final inversion result, namely the high-precision longitudinal and shear wave velocity structure of the strong scattering medium.

Example 1

The overall process of the present invention is shown in FIG1 .

Assume that the underground real P- and S-wave velocity models are shown in Figures 2a and 2b, respectively. In the real velocity model, the background velocity is low, and there is a high-speed, strongly scattering salt dome in the middle. Forward simulation on the real model can obtain the observed seismic record. The source wavelet waveform and spectrum are shown in Figures 3a and 3b, respectively. In order to simulate the lack of low-frequency information in actual seismic acquisition, the Ricker wavelet is high-pass filtered to cut off the low-frequency information below 3Hz. The main frequency of the source wavelet is about 9Hz. The background velocity analysis is performed using the observed records to obtain the approximate range of the background P- and S-wave velocities, and the initial P- and S-wave velocity models are established as shown in Figures 4a and 4b. The initial P- and S-wave velocity models do not contain prior information of strong scatterers.

In order to compare the inversion effect of the method of the present invention, the conventional elastic wave full waveform inversion is first performed on the initial P-wave and S-wave velocity models, and the P-wave and S-wave velocity inversion results are shown in Figures 5a and 5b, respectively. It can be seen that due to the lack of low-frequency information, the conventional elastic wave full waveform inversion can only obtain partial information on the top interface of the strong scatterer, and cannot restore the velocity information inside the strong scatterer and the morphological information of the strong scatterer.

Then, the unconstrained elastic wave direct envelope inversion method is performed, and the P-wave and S-wave velocity inversion results are shown in Figures 6a and 6b, respectively. It can be seen that although the morphology of the strong scatterer is restored to a certain extent, the velocity inside the salt dome is still deviated from the true value. Under the unconstrained condition, the results of Figures 6a and 6b are used as the initial model to perform conventional elastic wave full waveform inversion, and the final P-wave and S-wave velocity inversion results are shown in Figures 6c and 6d. It can be seen that the final result shows that the internal velocity uniformity of the salt dome is not good, and the lower boundary of the salt dome is not sufficiently characterized.

Then, the method of the present invention is used for inversion. Taking the model shown in Figure 4 as the initial model, according to steps (e) to (j), the P- and S-wave velocity models that can be obtained are shown in Figures 7a and 7b. It can be seen that the large-scale P- and S-wave velocity structure information of the strong scatterer has been successfully restored, and the velocity information in the salt dome is closer to the true value. Taking Figures 7a and 7b as the initial model, step (k) is executed, that is, anisotropic total variation constrained elastic wave full waveform inversion is performed, and the final P- and S-wave velocity inversion results are shown in Figures 7c and 7d. It can be seen that the boundary information and internal velocity of the strong scatterer are well inverted. The inversion effects of the interior, lower boundary and subsalt area of the salt dome are better than those in the unconstrained case. In short, the final inversion results of the method proposed in the present invention are more accurate for the internal velocity of the salt dome and the boundary of the salt dome, and its overall effect is significantly better than that of the conventional method (Figures 5 and 6).

Claims (1)

1. A method for elastic wave direct envelope inversion based on anisotropic total variation constraint, characterized in that, when low-frequency information of seismic data is missing and there is no model prior information, anisotropic total variation constraint is imposed on the elastic wave direct envelope inversion process to obtain high-quality large-scale P- and S-wave velocity structures of elastic strong scattering media; the obtained high-quality large-scale P- and S-wave velocity structures are used as the initial model of anisotropic total variation constraint elastic wave full waveform inversion to obtain high-precision P- and S-wave velocity structures of elastic strong scattering media; The following steps are involved: a. Install MATLAB software platform; b. Perform static correction and denoising preprocessing on the data to obtain high-quality elastic wave multi-component observation seismic data; c. Perform wavelet estimation on seismic data and extract the source wavelet of each shot data; d. By analyzing the background velocity, the approximate range of the background P-wave and S-wave velocities is obtained, and a background velocity model is generated, which does not contain any prior information of strong scatterers, and is used as the initial P-wave velocity model v _p0 and the initial S-wave velocity model v _s0 for inversion respectively; e. Calculate elastic wave multi-component simulated seismic data based on the initial longitudinal and transverse wave velocity models, and obtain simulated envelope data by taking the envelope of the simulated seismic data

The superscript i indicates the component in the i direction. For the two-dimensional case, it refers to the x horizontal direction and the z vertical direction. The observed seismic data is enveloping to obtain the observed envelope data.

The objective function σ _EDEI of elastic wave direct envelope inversion is calculated by formula (1):

In the formula, the subscript sr of the summation symbol indicates the integration of all sources and receivers, t indicates time, and T indicates the total recording time length; f. Calculate the simulated seismic wave field on the initial model, decompose the simulated seismic wave field to obtain the positive longitudinal wave field and the positive transverse wave field, and take the envelope to obtain the longitudinal wave positive envelope field

and the shear wave positive envelope field

g. Calculate the difference between the simulated envelope data and the observed envelope data to obtain the companion source, and then transmit the companion source back to obtain the companion envelope field; perform wave field mode decomposition on the companion envelope field to obtain the longitudinal wave companion envelope field

and the shear wave envelope field

h. The P-wave velocity gradient of direct envelope inversion of elastic waves is obtained by performing zero-delay cross-correlation between the P-wave forward envelope field and the P-wave accompanying envelope field, as shown in formula (2):

Where _vp represents the P-wave velocity, ρ represents the density, and the S-wave velocity gradient of the elastic wave direct envelope inversion can be calculated by formula (3):

In the formula, _vs represents the shear wave velocity, · represents the dot product, and μ represents the shear modulus; i. Select an appropriate step size and use the steepest descent method to update the P-wave and S-wave velocity models. Assuming that the current number of iterations is m, the P-wave velocity and S-wave velocity models obtained by the current iteration are respectively expressed as

and

It means that the process of applying anisotropic total variation constraints to the two is equivalent to solving the optimization problem shown in formulas (4) and (5):

Where J ₁ and J ₂ represent the objective functions for applying anisotropic total variation constraints to the P-wave and S-wave velocity models, α ₁ and α ₂ represent the update steps of the P-wave and S-wave velocities, respectively.

and

are the updated values of the P-wave and S-wave velocities at the mth iteration, _λ1 and _λ2 are the weight coefficients for applying anisotropic total variation constraints on the P-wave and S-wave velocities, ||·|| represents the binary norm, and ||·|| _ATV represents the anisotropic total variation norm. The specific expression for calculating the anisotropic total variation norm for the velocity model v is shown in formula (6):

Where nz and nx represent the number of grid points in the vertical and horizontal directions of the model, respectively. Solve the optimization problems shown in formulas (4) and (5) to obtain the P- and S-wave velocity models after anisotropic total variation constraints as the velocity model after the current iteration update. j. On the updated model, the iteration stop condition is judged; if the stop condition is not met, the updated P- and S-wave velocity models are used as the initial model, and the iterative calculation is continued in step e; if the stop condition is met, the output result is the large-scale P- and S-wave velocity structure v _pt and v _st of the strongly scattering medium after the anisotropic total variation constraint; k. Using v _pt and v _st as the initial models, anisotropic total variation constrained elastic wave full waveform inversion is performed to obtain the final inversion result, namely the high-precision longitudinal and shear wave velocity structure of the strong scattering medium.

Images