A high-precision joint pre-stack inversion method of P-wave and S-wave in parameter domain

Abstract

The present application relates to oil and gas geophysical exploration technology field, more particularly, to a kind of high-precision ray parameter domain P-S wave joint prestack inversion method, it includes the following steps: S1: obtaining the PP wave prestack ray parameter domain gather and PS wave prestack ray parameter domain gather for inversion.S2: establishing PP wave subwave matrix and PS wave subwave matrix.S3: obtaining the wellhead position of target study area's longitudinal wave velocity, transverse wave velocity and density, constructs inversion initial model.S4: establishing L 1‑2 Norm regularization sparse constraint target function.S5: solving target function.S6: output inversion result.The present application overcomes the problem that inversion result is not accurate enough in prior art by using L 1‑2 Norm regularization sparse constraint, and provides a kind of solving algorithm of target function using L 1‑2 Norm sparse constraint, more accurate inversion result can be obtained.

Description

Technical Field

[0001] This invention relates to the field of oil and gas geophysical exploration technology, and more specifically, to a high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method. Background Technology

[0002] In oil and gas exploration and development, information such as subsurface lithology, fluid properties, and porosity is crucial for reservoir prediction and fluid identification. This information can be used for prediction using elastic parameters such as P-wave and S-wave velocities and density. The primary method for obtaining these subsurface elastic parameters is pre-stack seismic inversion. Pre-stack inversion utilizes the characteristic of reflection amplitude varying with offset to invert subsurface elastic parameters from seismic records at different offsets, aiding in the identification of subsurface lithology and fluid characteristics and improving reservoir prediction accuracy.

[0003] In seismic inversion, noise in seismic records and inaccuracies in wavelet forward modeling operators can lead to multiple solutions and ill-posed problems. A common solution to this type of ill-posed problem is to use regularization methods for constraint solving. This involves adding appropriate prior constraints to the inversion problem, transforming it into a well-posed problem before solving, ensuring the inversion results are consistent with the prior constraint characteristics. In seismic inversion, the subsurface strata can be considered to exhibit a sparse distribution; therefore, sparsity regularization is frequently used to ensure that the inversion results are sparse and conform to stratigraphic patterns.

[0004] In terms of regularization constraints, sparsity norms are mainly used to obtain sparse solutions. The L0 norm is the best indicator of sparsity; however, optimization problems involving the L0 norm are NP-hard and difficult to solve. The L1 norm, as a convex approximation of the L0 norm, is widely used in seismic inversion. The objective function constrained by the L1 norm can be expressed as... ,in, It is a regularization parameter that controls the degree of influence of the constraint terms on the inversion. If the value is too small, the inversion results will not be sparse enough, which does not conform to the stratigraphic assumption; If the L1 norm is too large, the inversion results will be very sparse, resulting in the loss of a lot of information. As its application becomes more widespread, the L1 norm has gradually revealed its shortcomings, such as the difficulty in obtaining weak reflections, leading to inaccurate inversion results. Summary of the Invention

[0005] To address the problem that the inversion results are not accurate enough when using L1 norm for regularization constraints in the existing technology, this invention provides a high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method, which can effectively solve the ill-posedness in the inversion process and obtain more accurate inversion results.

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

[0007] A high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method includes the following steps:

[0008] S1: Collect pre-stack migration gathers for PP waves and PS waves in the target study area, and perform gather transformation on these gathers to obtain pre-stack ray parameter domain gathers for PP waves and PS waves used for inversion. The PP wave pre-stack ray parameter domain gathers represent P-wave pre-stack ray parameter domain gathers, and the PS wave pre-stack ray parameter domain gathers represent S-wave pre-stack ray parameter domain gathers; that is, transform the S-wave pre-stack ray parameter domain gathers. Generally, pre-stack inversion is performed in the angular domain, requiring the conversion of seismic data from pre-stack migration gathers to pre-stack angular gathers. The calculation of angular gathers requires prior knowledge of strata velocity information, making the entire inversion process dependent on an accurate velocity model. In contrast, the calculation of ray parameter domain gathers does not require prior velocity information, offering the advantage of being entirely data-driven, and ray parameter domain gathers more intuitively reflect the propagation path of seismic waves. Therefore, this step chooses to perform pre-stack inversion in the ray parameter domain.

[0009] S2: Based on the pre-stack ray parameter domain gathers of the PP wave and the pre-stack ray parameter domain gathers of the PS wave, extract the PP seismic wavelet and the PS seismic wavelet for inversion, and establish the PP wavelet matrix and the PS wavelet matrix respectively.

[0010] S3: Obtain the P-wave velocity, S-wave velocity, and density at the wellhead location of the target study area, and then calculate the logarithms of the P-wave velocity, S-wave velocity, and density respectively. Construct an initial inversion model based on the logarithms of the P-wave velocity, S-wave velocity, and density and the seismic horizon information.

[0011] S4: Establish L 1-2 The objective function of norm regularization sparse constraint is to substitute the PP wavelet matrix, the PS wavelet matrix, and the inversion initial model into the L... 1-2 The objective function of norm regularization sparse constraints.

[0012] S5: Solve for the objective function.

[0013] S6: Output the logarithmic inversion results of the P-wave velocity, S-wave velocity and density, and then perform exponential calculation processing on the logarithmic inversion results of the P-wave velocity, S-wave velocity and density of the target study area to obtain the P-wave velocity, S-wave velocity and density of the subsurface. Thus, the inversion results are obtained.

[0014] In the above technical solution, L 1-2The norm is represented by the difference between the L1 norm and the L2 norm. It possesses better sparsity properties than the L1 norm, and sparsity reflects the distribution of boundaries and anomalies in underground structures. In the inversion, the conventional L1 norm constraint is replaced with L... 1-2 Norm constraints can reduce the complexity of the inversion model, improve interpretability and noise resistance of the inversion results, and thus obtain more accurate inversion results.

[0015] Preferably, in step S1, the approximate formula is rewritten in ray parameter form, and joint inversion is performed in the ray parameter domain. The formula for calculating the longitudinal wave velocity reflection coefficient when the ray parameter is p is:

[0016] ,

[0017] The formula for calculating the reflection coefficient of the converted shear wave velocity when the ray parameter is p is:

[0018] ,

[0019] These two formulas reflect the relationship between the magnitude of the reflection coefficient and the ray parameters, where, These are ray parameters. , and It is the average value of the longitudinal wave velocity, transverse wave velocity, and density of the upper and lower media. It is the difference in longitudinal wave velocity between the upper and lower media. It is the difference in transverse wave velocity between the upper and lower media. It is the density difference between the upper and lower media.

[0020] Preferably, for ease of calculation, in step S1, the formulas for calculating the reflection coefficients of the PP wave and the PS wave are expressed as a linear sum of the reflection coefficients of three elastic parameters (longitudinal wave velocity reflection coefficient, transverse wave velocity reflection coefficient, and density reflection coefficient), which is expressed as:

[0021] ,

[0022] in, , and These are the reflection coefficient vectors for P-wave velocity, S-wave velocity, and density, calculated using P-wave velocity, S-wave velocity, and density, respectively. , , ;and , , , , It is related to the ray parameters The correlation coefficients are expressed as follows:

[0023] ,

[0024] Therefore, in the ray parameter domain, when the ray parameter is At that time, the earthquake occurred in the first... The longitudinal wave reflection coefficient and the converted transverse wave reflection coefficient at each sampling point are expressed as follows:

[0025] ,

[0026] The reflection coefficients of longitudinal and transverse waves under different ray parameters are represented by matrices as follows:

[0027] ,

[0028] in, The ray parameters are The longitudinal wave reflection coefficient vector at time, To convert the transverse wave reflection coefficient vector; It is a diagonal coefficient matrix. , , and It is also a diagonal coefficient matrix; , and It is a vector of longitudinal wave velocity reflection coefficient, transverse wave velocity reflection coefficient, and density reflection coefficient calculated using longitudinal wave velocity, transverse wave velocity, and density, respectively.

[0029] Preferably, in step S3, according to the convolution forward model, the seismic record can be represented as the convolution of the reflection coefficient and the seismic wavelet. For seismic gathers with different pre-stack ray parameters, the forward model is represented by a matrix as follows:

[0030] ,

[0031] in, The ray parameters are P-wave seismic data vector at time, It is a transformation of shear wave seismic data vectors; It corresponds to the first The longitudinal wavelet matrix with ray parameters, It is a transformation transverse wavelet matrix;

[0032] The noisy forward convolution model is represented as:

[0033] ,

[0034] in, These are P-wave and S-wave seismic data vectors under different ray parameters; It is a diagonal matrix of wavelets composed of longitudinal and transverse wavelet matrices under different ray parameters; It is a diagonal matrix , , , and The coefficient matrix formed; It is the overall elastic parameter reflection coefficient vector composed of the longitudinal wave velocity, transverse wave velocity, and density reflection coefficient vector; It is a random noise vector. This model can effectively link P-wave and S-wave seismic data, laying the foundation for joint inversion.

[0035] Based on L 1-2 The objective function of the norm constraint can be expressed as: ,in, It is also a regularization parameter. The parameter to be inverted in this formula... It is the reflection coefficient of the underground elastic parameters. The conventional inversion method is to obtain... Then, the underground elastic parameters are obtained by combining the initial model with trace integrals. However, this indirect method has certain problems; if some reflection coefficients are not accurately obtained, errors will accumulate. Therefore, to overcome this problem, this invention utilizes the idea of ​​TV regularization and adds initial model constraints to directly invert the underground elastic parameters. TV regularization uses a difference matrix to rewrite the terms to be inverted, so as to... For example:

[0036] ,

[0037] make ,vector It can then be expressed as ,in, It is a first-order difference matrix, represented as:

[0038] ,

[0039] Therefore, based on the above calculation formula and combined with the initial model constraints, the objective function is to solve for the underground elastic parameters (P-wave velocity, S-wave velocity, and density) and reflection coefficient. This becomes a logarithm for solving for the elastic parameters (P-wave velocity, S-wave velocity, and density) underground. ,right Exponential calculations can directly yield the elastic parameters of the subsurface (longitudinal wave velocity, transverse wave velocity, and density).

[0040] In summary, preferably, in step S4, the expression of the objective function is:

[0041] ,

[0042] in, , It is a regularization parameter that controls the sparsity of the inversion. It is a regularization parameter that controls the initial model weights; This is the initial model, which can be created using actual well logging data; m is the logarithm of the subsurface elastic parameter. It is a first-order difference matrix, represented as:

[0043] .

[0044] Preferably, in step S5, two auxiliary variables are introduced to construct an augmented Lagrangian form, and the objective function is solved using the Convexity Difference Algorithm (DCA) + Alternating Multiplier Algorithm (ADMM). The optimal solution and the expressions for the two auxiliary variables are as follows:

[0045] ,

[0046] in, , These are introduced auxiliary variables. It is a penalty parameter that controls the convergence speed of the iteration; sthresh represents the calculation of the soft threshold. This represents the optimal solution; after obtaining the optimal solution, an exponential operation is performed. Thus, the P-wave velocity, S-wave velocity, and density of the subsurface of the target study area can be obtained.

[0047] Preferably, in step S1, the following is adopted: The target study area is transformed using either the transformation processing method or the bending ray tracing method to convert the PP wave pre-stack offset gather and the PS wave pre-stack offset gather, thereby obtaining the PP wave pre-stack ray parameter domain gather and the PS wave pre-stack ray parameter domain gather.

[0048] Preferably, in step S3, the well logging data in the working area of ​​the PP wave pre-stack ray parameter domain gather and the PS wave pre-stack ray parameter domain gather are respectively subjected to low-pass filtering to obtain their low-frequency trends, and the inversion initial model is obtained by interpolation based on the seismic horizon information, which is beneficial to improving the accuracy and reliability of pre-stack inversion.

[0049] Preferably, in step S5, the model parameters of the initial inversion model are updated; before step S6, the maximum number of iterations, the maximum value of the termination update change, and the maximum value of the total error are set; in step S6, the changes in model parameters before and after the initial inversion model update are calculated, and the total error is calculated based on the pre-stack ray parameter domain gather of the PP wave, the synthetic record of the PP wave, the pre-stack ray parameter domain gather of the PS wave, and the synthetic record of the PS wave. Then, it is determined whether the number of iterations is greater than or equal to the maximum number of iterations, or whether the change is greater than or equal to the maximum value of the termination update change and the total error is greater than or equal to the maximum value of the total error. If so, the inversion result is output; otherwise, steps S5 to S6 are repeated. It can be understood that the number of times steps S5 to S6 are repeated is the number of iterations. Outputting the inversion result only after the iteration termination condition is met can make the model obtained by each update closer to the real underground model, thereby obtaining a more stable and accurate inversion result.

[0050] Preferably, the method further includes step S7: adding noise with different signal-to-noise ratios to the pre-stack ray parameter domain gathers of the PP wave and the pre-stack ray parameter domain gathers of the PS wave, then calculating the mean relative error (MRE) and correlation coefficient (CC) between the real model and the initial inversion model, and analyzing the stability and accuracy of the inversion results based on the mean relative error (MRE) and the correlation coefficient (CC).

[0051] The beneficial effects of this invention: Introducing L into the inversion process 1-2 Norm, replacing the regular L1 norm constraint with L 1-2 Norm constraints, and provide a method using L 1-2 The algorithm for solving the objective function of norm sparse constraints can obtain more accurate inversion results. This method has good noise resistance. Even with low signal-to-noise ratio, the average relative error and correlation coefficient between the inversion results and the real model can be maintained at a high level, which has high application value. Attached Figure Description

[0052] Figure 1 This is a flowchart of a high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method;

[0053] Figure 2 This is a schematic diagram of the wave transformation when a longitudinal wave is incident;

[0054] Figure 3 It is a plane contour map of the L0 norm;

[0055] Figure 4 It is a plane contour map of the L1 norm;

[0056] Figure 5 It is a plane contour map of the L2 norm;

[0057] Figure 6 It is L 1-2 Planar contour plot of norm;

[0058] Figure 7 It is a parameter information diagram of a multi-layer uniform model;

[0059] Figure 8 This is a schematic diagram of a multi-layered uniform model;

[0060] Figure 9 This is a schematic diagram of the longitudinal wave gather in the ray parameter domain;

[0061] Figure 10 This is a schematic diagram of the transverse wave gather converted from the ray parameter domain;

[0062] Figure 11 This is a schematic diagram of the initial model obtained by inversion from a multi-layer uniform model.

[0063] Figure 12 This is a schematic diagram of the inversion results under noise-free conditions;

[0064] Figure 13 This is a schematic diagram of the inversion results under the condition of a signal-to-noise ratio of 10;

[0065] Figure 14 This is a schematic diagram of the inversion results under the condition of a signal-to-noise ratio of 5;

[0066] Figure 15 This is a graph showing the average relative error and correlation coefficient between the inversion results and the actual model. Detailed Implementation

[0067] The technical solution of the present invention will be further described in detail below through specific embodiments and in conjunction with the accompanying drawings:

[0068] Example 1

[0069] like Figure 1 The high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method shown includes the following steps:

[0070] S1: Collect PP-wave pre-stack migration gathers and PS-wave pre-stack migration gathers for the target study area, and perform gather transformation on the PP-wave and PS-wave pre-stack migration gathers to obtain PP-wave and PS-wave pre-stack ray parameter domain gathers for inversion. Generally, pre-stack inversion is performed in the angular domain, requiring the conversion of seismic data from pre-stack migration gathers to pre-stack angular gathers. The calculation of angular gathers requires prior knowledge of formation velocity information, making the entire inversion process dependent on an accurate velocity model. However, extracting ray parameter domain gathers does not rely on a velocity model, thus avoiding this situation. Figure 2 This demonstrates the phenomenon of converted shear waves generated when incident P-waves encounter elastic interfaces during seismic exploration. and This represents the elastic interface of the formation. Wave represents the incident longitudinal wave. The wave represents a reflected longitudinal wave. This represents a transverse wave. From Figure 2 It can be observed that the reflection angles of reflected P-waves and reflected S-waves at the same underground location are different, but the ray parameters remain unchanged. Therefore, joint inversion in the ray parameter domain avoids angle conversion, and the ray parameter domain gather can more intuitively reflect the propagation path of seismic waves.

[0071] S2: Based on the pre-stack ray parameter domain gathers of PP waves and PS waves, the PP seismic wavelets and PS seismic wavelets used for inversion are extracted respectively, and the PP wavelet matrix and PS wavelet matrix are established respectively.

[0072] S3: Obtain the P-wave velocity, S-wave velocity, and density at the wellhead location in the target study area, then calculate the logarithms of the P-wave velocity, S-wave velocity, and density respectively, and construct the initial inversion model based on the logarithms of the P-wave velocity, S-wave velocity, and density and the seismic horizon information.

[0073] S4: Establish L 1-2 The objective function of norm regularization sparse constraint is to substitute the PP wavelet matrix, PS wavelet matrix, and inversion initial model into L. 1-2 The objective function of norm regularization sparse constraints. For example... Figures 3 to 6 As shown, these figures illustrate plane contour maps with different norm constraints. The contour values ​​decrease closer to the origin. The optimization process of the objective function can be understood as finding the solution where the regularization term is minimized while keeping the least squares term constant; that is, searching for the minimum contour value. Regarding the distribution of contour lines, if the overall value is closer to the origin... shaft and The more sparse the corresponding inversion solution, the more sparsy the axis. It can be observed that the L0 norm is the sparsest of these constraints, with its contour lines all distributed along the axis. shaft and On the axis. Compared to the L1 norm, L 1-2 The minimum contour line of the norm is closer to shaft and The axis, as a whole, is closer to the L0 norm, so L 1-2 The sparsity of norm constraints is better than that of L1 norm constraints.

[0074] S5: Solve for the objective function.

[0075] S6: Output the logarithmic inversion results of P-wave velocity, S-wave velocity, and density. Then, perform exponential calculation on the logarithmic inversion results of P-wave velocity, S-wave velocity, and density to obtain the P-wave velocity, S-wave velocity, and density in the subsurface of the target study area. This completes the inversion results.

[0076] Further, in step S1, the approximate formula is rewritten in ray parameter representation, and joint inversion is performed in the ray parameter domain. The formula for calculating the longitudinal wave velocity reflection coefficient when the ray parameter is p is:

[0077] ,

[0078] The formula for calculating the reflection coefficient of the transverse wave velocity when the ray parameter is p is:

[0079] ,

[0080] These two formulas reflect the relationship between the magnitude of the reflection coefficient and the ray parameters, where, These are ray parameters. , and It is the average value of the longitudinal wave velocity, transverse wave velocity, and density of the upper and lower media. It is the difference in longitudinal wave velocity between the upper and lower media. It is the difference in transverse wave velocity between the upper and lower media. It is the density difference between the upper and lower media.

[0081] Furthermore, for ease of calculation, in step S1, the formulas for calculating the reflection coefficients of PP waves and PS waves are expressed as a linear sum of the reflection coefficients of three elastic parameters (longitudinal wave velocity coefficient reflection, transverse wave velocity reflection coefficient, and density reflection coefficient), which is expressed as:

[0082] ,

[0083] in, , and These are the reflection coefficient vectors for P-wave velocity, S-wave velocity, and density, calculated using P-wave velocity, S-wave velocity, and density, respectively. , , ;and , , , , It is related to the ray parameters The correlation coefficients are expressed as follows:

[0084] ,

[0085] Therefore, in the ray parameter domain, when the ray parameter is At that time, the earthquake occurred in the first... The longitudinal wave reflection coefficient and the converted transverse wave reflection coefficient at each sampling point are expressed as follows:

[0086] ,

[0087] The reflection coefficients of longitudinal and transverse waves under different ray parameters are represented by matrices as follows:

[0088] ,

[0089] in, The ray parameters are The longitudinal wave reflection coefficient vector at time, To convert the transverse wave reflection coefficient vector; It is a diagonal coefficient matrix. , , and It is also a diagonal coefficient matrix; , and It is a vector of longitudinal wave velocity reflection coefficient, transverse wave velocity reflection coefficient, and density reflection coefficient calculated using longitudinal wave velocity, transverse wave velocity, and density, respectively.

[0090] Furthermore, in step S3, according to the convolution forward model, the seismic record can be represented as the convolution of the reflection coefficient and the seismic wavelet. For seismic gathers with different pre-stack ray parameters, the forward model is represented by a matrix as follows:

[0091] ,

[0092] in, The ray parameters are P-wave seismic data vector at time, It is a transformation of shear wave seismic data vectors; It corresponds to the first The longitudinal wavelet matrix with ray parameters, It is a transformation transverse wavelet matrix;

[0093] The noisy forward convolution model is represented as:

[0094] ,

[0095] in, These are P-wave and S-wave seismic data vectors under different ray parameters; It is a diagonal matrix of wavelets composed of longitudinal and transverse wavelet matrices under different ray parameters; It is a diagonal matrix , , , and The coefficient matrix formed; It is the overall elastic parameter reflection coefficient vector composed of the longitudinal wave velocity, transverse wave velocity, and density reflection coefficient vector; It is a random noise vector. This model can effectively link P-wave and S-wave seismic data, laying the foundation for joint inversion.

[0096] Based on L 1-2 The objective function of the norm constraint can be expressed as: ,in, It is also a regularization parameter. The parameter to be inverted in this formula... It is the reflection coefficient of the underground elastic parameters. The conventional inversion method is to obtain... Then, the underground elastic parameters are obtained by combining the initial model with trace integrals. However, this indirect method has certain problems; if some reflection coefficients are not accurately obtained, errors will accumulate. Therefore, to overcome this problem, this invention utilizes the idea of ​​TV regularization and adds initial model constraints to directly invert the underground elastic parameters. TV regularization uses a difference matrix to rewrite the terms to be inverted, so as to... For example:

[0097] ,

[0098] make ,vector It can then be expressed as ,in, It is a first-order difference matrix, represented as:

[0099] ,

[0100] Therefore, based on the above calculation formula and combined with the initial model constraints, the objective function is to solve for the underground elastic parameters (P-wave velocity, S-wave velocity, and density) and reflection coefficient. This becomes a logarithm for solving for the elastic parameters (P-wave velocity, S-wave velocity, and density) underground. ,right Exponential calculations can directly yield the elastic parameters of the subsurface (longitudinal wave velocity, transverse wave velocity, and density).

[0101] In summary, further, in step S4, the expression for the objective function is:

[0102] ,

[0103] in, , It is a regularization parameter that controls the sparsity of the inversion. It is a regularization parameter that controls the initial model weights; This is the initial model, which can be created using actual well logging data; m is the logarithm of the subsurface elastic parameter. It is a first-order difference matrix, represented as:

[0104] .

[0105] Furthermore, in step S5, two auxiliary variables are introduced to construct the augmented Lagrangian form, and the objective function is solved using the Convexity Difference Algorithm (DCA) + Alternating Multiplier Algorithm (ADMM). The optimal solution and the expressions for the two auxiliary variables are as follows:

[0106] ,

[0107] in, , These are introduced auxiliary variables. It is a penalty parameter that controls the convergence speed of the iteration; sthresh represents the calculation of the soft threshold. This represents the optimal solution; after obtaining the optimal solution, exponential operations are performed. This allows us to obtain the P-wave velocity, S-wave velocity, and density of the subsurface in the target study area.

[0108] The working principle or workflow of this embodiment is as follows: Select a target study area and design a multi-layered homogeneous model. It can be understood that a multi-layered homogeneous model refers to a subsurface model that is spatially divided into multiple horizontal layers, and each layer is assumed to have homogeneous physical properties. The model parameters of this multi-layered homogeneous model are as follows: Figure 7 As shown, a total of 361 sampling points were collected in terms of time depth, with a sampling interval of 2ms. Figure 8 The curves for this model are shown below. The curve on the left represents the longitudinal wave velocity, the curve in the middle represents the transverse wave velocity, and the curve on the right represents the density.

[0109] Forward modeling of the above model using ray parameter domain gathers is performed. Figure 9 and Figure 10The images show the noise-free P-wave gather and the converted S-wave gather, respectively, generated during forward modeling. The seismic wavelet is a Ricker wavelet with a dominant frequency of 30 Hz, and the ray parameter values ​​are set to 0.5 × 10⁻⁶. -4 s / m, 1×10 -4 s / m, 1.5×10 -4 s / m, 2×10 -4 s / m, 2.5×10 -4 s / m, 3×10 -4 s / m. By performing forward modeling of the ray parameter domain gathers of this model, a relatively obvious AVP (amplitude versus ray-parameter) phenomenon can be found. Therefore, pre-stack inversion in the ray parameter domain can be called AVP inversion.

[0110] The inversion initial model obtained in step S4 is as follows Figure 11 As shown. Figure 12 The inversion results are under noise-free conditions. The curve represents the initial inversion model, and both the real model and the inversion results are stepped broken lines. As can be seen from the figure, the two basically overlap, indicating that under noise-free conditions, the inversion results and the real model match very well. The inversion effect of the three elastic parameters is good, which proves the accuracy and feasibility of the inversion method.

[0111] The beneficial effect of this embodiment is that it introduces L into the inversion process. 1-2 Norm, replacing the regular L1 norm constraint with L 1-2 Norm constraints, and provide a method using L 1-2 Algorithms for solving the objective function of norm sparse constraints can yield more accurate inversion results.

[0112] Example 2

[0113] This embodiment further supplements steps S1 to S6 based on embodiment 1.

[0114] In step S1, the following is adopted: The target study area is transformed using either the transformation processing method or the bending ray tracing method to obtain the PP wave pre-stack offset gather and the PS wave pre-stack offset gather.

[0115] Furthermore, in step S3, the well logging data in the working area of ​​the PP wave pre-stack ray parameter domain gather and the PS wave pre-stack ray parameter domain gather are subjected to low-pass filtering to obtain their low-frequency trends. Interpolation is performed according to the seismic horizon to obtain the initial inversion model for inversion, which helps to improve the accuracy and reliability of pre-stack inversion.

[0116] Further, in step S5, the model parameters of the initial inversion model are updated. Before step S6, the maximum number of iterations, the maximum value of the termination update change, and the maximum value of the total error are set. In step S6, the changes in model parameters before and after the initial inversion model update are calculated. Additionally, the total error is calculated based on the PP wave pre-stack ray parameter domain gather, the PP wave composite record, the PS wave pre-stack ray parameter domain gather, and the PS wave composite record. Then, it is determined whether the number of iterations is greater than or equal to the maximum number of iterations, whether the change is greater than or equal to the maximum termination update change, and whether the total error is greater than or equal to the maximum total error. If so, the inversion result is output; otherwise, steps S5 to S6 are repeated. It can be understood that the number of times steps S5 to S6 are repeated is the number of iterations. Outputting the inversion result only after the iteration termination condition is met allows each updated model to more closely approximate the actual underground model, resulting in a more stable and accurate inversion result.

[0117] The beneficial effects of this embodiment are the same as those of Embodiment 1.

[0118] Example 3

[0119] Based on Example 2, this embodiment also includes step S7: adding noise with different signal-to-noise ratios to the PP wave pre-stack ray parameter domain gather and the PS wave pre-stack ray parameter domain gather, then calculating the mean relative error (MRE) and correlation coefficient (CC) between the real model and the initial inversion model, and analyzing the stability and accuracy of the inversion results based on the mean relative error (MRE) and correlation coefficient (CC).

[0120] Specifically, to test the noise resistance of the inversion method, different levels of noise were added to the pre-stack gathers, with signal-to-noise ratios of 10 and 5, respectively. The resulting initial inversion model and inversion results are as follows: Figure 13 and Figure 14 As shown. Similarly, Figure 13 and Figure 14 The curves in the diagram represent the initial inversion model, while both the actual model and the inversion results are stepped broken lines. Figure 13 and Figure 14 As can be seen from the results, in noisy conditions, although the inversion effect decreases as the signal-to-noise ratio decreases, it still tends to the real model, proving that the inversion method has good noise resistance.

[0121] Figure 15 The figure shows the average relative error and correlation coefficient between the inversion results and the true model. Figure 15Analysis of the average relative error and correlation coefficient reveals that as noise increases, the average relative error gradually increases while the correlation coefficient decreases. However, overall, both remain at a good level, indicating that the method is robust and can overcome the influence of noise. The inverted P-wave velocity, S-wave velocity, and density are all relatively accurate.

[0122] Other features, working principles, and beneficial effects of this embodiment are the same as those of Embodiment 2.

[0123] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art will recognize that other variations or modifications can be made based on the above description, and it is neither necessary nor possible to exhaustively describe all possible implementations here. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the claims of the present invention.

Claims

1. A high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method, characterized in that, Includes the following steps: S1: Collect PP wave pre-stack offset gathers and PS wave pre-stack offset gathers in the target study area, and perform gather transformation on the PP wave pre-stack offset gathers and PS wave pre-stack offset gathers to obtain PP wave pre-stack ray parameter domain gathers and PS wave pre-stack ray parameter domain gathers for inversion. S2: Based on the pre-stack ray parameter domain gather of the PP wave and the pre-stack ray parameter domain gather of the PS wave, extract the PP seismic wavelet and the PS seismic wavelet for inversion, and establish the PP wavelet matrix and the PS wavelet matrix respectively. S3: Obtain the P-wave velocity, S-wave velocity, and density at the wellhead location of the target study area, and then calculate the logarithms of the P-wave velocity, S-wave velocity, and density respectively. Construct an initial inversion model based on the logarithms of the P-wave velocity, S-wave velocity, and density and the seismic horizon information. S4: Establish L 1-2 The objective function of norm regularization sparse constraint is to substitute the PP wavelet matrix, the PS wavelet matrix, and the inversion initial model into the... L 1-2 The objective function of norm regularization sparse constraints; S5: Solve for the objective function; S6: Output the logarithmic inversion results of the P-wave velocity, S-wave velocity and density, and then perform exponential calculation processing on the logarithmic inversion results of the P-wave velocity, S-wave velocity and density to obtain the P-wave velocity, S-wave velocity and density of the subsurface of the target study area. In step S4, the expression for the objective function is: ， in, ; It is a regularization parameter that controls the sparsity of the inversion. It is a regularization parameter that controls the initial model weights; This is the initial model, which can be created using actual well logging data; m is the logarithm of the subsurface elastic parameters to be solved. It is a first-order difference matrix, represented as: ； In step S5, two auxiliary variables are introduced to construct the augmented Lagrangian form, and the objective function is solved using the Convexity Algorithm (DCA) + Alternating Multiplier Algorithm (ADMM). The optimal solution and the expressions for the two auxiliary variables are as follows: ， in, , These are introduced auxiliary variables. It is a penalty parameter that controls the convergence speed of the iteration; sthresh represents the calculation of the soft threshold. This represents the optimal solution; after obtaining the optimal solution, an exponential operation is performed. Thus, the P-wave velocity, S-wave velocity, and density of the subsurface of the target study area can be obtained.

2. The high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method according to claim 1, characterized in that, In step S1, inversion is performed in the ray parameter domain. The formula for calculating the longitudinal wave velocity reflection coefficient when the ray parameter is p is: ， The formula for calculating the reflection coefficient of the converted shear wave velocity when the ray parameter is p is: ， in, These are ray parameters. , and These are the average values ​​of the longitudinal wave velocity, transverse wave velocity, and density of the upper and lower media, respectively. It is the difference in longitudinal wave velocity between the upper and lower media. It is the difference in transverse wave velocity between the upper and lower media. It is the density difference between the upper and lower media.

3. The high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method according to claim 2, characterized in that, In step S1, the formulas for calculating the longitudinal wave velocity reflection coefficient and the converted transverse wave velocity reflection coefficient can be linearly expressed as follows: ， in, , and These are the reflection coefficient vectors for P-wave velocity, S-wave velocity, and density, calculated using P-wave velocity, S-wave velocity, and density, respectively. , , ;and , , , , It is related to the ray parameters The correlation coefficients are expressed as follows: ， In the ray parameter domain, when the ray parameter is At that time, the earthquake channel number The longitudinal wave velocity reflection coefficient and the converted transverse wave velocity reflection coefficient at each sampling point are expressed as follows: ， The reflection coefficients of longitudinal and transverse waves under different ray parameters are represented by matrices as follows: ， in, The ray parameters are The longitudinal wave reflection coefficient vector at time, To convert the transverse wave reflection coefficient vector; It is a diagonal coefficient matrix. , , and It is also a diagonal coefficient matrix; , and It is a vector of longitudinal wave velocity reflection coefficient, transverse wave velocity reflection coefficient, and density reflection coefficient calculated using longitudinal wave velocity, transverse wave velocity, and density, respectively.

4. The high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method according to claim 3, characterized in that, In step S3, for seismic gathers with different pre-stack ray parameters, their forward model is represented by a matrix as follows: ， in, The ray parameters are P-wave seismic data vector at time, It is a transformation of shear wave seismic data vectors; It corresponds to the first The longitudinal wavelet matrix with ray parameters, It is a transformation transverse wavelet matrix; The noisy forward convolution model is represented as: ， in, These are P-wave and S-wave seismic data vectors under different ray parameters; It is a diagonal matrix of wavelets composed of longitudinal and transverse wavelet matrices under different ray parameters; It is a diagonal matrix , , , and The coefficient matrix formed; It is the overall elastic parameter reflection coefficient vector composed of the longitudinal wave velocity, transverse wave velocity, and density reflection coefficient vector; It is a random noise vector.

5. The high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method according to claim 1, characterized in that, In step S1, the following is adopted: The target study area is transformed using either the transformation processing method or the bending ray tracing method to convert the PP wave pre-stack offset gather and the PS wave pre-stack offset gather, thereby obtaining the PP wave pre-stack ray parameter domain gather and the PS wave pre-stack ray parameter domain gather.

6. The high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method according to claim 1, characterized in that, In step S3, the well data within the working area of ​​the PP wave pre-stack ray parameter domain gather and the PS wave pre-stack ray parameter domain gather are subjected to low-pass filtering to obtain their low-frequency trends, and the inversion initial model is obtained by interpolation based on the seismic horizon information.

7. The high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method according to claim 1, characterized in that, In step S5, the model parameters of the initial inversion model are updated. Before step S6, the maximum number of iterations, the maximum value of the termination update change, and the maximum value of the total error are set. In step S6, the changes in model parameters before and after the initial inversion model update are calculated. The total error is then calculated based on the pre-stack ray parameter domain gather of the PP wave, the synthetic record of the PP wave, the pre-stack ray parameter domain gather of the PS wave, and the synthetic record of the PS wave. Then, it is determined whether the number of iterations is greater than or equal to the maximum number of iterations, or whether the change is greater than or equal to the maximum value of the termination update change and the total error is greater than or equal to the maximum value of the total error. If so, the inversion result is output; otherwise, steps S5 to S6 are repeated.

8. A high-precision ray parameter domain P-wave and S-wave joint pre-stack inversion method according to any one of claims 1 to 7, characterized in that, The method also includes step S7: adding noise with different signal-to-noise ratios to the pre-stack ray parameter domain gathers of the PP wave and the pre-stack ray parameter domain gathers of the PS wave, then calculating the average relative error and correlation coefficient between the real model and the initial inversion model, and analyzing the stability and accuracy of the inversion results based on the average relative error and the correlation coefficient.

Images