A prestack multi-wave joint inversion method constrained by ultra-differential Laplace blocks

Abstract

The invention provides a pre-stack multi-wave joint inversion method with ultra-differential Laplace block constraints. The method comprises the following steps: firstly, a statistical wavelet is extracted based on seismic data, and the distribution law satisfied by the spatial derivatives of the P-wave velocity, S-wave velocity and density of the work area is statistically analyzed based on the well logging data to determine the selection of parameters in the differential Laplace distribution. Then, an initial P-wave and S-wave velocity and density model that conforms to the sedimentary mode and the seismic phase is established. Then, initial synthetic PP waves and PS seismic data are calculated by convolution wavelets, and the residuals are compared with the real seismic data to calculate the derivatives of the forward operator to the model parameters. Then, based on the Bayesian principle, an ultra-differential Laplace distribution block constraint is added to construct an inversion objective function in the sense of maximum a posteriori probability to obtain an inversion result. Finally, according to the Gauss-Newton method iteration, a final elastic parameter inversion result is obtained to guide the formation interpretation and reservoir prediction, and effectively improve the accuracy and effectiveness of the seismic elastic parameter inversion.

Description

Technical Field

[0001] The invention belongs to the technical field of seismic attribute prediction in oil and gas resource exploration, and in particular relates to a pre-stack multi-wave joint inversion method constrained by an ultra-differential Laplace block. Background Art

[0002] Formation elastic parameters are important attributes that describe formation properties and are the bridge between seismic amplitude information and reservoir physical parameters. Accurate elastic parameters (such as P-wave velocity, S-wave velocity and density, etc.) can not only be used for formation interpretation, but also can be used to obtain physical parameters of the formation (such as porosity, permeability and water saturation, etc.). Therefore, prestack elastic parameter inversion has become a routine workflow in oil and gas exploration. The physical completeness of the direct problem is a necessary condition for the successful solution of the inverse problem. In the prior art, the seismic prestack forward modeling method is generally based on the Zoeppritz equation and its linear approximation, including the Aki-Richard approximation, the Hilterman approximation, the Simth-Gidlow approximation, the Fatti approximation and the Gray approximation. These linear approximations are all proposed under different assumptions, and the formation elastic parameters used in the equations are not the same. Therefore, different formation elastic parameters can be directly inverted using different approximations. Compared with the exact Zoeppritz equation, the inversion based on its linear approximation is mathematically simpler. Due to linear approximation, in all approximate formulas, there is a linear relationship between the reflection coefficient and the formation elastic parameters. Therefore, the inversion based on the linear approximate formula is a linear inversion method; while the exact Zoeppritz equation is a nonlinear equation, so the inversion based on the exact formula is obviously more complicated. At first, the prestack inversion technology was mainly based on the inversion of the Zoeppritz equation approximate formula. However, due to many assumptions of its approximate formula, such as small angle incidence (the incident angle does not exceed 30°), the elastic parameters on both sides of the interface do not change much, etc., the inversion based on the approximate formula is obviously insufficient in accuracy and easy to produce illusions. Therefore, some scholars later studied the prestack inversion method based on the exact Zoeppritz equation. However, the exact Zoeppritz equation and its approximate formula are based on the single interface assumption, without considering transmission loss, conversion waves and interlayer multiple waves, which is obviously inconsistent with the propagation of seismic waves in actual formations, and also leads to many illusions in the inversion results.

[0003] At present, the main prestack inversion techniques are based on the inversion of pure P-waves (PP waves). With the development of four-component acquisition technology, converted wave (PS wave) information has attracted more and more attention from geophysicists. Compared with PP waves, PS waves have a stronger ability to penetrate fluid-containing reservoirs, especially gas cloud areas, and can be better imaged on seismic profiles, so they are more conducive to reservoir inversion. In addition, PS waves are more sensitive to S-wave velocity and density information than PP waves. Therefore, the joint inversion based on PP and PS waves can not only obtain accurate P-wave velocity information, but also more accurately obtain S-wave velocity and density information. In addition, due to the pathological nature of the forward operator, regularization terms are usually required in the inversion to stabilize the inversion results. The selection of regularization terms has a great influence on the accuracy and resolution of the inversion results. The most common regularization methods are Gaussian constraints, Cauchy constraints, and differential Laplace constraints. Compared with Gaussian constraints, the latter two have greatly improved the resolution of the inversion results. Nevertheless, both the Cauchy constraint and the differential Laplace constraint have their own drawbacks, which affect the quality of the final inversion results. Therefore, the seismic pre-stack inversion method in the prior art mainly has the following problems: 1. The forward operator cannot correctly characterize the real propagation process of seismic waves underground, which leads to many illusions in the inversion results; 2. The PP wave is severely attenuated when passing through a fluid-containing reservoir, especially when there is a gas cloud area in the reservoir, so the imaging effect of the reservoir area on the seismic profile is poor, which affects the inversion results. In addition, the PP wave is less sensitive to shear wave velocity and density, resulting in low inversion accuracy of the two; 3. The inversion based on the conventional regularization term cannot well characterize the blocky nature of the formation, resulting in a decrease in the resolution of the inversion and difficulty in characterizing the reservoir boundary. Summary of the invention

[0004] The problem to be solved by the invention is to overcome the defects in the above-mentioned prior art and propose a pre-stack multi-wave joint inversion method constrained by ultra-differential Laplace blocks to effectively improve the accuracy and effectiveness of seismic elastic parameter inversion.

[0005] In order to solve the above technical problems, the technical solution created by the present invention is implemented as follows:

[0006] A prestack multi-wave joint inversion method constrained by ultra-differential Laplace blocks comprises the following steps:

[0007] 101. Extract statistical wavelets based on seismic data; Calculate the distribution law satisfied by the spatial derivatives of the longitudinal wave velocity, shear wave velocity and density of the work area based on the well logging data to determine the selection of parameters in the differential Laplace distribution; According to Bayesian theory, the prior information of the model characterizes certain characteristics of the variables to be inverted. Therefore, it is first necessary to perform statistical analysis on the variables to be inverted and obtain appropriate parameters to ensure that the distribution function under the parameters can contain all or most of the parameters to be inverted. When inverting actual data, the process often requires statistics on the well logging data, and the distribution law that satisfies the geological conditions of the area is obtained through the well data;

[0008] 102. Combine well logging data, seismic horizons and seismic structural interpretation data to establish an initial acoustic impedance model that conforms to the sedimentary model and seismic phase. Specifically, use seismic structural interpretation data to establish a geological model based on the sedimentary model, and interpolate and extrapolate the well logging data according to the structural model to obtain an initial smooth P- and S-wave velocity and density model;

[0009] 103. Calculate the formation reflection coefficient based on the analytical solution of the one-dimensional wave equation, calculate the initial synthetic PP wave and PS seismic data by convolution wavelet, and compare with the real seismic data to calculate the residual, so as to obtain the derivative of the forward operator to the model parameter; for the model data, the wavelet is artificially given, so it is sufficient to directly use the wavelet and the above-calculated reflection coefficient for convolution. For the actual seismic data, as described in step 101, it is necessary to count the wavelets from the original data, and it is usually necessary to count the data at different angles separately to obtain more accurate seismic wavelets, so as to ensure the accuracy of the inversion results;

[0010] 104. Based on the Bayesian principle, the ultra-differential Laplace distribution block constraint is added to construct the inversion objective function in the sense of maximum a posteriori probability to obtain the inversion result; the solution of the inverse problem is usually ill-posed, which causes the quality of the inversion result to deteriorate. Therefore, adding some prior information of the parameters to be inverted is conducive to updating the solution in the direction of the prior information. For the pre-stack three-parameter inversion, adding the ultra-differential Laplace block constraint can improve the resolution of the inversion result and make the characterization of the formation boundary clearer;

[0011] 105. According to the iteration of Gauss-Newton method, the final elastic parameter inversion result is obtained to guide the formation interpretation and reservoir prediction.

[0012] Further, in step 102, a spatial interpolation method is used to establish a model.

[0013] Further, the method of establishing the model in step 102 includes:

[0014] First, the scattered point interpolation method is used to interpolate the data of each layer to complete the geological layer modeling; then the P-wave and S-wave velocities and densities are transversely interpolated according to the geological layers to calculate the P-wave and S-wave velocities and density values ​​of each point underground to complete the initial elastic parameter modeling.

[0015] Further, in step 103, it is assumed that the underground medium is a horizontal layered medium composed of N isotropic layers, a point source and several detectors are arranged on the top interface of the first layer of the medium, and the total reflectivity It can be obtained by a 6-dimensional vector Recursively find

[0016] , (1)

[0017] Each vector represents the total response of all media below the nth layer.

[0018] , (2)

[0019] in is the transfer matrix of the nth layer medium, Describes the effect of the nth layer of medium on the seismic wave phase, and The effect of the nth interface on the amplitude of the downgoing and upgoing waves is described respectively;

[0020] Assume that there is only transmitted wave but no reflected wave in the last layer of medium, so the reflectivity of all waveforms is 0, that is,

[0021] , (3)

[0022] Recursively using formula (2) from layer N to layer 1, we can get the vector of layer 0:

[0023] , (4)

[0024] Then the total reflectivity of PP and PS waves of all strata can be obtained.

[0025] , (5)

[0026] Calculate the frequency and slowness range Finally, the seismic trace gathers in the time-space domain can be obtained by using the slowness integral and the frequency integral.

[0027] , (6)

[0028] in, Fourier transform of seismic wavelet, is the zero-order Bessel function, then p is sampled, and then the total reflectivity Perform frequency integration

[0029] , (7)

[0030] That is, the intercept time-ray parameter domain seismic record is obtained.

[0031] Furthermore, for equation (6), the slowness p needs to be fully sampled to avoid spatial aliasing effects.

[0032] Furthermore, in order to avoid oversampling of p, which would lead to a significant decrease in computational efficiency, it is allowed to directly calculate the total reflectivity Perform frequency integration.

[0033] Further, the step of obtaining the inversion objective function in step 104 includes:

[0034] First, the Bayesian method is applied to add the ultra-differential Laplace distribution block constraint that can better characterize the boundary information. The mathematical expression used is:

[0035] , (8)

[0036] in, is the posterior probability distribution; It is the likelihood function that maps from the model parameter space to the observed data space; is the prior distribution; is the marginal probability density. After omitting the marginal probability density term, the Bayesian equation can be expressed as follows

[0037] , (9)

[0038] For the ultradifferential Laplace block constraint, the prior distribution is

[0039] , (10)

[0040] in, It can be expressed as

[0041] , (11)

[0042] In the above formula, μ is the mean of different model parameters, usually set as the initial model. is the scaling factor of different model parameters, D is the first-order difference operator, which can be expressed as

[0043] , (12)

[0044] It can be expressed as

[0045] , (13)

[0046] The likelihood function can be expressed as

[0047] , (14)

[0048] According to formula (9), the posterior probability distribution based on the ultra-differential Laplace block constraint can be obtained:

[0049] , (15)

[0050] Therefore, the inversion objective function can be established as

[0051] , (16)

[0052] Finally, the objective function of prestack joint inversion based on ultradifferential Laplace block constraints can be expressed as

[0053] (17)

[0054] Furthermore, the marginal probability density is a constant whose definition rule is to make the sum of the posterior probability distribution equal to 1.

[0055] Furthermore, in the inversion, in addition to adding the ultra-differential Laplace distribution to characterize the boundary information, it is also necessary to add the Gaussian distribution to avoid severe jitter in the inversion result.

[0056] Furthermore, in the inversion, the objective function should include not only the data residual of the PP wave but also the data residual of the PS wave.

[0057] The advantages and positive effects of the invention are:

[0058] The method of the present invention can fully consider the transmission loss, conversion wave and interlayer multiple waves based on the analytical solution of the one-dimensional wave equation, so as to obtain a more accurate inversion result. PS conversion wave information is added on the basis of the traditional PP wave inversion, so as to obtain a more stable and accurate shear wave velocity and density inversion result. And based on the Bayesian theory, by adding the ultra-differential Laplace block constraint, the derivative of the model parameter is constrained, so as to obtain a more accurate characterization of the boundary information of the formation and obtain a higher resolution inversion result. Compared with the traditional Cauchy constraint and Laplace distribution, the new constraint method theoretically has a better effect on improving the resolution, thereby providing an accurate data body for subsequent reservoir interpretation and physical property parameter inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 It is a flow chart of the nonlinear prestack multi-wave joint inversion method based on ultra-differential Laplace block constraints of the present invention;

[0060] Figure 2 The comparison of the probability density functions of the ultra-differential Laplace distribution used in the embodiment of the present invention and the traditional Gaussian, Cauchy and Laplace distributions (2a), and the comparison of the probability density functions of the ultra-differential Laplace distribution with different parameters p (2b);

[0061] Figure 3 are noise-free pre-stack PP wave (3a) and PS wave (3b) seismic data input in the embodiment of the method of the present invention;

[0062] Figure 4 The embodiment of the present invention is based on Figure 3 The inversion results obtained from the noise-free prestack data, where the dashed line is the input initial model, the solid line is the true model, and the dotted line is the inversion result;

[0063] Figure 5 is a pre-stack seismic record with a signal-to-noise ratio of 10 dB input in an embodiment of the present invention;

[0064] Figure 6 The embodiment of the present invention is based on Figure 5 The inversion results obtained from the noise-free prestack data, where the dashed line is the input initial model, the solid line is the true model, and the dotted line is the inversion result. DETAILED DESCRIPTION

[0065] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0066] In the description of the present invention, it should be understood that the terms "center", "longitudinal", "lateral", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside" and the like indicate positions or positional relationships based on the positions or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation on the present invention. In addition, the terms "first", "second", etc. are only used for descriptive purposes, and cannot be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, features defined as "first", "second", etc. may explicitly or implicitly include one or more of the features. In the description of the present invention, unless otherwise specified, "multiple" means two or more.

[0067] In the description of the invention, it should be noted that, unless otherwise clearly specified and limited, the terms "installation", "connection" and "connection" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be indirectly connected through an intermediate medium, or it can be the internal communication of two components. For ordinary technicians in this field, the specific meanings of the above terms in the invention can be understood according to specific circumstances.

[0068] In view of the instability and incompleteness of the forward operator in the existing seismic pre-stack inversion method, the specific embodiments of the present invention are described in detail below. Figures 1 to 6 As shown, the present invention provides a pre-stack multi-wave joint inversion method constrained by an ultra-differential Laplace block:

[0069] The present invention assumes that the seismic wavelet is known before inversion, and therefore, it is necessary to extract the wavelet by using a statistical method based on the actual seismic pre-stack gathers and well logging data. In addition, the added ultra-differential Laplace block distribution needs to contain the model parameters to be inverted. Therefore, before inversion, it is necessary to perform a statistical analysis on the P- and S-wave velocities and density information through well logging data to determine appropriate parameters so that the obtained ultra-differential Laplace distribution satisfies the distribution law obeyed by the formation elastic parameters of the work area.

[0070] Combine logging data, seismic horizons and seismic structural interpretation data to establish an initial model that conforms to the sedimentary model and seismic phase: The elastic parameter model is mainly established by using the spatial interpolation method. The technical process is to first use the scattered point interpolation method to interpolate the data of each horizon to complete the geological horizon modeling, and then perform the longitudinal and transverse wave velocity and density lateral interpolation according to the geological horizon, that is, to perform the lateral interpolation of the logging information, calculate the longitudinal and transverse wave velocity and density values ​​of each point underground, and complete the task of initial elastic parameter modeling. For model data, low-frequency models (such as Figure 4 and Figure 6 The dashed line in the figure) is the true elastic parameter ( Figure 4 and Figure 6 The low-pass filtering result is shown in the solid line in .

[0071] The formation reflection coefficient is calculated based on the analytical solution of the one-dimensional wave equation, and the initial synthetic PP wave and PS seismic data are calculated by convolution wavelet. Assuming that the underground medium is a horizontal layered medium composed of N isotropic layers, a point source and several detectors are arranged on the top interface of the first layer of medium. Total reflectivity (i.e. the solution of the fully elastic wave equation in the frequency-slowness domain) can be obtained by a 6-dimensional vector Recursively find

[0072] , (1)

[0073] Each vector Represents the overall response of all media below the nth layer

[0074] , (2)

[0075] in is the transfer matrix of the nth layer medium, Describes the effect of the nth layer of medium on the seismic wave phase, and The influence of the nth interface on the amplitude of the downgoing wave and the upgoing wave is described respectively.

[0076] First, assume that there is only transmitted wave but no reflected wave in the last layer of medium, so the reflectivity of all waveforms is 0, that is,

[0077] , (3)

[0078] Recursively using formula (2) from layer N to layer 1, we can get the vector of layer 0:

[0079] , (4)

[0080] According to the above definition of vector, the total reflectivity of PP and PS waves of all strata can be obtained.

[0081] , (5)

[0082] Calculate the frequency and slowness range Finally, the seismic trace gathers in the time-space domain can be obtained by using the slowness integral and the frequency integral.

[0083] , (6)

[0084] in, Fourier transform of seismic wavelet, is the zero-order Bessel function. However, the above formula requires sufficient sampling of the slowness p, otherwise spatial aliasing will occur. Then oversampling p will lead to a significant decrease in computational efficiency. In order to save computational costs, the total reflectivity is directly Perform frequency integration

[0085] , (7)

[0086] That is, the intercept time-ray parameter domain seismic record.

[0087] (4) Based on the Bayesian principle, we add the ultra-differential Laplace distribution block constraint that can better describe the boundary information, and construct the inversion objective function in the sense of maximum a posteriori probability. The Bayesian method is a probabilistic statistical method used to calculate conditional probability. Mathematically, its expression is:

[0088] , (8)

[0089] in, is the posterior probability distribution; It is the likelihood function that maps from the model parameter space to the observed data space; It is the prior knowledge of the model parameters before observation, called the prior distribution; is a constant, called the marginal probability density, and its definition rule is to make the sum of the posterior probability distribution equal to 1. Generally, the marginal probability density term is omitted and the Bayesian equation is expressed as follows

[0090] , (9)

[0091] For the ultradifferential Laplace block constraint, the prior distribution is

[0092] , (10)

[0093] in, It can be expressed as

[0094] , (11)

[0095] In the above formula, μ is the mean of different model parameters, usually set as the initial model. is the scaling factor of different model parameters. D is the first-order difference operator, which can be expressed as

[0096] , (12)

[0097] It can be expressed as

[0098] , (13)

[0099] The likelihood function can be expressed as

[0100] , (14)

[0101] According to formula (9), the posterior probability distribution based on the ultra-differential Laplace block constraint can be obtained:

[0102] , (15)

[0103] Therefore, the inversion objective function can be established as

[0104] , (16)

[0105] Usually, in inversion, in addition to adding the ultra-differential Laplace distribution to characterize the boundary information, it is also necessary to add the Gaussian distribution to make the inversion result less jittery. In addition, for joint inversion, in addition to the PP wave data residual, the objective function should also include the PS wave data residual term. Finally, the objective function of prestack joint inversion based on ultra-differential Laplace block constraints can be expressed as

[0106] (17)

[0107] According to the Gauss-Newton method, we can obtain the final inversion results of P-wave and S-wave velocity and density. The iterative formula of the Gauss-Newton method is:

[0108] , (18)

[0109] in, is the first-order derivative (gradient) of the objective function with respect to the model parameters, is the second-order derivative of the objective function with respect to the model parameters, namely the Hessian matrix, and the two can be expressed as

[0110] , (19) , (20)

[0111] Among them, the second-order derivative terms of the forward operator with respect to the model parameters are ignored in the Hessian matrix.

[0112] After setting the number of iterations, the final elastic parameter inversion results can be obtained (such as Figure 4 and Figure 6 ). Figure 4 The embodiment of the present invention is based on Figure 3 The elastic parameter inversion results obtained by inverting the noise-free seismic data, where the dotted line is the input initial model, the solid line is the true model, and the dotted line is the inversion result. Figure 6 Based on Figure 5 The inversion results of noisy seismic data (signal-to-noise ratio is 10), where the dashed line is the input initial model, the solid line is the true model, and the dotted line is the inversion result.

[0113] The present invention is proposed on the basis of studying the following problems: (1) The forward modeling operator cannot correctly characterize the real propagation process of seismic waves underground, which leads to many false images in the inversion results; (2) The PP wave is severely attenuated when passing through a fluid-containing reservoir, especially when there is a gas cloud area in the reservoir, so the imaging effect of the reservoir area on the seismic profile is poor, affecting the inversion result. In addition, the PP wave is less sensitive to shear wave velocity and density, resulting in low inversion accuracy of the two; (3) The inversion based on the conventional regularization term cannot well characterize the blocky nature of the formation, resulting in a decrease in the resolution of the inversion and difficulty in characterizing the reservoir boundary. The present invention first introduces the analytical solution of the one-dimensional wave equation and points out the advantages of the forward modeling method compared with the conventional forward modeling method based on the Zoeppritz equation and its approximation. Then it is pointed out that the imaging effect of the converted wave data on the reservoir (especially the gas cloud area) is better than that of the longitudinal wave data, and then the PP and PS wave joint inversion method based on the analytical solution of the one-dimensional wave equation is introduced. Finally, in order to accurately characterize the reservoir boundary and improve the resolution and stability of the inversion results, the ultra-differential Laplace block constraint is introduced. The advantage of this constraint over the traditional Gaussian, Cauchy and Laplace distributions is demonstrated from the distribution image. Combined with Bayesian theory, the final high-resolution inversion result is obtained through the Gauss-Newton method, providing a data basis for formation interpretation and reservoir prediction.

[0114] In summary, the present invention has the following advantages:

[0115] 1. This technical solution can provide a high-resolution pre-stack elastic parameter inversion method. The traditional inversion method based on the Zoeppritz equation and its approximation cannot take into account the complex propagation law of seismic waves in the real strata, while this technical solution can fully consider the transmission loss, conversion wave and interlayer multiple waves based on the analytical solution of the one-dimensional wave equation, thereby obtaining a more accurate inversion result.

[0116] 2. This technical solution adds PS conversion wave information to the traditional PP wave inversion, thereby obtaining more stable and accurate shear wave velocity and density inversion results.

[0117] 3. Based on Bayesian theory, this technical solution constrains the derivatives of model parameters by adding ultra-differential Laplace block constraints, thereby obtaining a more accurate characterization of the boundary information of the formation and obtaining higher-resolution inversion results. Compared with traditional Cauchy constraints and Laplace distribution, the new constraint method theoretically has a better effect on improving resolution, thereby providing accurate data for subsequent reservoir interpretation and physical property parameter inversion.

[0118] It is obvious to those skilled in the art that the present invention is not limited to the details of the above exemplary embodiments, and that the present invention can be implemented in other specific forms without departing from the spirit or essential features of the present invention.

[0119] Therefore, no matter from which point of view, the embodiments should be regarded as exemplary and non-restrictive, and the scope of the invention is limited by the appended claims rather than the above description, so it is intended to include all changes within the meaning and scope of the equivalent elements of the claims in the invention. Any reference numerals in the claims should not be regarded as limiting the claims involved.

[0120] In addition, it should be understood that although the present specification is described according to implementation modes, not every implementation mode contains only one independent technical solution. This description of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment may also be appropriately combined to form other implementation modes that can be understood by those skilled in the art.

Claims

1. A pre-stack multi-wave joint inversion method constrained by ultra-differential Laplace blocks, characterized in that: The steps include:

101. Extract statistical wavelets based on seismic data; determine the distribution law satisfied by the spatial derivatives of the longitudinal wave velocity, shear wave velocity and density in the statistical work area based on the well logging data to determine the selection of parameters in the differential Laplace distribution; 102. Combine well logging data, seismic horizons and seismic structural interpretation data to establish initial P- and S-wave velocity and density models that conform to the sedimentary model and seismic phase; 103. Calculate the formation reflection coefficient based on the analytical solution of the one-dimensional wave equation, calculate the initial synthetic PP wave and PS seismic data through the convolution wavelet, and compare with the real seismic data to calculate the residual to obtain the derivative of the forward operator with respect to the model parameters; 104. Based on the Bayesian principle, the ultra-differential Laplace distribution block constraint is added to construct the inversion objective function in the sense of maximum a posteriori probability to obtain the inversion result; 105. According to the iteration of Gauss-Newton method, the final elastic parameter inversion result is obtained to guide the formation interpretation and reservoir prediction.

2. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 1, characterized in that: In step 102, a model is established using a spatial interpolation method.

3. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 1 or 2, characterized in that: The method of establishing the model in step 102 includes: First, the scattered point interpolation method is used to interpolate the data of each layer to complete the geological layer modeling; then the P-wave and S-wave velocities and densities are transversely interpolated according to the geological layers to calculate the P-wave and S-wave velocities and density values ​​of each point underground to complete the initial elastic parameter modeling.

4. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 1, characterized in that: In step 103, it is assumed that the underground medium is a horizontal layered medium composed of N isotropic layers. A point source and several detectors are arranged on the top interface of the first layer of the medium. The total reflectivity It can be obtained by a 6-dimensional vector Recursively find ，(1) Each vector Represents the overall response of all media below the nth layer ，(2) in is the transfer matrix of the nth layer medium, Describes the effect of the nth layer of medium on the seismic wave phase, and The effect of the nth interface on the amplitude of the downgoing and upgoing waves is described respectively; Assume that there is only transmitted wave but no reflected wave in the last layer of medium, so the reflectivity of all waveforms is 0, that is, ,(3) Recursively using formula (2) from layer N to layer 1, we can get the vector of layer 0: ,(4) Then the total reflectivity of PP and PS waves of all strata can be obtained. ，(5) Calculate the frequency and slowness range Finally, the seismic trace gathers in the time-space domain can be obtained by using the slowness integral and the frequency integral. ，(6) in, Fourier transform of seismic wavelet, is the zero-order Bessel function, then p is sampled, and then the total reflectivity Perform frequency integration ，(7) That is, the intercept time-ray parameter domain seismic record is obtained.

5. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 4, characterized in that: For equation (6), the slowness p needs to be fully sampled to avoid spatial aliasing effects.

6. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 4, characterized in that: In order to avoid oversampling of p, which will significantly reduce the computational efficiency, it is allowed to directly calculate the total reflectivity Perform frequency integration.

7. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 1, characterized in that: The steps of obtaining the inversion objective function in step 104 include: First, the Bayesian method is applied to add the ultra-differential Laplace distribution block constraint that can better characterize the boundary information. The mathematical expression used is: ，(8) in, is the posterior probability distribution; It is the likelihood function that maps from the model parameter space to the observed data space; is the prior distribution; is the marginal probability density. After omitting the marginal probability density term, the Bayesian equation can be expressed as follows ，(9) For the ultradifferential Laplace block constraint, the prior distribution is ，(10) in, It can be expressed as ，(11) In the above formula, μ is the mean of different model parameters, usually set as the initial model. is the scaling factor of different model parameters, j=1,2,3, D is the first-order difference operator, which can be expressed as ，(12) It can be expressed as ，(13) The likelihood function can be expressed as ，(14) According to formula (9), the posterior probability distribution based on the ultra-differential Laplace block constraint can be obtained: ，(15) Therefore, the inversion objective function can be established as ，(16) Finally, the objective function of prestack joint inversion based on ultradifferential Laplace block constraints can be expressed as (17)。 8. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 7, characterized in that: The marginal probability density is a constant whose definition rule makes the sum of the posterior probability distribution equal to 1.

9. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 7, characterized in that: In the inversion, in addition to adding the ultra-differential Laplace distribution to characterize the boundary information, it is also necessary to add the Gaussian distribution to avoid severe jitter in the inversion result.

10. The method for prestack multi-wave joint inversion constrained by ultra-differential Laplace block according to claim 7, characterized in that: In the inversion, the objective function should include not only the PP wave data residual but also the PS wave data residual.

Images