A method for constructing constrained attenuation compensation velocity model of deep water shallow seismic data

Abstract

The application discloses a deep-water shallow seismic data structure-constrained attenuation compensation velocity modeling method, which comprises the following steps: obtaining an initial velocity model by using a tomographic inversion method and performing reverse-time migration imaging to obtain an initial imaging result; and obtaining a seismic data dip angle field by using a dip angle prediction technology; a structure-constrained velocity modeling target functional is established by using the dip angle field and multi-frequency band seismic data; wave field continuation for Q compensation is performed by using a decoupled viscous wave equation, and the gradient of the functional is calculated to update the velocity model; a multi-scale multi-frequency band inversion strategy is used to solve a high-precision velocity inversion optimization problem to obtain a high-precision deep-water shallow velocity field. By applying attenuation compensation to full waveform inversion, a structure guide constraint is introduced into a full waveform inversion target function, and a velocity parameter model is constrained by using seismic imaging structure information, so that the precision and efficiency of inversion can be effectively improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of deep water shallow seismic exploration, and more particularly to a deep water shallow seismic data structure constraint attenuation compensation velocity modeling method. BACKGROUND

[0002] With the increasing level and connotation of geosciences research and petroleum exploration, the demand for seismic exploration is gradually increasing, and the precision of seismic exploration is continuously improving. In order to meet the needs of dynamic prediction of oil and gas reservoirs, lithology analysis and structure analysis, higher quality seismic data processing is carried out. Seismic imaging technology can intuitively and clearly display the properties of underground material structure in the form of images. The iterative inversion of the initial model is the basis of seismic imaging technology, and the velocity modeling is a key step in the process of seismic imaging processing. The quality of the velocity model built will directly affect the final results of seismic data processing, that is, the quality of the seismic profile. That is, obtaining an accurate velocity model helps to improve the precision of seismic imaging. However, due to the strong nonlinearity of full waveform inversion (FWI), the solving process of full waveform inversion requires a large amount of computing cost. In addition, the actual complex medium of stratum has an absorption effect on seismic waves, which causes changes in amplitude and frequency of seismic waves during propagation. Therefore, when calculating the gradient of full waveform inversion, energy attenuation and phase distortion will be caused, which leads to the weakening of deep gradient energy and the reduction of deep inversion precision. Therefore, additional iterations are usually required to make up for this deficiency.

[0003] The existing velocity model establishment method mainly includes the following steps: 1. determining the objective function of the seismic data velocity model; 2. calculating the gradient of the objective function based on the horizon constraint, that is, increasing the constraint condition by introducing the horizon constraint in the process of calculating the gradient of the objective function; 3. adjusting the objective function according to the calculated gradient to obtain an adjusted velocity model. However, this velocity modeling method does not fully consider the absorption effect of actual stratum medium on seismic waves, so the inversion result will cause energy attenuation and phase distortion, affecting the inversion precision.

[0004] Full waveform inversion is a strong nonlinear and initial model dependent method, and other prior constraints can effectively reduce the ill-conditioned inversion and avoid the cycle skipping phenomenon of inversion. However, the existing velocity modeling method does not fully apply the concept of attenuation compensation, but simply increases the number of iterations to improve the inversion precision, which wastes a lot of computing power. SUMMARY

[0005] The purpose of this invention is to overcome the shortcomings of existing velocity modeling methods in terms of poor inversion accuracy, and to provide a structural constraint attenuation compensation velocity modeling method for deep-water shallow seismic data. Attenuation compensation is applied to the full waveform inversion to improve the inversion convergence rate, which can quickly obtain high-precision shallow velocity information, providing a foundation for subsequent seismic imaging and inversion.

[0006] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0007] This invention provides a method for modeling the velocity of structurally constrained attenuation compensation in deep-water shallow seismic data, comprising the following steps:

[0008] Step S1: Obtain the initial velocity model using the tomographic inversion method and perform reverse time migration imaging to obtain the initial imaging results; use dip prediction technology to obtain the dip field of the seismic data.

[0009] Step S2: Establish the objective functional for tectonic constraint velocity modeling using the dip field and multi-band seismic data;

[0010] Step S3: Perform Q-compensated wavefield extension using the decoupled viscous acoustic wave equation, calculate the gradient of the functional, and update the velocity model.

[0011] Step S4: Solve the high-precision velocity inversion optimization problem using a multi-scale, multi-band inversion strategy to obtain a high-precision deep-water shallow-layer velocity field.

[0012] Further, in step S1, firstly, migration imaging is performed using the initial velocity model to obtain initial imaging results, and then the initial imaging results are used to establish the mathematical relationship between travel time and tilt field:

[0013]

[0014] Obtain the dip field information of the initial earthquake model, where The zero offset travel time curve is represented by t, where t represents the travel time, t0 represents the first travel time, x represents the offset, x0 represents the zero offset, and p represents the tilt field.

[0015] Furthermore, in step S2, a constraint term on the inversion model is introduced into the objective function of the velocity modeling inversion, effectively improving the stability and convergence rate of the inverse problem. The dip field information and adaptive total variation constraints are introduced into the objective function to establish a constrained attenuation-compensated velocity modeling objective functional.

[0016]

[0017] Among them, J Q It is the objective function, u(v) c Q) represents the simulation data considering attenuation, f(v) rQ) represents the observed data with attenuation, v c represents the current velocity model, v r represents the final velocity model, is the structural constraint regularization term, and are the gradient operators parallel and perpendicular to the structural direction of the velocity model, respectively.

[0018] Further, the structural constraint regularization term is λ1, λ2, respectively, the structural constraint regularization weight factors corresponding to and The expressions of the gradient operators and can be extended as follows:

[0019]

[0020] wherein ε1 and ε2 are the weighting coefficients corresponding to and and and are the gradient operators along the horizontal and vertical directions, respectively, and v(x, z) represents the velocity model; using the dip prediction technique, the local dip information corresponding to each point of the velocity model is obtained, and the local dip information is combined with the horizontal and vertical gradient operators, so that the gradient constraint operators parallel and perpendicular to the structural direction are obtained. Compared with the horizontal and vertical gradient operators commonly used, the new structural constraint operators and can more effectively protect the boundary information of the inversion model, and only a small amount of calculation is required relative to the full waveform inversion; after detecting the boundary information of the model, the obtained boundary information is used to make certain modifications to the model update gradient, so that the model gradient is more sensitive to the interface.

[0021] In addition, the parameters ε1 and ε2 are set as parameters related to the dip p, i.e. ε1(p) and ε2(p), the larger the dip p, the smaller the ε1(p), and the larger the ε2(p), and vice versa. Such constraints ensure that the inversion constraint process can be automatically adjusted in the case of large dip estimation error, and reliable and stable results can be obtained.

[0022] Further, the structural constraint regularization term includes the constraint parallel to the structural direction using the L2 norm and the constraint perpendicular to the structural direction using the L1 norm . The constraint parallel to the structural direction The constraint employs L2 norm The purpose is to make the inverted velocity model more smooth along the structural direction, while perpendicular to the structural direction The constraint employs L1 norm The purpose is to make the inverted velocity model more sparse perpendicular to the structural direction, which ensures the inverted velocity model has the function of blockiness, and is more consistent with the assumption of the actual underground velocity model.

[0023] Further, the wave field is obtained by the attenuation compensation wave field continuation algorithm, but since the L1 norm regularization is not convenient for direct derivation, the Split-Bregman iteration method is introduced to convert the L1 norm constraint problem into an equivalent L2 norm constraint, so as to uniformly solve the derivative of the structural constraint on the velocity model. After conversion, the structural constraint has the following L2 norm form:

[0024]

[0025]

[0026] Wherein, m1, m2, n1 and n2 are iteration variables.

[0027] Further, in the step S3, when the seismic wave propagates underground, it is affected by absorption and attenuation, resulting in weakening of the amplitude of the seismic wave. The decoupled viscous wave equation is used to describe this process. The continuation formula of the decoupled viscous wave equation in the time direction is:

[0028]

[0029] Wherein, u(x, t+NΔt) is the wave field at the spatial position x at the time t+NΔt, k is the wave number, i is the imaginary unit, is the wave field in the Fourier domain at the time (N-1)Δt, is a phase shift function of the compensation wave field; the wave field continuation operator is decomposed by using a low-rank decomposition algorithm, and part of the elements in the partial wave field continuation operator is selected to approximate the propagation operator, and the calculation efficiency of the method is improved by fast Fourier transform.

[0030] Further, by means of the adjoint state method, the gradient of the objective function with respect to the velocity can be equivalent to the correlation of the forward wave field and the backward wave field and the product of the velocity derivative with respect to the forward operator. The velocity field update gradient expression can be written in the following form:

[0031]

[0032] Wherein, the superscript T represents the conjugate transpose operator, J Qis the objective function, Q is the quality factor used to compensate for the energy loss during wavefield propagation, u c (v c is the compensated forward wavefield, A(Q) is the forward operator, is the compensated forward operator, R is the sampling operator, d c is the data residual, m1, m2, n1 and n2 are iteration variables, λ1, λ2 are the corresponding and is the construction of the constraint regularization weight factor. Where the earth is viscoelastic, it will distort the amplitude and phase of the propagating seismic wave, and the attenuation of the seismic wave can be quantified by the quality factor Q, which describes the relationship between the amplitude loss and phase distortion of the seismic wave during propagation and the propagation distance. A lower Q value means that the energy loss or attenuation of each periodic wave is larger. The objective functional is solved by using the optimization iteration method, which combines the advantages of the split Bregman iteration method and the adjoint state method. The L1 norm optimization problem which is not easy to derive directly is converted into an equivalent L2 norm optimization problem to obtain reliable and effective velocity modeling results. Where the gradient of the constraint attenuation compensation velocity modeling is divided into two parts of the weighted sum, the first part is i.e. the conventional gradient, which is the autocorrelation of the forward wavefield and the backward residual wavefield multiplied by the derivative of the forward operator with respect to the velocity; the second part is which is the derivative of the adaptive direction operator with respect to the model parameter. The additional adaptive direction operator provides a more accurate update gradient, accelerates the convergence rate of the full waveform inversion, and improves the accuracy of the inversion result.

[0033] Further, the multi-scale multi-band inversion strategy in the step S4 comprises the following steps:

[0034] Step S41, first use low-frequency data (frequency range between 5-9Hz) to recover the background velocity, and continuously increase the frequency band width of the data used in the inversion process (the frequency range is usually up to 5-20Hz), and apply additional high-frequency data information to the full waveform inversion to obtain a high-resolution inversion result.

[0035] Step S42, the dip angle field information is updated with the increase of the number of iterations, and the dip angle field is calculated by using the imaging result obtained by using the iteration velocity after each iteration, and the new dip angle field is used for constraint in the subsequent iteration.

[0036] Further, the iteration number is adjusted according to the accuracy of the initial model used, and the iteration number can be appropriately reduced when the accuracy of the initial model used is high. According to the inversion simulation result, the iteration number is 5-10 times, which is relatively stable, and the local structure information can effectively depict the detailed information in the velocity model. Compared with the conventional FWI, the structure constraint can provide more accurate update information, and effectively improve the inversion accuracy.

[0037] Compared with the prior art, the beneficial effects of the present application are that on the basis of the conventional full waveform inversion theory, the attenuation compensation is applied to the full waveform inversion to improve the inversion convergence rate by using the fractional-order decoupled viscoacoustic wave equation, the structure-guided constraint is introduced into the full waveform inversion objective function, the geological structure dip information is obtained through plane wave decomposition, and the velocity parameter model is constrained by the seismic imaging structure information, so that the inversion accuracy and efficiency can be effectively improved. BRIEF DESCRIPTION OF DRAWINGS

[0038] Figure 1 The flowchart of the attenuation compensation velocity modeling method of the present application is shown in the figure.

[0039] Figure 2 The deep water actual seismic shot gather data is shown in the figure.

[0040] Figure 3 The initial velocity model is shown in the figure.

[0041] Figure 4 The final velocity model obtained by the attenuation compensation velocity modeling method of the present application is shown in the figure. DETAILED DESCRIPTION

[0042] The present application will be further described below in combination with specific embodiments. The accompanying drawings are only used for illustrative description, and cannot be understood as the limitation of the present patent; in order to better illustrate the present embodiments, some components in the drawings may be omitted, enlarged or reduced, and do not represent the size of the actual product; for those skilled in the art, it is understandable that some well-known structures and their descriptions in the drawings may be omitted.

[0043] The same or similar reference numerals in the drawings of the embodiments of the present application correspond to the same or similar components; in the description of the present application, it is understood that if the orientations or positional relationships indicated by the terms "front", "back", "left", "right" and the like are based on the orientations or positional relationships shown in the drawings, they are only for the convenience of describing the present application and simplifying the description, and do not indicate or imply that the devices or elements referred to must have a particular orientation, be constructed and operated in a particular orientation, therefore the terms describing the positional relationship in the drawings are only used for exemplary illustration, and cannot be understood as a limitation on the present patent. For those of ordinary skill in the art, the specific meanings of the above terms can be understood according to the specific circumstances. In addition, the descriptions such as "first", "second" and the like in the present application are only for the purpose of description, and cannot be understood as indicating or implying the relative importance of the technical features indicated or implying the number of technical features indicated. Therefore, the features limited by "first", "second" can be explicitly or implicitly included at least one of the features.

[0044] Embodiment one:

[0045] Referring to Figure 1 The embodiment provides a deep water shallow layer seismic data structure constraint attenuation compensation velocity modeling method, comprising the following steps:

[0046] Step S1, using a tomographic inversion method to obtain an initial velocity model and performing reverse-time migration imaging to obtain an initial imaging result, and using an inclination prediction technology to obtain a seismic data inclination field;

[0047] Step S2, using the inclination field and multi-frequency band seismic data to establish a structure constraint velocity modeling target functional;

[0048] Step S3, using decoupled viscous wave equation to carry out Q compensation wave field continuation, and calculating the gradient of the functional to update the velocity model;

[0049] Step S4, using a multi-scale multi-frequency band inversion strategy to solve a high-precision velocity inversion optimization problem to obtain a high-precision deep water shallow layer velocity field.

[0050] Referring to Figures 2 to 4 , Figure 2 For 3-shot data in deep water actual seismic shot data, Figure 3 is an initial velocity model, and through the steps S1 to S4, a final velocity model as shown in Figure 4 is obtained by using a different scale multi-frequency band iteration method. Wherein, attenuation compensation (Q compensation) is applied to full waveform inversion to improve the convergence rate of inversion, a structure guiding constraint can be introduced in the full waveform inversion target function, the inclination information of geological structure is obtained by plane wave decomposition, and the velocity parameter model is constrained by seismic imaging structure information, which can effectively improve the precision and efficiency of inversion.

[0051] Embodiment two:

[0052] On the basis of the embodiment one, the step S2 further comprises the following steps:

[0053] Step S21, introducing a constraint term of the inversion model in the velocity modeling inversion objective function, effectively improving the stability and convergence rate of the inverse problem, introducing the dip angle field information and the adaptive total variation constraint into the objective function, and establishing a constrained attenuation compensation velocity modeling objective functional:

[0054]

[0055] Wherein, J Q is the objective function, u(v c , Q) represents the simulation data considering attenuation, f(v r , Q) represents the observation data containing attenuation, v c represents the current velocity model, v r represents the final velocity model, is a constrained regularization term, and are gradient operators parallel and perpendicular to the velocity model construction direction respectively.

[0056] Step S22, constructing the constrained regularization term in the equation, λ1 and λ2 are respectively the construction constraint regularization weight factors corresponding to and The expressions of the gradient operators and can be respectively expanded as:

[0057]

[0058]

[0059] Wherein, ε1 and ε2 are respectively the weighting coefficients corresponding to and , and and are gradient operators along the horizontal direction and the vertical direction respectively, and v(x, z) represents the velocity model; using the dip angle prediction technology, the local dip angle information corresponding to each point of the velocity model is obtained, and the local dip angle information is combined with the horizontal and vertical direction gradient operators, so that the gradient constraint operators parallel and perpendicular to the construction direction are obtained. Compared with the commonly used horizontal and vertical direction gradient operators, the new construction constraint operators are and The boundary information of the inversion model can be protected more effectively, and this process only needs a small amount of calculation relative to full waveform inversion; after the boundary information of the model is detected, the obtained boundary information is used to correct the gradient of the data residual, so that the model gradient is more sensitive to the interface.

[0060] In step S23, the parameters ε1 and ε2 are set as parameters related to the dip angle p, that is, ε1(p) and ε2(p), the larger the dip angle p, the smaller ε1(p) and the larger ε2(p), and vice versa. Such constraints ensure that the inversion constraint process can be automatically adjusted in the case of large dip angle estimation error, and reliable and stable results can be obtained.

[0061] In step S24, the constraint regularization term is constructed including the constraint parallel to the structure direction using the L2 norm , and the constraint perpendicular to the structure direction using the L1 norm . The constraint parallel to the structure direction uses the L2 norm , and the purpose is to make the inverted velocity model more smooth along the structure direction, and the constraint perpendicular to the structure direction uses the L1 norm , and the purpose is to make the inverted velocity model more sparse perpendicular to the structure direction, so as to ensure that the inverted velocity model has a block function and is more consistent with the assumption of the actual underground velocity model.

[0062] In step S25, the wave field is obtained by the attenuation compensation wave field continuation algorithm, but the L1 norm regularization is not convenient for direct derivation, so the Split-Bregman iteration method is introduced to convert the L1 norm constraint problem into an equivalent L2 norm constraint, so as to uniformly solve the derivative of the structure constraint to the velocity model. After conversion, the structure constraint has the following L2 norm form:

[0063]

[0064]

[0065] Wherein, m1, m2, n1 and n2 are iteration variables.

[0066] Embodiment three:

[0067] On the basis of embodiment two, in step S3, when the seismic wave propagates underground, the amplitude of the seismic wave is weakened due to the effect of absorption and attenuation, and a decoupled viscous wave equation is used to describe this process. The continuation formula in the time direction of the decoupled viscous wave equation is:

[0068]

[0069] where u(x, t+NΔt) is the wavefield at spatial position x and time t+NΔt, k is the wave number, is the wavefield in Fourier domain at time (N-1)Δt, is the phase shift function of the compensated wavefield; the wavefield extrapolation operator is decomposed by low-rank decomposition algorithm, and part of the elements in the wavefield extrapolation operator is selected to approximate the propagation operator, and the calculation efficiency of the method is improved by fast Fourier transform.

[0070] Further, by the adjoint state method, the gradient of the objective function with respect to the velocity can be equivalent to the correlation of the forward wavefield and the backward residual wavefield multiplied by the derivative of the forward operator with respect to the velocity, and the velocity field update gradient expression can be written as follows:

[0071]

[0072] where the superscript T represents the conjugate transpose operator, J Q is the objective function, Q is the quality factor used to compensate for the energy loss in the wavefield propagation process, u c (v c , Q) is the compensated forward wavefield, A(Q) is the forward operator, is the compensated forward operator, R is the sampling operator, d c is the data residual, m1, m2, n1 and n2 are iteration variables, λ1 and λ2 are the regularization weight factors constructed according to and respectively. Where the earth is viscoelastic, which distorts the amplitude and phase of the propagating seismic wave, and the attenuation of the seismic wave can be quantified by the quality factor Q, which describes the relationship between the amplitude loss and the phase distortion of the seismic wave during propagation and the propagation distance. A lower Q value means greater energy loss or greater attenuation of each periodic wave. The objective functional is solved by using the optimization iteration method, which combines the advantages of the split Bregman iteration method and the adjoint state method, and converts the L1 norm optimization problem which is not easy to derive directly into an equivalent L2 norm optimization problem to obtain reliable and effective velocity modeling results.

[0073] where the gradient of the constraint attenuation compensated velocity modeling is divided into the weighted sum of two parts, the first part is i.e. the conventional gradient, which is the autocorrelation of the forward wavefield and the backward residual wavefield multiplied by the derivative of the forward operator with respect to the velocity; the second part is which is the derivative of the adaptive direction operator with respect to the model parameter. The additional adaptive direction operator provides a more accurate update gradient, accelerates the convergence rate of the full waveform inversion, and improves the accuracy of the inversion result.

[0074] In the specific contents of the foregoing specific embodiments, each technical feature can be combined arbitrarily without contradiction. In order to make the description simple, all possible combinations of the foregoing technical features are not described, however, as long as the combinations of the technical features do not contradict, they should be considered as the scope of the present disclosure.

[0075] Obviously, the above embodiments of the present application are merely exemplary and are not intended to limit the implementation modes of the present application. Based on the above description, other different forms of changes or variations can be made by those skilled in the art. Here, it is not necessary and also impossible to exhaust all the implementation modes. Any modification, equivalent replacement and improvement, etc. made within the spirit and principle of the present application should be included in the protection scope of the claims of the present application.

Claims

1. A method for constructing a constrained velocity model for deep water shallow seismic data with attenuation compensation, characterized in that, The method comprises the following steps: Step S1, obtaining an initial velocity model by using a tomographic inversion method and performing reverse-time migration imaging to obtain an initial imaging result, and using an angle prediction technique to obtain a seismic data dip angle field; Step S2, establishing a structure-constrained velocity modeling target functional based on the dip angle field and multi-band seismic data, wherein in the step S2, a constraint term for an inversion model is introduced into a velocity modeling inversion optimization target function, thereby effectively improving the stability and convergence rate of the inverse problem, the dip angle field information and an adaptive total variation constraint are introduced into the target function, and a structure-constrained attenuation compensation velocity modeling target functional is established: wherein, is the objective function, represents simulated data taking into account attenuation, represents observed data containing attenuation, is a construction constraint regularization term; wherein, represents a current velocity model, represents a final velocity model, and are gradient operators parallel and perpendicular to the velocity model construction direction, respectively, , are construction constraint regularization weight factors corresponding to and , respectively. Step S3, performing wave field continuation by using a decoupled viscoacoustic wave equation for Q compensation, and calculating the gradient of the functional to update the velocity model; Step S4, solving a high-precision velocity inversion optimization problem by using a multi-scale multi-band inversion strategy to obtain a high-precision deep water shallow layer velocity field.

2. The deep water shallow seismic data structure-constrained attenuation-compensated velocity modeling method according to claim 1, characterized in that, In the step S1, first, initial imaging is performed by using an initial velocity model to obtain an initial imaging result, and then a mathematical relationship between travel time and the dip angle field is used to obtain a travel time difference: （1） to obtain the dip field information of the initial model of the earthquake, wherein represents a zero-offset travel time curve, represents a travel time, represents a first trace travel time, represents a offset distance, represents a zero-offset, represents a dip field.

3. The method according to claim 1, wherein, is a construction constraint regularization term, where, , are construction constraint regularization weight factors corresponding to and respectively; the expressions of the gradient operators and can be extended to: wherein, and are weighting coefficients corresponding to and respectively, and are gradient operators along the horizontal and vertical directions respectively, represents a velocity model; using the dip prediction technique, local dip information corresponding to each point of the velocity model is obtained, and the local dip information is combined with the horizontal and vertical direction gradient operators, so that gradient constraint operators parallel and perpendicular to the structural direction are obtained. In addition, parameters and Set with tilt angle The relevant parameters, namely and ,inclination The larger, The smaller, The larger the angle, the lower the inclination. The smaller, The larger, The smaller.

4. The method according to claim 3, wherein, The construction constraint regularizing term includes adoption of L2 norm parallel to the construction direction and adoption of L1 norm perpendicular to the construction direction .

5. The method according to claim 4, wherein, The split Bregman iteration method is introduced to convert the L1 norm constraint problem which is not convenient for direct derivation into an equivalent L2 norm constraint, so as to uniformly solve the derivative of the structure constraint with respect to the velocity model, and the transformed structure constraint is in the following L2 norm form: wherein , , and are iteration variables.

6. The method according to claim 1, wherein, In the step S3, the time direction continuation formula of the decoupled viscoacoustic wave equation is: wherein, is the spatial position at the time instant, is the wave number, i is the imaginary unit, is the wave field in the Fourier domain at the time instant, is the phase shift function compensating the wave field; the wave field extrapolation operator is decomposed using a low-rank decomposition algorithm, and some elements in some wave field extrapolation operators are selected to approximate the propagation operator, and the calculation efficiency of the method is improved by fast Fourier transform.

7. The method according to claim 6, wherein, By using the adjoint state method, the gradient of the target function with respect to the velocity can be equivalent to the correlation of the forward wave field and the backward wave field and the product of the forward operator derivative with respect to the velocity, and the velocity field update gradient expression can be written in the following form: where the superscript T denotes the conjugate transpose operator, is the objective function, is the quality factor used to compensate for energy loss during wavefield propagation, is the forward wavefield with compensation, is the forward operator, is the forward operator with compensation, is the sampling operator, is the data residual, , , and are the iteration variables, , are the constructed constraint regularization weight factors corresponding to and , respectively.

8. The method according to claim 1, wherein, The multi-scale multi-band inversion strategy in the step S4 comprises the following steps: Step S41, using low-frequency data in a frequency range of 5-9 Hz to restore the background velocity, and then adding high-frequency data in a frequency range of 5-20 Hz to the full waveform inversion to obtain a high-resolution velocity inversion result; Step S42, updating the dip angle field information as the number of iterations increases, and obtaining an imaging result by using the iterated velocity every several iterations to calculate the dip angle field, and using the new dip angle field for constraint in subsequent iterations.

9. The method according to claim 8, wherein, The number of iterations is 5-10.

Images