Porosity and fluid factor seismic synchronous inversion method, device, medium and equipment

Abstract

The invention relates to a porosity and fluid factor seismic synchronous inversion method, device, medium and equipment, comprising the following steps: deriving an isotropic saturated medium precise Zeoppritz equation; forward modeling seismic angle gathers based on well logging data and the isotropic saturated medium precise Zeoppritz equation, and determining amplitude scaling factors in combination with actual wellside seismic data; statistically analyzing the mean of each model parameter based on all well logging data in a work area, obtaining the autocorrelation coefficient and the mutual correlation coefficient of each model parameter, constructing a parameter-related covariance matrix, and forming a model parameter prior distribution function that conforms to the work area; smoothing the well logging data and extrapolating along the interpretation layer to obtain an initial parameter model; constructing an inversion objective function in the sense of maximum a posteriori probability, and solving the inversion objective function to obtain updated and iterative porosity and fluid factor model parameter formulas.

Description

Technical Field

[0001] The invention relates to a porosity and fluid factor seismic synchronous inversion method, device, medium and equipment, belonging to the field of conventional or unconventional oil and gas reservoir seismic exploration and fluid prediction. Background Art

[0002] Seismic inversion is an effective method to obtain internal images of underground media and to describe reservoirs in detail. It is also the final manifestation of high-resolution seismic exploration. Seismic data inversion has greatly improved the value of reservoir characterization. As the focus of seismic research gradually shifts from exploration to development, and from conventional oil and gas to unconventional oil and gas, especially shale gas, the technology of revealing the fine distribution characteristics of underground oil and gas reservoirs through means such as seismic inversion and carving oil and gas reservoirs in detail has also received more and more attention. Under the guidance of rock physics theory, elastic parameters that are sensitive to abnormal characteristics of reservoir fluids can be characterized as longitudinal and transverse wave fluid parameters. From the current research, seismic prestack inversion research mainly involves AVO inversion methods based on approximate and exact Zoeppritz equations. The linear inversion method based on the approximate Zoeppritz equation can stably solve the ill-conditioned inversion problem or the local extreme value problem, but the various assumptions in the derivation process cause large errors in the reflection coefficient in the medium and large angle range, fail to make full use of the rich prestack information, and cannot meet the requirements of fine characterization of actual reservoirs in terms of inversion accuracy, especially for actual three-dimensional large offset seismic data. For isotropic saturated media, the traditional approximate reflection coefficient based on the Zoeppritz equation is an expression that includes elastic parameters and fluid characteristic terms. However, the elastic parameters or fluid characteristic terms obtained by conventional inversion are difficult to directly characterize the reservoir porosity distribution characteristics, and step-by-step inversion or indirect calculation will introduce cumulative errors. Summary of the invention

[0003] In view of the above technical problems, the present invention provides a porosity and fluid factor seismic synchronous inversion method, device, medium and equipment, which can meet the requirements of seismic pre-stack inversion for quantitative characterization and accurate prediction of fluid characteristics and porosity distribution laws of isotropic saturated media.

[0004] To achieve the above object, the present invention adopts the following technical solutions:

[0005] A porosity and fluid factor seismic synchronous inversion method comprises the following steps:

[0006] Derive the exact Zeoppritz equation for isotropic saturated media;

[0007] Extract seismic wavelets that depend on the incident angle based on seismic data of oil and gas reservoirs, simulate seismic angle gathers based on well logging data and accurate Zeoppritz equations for isotropic saturated media, and determine amplitude scaling factors in combination with actual wellside seismic data;

[0008] Based on all the logging data in the work area, the mean of each model parameter is calculated, the autocorrelation coefficient and the mutual correlation coefficient of each model parameter are obtained, and the covariance matrix related to the parameters is constructed to form a priori distribution function of the model parameters that conforms to the work area;

[0009] The logging data are smoothed and extrapolated along the interpreted horizon to obtain the initial parameter model;

[0010] Based on Bayesian theory, model parameter prior distribution function, amplitude scaling factor and initial parameter model, an inversion objective function in the sense of maximum a posteriori probability is constructed, and the inversion objective function is solved to obtain updated iterative porosity and fluid factor model parameter formulas.

[0011] The porosity and fluid factor seismic synchronous inversion method preferably extracts seismic wavelets that depend on the incident angle based on seismic data of oil and gas reservoirs, forward models seismic angle gathers based on well logging data and accurate Zeoppritz equations for isotropic saturated media, and determines amplitude scaling factors in combination with actual wellside seismic data, specifically including:

[0012] Assuming that the seismic wavelet is known before inversion, it is necessary to extract the wavelet using statistical methods based on the actual pre-stack seismic gathers and logging data. The wavelet will change in waveform or frequency due to the influence of the formation during propagation. Extracting the seismic wavelet that depends on the incident angle can effectively improve the amplitude matching degree.

[0013] Taking logging data as input, the model uses the accurate Zeoppritz equation of isotropic saturated medium to forward simulate the PP wave gathers in the angle domain. 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.

[0014] In the porosity and fluid factor seismic synchronous inversion method, preferably, the model parameter prior distribution function is:

[0015]

[0016] in, is the parameter model, ρ is the density, φ is the porosity, F is the improved fluid coefficient, is the longitudinal wave velocity of dry rock, Dry rock shear wave velocity; μ and C m are the mean of the model parameters and the covariance matrix of the model parameters respectively.

[0017] The porosity and fluid factor seismic synchronous inversion method preferably smoothes the logging data and extrapolates along the interpretation horizon to obtain an initial parameter model, as follows:

[0018] First, the scattered point interpolation method is used to interpolate the data of each layer to complete the geological layer modeling. Then, the model parameters are transversely interpolated according to the geological layer, that is, the logging information is transversely interpolated to calculate the model parameter value at each point underground, and the establishment of the initial new attribute parameter model is completed.

[0019] In the porosity and fluid factor seismic synchronous inversion method, preferably, the specific solution process of the inversion objective function is as follows:

[0020] By taking the derivative of the inversion objective function with respect to the model parameters and solving the minimum value of the inversion objective function, which is equivalent to setting its derivative equal to zero, the updated iterative porosity and fluid factor model parameter formulas are obtained by combining the iterative least squares algorithm.

[0021] The second aspect of the present invention provides a porosity and fluid factor seismic synchronous inversion device, comprising:

[0022] The first processing unit is used to derive the exact Zeoppritz equation for isotropic saturated media;

[0023] The second processing unit is used to extract the seismic wavelet that depends on the incident angle based on the seismic data of the oil and gas reservoir, to simulate the seismic angle gather based on the well logging data and the accurate Zeoppritz equation of isotropic saturated medium, and to determine the amplitude scaling factor in combination with the actual wellside seismic data;

[0024] The third processing unit is used to calculate the mean of each model parameter based on all the logging data in the work area, obtain the autocorrelation coefficient and the mutual correlation coefficient of each model parameter, construct the covariance matrix related to the parameters, and form a priori distribution function of the model parameters that conforms to the work area;

[0025] The fourth processing unit is used to perform smoothing processing on the logging data and extrapolate along the interpretation layer to obtain an initial parameter model;

[0026] The fifth processing unit is used to construct an inversion objective function in the sense of maximum a posteriori probability based on Bayesian theory, model parameter prior distribution function, amplitude scaling factor and initial parameter model, and solve the inversion objective function to obtain updated iterative porosity and fluid factor model parameter formulas.

[0027] A third aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-mentioned porosity and fluid factor seismic synchronous inversion method.

[0028] A fourth aspect of the present invention provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the above-mentioned porosity and fluid factor seismic synchronous inversion method when executing the computer program.

[0029] The present invention adopts the above technical solution, which has the following advantages:

[0030] 1. Compared with the linear AVO inversion of the approximate Zeoppritz equation, the seismic prestack saturated medium AVO precise Zeoppritz equation inversion adopted in the present invention has high inversion calculation accuracy based on the precise equation, strong purposefulness of inversion parameters, avoids the cumulative error caused by step-by-step inversion or indirect calculation, and has good industrial production application value.

[0031] 2. The precise Zeoppritz equation of multi-parameter reflection coefficient of saturated medium established in the present invention can directly estimate the fluid factor and porosity without reducing the inversion accuracy, and guide the prediction and evaluation of reservoir fluid characteristics, compared with the precise Zeoppritz equation of elastic medium reflection coefficient. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is a processing flow chart of a porosity and fluid factor AVO inversion method based on an accurate Zeoppritz equation provided in one embodiment of the present invention;

[0033] Figure 2 is a CDP stacked cross-sectional diagram provided by this embodiment of the present invention, wherein the vertical line in the diagram is the well logging position;

[0034] Figure 3 a-3b is the elastic parameter (P-wave velocity, S-wave velocity, density) inversion profile provided by this embodiment of the present invention, and the vertical curve is the well logging data, wherein Figure 3 a is the P-wave velocity inversion profile, Figure 3 b is the shear wave velocity inversion section, Figure 3 c is the density inversion profile;

[0035] Figure 4 a-4b is the inversion section of porosity and fluid factor, and the vertical curve is the logging data. Figure 4 a is the porosity inversion section, 4b is the fluid factor inversion section;

[0036] Figure 5 a-5c is the synthetic pseudo PS wave gather of the logging position inversion curve ( Figure 5 c) and the actual PP channel gather ( Figure 5 a) contrast, Figure 5b) is the flattened gather after KL transformation;

[0037] Figure 6 a-6b is a comparison chart of the inversion curve of the logging position. Figure 6 a is the inversion result of PP wave alone, Figure 6 b is the result of multi-wave joint inversion; the solid line is the inversion curve, the dotted solid line is the actual logging curve, and the dotted line is the initial smoothing model. DETAILED DESCRIPTION

[0038] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention is described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary persons in the art without creative work are within the scope of protection of the present invention.

[0039] The present invention aims at predicting reservoir fluid characteristics of isotropic saturated media. The traditional Zoeppritz equation-based approximate reflection coefficient has too many unknown parameters, resulting in inversion pathology, and the fluid parameters and porosity cannot be directly estimated. A porosity and fluid factor AVO inversion method based on the precise Zeoppritz equation is proposed to solve the problem. First, the saturated seismic wave velocity containing fluid parameters and porosity is introduced into the elastic medium Zoeppritz equation to characterize the isotropic saturated medium. Then, using perturbation theory and Taylor expansion method, the problem of solving the objective function is transformed into obtaining the gradient or partial derivative of the reflection coefficient with respect to the model parameters. Combining the iterative least squares algorithm and Bayesian theory, an iterative update expression of the model parameters is established, and a porosity and fluid factor seismic AVO synchronous inversion method based on the precise Zeoppritz equation is proposed, which finally realizes the quantitative characterization and accurate prediction of fluid characteristics and porosity distribution law.

[0040] The technical solution of the present invention is explained in detail below in conjunction with specific implementation methods.

[0041] like Figure 1 As shown, the porosity and fluid factor AVO inversion method based on the accurate Zeoppritz equation provided by the present invention has the following specific steps:

[0042] Step 1: Extract angle-dependent wavelets based on seismic data; forward model the seismic angle gathers based on well logging data and the accurate Zeoppritz equation for isotropic saturated media, and determine the amplitude scaling factor in combination with actual wellside seismic data.

[0043] The specific details are as follows:

[0044] The present invention assumes that the seismic wavelet is known before inversion, so it is necessary to extract the wavelet by using statistical methods based on the actual pre-stack seismic gathers and logging data. The wavelet will change in waveform or frequency due to the influence of the stratum during propagation. Extracting the seismic wavelet that depends on the incident angle can effectively improve the amplitude matching degree. Figure 2 As shown, from Figure 2 It can be seen that there is a strong reflection axis at the position of 2210-2255ms.

[0045] The actual seismic amplitude is often a relative value. There is a certain numerical difference between the amplitude of seismic data simulated by the accurate Zeoppritz equation of isotropic saturated media and the actual amplitude. The well logging data is used as the input model to simulate the PP wave gather in the angle domain using the accurate Zeoppritz equation of isotropic saturated media. 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. When the signal-to-noise ratio of seismic data is high, a unified amplitude scaling factor is used for each channel of the angle gather to ensure the relationship between the amplitude and the offset; when the signal-to-noise ratio is low, the amplitude scaling factor can be calculated separately for near, medium and far offsets to ensure the best match between the simulated record and the actual record and reduce the impact of noise on the inversion process.

[0046] Step 2: Calculate the prior information of the model parameters based on all the logging data in the work area, including the mean of the three AVO elastic parameters (P-wave velocity, S-wave velocity, density), porosity and fluid parameters, and the covariance matrix related to a total of five model parameters. The details are as follows:

[0047] The present invention adopts the Cauchy distribution function as the prior distribution function. Based on the statistical means of all the logging data in the work area, the autocorrelation coefficients and the mutual correlation coefficients of the parameters are obtained, and the covariance matrix related to the parameters is constructed to form the model parameter prior distribution function that conforms to the work area. In the subsequent inversion objective function, the corresponding regularized expression of the Cauchy distribution function is:

[0048]

[0049] in, is a parameter model, ρ is density, φ is porosity, F is the improved fluid coefficient, is the longitudinal wave velocity of dry rock, Dry rock shear wave velocity; μ and C m are the mean of the statistical model parameters and the covariance matrix of the model parameters. Since the Cauchy distribution integrates the correlation between the five parameters through the covariance matrix, the uncertainty between the model parameters is reduced.

[0050] Step 3: Smooth the logging data and extrapolate along the interpretation horizon to obtain the initial parameter model. That is, the initial parameter model is mainly established by using the three-dimensional space interpolation method. The technical process is to first interpolate the data of each horizon using the scattered point interpolation method to complete the geological horizon modeling, and then perform horizontal interpolation of model parameters according to the geological horizon, that is, horizontal interpolation of logging information, calculate the model parameter value at each point underground, and complete the task of initial new attribute parameter modeling.

[0051] Step 4: Derive the exact Zeoppritz equation for isotropic saturated media:

[0052] In order to introduce fluid characteristics, the longitudinal and transverse wave velocities (V P and V S ) is replaced by the saturated longitudinal and transverse wave velocities, so the exact Zeoppritz equation for the saturated medium can be expressed as:

[0053] S·R=b,

[0054] Where R represents the saturated medium reflection coefficient matrix, and denote the longitudinal wave reflection and transmission coefficients of saturated media, respectively. and represents the shear wave reflection and transmission coefficient of saturated medium; e 11 ~e 44 are the coefficient equations related to the incident angle and the P-wave, S-wave velocity and density of the upper and lower interfaces of the saturated medium;

[0055]

[0056] in,

[0057]

[0058] S 12 is the square value of the ratio of the shear wave velocity to the longitudinal wave velocity of the dry rock on the upper interface, S 13 is the ratio of the longitudinal wave velocity of the lower interface to the upper interface of the saturated medium, S 14 S is the square value of the ratio of the shear wave velocity of the lower interface to the longitudinal wave velocity of the upper interface of the saturated medium. 22 is the ratio of the shear wave velocity to the longitudinal wave velocity on the interface of the saturated medium, S 23 S is the square value of the ratio of the longitudinal wave velocity between the lower interface and the upper interface of the saturated medium. 24 is the ratio of the shear wave velocity of the lower interface to the longitudinal wave velocity of the upper interface of the saturated medium, S 32 is the ratio of the longitudinal wave velocity to the shear wave velocity on the interface of the saturated medium, S 33S is the square value of the ratio of the shear wave velocity between the lower interface and the upper interface of the saturated medium. 34 is the ratio of the product of the shear wave velocity of the lower interface of the saturated medium and the longitudinal wave velocity of the upper interface to the square of the longitudinal wave velocity of the upper interface, ρ is the density, θ1 is the incident angle, represents the P-wave velocity on the interface of saturated rock, and denote the shear wave velocities on the interface of dry rock and saturated rock, respectively. and represent the P-wave velocity and S-wave velocity of the interface under saturated rock, respectively. and They can be expressed as:

[0059]

[0060] Where, f is the fluid factor. According to Gassmann (1951), the Russell fluid factor can be expressed as

[0061]

[0062] Where φ represents the porosity, K dry is the bulk modulus of dry rock, K m represents the bulk modulus of solid minerals, K fl represents the bulk modulus of the saturated fluid in the rock. Assume that K fl Much smaller than K m , ignoring the above equations, the improved fluid coefficient F can be expressed as:

[0063] F=f / φ=K fl (1-K dry / K m ) 2 / φ 2 (7)

[0064] Combining equations (1) and (5, 6, 7) produces the exact Zeoppritz equation for saturated media, so equation (4) can be written as:

[0065]

[0066] Among them, F1 and F2 represent the fluid coefficient F of the upper and lower interfaces respectively, φ1 and φ2 represent the porosity φ of the upper and lower interfaces respectively, is the P-wave velocity of dry rock on the upper interface, is the P-wave velocity of dry rock at the lower interface, is the shear wave velocity of dry rock at the lower interface.

[0067] Step 5: Based on the Bayesian theory, the prior distribution function of the model parameters, the amplitude scaling factor and the initial parameter model, the inversion objective function in the sense of maximum a posteriori probability is constructed. By taking the derivative of the objective function with respect to the model parameters and solving the minimum value of the objective function, which is equivalent to setting its derivative equal to zero, the updated iterative porosity and fluid factor model parameter formulas are obtained by combining the iterative least squares algorithm.

[0068] The Bayesian method is a probability statistical method used to calculate conditional probability. Mathematically, its expression is:

[0069] P(m|d)∝P(d|m)P(m) (9)

[0070] Among them, P(m|d) is the posterior probability distribution; P(d|m) is the likelihood function that maps from the model parameter space to the observation data space; P(m) is the prior understanding of the model parameters before observation, called the prior distribution.

[0071] Assuming that the noises are independent of each other and the model parameters follow the Cauchy distribution, the posterior probability distribution of the multi-wave observation data can be expressed as

[0072]

[0073] Where superscript T represents the transpose of the matrix, Γ represents the nonlinear forward operator, and d pp is the actual observation data of longitudinal waves, Γ pp is the longitudinal wave synthetic data, d pS To convert the actual observed data of shear waves, Γ pS is to convert the shear wave synthetic data, m is the model parameter to be inverted, and the diagonal matrix χ = D T ξ -1 D is the statistical information related to the model parameters, ξ is the correlation matrix of the observed data, D is the difference matrix, ω and N m Respectively represent the mean vector of unknown model parameters and the length of unknown model parameters, τ=σ pp / σ ps is the weight associated with the PS wave, and υ=σ pp is the weight associated with the prior distribution, σ pp and σ ps are the noise variances of PP and PS waves respectively.

[0074] Based on perturbation theory and Taylor expansion, the generalized linear inversion strategy can transform the nonlinear problem into a linear problem by omitting the second-order terms and higher-order terms:

[0075]

[0076] Where d is the multi-wave observation data, Δd is the disturbance value of the observation data, Γ represents the nonlinear forward operator, and m is the model parameter to be inverted. is the first-order Taylor expansion of the model parameter perturbation value Δm to be inverted, O(Δm 2 ) are the higher-order infinitesimal terms of Taylor expansion that can be omitted.

[0077] At this point, the estimation of unknown model parameters can be transformed into solving the perturbation of model parameters:

[0078]

[0079] The inversion problem is converted into a direct solution to the maximum a posteriori probability, that is, to solve the minimum value of the posterior probability distribution or the objective function, that is, to find the solution where the derivative of the objective function with respect to the model parameters is zero. The simultaneous iterative least squares algorithm is simplified to obtain:

[0080] m k+1 =m k +η k (H k ) -1 (γ k ),k=0,1,2... (13)

[0081] Among them, m k+1 is the model parameter to be inverted at the k+1th iteration, m k is the model parameter to be inverted at the kth iteration, η k represents the step weight of the kth iteration, represents the Hessian matrix of the objective function,

[0082] The following uses actual data and well logging data to test the technical effect of the method of the present invention:

[0083] Figure 2 This is a CDP superimposed profile diagram. It can be seen from the figure that there is a strong reflection axis at the position of 2210-2255ms. Figure 3 is the elastic parameter inversion profile, from Figure 3 It can be seen that the elasticity parameter has relatively low values ​​at 2210-2255ms; Figure 4 is the inversion profile of porosity and fluid factor, from Figure 4 It can be seen that the porosity has a relatively low value at 2210-2255ms, while the fluid parameters have a relatively high value at 2210-2255ms; Figure 2 The strong reflection axis of the CDP stacked section has a good correspondence. At the same time, it can be seen that the inverted model parameter curve is very consistent with the corresponding actual logging curve, and the gas-bearing reservoir development section is 2210-2255ms.

[0084] In order to verify the advantages of multi-wave joint inversion over PP wave single inversion, PS wave pseudo gathers were synthesized using the saturated medium exact equation (Formulas 1 and 4) and logging data. Figure 5 Partial angle gathers from 8 to 20 degrees are shown. Figure 6 Comparison of the inversion curves of the logging position. It can be seen from the figure that compared with the PP wave single inversion result, the multi-wave joint inversion curve has better consistency with the logging curve, which proves that the multi-wave joint inversion method can effectively improve the inversion accuracy of the estimated parameters.

[0085] The second aspect of the present invention provides a porosity and fluid factor seismic synchronous inversion device, comprising:

[0086] The first processing unit is used to derive the exact Zeoppritz equation for isotropic saturated media;

[0087] The second processing unit is used to extract the seismic wavelet that depends on the incident angle based on the seismic data of the oil and gas reservoir, to simulate the seismic angle gather based on the well logging data and the accurate Zeoppritz equation of isotropic saturated medium, and to determine the amplitude scaling factor in combination with the actual wellside seismic data;

[0088] The third processing unit is used to calculate the mean of each model parameter based on all the logging data in the work area, obtain the autocorrelation coefficient and the mutual correlation coefficient of each model parameter, construct the covariance matrix related to the parameters, and form a priori distribution function of the model parameters that conforms to the work area;

[0089] The fourth processing unit is used to perform smoothing processing on the logging data and extrapolate along the interpretation layer to obtain an initial parameter model;

[0090] The fifth processing unit is used to construct an inversion objective function in the sense of maximum a posteriori probability based on Bayesian theory, model parameter prior distribution function, amplitude scaling factor and initial parameter model, and solve the inversion objective function to obtain updated iterative porosity and fluid factor model parameter formulas.

[0091] A third aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-mentioned porosity and fluid factor seismic synchronous inversion method.

[0092] A fourth aspect of the present invention provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the above-mentioned porosity and fluid factor seismic synchronous inversion method when executing the computer program.

[0093] The present invention is described in terms of flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to specific embodiments. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0094] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0095] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0096] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A porosity and fluid factor seismic synchronous inversion method, characterized in that: The steps include: Derive the exact Zeoppritz equation for isotropic saturated media; Seismic data from oil and gas reservoirs are used to extract seismic wavelets that depend on the incident angle, and to simulate seismic angle gathers based on well logging data and the accurate Zeoppritz equation for isotropic saturated media, and the amplitude scaling factor is determined in combination with actual wellside seismic data. Based on all the logging data in the work area, the mean of each model parameter is calculated, the autocorrelation coefficient and the mutual correlation coefficient of each model parameter are obtained, and the covariance matrix related to the parameters is constructed to form a priori distribution function of the model parameters that conforms to the work area; The logging data are smoothed and extrapolated along the interpreted horizon to obtain the initial parameter model; Based on Bayesian theory, model parameter prior distribution function, amplitude scaling factor and initial parameter model, the inversion objective function in the sense of maximum a posteriori probability is constructed, and the inversion objective function is solved to obtain updated iterative porosity and fluid factor model parameter formulas. The exact Zeoppritz equation for an isotropic saturated medium is given by: Among them, S 12 S is the square value of the ratio of the shear wave velocity to the longitudinal wave velocity of the dry rock on the upper interface; 13 S is the ratio of the longitudinal wave velocity between the lower interface and the upper interface of the saturated medium; 14 S is the square value of the ratio of the shear wave velocity of the lower interface to the longitudinal wave velocity of the upper interface of the saturated medium; 22 is the ratio of the shear wave velocity to the longitudinal wave velocity on the interface of the saturated medium; S 23 S is the square value of the ratio of the longitudinal wave velocity between the lower interface and the upper interface of the saturated medium; 24 S is the ratio of the shear wave velocity of the lower interface to the longitudinal wave velocity of the upper interface of the saturated medium; 32 is the ratio of the longitudinal wave velocity to the shear wave velocity on the interface of the saturated medium; S 33 S is the square value of the ratio of the shear wave velocity between the lower interface and the upper interface of the saturated medium; 34 It is the ratio of the product of the shear wave velocity of the lower interface of the saturated medium and the longitudinal wave velocity of the upper interface to the square of the longitudinal wave velocity of the upper interface; represents the shear wave velocity of the upper interface of dry rock; F1 and F2 represent the fluid coefficient F of the upper and lower interfaces respectively; φ1 and φ2 represent the porosity of the upper and lower interfaces respectively. is the P-wave velocity of dry rock on the upper interface, is the P-wave velocity of dry rock at the lower interface, is the shear wave velocity of dry rock at the lower interface; ρ1 is the density.

2. The porosity and fluid factor seismic synchronous inversion method according to claim 1, characterized in that: Seismic data from oil and gas reservoirs are used to extract seismic wavelets that depend on the incident angle, and the seismic angle gathers are simulated based on the accurate Zeoppritz equation of well logging data and isotropic saturated media. The amplitude scaling factor is determined in combination with the actual wellside seismic data. Specifically, the following are included: Assuming that the seismic wavelet is known before inversion, it is necessary to extract the wavelet using statistical methods based on the actual pre-stack seismic gathers and logging data. The wavelet will change in waveform or frequency due to the influence of the formation during propagation. Extracting the seismic wavelet that depends on the incident angle can effectively improve the amplitude matching degree. Taking logging data as input, the model uses the accurate Zeoppritz equation of isotropic saturated medium to forward simulate the PP wave gathers in the angle domain. 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.

3. The porosity and fluid factor seismic synchronous inversion method according to claim 1, characterized in that: The prior distribution function of the model parameters is: in, is the parameter model, ρ is the density, φ is the porosity, F is the improved fluid coefficient, is the longitudinal wave velocity of dry rock, Dry rock shear wave velocity; μ and C m are the mean of the model parameters and the covariance matrix of the model parameters respectively.

4. The porosity and fluid factor seismic synchronous inversion method according to claim 1, characterized in that: The logging data is smoothed and extrapolated along the interpretation horizon to obtain the initial parameter model, as follows: First, the scattered point interpolation method is used to interpolate the data of each layer to complete the geological layer modeling. Then, the model parameters are transversely interpolated according to the geological layer, that is, the logging information is transversely interpolated to calculate the model parameter value at each point underground, and the establishment of the initial new attribute parameter model is completed.

5. The porosity and fluid factor seismic synchronous inversion method according to claim 1, characterized in that: The specific solution process of the inversion objective function is as follows: By taking the derivative of the inversion objective function with respect to the model parameters and solving the minimum value of the inversion objective function, which is equivalent to setting its derivative equal to zero, the updated iterative porosity and fluid factor model parameter formulas are obtained by combining the iterative least squares algorithm.

6. A porosity and fluid factor seismic synchronous inversion device, characterized in that: include: The first processing unit is used to derive the accurate Zeoppritz equation of isotropic saturated medium. The accurate Zeoppritz equation of isotropic saturated medium is as follows: Among them, S 12 S is the square value of the ratio of the shear wave velocity to the longitudinal wave velocity of the dry rock on the upper interface; 13 S is the ratio of the longitudinal wave velocity between the lower interface and the upper interface of the saturated medium; 14 S is the square value of the ratio of the shear wave velocity of the lower interface to the longitudinal wave velocity of the upper interface of the saturated medium; 22 is the ratio of the shear wave velocity to the longitudinal wave velocity on the interface of the saturated medium; S 23 S is the square value of the ratio of the longitudinal wave velocity between the lower interface and the upper interface of the saturated medium; 24 S is the ratio of the shear wave velocity of the lower interface to the longitudinal wave velocity of the upper interface of the saturated medium; 32 is the ratio of the longitudinal wave velocity to the shear wave velocity on the interface of the saturated medium; S 33 S is the square value of the ratio of the shear wave velocity between the lower interface and the upper interface of the saturated medium; 34 It is the ratio of the product of the shear wave velocity of the lower interface of the saturated medium and the longitudinal wave velocity of the upper interface to the square of the longitudinal wave velocity of the upper interface; represents the shear wave velocity of the upper interface of dry rock; F1 and F2 represent the fluid coefficient F of the upper and lower interfaces respectively; φ1 and φ2 represent the porosity of the upper and lower interfaces respectively. is the P-wave velocity of dry rock on the upper interface, is the P-wave velocity of dry rock at the lower interface, is the shear wave velocity of dry rock at the lower interface; ρ1 is the density; The second processing unit is used to extract the seismic wavelet that depends on the incident angle based on the seismic data of the oil and gas reservoir, to simulate the seismic angle gather based on the well logging data and the accurate Zeoppritz equation of isotropic saturated medium, and to determine the amplitude scaling factor in combination with the actual wellside seismic data; The third processing unit is used to calculate the mean value of each model parameter based on all the logging data in the work area, obtain the autocorrelation coefficient and the mutual correlation coefficient of each model parameter, construct a covariance matrix related to the parameters, and form a priori distribution function of the model parameters that conforms to the work area; The fourth processing unit is used to perform smoothing processing on the logging data and extrapolate along the interpretation layer to obtain an initial parameter model; The fifth processing unit is used to construct an inversion objective function in the sense of maximum a posteriori probability based on Bayesian theory, model parameter prior distribution function, amplitude scaling factor and initial parameter model, and solve the inversion objective function to obtain updated iterative porosity and fluid factor model parameter formulas.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the porosity and fluid factor seismic synchronous inversion method described in any one of claims 1 to 5 are implemented.

8. A computer device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of the porosity and fluid factor seismic synchronous inversion method described in any one of claims 1 to 5 are implemented.

Images