Prestack earthquake prediction method and device

Abstract

The invention discloses a pre-stack earthquake prediction method and device. The method comprises the following steps: determining reservoir logging sensitive elastic parameters and sedimentary facies plane attributes of a target interval of a research area; according to the reservoir logging sensitive elastic parameters and the sedimentary facies plane attributes, phase control sensitive parameter plane constraint information is established by adopting an in-situ co-Kriging method; establishing a three-dimensional prior constraint model by using a direct sequential simulation method according to the phase control sensitive parameter plane constraint information and the reservoir logging sensitive elastic parameters; and according to the three-dimensional prior constraint model and the pre-stack seismic data of the research area, a Bayesian pre-stack probability inversion method is adopted to predict and obtain pre-stack elastic parameters of a target layer section of the research area.

Description

Technical Field

[0001] The present invention relates to the field of oil and gas exploration technology, and in particular to a method and device for pre-stack seismic prediction. Background Art

[0002] Prestack seismic random inversion is an important means to identify lithology and fluids in complex reservoirs. Prestack random inversion has good fidelity and can obtain multiple parameters such as shear wave impedance, longitudinal and transverse wave velocity ratio, Poisson's ratio, shear modulus, etc., providing richer inversion results for lithology and fluid identification. In the process of complex reservoir prediction, the seismic random inversion method that integrates logging data can broaden the spectrum of the inversion results and improve the resolution of the inversion profile under the constraints of seismic data.

[0003] Conventional seismic random inversion methods have been widely used in the field of seismic reservoir prediction. Commonly used seismic random inversion methods are mainly divided into geostatistical inversion methods and Bayesian probabilistic inversion methods. The geostatistical-based random inversion method uses random simulation methods as a means of model perturbation, uses global optimization algorithms as model optimization strategies, and makes full use of the advantages of logging data in vertical high resolution to improve the vertical resolution of seismic inversion. However, the geostatistical inversion process is relatively complex, with low computational efficiency, and is based on logging data. The prediction effect is affected by the number and distribution of logging data, and the lateral continuity is poor. These problems restrict the industrial application of geostatistical inversion methods.

[0004] In the framework of Bayesian theory, the probabilistic inversion method constructs the prior probability density function of inversion based on the prior constraint model, uses a random sampling algorithm to sample the prior probability density function of the inversion parameters, and matches the samples with the actual seismic records through the likelihood function. Through continuous sample iteration, the effective prediction of the inversion parameters is achieved. Probabilistic inversion relies on the prior constraint model. When the prior constraint model differs greatly from the true model, the reliability of the inversion result is poor. Summary of the invention

[0005] The inventors of the present application have found that the prior constraint model of traditional Bayesian probabilistic inversion is usually obtained by interpolation of well logging data or geostatistical simulation using well logging data as hard data, without considering the impact of geological knowledge on prior information.

[0006] To this end, the present invention proposes a method and device for prestack seismic prediction. The method introduces sedimentary phase plane attribute constraint information into the process of constructing a priori constraint model, combines it with logging data to construct a phase-controlled priori constraint model, and applies it to Bayesian probability inversion, thereby effectively improving the accuracy and rationality of prestack seismic inversion.

[0007] In a first aspect, the present application provides a method for pre-stack earthquake prediction, the method comprising:

[0008] Determine the reservoir logging sensitive elastic parameters and sedimentary facies plane attributes of the target layer in the study area;

[0009] According to the reservoir logging sensitive elastic parameters and the sedimentary phase plane attributes, using the isotopic co-kriging method to establish the phase-controlled sensitive parameter plane constraint information;

[0010] According to the phase-controlled sensitive parameter plane constraint information and the reservoir logging sensitive elastic parameters, a three-dimensional priori constraint model is established by using a direct sequential simulation method;

[0011] According to the three-dimensional prior constraint model and the prestack seismic data of the study area, the prestack elastic parameters of the target layer in the study area are predicted by using the Bayesian prestack probability inversion method.

[0012] In a second aspect, an embodiment of the present invention provides a device for prestack seismic prediction, the device comprising: a memory and a processor; the memory is used to store a program for performing a prestack seismic prediction method, and the processor is used to read and execute the program for performing a prestack seismic prediction method, and execute any one of the methods described in the above embodiments.

[0013] In a third aspect, an embodiment of the present invention provides a computer-readable storage medium, on which a data processing program is stored, and the data processing program is executed by a processor to implement any one of the pre-stack seismic prediction methods.

[0014] Compared with the related art, a method and device for prestack seismic prediction includes: determining the reservoir logging sensitive elastic parameters and sedimentary phase plane attributes of the target layer section in the study area; using the co-positioned co-kriging method to establish the phase-controlled sensitive parameter plane constraint information according to the reservoir logging sensitive elastic parameters and the sedimentary phase plane attributes; using the direct sequential simulation method to establish a three-dimensional priori constraint model according to the phase-controlled sensitive parameter plane constraint information and the reservoir logging sensitive elastic parameters; using the Bayesian prestack probability inversion method to predict the prestack elastic parameters of the target layer section in the study area according to the three-dimensional priori constraint model and the prestack seismic data of the study area. The present application introduces the sedimentary phase plane attribute constraint information in the priori constraint model construction process, combines it with the logging data to construct the phase-controlled priori constraint model, and applies it to the Bayesian probability inversion, thereby effectively improving the accuracy and rationality of prestack seismic inversion.

[0015] Other features and advantages of the present application will be described in the following description, and partly become apparent from the description, or be understood by implementing the present application. Other advantages of the present application can be realized and obtained by the schemes described in the description and the drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] The accompanying drawings are used to provide an understanding of the technical solution of the present application and constitute a part of the specification. Together with the embodiments of the present application, they are used to explain the technical solution of the present application and do not constitute a limitation on the technical solution of the present application.

[0017] Figure 1 is a flow chart of a pre-stack seismic inversion method according to an embodiment of the present application;

[0018] Figure 2 is a schematic diagram of a pre-stack seismic inversion device according to an embodiment of the present application;

[0019] Figure 3 is a flow chart of a prestack seismic inversion method in some exemplary embodiments;

[0020] Figure 4 The invention relates to elastic sensitive parameter cross-section analysis of well logging data in some exemplary embodiments.

[0021] Figure 5 is a plan view of the deposition phase of the K1 layer in some exemplary embodiments.

[0022] Figure 6 a is a schematic diagram of a hard data plane of the K1 layer under a large-scale grid in some exemplary embodiments;

[0023] Figure 6 b is a schematic diagram of discrete attributes of sedimentary facies sensitive parameters on a large-scale grid in some exemplary embodiments;

[0024] Figure 6 c is a schematic diagram of a coarsened plane sensitive parameter constraint attribute plane under a large-scale grid in some exemplary embodiments;

[0025] Figure 7 a is a schematic diagram of a hard data plane of the K1 layer under a small-scale grid in some exemplary embodiments;

[0026] Figure 7 b is a schematic diagram of the plane attributes of the seismic RMS plane under a small-scale grid in some exemplary embodiments;

[0027] Figure 7 c is a schematic diagram of a plane sensitive parameter constraint attribute plane under a small-scale grid in some exemplary embodiments;

[0028] Figure 8 a is a cross-sectional view of the longitudinal wave modulus through the well in some exemplary embodiments;

[0029] Figure 8 b is a cross-sectional view of the longitudinal and transverse wave velocity ratios through the well in some exemplary embodiments;

[0030] Figure 8c is a cross-sectional view of a density three-dimensional priori constraint model through a well in some exemplary embodiments;

[0031] Fig. 9 is a comparative schematic diagram of a facies-controlled prior constraint model plane diagram, a traditional geostatistical prior model plane diagram, and a plane along-layer attribute plane diagram of a well interpolation prior model in some exemplary embodiments;

[0032] Fig.10 is a schematic diagram of the cross-well profile comparison of the longitudinal and transverse wave velocity ratios and longitudinal wave modulus inversion results and the sand body estimation results in some exemplary embodiments;

[0033] Fig.11 It is a schematic diagram comparing the plane attribute diagram of the K1 layer VP / VS phase-controlled Bayesian probabilistic inversion results and the traditional seismic random inversion results in some exemplary embodiments. DETAILED DESCRIPTION

[0034] The present application describes multiple embodiments, but the description is exemplary rather than restrictive, and it is obvious to those skilled in the art that there may be more embodiments and implementations within the scope of the embodiments described in the present application. Although many possible feature combinations are shown in the drawings and discussed in the specific embodiments, many other combinations of the disclosed features are also possible. Unless specifically limited, any feature or element of any embodiment may be used in combination with any other feature or element in any other embodiment, or may replace any other feature or element in any other embodiment.

[0035] The present application includes and contemplates combinations of features and elements known to those of ordinary skill in the art. The embodiments, features and elements disclosed in the present application may also be combined with any conventional features or elements to form a unique invention scheme defined by the claims. Any features or elements of any embodiment may also be combined with features or elements from other invention schemes to form another unique invention scheme defined by the claims. Therefore, it should be understood that any feature shown and / or discussed in the present application may be implemented individually or in any appropriate combination. Therefore, except for the limitations made according to the attached claims and their equivalents, the embodiments are not subject to other restrictions. In addition, various modifications and changes may be made within the scope of protection of the attached claims.

[0036] In addition, when describing representative embodiments, the specification may have presented the method and / or process as a specific sequence of steps. However, to the extent that the method or process does not rely on the specific order of the steps described herein, the method or process should not be limited to the steps of the specific order described. As will be understood by those of ordinary skill in the art, other sequences of steps are also possible. Therefore, the specific sequence of the steps set forth in the specification should not be interpreted as a limitation to the claims. In addition, the claims for the method and / or process should not be limited to the steps of performing them in the order written, and those skilled in the art can easily understand that these sequences can be changed and still remain within the spirit and scope of the embodiments of the present application.

[0037] The embodiment of the present invention provides a method for pre-stack earthquake prediction, such as Figure 1 As shown, the method includes steps S100-S130, which are specifically as follows:

[0038] S100: Determine the reservoir logging sensitive elastic parameters and sedimentary facies plane attributes of the target layer in the study area;

[0039] S110: establishing the phase-controlled sensitive parameter plane constraint information by using the isotopic co-kriging method according to the reservoir logging sensitive elastic parameters and the sedimentary phase plane attributes;

[0040] S120: establishing a three-dimensional priori constraint model using a direct sequential simulation method according to the phase-controlled sensitive parameter plane constraint information and the reservoir logging sensitive elastic parameters;

[0041] S130: According to the three-dimensional prior constraint model and the prestack seismic data of the study area, a Bayesian prestack probability inversion method is used to predict the prestack elastic parameters of the target layer segment in the study area.

[0042] In this embodiment, the plane attributes of the sedimentary phase can be determined by conventional means in the art. Since the seismic sedimentary phase belongs to the concept of a phase, it is difficult to directly apply it to the construction of plane sensitive parameter constraint information. Therefore, through logging data, the distribution characteristics of sensitive parameters of elastic parameters of different sedimentary phases are statistically analyzed, and they are filled into the sedimentary phase map to establish a discrete attribute map of sedimentary phase sensitive parameters.

[0043] In an exemplary embodiment, the reservoir logging sensitive elastic parameters of the target layer section in the study area are determined by performing rock physics analysis on the logging data of the target layer section in the study area; the elastic parameters include: P-wave modulus, P-wave velocity ratio and density variation. For example, the P-wave velocity ratio is determined as a sensitive elastic parameter through rock physics analysis.

[0044] In an exemplary embodiment, a seismic forward modeling relationship is determined.

[0045] The basis of prestack inversion is the Zeoppritz formula and its approximate formula. Through the analysis of the actual work area, it is found that the lithology sensitive elastic parameters are the ratio of longitudinal and transverse wave velocities. In order to reduce the accumulation of errors, the Zeoppritz approximate formula of longitudinal wave velocity M, longitudinal and transverse wave velocity ratio ν and density ρ is derived:

[0046]

[0047] Where θ is the incident angle, and R(θ) represents the PP wave reflection coefficient; They represent the average values ​​of underground P-wave modulus, P-wave velocity ratio and density respectively; ΔM, Δν and Δρ represent the changes of P-wave modulus, P-wave velocity ratio and density on both sides of the stratum interface respectively.

[0048] Given a wavelet W, seismic data can be characterized using a convolution model:

[0049]

[0050] Where d(t,θ) is the seismic amplitude at time t, t is time, θ is the incident angle, R is the reflection sparseness, and τ is the integral variable, ranging from 0 to t.

[0051] Using formula (1) and formula (2), the forward relationship between sensitive parameters and seismic records can be obtained.

[0052] In an exemplary embodiment, based on the reservoir logging sensitive elastic parameters and the sedimentary phase plane attributes, the co-position co-kriging method is used to establish the phase-controlled sensitive parameter plane constraint information, including:

[0053] Establishing large-scale phase-controlled plane sensitive parameter constraint information according to the sedimentary phase plane attributes;

[0054] According to the reservoir logging sensitive elastic parameters and the large-scale phase-controlled plane sensitive parameter constraint information, the phase-controlled plane sensitive parameter constraint information is established by using the co-positioned cokriging method.

[0055] In an exemplary embodiment, establishing large-scale phase-controlled plane sensitive parameter constraint information according to the depositional phase plane attributes includes:

[0056] The first step is to process the original seismic grid data to obtain large-scale seismic grid data;

[0057] In this step, the original seismic grid is a small-scale seismic grid, generally a grid with a line number of 1X1, and the grid is coarsened to obtain a large-scale seismic grid, for example: the line number is 10X10, or it can be 20X20, which can be adjusted according to the specific situation.

[0058] The second step is to determine the plane variogram of the large-scale seismic grid data using the seismic RMS plane attributes of the target layer in the study area;

[0059] In an exemplary embodiment, the theoretical model of the variogram is:

[0060] In the x direction, the variance function γ(x,h) represents half of the variance of the value of the variable Z(x) at the two points x and x+h in space, which can be specifically represented as:

[0061]

[0062] In the above formula, γ(x,h) is the variogram, Z(x) represents the variable value of Z at x, x and x+h are two points in space, and Var represents the variance.

[0063] The experimental variogram is actually a series of discrete points. Therefore, it is necessary to fit the experimental variogram to obtain its specific expression, that is, the theoretical model of the variogram. The spherical model is the most commonly used theoretical model of the variogram. In addition, there are Gaussian model, exponential model and nugget model.

[0064] In this embodiment, the plane variogram of the large-scale seismic grid data is determined by using the seismic root mean square plane attributes of the target layer segment in the study area; the variogram is:

[0065]

[0066] In the formula, γ * (h) represents the variogram value, x i represents the i-th sampling position, x j represents the jth sampling position, N(h) represents the number of sampling point pairs with a distance of h, and h ij It represents the distance between the i-th sampling point and the j-th sampling point, and Z(x) represents the variable value of Z at x.

[0067] According to formula (4), the longitudinal and plane variograms of the model can be calculated based on the well logging P-wave modulus, P-wave velocity ratio and density curve and the RMS plane attributes of the seismic along the layer.

[0068] The third step is to establish the constraint information of large-scale phase-controlled plane sensitive parameters using the co-kriging estimation method based on the determined plane variogram, reservoir logging sensitive elastic parameters, and sedimentary phase plane attributes. The co-kriging method is:

[0069]

[0070]

[0071] In the formula, Z1 represents the discrete attributes of sedimentary facies sensitive parameters, Z 2 Represents the sensitive parameter attributes of the phased plane; Z 2 (x 0 ) * Represents Z 2 In x 0 The kriging mean at i is the kriging weight, Z 2 (x i ) represents x i Z 2 The variable value of Respectively represented in x i Z 2 The Kriging mean of 0 is the cokriging weight, Z 1 (x 0 ) represents x 0 Z 1 The variable value of Characterization x 0 Z 1 The Kriging mean of Characterization x 0 Z 2 The Kriging mean of ; Represents Z 2 In x 0 The variance at Z 2 (x 0 ) represents x 0 Z 2 The variable value of Z 2 (x i ) represents x i Z 2 The variable value of Z can be directly obtained based on the plane properties. 1 (x 0 ) represents x 0 Z 1 The variable value of Z can be directly obtained based on the plane properties. 2 (x 0 ) represents x 0 Z 2 The variable value can be directly obtained based on the plane properties.

[0072] Kriging weight λ i (i=1,2,3,...,N) can be calculated by the following formula:

[0073]

[0074] C(x i ,x j)(i,j=1,2,...,N) represents the point x i At point x j The covariance between the simulated data, C(x i ,x 0 ) represents the point x to be estimated 0 and the adjacent point x i is the covariance between the data to be simulated, and N is the number of adjacent points involved in the estimation.

[0075] In an exemplary embodiment, based on the reservoir logging sensitive elastic parameters and the large-scale phase-controlled plane sensitive parameter constraint information, the phase-controlled plane sensitive parameter constraint information is established using the co-location co-kriging method, including: refining the large-scale phase-controlled plane sensitive parameter constraint information to the original seismic grid through the plane variation function; based on the data phase-controlled plane sensitive parameter constraint information on the original seismic grid and the seismic root mean square plane attribute, the phase-controlled plane sensitive parameter constraint information is established using the co-location co-kriging method. In this embodiment, the elastic parameter estimation values ​​on the large-scale grid points are redistributed to the small grid points and used as hard data, and the plane seismic root mean square amplitude attribute is used as soft data, and the sensitive parameter values ​​on the small-scale seismic grid points are calculated using co-location co-kriging.

[0076] In an exemplary embodiment, according to the plane constraint information of the phase-controlled sensitive parameter and the sensitive elastic parameter of the reservoir logging, a three-dimensional prior constraint model is established by a direct sequential simulation method, and then the correlation coefficient between the sensitive elastic parameter and the other two elastic parameters is determined; according to the sensitive elastic parameter of the reservoir logging, the three-dimensional prior constraint model and the correlation coefficient, a three-dimensional prior constraint model of the other two elastic parameters is established by direct sequential simulation. In this embodiment, prestack inversion requires three elastic parameters, namely, P-wave modulus, P-wave velocity ratio and density. When the P-wave velocity ratio is determined as a sensitive elastic parameter based on rock physics analysis, the correlation coefficient between the P-wave velocity ratio and the other two elastic parameters (i.e., P-wave modulus and density) is determined, and a three-dimensional prior constraint model of the other two elastic parameters is established by direct sequential co-simulation.

[0077] The basis of direct sequential simulation is the simple cokriging method. When some parameter values ​​and covariance functions of the simulation space are known, the following formula can be used to calculate the kriging mean Z(x 0 ) * and Kriging variance Characterization:

[0078]

[0079]

[0080] In the formula, Represents the point x to be estimated0 The mean of the data to be simulated, λ i Characterize adjacent points x i The weight coefficient of the data value to be simulated, Z(x i ) represents the adjacent point x i The variable value of the data to be simulated, For adjacent points x i The mean of the data to be simulated, C(0) is the point x to be estimated 0 The variance of the data to be simulated, C(x i ,x 0 ) represents the point x to be estimated 0 and the adjacent point x i is the covariance between the data to be simulated, and N is the number of adjacent points involved in the estimation.

[0081] In an exemplary embodiment, the three-dimensional prior constraint model corresponding to the three elastic parameters is used as the mean, and the P-wave modulus, P-wave velocity ratio and density curve statistical variance in the logging data are used as the variance to establish the prior probability distribution of the inversion parameters, and sampling is performed from the prior probability distribution to obtain a new three-parameter inversion model. In this embodiment, after obtaining the three-dimensional prior constraint model in x 0 The kriging mean Z(x 0 ) * and the kriging variance σ E (x 0 ), then establish the local sampling space:

[0082]

[0083]

[0084] In the formula, Z i (x 0 ) is the local sampling space N(Z(x 0 ) * ,σ E (x 0 )), where n is the number of sampling values.

[0085] After establishing the local sampling space, define a global cumulative probability density function κ(x 0 ) * is Z(x 0 ) * The corresponding Gaussian value. Generate a random Gaussian number κ G , and the corresponding inversion parameter value is the simulation value Z s (x 0 ).

[0086] In an exemplary embodiment, the synthetic seismic record corresponding to the new three-parameter inversion model is determined through the pre-stack seismic forward modeling relationship; and the likelihood function value of the inversion parameter is calculated based on the synthetic seismic record.

[0087] The inversion can be performed using the Bayesian probability inversion method. Since the inversion has multiple solutions, the main goal of the probability inversion is to obtain the possible values ​​of the inversion parameters, that is, the posterior probability distribution. According to the Bayesian theorem, the posterior probability distribution p(m|d) mainly consists of the prior probability distribution p(m) and the likelihood function p(d|m), which can be characterized as:

[0088] p(m|d)∝p(d|m)·p(m) (11)

[0089] In the above formula, m is the inversion parameter, m = (M, ν, ρ) T , M is the longitudinal wave modulus, V is the ratio of longitudinal and transverse wave velocity, ρ is the density, d is the observed seismic record, p(m) represents the prior probability distribution, p(d|m) represents the likelihood function value, and p(m|d) represents the posterior probability distribution.

[0090] In an exemplary embodiment, the prior probability distribution can be set to a priori constraint model is the mean, which refers to the prior constraint model established in the previous step. Each sample point value of this model is taken as the mean of the prior probability distribution, and the variance of the prior constraint model is taken as the Gaussian distribution with variance;

[0091] Among them, the prior probability distribution is:

[0092]

[0093] Where m is the inversion parameter, m = (M, ν, ρ) T , M is the longitudinal wave modulus, V is the longitudinal and transverse wave velocity ratio, ρ is the density; The variance diagonal matrix of the prior constraint model can be directly obtained based on the prior constraint model; ξ is the number of sample points of the inversion parameters.

[0094] Sampling from the prior probability distribution is done according to the following formula:

[0095]

[0096] In the formula, ζ t is the Gaussian transformation value corresponding to the porosity model value of the current iteration (which can be obtained after Gaussian transformation), δ is a random number between [-1,1], τ is the step size, Characterize the inverse Gaussian transform.

[0097] In an exemplary embodiment, the likelihood function p(d|m) represents the similarity between the seismic record synthesized by the inversion parameters and the actual seismic record, assuming that the seismic noise satisfies the mean of 0 and the variance of Gaussian distribution, it can be characterized as:

[0098]

[0099] In the formula, C e is the variance diagonal matrix of seismic noise (the variance of seismic noise is written in the form of a diagonal matrix), (a quantity that can be directly calculated), ξ represents the unit diagonal matrix of rank N (the unit diagonal matrix is ​​fixed), f represents the forward operator (the forward relationship between the sensitive parameters and the seismic records obtained using formula (1) and formula (2)), and Gaussian distribution (this is based on empirical assumptions).

[0100] Substituting equations (12) and (14) into equation (11), we can obtain the posterior probability distribution expression:

[0101]

[0102] In the above formula, f(m) is the forward function.

[0103] Therefore, the objective function of the inversion can be characterized as:

[0104]

[0105] In the above formula, J(m) is the inversion objective function, m is the inversion parameter, m = (M, ν, ρ) T , M is the longitudinal wave modulus, V is the longitudinal and transverse wave velocity ratio, ρ is the density; The variance diagonal matrix of the prior constraint model can be directly obtained based on the prior constraint model; d is the observed seismic record, f(m) is the forward function; C e is the diagonal matrix of variance of the seismic noise.

[0106] In an exemplary embodiment, the final model of the inversion is determined using the Metropolis criterion;

[0107] After the final model is determined, the objective function of the inversion is established and repeated iterations are performed. When the objective function value is less than the threshold value or the number of iterations reaches the maximum value, a statistical analysis is performed on multiple inversion results to determine the inversion result.

[0108] 6) Metropolis criteria judgment

[0109] The Metropolis criterion is to compare the sampled values ​​to determine whether to keep them. The acceptance probability can be expressed as the following formula:

[0110]

[0111] Among them, J is the objective function of the inversion, k and T are adjustable parameters, and m t is the current model, m * is the new model generated after the disturbance. At this time, a random number between [0,1] is randomly generated, and the random number is compared with the acceptance probability to determine whether to update the model; that is, if the random number is less than or equal to the acceptance probability value, the model is updated, and if the random number is greater than the acceptance probability value, the original value is retained.

[0112] Through a series of iterations, the final inversion result can be obtained by statistically analyzing the inversion parameter values ​​under steady state.

[0113] An embodiment of the present invention provides a device for prestack seismic prediction, the device comprising: a memory and a processor; the memory is used to store a program for performing a prestack seismic prediction method, and the processor is used to read and execute the program for performing a prestack seismic prediction method, and execute any one of the methods in the above embodiments.

[0114] An embodiment of the present invention provides a computer-readable storage medium, on which a data processing program is stored, and the data processing program is executed by a processor to implement any one of the pre-stack seismic prediction methods.

[0115] Example 1

[0116] A pre-stack earthquake prediction method based on phase-controlled Bayesian probabilistic inversion, the main process of which is as follows:

[0117] (1) Optimize the reservoir sensitive elastic parameters of the target layer and establish the forward relationship between the three parameters including the sensitive parameters and the prestack seismic records.

[0118] (2) Based on the planar geological sedimentary facies map and the sensitive elastic parameter values ​​of the well points along the layer, a discrete attribute map of the sedimentary facies sensitive parameters is established.

[0119] (3) The small-scale seismic grid (channel number is 1x1) is coarsened to establish a large-scale seismic grid (channel number is 10x10). The plane variogram of the parameters of large and small scales is calculated using the seismic root mean square attribute plane, and the vertical variogram is calculated using the logging data.

[0120] (4) Based on the large-scale plane variation function, the sensitive elastic parameter values ​​of the well point positions along the layer are used as hard data, and the discrete attribute map of the sedimentary facies sensitive parameters is used as soft data. The co-positioned cokriging method is used to establish the constraint information of the large-scale facies-controlled plane sensitive parameters.

[0121] (5) The large-scale plane sensitive parameter constraint information is assigned to the small-scale seismic grid and used as hard data, and the seismic root mean square attributes are used as soft data. The co-location kriging method is used again to establish the small-scale phase-controlled plane sensitive parameter constraint information.

[0122] (6) Based on the planar and longitudinal variograms, the logging data and planar sensitive parameter constraint information are used as hard data, and a three-dimensional prior constraint model of sensitive parameters is established using direct sequential simulation.

[0123] (7) Taking the logging data as hard data and the three-dimensional prior constraint model of sensitive parameters as soft data, and introducing the correlation coefficient between parameters as weights, direct sequential co-simulation is used to establish the three-dimensional prior constraint model of the other two inversion parameters.

[0124] (8) Taking the three-dimensional prior constraint model as the mean and the statistical variance of the logging data as the variance, the prior probability distribution of the inversion parameters is established, and sampling is performed from the prior probability distribution to obtain a new three-parameter inversion model.

[0125] (9) The synthetic seismic records corresponding to the new model are calculated through the pre-stack seismic forward modeling relationship, and combined with the actual seismic records to calculate the likelihood function value of the inversion parameters.

[0126] (10) The objective function of the inversion is established within the framework of Bayesian theory, and the Metropolis criterion is used to determine whether to retain the new model.

[0127] (11) Repeat steps (8) to (10) until the objective function value is less than the threshold value or the number of iterations reaches the maximum value.

[0128] (12) Perform statistical analysis on multiple inversion results and select the maximum probability under the steady state as the inversion result.

[0129] Example 2

[0130] The actual work area of ​​tight sandstone in northern China was used for verification. The southern part of the work area mainly develops the distributary channel of the river delta and the delta front phase, with thicker sand layers, and the sand bodies in the north are thinner. The sweet spots inside the sand bodies are mainly concentrated inside the high-energy thick-layer river channel sand bodies, which have the characteristics of low elastic modulus and relatively high porosity and permeability. The details are as follows:

[0131] Step 1: Determine the reservoir logging sensitive elastic parameters of the target layer in the study area

[0132] According to the logging data of the target layer in the study area, the rock physical analysis is carried out to determine the reservoir logging sensitive elastic parameters of the target layer. For example, by conducting rock physical analysis on multiple wells in the study area, such as Figure 4As shown, it can be found that the longitudinal and transverse wave velocity ratio is more sensitive to the sand body, and the longitudinal wave modulus is more sensitive to the sweet spot inside the sand body.

[0133] Step 2, establishing sedimentary phase plane attributes;

[0134] Through logging data, the distribution characteristics of different sedimentary facies sensitive parameters are statistically analyzed and filled into the sedimentary facies map to establish a discrete attribute plane map of sedimentary facies sensitive parameters, namely, sedimentary facies plane attributes. Taking the K1 layer as an example, Figure 5 Represents the sedimentary phase plane of the K1 layer.

[0135] Step 3: Based on the reservoir logging sensitive elastic parameters and sedimentary phase plane attributes, the isotopic co-kriging method is used to establish the phase-controlled sensitive parameter plane constraint information. The specific implementation process includes:

[0136] Step 31. Process the small-scale seismic grid data to obtain large-scale seismic grid data;

[0137] The small-scale seismic grid is generally the original seismic data grid, that is, the line number 1X1. The small-scale seismic grid is coarsened to establish a large-scale seismic grid. For example, the large-scale seismic grid is the line number 10X10. The coarsening process can reduce the number of points to be estimated.

[0138] For example: Figure 6 a is a schematic diagram of the hard data plane of the K1 layer under the large-scale grid; Figure 6 b is a schematic diagram of the discrete attribute plane of sedimentary facies sensitive parameters under large-scale grid; Figure 6 c is a schematic diagram of the coarsened plane sensitive parameter constraint attribute plane under large-scale grid.

[0139] Step 32. Determine the horizontal variogram based on the root mean square attributes and the planar distribution law of the sedimentary facies, and determine the vertical variogram based on the logging information. On this basis, carry out geostatistical random inversion.

[0140] The discrete sedimentary facies values ​​of the well point locations in the study area are assigned to the nearest grid points to obtain the large-scale plane variogram values.

[0141] Wherein, the variogram is:

[0142]

[0143] In the formula, γ * (h) represents the large-scale plane variation function value, x i , x j They represent different sampling positions respectively, and N(h) is the number of sampling point pairs with a distance of h.

[0144] Step 33. Based on the large-scale plane variation function, the reservoir logging sensitive elastic parameters are used as hard data, and the discrete attributes of the sedimentary phase sensitive parameters are used as soft data. The large-scale phase-controlled plane sensitive parameter constraint information is established by using the isotopic cokriging estimation method;

[0145] The theoretical model of the co-kriging method is as follows:

[0146]

[0147]

[0148] In the formula, Z 1 represents the discrete attributes of sedimentary facies sensitive parameters, Z 2 Represents the large-scale elastic sensitive parameter properties; Z 2 (x 0 ) * Represents Z 2 At the point x to be estimated 0 The Kriging mean at , Represents Z 2 At the point x to be estimated 0 The variance at ; Characterization x 0 Z 2 The mean of Respectively represented in x i Z 2 The mean of Characterization x 0 Z 1 The mean value, Z 2 (x i ) represents x i Z 2 With Z 1 The variable value, Z 1 (x 0 ) represent x 0 Z 2 With Z 1 The variable value, λ i is the Kriging weight, η 0 is the cokriging weight, Z 1 (x 0 ) represents the point x to be estimated 0 The known variable Z 1 For example: Figure 7 a is a schematic diagram of the hard data plane of the K1 layer under a small-scale grid; Figure 7 b is a schematic diagram of the plane attributes of the earthquake RMS plane under the small-scale grid; Figure 7 c is a schematic diagram of the plane sensitive parameter constraint attribute under small-scale grid.

[0149] Step 34. Assign the large-scale phase-controlled plane sensitive parameter constraint information to the small-scale seismic grid and use it as hard data, and the seismic root mean square attributes as soft data, and use the co-location kriging method to establish the small-scale phase-controlled plane sensitive parameter constraint information.

[0150] Step 4: establishing a three-dimensional priori constraint model using a direct sequential simulation method according to the phase-controlled sensitive parameter plane constraint information and the reservoir logging sensitive elastic parameters; Figure 8 a is a cross-sectional view of the longitudinal wave modulus through the well in some exemplary embodiments; Figure 8 b is a cross-sectional view of the longitudinal and transverse wave velocity ratios through the well in some exemplary embodiments; Figure 8 c is a cross-sectional view of a density three-dimensional prior constraint model through a well in some exemplary embodiments.

[0151] Step 5: Use the logging data as hard data, the three-dimensional prior constraint model of sensitive parameters as soft data, and introduce the correlation coefficient between parameters as weights. Use direct sequential simulation to establish the three-dimensional prior constraint model of the other two inversion parameters to obtain a new three-parameter inversion model. Fig. 9 is a comparative schematic diagram of a facies-controlled prior constraint model plane diagram, a traditional geostatistical prior model plane diagram, and a plane along-layer attribute plane diagram of a well interpolation prior model in some exemplary embodiments;

[0152] Step 6: Determine the synthetic seismic record corresponding to the new three-parameter inversion model through the pre-stack seismic forward modeling relationship; and calculate the likelihood function value of the inversion parameter based on the synthetic seismic record.

[0153] Step 7: Bayesian Probability Inversion

[0154] Since inversion has multiple solutions, the main goal of probabilistic inversion is to obtain the possible values ​​of the inversion parameters, that is, the posterior probability distribution.

[0155] Step 8: Prestack elastic parameter inversion, the process includes:

[0156] Step 81. Determine the final inversion model by using the Metropolis criterion;

[0157] The Metropolis criterion is to compare the sampled values ​​to determine whether to keep them. The acceptance probability can be expressed as the following formula:

[0158]

[0159] Among them, J is the objective function of the inversion, k and T are adjustable parameters, and m t is the current model, m *is the new model generated after the disturbance. At this time, a random number between [0,1] is randomly generated, and the random number is compared with the acceptance probability to determine whether to update the model.

[0160] Through a series of iterations, the final inversion result can be obtained by statistically analyzing the inversion parameter values ​​under steady state.

[0161] Step 82. After the final model is determined, the objective function of the inversion is established and repeated iterations are performed. When the objective function value is less than the threshold value or the number of iterations reaches the maximum value, a statistical analysis is performed on multiple inversion results to determine the inversion result.

[0162] Using the above process, the well section of the prior constraint model established based on the logging data and the phase control plane attribute constraint information is Figure 8 It can be found that the model has a high resolution and is consistent with the logging data. Therefore, this model is used as a priori constraint model to directly estimate the P-wave modulus and P-wave velocity ratio through Bayesian probability inversion. Fig. 9 In the figure, the K1 layer facies-controlled prior constraint model (left) is compared with the traditional geostatistical prior model (middle) and the well interpolation prior model (right) in terms of plane layer attributes. It can be found that the facies-controlled prior constraint model is consistent with geological knowledge and has better rationality. Fig.10 It is a cross-well profile diagram of the inversion results of the P-wave velocity ratio and the P-wave modulus and the sand body estimation results. It can be found that the inversion results of sensitive parameters are consistent with the logging data, and the sand body results are well matched with the logging interpretation results. Fig.11 The figure shows the comparison of the plane attribute map of the VP / VS phase-controlled Bayesian probabilistic inversion results of the K1 layer (left) and the traditional seismic random inversion results (right). The figure shows the comparison of the average attribute maps along the layer of the VP / VS inversion results extracted by different inversion methods in the K1 layer. It can be found that the prediction results of the phase-controlled Bayesian probabilistic inversion are consistent with the geological understanding, which verifies the rationality of the inversion method. The traditional Bayesian probabilistic inversion results have certain differences from the geological understanding. Therefore, the Bayesian probabilistic inversion method based on the phase-controlled prior constraint model has certain advantages.

[0163] It will be appreciated by those skilled in the art that all or some of the steps, systems, and functional modules / units in the methods disclosed above may be implemented as software, firmware, hardware, and appropriate combinations thereof. In hardware implementations, the division between the functional modules / units mentioned in the above description does not necessarily correspond to the division of physical components; for example, a physical component may have multiple functions, or a function or step may be performed by several physical components in cooperation. Some or all components may be implemented as software executed by a processor, such as a digital signal processor or a microprocessor, or implemented as hardware, or implemented as an integrated circuit, such as an application-specific integrated circuit. Such software may be distributed on a computer-readable medium, which may include a computer storage medium (or non-transitory medium) and a communication medium (or temporary medium). As known to those skilled in the art, the term computer storage medium includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storing information (such as computer-readable instructions, data structures, program modules, or other data). Computer storage media include, but are not limited to, RAM, ROM, EEPROM, flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tapes, magnetic disk storage or other magnetic storage devices, or any other medium that can be used to store the desired information and can be accessed by a computer. In addition, it is well known to those of ordinary skill in the art that communication media typically contain computer-readable instructions, data structures, program modules, or other data in a modulated data signal such as a carrier wave or other transport mechanism, and may include any information delivery media.

Claims

1. A method for pre-stack earthquake prediction, characterized in that: The method comprises: Determine the reservoir logging sensitive elastic parameters and sedimentary facies plane attributes of the target layer in the study area; According to the reservoir logging sensitive elastic parameters and the sedimentary phase plane attributes, using the isotopic co-kriging method to establish the phase-controlled sensitive parameter plane constraint information; According to the phase-controlled sensitive parameter plane constraint information and the reservoir logging sensitive elastic parameters, a three-dimensional priori constraint model is established by using a direct sequential simulation method; According to the three-dimensional prior constraint model and the prestack seismic data of the study area, the prestack elastic parameters of the target layer in the study area are predicted by using the Bayesian prestack probability inversion method.

2. The method for pre-stack earthquake prediction according to claim 1, characterized in that: The reservoir logging sensitive elastic parameters of the target layer section in the study area are determined by performing rock physical analysis on the logging data of the target layer section in the study area; The elastic parameters include: longitudinal wave modulus, longitudinal and transverse wave velocity ratio and density variation.

3. The method for pre-stack earthquake prediction according to claim 1, characterized in that: The method of establishing the phase-controlled sensitive parameter plane constraint information by using the co-location co-kriging method according to the reservoir logging sensitive elastic parameters and the sedimentary phase plane attributes includes: Establishing large-scale phase-controlled plane sensitive parameter constraint information according to the sedimentary phase plane attributes; According to the reservoir logging sensitive elastic parameters and the large-scale phase-controlled plane sensitive parameter constraint information, the phase-controlled plane sensitive parameter constraint information is established by using the co-positioned cokriging method.

4. The method for pre-stack earthquake prediction according to claim 3, characterized in that: The step of establishing large-scale phase-controlled plane sensitive parameter constraint information according to the depositional phase plane attributes includes: Processing the original seismic grid data to obtain large-scale seismic grid data; The plane variogram of large-scale seismic grid data is determined using the seismic RMS plane attributes of the target layer in the study area. Based on the determined plane variogram, according to the reservoir logging sensitive elastic parameters and sedimentary phase plane attributes, the co-positioned cokriging method is used to establish large-scale phase-controlled plane sensitive parameter constraint information.

5. The method for pre-stack earthquake prediction according to claim 4, characterized in that: The variogram is: In the formula, γ * (h) represents the variogram value, x i represents the i-th sampling position, x j represents the jth sampling position, N(h) represents the number of sampling point pairs with a distance of h, and h ij Represents the distance between the i-th sampling point and the j-th sampling point.

6. The method for pre-stack earthquake prediction according to claim 5, characterized in that: The method of establishing the phase-controlled plane sensitive parameter constraint information by using the co-location cokriging method according to the reservoir logging sensitive elastic parameters and the large-scale phase-controlled plane sensitive parameter constraint information includes: The large-scale phase-controlled plane sensitive parameter constraint information is combined with the original seismic grid through a plane variation function to obtain the data phase-controlled plane sensitive parameter constraint information on the original seismic grid; The phase control plane sensitive parameter constraint information is established based on the data phase control plane sensitive parameter constraint information on the original seismic grid and the seismic root mean square plane attributes by using the co-location kriging method.

7. The method for pre-stack earthquake prediction according to claim 6, characterized in that: The isotopic cokriging method is: Where Z1 represents the discrete attribute of the sedimentary phase sensitive parameter, Z2 represents the attribute of the phase-controlled plane sensitive parameter; Z2(x0) * represents the Kriging mean of Z2 at x0, λ i is the Kriging weight, Z2(x i ) represents x i The variable value of Z2 at Respectively represented in x i The kriging mean of Z2 at ; η0 is the cokriging weight, Z1(x0) represents the variable value of Z1 at x0, represents the kriging mean of Z1 at x0, Characterizes the Kriging mean of Z2 at x0; Represents the variance of Z2 at x0; Z2(x0) represents the variable value of Z2 at x0.

8. The method for pre-stack earthquake prediction according to claim 2, characterized in that: After establishing the three-dimensional priori constraint model by using the direct sequential simulation method according to the phase-controlled sensitive parameter plane constraint information and the reservoir logging sensitive elastic parameters, the method further includes: Determine the correlation coefficient between the sensitive elastic parameter and the other two elastic parameters in the elastic parameters; According to the reservoir logging sensitive elastic parameters, the three-dimensional priori constraint model and the correlation coefficient, direct sequential co-simulation is used to establish the three-dimensional priori constraint models of the remaining two elastic parameters.

9. The method for pre-stack earthquake prediction according to claim 8, characterized in that: The direct sequential simulation method is: Use the following formula to calculate the Kriging mean Z(x0): * and Kriging variance Characterization: Where Z(x0) * represents the Kriging mean, represents the mean of the simulated data at x0, λ i Characterize adjacent points x i The weight coefficient of the data value to be simulated, Z(x i ) represents the adjacent point x i The variable value of the data to be simulated, For adjacent points x i is the mean of the data to be simulated, is the Kriging variance, C(0) is the variance of the simulated data at x0, C(x i ,x0) represents the distance between the estimated point x0 and the adjacent point x i is the covariance between the data to be simulated, and N is the total number of adjacent points involved in the estimation.

10. The method for pre-stack earthquake prediction according to claim 8, characterized in that: The method of predicting the prestack elastic parameters of the target layer in the study area by using the Bayesian prestack probability inversion method according to the three-dimensional prior constraint model and the prestack seismic data of the study area includes: The three-dimensional prior constraint model corresponding to the three elastic parameters established is used as the mean of the prior probability distribution; The prior probability distribution of inversion parameters is established by taking the logging P-wave modulus, P-wave velocity ratio and density curve statistical variance as the variance; Sampling from the prior probability distribution to obtain a new three-parameter inversion model; The prestack elastic parameters of the target layer in the study area are predicted based on the new three-parameter inversion model.

11. The method for pre-stack earthquake prediction according to claim 10, characterized in that: After obtaining the new three-parameter inversion model, the method further includes: Through the pre-stack seismic forward modeling relationship, the synthetic seismic record corresponding to the new three-parameter inversion model is determined; Based on the synthetic seismic record, the likelihood function value of the inversion parameter is calculated.

12. The method for pre-stack earthquake prediction according to claim 11, characterized in that: After calculating the likelihood function value of the inversion parameter, the method further includes: Use the Metropolis criterion to determine the final inversion model; After the final model is determined, the objective function of the inversion is established and repeated iterations are performed. When the objective function value is less than the threshold value or the number of iterations reaches the maximum value, a statistical analysis is performed on multiple inversion results to determine the inversion result.

13. A device for pre-stack earthquake prediction, characterized in that: The device comprises: a memory and a processor; the memory is used to store a program for performing a prestack seismic prediction method, and the processor is used to read and execute the program for performing a prestack seismic prediction method, and execute the method described in any one of claims 1-12.

14. A computer-readable storage medium having a data processing program stored thereon, wherein the data processing program is executed by a processor to implement the pre-stack seismic prediction method according to any one of claims 1 to 12.

Images