Pre-stack earthquake wide-angle AVO inversion method based on second-order approximation

Abstract

The invention relates to the technical field of pre-stack seismic inversion, and discloses a pre-stack seismic wide-angle AVO inversion method based on a second-order approximation, and the method comprises the steps: constructing an inversion objective function through a seismic plane wave second-order approximation on a solid / solid medium horizontal interface; and solving the inversion objective function by adopting an alternating direction multiplier method to realize high-resolution inversion of the longitudinal wave velocity, the transverse wave velocity and the density. According to the method, an inversion objective function is constructed by using a second-order approximation expression, the complex inversion objective function is converted into a plurality of single-parameter linear sub-problems easy to solve by using an alternating direction multiplier algorithm, and meanwhile, a Hadamard product operator is introduced to decompose and solve a high-order sub-problem. Model test and practical application show that the longitudinal wave velocity, the transverse wave velocity and the density resolution ratio inverted by the method are higher, the thin layer identification capability is stronger, the pseudo layer phenomenon can be greatly weakened, and a certain support can be provided for improving the prediction precision of the complex oil and gas reservoir.

Description

Technical Field

[0001] The present invention relates to the technical field of pre-stack seismic inversion, and in particular to a pre-stack seismic wide-angle AVO inversion method based on a second-order approximation. Background Art

[0002] Prestack seismic inversion is a key technique for obtaining subsurface rock elastic parameters (such as P-wave velocity, S-wave velocity, and density) from observed prestack seismic data. These elastic parameters are crucial for reservoir prediction. The variation of seismic wave reflection amplitude with incident angle (AVO) or offset (AVA) is the foundation of prestack seismic inversion. This variation is primarily influenced by contrasting rock properties. The Zoeppritz equation (Zoeppritz, 1919) establishes a relationship between the reflection amplitude of a P-wave incident on an adjacent isotropic, homogeneous, single-phase medium and its variation with offset. However, due to its complexity, this equation is difficult to directly apply to actual prestack seismic inversion and interpretation. Consequently, numerous researchers have simplified the Zoeppritz equation in recent decades. Bortfeld's (1961) approximation greatly simplified the calculation process, making it more efficient in practical applications; Aki and Richards (1980) derived an approximate equation for the P-wave reflection coefficient based on P-wave and S-wave velocities and density under the assumptions of isotropy and weak contrast in medium parameters on both sides of the interface; Shuey (1985) established an approximate equation for the reflection coefficient based on Poisson's ratio, S-wave velocity and density, which can effectively highlight the characteristics of oil and gas; Fatti (1994) established an approximate equation for the reflection coefficient of P-wave impedance, S-wave impedance and density, which realizes the efficient inversion of P-wave and S-wave impedance parameters from seismic data; in order to better identify fluids, Gray and Verttas (2002) proposed an approximate equation for the reflection coefficient to directly obtain the Lame constant and shear modulus, realizing the direct extraction of Lame parameters from pre-stack seismic data.

[0003] Improving inversion accuracy has long been a concern for many scholars both domestically and internationally. However, due to limitations in seismic data quality and inversion theory, the solution to the inverse problem is often ill-posed, meaning that more than one reflection coefficient sequence satisfies the inversion equation. With the continuous advancement of exploration and development, oil and gas reservoirs with simple geological structures and easy access are becoming increasingly rare, while interpretation theories and methods for reservoirs with complex geological structures and deeper burials are limited. To overcome the increasingly complex challenges of oil and gas exploration and development, some studies have considered implementing three-parameter (P-wave velocity, S-wave velocity, and density) generalized linear or nonlinear inversion based on the Zoeppritz equation (Zhi et al., 2018; Gholami et al., 2018; Zhang et al., 2013; Liu et al., 2012) and Bayesian inversion (Zhou et al., 2017). However, the Zoeppritz equation, in addition to P-wave reflection, also couples scattered waves such as S-wave reflection, P-wave transmission, and S-wave transmission. Therefore, some researchers have used the Zoeppritz equation to achieve a joint PP-PS wave inversion (Lu et al., 2015). Others have considered using the exact analytical expression for the P-wave reflection coefficient derived by Aki-Richards (1980) and combining it with nonlinear or generalized linear inversion algorithms to solve it (Zhi and Gu, 2018; Zong et al., 2013). Others have derived a third-order approximation to this exact analytical expression, enabling three-parameter nonlinear inversion (Cheng et al., 2018). Yang et al. (2023) derived a quadratic PP reflection approximation for the Graebner equation and solved it using the linear alternating direction multiplier method, improving the accuracy of the density inversion. Furthermore, some researchers have proposed a method for inverting three parameters based on the wave equation to achieve higher-precision inversion (Sun Chengyu et al., 2019).

[0004] Prestack AVO inversion is a key technique for obtaining subsurface rock elastic parameters from observed seismic data. These elastic parameters are crucial for reservoir prediction. Prestack inversion methods, based on the Aki and Richards approximation and its variations, are only applicable to low-angle incidence scenarios. The lack of high-angle gather information for inversion results in low shear wave velocity and density inversion resolution. Summary of the Invention

[0005] The purpose of the present invention is to provide a pre-stack seismic wide-angle AVO inversion method based on a second-order approximation to address the problem of low resolution of shear wave velocity and density inversion due to the lack of large-angle gather information in the inversion. To solve this problem, the present invention proposes a three-parameter wide-angle AVO inversion using a second-order approximation. Analysis of a sandstone medium model shows that the second-order approximation has higher reflection coefficient accuracy than the first-order approximation in a large incident angle range (>30°). The inversion objective function is constructed using the second-order approximation, and the alternating direction multiplier algorithm (ADMM) is used to transform the complex inversion objective function into multiple easily solvable single-parameter linear subproblems. The Hadamard product operator is introduced to decompose and solve higher-order subproblems. Model testing and practical applications show that the proposed method has higher resolution of P-wave velocity, S-wave velocity, and density inverted, has stronger thin-layer recognition capabilities, and can greatly reduce the pseudo-layer phenomenon, which can provide certain support for improving the prediction accuracy of complex oil and gas reservoirs.

[0006] In order to achieve the above purpose, the technical solutions adopted are as follows:

[0007] The present invention provides a pre-stack seismic wide-angle AVO inversion method based on a second-order approximation, the method comprising:

[0008] The inversion objective function is constructed using the second-order approximation of seismic plane waves on the solid / solid medium horizontal interface;

[0009] The alternating direction multiplier method is used to solve the inversion objective function to achieve high-resolution inversion of P-wave velocity, S-wave velocity and density. Furthermore, the second-order approximation of the seismic plane wave on the solid / solid medium horizontal interface is expressed as:

[0010] R(θ)=R1(θ)+R2(θ) (1)

[0011]

[0012] where R(θ) is the second-order approximation of the seismic plane wave at the solid / solid medium horizontal interface, R1(θ) is the first-order term of the reflection coefficient, and R2(θ) is the second-order term of the reflection coefficient; θ is the average incident angle of the upper and lower layers; α, β, and ρ are the average P-wave velocity, S-wave velocity, and density of the upper and lower layers, respectively; Δα, Δβ, and Δρ represent the differences in P-wave velocity, S-wave velocity, and density between the lower and upper layers, respectively; c1(θ) is the controlling coefficient of the first-order P-wave velocity reflection coefficient term, c2(θ) is the controlling coefficient of the first-order S-wave velocity reflection coefficient term, c3(θ) is the controlling coefficient of the first-order density reflection coefficient term, x1(θ) is the controlling coefficient of the second-order S-wave velocity reflection coefficient term, x2(θ) is the controlling coefficient of the combined second-order S-wave velocity and density reflection coefficient term, and x3(θ) is the controlling coefficient of the second-order density reflection coefficient term. Y(θ) and H(θ) are the influencing factors of the second-order reflection coefficient.

[0013] Furthermore, using the second-order approximation of seismic plane waves on the solid / solid medium horizontal interface, the inversion objective function is constructed as follows:

[0014] Describing the angle gather as a form of convolution of an angle wavelet and an angle reflection coefficient sequence; wherein the form of convolution of the angle wavelet and the angle reflection coefficient sequence includes seismic responses generated by first-order and second-order terms of the reflection coefficient;

[0015] Based on the seismic response generated by the first-order term of the reflection coefficient, a set of equations is established by stacking data from several partial angle gathers.

[0016] Based on the seismic response generated by the second-order term of the reflection coefficient, a relationship is established by stacking data from several partial angle gathers.

[0017] According to the equations and relations, a three-parameter inversion equation based on a second-order approximation of the reflection coefficient is obtained;

[0018] Under the weak elasticity assumption, a low-frequency model constraint equation is obtained according to the relationship between the three-parameter reflection terms of the three-parameter inversion equation and the logarithms of the corresponding parameters;

[0019] An inversion objective function is constructed according to the three-parameter inversion equation and the low-frequency model constraint equation.

[0020] Furthermore, the convolution of the angle wavelet and the angle reflection coefficient sequence is expressed as:

[0021] d(θ)=d1(θ)+d2(θ)=w(θ)*R1(θ)+w(θ)*R2(θ) (10)

[0022] d1(θ)=c1(θ)w(θ)*r P +c2(θ)w(θ)*r S+c3(θ)w(θ)*r ρ (11)

[0023] d2(θ)=x1(θ)w(θ)*r S ·r S +x2(θ)w(θ)*r S ·r ρ +x3(θ)w(θ)*r ρ ·r ρ (12)

[0024] Where d1(θ) and d2(θ) represent the seismic responses generated by the first-order and second-order terms of the reflection coefficient, respectively; r P 、r S and r ρ Represent the reflection coefficient sequences of P-wave velocity, S-wave velocity and density respectively; * represents the convolution operator; . represents the Hadamard product operator; w(θ) represents the angle wavelet; d(θ) represents the angle gather.

[0025] Furthermore, the system of equations is expressed as:

[0026]

[0027] d1=G1r(14)

[0028] Where, d1(θ i ) represents the seismic response generated by the first-order term of the i-th (i=1, 2, ..., h) angular reflection coefficient, c1(θ i ) represents the control coefficient of the first-order longitudinal wave velocity reflection coefficient term at the i-th (i=1, 2, ..., h) angle, c2(θ i ) represents the control coefficient of the first-order shear wave velocity reflection coefficient term at the i-th (i=1, 2, ..., h) angle, c3(θ i ) represents the control coefficient of the first-order density reflection coefficient term at the i-th (i=1, 2, ..., h) angle, W(θ i ) represents the wavelet matrix of the i-th (i=1, 2, ..., h) angle, d1 represents the seismic response generated by the first-order term of the reflection coefficient, G1 represents the wavelet kernel matrix of the first-order term of the reflection coefficient, and r represents the quantity to be solved;

[0029] The relationship is expressed as:

[0030]

[0031] Where, W(θ i ) represents the wavelet matrix of the i-th (i=1, 2, ..., h) angle, I is the same as r S The identity matrix of the same dimension, x1(θ i) represents the control coefficient of the second-order shear wave velocity reflection coefficient term at the i-th (i=1, 2, ..., h) angle, x2(θ i ) represents the control coefficient of the joint term of the second-order shear wave velocity reflection coefficient and density reflection coefficient at the i-th (i=1, 2, ..., h) angle, x3(θ i ) represents the control coefficient of the i-th (i=1, 2, ..., h) angle second-order density reflection coefficient term, K1 and K2 represent elementary transformation matrices, d2 represents the seismic response generated by the second-order reflection coefficient term, and G2 represents the wavelet kernel matrix of the second-order reflection coefficient term;

[0032] Substituting equation (16) into equation (15), we can obtain the simplified form of the relationship, which is expressed as:

[0033] d2=G2(K1r)·(K2r) (17)

[0034] Combining equations (14) and (17), we obtain the three-parameter inversion equation based on the second-order approximation of the reflection coefficient, which is expressed as:

[0035] d=d1+d2=G1r+G2(K1r)·(K2r) (18)

[0036] Where d represents the earthquake record.

[0037] Furthermore, under the weak elasticity assumption, a low-frequency model constraint equation is obtained based on the relationship between the three-parameter reflection term and the logarithm of the corresponding parameter in the three-parameter inversion equation, and an inversion objective function is constructed based on the three-parameter inversion equation and the low-frequency model constraint equation, including:

[0038] Under the weak elasticity assumption, the relationship between the three-parameter reflection term and the logarithm of the corresponding parameter is expressed as:

[0039] ε P =C′r P , ε S =C′r S , ε ρ =C′r ρ (19)

[0040]

[0041] Where α0, β0, and ρ0 represent the initial longitudinal wave velocity, shear wave velocity, and density, respectively. i , β i , ρ i are the P-wave velocity, S-wave velocity and density of the i-th layer respectively, C′ is the integration matrix;

[0042] Combining equations (19) to (22), we get the low-frequency model constraint equation:

[0043]

[0044] ε=Cr( 24 )

[0045] Where η P ,η S and η ρ Respectively represent the model constraint parameters of P-wave velocity, S-wave velocity and density, which are used to adjust the model constraint proportions of P-wave velocity, S-wave velocity and density; ε represents the low-frequency model;

[0046] Combining equations (18) and (24), we construct the following inversion objective function f(r):

[0047]

[0048] Where η>0 is the low-frequency model constraint parameter, which determines the low-frequency model term and fidelity Relative contribution to the inversion results.

[0049] Furthermore, the alternating direction multiplier method is used to solve the inversion objective function to achieve high-resolution inversion of P-wave velocity, S-wave velocity, and density, including:

[0050] By introducing the first variable and the second variable into the inversion objective function, the process of optimizing K1r and K2r is transformed into optimizing the first variable and the second variable, and the inversion objective function is transformed into a multi-constrained optimization problem;

[0051] The Lagrangian parameters λ1 and λ2 are introduced by using the Lagrangian multiplier method, and the constraints are added to the inversion objective function to obtain the unconstrained augmented Lagrangian inversion objective function.

[0052] Define the first dual variable and the second dual variable to obtain the unconstrained inversion objective function;

[0053] The variable to be solved, the first variable, the second variable, the first dual variable and the second dual variable are solved iteratively in sequence until the convergence condition is met, and the longitudinal wave velocity, shear wave velocity and density are output.

[0054] Furthermore, the multi-constraint optimization problem is expressed as:

[0055]

[0056] Where, f η (r) represents a multi-constrained optimization problem, X1 and X2 represent the first and second variables, respectively;

[0057] The inversion objective function of the unconstrained augmented Lagrangian form is expressed as:

[0058]

[0059] Where, L η，κ (r, X1, X2, λ1, λ2) represents the inversion objective function of the unconstrained augmented Lagrangian form;

[0060] Define the first dual variable and the second dual variable, expressed as The obtained unconstrained inversion objective function is expressed as:

[0061]

[0062] Where, L η，κ (r, X1, X2, Z1, Z2) represents the unconstrained inversion objective function; κ is a penalty parameter used to control the convergence speed.

[0063] Furthermore, in the process of iteratively solving the quantity to be solved, the first variable, the second variable, the first dual variable, and the second dual variable in the i-th loop:

[0064] The quantity to be solved r is optimized by the following method:

[0065] Taking r as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0066]

[0067] Where, represents the first variable solved in the previous iteration, represents the second variable solved in the previous iteration, represents the first dual variable solved in the previous iteration, represents the second dual variable solved in the previous iteration;

[0068] The process of solving the single variable optimization problem in equation (29) using the least squares method is as follows:

[0069]

[0070] Where T represents matrix transpose; r i+1 Represents the variable to be solved after the i-th iteration update;

[0071] The first variable X1 is solved by optimizing as follows:

[0072] Taking X1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0073]

[0074] The process of solving the single variable optimization problem in equation (31) using the least squares method is as follows:

[0075]

[0076] Where I is the identity matrix of the same dimension as X1;

[0077] The second variable X2 is solved by optimizing as follows:

[0078] Taking X2 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0079]

[0080] The process of solving the single variable optimization problem in equation (33) using the least squares method is as follows:

[0081]

[0082] The first dual variable Z1 is solved by optimizing as follows:

[0083] Taking Z1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0084]

[0085] The process of solving the single variable optimization problem in equation (35) using the least squares method is as follows:

[0086] Z1=Z1+κ(X1-K1r i+1 ) (36)

[0087] The first dual variable Z2 is solved by optimizing as follows:

[0088] Taking Z1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0089]

[0090] The process of solving the single variable optimization problem in equation (37) using the least squares method is as follows:

[0091] Z2=Z2+κ(X2-K2r i+1 ) (38).

[0092] Furthermore, the convergence condition is Here, tol represents the convergence error.

[0093] The beneficial effects of the present invention are:

[0094] Based on the derived second-order approximate analytical expression for the P-wave reflection coefficient, this paper establishes a second-order constrained AVO inversion equation. The Hadamard product operator is used to simplify the objective function based on the quadratic reflection coefficient. Combined with low-frequency model constraints, a convex optimization inversion objective function is constructed. Taking into account the complexity of the objective function, the alternating direction method of multipliers (ADMM) is employed to iteratively solve the objective function. Finally, model testing and application to real data yield high-resolution P-wave velocity, S-wave velocity, and density parameters, demonstrating the feasibility and effectiveness of this method. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] Figure 1 A flowchart of a pre-stack seismic wide-angle AVO inversion method based on a second-order approximation is provided in an embodiment of the present invention.

[0096] Figure 2 A schematic diagram of the variation of the longitudinal wave reflection coefficient with the incident angle provided in an embodiment of the present invention. In the figure, the three line segments represent the longitudinal wave reflection coefficient calculated based on the exact Zoeppritz equation, the Aki and Richards approximation, and the second-order approximation, respectively.

[0097] Figure 3 A flowchart for constructing the inversion objective function provided in an embodiment of the present invention.

[0098] Figure 4 A flowchart of solving an inversion objective function using an alternating direction multiplier method is provided in an embodiment of the present invention.

[0099] Figure 5 Schematic diagram of the original theoretical model provided by an embodiment of the present invention, wherein (a) longitudinal wave velocity (b) shear wave velocity (c) density.

[0100] Figure 6 Schematic diagram of a low-frequency model provided by an embodiment of the present invention, wherein (a) longitudinal wave velocity (b) shear wave velocity (c) density.

[0101] Figure 7 Schematic diagram of the angle gathers required for conventional AVO inversion provided in an embodiment of the present invention, where (a), (b), and (c) are synthetic seismic records with incident angles of 5°, 15°, and 25°, respectively.

[0102] Figure 8 Schematic diagram of the angle gathers required for wide-angle AVO inversion provided in an embodiment of the present invention, wherein (a), (b), and (c) are synthetic seismic records with incident angles of 8°, 24°, and 40°, respectively.

[0103] Figure 9 A conventional AVO inversion result diagram provided in an embodiment of the present invention, wherein (a) P-wave velocity (b) S-wave velocity (c) density.

[0104] Figure 10 AVO inversion result diagram of the new method provided by an embodiment of the present invention, including (a) longitudinal wave velocity, (b) shear wave velocity, and (c) density.

[0105] Figure 11 Schematic diagram of the 40th three-parameter generated angle gather provided in an embodiment of the present invention, where (a) is a noise-free single-channel angle gather, (b) is a single-channel angle gather with 5% noise, and (c) is a single-channel angle gather with 10% noise.

[0106] Figure 12 A comparison of the pre-stack AVO inversion results of the traditional method and the new method provided in an embodiment of the present invention, wherein (a) the noise-free single-channel inversion result, (b) the single-channel inversion result with 5% noise, and (c) the single-channel inversion result with 10% noise.

[0107] Figure 13 Schematic diagram of partial angle gather stacking data required for 10 traditional pre-stack AVO inversion provided in an embodiment of the present invention, where (a), (b), and (c) are partial angle gather stacking data with incident angles of 5°, 15°, and 25°, respectively.

[0108] Figure 14 Schematic diagram of partial angle gather stacking data required for pre-stack AVO inversion of the new method provided in an embodiment of the present invention, wherein (a), (b), and (c) are partial angle gather stacking data with incident angles of 8°, 24°, and 40°, respectively.

[0109] Figure 15 Schematic diagram of the traditional pre-stack AVO inversion method provided in an embodiment of the present invention, where (a), (b), and (c) are the inversion results of the P-wave velocity, S-wave velocity, and density, respectively.

[0110] Figure 16 Schematic diagram of pre-stack AVO inversion of the new method provided in an embodiment of the present invention, where (a), (b), and (c) are the inversion results of P-wave velocity, S-wave velocity, and density, respectively. DETAILED DESCRIPTION

[0111] The following describes the embodiments of the present invention through specific examples. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments. The details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the following embodiments and features in the embodiments can be combined with each other unless they conflict.

[0112] The specific implementation of the present invention is further described in detail below with reference to the accompanying drawings and examples.

[0113] The present invention provides a pre-stack seismic wide-angle AVO inversion method based on a second-order approximation. Figure 1 As shown, the method includes the following steps:

[0114] S100, constructing the inversion objective function using the second-order approximation of seismic plane waves on the solid / solid medium horizontal interface;

[0115] S200, using the alternating direction multiplier method to solve the inversion objective function and achieve high-resolution inversion of P-wave velocity, S-wave velocity and density.

[0116] In some embodiments, in step S100, the scattering equation of seismic plane waves at the solid / solid horizontal medium interface can be described by Zoeppritz (Zoeppritz, 1916). Furthermore, Aki and Richards (1980) derived its exact analytical expression and first-order linear approximation, laying the theoretical foundation for pre-stack seismic AVO inversion. Based on previous research results, Zhou et al. (2022) derived a second-order approximation for the longitudinal wave reflection coefficient at the horizontal medium interface, R(θ).

[0117] R(θ)=R1(θ)+R2(θ) (1)

[0118]

[0119]

[0120] Where: R1(θ) is the first-order term of the reflection coefficient, i.e., the Aki and Richards approximation, and R2(θ) is the second-order term of the reflection coefficient; θ is the average incident angle of the upper and lower media; α, β, and ρ are the average P-wave velocity, S-wave velocity, and density of the upper and lower layers, respectively. The subscripts 1 and 2 represent the parameters of the upper and lower media, respectively. Δα, Δβ, and Δρ represent the differences in P-wave velocity, S-wave velocity, and density between the lower and upper layers, respectively. c1(θ) is the controlling coefficient of the first-order P-wave velocity reflection coefficient term, c2(θ) is the controlling coefficient of the first-order S-wave velocity reflection coefficient term, c3(θ) is the controlling coefficient of the first-order density reflection coefficient term, x1(θ) is the controlling coefficient of the second-order S-wave velocity reflection coefficient term, x2(θ) is the controlling coefficient of the combined term of the second-order S-wave velocity and density reflection coefficients, and x3(θ) is the controlling coefficient of the second-order density reflection coefficient term. Y(θ) and H(θ) are the influencing factors of the second-order reflection coefficient.

[0121] The analysis of the variation of reflection coefficient with incident angle is the basis of pre-stack AVO inversion. In order to compare and analyze the variation of Aki and Richards approximation and second-order approximation with incident angle, a single reflection interface sandstone medium model (Table 1) was established in this example. The analysis results are shown in Table 1. Figure 2 shown.

[0122] Table 1 Sandstone model parameters

[0123]

[0124] analyze Figure 2 It can be seen that the Aki and Richards approximation is effective within the incident angle range of less than 30°. As the incident angle increases, the accuracy of the Aki and Richards approximation decreases. However, the second-order approximation is effective within the incident angle range of less than 45°, and it still has higher accuracy when the incident angle is greater than 26°. This shows that the Aki and Richards approximation is only applicable to prestack AVO inversion within a small angle range, while the second-order approximation is applicable to prestack AVO inversion within a larger angle range. Moreover, wide-angle seismic data can contain more subsurface medium information, thereby improving the resolution of prestack inversion.

[0125] In some embodiments, as Figure 3 As shown in Figure 2, using the second-order approximation of seismic plane waves on the solid / solid medium horizontal interface, the inversion objective function is constructed through the following steps:

[0126] S101. Describe the angle gather as a form of convolution of an angle wavelet and an angle reflection coefficient sequence; wherein the form of convolution of the angle wavelet and the angle reflection coefficient sequence includes seismic responses generated by first-order and second-order terms of the reflection coefficient.

[0127] Exemplarily, the convolution of the angle wavelet and the angle reflection coefficient sequence is expressed as:

[0128] d(θ)=d1(θ)+d2(θ)=w(θ)*R1(θ)+w(θ)*R2(θ) (10)

[0129] d1(θ)=c1(θ)w(θ)*r P +c2(θ)w(θ)*r S +c3(θ)w(θ)*r ρ (11)

[0130] d2(θ)=x1(θ)w(θ)*r S ·r S +x2(θ)w(θ)*r S ·r ρ +x3(θ)w(θ)*r ρ·r ρ (12)

[0131] Where d1(θ) and d2(θ) represent the seismic responses generated by the first-order and second-order terms of the reflection coefficient, respectively; r P 、r S and r ρ Represent the reflection coefficient sequences of P-wave velocity, S-wave velocity and density respectively; * represents the convolution operator; . represents the Hadamard product operator; w(θ) represents the angle wavelet; d(θ) represents the angle gather.

[0132] S102. Based on the seismic response generated by the first-order term of the reflection coefficient, a set of equations is established using the superimposed data of several partial angle gathers.

[0133] Exemplarily, the system of equations is expressed as:

[0134]

[0135] d1=G1r (14)

[0136] Where, d1(θ i ) represents the seismic response generated by the first-order term of the i-th (i=1, 2, ..., h) angular reflection coefficient, c1(θ i ) represents the control coefficient of the first-order longitudinal wave velocity reflection coefficient term at the i-th (i=1, 2, ..., h) angle, c2(θ i ) represents the control coefficient of the first-order shear wave velocity reflection coefficient term at the i-th (i=1, 2, ..., h) angle, c3(θ i ) represents the control coefficient of the first-order density reflection coefficient term at the i-th (i=1, 2, ..., h) angle, W(θ i ) represents the wavelet matrix of the i-th (i=1, 2, ..., h) angle, d1 represents the seismic response generated by the first-order term of the reflection coefficient, G1 represents the wavelet kernel matrix of the first-order term of the reflection coefficient, and r represents the quantity to be solved.

[0137] S103. Based on the seismic response generated by the second-order term of the reflection coefficient, a relationship is established by using the superimposed data of several partial angle gathers.

[0138] Exemplarily, the relationship is expressed as:

[0139]

[0140]

[0141] Where, W(θ i ) represents the wavelet matrix of the i-th (i=1, 2, ..., h) angle, I is the unit matrix of the same dimension as rS, x1(θ i) represents the control coefficient of the second-order shear wave velocity reflection coefficient term at the i-th (i=1, 2, ..., h) angle, x2(θ i ) represents the control coefficient of the joint term of the second-order shear wave velocity reflection coefficient and density reflection coefficient at the i-th (i=1, 2, ..., h) angle, x3(θ i ) represents the control coefficient of the i-th (i=1, 2, ..., h) angle second-order density reflection coefficient term, K1 and K2 represent elementary transformation matrices, d2 represents the seismic response generated by the second-order reflection coefficient term, and G2 represents the wavelet kernel matrix of the second-order reflection coefficient term;

[0142] Substituting equation (16) into equation (15), we can obtain the simplified form of the relationship, which is expressed as:

[0143] d2=G2(K1r)·(K2r) (17)

[0144] S104. According to the equation group and the relationship, a three-parameter inversion equation based on the second-order approximation of the reflection coefficient is obtained.

[0145] For example, by combining equations (14) and (17), a three-parameter inversion equation based on the second-order approximation of the reflection coefficient is obtained, which is expressed as:

[0146] d=d1+d2=G1r+G2(K1r)·(K2r) (18)

[0147] Where d represents the earthquake record.

[0148] S105. Under the weak elasticity assumption, a low-frequency model constraint equation is obtained based on the relationship between the three-parameter reflection terms of the three-parameter inversion equation and the logarithms of the corresponding parameters.

[0149] For example, the actual seismic data usually received lacks low-frequency components, which increases the uncertainty of the pre-stack inversion results. Therefore, it is necessary to combine logging data, layer data and prior geological knowledge to compensate for the missing low-frequency components. Under the weak elasticity assumption, there is the following relationship between the three-parameter reflection term and the logarithm of the corresponding parameter:

[0150] ε P =C′r P , ε S =C′r S , ε ρ =C′r ρ (19)

[0151]

[0152] Where α0, β0, and ρ0 represent the initial longitudinal wave velocity, shear wave velocity, and density, respectively. i , β i , ρi are the longitudinal wave velocity, shear wave velocity and density of the i-th layer respectively, and C′ is the integral matrix.

[0153] By combining equations (19) to (22), we can obtain the low-frequency model constraint equation as follows:

[0154]

[0155] ε=Cr (24)

[0156] Where η P , η S , η ρ They represent the model constraint parameters of P-wave velocity, S-wave velocity, and density, respectively, and are used to adjust the model constraint proportions of P-wave velocity, S-wave velocity, and density; ε represents the low-frequency model.

[0157] S106. Construct an inversion objective function based on the three-parameter inversion equation and the low-frequency model constraint equation.

[0158] For example, by combining equations (18) and (24), the following inversion objective function f(r) is constructed:

[0159]

[0160] Where η>0 is the low-frequency model constraint parameter, which determines the low-frequency model term and fidelity The relative contribution to the inversion result. The larger the value, the greater the contribution of the low-frequency model term to the inversion, and the more the inversion result tends to the low-frequency model.

[0161] In some embodiments, the inversion objective function formula (25) is an L2 norm optimization problem in terms of expression. However, this expression contains third- and fourth-order terms of the quantity to be solved r, which is difficult to solve using the conjugate gradient method. This embodiment adopts the ADMM algorithm, which introduces dual variables to decompose the original complex function into multiple simple sub-functions for iterative solution. Figure 4 As shown in Figure 2, the specific solution process includes:

[0162] S201 , by introducing the first variable and the second variable into the inversion objective function, the process of optimizing K1r and K2r is transformed into the process of optimizing the first variable and the second variable, and the inversion objective function is transformed into a multi-constrained optimization problem.

[0163] For example, by introducing the first variable X1 and the second variable X2 into the inversion objective function (Equation 25), the process of optimizing K1r and K2r is transformed into the first variable X1 and the second variable X2, and the original inversion objective function is transformed into a multi-constrained optimization problem:

[0164]

[0165] Where, f η (r) represents a multi-constrained optimization problem.

[0166] S202. Using the Lagrange multiplier method to introduce Lagrange parameters λ1 and λ2, adding constraints to the inversion objective function, and obtaining an unconstrained augmented Lagrangian inversion objective function.

[0167] For example, by using the Lagrange multiplier method to introduce the Lagrange parameters λ1 and λ2, the constraints can be added to the inversion objective function, and the inversion objective function in the unconstrained augmented Lagrangian form is obtained as follows:

[0168]

[0169] Where, L η，κ (r, X1, X2, λ1, λ2) represents the inversion objective function in the unconstrained augmented Lagrangian form.

[0170] S203: Define a first dual variable and a second dual variable to obtain an unconstrained inversion objective function.

[0171] For example, the first dual variable and the second dual variable are defined as Finally, the unconstrained inversion objective function is obtained as:

[0172]

[0173] Where, L η，κ (r, X1, X2, Z1, Z2) represents the unconstrained inversion objective function; κ is a penalty parameter used to control the convergence speed.

[0174] S204. Iteratively solve the variable to be solved, the first variable, the second variable, the first dual variable, and the second dual variable in sequence until the convergence condition is met, and output the longitudinal wave velocity, the shear wave velocity, and the density.

[0175] After the changes in steps S201-S203, the original single-variable complex optimization problem of the objective function is transformed into a simple multi-variable optimization problem. The optimal solution can be obtained by iteratively solving the variables r, X1, X2, Z1, and Z2 in sequence. In the i-th loop iteration, the optimization solution for each variable is as follows:

[0176] The quantity to be solved r is optimized by the following method:

[0177] Taking r as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0178]

[0179] Where, represents the first variable solved in the previous iteration, represents the second variable solved in the previous iteration, represents the first dual variable solved in the previous iteration, represents the second dual variable solved in the previous iteration;

[0180] The process of solving the single variable optimization problem in equation (29) using the least squares method is as follows:

[0181]

[0182] Where T represents matrix transpose; r i+1 Represents the variable to be solved after the i-th iterative update;

[0183] The first variable X1 is solved by optimizing as follows:

[0184] Taking X1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0185]

[0186] The process of solving the single variable optimization problem in equation (31) using the least squares method is as follows:

[0187]

[0188] Where I is the identity matrix of the same dimension as X1;

[0189] The second variable X2 is solved by optimizing as follows:

[0190] Taking X2 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0191]

[0192] The process of solving the single variable optimization problem in equation (33) using the least squares method is as follows:

[0193]

[0194] The first dual variable Z1 is solved by optimizing as follows:

[0195] Taking Z1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0196]

[0197] The process of solving the single variable optimization problem in equation (35) using the least squares method is as follows:

[0198] Z1=Z1+κ(X1-K1r i+1 ) (36)

[0199] The first dual variable Z2 is solved by optimizing as follows:

[0200] Taking Z1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as:

[0201]

[0202] The process of solving the single variable optimization problem in equation (37) using the least squares method is as follows:

[0203] Z2=Z2+κ(X2-K2r i+1 ) (38).

[0204] Given the initial data and setting the constraint parameter values, repeat the above iterative solution process and determine whether the convergence conditions are met. If the conditions are met, the loop is exited and the optimal inversion result of the three parameters is obtained.

[0205] In an exemplary embodiment, the pre-stack seismic wide-angle AVO inversion method based on the second-order approximation can be implemented by the following steps:

[0206] Step 1: Input / acquire angle gather data D(θ), angle wavelet w(θ), incident angle θ; low-frequency model data L, model constraint parameters η, η P , η S , η ρ ; Penalty parameter κ, convergence error tol, maximum number of iterations M.

[0207] Step 2: Use the angle wavelet w(θ) to construct the wavelet matrix W(θ), combine the incident angle θ to construct the observation matrices G1 and G2, and construct the integral matrix C, with seismic trace = 1.

[0208] Step 3: If trace<t_Tota, obtain single-channel data d(θ)=D(θ, trace) and single-channel low-frequency model ε=L(trace).

[0209] Step 4, initialization: i = 0, r 0 =ε、X1 0 =K1r, X2 0 =K2r, Z1 0 =0, Z2 0 =0.

[0210] Step 5, if \(i < M\), then sequentially execute Step 6 to Step 9.

[0211] Step 6, update the reflection coefficient \(r\) using Equation (30) i+1

[0212] Step 7, sequentially update the iteration parameters using Equation (32) and Equation (34)

[0213] Step 8, sequentially update the dual variables using Equation (36) and Equation (38)

[0214] Step 9, determine whether the convergence condition is satisfied

[0215] 1) If the condition is satisfied, record the result \(r = r\) i+1 , and exit the loop;

[0216] 2) If the condition is not satisfied, then \(i = i + 1\), and return to Step 5.

[0217] Step 10, calculate the P-wave velocity \(\alpha\), S-wave velocity \(\beta\), and density \(\rho\) using the final inversion result \(r\).

[0218] To verify the correctness and effectiveness of the method proposed in this paper, this embodiment intercepted part of the data ([[]] Figure 5 ) in the Marmousii model, and carried out a comparative analysis of wide-angle pre-stack AVO inversion and traditional pre-stack AVO inversion. Figure 5 In, (a), (b), and (c) are the P-wave velocity, S-wave velocity, and density models respectively; the target layer in the model is a sandstone reservoir (at the black arrow), and the P-wave velocity and density of this gas layer both show obvious low values, while the S-wave velocity remains basically unchanged. The low-frequency model required for inversion is as Figure 6 shown, where (a), (b), and (c) in the figure are the low-frequency models of the P-wave velocity, S-wave velocity, and density respectively.

[0219] Bringing the original model data into the Zoeppritz equation can obtain the angle-dependent P-wave reflection coefficient, and then convolving it with a 20 Hz Ricker wavelet can obtain the angle gather data required for pre-stack inversion ([[]] Figure 7 and 8 ). Figure 7 is the angle gather required for conventional AVO inversion. In the figure, (a), (b), and (c) are the synthetic seismograms with incident angles of 5°, 15°, and 25° respectively; Figure 8 is the angle gather required for wide-angle AVO inversion. In the figure, (a), (b), and (c) are the synthetic seismograms with incident angles of 8°, 24°, and 40° respectively.

[0220] The prestack AVO three-parameter inversion was carried out using the traditional method and the new method. The inversion results are shown in the following table. Figure 9 and Figure 10 shown. Figure 9 This is the traditional pre-stack AVO inversion. Figures (a), (b), and (c) are the inversion results of P-wave velocity, S-wave velocity, and density, respectively. Figure 10 This is the new method for pre-stack AVO inversion. Figures (a), (b), and (c) are the inversion results of P-wave velocity, S-wave velocity, and density, respectively.

[0221] Comparative Analysis Figure 9 and 10 It can be seen that: 1) The traditional AVO inversion method produces a ghost reflection at the bottom of the gas sandstone reservoir ( Figure 9 Black arrow), the inversion method proposed in this paper can better identify the reservoir; 2) The shear wave velocity of the original model in the gas sandstone reservoir is about 540m / s ( Figure 5 (b) black arrow), the inversion results of the proposed method are basically consistent with the model ( Figure 10 (b) black arrow), while the inversion results of the traditional method are quite different from it ( Figure 9 3) Compared with the traditional AVO inversion results, the proposed inversion method also has better inversion effect for thin layers. In summary, the inversion accuracy of the proposed method is better than that of the traditional AVO inversion results.

[0222] In order to further verify the noise resistance of the proposed method, this paper extracts the 40th three-parameter curve of the model ( Figure 12 The black solid line in the middle) and convolving it with the 20Hz Ricker wavelet, we get the single-channel angle gather, as shown in Figure 11 (a), and add 5% ( Figure 11 (b)) and 10% ( Figure 11 The noise in (b) is used to obtain the noisy single-channel angle gather.

[0223] Blue Earthquake Record ( Figure 11 (a) is used for traditional AVO inversion to obtain the three-parameter inversion results ( Figure 12 (a) blue dashed line); similarly, the red earthquake record ( Figure 11 The three-parameter inversion results obtained by the new AVO inversion method in (a) are as follows Figure 12 (a) The red solid line. Similarly, we can get the noise 5% ( Figure 12 (b)) and 10% ( Figure 12 Inversion results obtained by two inversion methods in (c).

[0224] from Figure 12As can be seen from (a), the prestack inversion results of the new method are generally better than those of the traditional method; for the inversion results of the angle gathers containing noise ( Figure 12 Figures (b) and (c) show that the inversion results of the proposed method are stable and have strong noise resistance.

[0225] In order to further verify the effectiveness and practicality of the proposed method, this paper selected the actual seismic data of an oil field as the test data ( Figure 13 and 14 ). Figure 13 The partial angle gather stacking data required for traditional pre-stack AVO inversion. (a), (b), and (c) are partial angle gather stacking data with incident angles of 3° to 8° (center angle is 5°), 15°, and 25°, respectively. Figure 14 These are the partial angle gather stacking data required for the pre-stack AVO inversion method proposed in this paper. In the figure, (a), (b), and (c) are the partial angle gather stacking data with incident angles of 8°, 24°, and 40°, respectively; the black straight line in the figure represents the specific position of well A in the seismic profile.

[0226] use Figure 13 and 14 The angle gather data in the pre-stack AVO inversion were carried out using the traditional method and the new method, and the inversion results correspond to Figure 15 and 16 . Figure 15 This is the traditional pre-stack AVO inversion method. Figures (a), (b), and (c) are the inversion results of P-wave velocity, S-wave velocity, and density, respectively. Figure 16 This is the new method for pre-stack AVO inversion. Figures (a), (b), and (c) are the inversion results of P-wave velocity, S-wave velocity, and density, respectively.

[0227] Comparative Analysis Figure 15 and 16 It can be seen that: 1) The inversion results of the new method are in good agreement with the actual logging interpretation results ( Figure 15 (b), (c) and Figure 16 (b) and (c) black arrows); 2) Compared with the traditional AVO inversion results, the newly proposed inversion method has better inversion effect for thin layers ( Figure 15 3) Comparing the density inversion results, it can be seen that the traditional AVO inversion method will produce strong ghost reflections, while the new method has higher vertical and horizontal resolution ( Figure 15 (c) and Figure 16 In summary, compared with the traditional AVO inversion results, the inversion results of the inversion method proposed in this paper have higher accuracy and resolution.

[0228] In summary, rock elastic parameters are important parameters reflecting subsurface reservoir information, and prestack seismic AVO inversion is one of the primary methods for obtaining these parameters. This paper constructs a prestack inversion objective function using a derived second-order approximation for seismic plane waves at the solid / solid interface and solves it using the alternating direction method of multipliers (ADMM). Ultimately, high-resolution inversion of P-wave velocity, S-wave velocity, and density is achieved. Compared with traditional prestack AVO inversion methods, the second-order approximation employed in this paper achieves higher precision in reflection coefficients at high angles of incidence. Therefore, prestack AVO inversion can more fully utilize high-angle gather information, thereby improving the inversion resolution of the three parameters (particularly S-wave velocity and density). The proposed method provides theoretical and methodological support for improving the prediction accuracy of complex oil and gas reservoirs and more effectively advancing their exploration and development.

[0229] The above embodiments are only used to illustrate the present invention, and are not intended to limit the present invention. Ordinary technicians in the relevant technical field may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent technical solutions also fall within the scope of the present invention. The scope of patent protection of the present invention should be defined by the claims.

Claims

1. A pre-stack seismic wide-angle AVO inversion method based on a second-order approximation, characterized in that: The method comprises: The inversion objective function is constructed using the second-order approximation of seismic plane waves on the solid / solid medium horizontal interface; The alternating direction multiplier method is used to solve the inversion objective function to achieve high-resolution inversion of P-wave velocity, S-wave velocity and density.

2. The method according to claim 1, wherein The second-order approximate formula of the seismic plane wave on the solid / solid medium horizontal interface is expressed as: R(θ)=R1(θ)+R2(θ) (1) Where: R(θ) is the second-order approximation of the seismic plane wave on the solid / solid medium horizontal interface, R1(θ) is the first-order term of the reflection coefficient, R2(θ) is the second-order term of the reflection coefficient; θ is the average incident angle of the upper and lower layers; α, β, ρ are the average P-wave velocity, S-wave velocity and density of the upper and lower layers, respectively; Δα, Δβ and Δρ represent the differences in P-wave velocity, S-wave velocity and density between the lower and upper layers, respectively; c1(θ) is the control coefficient of the first-order P-wave velocity reflection coefficient term, c2(θ) is the control coefficient of the first-order S-wave velocity reflection coefficient term, c3(θ) is the control coefficient of the first-order density reflection coefficient term, x1(θ) is the control coefficient of the second-order S-wave velocity reflection coefficient term, x2(θ) is the control coefficient of the combined term of the second-order S-wave velocity reflection coefficient and density reflection coefficient, x3(θ) is the control coefficient of the second-order density reflection coefficient term, and Y(θ) and H(θ) are the influencing factors of the second-order reflection coefficient.

3. The method according to claim 2, wherein Using the second-order approximation of seismic plane waves on the solid / solid medium horizontal interface, the inversion objective function is constructed as follows: Describing the angle gather as a form of convolution of an angle wavelet and an angle reflection coefficient sequence; wherein the form of convolution of the angle wavelet and the angle reflection coefficient sequence includes seismic responses generated by first-order and second-order terms of the reflection coefficient; Based on the seismic response generated by the first-order term of the reflection coefficient, a set of equations is established by stacking data from several partial angle gathers. Based on the seismic response generated by the second-order term of the reflection coefficient, a relationship is established by stacking data from several partial angle gathers. According to the equations and relations, a three-parameter inversion equation based on a second-order approximation of the reflection coefficient is obtained; Under the weak elasticity assumption, a low-frequency model constraint equation is obtained according to the relationship between the three-parameter reflection terms of the three-parameter inversion equation and the logarithms of the corresponding parameters; An inversion objective function is constructed according to the three-parameter inversion equation and the low-frequency model constraint equation.

4. The method according to claim 3, wherein The convolution of the angle wavelet and the angle reflection coefficient sequence is expressed as: d(θ)=d1(θ)+d2(θ)=w(θ)*R1(θ)+w(θ)*R2(θ) (10) d1(θ)=c1(θ)w(θ)*r P +c2(θ)w(θ)*r S +c3(θ)w(θ)*r ρ (11) d2(θ)=x1(θ)w(θ)*r S ·r S +x2(θ)w(θ)*r S ·r ρ +x3(θ)w(θ)*r ρ ·r ρ (12) Where d1(θ) and d2(θ) represent the seismic responses generated by the first-order and second-order terms of the reflection coefficient, respectively; r P 、r S and r ρ Represent the reflection coefficient sequences of P-wave velocity, S-wave velocity and density respectively; * represents the convolution operator; . represents the Hadamard product operator; w(θ) represents the angle wavelet; d(θ) represents the angle gather.

5. The method according to claim 4, wherein The system of equations is expressed as: d1=G1r(14) Where, d1(θ i ) represents the seismic response generated by the first-order term of the i-th angle reflection coefficient, c1(θ i ) represents the control coefficient of the first-order longitudinal wave velocity reflection coefficient term at the i-th angle, c2(θ i ) represents the control coefficient of the first-order shear wave velocity reflection coefficient term at the i-th angle, c3(θ i ) represents the control coefficient of the first-order density reflection coefficient term at the i-th angle, W(θ i ) represents the wavelet matrix of the i-th angle, d1 represents the seismic response generated by the first-order term of the reflection coefficient, G1 represents the wavelet kernel matrix of the first-order term of the reflection coefficient, and r represents the quantity to be solved; The relationship is expressed as: Where, W(θ i ) represents the wavelet matrix of the i-th angle, I is related to r S The identity matrix of the same dimension, x1(θ i ) represents the control coefficient of the second-order shear wave velocity reflection coefficient term at the i-th angle, x2(θ i ) represents the control coefficient of the joint term of the second-order shear wave velocity reflection coefficient and density reflection coefficient at the i-th angle, x3(θ i ) represents the control coefficient of the second-order density reflection coefficient term of the i-th angle, K1 and K2 represent the elementary transformation matrices, d2 represents the seismic response generated by the second-order reflection coefficient term, and G2 represents the wavelet kernel matrix of the second-order reflection coefficient term; Substituting equation (16) into equation (15), we can obtain the simplified form of the relationship, which is expressed as: d2=G2(K1r)·(K2r) (17) Combining equations (14) and (17), we obtain the three-parameter inversion equation based on the second-order approximation of the reflection coefficient, which is expressed as: d=d1+d2=G1r+G2(K1r)·(K2r) (18) Where d represents the earthquake record.

6. The method according to claim 5, characterized in that, under the weak elasticity assumption, a low-frequency model constraint equation is obtained according to the relationship between the three-parameter reflection terms of the three-parameter inversion equation and the logarithms of the corresponding parameters, and an inversion objective function is constructed according to the three-parameter inversion equation and the low-frequency model constraint equation, comprising: Under the weak elasticity assumption, the relationship between the three-parameter reflection term and the logarithm of the corresponding parameter is expressed as: e P =C′r P ,he S =C′r S ,he ρ =C′r ρ (19) Where α0, β0, and ρ0 represent the initial longitudinal wave velocity, shear wave velocity, and density, respectively. i , β i , ρ i are the P-wave velocity, S-wave velocity and density of the i-th layer respectively, C′ is the integration matrix; Combining equations (19) to (22), we get the low-frequency model constraint equation: ε=Cr (24) Where η P ,η S and η ρ Respectively represent the model constraint parameters of P-wave velocity, S-wave velocity and density, which are used to adjust the model constraint proportions of P-wave velocity, S-wave velocity and density; ε represents the low-frequency model; Combining equations (18) and (24), we construct the following inversion objective function f(r): Where η>0 is the low-frequency model constraint parameter, which determines the low-frequency model term and fidelity Relative contribution to the inversion results.

7. The method according to claim 6, wherein The alternating direction multiplier method is used to solve the inversion objective function to achieve high-resolution inversion of P-wave velocity, S-wave velocity, and density, including: By introducing the first variable and the second variable into the inversion objective function, the process of optimizing K1r and K2r is transformed into optimizing the first variable and the second variable, and the inversion objective function is transformed into a multi-constrained optimization problem; The Lagrangian parameters λ1 and λ2 are introduced by using the Lagrangian multiplier method, and the constraints are added to the inversion objective function to obtain the inversion objective function in the form of unconstrained augmented Lagrangian. Define the first dual variable and the second dual variable to obtain the unconstrained inversion objective function; The variable to be solved, the first variable, the second variable, the first dual variable and the second dual variable are solved iteratively in sequence until the convergence condition is met, and the longitudinal wave velocity, shear wave velocity and density are output.

8. The method according to claim 7, wherein The multi-constrained optimization problem is expressed as: Where, f η (r) represents a multi-constrained optimization problem, X1 and X2 represent the first and second variables, respectively; The inversion objective function of the unconstrained augmented Lagrangian form is expressed as: Where, L η，κ (r, X1, X2, λ1, λ2) represents the inversion objective function of the unconstrained augmented Lagrangian form; Define the first dual variable and the second dual variable, expressed as The obtained unconstrained inversion objective function is expressed as: Where, L η，κ (r, X1, X2, Z1, Z2) represents the unconstrained inversion objective function; κ is a penalty parameter used to control the convergence speed.

9. The method according to claim 8, wherein In the process of iteratively solving the quantity to be solved, the first variable, the second variable, the first dual variable and the second dual variable in the i-th loop: The quantity to be solved r is optimized by the following method: Taking r as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as: Where, represents the first variable solved in the previous iteration, represents the second variable solved in the previous iteration, represents the first dual variable solved in the previous iteration, represents the second dual variable solved in the previous iteration; The process of solving the single variable optimization problem in equation (29) using the least squares method is as follows: Where T represents matrix transpose; r i+1 Represents the variable to be solved after the i-th iterative update; The first variable X1 is solved by optimizing as follows: Taking X1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as: The process of solving the single variable optimization problem in equation (31) using the least squares method is as follows: Where, is the identity matrix of the same dimension as X1; The second variable X2 is solved by optimizing as follows: Taking X2 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as: The process of solving the single variable optimization problem in equation (33) using the least squares method is as follows: The first dual variable Z1 is optimized and solved as follows: Taking Z1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as: The process of solving the single variable optimization problem in equation (35) using the least squares method is as follows: Z1=Z1+κ(X1-K1r i+1 ) (36) The first dual variable Z2 is solved by optimizing as follows: Taking Z1 as a variable and other parameters as known quantities, the unconstrained inversion objective function is simplified as: The process of solving the single variable optimization problem in equation (37) using the least squares method is as follows: Z2=Z2+κ(X2-K2r i+1 ) (38)。 10. The method according to claim 9, wherein The convergence condition is Here, tol represents the convergence error.

Images