Layer-constrained multi-gaussian fast prestack stochastic inversion method
The multivariate Gaussian fast pre-stack stochastic inversion method with stratigraphic constraints solves the problem of large computational cost of stochastic inversion, improves computational efficiency and inversion results, and is particularly suitable for high-resolution reservoir prediction of complex oil and gas reservoirs.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- NORTHEAST GASOLINEEUM UNIV
- Filing Date
- 2023-01-29
- Publication Date
- 2026-04-28
AI Technical Summary
Existing stochastic inversion methods are computationally intensive and difficult to solve effectively in practical applications. In particular, the inversion process of large kernel matrices in Bayesian linearized stochastic inversion is difficult, which affects the efficiency and effectiveness of high-resolution reservoir prediction.
A multivariate Gaussian fast pre-stack stochastic inversion method with horizon constraints is adopted. By introducing seismic interpretation horizon information to improve the covariance matrix and constructing a complex horizon control operator, and combining the independence assumption of well logging data and seismic data, large matrix inversion is avoided. Bayesian linearized stochastic inversion is used to obtain high-resolution model parameter distribution.
It improves computational efficiency by more than 30%, reduces computer memory requirements, effectively characterizes complex geological structures, and improves the resolution and accuracy of inversion results, making it suitable for the development of remaining oil in complex oil and gas reservoirs.
Smart Images
Figure CN115951410B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of petroleum exploration, specifically relating to a stratigraphically constrained multivariate Gaussian fast pre-stack stochastic inversion method. Background Technology
[0002] my country's demand for oil and gas resources has shown a significant upward trend with economic development. Seismic exploration is one of the most critical technologies in the field of oil and gas exploration, and it has profound significance for increasing the exploration and development efforts of oil and gas resources. As the exploitation level of oil and gas fields in my country increases, the oil and gas resources provided by traditional structural oil and gas reservoirs are nearing depletion. Exploration targets are becoming increasingly complex, the quality of remaining resources is gradually deteriorating, and the difficulty of exploration is constantly increasing, while the production of oil and natural gas resources has also been gradually declining recently. Currently, most oil and gas reservoirs in China have entered a high water-cut exploration stage, making oil extraction in various oil fields quite difficult. Therefore, the exploration of remaining oil has gradually become a focus of exploration for major oil fields. The key to improving the recovery rate of remaining oil lies in high-resolution reservoir prediction technology, which is also one of the technologies most concerned by major oil fields in China.
[0003] High-resolution seismic inversion is a key technology for reservoir prediction and enhanced oil recovery. Seismic inversion integrates and utilizes multi-source information from geology, seismic exploration, and geophysical logging to estimate the distribution of various subsurface parameters. The goal of seismic inversion is to estimate subsurface parameters using surface observation data, thereby providing a basis for reservoir or mineral development. Generally, seismic inversion methods can be divided into two main categories: deterministic inversion methods and stochastic inversion methods. In seismic exploration, deterministic inversion methods estimate the optimal solutions for subsurface elastic or physical parameters using optimization algorithms. To date, there has been much research on deterministic inversion methods. Since most geophysical inverse problems are usually ill-posed, they are not only difficult to solve but also have strong multiple solutions. Therefore, in inversion, exploring the solution space of model parameters is more reasonable than solving for the optimal solution. To achieve this, stochastic inversion was proposed and has been continuously improved in previous studies. Currently, stochastic inversion can comprehensively utilize multi-scale information such as seismic and well logging data to estimate the posterior probability distribution of subsurface parameters. Studies have shown that stochastic inversion results have higher resolution than conventional deterministic inversion.
[0004] Currently, stochastic inversion can be divided into two main categories: iterative geostatistical stochastic inversion and linear Bayesian stochastic inversion. Iterative geostatistical stochastic inversion uses geostatistical stochastic simulation methods to update model parameters through perturbation, obtaining multiple equally probable stochastic solutions under the constraint of seismic data. Although iterative geostatistical stochastic inversion has been extensively studied, the large amount of seismic forward modeling involved makes its computational cost very high, limiting its practical application. Compared to iterative methods, Bayesian linearized stochastic inversion requires less iterative computation and is more efficient. This method assumes that both model parameters and observational data (such as seismic or well data) follow a Gaussian or multivariate Gaussian distribution. Based on linear Gaussian theory, the posterior probability distribution of the model parameters can be directly estimated. This distribution includes a posterior mean and a posterior covariance matrix, representing the smoothed estimate and the uncertainty estimate of the model parameters, respectively. Using stochastic simulation techniques, a large number of equally probable stochastic solutions can be obtained from this posterior probability distribution. However, a key and unavoidable problem in traditional Bayesian linearized stochastic inversion is the inversion of a large kernel matrix. In practical applications, the dimension of this kernel matrix may exceed ten thousand or even one hundred thousand, making its inversion process often difficult to implement. Therefore, the application of Bayesian linearized stochastic inversion poses a significant challenge to both computer performance and matrix inversion algorithms.
[0005] Both of the aforementioned stochastic inversion methods have their own drawbacks. Iterative geostatistical inversion methods are limited by countless seismic forward modeling calculations, while Bayesian linearized stochastic inversion is limited by the inversion of large kernel matrices. Therefore, high computational cost seems to be an unavoidable problem in stochastic inversion. Although stochastic inversion is the most effective high-resolution reservoir characterization technique, its high computational cost severely hinders its practical application. Most current research fails to effectively address the issue of the large computational cost of stochastic inversion. Summary of the Invention
[0006] To address the problems described in the background art, this invention provides a layer-constrained multivariate Gaussian fast pre-stack stochastic inversion method to solve the problems of low computational efficiency and low ability to characterize complex geological bodies in existing seismic stochastic inversion techniques.
[0007] The technical solution adopted in this invention is: a layer-constrained multivariate Gaussian fast pre-stack stochastic inversion method, characterized in that: the layer-constrained multivariate Gaussian fast pre-stack stochastic inversion method includes the following steps:
[0008] Step 1: Determine the horizontal and vertical lateral line ranges and the longitudinal time range of the inversion data; select comprehensive and high-quality logging data as experimental wells for inversion, select inversion quality control wells, and perform well-seismic calibration and wavelet extraction; determine the horizontal and vertical lateral line ranges and the longitudinal time range of the inversion data; select experimental wells for inversion, select inversion control wells, and perform well-seismic calibration and wavelet extraction.
[0009] The selected data range should include the target area of the oilfield exploration focus in the plane, and in the time direction, it should include the target layer and the interval of about 50 milliseconds above and below the target layer. A time range that is too large will increase the computational load of the inversion, while a range that is too small will make it difficult to eliminate the sidelobe effect of the wavelet during inversion. Regarding the selection of well data, firstly, well locations that have passed through the exploration target should be selected in the work area map; secondly, the well data should penetrate the target layer vertically; thirdly, the well data should include parameters such as P-wave transit time, density, and S-wave transit time. If S-wave data is unavailable, S-wave prediction is required; finally, while ensuring the quality of the well data, as many wells as possible should be used to ensure high resolution of the inversion results; one to two wells should be selected from the wells within the target area as quality control wells, which will not participate in the inversion but will only be used to verify the inversion effect.
[0010] Step 2: Based on the single-channel inversion strategy (to avoid large matrix inversion), use Bayesian linearized stochastic inversion to obtain the posterior probability distribution of P-wave and S-wave velocity and density model parameters under seismic data constraints.
[0011] Assume the model parameters m (P-wave velocity, S-wave velocity, density) follow a multivariate Gaussian distribution:
[0012] m~N(μ m , ∑ m (1)
[0013] In the inversion v p v s ρ represents the P-wave velocity, S-wave velocity, and density of channel 1, respectively, and T is the transpose of the matrix; μ m and ∑ m Let m be the mean of the posterior probability distribution and the multivariate covariance matrix.
[0014] In the traditional method, ∑ m It can be obtained from the variation function of the model parameters. This matrix also measures the correlation between different sampling points based on the distance between them. However, in practice, due to the complexity of geological structures, the correlation between two closely spaced sampling points may not be high. It is necessary to introduce information such as the stratigraphic position of seismic interpretation to constrain the correlation.
[0015] By interpreting horizon information through seismic data, the covariance matrix is improved, and a covariance matrix with complex horizon constraints is proposed. When calculating the improved matrix, the horizon information is used as a reference surface to calculate the "true" distance between different sampling points under horizon constraints, and then the covariance function between different sampling points is calculated. Let the improved matrix be ∑ hor Then, the posterior probability distribution of the model parameters under the hierarchical constraints is:
[0016] m~N(μ m , ∑ hor (2)
[0017] The seismic forward modeling expression based on the convolution model is:
[0018] d = Gm + e, (3)
[0019] Where G is the pre-stack forward modeling operator, which can be obtained through the Zoeppritz equation or its approximation, and e is the error term, which is assumed to follow a Gaussian distribution with zero mean:
[0020] e~N(0,∑ e (4)
[0021] In the formula, ∑e is the covariance matrix of e. According to formulas (2) and (3), the probability distribution of the earthquake record can be obtained:
[0022] d~N(Gμ m , ∑ d (5)
[0023] In the formula, Gμ m Let ∑ be the mean of the earthquake records. d =G∑ hor G T +∑ e Let be the covariance matrix of the seismic signal. Based on multivariate linear Gaussian theory, given the seismic record d, the conditional probability distribution of the model parameter m is still a Gaussian distribution, which can be expressed as:
[0024] m|d~N(μ m|d , ∑ m|d (6)
[0025] In the formula, the posterior mean μ m|d and the posterior multivariate covariance matrix ∑ m|d It can be represented as:
[0026]
[0027] Where, ∑ md =∑ hor G TThe covariance matrix between seismic data and model parameters under the constraint of horizon information. This refers to information related to the model parameter m derived from seismic data. I is defined as the degree to which the uncertainty of seismic data with respect to m is reduced. d and V d For earthquake information; in formula (7), μ m|d Considered as the result of deterministic inversion, ∑ m|d These two parameters are used to characterize the strength of the uncertainty in the inversion result. Based on these two parameters, a stochastic inversion solution is obtained by using LU decomposition simulation or Gibbs simulation.
[0028] Step 3: Based on the stratigraphic information interpreted from the seismic data, construct a complex stratigraphic control operator. Based on geostatistics, under the constraints of the stratigraphic control operator, establish a forward modeling expression between the model parameters and well logging data, and finally obtain the posterior probability distribution of the model parameters under the constraints of well logging data.
[0029] In the inversion, m contains the values of multiple model parameters, and this multiple data includes well-side channels; then the model parameters of the well-side channels can be represented as Hm, where the operator H is a sparse matrix, which is used to select the well-side channels from m. Therefore, the well-seismic relationship in traditional stochastic inversion can be expressed as:
[0030] h = Hm, (8)
[0031] Where h represents well logging data, which is sparsely distributed; however, traditional Bayesian linearization inversion does not consider geological structure. When the geological structure is complex, it is unreasonable to represent the relationship between well logging data and model parameters solely through the linear operator H. Here, seismic interpretation stratigraphic information is introduced to improve the traditional operator, making the relationship between the improved model parameters and well logging data as follows:
[0032] h = H hor m, (9)
[0033] Wherein, matrix H hor Constrained by seismic interpretation results, it is defined as a complex horizon control operator, whose number of rows and columns is controlled by the horizon. When the horizon information is complex, H is no longer a square matrix.
[0034] Based on the above formula, the probability distribution of well logging data is as follows:
[0035]
[0036] Step 4: Based on the assumption that well logging data and seismic data are statistically independent, the posterior probability distributions of model parameters constrained by seismic and well logging data are fused to obtain the probability distribution of model parameters constrained by both well logging and seismic data. Different inversion solutions are obtained through stochastic simulation.
[0037] Based on linear Gaussian theory, under the joint constraints of seismic data d and well logging data h, the conditional probability distribution of the model parameter m is still a Gaussian distribution:
[0038] m|(h,d)~N(μ) m|(h，d) , ∑ m|(h，d) (11)
[0039] Among them, the posterior mean μ of the model parameters under simultaneous well-seismic constraints m|(h，d) and the posterior multivariate covariance matrix ∑ m|(h，d) Represented as:
[0040]
[0041] Where, ∑ hh Let ∑ be the covariance matrix of the well data. hd It is a covariance matrix representing the correlation between well logging data and seismic data, ∑ hh and ∑ hd The specific expression is:
[0042]
[0043] Where, μ m|(h，d) Considered as the result of deterministic inversion under well seismic constraints, ∑ m|(h，d) This method is used to characterize the strength of uncertainty in the inversion result; through stochastic simulation, stochastic inversion solutions with different and equally probable well-seismic simultaneous constraints can be obtained; however, this method has a serious drawback, namely, when solving for the posterior probability distribution, the kernel matrix C must be considered. hd The inversion process, which cannot be avoided through approximation in traditional linear Bayesian stochastic inversion methods, involves C. hd The seismic forward modeling operator G is complex and difficult to invert, and its inversion process is both difficult and time-consuming, which greatly limits the practical application of the traditional BLI method.
[0044] Fast stochastic inversion is based on the assumption that well logging data and seismic data are independent, ignoring the statistical correlation between them. This assumption is based on several reasons: First, seismic and well logging data are two different types of data, acquired using different methods and operating in different frequency bands, resulting in low statistical correlation. Second, seismic data reflects changes in subsurface reflection coefficients, i.e., relative changes in elastic parameters such as velocity; while well logging data represents absolute values of elastic parameters, leading to significant physical differences between the two. Based on this assumption, it can be concluded that:
[0045] ∑ hd =0 (14)
[0046] Therefore, the kernel matrix in formula (12) can be transformed into ：
[0047]
[0048] The above formula can be further transformed into:
[0049]
[0050] Further derivation shows that the above formula can be transformed into the following form:
[0051]
[0052] The above expression can be represented in a concise form as follows:
[0053]
[0054] Among them I b It is information related to the absolute value of the model parameter m from well logging data, V b I represents the degree to which the uncertainty of well logging data is reduced for m; b and V b For well logging information, its expression is:
[0055]
[0056] According to formula (19), by inputting seismic data, well logging data, and stratigraphic data, high-resolution well-seismic joint inversion results under stratigraphic information constraints can be obtained.
[0057] The beneficial effects of this invention are as follows: By introducing the assumption of statistical independence of well and seismic information in Bayesian stochastic inversion, the inversion of large matrices is avoided, greatly improving the computational efficiency of stochastic inversion and reducing the computer memory requirements. Actual data testing verifies that, while maintaining inversion quality, the proposed method improves computational efficiency by more than 30% compared to traditional Bayesian stochastic inversion methods, while maintaining consistent inversion results. Furthermore, the introduction of a stratigraphic control operator effectively addresses the high requirement for geological body stability in traditional geostatistical stochastic inversion, which is beneficial for characterizing complex geological structures and is of great significance for the development of remaining oil in complex oil and gas reservoirs. Attached Figure Description
[0058] Figure 1 This is a diagram of pre-stack gather data with an incident angle of 5 degrees in Example 1;
[0059] Figure 2 This is a pre-stack gather data diagram with an incident angle of 15 degrees in Example 1;
[0060] Figure 3 This is a pre-stack gather data diagram with an incident angle of 25 degrees in Example 1;
[0061] Figure 4 This is a diagram showing the Bayesian linearization inversion result of the longitudinal wave velocity in Example 1;
[0062] Figure 5 This is a diagram showing the Bayesian linearization inversion result of the shear wave velocity in Example 1;
[0063] Figure 6 This is a diagram showing the Bayesian linearization inversion result of the density in Example 1;
[0064] Figure 7 To test the target horizon shown in the seismic record, and the data from two wells passing through the horizon that participated in the inversion, a geological structure model was established for the inversion.
[0065] Figure 8 The result of the stochastic inversion of P-wave velocity is obtained by efficient multivariate Gaussian well seismic inversion based on layer constraints in Example 1.
[0066] Figure 9 The result is the stochastic inversion of shear wave velocity obtained by efficient multivariate Gaussian well seismic inversion based on stratigraphic constraints in Example 1.
[0067] Figure 10 This is a diagram showing the density stochastic inversion results obtained from the efficient multivariate Gaussian well seismic stochastic inversion based on stratigraphic constraints in Example 1.
[0068] Figure 11 This is a comparison chart of the inversion results of P-wave and S-wave velocities and densities at well W2 in Example 1 with the well curve, where the dashed line represents the well curve and the solid line represents the inversion results.
[0069] Figure 12 This is a histogram comparison of the P-wave velocity inversion result (left) and the P-wave velocity well curve (right) at well W2 in Example 1.
[0070] Figure 13 This is a histogram comparison of the shear wave velocity inversion result (left) and the shear wave velocity well curve (right) at well W2 in Example 1.
[0071] Figure 14 This is a histogram comparison of the density inversion results (left) and the density well curve (right) at the verification well W2 in Example 1. Detailed Implementation
[0072] Example 1
[0073] Referring to the figures,
[0074] A hierarchical-constrained multivariate Gaussian fast pre-stack stochastic inversion method, characterized by the following steps:
[0075] Step 1: Determine the horizontal and vertical lateral line ranges and the longitudinal time range of the inversion data; select comprehensive and high-quality logging data as experimental wells for inversion, select inversion quality control wells, and perform well-seismic calibration and wavelet extraction; determine the horizontal and vertical lateral line ranges and the longitudinal time range of the inversion data; select experimental wells for inversion, select inversion control wells, and perform well-seismic calibration and wavelet extraction.
[0076] The selected data range should include the target area of the oilfield exploration focus in the plane, and in the time direction, it should include the target layer and the interval of about 50 milliseconds above and below the target layer. A time range that is too large will increase the computational load of the inversion, while a range that is too small will make it difficult to eliminate the sidelobe effect of the wavelet during inversion. Regarding the selection of well data, firstly, well locations that have passed through the exploration target should be selected in the work area map; secondly, the well data should penetrate the target layer vertically; thirdly, the well data should include parameters such as P-wave transit time, density, and S-wave transit time. If S-wave data is unavailable, S-wave prediction is required; finally, while ensuring the quality of the well data, as many wells as possible should be used to ensure high resolution of the inversion results; one to two wells should be selected from the wells within the target area as quality control wells, which will not participate in the inversion but will only be used to verify the inversion effect.
[0077] Step 2: Based on the single-channel inversion strategy (to avoid large matrix inversion), use Bayesian linearized stochastic inversion to obtain the posterior probability distribution of P-wave and S-wave velocity and density model parameters under seismic data constraints.
[0078] Assume the model parameters m (P-wave velocity, S-wave velocity, density) follow a multivariate Gaussian distribution:
[0079] m~N(μ m , ∑ m (1)
[0080] In the inversion v p v s ρ represents the P-wave velocity, S-wave velocity, and density of channel 1, respectively, and T is the transpose of the matrix; μ m and ∑ m Let m be the mean of the posterior probability distribution and the multivariate covariance matrix.
[0081] In the traditional method, ∑ m It can be obtained from the variation function of the model parameters. This matrix also measures the correlation between different sampling points based on the distance between them. However, in practice, due to the complexity of geological structures, the correlation between two closely spaced sampling points may not be high. It is necessary to introduce information such as the stratigraphic position of seismic interpretation to constrain the correlation.
[0082] By interpreting horizon information through seismic data, the covariance matrix is improved, and a covariance matrix with complex horizon constraints is proposed. When calculating the improved matrix, the horizon information is used as a reference surface to calculate the "true" distance between different sampling points under horizon constraints, and then the covariance function between different sampling points is calculated. Let the improved matrix be ∑ hor Then, the posterior probability distribution of the model parameters under the hierarchical constraints is:
[0083] m~N(μ m , ∑ hor (2)
[0084] The seismic forward modeling expression based on the convolution model is:
[0085] d = Gm + e, (3)
[0086] Where G is the pre-stack forward modeling operator, which can be obtained through the Zoeppritz equation or its approximation, and e is the error term, which is assumed to follow a Gaussian distribution with zero mean:
[0087] e~N(0,∑ e (4)
[0088] In the formula ∑ e Let e be the covariance matrix. According to formulas (2) and (3), the probability distribution of the earthquake record can be obtained:
[0089] d~N(Gμ m , ∑ d (5)
[0090] In the formula, Gμ m Let ∑ be the mean of the earthquake records. d =G∑ hor G T +∑ e Let be the covariance matrix of the seismic signal. Based on multivariate linear Gaussian theory, given the seismic record d, the conditional probability distribution of the model parameter m is still a Gaussian distribution, which can be expressed as:
[0091] m|d~N(μ m|d , ∑ m|d (6)
[0092] In the formula, the posterior mean μ m|d and the posterior multivariate covariance matrix ∑ m|d It can be represented as:
[0093]
[0094] Where, ∑ md =∑ hor G TThe covariance matrix between seismic data and model parameters under the constraint of horizon information. This refers to information related to the model parameter m derived from seismic data. I is defined as the degree to which the uncertainty of seismic data with respect to m is reduced. d and V d For earthquake information; in formula (7), μ m|d Considered as the result of deterministic inversion, ∑ m|d These two parameters are used to characterize the strength of the uncertainty in the inversion result. Based on these two parameters, a stochastic inversion solution is obtained by using LU decomposition simulation or Gibbs simulation.
[0095] Step 3: Based on the stratigraphic information interpreted from the seismic data, construct a complex stratigraphic control operator. Based on geostatistics, under the constraints of the stratigraphic control operator, establish a forward modeling expression between the model parameters and well logging data, and finally obtain the posterior probability distribution of the model parameters under the constraints of well logging data.
[0096] In the inversion, m contains the values of multiple model parameters, and this multiple data includes well-side channels; then the model parameters of the well-side channels can be represented as Hm, where the operator H is a sparse matrix, which is used to select the well-side channels from m. Therefore, the well-seismic relationship in traditional stochastic inversion can be expressed as:
[0097] h = Hm, (8)
[0098] Where h represents well logging data, which is sparsely distributed; however, traditional Bayesian linearization inversion does not consider geological structure. When the geological structure is complex, it is unreasonable to represent the relationship between well logging data and model parameters solely through the linear operator H. Here, seismic interpretation stratigraphic information is introduced to improve the traditional operator, making the relationship between the improved model parameters and well logging data as follows:
[0099] h = H hor m, (9)
[0100] Wherein, matrix H hor Constrained by seismic interpretation results, it is defined as a complex horizon control operator, whose number of rows and columns is controlled by the horizon. When the horizon information is complex, H is no longer a square matrix.
[0101] Based on the above formula, the probability distribution of well logging data is as follows:
[0102]
[0103] Step 4: Based on the assumption that well logging data and seismic data are statistically independent, the posterior probability distributions of model parameters constrained by seismic and well logging data are fused to obtain the probability distribution of model parameters constrained by both well logging and seismic data. Different inversion solutions are obtained through stochastic simulation.
[0104] Based on linear Gaussian theory, under the joint constraints of seismic data d and well logging data h, the conditional probability distribution of the model parameter m is still a Gaussian distribution:
[0105] m|(h,d)~N(μ) m|(h，d) , ∑ m|(h，d) (11)
[0106] Among them, the posterior mean μ of the model parameters under simultaneous well-seismic constraints m|(h，d) and the posterior multivariate covariance matrix ∑ m|(h，d) Represented as:
[0107]
[0108] Where, ∑ hh Let ∑ be the covariance matrix of the well data. hd It is a covariance matrix representing the correlation between well logging data and seismic data, ∑ hh and ∑ hd The specific expression is:
[0109]
[0110] Where, μ m|(h，d) Considered as the result of deterministic inversion under well seismic constraints, ∑ m|(h，d) This method is used to characterize the strength of uncertainty in the inversion result; through stochastic simulation, stochastic inversion solutions with different and equally probable well-seismic simultaneous constraints can be obtained; however, this method has a serious drawback, namely, when solving for the posterior probability distribution, the kernel matrix C must be considered. hd The inversion process, which cannot be avoided through approximation in traditional linear Bayesian stochastic inversion methods, involves C. hd The seismic forward modeling operator G is complex and difficult to invert, and its inversion process is both difficult and time-consuming, which greatly limits the practical application of the traditional BLI method.
[0111] Fast stochastic inversion is based on the assumption that well logging data and seismic data are independent, ignoring the statistical correlation between them. This assumption is based on several reasons: First, seismic and well logging data are two different types of data, acquired using different methods and operating in different frequency bands, resulting in low statistical correlation. Second, seismic data reflects changes in subsurface reflection coefficients, i.e., relative changes in elastic parameters such as velocity; while well logging data represents absolute values of elastic parameters, leading to significant physical differences between the two. Based on this assumption, it can be concluded that:
[0112] ∑ hd =0 (14) Therefore, the kernel matrix in formula (12) can be transformed into:
[0113]
[0114] The above formula can be further transformed into:
[0115]
[0116] Further derivation shows that the above formula can be transformed into the following form:
[0117]
[0118] The above expression can be represented in a concise form as follows:
[0119]
[0120] Among them I b It is information related to the absolute value of the model parameter m from well logging data, V b I represents the degree to which the uncertainty of well logging data is reduced for m; b and V b For well logging information, its expression is:
[0121]
[0122] According to formula (19), by inputting seismic data, well logging data, and stratigraphic data, high-resolution well-seismic joint inversion results under stratigraphic information constraints can be obtained.
[0123] Figure 1 , Figure 2 and Figure 3 The two-dimensional angle gathers used in this embodiment with incident angles of 5, 15, and 25 degrees are shown respectively. This two-dimensional profile contains 600 transverse data points at a distance of 6 kilometers and 400 longitudinal time sampling points with a time range of 600-1000 milliseconds. Figure 1 The angle channel set shown is marked with dashed lines, indicating the three channels W1, W2, and W3 through which the profile passes. W1 and W3 are used for inversion, while W2 is used to verify the inversion effect.
[0124] Figure 4 , Figure 5 and Figure 6 The posterior probability mean values corresponding to the three parameters of P-wave and S-wave velocity and density are displayed. Because only seismic data is used here, the inversion results are relatively smooth and have low resolution, which is completely unacceptable for most oilfield development needs at this stage. Figure 7 Showing Figures 1-3 The target horizon shown in the seismic record, and the data from two wells that participated in the inversion, passing through the horizon.
[0125] Figure 8 , Figure 9 ,Figure 10 The single-shot stochastic inversion results for P-wave and S-wave velocities and densities are presented respectively. It can be seen that the inversion results have high resolution, sufficient to meet the needs of the oilfield during the development phase. At the inversion quality control wells, a one-dimensional comparison of the inversion results is first performed. Figure 11 The comparison results are shown. The solid line represents the well curves of the three pre-stack parameters, while the dashed line represents the random inversion results. It can be seen that the inversion results not only have high resolution, but also have a high degree of agreement with the well data, which is something that conventional deterministic inversion cannot achieve.
[0126] In addition, the distribution histograms of the inversion results well data are compared to determine whether the stochastic inversion method can reasonably characterize the distribution of model parameters. Figure 12 , Figure 13 and Figure 14 The comparison results of P-wave and S-wave velocities and densities are shown separately. The light-colored histogram on the left is the distribution histogram of the inversion results at the verification well location, while the dark-colored histogram on the right represents the distribution histogram of the actual well logging curves. The comparison reveals that the inversion results and the well logging curve histograms are almost identical, which is sufficient to demonstrate that the inversion strategy proposed in this patent can accurately depict the true distribution of subsurface parameters.
[0127] Finally, the method of this invention is compared with the traditional Bayesian stochastic inversion method, and the comparison results are shown in Table 1. The first row of Table 1 represents the number of channels simultaneously inverted during the above two-dimensional inversion test, which are 1, 5, 10, and 15 respectively. The second row shows the time consumed by the method of this invention for two-dimensional inversion. The third row shows the time consumed by the traditional Bayesian linearized stochastic inversion. The fourth row shows the percentage improvement in computational efficiency. The comparison shows that although the method of this invention is based on the traditional Bayesian stochastic inversion, its computational efficiency is improved by about 30%. The essential reason for this improvement is that this method is based on the assumption of statistical independence of well and seismic data, thus eliminating the need for large-scale matrix inverse operations.
[0128] Table 1
[0129]
[0130]
[0131] To address the aforementioned issues in stochastic inversion, this method, based on the assumption of independence between seismic signals and well logging records, improves Bayesian linearized stochastic inversion and proposes an accelerated method. This method separates seismic and well logging information in traditional Bayesian linearized stochastic inversion, retrieving them separately. In the calculation of seismic information, a single-channel stochastic inversion method or conventional deterministic inversion is employed, significantly reducing the dimensionality of the required inverse matrix or directly avoiding matrix inversion using optimization algorithms. This substantially improves the computational efficiency of stochastic inversion. Furthermore, to enhance the ability of stochastic inversion to characterize complex geological structures, complex horizon control operators are constructed using seismically interpreted stratigraphic information and introduced into the fast Bayesian linearized stochastic inversion, ultimately forming a horizon-constrained multivariate Gaussian well-seismic fast stochastic inversion method. Actual data testing shows that the improved fast method is several times more computationally efficient than traditional Bayesian linearized stochastic inversion, while achieving the same inversion results.
Claims
1. A hierarchical constrained multivariate Gaussian fast pre-stack stochastic inversion method, characterized in that: The hierarchical constrained multivariate Gaussian fast pre-stack stochastic inversion method includes the following steps: Step 1: Determine the horizontal and vertical lateral line ranges and the vertical time range of the inversion data; select the experimental wells to participate in the inversion, select the control wells for inversion, and perform well-seismic calibration and wavelet extraction; Step 2: Based on the single-channel inversion strategy, use Bayesian linearized stochastic inversion to obtain the posterior probability distribution of P-wave and S-wave velocity and density model parameters under seismic data constraints. By interpreting horizon information through seismic data, the covariance matrix is improved, and a covariance matrix with complex horizon constraints is proposed. Let the improved matrix be ∑ hor Then, the posterior probability distribution of the model parameters under the hierarchical constraints is: m~N(μ m ,∑ hor ), (2) The seismic forward modeling expression based on the convolution model is: d = Gm + e, (3) Where G is the pre-stack forward operator, and e is the error term, assumed to follow a Gaussian distribution with zero mean: e~N(0,∑ e ), (4) In the formula ∑ e Let e be the covariance matrix. According to formulas (2) and (3), the probability distribution of the earthquake record can be obtained: d~N(Gμ m ,∑ d ), (5) In the formula, Gμ m Let ∑ be the mean of the earthquake records. d =G∑ hor G T +∑ e Let be the covariance matrix of the seismic signal. Based on multivariate linear Gaussian theory, given the seismic record d, the conditional probability distribution of the model parameter m is still a Gaussian distribution, which can be expressed as: m|d~N(μ m|d ,∑ m|d ), (6) In the formula, the posterior mean μ m|d and the posterior multivariate covariance matrix ∑ m|d It can be represented as: Where, ∑ md =∑ hor G T The covariance matrix between seismic data and model parameters under the constraint of horizon information. This refers to information related to the model parameter m derived from seismic data. I is defined as the degree to which the uncertainty of seismic data with respect to m is reduced. d and V d For earthquake information; in formula (7), μ m|d Considered as the result of deterministic inversion, ∑ m|d These two parameters are used to characterize the strength of the uncertainty in the inversion result. Based on these two parameters, a stochastic inversion solution is obtained by using LU decomposition simulation or Gibbs simulation. Step 3: Based on the stratigraphic information interpreted from the seismic data, construct a complex stratigraphic control operator. Under the constraints of the stratigraphic control operator, establish a forward modeling expression between the model parameters and the well logging data, and finally obtain the posterior probability distribution of the model parameters under the constraints of the well logging data. By incorporating seismic interpretation horizon information, the traditional operator is improved, and the relationship between the improved model parameters and well logging data is defined as follows: h=H hor m, (9) Wherein, matrix H hor Constrained by seismic interpretation results, it is defined as a complex horizon control operator, whose number of rows and columns is controlled by the horizon. Based on the above formula, the probability distribution of well logging data is as follows: Step 4: Based on the assumption that well logging data and seismic data are statistically independent, the posterior probability distributions of model parameters constrained by seismic and well logging data are fused to obtain the probability distribution of model parameters constrained by both well logging and seismic data. Different inversion solutions are obtained through stochastic simulation. Based on linear Gaussian theory, under the joint constraints of seismic data d and well logging data h, the conditional probability distribution of the model parameter m is still a Gaussian distribution: m|(h,d)~N(μ m|(h，d) ,∑ m|(h，d) ) (11) Among them, the posterior mean μ of the model parameters under simultaneous well-seismic constraints m|(h，d) and the posterior multivariate covariance matrix ∑ m|(h，d) Represented as: Where, ∑ hh Let ∑ be the covariance matrix of the well data. hd It is a covariance matrix representing the correlation between well logging data and seismic data, ∑ hh and ∑ hd The specific expression is: Where, μ m|(h，d) Considered as the result of deterministic inversion under well seismic constraints, ∑ m|(h，d) This is used to characterize the strength of the uncertainty in the inversion result; through stochastic simulation, stochastic inversion solutions with different and equally probable well-seismic simultaneous constraints can be obtained. Fast stochastic inversion is based on the assumption that well logging data and seismic data are independent, ignoring the statistical correlation between them; based on this assumption, we can conclude that: ∑ hd ＝0 (14) Therefore, the kernel matrix in formula (12) can be transformed into: The above formula can be further transformed into: Further derivation shows that the above formula can be transformed into the following form: The above expression can be represented in a concise form as follows: Among them I b It is information related to the absolute value of the model parameter m from well logging data, V b I represents the degree to which the uncertainty of well logging data is reduced for m; b and V b For well logging information, its expression is: According to formula (19), by inputting seismic data, well logging data, and stratigraphic data, high-resolution well-seismic joint inversion results under stratigraphic information constraints can be obtained.