CN108845351A - A kind of VSP seismic data converted wave full waveform inversion method - Google Patents

Abstract

The invention discloses a full waveform inversion method of converted wave of VSP seismic data. Separation of upgoing PS converted waves by global optimization parameter inversion method; calculation of shear wave velocity according to empirical formula and smoothing to obtain initial shear wave velocity model; solution of elastic wave equation for separation of longitudinal and transverse waves to obtain forward longitudinal and transverse waves at each time step Wave field and simulated VSP PS wave records; by solving the adjoint elastic wave equation of the separation of longitudinal and transverse waves, the reverse longitudinal and transverse wave field at each time step is obtained; the forward longitudinal wave field and the backward transverse wave field are used to calculate the shear wave velocity gradient; The three-point parabola fitting method is used to obtain the iteration step size; the shear wave velocity model is updated by the L-BFGS method; the Wiener filter is used to perform multi-scale inversion from low frequency to high frequency. The beneficial effect of the invention is to improve the accuracy of full waveform inversion of converted wave of VSP seismic data through the separation strategy of longitudinal and transverse waves, and improve the inversion effect of complex structures.

Description

A Converted Wave Full Waveform Inversion Method for VSP Seismic Data

technical field

The invention belongs to the technical field of exploration geophysics, and relates to a new full waveform inversion method, which can be used for VSP seismic data inversion to obtain a high-precision converted wave velocity model.

Background technique

With the deepening of oil and gas exploration and development and the continuous improvement of borehole seismic technology, the research and application of VSP (vertical seismic profiling) exploration technology has also been developed rapidly. In areas with complex structures or complex media, VSP technology can effectively understand the geological structure, petrophysical properties and other characteristics near the wellbore. Through VSP imaging combined with ground three-dimensional seismic data and lithological data, fine seismic attribute analysis can be carried out, and key layers can be further verified, the contact relationship between rock layers can be clarified, an accurate geological model can be established, and stratum heterogeneity Sexual research, providing technical support for the adjustment of the development plan. Compared with surface seismic data, VSP data contain richer wave field information, high resolution, high signal-to-noise ratio, and the seismic records received in the well only pass through a low deceleration zone once, so the attenuation of seismic wave energy compared with surface data Smaller, so it has its unique advantages for deep imaging. Secondly, due to observation limitations in ground seismic exploration, it is difficult to observe the salt dome body and its flanks and sub-salt information in the target layer. VSP technology can make up for these deficiencies in ground data. Therefore, no matter in the field of exploration or development, VSP technology has its unique advantages and functions.

Velocity is the core parameter in the data processing flow of seismic data. High-quality velocity modeling methods are decisive for subsequent migration imaging, lithology, lithofacies interpretation, and reservoir prediction. So far, the data processing flow of VSP P-wave data has been relatively perfect, and the inversion accuracy of P-wave velocity is still acceptable. However, due to the influence of many factors (such as the extreme asymmetry of shear wave propagation, the coupling crosstalk of longitudinal and transverse waves, and the anisotropy of shear waves, etc.), the accuracy of converted wave velocity inversion is still low. Elastic full waveform inversion (EFWI) has become increasingly popular in exploration geophysics for its potential to obtain accurate P-wave velocities, S-wave velocities, and densities. However, due to multi-parameter coupling, elastic wave full waveform inversion has serious nonlinearity and strong multi-solution, so it has not been widely promoted in industry. In order to eliminate the influence of P-wave velocity on the accuracy of S-wave velocity inversion, it is necessary to further study the converted wave full waveform inversion method.

The invention provides a new method for full waveform inversion of converted wave of VSP seismic data. For the VSP observation mode (surface blasting, well reception), the PS wave (converted wave) record is separated by the global optimization parameter inversion method, and the separated converted wave is used for inversion to alleviate nonlinearity and multi-solution; by solving the longitudinal and transverse waves Separate the elastic wave equation and the adjoint equation to obtain the longitudinal and transverse wave fields at different times, and use the forward longitudinal wave field and the accompanying shear wave field to calculate the gradient of the objective function with respect to the shear wave velocity, weaken the crosstalk effect between different wave modes, and increase the shear wave velocity inversion accuracy.

Contents of the invention

The purpose of the present invention is to provide a full waveform inversion method of converted wave of VSP seismic data.

The technical scheme adopted in the present invention is to carry out according to the following steps:

(1) PS wave (converted wave) separation of original VSP (vertical seismic profile) seismic data;

Separation of uplink PS converted waves by global optimization parameter inversion method;

(2) Establishment of the initial model of shear wave velocity;

Input the longitudinal wave velocity model, obtain the shear wave velocity model according to the empirical formula, and smooth it to obtain the initial shear wave velocity model;

(3) The source wave field is extended forward to separate the longitudinal and transverse wave fields;

The forward P-S wave field and the simulated VSP PS-wave record are obtained at each time step by solving the elastic wave equation for the separation of P- and S-waves.

(4) Reverse continuation of the converted wave residual to separate the longitudinal and transverse wave fields;

By solving the adjoint elastic wave equations of the separation of the longitudinal and transverse waves, the backpropagating longitudinal and transverse wave fields are obtained at each time step;

(5) Calculate the shear wave velocity gradient.

Calculate the shear wave velocity gradient using the forward longitudinal wave field and the reverse propagating shear wave field obtained in steps (3) and (4);

(6) Find the iteration step size;

Use the three-point parabola fitting method to find the iteration step size;

(7) Update the shear wave velocity model;

When the convergence condition is not met, repeat step (3) to step (6), when the convergence condition is met, the iteration ends;

(8) Multi-scale inversion;

Wiener filtering is used to divide the VSP PS converted wave records, and the inversion is carried out step by step from low frequency to high frequency.

Specifically, the PS converted wave separation method of the original VSP seismic data in step (1) of the converted wave full waveform inversion of VSP seismic data is as follows:

The parameter inversion method is a method for separating wave fields of VSP data based on seismic wave apparent velocity and polarization characteristics. In the present invention, the global optimized parameter inversion method is used to separate PS converted waves.

For each depth point, the minimized target functional of observed data and simulated records is established:

Among them, Ψ _Pdown , Ψ _Pup , Ψ _Sdown and Ψ _Sup are the incident angles of downgoing P-wave, upgoing P-wave, downgoing S-wave and upgoing S-wave, respectively; v _P and v _S are the P- and S-wave velocities; u and u _obs are the simulated and observed wave fields , ω is the angular frequency, [ω ₁ ,ω ₁ ] is the frequency range.

The optimal propagation speed and propagation angle can be obtained through the iteration of the global optimization genetic algorithm. Solving equation (2) realizes the separation of downgoing longitudinal wave, upgoing longitudinal wave, downgoing shear wave and upgoing shear wave.

u(ω)=[u(-M ₁ ,ω)...u(0,ω)...u(M ₂ ,ω)] ^T

＝a _Pdown u _Pdown (ω)+a _Pup u _Pup (ω)+a _Sdown u _Sdown (ω)+a _Sup u _Sup (ω), (2)

Among them, [u _Pdown , u _Pup , u _Sdown , u _Sup ] are the wave fields of downgoing longitudinal wave, upgoing longitudinal wave, downgoing shear wave and upgoing shear wave respectively, and [a _Pdown , a _Pup , a _Sdown , a _up ] are the corresponding polarization vectors, [-M ₁ ,…,0,…,M ₂ ] are the numbers of the receiver points in the spatial window.

When the seismic source is a longitudinal wave source, the separated upward shear wave u _Sup (frequency domain) is inversely Fourier transformed to obtain the observed PS converted wave record (decomposed into horizontal and vertical components according to the propagation angle).

The establishment method of the initial model of shear wave velocity in step (2) is as follows:

The empirical formula for compressional and shear wave velocity:

According to the input P-wave velocity model, the empirical formula is used to obtain the shear-wave velocity model, which is smoothed by Gaussian window to obtain the initial shear-wave velocity model.

In step (3), the source wave field is extended forward, and the method of separating the longitudinal and transverse wave fields is as follows:

The first-order elastic wave velocity-stress equation:

Among them, λ and μ are Lame constants, ρ is density, (v _x , v _z ) is particle vibration velocity vector, (τ _xx ,τ _zz ,τ _xz ) is stress vector.

Adopt the velocity-stress elastic wave equation of longitudinal and transverse wave separation in the present invention:

Among them, (v _xp , v _zp ) is the particle vibration velocity vector of longitudinal wave, (v _xs , v _zs ) is the particle vibration velocity vector of shear wave, (τ _xxp ,τ _zzp ) is the stress vector of longitudinal wave, (τ _xxs ,τ _zzs ,τ _xzs ) is the shear wave stress vector.

The high-order staggered grid finite difference method is used to numerically solve Equation (6) to obtain the forward P- and S-wave field and simulated VSP PS wave records at each time step.

Step (4) Reverse continuation of the converted wave residual, and the method of separating the longitudinal and transverse wave fields is as follows:

In the adjoint state method, the adjoint operator satisfies:

<Lu,u ^* >＝<u,L ^* u ^* >, (7)

Among them: L is the forward operator, u is the particle velocity and stress vector, "*" is its adjoint state. Substituting Equation (4) into Equation (7) yields the conventional adjoint equation:

Similarly, based on equation (6), the adjoint equation for the separation of P- and S-waves can be obtained:

Among them, "*" indicates the accompanying wave field of the corresponding wave field, For the observation and simulation of VSP seismic records coupled with compressional and shear waves, Separated observations and simulated VSP converted wave seismic records.

Taking the simulated and observed VSP PS wave recording residuals as the seismic source, equations (10) and (11) are numerically solved using high-order staggered finite differences to obtain the backpropagating P-S wave field at each time step.

Step (5) The method of calculating the shear wave velocity gradient is as follows:

The shear wave velocity gradient in traditional full waveform inversion is:

Among them, E is the target functional.

The present invention calculates the shear wave velocity gradient by using the forward longitudinal wave field and the backpropagation shear wave field:

in,

By adopting the new gradient formula (Equation 13), the interference of P-wave velocity can be alleviated, and the inversion accuracy of S-wave velocity can be improved.

Step (6) The method of calculating the iteration step size is as follows:

In the present invention, a three-point parabolic interpolation method is used to calculate the step size. Select three trial points, and use these three points to construct a parabola. The extreme point of this parabola is the optimal step size.

The target functionals corresponding to the three trial steps a ₀ , a ₁ , and a ₂ are E _s0 , E _s1 , E _s2 , and the optimal step size is:

Step (7) The method of updating the shear wave velocity model is as follows:

Update the shear wave velocity by the L-BFGS method:

Among them, m _k and m _k+1 are the shear wave velocity models of the current iteration and the next iteration, respectively, and H _k is the approximate matrix of the inverse of the Hessian matrix.

The multi-scale inversion method in step 8 is as follows:

Wiener filtering is used to divide the VSP PS converted wave records, and the inversion is carried out step by step from low frequency to high frequency. The Wiener filter is:

Among them, f _wiener is the wiener filter, W _original is the original main frequency wavelet, W _target target main frequency wavelet, ω is the angular frequency, ε is the stability coefficient, and * indicates the complex conjugate.

The large-scale inversion result is used as the initial model of the adjacent small-scale inversion, and steps (3) to (7) are repeated until the convergence condition is satisfied.

The beneficial effect of the invention is to improve the precision of full waveform inversion of converted wave of VSP seismic data through the separation strategy of longitudinal and transverse waves. The purpose of the invention is to make up for the shortcomings of the lighting limitation of the ground observation system, improve the modeling accuracy of the shear wave velocity, especially the inversion effect of the flank of the salt dome and the target layer with a steep dip angle.

Description of drawings

Fig. 1 Flowchart of converted wave full waveform inversion method of VSP seismic data;

Figure 2 Complex structure model;

Fig. 3 The shear wave velocity gradient of the first iteration at the first scale;

Figure 4 S-wave velocity inversion results;

Fig. 5 Comparison diagram of shear wave velocity curves of the longitudinal sideline beside the well;

Fig. 6 Convergence curve of converted shear wave target functional.

Detailed ways

The present invention will be described in detail below in combination with specific embodiments.

As shown in Figure 1, it is a flow chart of implementing VSP seismic data converted wave full waveform inversion for the present invention, specifically comprising:

(1) PS converted wave separation of original VSP seismic data;

Separation of uplink PS converted waves by global optimization parameter inversion method;

(2) Establishment of the initial model of shear wave velocity;

Input the longitudinal wave velocity model, obtain the shear wave velocity model according to the empirical formula, and smooth it to obtain the initial shear wave velocity model;

(3) The source wave field is extended forward to separate the longitudinal and transverse wave fields;

The forward P-S wave field and the simulated VSP PS-wave record are obtained at each time step by solving the elastic wave equation for the separation of P- and S-waves.

(4) Reverse continuation of the converted wave residual to separate the longitudinal and transverse wave fields;

By solving the adjoint elastic wave equations of the separation of the longitudinal and transverse waves, the backpropagating longitudinal and transverse wave fields are obtained at each time step;

(5) Calculate the shear wave velocity gradient.

Calculate the shear wave velocity gradient using the forward longitudinal wave field and the reverse propagating shear wave field obtained in steps (3) and (4);

(6) Find the iteration step size;

Use the three-point parabola fitting method to find the iteration step size;

(7) Update the shear wave velocity model;

When the convergence condition is not met, repeat step (3) to step (6), when the convergence condition is met, the iteration ends;

(8) Multi-scale inversion;

Wiener filtering is used to divide the VSP PS converted wave records, and the inversion is carried out step by step from low frequency to high frequency.

In step (1), the PS converted wave separation method of the original VSP seismic data is as follows:

In the present invention, a global optimization parameter inversion method is used to separate PS converted waves.

For each depth point, the minimized target functional of observed data and simulated records is established:

Among them, Ψ _Pdown , Ψ _Pup , Ψ _Sdown and Ψ _Sup are the incident angles of downgoing P-wave, upgoing P-wave, downgoing S-wave and upgoing S-wave, respectively; v _P and v _S are the P- and S-wave velocities; u and u _obs are the simulated and observed wave fields , ω is the angular frequency, [ω ₁ ,ω ₁ ] is the frequency range.

The optimal propagation speed and propagation angle can be obtained through the iteration of the global optimization genetic algorithm. Solving equation (2) realizes the separation of downgoing longitudinal wave, upgoing longitudinal wave, downgoing shear wave and upgoing shear wave.

u(ω)=[u(-M ₁ ,ω)...u(0,ω)...u(M ₂ ,ω)] ^T

＝a _Pdown u _Pdown (ω)+a _Pup u _Pup (ω)+a _Sdown u _Sdown (ω)+a _Sup u _Sup (ω), (2)

Among them, [u _Pdown , u _Pup , u _{Sdown ,} u _Sup ] are respectively the wave fields of downgoing longitudinal wave, upgoing longitudinal wave, downgoing shear wave and upgoing shear wave, and [a _Pdown , a _Pup , a _Sdown , a _up ] are corresponding polarization vectors, [-M ₁ ,…,0,…,M ₂ ] are the numbers of the receiver points in the spatial window.

When the seismic source is a longitudinal wave source, the separated upward shear wave u _Sup (frequency domain) is inversely Fourier transformed to obtain the observed PS converted wave record (decomposed into horizontal and vertical components according to the propagation angle).

The establishment method of the initial model of shear wave velocity in step (2) is as follows:

The empirical formula for compressional and shear wave velocity:

According to the input P-wave velocity model, the empirical formula is used to obtain the shear-wave velocity model, which is smoothed by Gaussian window to obtain the initial shear-wave velocity model.

In step (3), the source wave field is extended forward, and the method of separating the longitudinal and transverse wave fields is as follows:

The first-order elastic wave velocity-stress equation:

Among them, λ and μ are Lame constants, ρ is density, (v _x , v _z ) is particle vibration velocity vector, (τ _xx ,τ _zz ,τ _xz ) is stress vector.

Adopt the velocity-stress elastic wave equation of longitudinal and transverse wave separation in the present invention:

Among them, (v _xp , v _zp ) is the particle vibration velocity vector of longitudinal wave, (v _xs , v _zs ) is the particle vibration velocity vector of shear wave, (τ _xxp ,τ _zzp ) is the stress vector of longitudinal wave, (τ _xxs ,τ _zzs ,τ _xzs ) is the shear wave stress vector.

The high-order staggered grid finite difference method is used to numerically solve Equation (6) to obtain the forward P- and S-wave field and simulated VSP PS wave records at each time step.

Step (4) Reverse continuation of the converted wave residual, and the method of separating the longitudinal and transverse wave fields is as follows:

In the adjoint state method, the adjoint operator satisfies:

<Lu,u ^* >＝<u,L ^* u ^* >, (7)

Among them: L is the forward operator, u is the particle velocity and stress vector, "*" is its adjoint state. Substituting Equation (4) into Equation (7) yields the conventional adjoint equation:

Similarly, based on equation (6), the adjoint equation for the separation of P- and S-waves can be obtained:

Among them, "*" indicates the accompanying wave field of the corresponding wave field, For the observation and simulation of VSP seismic records coupled with compressional and shear waves, Separated observations and simulated VSP converted wave seismic records.

Taking the simulated and observed VSP PS wave recording residuals as the seismic source, equations (10) and (11) are numerically solved using high-order staggered finite differences to obtain the backpropagating P-S wave field at each time step.

Step (5) The method of calculating the shear wave velocity gradient is as follows:

The shear wave velocity gradient in traditional full waveform inversion is:

Among them, E is the target functional.

The present invention calculates the shear wave velocity gradient by using the forward longitudinal wave field and the backpropagation shear wave field:

in,

By adopting the new gradient formula (Equation 13), the interference of P-wave velocity can be alleviated, and the inversion accuracy of S-wave velocity can be improved.

Step (6) The method of calculating the iteration step size is as follows:

In the present invention, a three-point parabolic interpolation method is used to calculate the step size. Select three trial points, and use these three points to construct a parabola. The extreme point of this parabola is the optimal step size.

The target functionals corresponding to the three trial steps a ₀ , a ₁ , and a ₂ are E _s0 , E _s1 , E _s2 , and the optimal step size is:

Step (7) The method of updating the shear wave velocity model is as follows:

Update the shear wave velocity by the L-BFGS method:

Among them, m _k and m _k+1 are the shear wave velocity models of the current iteration and the next iteration, respectively, and H _k is the approximate matrix of the inverse of the Hessian matrix.

The multi-scale inversion method in step (8) is as follows:

Wiener filtering is used to divide the VSP PS converted wave records, and the inversion is carried out step by step from low frequency to high frequency. The Wiener filter is:

Among them, f _wiener is the wiener filter, W _original is the original main frequency wavelet, W _target target main frequency wavelet, ω is the angular frequency, ε is the stability coefficient, and * indicates the complex conjugate.

The large-scale inversion result is used as the initial model of the adjacent small-scale inversion, and steps (3) to (7) are repeated until the convergence condition is met.

Due to the adoption of the above technical solutions, the present invention has the following advantages: 1. The modeling accuracy of the shear wave velocity model can be improved. 2. It can improve the inversion accuracy of complex structures (such as salt dome flanks and steep dip angle target layers).

Next, an example is used to analyze the effectiveness of the converted wave full waveform inversion method of VSP seismic data proposed in the present invention.

(a)-(c) in Fig. 2 are the true P-wave velocity, true S-wave velocity and initial S-wave velocity model of a complex structural model, respectively. The calculation area is 200×200 grid points, the spatial interval is 10m, the time step is 1ms, and the total recording time is 1.7s. The source is the 15Hz Reich wavelet, which is added to the normal stress (longitudinal wave source). There is a high-speed body at the lower right of the model, and the shape of the longitudinal wave velocity and the shear wave velocity are different at the high-speed body. The well is located at x=1980m, and 180 geophones are arranged in the well, one every 10 meters from z=100m. The shot point is located on the surface, and a shot is fired every 70m from x=150m, a total of 26 shots. The second-order time and tenth-order precision finite difference method is used for numerical simulation, and the PML absorbs the boundary to eliminate the boundary reflection. Two scales are used for multi-scale inversion, the first scale is 15Hz and the second scale is 30Hz. Figure 3 shows the shear wave velocity gradient for the first iteration of the first scale. It can be seen that the shallow structure beside the well is clear, and the deep structure and the flank of the salt dome are also accurately depicted. (a) and (b) in Fig. 4 are the inversion results of the first-scale and second-scale shear wave velocities, respectively. It can be seen from the figure that the VSP seismic data converted wave full waveform inversion method proposed by the present invention has high inversion accuracy, the wellside and salt dome flank structures are clearly visible, and the high-speed anomaly body is also well recovered. Fig. 5 is a comparative graph of true shear wave velocity, initial velocity and inversion results at x=1650m. The inversion results are in good agreement with the real results. Fig. 6 shows the convergence curve of the objective function Es, which further shows the validity of the converted wave full waveform inversion method for VSP seismic data of the present invention.

The invention is a new converted wave full waveform inversion method of VSP seismic data, which can obtain a high-precision shear wave velocity model; and can greatly improve the inversion accuracy of complex geological structures (such as salt dome flanks and steep interfaces).

The above description is only a preferred embodiment of the present invention, and does not limit the present invention in any form. Any simple modifications made to the above embodiments according to the technical essence of the present invention, equivalent changes and modifications, all belong to this invention. within the scope of the technical solution of the invention.

Claims (9)

1. A full waveform inversion method for VSP seismic data converted waves is characterized by comprising the following steps:

(1) separating PS converted waves of original VSP seismic data;

separating the uplink PS converted wave by a global optimization parameter inversion method;

(2) establishing a transverse wave velocity initial model;

inputting a longitudinal wave velocity model to obtain a transverse wave velocity model, and smoothing the transverse wave velocity model to obtain an initial transverse wave velocity model;

(3) carrying out forward continuation on a seismic source wave field, and separating a longitudinal and transverse wave field;

obtaining a positive transmission longitudinal and transverse wave field of each time step and a simulated VSPPS wave record by solving an elastic wave equation of longitudinal and transverse wave separation;

(4) converting wave residual error reverse continuation, and separating longitudinal and transverse wave fields;

obtaining a reverse transmission longitudinal and transverse wave field of each time step by solving an adjoint elastic wave equation of longitudinal and transverse wave separation;

(5) calculating the velocity gradient of the transverse wave;

calculating the shear wave velocity gradient by adopting the forward-transmission longitudinal wave field obtained in the step (3) and the backward-transmission shear wave field obtained in the step (4);

(6) calculating an iteration step length;

solving an iteration step length by adopting a three-point parabola fitting method;

(7) updating the transverse wave velocity model;

when the convergence condition is not met, repeating the steps (3) to (6), and when the convergence condition is met, finishing the iteration;

(8) performing multi-scale inversion;

and performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency.

2. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the PS conversion wave separation method of the original VSP seismic data in the step (1) is as follows:

separating PS converted waves by adopting a global optimization parameter inversion method, and establishing a minimized target functional of observation data and simulation record for each depth point:

therein, Ψ_Pdown,Ψ_Pup,Ψ_SdownAnd Ψ_SupIncident angles v of down longitudinal wave, up longitudinal wave, down transverse wave and up transverse wave respectively_PAnd v_SIs the longitudinal and transverse wave velocity, u and u_obsFor simulating and observing the wavefield, ω is the angular frequency, [ ω ]₁,ω₁]For the frequency range, the optimal propagation speed and propagation angle can be obtained through the iteration of a global optimization genetic algorithm; solving an equation (2) to realize separation of downlink longitudinal waves, uplink longitudinal waves, downlink transverse waves and uplink transverse waves;

u(ω)＝[u(-M₁,ω)…u(0,ω)…u(M₂,ω)]^T

＝a_Pdownu_Pdown(ω)+a_Pupu_Pup(ω)+a_Sdownu_Sdown(ω)+a_Supu_Sup(ω), (2)

wherein [ u ]_Pdown,u_Pup,u_Sdown,u_Sup]Respectively, a down-going longitudinal wave, an up-going longitudinal wave, a down-going transverse wave and an up-going transverse wave, [ a ]_Pdown,a_Pup,a_Sdown,a_up]For the corresponding polarization vector, [ -M₁,…,0,…,M₂]The number of the detection point in the space window is the number, when the seismic source is a longitudinal wave source, the separated up-going transverse wave u_SupAnd performing Fourier inversion to obtain an observation PS converted wave record.

3. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for establishing the shear wave velocity initial model in the step (2) is as follows:

longitudinal and transverse wave velocity empirical formula:

and obtaining a transverse wave velocity model by adopting an empirical formula according to the input longitudinal wave velocity model, and smoothing the transverse wave velocity model through a Gaussian window to obtain an initial transverse wave velocity model.

4. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for separating the longitudinal wave field from the transverse wave field in the step (3) comprises the following steps:

first order elastic wave velocity-stress equation:

wherein λ and μ are Lame constants, ρ is density, (v)_x,v_z) Is the particle vibration velocity vector (tau)_xx,τ_zz,τ_xz) Is a stress vector;

adopting a speed-stress elastic wave equation of longitudinal and transverse wave separation:

wherein (v)_xp,v_zp) Is the vibration velocity vector of longitudinal wave particle, (v)_xs,v_zs) Is the transverse wave particle vibration velocity vector (tau)_xxp,τ_zzp) Is a longitudinal wave stress vector (tau)_xxs,τ_zzs,τ_xzs) Is the transverse wave stress vector;

and (3) solving the equation (6) by adopting a high-order staggered grid finite difference numerical value to obtain a positive transmission longitudinal and transverse wave field of each time step and simulating VSP PS wave record.

5. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the converted wave residual error reverse continuation in the step (4) and the method for separating the longitudinal wave field and the transverse wave field are as follows:

in the adjoint state method, the adjoint operator satisfies:

<Lu,u^*>＝<u,L^*u^*>, (7)

wherein: l is the positive operator, u is the particle velocity and stress vector, "# is its companion state, substituting equation (4) into equation (7) yields the conventional companion equation:

obtaining a concomitant equation of longitudinal and transverse wave separation based on equation (6):

wherein, represents the adjoint wave field of the corresponding wave field,for observation of the coupling of longitudinal and transverse waves and for simulation of VSP seismic recordings,for separated observation and simulation VSP converted wave seismic records, simulating and observing VSP PS wave record residual errors are used as seismic sources, and equations (10) and (11) are solved by adopting high-order interleaved finite difference numerical values to obtain a reverse propagation longitudinal and transverse wave field of each time step.

6. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for calculating the shear wave velocity gradient in the step (5) comprises the following steps:

the shear wave velocity gradient in the conventional full waveform inversion is as follows:

wherein E is a target functional, and the shear wave velocity gradient is calculated by adopting a forward-transmission longitudinal wave field and a backward-transmission shear wave field:

wherein,

the new gradient formula is adopted to relieve the interference of longitudinal wave velocity, and further the inversion accuracy of transverse wave velocity is improved.

7. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the iterative step length calculating method in the step (6) is as follows:

selecting three trial points, and constructing a parabola by using the three points, wherein the extreme point of the parabola is the optimal step length, and the three trial step lengths a₀,a₁,a₂Corresponding target functional is E_s0,E_s1,E_s2Then the optimal step length is:

8. the VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the method for updating the shear wave velocity model in the step (7) is as follows:

updating the transverse wave velocity by an L-BFGS method:

wherein m is_kAnd m_k+1Shear wave velocity models, H, for the current iteration and the next iteration, respectively_kAn approximation matrix that is the inverse of the hessian matrix.

9. The VSP seismic data converted waveform full waveform inversion method of claim 1, wherein: the multi-scale inversion method in the step (8) is as follows:

performing frequency division processing on the VSP PS converted wave record by using wiener filtering, and gradually performing inversion from low frequency to high frequency, wherein the wiener filter is as follows:

wherein f is_wienerIs a wiener filter, W_originalIs a primary dominant frequency wavelet, W_targetAnd (3) repeating the steps (3) to (7) until a convergence condition is met, wherein omega is angular frequency, epsilon is a stable coefficient, represents complex conjugate, and the inversion result with a large scale is used as an initial model of adjacent small-scale inversion.