A seismic inversion method jointly constrained by a physical model and prior information

Abstract

The present invention relates to the technical field of seismic exploration, and specifically discloses a seismic inversion method jointly constrained by a physical model and prior information, which includes extracting a seismic wavelet from seismic data to determine a wavelet amplitude scaling factor; statistically analyzing the prior information of impedance parameters; establishing an initial impedance parameter model by using seismic structure interpretation data and logging data; obtaining a simplified approximate equation based on the assumption of weak elastic differences at interfaces, and using this simplified equation to forward model a seismic trace gather and calculate an inversion residual; rewriting the objective function as a function of impedance parameters by using the idea of generalized linear inversion, solving the impedance parameters by using an iterative reweighted least squares algorithm, and updating the impedance parameters; repeating the above steps until the inversion residual meets the requirements or reaches the maximum number of iterations, and outputting the final processing result. The advantage of the present invention is that it can achieve stable and high-precision prediction of impedance parameters.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration, and particularly to a seismic inversion method jointly constrained by a physical model and prior information. Background Art

[0002] Seismic exploration is a method of inferring the underground geological structure by artificially exciting seismic waves and recording their propagation characteristics underground. The propagation velocity, reflection, and transmission characteristics of seismic waves in different media carry rich geological information, and seismic impedance is one of the key parameters describing this information. Seismic impedance is defined as the product of the density of a rock and its P-wave velocity, which directly reflects the physical properties of the rock, such as porosity, saturation, rock type, etc. These properties are crucial for identifying oil and gas reservoirs and other economic minerals. With the continuous growth of global energy demand, the search for and development of new oil and gas resources has become increasingly urgent. The depletion of traditional oil and gas fields and the scarcity of new resources require exploration technologies to be more efficient and accurate. Seismic impedance inversion can provide high-resolution underground images, helping geologists accurately identify the boundaries, fluid properties, reservoir quality, and oil and gas content of oil and gas reservoirs, thereby increasing the exploration success rate and reducing drilling risks and costs. The core objective of wave impedance inversion is to convert seismic data into parameters that can directly reflect the physical properties of rocks, which is an indispensable step in reservoir prediction and reservoir description. As a comprehensive parameter, wave impedance can bridge the gap between the physical propagation characteristics of seismic waves and the actual geological properties of rocks. Through wave impedance inversion, geologists can indirectly estimate reservoir parameters such as porosity and permeability of formations, which is crucial for evaluating the productivity of reservoirs and formulating development strategies. With the in-depth research, various wave impedance inversion methods have emerged, including but not limited to trace integration inversion, generalized linear inversion, iterative inversion, non-linear inversion, etc. Each method has its advantages and limitations. For example, trace integration inversion is simple and direct but has limited accuracy, while generalized linear inversion has higher accuracy but is vulnerable to high-frequency noise. Model-based inversion technology breaks through the limitations of traditional seismic resolution and can theoretically obtain the same resolution as well logging data. However, the well logging information and high- and low-frequency components provided by the model in the inversion results lead to the non-uniqueness of inversion, and it is restricted by the number of wells and the distribution of well patterns. To reduce the non-uniqueness of inversion, some prior information is usually introduced during the inversion process. However, the existing methods for imposing constraints on prior information (such as the initial model) directly calculate the error with the inversion results, resulting in a decrease in the accuracy of the inversion results.

[0003] In summary, the current research on seismic impedance inversion methods has the following problems: 1. In the seismic impedance inversion by trace integration, due to the influence of the band-limited seismic wavelet, the resolution of the inversion result is low, and it is greatly affected by the initial impedance value, with a large cumulative error; 2. For the method of first inverting the reflection coefficient and then recursively calculating the impedance, it is greatly affected by seismic data noise and has a large cumulative error; 3. In order to reduce the non-uniqueness of the inversion result in the existing model-based seismic impedance inversion methods, prior constraint information such as the initial model is introduced, but directly using this as a constraint reduces the resolution of the inversion result; 4. The deterministic seismic impedance inversion method can only give a unique inversion result and cannot evaluate the reliability of the result, increasing the risk of subsequent interpretation. Summary of the Invention

[0004] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a seismic inversion method jointly constrained by a physical model and prior information. By adopting a model-based inversion strategy, a physical model is introduced to ensure that the prediction result satisfies the observed seismic data. The prior information constraint term of the initial model is designed to calculate the error between the low frequency of the inversion result and the prior of the initial model, reducing the non-uniqueness of the inversion result while improving the resolution of the inversion result. Based on the Bayesian framework, while predicting the inversion result, the uncertainty of the inversion result is given, meeting the requirements of high-precision seismic exploration and fine reservoir characterization.

[0005] The purpose of the present invention is realized through the following technical solutions: A seismic inversion method jointly constrained by a physical model and prior information, comprising the following steps:

[0006] Step 110: Extract the wavelet using the actual seismic data, forward simulate the seismic trace gather based on the well logging data and the post-stack seismic reflection coefficient equation, and determine the amplitude scaling factor in combination with the actual seismic data beside the well.

[0007] Step 120: Extract the impedance parameters and their means from all the well logging data in the work area, and statistically calculate their variances and vertical variogram matrices.

[0008] Step 130: Use the seismic data horizon interpretation data and well logging data to establish an initial impedance parameter model in the time domain.

[0009] Step 140: Starting from the accurate post-stack reflection coefficient equation, based on the assumption of weak elastic differences at the interface, a linear approximation equation is derived. Based on the initial impedance parameter model in the time domain and the simplified equation, forward simulate the post-stack seismic trace gather, and directly calculate the inversion residual from the forward simulation record and the actual record.

[0010] Step 150: Based on the Bayesian principle, introduce a sparse prior constraint term and an initial model prior information constraint term to construct an inversion objective function in the sense of maximum a posteriori probability, and use the generalized linear inversion idea to solve the inversion objective function to obtain the solution expression of the impedance parameters.

[0011] Step 160: Calculate the model parameters by using the iterative reweighted least squares algorithm according to the solution expression formula of the model parameters and the inversion residual.

[0012] Step 170: Repeat the above steps 140, 150, and 160 iteratively, and control the maximum number of iterations through the inversion residual to obtain the optimal inversion result of the impedance parameters.

[0013] Specifically, in step 110, using the logging data as the input model, the reflection coefficient is calculated by using the post-stack seismic reflection coefficient equation:

[0014] In the formula, m represents the logging impedance parameter data, and r represents the calculated reflection coefficient.

[0015] Then, the reflection coefficient is convolved with the extracted seismic wavelet to obtain a seismic trace gather, which is compared with the actual seismic trace gather beside the well, the amplitude scaling factor is calculated, and applied to the extracted seismic wavelet to achieve the amplitude matching between the simulated record and the actual record:

[0016]

[0017] In the formula, x represents the synthetic seismic record, represents the convolution operator, and w represents the extracted wavelet.

[0018] Specifically, in step 120, the required impedance parameters are obtained through logging data analysis, the autocorrelation coefficient of the impedance parameters is calculated, the variance matrix is constructed, and the vertical variogram is calculated based on the impedance parameters to form a prior distribution function of the impedance parameters that conforms to this work area.

[0019] The calculation formula of the vertical variogram is as follows:

[0020]

[0021] Among them, v represents the variogram value calculated based on well data, h represents the lag distance, N represents the length of the logging data, z(x i ) represents the impedance parameter value at a certain position i.

[0022] Specifically, in step 130, using the seismic structure interpretation data, a geological model is established based on the sedimentation model, and the logging data is interpolated and extrapolated according to the structural model to obtain the initial impedance parameter model of each survey line; the impedance parameter model is established by using the spatial interpolation method. First, the scattered point interpolation method is used to interpolate the data of each horizon to complete the geological horizon modeling, and then the impedance parameter is interpolated horizontally according to the geological horizon to calculate the impedance parameter value at each point underground to complete the task of initial impedance parameter modeling.

[0023] Specifically, in step 140, based on the weak elastic difference hypothesis, a linear approximation equation is derived:

[0024]

[0025] In the formula, m represents impedance parameter data, r represents the calculated reflection coefficient, and ln represents taking the logarithm of the data;

[0026] Then, using the initial impedance parameter model in the time domain as the input, the reflection coefficient vector is directly calculated using the simplified equation. The seismic wavelet is convolved with the reflection coefficient to obtain the post-stack seismic trace gather, and the difference is taken with the actual seismic trace gather to obtain the inversion residual:

[0027] In the formula, K represents the wavelet matrix constructed based on w, D represents the difference operator, and m represents impedance parameter data.

[0028] Specifically, in step 150, the inversion result is constrained using the initial model, and it is assumed that the seismic data noise follows a Gaussian distribution, the prior information constraint term of the initial model follows a Gaussian distribution, and the model follows a sparse distribution. Then, the inversion likelihood function and the prior probability distribution respectively satisfy a Gaussian distribution and a joint distribution of Gaussian and sparse. According to Bayes' principle, the posterior probability distribution function is obtained by synthesizing the inversion likelihood function and the prior distribution function. The objective function of the inversion is determined based on the posterior probability, and the derivative of the objective function with respect to the model parameters is taken to obtain the iterative solution formula;

[0029] Let the impedance model parameter be z T =(z 1 , z 2 ,..., z n ) T , and the observed seismic data be x T =(x 1 , x 2 ,..., x n ) T , according to Bayes' theory, in the case of known post-stack seismic data, the problem of inverting the impedance parameters of the underground medium can be reduced to solving a posterior probability function:

[0030]

[0031] where P(x)=∫p(x|z)p(z)dz is the normalization factor, which can be regarded as a constant, P(x|z) is the likelihood function, and P(z) is the prior probability distribution. Assume that the likelihood function satisfies a Gaussian distribution:

[0032]

[0033] In the above formula, Γ is the forward operator, CX is the covariance matrix of noise, N x is the length of the observed data. Based on Bayesian theory, combined with the prior model, after canceling out the normalization constant, the posterior probability distribution of post-stack seismic data inversion can be expressed as:

[0034]

[0035] In the formula, Γ is the forward operator, Φ(z) represents the sparse constraint term, and z 0 represents the initial impedance parameter, represents the smoothing operator, C X is the covariance matrix of noise, C z represents the model covariance matrix, which can be obtained by the variance α and the variogram matrix υ through the Kronecker product. The specific expressions of Φ(z) and are as follows:

[0036] Φ(z) = ||Dz|| 0 (6)

[0037] In the formula, D represents the difference operator. If three-point smoothing is adopted, the expression of can be written as:

[0038]

[0039] According to different specific smoothing degrees, the expression of can be modified according to Equation (10), and

[0040] The optimal solution can be obtained by solving the maximum value of Equation (8), which is equivalent to solving the solution corresponding to the minimum value of the following objective function,

[0041]

[0042] Taking the derivative of the above-obtained objective function with respect to the impedance parameter z and setting the derivative equal to zero. Since the 0-norm is not differentiable, the iteratively reweighted least squares algorithm is used for solution, and the updated iteration formula of the model parameter can be obtained:

[0043]

[0044] In the formula, L = 0.5·W·D represents the forward operator; y = ln(z) represents the logarithm of the impedance parameter; W represents the wavelet matrix, X represents the seismic record, I is the N x ×N x identity matrix, N x is the length of the observed data; diag represents the operator for constructing the diagonal matrix; λ 1 and λ2 respectively represent the weights of the initial model prior information constraint term and the sparse constraint term; μ is a constant to ensure diagonal dominance of the matrix. Specifically, in step 170, based on the Bayesian framework, while calculating the optimal impedance parameter inversion result, the uncertainty information of the inversion result is calculated. The seismic inversion method based on the joint constraint of the physical model and prior information can predict the inversion result and give the uncertainty of the inversion result at the same time, realizing the quantitative reliability evaluation of the inversion result;

[0045] ∑ = C z -((LC z ) T / (LC z L T +μI))·LC z (10)

[0046] In the formula, C z represents the model covariance matrix; L represents the forward operator; μ is a constant to ensure diagonal dominance of the matrix; I is the N x ×N x identity matrix, and N x is the length of the observed data.

[0047] A computer device includes a memory and a processor. The memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions. When the computer-executable instructions are executed by the processor, the steps of the above method are implemented.

[0048] A computer-readable storage medium stores computer-executable instructions. When the computer-executable instructions are executed by a processor, the steps of the above method are implemented.

[0049] The present invention has the following advantages:

[0050] 1. The present invention directly estimates impedance from seismic data, avoiding the cumulative error caused by the trace integration and recursive inversion methods;

[0051] 2. The present invention introduces a physical model to ensure that the prediction result meets the observed seismic data, and designs an initial model prior information constraint term to reduce the non-uniqueness of the inversion result and improve the accuracy of the inversion result while;

[0052] 3. The present invention introduces a sparse prior constraint term and a vertical constraint term, which can well highlight the boundary characteristics of the inversion result and ensure the vertical continuity of the inversion result;

[0053] 4. Based on the Bayesian framework, the uncertainty of the inversion result is given while predicting the inversion result, and the quantitative reliability evaluation of the inversion result can be realized. Description of the Drawings

[0054] Figure 1 Flow chart of seismic inversion method jointly constrained by physical model and prior information of the present invention;

[0055] Figure 2 Model diagram of noise-free post-stack seismic data and initial impedance parameter input for the present invention;

[0056] Figure 3 Model diagram of noisy post-stack seismic data and initial impedance parameter input for the present invention;

[0057] Figure 4 Impedance and its uncertainty information diagram obtained by seismic inversion method jointly constrained by physical model and prior information without noise for the present invention;

[0058] Figure 5 Impedance and its uncertainty information diagram obtained by seismic inversion method jointly constrained by physical model and prior information with noise for the present invention. Detailed implementation manners

[0059] The present invention will be further described below in conjunction with the accompanying drawings, but the protection scope of the present invention is not limited to the following. As Figures 1 to 5 shown, a seismic inversion method jointly constrained by physical model and prior information includes the following steps:

[0060] Step 110: Assume that the seismic wavelet is known before inversion. Therefore, it is necessary to extract the wavelet by statistical method based on the actual seismic trace gather and well logging data. The actual seismic amplitude is often a relative value, and there is a certain numerical difference between the amplitude of the seismic data forward-simulated by the post-stack seismic reflection coefficient equation and the actual amplitude. Using the well logging data as the input model, calculate the reflection coefficient by the post-stack seismic reflection coefficient equation:

[0061]

[0062] In the formula, m represents the well logging impedance parameter data, and r represents the calculated reflection coefficient;

[0063] Then convolve the reflection coefficient with the extracted seismic wavelet to obtain a seismic trace gather, compare it with the actual seismic trace gather beside the well, calculate the amplitude scaling factor, and apply it to the extracted seismic wavelet to achieve the amplitude matching between the simulated record and the actual record;

[0064]

[0065] In the formula, x represents the synthetic seismic record, represents the convolution operator, and w represents the extracted wavelet;

[0066] Step 120: Extract impedance parameters and their means from all well logging data in the work area, and statistically calculate their variances and vertical variogram matrices; obtain the required impedance parameters through well logging data analysis, calculate the autocorrelation coefficient of the impedance parameters, construct a variance matrix, and calculate the vertical variogram based on the impedance parameters to form a prior distribution function of the impedance parameters that conforms to this work area; the calculation formula of the vertical variogram is as follows:

[0067]

[0068] where υ represents the variogram value calculated based on well data, h represents the lag distance, N represents the length of well logging data, and z(x i ) represents the impedance parameter value at a certain position i

[0069] Step 130: Use seismic data horizon interpretation data and well logging data to establish an initial impedance parameter model in the time domain; use seismic structural interpretation data to establish a geological model based on the sedimentation pattern, and interpolate and extrapolate the well logging data according to the structural pattern to obtain the initial impedance parameter model of each survey line; establish an impedance parameter model using spatial interpolation methods. First, use the scattered point interpolation method to interpolate the data of each horizon to complete the geological horizon modeling, and then perform lateral interpolation of the impedance parameters according to the geological horizons to calculate the impedance parameter values at each point underground to complete the initial impedance parameter modeling;

[0070] Step 140: Starting from the accurate post-stack reflection coefficient equation, based on the assumption of weak elastic differences at the interface, derive a linear approximation equation:

[0071]

[0072] In the formula, m represents impedance parameter data, r represents the calculated reflection coefficient, and ln represents taking the logarithm of the data;

[0073] Based on the initial impedance parameter model in the time domain and the simplified equation, forward simulate the post-stack seismic gather. Using the initial impedance parameter model in the time domain as the input, directly calculate the reflection coefficient vector using the simplified equation, convolve the seismic wavelet with the reflection coefficient to obtain the post-stack seismic gather, and subtract it from the actual seismic gather to obtain the inversion residual:

[0074]

[0075] In the formula, K represents the wavelet matrix constructed based on w, D represents the difference operator, and m represents impedance parameter data.

[0076] Step 150: Based on the Bayesian principle, introduce a sparse prior constraint term and an initial model prior information constraint term to construct an inversion objective function in the sense of maximum a posteriori probability. Use the generalized linear inversion idea to solve the inversion objective function to obtain the solution expression of the impedance parameters. Constrain the inversion result using the initial model, and assume that the seismic data noise follows a Gaussian distribution, the initial model prior information constraint term follows a Gaussian distribution, and the model follows a sparse distribution. Then the inversion likelihood function and the prior probability distribution respectively satisfy a Gaussian distribution and a joint distribution of Gaussian and sparse. According to the Bayesian principle, synthesize the inversion likelihood function and the prior distribution function to obtain the posterior probability distribution function. Determine the inversion objective function according to the posterior probability, and take the derivative of the objective function with respect to the model parameters to obtain the iterative solution formula;

[0077] Let the impedance model parameter be z T =(z 1 ,z 2 ,L,z n ) T , and the observed seismic data be x T =(x 1 ,x 2 ,L,x n ) T . According to the Bayesian theory, it can be obtained that, given the post-stack seismic data, the problem of inverting the impedance parameters of the underground medium can be reduced to solving a posterior probability function:

[0078]

[0079] where P(x)=∫p(x|z)p(z)dz is the normalization factor, which can be regarded as a constant, P(x|z) is the likelihood function, and P(z) is the prior probability distribution. Assume that the likelihood function satisfies a Gaussian distribution:

[0080]

[0081] In the above formula, Γ is the forward operator, C X is the covariance matrix of the noise, N x is the length of the observed data. Based on the Bayesian theory, combined with the prior model, cancel out the normalization constant, then the post-stack seismic data inversion posterior probability distribution can be expressed as:

[0082]

[0083] In the formula, Γ is the forward operator, Φ(z) represents the sparse constraint term, z 0 represents the initial impedance parameter, represents the smoothing operator, C Xis the covariance matrix of noise, C z represents the model covariance matrix, which can be the variance statistically calculated from well logging data α and the variogram matrix υ obtained through the Kronecker product. The specific expressions of Φ(z) and are as follows:

[0084] Φ(z) = ||Dz|| 0 (14)

[0085] In the formula, D represents the difference operator. If three-point smoothing is adopted, the expression of

[0086]

[0087] According to different specific smoothing degrees, the expression of can be modified according to formula (10).

[0088] The optimal solution can be obtained by solving the maximum value of formula (8), which is equivalent to solving the solution corresponding to the minimum value of the following objective function:

[0089]

[0090] Taking the derivative of the above-obtained objective function with respect to the impedance parameter z and setting the derivative equal to zero. Since the 0-norm is not differentiable, the iteratively reweighted least squares algorithm is used for solving, and the updated iteration formula of the model parameters can be obtained:

[0091]

[0092] In the formula, L = 0.5·W·D represents the forward operator; y = ln(z) represents the logarithm of the impedance parameter; W represents the wavelet matrix, X represents the seismic record, I is the N x ×N x identity matrix, N x is the length of the observed data; diag represents the operator for constructing a diagonal matrix; λ 1 and λ 2 represent the weights of the initial model prior information constraint term and the sparse constraint term respectively; μ is a constant to ensure diagonal dominance of the matrix.

[0093] Step 160: Calculate the model parameters using the iteratively reweighted least squares algorithm according to the solution expression formula of the model parameters and the inversion residual;

[0094] Step 170: Repeat the above steps 140, 150, and 160 iteratively, and control the maximum number of iterations through the inversion residual to obtain the optimal impedance parameter inversion result; such as Figure 4 andFigure 5 As shown, it is the impedance parameter obtained by the seismic inversion method based on the joint constraint of the physical model and prior information in the embodiment of the present invention. In the figure, the vertical axis represents time in seconds, and the horizontal axis represents impedance (unit: g / cm 3 ·km / s)). The present invention can predict the impedance parameter information with high accuracy, such as Figure 5 As shown is the inversion result when adding random noise with a signal-to-noise ratio of 4. The introduction of the prior model plays a key role in maintaining the stability of the inversion process and improving the accuracy of the inversion result. Based on the Bayesian framework, while calculating the inversion result, the uncertainty information of the inversion result can be calculated (Equation 13), and the uncertainty is as Figure 4 shown on the right, Figure 5 as shown by the color of the bottom right figure. The present invention can give the uncertainty of the inversion result while predicting the inversion result, and can realize the quantitative reliability evaluation of the inversion result.

[0095] ∑ = C z -((LC z ) T / (LC z L T + μI))·LC z (18)

[0096] In the formula, C z represents the model covariance matrix; L represents the forward operator; μ is a constant to ensure diagonal dominance of the matrix; I is the N x ×N x identity matrix, and N x is the length of the observed data.

[0097] A computer device includes a memory and a processor. The memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions. When the computer-executable instructions are executed by the processor, the steps of the above method are implemented.

[0098] A computer-readable storage medium stores computer-executable instructions. When the computer-executable instructions are executed by the processor, the steps of the above method are implemented.

[0099] Through the above method, the present invention can achieve stable and high-precision prediction of impedance parameters.

Claims

1. A seismic inversion method jointly constrained by a physical model and prior information, characterized in that: The following steps are involved: Step 110: extracting wavelets using actual seismic data, forward modeling seismic gathers based on well logging data and post-stack seismic reflection coefficient equations, and determining amplitude scaling factors in combination with actual wellside seismic data; Step 120: extracting impedance parameters and their mean values ​​based on all logging data in the work area, and statistically calculating their variances and vertical variogram matrices; Step 130: using the seismic data layer interpretation information and the logging data, establish an initial impedance parameter model in the time domain; Step 140: using the accurate post-stack reflection coefficient equation and based on the assumption of weak elastic difference of the interface, a linear approximate equation is derived, and the post-stack seismic gathers are forward modeled based on the initial impedance parameter model in the time domain and the simplified equation, and the inversion residual is directly calculated from the forward simulation record and the actual record; Step 150: Based on the Bayesian principle, a sparse prior constraint term and an initial model prior information constraint term are introduced to construct an inversion objective function in the sense of maximum a posteriori probability, and the inversion objective function is solved using the generalized linear inversion idea to obtain a solution expression for the impedance parameter; Step 160: Calculate the model parameters using an iterative reweighted least squares algorithm based on the solution expression formula of the model parameters and the inversion residual; Step 170: Repeat the iteration steps 140, 150 and 160, and control the maximum number of iterations by the inversion residual, and output the optimal impedance parameter inversion result.

2. The seismic inversion method with joint constraints of physical model and prior information according to claim 1, characterized in that: In step 170, based on the Bayesian framework, while calculating the optimal impedance parameter inversion result, the uncertainty information of the inversion result is calculated to achieve quantitative reliability evaluation of the inversion result: Σ=C z -((LC z ) T / (LC z L T +μI))·LC z (13) In the formula, C z represents the model covariance matrix; L represents the forward operator; μ is a constant to ensure that the matrix is ​​diagonally dominant; I is N x ×N x The unit matrix, N x is the length of the observation data.

3. The seismic inversion method with joint constraints of physical model and prior information according to claim 1, characterized in that: In step 110, the reflection coefficient is calculated using the post-stack seismic reflection coefficient equation with the logging data as the input model: In the formula, m represents the logging impedance parameter data, and r represents the calculated reflection coefficient; Then, the reflection coefficient is convolved with the extracted seismic wavelet to obtain a seismic gather, which is compared with the actual well-side seismic gather, and the amplitude scaling factor is calculated and applied to the extracted seismic wavelet to achieve amplitude matching between the simulated record and the actual record: Where x represents the synthetic seismic record, represents the convolution operator, and w represents the extracted wavelet.

4. A seismic inversion method with joint constraints of physical model and prior information according to claim 3, characterized in that: In step 120, the required impedance parameters are obtained through logging data analysis, the autocorrelation coefficient of the impedance parameters is obtained, a variance matrix is ​​constructed, and a vertical variogram matrix is ​​calculated based on the impedance parameters to form an impedance parameter prior distribution function that meets the work area; The calculation formula of the vertical variation function is as follows: Among them, υ represents the variogram value calculated based on well data, h represents the lag distance, N represents the length of logging data, z(x i ) represents the impedance parameter value at a certain position i.

5. A seismic inversion method with joint constraints of physical model and prior information according to claim 4, characterized in that: In step 130, the seismic structural interpretation data is used to establish a geological model based on the sedimentation model, and the logging data is interpolated and extrapolated according to the structural model to obtain an initial impedance parameter model for each survey line.

6. A seismic inversion method with joint constraints of physical model and prior information according to claim 5, characterized in that: In step 130, the impedance parameter model is established using a spatial interpolation method, which first uses a scattered point interpolation method to interpolate the data of each layer to complete the geological layer modeling, and then performs lateral interpolation of the impedance parameters according to the geological layer to calculate the impedance parameter value at each point underground to complete the initial impedance parameter modeling.

7. A seismic inversion method with joint constraints of physical model and prior information according to claim 6, characterized in that: In step 140, based on the weak elastic difference assumption, a linear approximate equation is derived: In the formula, m represents the impedance parameter data, r represents the calculated reflection coefficient, and ln represents the logarithm of the data; Then, the initial impedance parameter model in the time domain is used as input, and the reflection coefficient vector is directly calculated using the simplified equation. The seismic wavelet is convolved with the reflection coefficient to obtain the post-stack seismic gathers, which are then subtracted from the actual seismic gathers to obtain the inversion residual: Where K represents the wavelet matrix constructed based on w, D represents the differential operator, and m represents the impedance parameter data.

8. A seismic inversion method with joint constraints of physical model and prior information according to claim 7, characterized in that: In step 150, the initial model is used to constrain the inversion results, and it is assumed that the seismic data noise obeys the Gaussian distribution, the initial model prior information constraint item obeys the Gaussian distribution and the model obeys the sparse distribution, then the inversion likelihood function and the prior probability distribution satisfy the Gaussian distribution and the joint distribution of Gaussian and sparse respectively. According to the Bayesian principle, the inversion likelihood function and the prior distribution function are integrated to obtain the posterior probability distribution function, and the objective function of the inversion is determined according to the posterior probability distribution function. The objective function is derived from the model parameters to obtain the iterative solution formula: Where L = 0.5·W·D represents the forward operator; D represents the difference operator; W represents the wavelet matrix, X represents the seismic record, represents the smoothing operator, y=ln(z) represents the logarithm of the impedance parameter; I is N x ×N x The unit matrix, N x is the length of the observed data; diag represents the operator for constructing a diagonal matrix; λ1 and λ2 represent the weights of the initial model prior information constraint and the sparse constraint, respectively; μ is a constant that ensures that the matrix is ​​diagonally dominant.

9. A computer device, characterized in that: The invention comprises a memory and a processor, wherein the memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions. When the computer-executable instructions are executed by the processor, the steps of the method according to any one of claims 1 to 8 are implemented.

10. A computer-readable storage medium storing computer-executable instructions, characterized in that: When the computer executable instructions are executed by a processor, the steps of the method described in any one of claims 1 to 8 are implemented.

Images