Reservoir physical property parameter random inversion method and system based on mixed prior information

Abstract

The invention relates to a reservoir physical property parameter random inversion method and system based on mixed prior information, and the method comprises the steps: building a seismic wave constrained low-frequency prior information model based on the similarity simulation of seismic waveforms; based on a fractal theory, constructing a high-frequency prior information model of seismic waves; combining the low-frequency prior information model and the high-frequency prior information model to construct a mixed prior information model; establishing a relationship between the reservoir physical property parameters and the elastic parameters, and combining the relationship with the mixed prior information model to obtain a reservoir physical property parameter inversion model of the seismic data; the reservoir physical property parameter inversion model is disturbed and iterated through a QA-MCMC algorithm, and when the least square error of the synthetic seismic record and the well bypass seismic record is minimum, a final reservoir physical property parameter inversion model is generated; and inputting seismic data to be measured into the final reservoir physical property parameter inversion model to obtain a reservoir physical property parameter inversion result. The problem of error accumulation during reservoir physical property parameter prediction by an indirect inversion method is solved.

Description

Technical Field

[0001] The invention relates to a reservoir physical property parameter random inversion method and system based on mixed prior information, belonging to the technical field of seismic wave inversion in geological exploration. Background Art

[0002] Reservoir physical parameters, such as permeability, porosity, fluid saturation and mud content, are important indicators for describing reservoir characteristics, fluid patterns and reservoir modeling. Their accurate estimation can provide a strong reference for reservoir prediction, improve the description accuracy of reservoir prediction, and reduce the cost of oil and gas exploration. It is a key research issue in current oil and gas exploration. Usually, reservoir physical parameters can only be calculated from logging data or obtained from core testing. Although logging data can provide information about parameters such as reservoir thickness, lithology, porosity and water saturation, it lacks spatial continuity.

[0003] Seismic wave impedance inversion methods can be divided into direct inversion and indirect inversion. Direct inversion refers to the method of inverting underground structures directly based on seismic wave observation data. Common direct inversion methods include full waveform inversion. Indirect inversion refers to the method of inverting underground structures by building models and calculating seismic wave propagation paths. Common indirect methods include tomography, regularized inversion, and genetic algorithms. At present, the indirect inversion method is mainly used to predict reservoir physical properties. However, there is a problem of error accumulation in predicting reservoir physical properties using the indirect inversion method, which causes the prediction results to be inconsistent with the actual situation. Summary of the invention

[0004] In view of the above problems, the purpose of the present invention is to provide a random inversion method and system for reservoir physical property parameters based on mixed prior information, which solves the error accumulation problem when predicting reservoir physical property parameters by indirect inversion method.

[0005] To achieve the above purpose, the present invention proposes the following technical solution: a random inversion method for reservoir physical property parameters based on mixed prior information. It includes the following steps:

[0006] Low-frequency prior information modeling: Based on the similarity simulation of seismic waveforms, a low-frequency prior information model constrained by seismic waves is constructed;

[0007] High-frequency prior information modeling: Based on fractal theory, a high-frequency prior information model of seismic waves is constructed;

[0008] Hybrid prior information modeling: combining low-frequency and high-frequency prior information models to construct a hybrid prior information model that comprehensively reflects reservoir parameter characteristics;

[0009] Reservoir physical property parameter inversion modeling: establishing the relationship between reservoir physical property parameters and elastic parameters, combining it with the hybrid prior information model, and obtaining a reservoir physical property parameter inversion model of seismic data;

[0010] Model optimization iteration: The reservoir physical property parameter inversion model is disturbed and iterated by using the QA-MCMC algorithm, and when the minimum square error between the synthetic seismic record and the well bypass seismic record is the minimum, the final reservoir physical property parameter inversion model is generated;

[0011] Result output: input the seismic data to be measured into the final reservoir physical property parameter inversion model to obtain the reservoir physical property parameter inversion result.

[0012] Furthermore, based on the similarity simulation of seismic waveforms, a method for constructing low-frequency prior information constrained by seismic waves is as follows: segmenting the well curve, wellside seismic data and measured seismic records, and establishing a well sample set, a wellside seismic trace sample set and a seismic data set to be inverted respectively; calculating the correlation coefficient between the known samples and the seismic data at the location to be inverted; screening all well logging physical property parameter samples from high to low according to the correlation coefficient; selecting the well logging physical property parameters with the correlation coefficient higher than a threshold; assigning different weights to the selected well logging physical property parameters, bringing them into the Kriging interpolation formula, and obtaining the Kriging parameters of the location to be inverted; using the Kriging parameters, combined with sequential Gaussian simulation, performing random sampling at each of the locations to be inverted to complete a similarity simulation; performing multiple similarity simulations, averaging the results of multiple similarity simulations, and obtaining prior information based on seismic waveform constraints.

[0013] Further, the calculation formula of the correlation coefficient is:

[0014]

[0015] Among them, r ij is x i The earthquake waveform and x j The similarity coefficient of the earthquake waveform at i ) is x i The variance of the variable at j ) is x j The variance of the variable, C ij is x l Place and x j The spatial covariance value at ;

[0016] After substituting the correlation coefficient into the Kriging variance formula, the variance within the region is obtained. The formula is:

[0017]

[0018] Among them, ω i and ω j are the weight coefficients in the Kriging interpolation formula, σ is the standard deviation of the measured value, and n is the number of sample points;

[0019] Find the variance within the region The minimum x j Solution: Combined with the constraints of unbiased estimation, construct the objective function, solve the coefficient matrix of the objective function, and obtain the weight coefficient ω based on the seismic waveform similarity coefficient i ; Substitute it into the ordinary Kriging interpolation formula and the Kriging variance formula to obtain the Kriging parameters of sequential Gaussian simulation. According to the random path, calculate the probability density cumulative distribution function of each position to be inverted in turn and perform random sampling to obtain the simulated values ​​of the physical property parameters.

[0020] Furthermore, based on fractal theory, a method for constructing high-frequency prior information of seismic waves is as follows: using a fractional Gaussian distribution model to simulate the frequency range of logging data and construct an initial high-frequency prior information model; using an autocovariance function as a fractional Gaussian distribution to describe logging curves based on logging data; performing Fourier transform on the autocovariance function to obtain a power spectral density function; generating a random number sequence having the same distribution characteristics as the standardized logging data, and convolving the random number sequence with the power spectral density function to obtain a high-frequency prior information model of seismic waves.

[0021] Furthermore, the high-frequency prior information model of seismic waves is:

[0022]

[0023] Among them, Y(k) is the spectrum value after Fourier transformation, M is the maximum number or range of random number generation, c is a constant, N is the total amount of data, and V(k) is the standardized random number.

[0024] Furthermore, the hybrid prior information model is:

[0025] m=αm 1 +(1-α)m 2

[0026] Among them, m 1 is the low-frequency prior information; m 2 is the high-frequency prior information; α is the weight coefficient of the prior initial model obtained based on the seismic waveform similarity simulation; m is the mixed prior information.

[0027] Furthermore, the relationship between the reservoir physical property parameters and the elastic parameters is:

[0028] d=f(m)+ε=f(g(r))+e

[0029] Among them, d is the reservoir physical parameter, ε is the error term, m is the elastic parameter, r is the rock physical parameter, g is the rock physics model, f is the seismic forward model, e is the mean value of 0, and the covariance is ∑ eThe data error of Gaussian distribution is e~N(0,∑ e ).

[0030] Furthermore, the method of perturbing and iterating the reservoir physical property parameter inversion model by using the QA-MCMC algorithm is as follows: according to the parameter distribution characteristics at the current moment in the Markov chain process that are only related to the previous moment, a physical property parameter distribution model is established; the relationship between the physical property parameter distribution model at the next moment and the initial parameter distribution model is obtained, which is used as a physical property parameter disturbance, and whether to accept the physical property parameter distribution model at the next moment is judged by the Metropolis sampling criterion, and the distribution model is optimized by the likelihood function to obtain a physical property parameter state model after the MCMC inversion algorithm; the QA algorithm is introduced into the physical property parameter state model to limit the range of the solution space and obtain the final reservoir physical property parameter inversion model.

[0031] Furthermore, the likelihood function uses low-frequency prior information and high-frequency prior information as constraints, and assigns different weights to different prior information items; the final reservoir physical property parameter inversion model is:

[0032]

[0033] in, is the sample point x t and target point The similarity weight between t ) is the sample point x t The probability density function of is the sample point The probability density function of x t is the sample point, is the target point, ex is the residual of the data, and φ(v) is the probability density function.

[0034] The present invention also discloses a reservoir physical property parameter random inversion system based on hybrid prior information, comprising: a low-frequency prior information model establishment module, used for constructing a low-frequency prior information model constrained by seismic waves according to similarity simulation of seismic waveforms; a high-frequency prior information model establishment module, used for constructing a high-frequency prior information model of seismic waves according to fractal theory; a hybrid prior information model establishment module, used for combining the low-frequency prior information model with the high-frequency prior information model to construct a hybrid prior information model; a physical property elastic parameter relationship module, used for establishing a relationship between reservoir physical property parameters and elastic parameters, combining the relationship with the hybrid prior information model to obtain a reservoir physical property parameter inversion model of seismic data; a QA-MCMC algorithm module, used for perturbing and iterating the reservoir physical property parameter inversion model through a QA-MCMC algorithm, and generating a final reservoir physical property parameter inversion model when the minimum square error between the synthetic seismic record and the well bypass seismic record is the smallest; and a result output module, used for inputting the seismic data to be measured into the final reservoir physical property parameter inversion model to obtain a reservoir physical property parameter inversion result.

[0035] The technical solution of the present invention has at least the following technical effects or advantages:

[0036] 1. The present invention successfully realizes direct seismic inversion of reservoir physical parameters. When constructing the initial model, the constraints of seismic waveform characteristics are fully considered, and high-frequency components are added to improve the resolution of the model.

[0037] 2. The present invention replaces the traditional variogram with the similarity coefficient of the seismic waveform, thereby achieving more accurate underground reservoir parameter interpolation and random simulation. It fully integrates the information carried by the seismic waveform into the interpolation and simulation results, thereby improving the accuracy of the results. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 is a flow chart of a random inversion method for reservoir physical property parameters based on mixed prior information in one embodiment of the present invention;

[0039] Figure 2 is a flow chart of a seismic waveform similarity simulation method in one embodiment of the present invention;

[0040] Figure 3 is a schematic diagram of similar simulation of earthquake waveforms in one embodiment of the present invention;

[0041] Figure 4 is a diagram of reservoir porosity inversion results in one embodiment of the present invention;

[0042] Figure 5 is a diagram of reservoir water saturation inversion results in one embodiment of the present invention;

[0043] Figure 6This is a diagram of reservoir mud content inversion results in one embodiment of the present invention. DETAILED DESCRIPTION

[0044] In order to enable those skilled in the art to better understand the technical solution of the present invention, the present invention is described in detail through specific embodiments. However, it should be understood that the provision of specific embodiments is only for a better understanding of the present invention, and they should not be understood as limitations of the present invention. In the description of the present invention, it should be understood that the terms used are only for the purpose of description and cannot be understood as indicating or implying relative importance.

[0045] In order to solve the problem of error accumulation in the indirect inversion of reservoir physical property parameters, the present invention proposes a random inversion method and system for reservoir physical property parameters based on mixed prior information, a formula for constructing a prior model of physical property parameters based on seismic waveform constraints, and a high-frequency prior model is obtained by using fractal theory, and the two are linearly weighted to obtain a mixed prior model of physical property parameters that can better reflect the spatial structure of the strata and the details of the geological body; the law of changes in seismic response characteristics with physical property parameters is analyzed, a linear or rock physics model is established to characterize the relationship between elasticity and physical property parameters, and a reservoir physical property parameter inversion formula based on seismic data is derived, which successfully realizes direct seismic inversion of reservoir physical property parameters. In view of the dependence of random inversion on the initial model, the method for constructing the initial model of physical property parameters is improved. When constructing the initial model, the constraints of seismic waveform characteristics are fully considered, and high-frequency components are added to improve the resolution of the model. The scheme of the present invention is described in detail below through embodiments in conjunction with the accompanying drawings.

[0046] Embodiment 1

[0047] This embodiment discloses a random inversion method for reservoir physical property parameters based on mixed prior information, such as Figure 1 As shown, the following steps are included:

[0048] S1: Based on the similarity simulation of seismic waveforms, a low-frequency prior information model constrained by seismic waves is constructed.

[0049] like Figure 2 As shown in the figure, first, the well curve, wellside seismic data and measured seismic records are segmented, and the well sample set, wellside seismic trace sample set and seismic data set to be inverted are established respectively. The well sample set records the reservoir physical property parameters corresponding to the well logging data; the wellside seismic trace sample set extracts the seismic waveform of the wellside seismic trace; the seismic data set to be inverted: extracts the seismic waveform of the target reservoir.

[0050] Subsequently, the correlation coefficient between the known sample and the seismic data at the location to be inverted is calculated, and all logging physical parameter samples are screened from high to low according to the correlation coefficient. In the prior art, the correlation coefficient between the sample and the seismic waveform at the location to be inverted is:

[0051]

[0052] Among them, r ij is the correlation coefficient between the sample and the seismic waveform at the location to be inverted; D(x i ) is the seismic waveform of the sample point; D(x j ) is the seismic waveform of the point to be inverted.

[0053] The well logging physical property parameters with correlation coefficients higher than the threshold are selected; different weights are assigned to the selected well logging physical property parameters, and the parameters are introduced into the Kriging interpolation formula to obtain the Kriging parameters of the position to be inverted, namely the Kriging estimate and the Kriging covariance.

[0054] The method to obtain Kriging parameters is:

[0055] The existing Kriging interpolation formula is:

[0056]

[0057] Among them, σ 2 is the variance within the region; C ij =C(x i ,x j ) is x i Place and x j The spatial covariance value at is usually obtained through the variogram, where Var is the variance of the data, x 0 is the parameter value at the initial position, is the parameter value after Kriging interpolation, ω i and ω j are weight coefficients in the Kriging interpolation formula, and n is the number of sample points.

[0058] In this embodiment, a method based on seismic waveform similarity simulation is proposed. By calculating the similarity coefficient of the seismic waveforms corresponding to two sets of physical parameters, C is obtained by replacing the variation function with the similarity coefficient. ij , the calculation formula of the similarity coefficient of the seismic waveform is as follows:

[0059]

[0060] Among them, r ij is x i The earthquake waveform and x j The similarity coefficient of the earthquake waveform at i ) is x i The variance of the variable at j ) is x j The variance of the variable, C ij is x i Place and x j The spatial covariance value at .

[0061] Substituting the similarity coefficient of the seismic waveform into the existing Kriging interpolation formula, the Kriging interpolation formula in this embodiment is obtained as follows:

[0062]

[0063] Find the x where the variance is the smallest in the Kriging interpolation formula j The solution; Combined with the constraints of unbiased estimation, construct the objective function;

[0064]

[0065] Where λ is the Lagrange multiplier, ω i is the weight coefficient of the similarity coefficient of the seismic waveform.

[0066] Solve the coefficient matrix of the objective function and obtain the weight coefficient ω based on the seismic waveform similarity coefficient i ; Substitute it into the ordinary Kriging interpolation formula and the Kriging variance formula to obtain the Kriging parameters of sequential Gaussian simulation.

[0067] Through the Kriging parameters, combined with sequential Gaussian simulation, random sampling is performed at each location to be inverted to complete a similarity simulation, that is, according to the random path, the probability density cumulative distribution function of each location to be inverted is calculated in turn and randomly sampled to obtain the simulated value of the physical property parameters. Multiple similarity simulations are performed, and the results of multiple similarity simulations are averaged to obtain prior information based on seismic waveform constraints.

[0068] S2 constructs a high-frequency prior information model of seismic waves based on fractal theory.

[0069] Due to the complexity and heterogeneity of underground media, the assumption of Gaussian distribution usually cannot accurately describe the distribution of underground media. Therefore, fractional Gaussian distribution is used to construct an initial model that conforms to the logging frequency range to optimize the underground media model.

[0070] In view of the complexity of underground media, fractional Gaussian distribution is used to describe the detailed changes of reservoir physical parameters. The fractional Gaussian distribution consists of mean μ, variance σ 2 and Hurst coefficient h. The power exponent of the covariance function and power spectrum density function of the logging data is the Hurst coefficient. Therefore, the autocovariance function R(t) with a mean of 0 can be used as a fractional Gaussian distribution to describe the logging curve; the formula of the autocovariance function is:

[0071] R(t)=0.5σ 2 (|t+1| 2h -2|t| 2h +|t-1| 2h )

[0072] According to the Wiener-Khinchin theorem, the autocovariance function is transformed by Fourier transform to obtain the power spectrum density function S(k). The formula of the power spectrum density function S(k) is:

[0073]

[0074] Generate a random number sequence V(k) with the same distribution characteristics as the standardized logging data. The calculation formula of the random number sequence V(k) is:

[0075]

[0076] Where W(k) is an independent Gaussian distributed random number.

[0077] If the mean of the standardized logging data is μ, the Gaussian random number independently distributed with the power spectrum density function S(k) is W(k)k=0,1…,M-1, and it is standardized to obtain V(k):

[0078]

[0079] Finally, the first N elements of the discrete Fourier transform of V(k) are calculated to obtain the simulated sequence Y(k), which is the high-frequency prior information model:

[0080]

[0081] Among them, Y(k) is the Fourier spectrum value of the seismic wave, M is the size of the data set, c is a constant, N is the total amount of data, V(k) is the standardized random number, V * It is the weight factor finally generated in the mixed prior information model.

[0082] S3 combines the low-frequency prior information model and the high-frequency prior information model to construct a hybrid prior information model.

[0083] Combining seismic waveform similarity simulation with fractal theory, we can obtain mixed prior information that can reflect the overall trend of the stratum and precisely characterize the stratum. The mixed prior information model m is:

[0084] m=αm 1 +(1-α)m 2

[0085] Among them, m 1 is the low-frequency prior information; m 2 is the high-frequency prior information; α is the weight coefficient of the prior initial model obtained based on the seismic waveform similarity simulation. The synthetic seismic records of the mixed prior model under different weight coefficients are obtained respectively, and the weight coefficient with the highest correlation between the synthetic seismic record and the measured seismic record is selected.

[0086] S4 establishes the relationship between reservoir physical property parameters and elastic parameters, combines it with the hybrid prior information model, and obtains the reservoir physical property parameter inversion model of seismic data.

[0087] The relationship between physical parameters and elastic parameters is characterized by rock physics model or logging curve analysis, and a direct relationship between reservoir physical parameters and seismic response is established, thereby achieving direct inversion of reservoir physical parameters. The relationship between reservoir physical parameters and elastic parameters is:

[0088] d=f(m)+ε=f(g(r))+e

[0089] Where d is the reservoir physical parameter, ε is the error term, m is the elastic parameter, r is the rock physical parameter, generally porosity, shale content and water saturation; g is the rock physics model, which is used to describe the relationship between rock elastic parameters and physical parameters; f is the seismic forward model, which is used to describe the relationship between rock elastic parameters and synthetic seismic records; e is the service mean of 0 and the covariance of ∑ e The Gaussian distribution of the data error eN(0,∑ e ).

[0090] If the relationship between reservoir physical properties and elastic parameters can be approximated as a linear relationship, it can be expressed in the form of a multivariate linear equation as follows:

[0091] V P =α 1 ·φ+α 2 ·V sh +α 3 ·S w +δ 1

[0092] V S =β 1 ·φ+β 2 ·V sh +β 3 ·S w +δ 2

[0093] ρ=γ 1 ·φ+γ 2 ·V sh +γ 3 ·S w +δ 3

[0094] Among them, V P is the longitudinal wave velocity, V S is the shear wave velocity, ρ is the density, φ is the porosity, V sh is the clay content, S w is the water saturation, α 1, α 2 and α 3 is the linear regression coefficient of porosity, β 1 , β 2 and β 3 is the linear regression coefficient of water saturation, γ 1 , γ 2 and γ 3 is the linear regression coefficient of mud content, δ 1 , δ 2 and δ 3 is the adjustment factor describing the parameter weight.

[0095] The formulas for rock physical parameters and seismic records derived based on statistical rock physics relationships and the Aki-Richard approximation formula are as follows:

[0096]

[0097] Among them, WAD and WADP are statistical parameters that characterize reservoir physical properties, θ is the incident angle, and n t and n e is the sample point, d and e are the description function and error function of reservoir physical property parameters respectively.

[0098] For example, if the P-wave velocity is approximately a linear combination of the three physical parameters of porosity, water saturation and shale content, the relationship between the S-wave velocity, density and the three physical parameters can be obtained. Therefore, the relationship between reservoir physical parameters and elastic parameters can be expressed as a simple linear relationship based on the rock physics modeling method of data fitting, and then the relationship between physical parameters and seismic response can be established to achieve direct inversion of reservoir physical parameters based on seismic data. The multivariate linear equation group of physical parameters and elastic parameters is as follows:

[0099] V P =-2618·φ-1709·V sh +1143·S w +3588

[0100] V S =-1767·φ-1072·V sh +518·S w +2080

[0101] ρ=-1620·φ-70·V sh +280·S w +2371.

[0102] S5 perturbs and iterates the reservoir physical property parameter inversion model through the QA-MCMC algorithm. When the minimum square error between the synthetic seismic record and the well bypass seismic record is the smallest, the final reservoir physical property parameter inversion model is generated.

[0103] The parameter distribution characteristics of the current moment in the Markov chain process are only related to the previous moment, which can be expressed as:

[0104] P(x t |x t-1 )∝L(x t )·π(x t )

[0105] Among them, L(x t ) is the likelihood function, π(x t ) is the prior distribution.

[0106] According to the fact that the parameter distribution characteristics of the current moment in the Markov chain process are only related to the previous moment, a physical property parameter distribution model is established, and its calculation formula is:

[0107] P(X t =s j |X t-1 =s j-1 ,X t-2 =s j-2 ,...,X 0 =s 0 )=P(X t =s j |X t-1 =s j-1 )

[0108] Where P(*) is a Markov chain, X i is a random variable, s j is the parameter distribution characteristic at state j.

[0109] The relationship between the physical property parameter distribution model at the next moment and the initial parameter distribution model is obtained and used as the physical property parameter disturbance. Specifically, assuming that the state space of the parameters at time t is π(t), the state space of the parameters at time t+1 can be expressed by the state space at time t, and by analogy, the relationship between the physical property parameter distribution model at the next moment and the initial parameter distribution model is:

[0110] π(t+1)=π(t)P t =π(t-1)P t-1 P t =L =π(0)P 0 LP t-1 P t

[0111] Among them, L is the objective function.

[0112] The Metroplis sampling criterion is used to judge whether to accept the physical property parameter distribution model at the next moment, and the acceptance probability formula after the physical property parameter disturbance is obtained. The proposed distribution given a symmetric distribution can simplify the acceptance probability formula, and the distribution model is optimized through the likelihood function. When constructing the likelihood function, low-frequency prior information and high-frequency prior information are used as constraints, and different prior information items are given different weights to improve the reliability of the inversion results.

[0113] Based on the Metropolis criterion, the probability of accepting the physical parameter distribution model at the next moment is:

[0114]

[0115] Among them, P accept is the probability of accepting the physical parameter distribution model at the next moment, L(x new ) is the likelihood function of the physical parameter distribution model at the next moment, L(x old ) is the likelihood function of the physical property parameter distribution model at this time.

[0116] The likelihood function expression for reservoir physical property inversion is as follows, in which the physical property parameters to be inverted and the logging data are added as constraints.

[0117]

[0118] Where L(s) is the likelihood estimate of the current state s, k is the total number of data points, N is the regularization factor used to control the complexity of the model, and d i It is a synthetic seismic record obtained by two-step forward modeling of physical parameters; is the measured earthquake record, R i is the a priori initial model of the physical parameters to be inverted, is the a priori initial model established based on well logging seismic data, μ i is the mean value of the physical property parameter to be inverted, is the mean of the measured earthquake data, σ i is the standard deviation of the physical property parameter to be inverted, is the standard deviation of the measured seismic data, h i is the Hurst coefficient of the physical property parameter to be inverted, is the Hurst coefficient of the measured seismic data, α 1 , α 2 , α 3 and α 4 are the weight coefficients of constraint information, which are determined according to the relative position of the inversion parameters and the wells, σ 2In this embodiment, when constructing the likelihood function, the constraints of prior information and logging data are added, which is conducive to accelerating the convergence speed of the inversion algorithm.

[0119] Obtain the physical property parameter state model after the MCMC inversion algorithm; introduce the QA algorithm into the physical property parameter state model to improve the initial convergence speed and limit the scope of the solution space to obtain the final reservoir physical property parameter inversion model. The final reservoir physical property parameter inversion model is:

[0120]

[0121] in, is the sample point x t and target point The similarity weight between t ) is the sample point x t The probability density function of is the sample point The probability density function of x t is the sample point, is the target point, ex is the residual of the data, and φ(v) is the probability density function.

[0122] Based on the prior model of physical property parameters and partial angle stacking gathers, the inversion objective function of physical property parameters is constructed, and the mixed prior model of physical property parameters is perturbed and iterated by using the QA-MCMC optimization algorithm. The initial temperature of the iteration is set to 200, and the annealing coefficient is set to 0.999. The initial step sizes of the three physical property parameter iterations are 0.01, 0.032, and 0.08, respectively, and the step size attenuation factor is 0.99. The minimum square error between the synthetic seismic record and the well bypass seismic record, as well as the relevant parameters of the mixed prior information of physical property parameters (such as mean, variance, Hurst coefficient, and low-frequency model, etc.) are obtained as constraints to construct the inversion objective function. The iteration ends when the temperature drops to 1e-5 or the objective function is less than 1e-6.

[0123] Figure 4 The black box is the target reservoir, and the well logging lithology interpretation is a 9m sand body, which shows high porosity. Comparing the well logging porosity curve with the inversion results at the well location, the inversion results of the high porosity layer in the well logging also show high porosity, and the low porosity layer corresponds to the low porosity value of the inversion result. Figure 5 This is the water saturation inversion result, which is also consistent with the well logging water saturation curve. Figure 6The inversion result of shale content. The high shale content layer corresponds to the high shale content value of the inversion result, and the logging interpretation is mudstone; the low shale content layer corresponds to the low shale content value of the inversion result, and the logging interpretation is a sandstone reservoir. Combining the logging interpretation curve and the inversion results of physical parameters, the reservoir section in the black box shows high porosity, low water saturation and low shale content, and it is inferred that the reservoir may be an oil-water mixed layer.

[0124] S6 inputs the seismic data to be measured into the final reservoir physical property parameter inversion model to obtain the reservoir physical property parameter inversion results.

[0125] Based on the statistical rock physics relationship and the AVO approximate formula, this embodiment derives the expression between the reservoir physical property parameters and the seismic response characteristics, and uses the random inversion method to achieve the direct seismic inversion of reservoir physical property parameters. Compared with the indirect inversion method of physical property parameters using elastic impedance or elastic parameters, the former can effectively avoid the cumulative error in the indirect inversion process and improve the accuracy of the inversion results of physical property parameters.

[0126] Embodiment 2

[0127] Based on the same inventive concept, this embodiment discloses a reservoir physical property parameter stochastic inversion system based on mixed prior information, comprising:

[0128] A low-frequency prior information model building module is used to build a low-frequency prior information model constrained by seismic waves based on similarity simulation of seismic waveforms;

[0129] A high-frequency prior information model building module is used to build a high-frequency prior information model of seismic waves based on fractal theory;

[0130] A hybrid prior information model building module is used to combine the low-frequency prior information model and the high-frequency prior information model to build a hybrid prior information model;

[0131] The physical property elastic parameter relationship module is used to establish the relationship between reservoir physical property parameters and elastic parameters, and combine it with the hybrid prior information model to obtain the reservoir physical property parameter inversion model of seismic data;

[0132] QA-MCMC algorithm module, used to perturb and iterate the reservoir physical property parameter inversion model through the QA-MCMC algorithm, and generate the final reservoir physical property parameter inversion model when the minimum square error between the synthetic seismic record and the well bypass seismic record is the smallest;

[0133] The result output module is used to input the measured seismic data into the final reservoir physical property parameter inversion model to obtain the reservoir physical property parameters. It should be understood by those skilled in the art that the embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program codes.

[0134] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. 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 produce 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.

[0135] 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.

[0136] 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.

[0137] 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 them. Although the present invention has been described in detail with reference to the above embodiments, ordinary technicians in the relevant field should understand that the specific implementation of the present invention can still be modified or replaced by equivalents, and any modification or equivalent replacement that does not deviate from the spirit and scope of the present invention should be included in the protection scope of the claims of the present invention. The above content is only a specific implementation of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of changes or replacements within the technical scope disclosed by the present invention, which should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A random inversion method for reservoir physical property parameters based on mixed prior information, characterized in that: The following steps are involved: Based on the similarity simulation of seismic waveforms, a low-frequency prior information model constrained by seismic waves is constructed; Based on fractal theory, a high-frequency prior information model of seismic waves is constructed; Combining the low-frequency prior information model and the high-frequency prior information model to construct a hybrid prior information model; Establishing the relationship between reservoir physical property parameters and elastic parameters, combining it with the hybrid prior information model, and obtaining a reservoir physical property parameter inversion model of seismic data; The reservoir physical property parameter inversion model is disturbed and iterated by using a QA-MCMC algorithm, and when the minimum square error between the synthetic seismic record and the well bypass seismic record is the minimum, a final reservoir physical property parameter inversion model is generated; The seismic data to be measured is input into the final reservoir physical property parameter inversion model to obtain the reservoir physical property parameter inversion result.

2. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 1, characterized in that: Based on the similarity simulation of seismic waveforms, the method of constructing low-frequency prior information constrained by seismic waves is as follows: segmenting the well curve, wellside seismic data and measured seismic records, and establishing the well sample set, wellside seismic trace sample set and seismic data set to be inverted respectively; Calculate the correlation coefficient between the known sample and the seismic data at the position to be inverted; screen all the logging physical property parameter samples from high to low according to the correlation coefficient; select the logging physical property parameters with the correlation coefficient higher than the threshold; assign different weights to the selected logging physical property parameters, bring them into the Kriging interpolation formula, and obtain the Kriging parameters of the position to be inverted; use the Kriging parameters, combined with sequential Gaussian simulation, to perform random sampling at each of the positions to be inverted to complete a similarity simulation; perform multiple similarity simulations, average the results of multiple similarity simulations, and obtain prior information based on seismic waveform constraints.

3. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 2, characterized in that: The calculation formula of the correlation coefficient is: Among them, r ij is x i The earthquake waveform and x j The similarity coefficient of the earthquake waveform at i ) is x i The variance of the variable at j ) is x j The variance of the variable, C ij is x i Place and x j The spatial covariance value at ; After substituting the correlation coefficient into the Kriging variance formula, the variance within the region is obtained. The formula is: Among them, ω i and ω j are the weight coefficients in the Kriging interpolation formula, σ is the standard deviation of the measured value, and n is the number of sample points; Find the variance within the region The minimum x j Solution: Combined with the constraints of unbiased estimation, construct the objective function, solve the coefficient matrix of the objective function, and obtain the weight coefficient ω based on the seismic waveform similarity coefficient i ; Substitute it into the ordinary Kriging interpolation formula and the Kriging variance formula to obtain the Kriging parameters of sequential Gaussian simulation. According to the random path, calculate the probability density cumulative distribution function of each position to be inverted in turn and perform random sampling to obtain the simulated values ​​of the physical property parameters.

4. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 1, characterized in that: Based on fractal theory, the method for constructing high-frequency prior information of seismic waves is as follows: using a fractional Gaussian distribution model to simulate the frequency range of logging data and construct an initial high-frequency prior information model; using an autocovariance function as a fractional Gaussian distribution to describe a logging curve based on logging data; performing Fourier transform on the autocovariance function to obtain a power spectral density function; generating a random number sequence having the same distribution characteristics as the standardized logging data, convolving the random number sequence with the power spectral density function, and obtaining a high-frequency prior information model of seismic waves.

5. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 4, characterized in that: The high-frequency prior information model of the seismic wave is: Among them, Y(k) is the spectrum value after Fourier transformation, M is the maximum number or range of random number generation, c is a constant, N is the total amount of data, and V(k) is the standardized random number.

6. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 1, characterized in that: The mixed prior information model is: m=αm1+(1-α)m2 Among them, m1 is the low-frequency prior information; m2 is the high-frequency prior information; α is the weight coefficient of the prior initial model obtained based on the seismic waveform similarity simulation; and m is the mixed prior information.

7. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 1, characterized in that: The relationship between the reservoir physical property parameters and the elastic parameters is: d=f(m)+ε=f(g(r))+e Among them, d is the reservoir physical parameter, ε is the error term, m is the elastic parameter, r is the rock physical parameter, g is the rock physics model, f is the seismic forward model, e is the mean, and the covariance is ∑ e Gaussian distribution of data errors.

8. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 1, characterized in that: The method for perturbing and iterating the reservoir physical property parameter inversion model by using the QA-MCMC algorithm is as follows: according to the parameter distribution characteristics at the current moment in the Markov chain process that are only related to the previous moment, a physical property parameter distribution model is established; the relationship between the physical property parameter distribution model at the next moment and the initial parameter distribution model is obtained, which is used as a physical property parameter perturbation, and whether to accept the physical property parameter distribution model at the next moment is judged by the Metropolis sampling criterion, and the distribution model is optimized by the likelihood function to obtain a physical property parameter state model after the MCMC inversion algorithm; the QA algorithm is introduced into the physical property parameter state model to limit the scope of the solution space, and the final reservoir physical property parameter inversion model is obtained.

9. The method for stochastic inversion of reservoir physical property parameters based on mixed prior information according to claim 8, characterized in that: The likelihood function uses low-frequency prior information and high-frequency prior information as constraints, and assigns different weights to different prior information items; the final reservoir physical property parameter inversion model is: in, is the sample point x t and target point The similarity weight between t ) is the sample point x t The probability density function of is the sample point The probability density function of x t is the sample point, is the target point, ex is the residual of the data, and φ(v) is the probability density function.

10. A stochastic inversion system for reservoir physical property parameters based on mixed prior information, characterized in that: include: A low-frequency prior information model building module is used to build a low-frequency prior information model constrained by seismic waves based on similarity simulation of seismic waveforms; A high-frequency prior information model building module is used to build a high-frequency prior information model of seismic waves based on fractal theory; A hybrid prior information model building module, used to combine the low-frequency prior information model and the high-frequency prior information model to build a hybrid prior information model; A physical property elastic parameter relationship module is used to establish a relationship between reservoir physical property parameters and elastic parameters, and combine it with the hybrid prior information model to obtain a reservoir physical property parameter inversion model of seismic data; A QA-MCMC algorithm module is used to perturb and iterate the reservoir physical property parameter inversion model through the QA-MCMC algorithm, and generate a final reservoir physical property parameter inversion model when the minimum square error between the synthetic seismic record and the well bypass seismic record is the smallest; The result output module is used to input the seismic data to be measured into the final reservoir physical property parameter inversion model to obtain the reservoir physical property parameter inversion results.

Images