A seismic structural feature oriented geostatistical modeling method

Abstract

The application discloses a kind of seismic structural feature oriented geostatistical modeling methods, belong to exploration geophysics field;The method comprises: based on the prior information of geostatistics, the multivariate Gaussian prior distribution of underground parameter model is established;From seismic data, structural guidance information is extracted, including calculating gradient, constructing structural tensor and carrying out eigenvalue decomposition to obtain local dominant direction;Accordingly, a structure-guided constraint operator is constructed;Well logging data constraints are combined with structural constraints to form an observation system;The prior distribution is updated using Bayes' theorem to obtain a posterior distribution containing posterior mean and covariance.The application realizes reasonable interpolation along the trend of strata, improves the accuracy and structural consistency of geological modeling, and can quantitatively predict uncertainty.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of exploration geophysics, and particularly relates to a seismic structural feature oriented geostatistical modeling method. BACKGROUND

[0002] The prediction of subsurface parameter distribution is of great significance in geological modeling and oil and gas reservoir characterization. The core is to infer the probability distribution of subsurface parameters using seismic, geological, logging and other multi-source data, and finally to establish a model set of subsurface elastic and physical parameters. In this technical field, the two-point geostatistical method has been widely used for a long time. This method mainly uses the variogram to describe the spatial correlation of model parameters, and can reasonably integrate logging data into the geological model. However, the two-point geostatistical method highly depends on the stationarity assumption of geological bodies, that is, it is assumed that the variogram and covariance of parameters do not change with the absolute position of the sampling point. This assumption makes it difficult to truly reflect the subsurface structural features when facing nonlinear geological structures such as channels, faults and irregular interlayer surfaces, resulting in serious deviations in the established model in complex geological structure areas.

[0003] To overcome the limitations of the two-point statistical method, the multi-point geostatistical method is proposed. This method breaks through the limitations of the two-point method in describing nonlinear geological bodies by introducing training images to capture complex spatial structures. However, the application of multi-point geostatistical method in complex reservoirs still largely depends on the quality and scope of the training images. When the training images cannot accurately represent the actual geological structure, systematic errors will occur in the established model, which constitutes a significant limitation in actual work area applications.

[0004] In recent years, introducing seismic data as local directional constraints has become an important way to enhance the structural consistency of geostatistical modeling. This type of method extracts the local principal direction, dip angle and structural plane orientation of seismic bodies through structure tensor analysis, and uses them to guide the interpolation process. However, most of these methods belong to deterministic interpolation methods, which can improve the geological rationality of the interpolation results, but lack quantitative characterization of the uncertainty of the prediction results. In oil and gas exploration and development decision-making, understanding the uncertainty range of the prediction results is directly related to risk assessment. This defect seriously limits the application value of this type of method in actual production.

[0005] Therefore, the exploration geophysical field urgently needs to develop a modeling method that can simultaneously consider geological structure rationality and uncertainty quantification, which not only overcomes the dependence of two-point statistics on stationarity, but also avoids the serious dependence of multi-point statistics on training images, while also achieving random modeling based on maintaining structural guidance advantages. This technical problem has become the direction of long-term efforts of technical personnel in this field. SUMMARY

[0006] To solve the above technical problems, the application provides a seismic structural feature oriented geostatistical modeling method, which realizes reasonable interpolation along the stratum trend, improves the accuracy and structural consistency of geological modeling, and can quantitatively predict uncertainty.

[0007] To achieve the above object, the application provides a seismic structural feature oriented geostatistical modeling method, which comprises the following steps:

[0008] Obtaining seismic data and logging data of an underground area;

[0009] Based on the geostatistical prior information, a prior probability distribution of the underground parameter model is established, and the prior probability distribution is a multivariate Gaussian distribution;

[0010] Extracting structural oriented information from the seismic data, wherein the extraction of the structural oriented information comprises calculating the gradient of the seismic data, constructing a structural tensor, and performing eigenvalue decomposition on the structural tensor to obtain a local dominant direction;

[0011] Based on the local dominant direction, a structural oriented constraint operator is constructed;

[0012] The logging data and the structural oriented constraint operator are combined to form a joint observation system;

[0013] Based on the Bayes theorem, the prior probability distribution is updated by using the joint observation system to obtain a posterior probability distribution, and the posterior probability distribution comprises a posterior mean and a posterior covariance.

[0014] Optionally, obtaining seismic data and logging data of an underground area comprises:

[0015] The seismic data is a post-stack seismic profile, which is generated by convolution of the reflection coefficient of the underground parameter model and a Ricker wavelet, and random noise defined by signal-to-noise ratio can be added to simulate actual observation conditions;

[0016] The logging data is P-wave velocity data measured at multiple well sites in the underground area, and the well sites are distributed at different trace positions of the seismic profile.

[0017] Optionally, establishing a prior probability distribution of an underground parameter model comprises:

[0018] The underground parameter model is a P-wave velocity model or a wave impedance model, and the underground parameter model is flattened into a one-dimensional vector;

[0019] The mean of the prior probability distribution reflects the prior cognition of the overall level of the underground parameter;

[0020] The covariance matrix of the prior probability distribution is calculated by a variogram, and is used to describe the spatial correlation of the model parameters.

[0021] Optionally, the extracting the structural orientation information comprises:

[0022] calculating a horizontal gradient and a vertical gradient of the seismic data;

[0023] constructing a structure tensor based on the horizontal gradient and the vertical gradient;

[0024] performing eigenvalue decomposition on the structure tensor to obtain a large eigenvalue and a small eigenvalue, and a corresponding eigenvector;

[0025] calculating an angle of a local dominant direction based on the eigenvector, the angle reflecting an orientation of a formation structure in the seismic profile.

[0026] Optionally, the constructing the structural orientation constrained operator comprises:

[0027] rotating a Cartesian coordinate system to a direction aligned with the eigenvector based on the angle of the local dominant direction;

[0028] constructing a first-order difference operator along a formation direction and a first-order difference operator perpendicular to the formation direction;

[0029] combining the first-order difference operators to form the structural orientation constrained operator for smoothing the model along the formation direction.

[0030] Optionally, the forming the joint observation system comprises:

[0031] associating the well logging data with the subsurface parameter model through a sparse sampling matrix to establish a linear observation relationship;

[0032] applying the structural orientation constrained operator to the subsurface parameter model to obtain a structural constraint term;

[0033] combining the well logging data and the structural constraint term into a joint observation vector and a joint observation matrix.

[0034] Optionally, the updating the prior probability distribution comprises:

[0035] calculating a posterior mean and a posterior covariance through a Bayesian linear Gaussian update formula based on the joint observation system;

[0036] the posterior mean representing an optimal model estimation fused with the well logging data and the structural constraint;

[0037] the posterior covariance representing a measure of uncertainty of the model estimation.

[0038] Optionally, the posterior probability distribution is used to generate a predictive model of the subsurface parameter distribution and an uncertainty quantification, the uncertainty quantification being represented by a confidence interval.

[0039] The present application discloses a kind of seismic structural feature oriented geostatistical modeling methods, by fusing seismic structural feature and geostatistical prior, effectively improve the accuracy and practicality of underground parameter modeling.Structure tensor constraint is introduced to accurately depict the local structural direction in seismic data, so that the interpolation process can be extended along the true stratigraphic strike, significantly improve the geological rationality of interpolation result, overcome the distortion problem of traditional method in complex structure area.Geostatistical framework based on Bayesian theory realizes the organic fusion of well logging data and seismic structure information, not only improves the accuracy and stability of interwell attribute prediction, but also can quantify the uncertainty of prediction result, provides reliable basis for reservoir risk assessment.The method reduces the dependence on the quality of training image, enhances its applicability under real geological conditions.At the same time, under the condition of actual data containing noise, good interpolation stability can be maintained by parameter adjustment, and strong structural recognition ability is also shown under the condition of less well logging number, with strong engineering practical value. BRIEF DESCRIPTION OF DRAWINGS

[0040] The accompanying drawings, which form a part of this application, are included to provide a further understanding of the application and are incorporated in and constitute a part of this application. The embodiments of this application, and of the accompanying drawings, are intended to explain the present application, and are not limiting of the present application. In the drawings:

[0041] Figure 1 It is a schematic diagram of the relationship between variogram and covariance of the embodiment of the present application;

[0042] Figure 2 It is a schematic diagram of the relationship between Cartesian coordinate system and stratigraphic coordinate of the embodiment of the present application;

[0043] Figure 3 It is a schematic diagram of the flow of a kind of seismic structural feature oriented geostatistical modeling method of the embodiment of the present application;

[0044] Figure 4 It is a schematic diagram of synthetic model and synthetic seismic section of the embodiment of the present application, wherein (a) is P-wave velocity;(b) is post-stack seismic section;(c) is post-stack seismic section containing noise;

[0045] Figure 5 It is a schematic diagram of interpolation result corresponding to different well number under noiseless condition of the embodiment of the present application, wherein (a) is seven well positions;(b) is the interpolation result corresponding to seven wells;(c) is five well positions;(d) is the interpolation result corresponding to five wells;(e) is three well positions;(f) is the interpolation result corresponding to three wells;

[0046] Figure 6Fig. 1 is a schematic diagram of interpolation results corresponding to different well numbers under a noise condition according to an embodiment of the present application, wherein (a) is the interpolation result corresponding to seven well positions; (b) is the interpolation result corresponding to five well positions; and (c) is the interpolation result corresponding to three well positions;

[0047] Figure 7 Fig. 2 is a one-dimensional schematic diagram of interpolation results and uncertainty estimation results corresponding to different well numbers under a noiseless condition according to an embodiment of the present application, wherein (a) is the interpolation result corresponding to the post-stack seismic profile; and (b) is the 95% confidence interval;

[0048] Figure 8 Fig. 3 is a comparison diagram of interpolation results when the parameter Var takes different values under a noiseless condition according to an embodiment of the present application, wherein (a) is the interpolation result corresponding to the noiseless condition when Var = 10 -3 ; (b) is the interpolation result corresponding to the noiseless condition when Var = 10 -4 ; (c) is the interpolation result corresponding to the noiseless condition when Var = 10 -5 ; (d) is the interpolation result corresponding to the noiseless condition when Var = 10 -6 ; (e) is the interpolation result corresponding to the noiseless condition when Var = 10 -7 ; (f) is the interpolation result corresponding to the noiseless condition when Var = 10 -8 ;

[0049] Figure 9 Fig. 4 is a comparison diagram of interpolation results when the parameter Var takes different values under a noise condition according to an embodiment of the present application, wherein (a) is the interpolation result corresponding to the noise condition with 8 dB signal-to-noise ratio when Var = 10 -3 ; (b) is the interpolation result corresponding to the noise condition with 8 dB signal-to-noise ratio when Var = 10 -4 ; (c) is the interpolation result corresponding to the noise condition with 8 dB signal-to-noise ratio when Var = 10 -5 ; (d) is the interpolation result corresponding to the noise condition with 8 dB signal-to-noise ratio when Var = 10 -6 ; (e) is the interpolation result corresponding to the noise condition with 8 dB signal-to-noise ratio when Var = 10 -7 ; (f) is the interpolation result corresponding to the noise condition with 8 dB signal-to-noise ratio when Var = 10 -8 ;

[0050] Figure 10 Fig. 5 is a schematic diagram of four well positions and a real seismic profile according to an embodiment of the present application, wherein (a) is the four well positions; and (b) is the post-stack seismic profile;

[0051] Figure 11Figures showing interpolation results of real data corresponding to different well numbers in embodiments of the present application, wherein (a) is the location of three wells; (b) is the interpolation result corresponding to the three wells; (c) is the location of two wells; (d) is the interpolation result corresponding to the two wells; (e) is the location of another group of two wells; and (f) is the corresponding interpolation result;

[0052] Figure 12 Figures showing the interpolation results of Var = 10 -3 Figures showing the 95% confidence interval of the 350th trace corresponding to the interpolation results when Var = 10. DETAILED DESCRIPTION

[0053] It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict. The present application will be described in detail below with reference to the accompanying drawings and in combination with embodiments.

[0054] It should be noted that the steps shown in the flowchart of the accompanying drawings can be executed in a computer system such as a group of computer executable instructions, and although the logical order is shown in the flowchart, in some cases, the steps shown or described herein can be executed in a different order.

[0055] As Figure 3 shown, the present embodiment provides a seismic structure feature oriented geostatistical modeling method, comprising:

[0056] obtaining seismic data and logging data of an underground region;

[0057] establishing a prior probability distribution of an underground parameter model based on geostatistical prior information, the prior probability distribution being a multivariate Gaussian distribution;

[0058] extracting structure oriented information from the seismic data, the extraction of structure oriented information comprising calculating the gradient of the seismic data, constructing a structure tensor, and performing eigenvalue decomposition on the structure tensor to obtain a local dominant direction;

[0059] constructing a structure oriented constraint operator based on the local dominant direction;

[0060] combining the logging data and the structure oriented constraint operator to form a joint observation system;

[0061] updating the prior probability distribution based on Bayes' theorem using the joint observation system to obtain a posterior probability distribution, the posterior probability distribution comprising a posterior mean and a posterior covariance.

[0062] Further, obtaining seismic data and logging data of an underground region comprises:

[0063] The seismic data is a post-stack seismic profile, which is generated by convolving the reflection coefficient of the subsurface parameter model with a Ricker wavelet, and adding random noise defined by the signal-to-noise ratio to simulate the actual observation conditions;

[0064] The logging data is the P-wave velocity data measured at multiple well positions in the subsurface area, and the well positions are distributed at different trace positions of the seismic profile.

[0065] Specifically, the implementation process of this embodiment includes:

[0066] Interpolation under the constraint of the structure tensor:

[0067] By introducing the structure tensor to characterize the directional features in the seismic data, extracting the local gradient information, and establishing a constraint operator that can describe the main tectonic direction and anisotropic features. This constraint enables the interpolation process to extend along the actual formation strike, thereby improving the consistency between the interpolation result and the true geological structure while maintaining the tectonic continuity.

[0068] Well interpolation based on geostatistical prior information:

[0069] Using Bayesian theory to introduce geostatistical prior information, generating the posterior mean and covariance, and realizing the organic combination of prediction results and uncertainty analysis. This method not only improves the accuracy and stability of well-to-well property prediction, but also can quantify the uncertainty of the prediction results, providing more reliable technical support for reservoir modeling and risk assessment.

[0070] Synthetic data testing:

[0071] Use a two-dimensional P-wave velocity model as the benchmark model, and its model is as shown in Figure 4 (a). The scale of this model is 187 nodes in the vertical (time) direction and 801 nodes in the horizontal (CDP) direction. The distance between every two CDPs is 10 meters. Its overall structure is complex, with obvious horizontal interlayer velocity variations and local fracture features, and can better simulate common geological structure situations. Figure 4 (b) of this model is the noiseless post-stack seismic profile corresponding to this model, and this post-stack seismic profile is obtained by convolving the model reflection coefficient with a Ricker wavelet. Adding random noise with a signal-to-noise ratio of 8 dB to the trace gather of this post-stack seismic profile gives Figure 4 (c), so as to simulate more realistic observation data for comparative testing.

[0072] To test the influence of the number of wells on the accuracy of the interpolation model, this embodiment sets three cases of the number of well constraints and compares their interpolation results. As shown in Figure 5(a) shows the position of well data, seven pseudo wells are used, corresponding to the real P-wave velocity model of 100th, 200th, 300th, 400th, 500th, 600th and 700th traces respectively, Figure 5 (b) is the interpolation result of the seven pseudo wells, and it is found that the result can more completely restore the main structural features in the velocity model and is close to the real model. Figure 5 (c) shows that the well positions are the 100th, 200th, 400th, 600th and 700th traces, and the corresponding interpolation results are as shown in (d), Figure 5 (d) shows the result, and the details of the result are weakened, especially in the complex structure area. Figure 5 (e) shows that the well positions are the 150th, 350th and 600th traces, and the corresponding interpolation results are as shown in (f), Figure 5 (f) shows that the resolution of the complex structure area is further weakened. The interpolation method has a certain degree of decline in the accuracy of the interpolation result with the decrease of the number of wells, but it can still be seen that even in the case of a small number of wells, the method still has strong structural identification ability and can accurately restore the underground geological structure features. The embodiment can restore the geological structure under the condition of limited well data, and has strong structural identification ability.

[0073] In the above three well constraint cases, random noise with a signal-to-noise ratio of 8 dB is added to test the performance. Figure 6 Fig. 8 is a schematic diagram of the interpolation results corresponding to different numbers of wells under the condition of noise according to the embodiment of the present application, wherein (a) is the interpolation result corresponding to the position of seven wells; (b) is the interpolation result corresponding to the position of five wells; (c) is the interpolation result corresponding to the position of three wells. The results show that the overall interpolation result has a certain degree of decline in accuracy under the condition of noise and the decrease of the number of wells, and the amplitude is related to the number of wells. Under the condition of seven wells, the influence of noise is small, and the interpolation result is close to the condition without noise, which indicates that the more the number of wells is, the stronger the noise resistance performance of the interpolation is.

[0074] Further, the prior probability distribution of the underground parameter model comprises:

[0075] The underground parameter model is a P-wave velocity model or a wave impedance model, and the underground parameter model is flattened into a one-dimensional vector;

[0076] The mean of the prior probability distribution reflects the prior cognition of the overall level of the underground parameter;

[0077] The covariance matrix of the prior probability distribution is calculated by a variation function and is used to describe the spatial correlation of the model parameters.

[0078] Specifically, the implementation process of the embodiment comprises:

[0079] Well interpolation method based on prior information from geostatistics: Let... This is a parametric model, obtained by flattening a two-dimensional P-wave velocity model and arranging it into a one-dimensional vector in column-major order. This arrangement is universal and also applicable to subsurface parameters such as wave impedance and S-wave velocity. (The parametric model is assumed to be a one-dimensional vector obtained by flattening a two-dimensional P-wave velocity model in column-major order.) Follows a multivariate Gaussian distribution:

[0080] (1);

[0081] in, Parametric model (a one-dimensional vector obtained by flattening the two-dimensional P-wave velocity model column-first). This is represented as the prior mean of the model, reflecting the prior knowledge of the overall level of this parameter. Represented as the prior covariance matrix, A multivariate Gaussian distribution can describe the spatial correlation of a model at different locations.

[0082] In practical applications, well data provides some known information, and geostatistical methods are used to establish a linear relationship to estimate subsurface parameters. Let... For locations with partial well logging data vectors, the locations without logging data are: . For the first The well data is then compared with the parametric model (a one-dimensional vector obtained by flattening the two-dimensional P-wave velocity model column-first). The following linear relationship exists between them:

[0083] (2);

[0084] in, For locations with partial logging data vectors, the locations without logging data are... ; This represents the logging error term, primarily stemming from logging instrument errors and the discrepancy between well data and actual data. It includes logging data vectors for some well locations. for:

[0085] (3);

[0086] matrix Defined as:

[0087] (4);

[0088] in, For the first Well data for the orifice; T is the matrix transpose; It is the identity matrix. The purpose of sparse sampling matrix is to pick up the model data at well locations for fitting well data. This way of construction makes a clear relationship between well data and model.

[0089] In geostatistics, variogram and covariance function are the core tools to characterize spatial correlation. Spatial characteristics of subsurface media are measured by variogram, while covariance function has the following relationship with variogram as shown in Figure 1

[0090] (5);

[0091] where, is the maximum value of covariance at zero distance, i.e. variance, is the covariance at any distance h, is the corresponding variogram, is the sill value.

[0092] Therefore, based on the linear Gaussian principle, its distribution function is obtained as:

[0093] (6);

[0094] where, and represent the prior prediction value and joint covariance matrix of model at well locations, respectively. The joint covariance is composed of the prediction variance matrix and the well error covariance matrix , which comprehensively considers the errors of both and determines the degree of response of model adjustment to data.

[0095] Since the parameter model (a one-dimensional vector obtained by column-first flattening of two-dimensional P-wave velocity model) and the well data vector with partial well locations are both linear Gaussian systems, their joint distribution is also Gaussian, so the following joint Gaussian form is obtained:

[0096] (7);

[0097] According to the conditional distribution property of Gaussian joint distribution, the posterior distribution of parameter model (a one-dimensional vector obtained by column-first flattening of two-dimensional P-wave velocity model) under the condition of well data vector with partial well locations is obtained as:

[0098] (8);

[0099] where,​​ is a parameter model (a one-dimensional vector obtained by column-wise flattening a two-dimensional P-wave velocity model); is a vector of logging data at the owning wellsite; represents a multivariate Gaussian distribution; is a vector of logging data at the owning wellsite is the model posterior mean under the condition is a vector of logging data at the owning wellsite is the model posterior covariance matrix under the condition

[0100] is a vector of logging data at the owning wellsite is the model posterior mean under the condition and the posterior covariance The expression of the posterior mean is:

[0101] (9);

[0102] (10);

[0103] wherein, is a vector of logging data at the owning wellsite is the model posterior mean under the condition is a vector of logging data at the owning wellsite is the model posterior covariance matrix under the condition represents a prior mean of the model; represents a prior covariance matrix; is a sparse sampling matrix; T is a matrix transpose; is a logging error covariance matrix; is a vector of logging data at the owning wellsite.

[0104] Further, the structural guiding information includes:

[0105] calculating a horizontal gradient and a vertical gradient of the seismic data;

[0106] constructing a structure tensor based on the horizontal gradient and the vertical gradient;

[0107] performing eigenvalue decomposition on the structure tensor to obtain a large eigenvalue and a small eigenvalue, and a corresponding eigenvector;

[0108] calculating an angle of a local dominant direction based on the eigenvector, the angle reflecting an orientation of a formation structure in the seismic profile.

[0109] Specifically, the implementation process of the embodiment includes:

[0110] Interpolation under structural tensor constraint: In order to further improve the geological rationality of the interpolation result, the core of the embodiment is to combine the structural information of the seismic data and the quantitative constraint of the logging data as the constraint condition of the interpolation model. Specifically, the gradient of the seismic profile is calculated to obtain the structural tensor matrix in the horizontal and vertical directions, and then the eigenvalue decomposition of the structural tensor is performed to extract the local dominant direction as the guiding information of the geological structure. The structural guiding information is integrated into the interpolation model to construct a joint Gaussian model.

[0111] The gradient components in the horizontal and vertical directions of the known seismic profile are represented as and The structural tensor is constructed according to the gradient, and the structural tensor is represented as:

[0112] (11);

[0113] wherein, is the structural tensor; is the gradient component in the horizontal direction of the seismic profile; is the gradient component in the vertical direction of the seismic profile;

[0114] The eigenvalue decomposition is performed on the structural tensor, and the diagonal line is arranged in descending order to obtain

[0115] (12);

[0116] wherein, is the eigenvalue matrix of the diagonalized structural tensor; is the local structural vertical eigenvalue of the seismic profile; is the eigenvalue in the parallel direction of the local structure of the seismic profile; and the corresponding eigenvector is obtained as

[0117] (13);

[0118] (14);

[0119] wherein, is the eigenvector parallel to the local structural direction of the seismic profile; is the eigenvector perpendicular to the local structural direction of the seismic profile; is the angle between the local structural direction of the seismic profile and the Cartesian coordinate system; and T is the matrix transpose.

[0120] Further, the structural guiding constraint operator includes:

[0121] Rotating the Cartesian coordinate system to the direction aligned with the eigenvector based on the angle of the local dominant direction;

[0122] Constructing the first-order difference operator along the strata direction and the first-order difference operator along the vertical strata direction;

[0123] Combining the first-order difference operators to form a structure-oriented constraint operator for smoothing the model along the strata direction.

[0124] Specifically, the implementation process of the embodiment includes:

[0125] By rotating the Cartesian coordinate system to the direction aligned with the eigenvector, a new difference operator is constructed. The operator enables the model to be smoothed along the local structure direction of the subsurface, and the rotation is as follows Figure 2 As shown in the figure, the black arrow represents and the purple dashed arrow represents and the purple dashed arrow represents and the purple dashed arrow represents the direction of the converted coordinate system, which is consistent with the direction of the eigenvector perpendicular and parallel to the local structure of the seismic profile; the green arrow represents the vector and ; is the angle between the Cartesian coordinate system and the direction of the converted coordinate system; in Figure 2 the black arrow represents and the black arrow represents the direction of the Cartesian coordinate system, and the Cartesian coordinate system is converted into a coordinate system with the same direction as the eigenvector perpendicular and parallel to the local structure of the seismic profile according to the angle information, in Figure 2 the coordinate system represented by the purple dashed arrow, and the first-order derivative operator in the horizontal direction and the first-order derivative operator in the vertical direction are converted into the structure operators and . The relationship between the structure operators and the gradient operator is represented as:

[0126] (15);

[0127] wherein, is the structure operator in the horizontal direction; is the structure operator in the vertical direction; is the eigenvector parallel to the local structure direction of the seismic profile; is the eigenvector perpendicular to the local structure direction of the seismic profile; is the first-order difference operator in the horizontal direction of the Cartesian coordinate system; is the first-order difference operator in the vertical direction of the Cartesian coordinate system; T is the matrix transpose;

[0128] When the model size is M×N, the operator is further written in matrix form:

[0129] (16) ;

[0130] (17) ;

[0131] (18) ;

[0132] (19) ;

[0133] (20) ;

[0134] (21) ;

[0135] wherein, is a diagonal matrix composed of different ; is a diagonal matrix composed of different ; M x N represents model size;

[0136] The isotropic first-order difference is weighted according to the local structure direction, so that a more actual structure conforming to local characteristics is obtained.

[0137] Further, forming the joint observation system comprises:

[0138] associating the well logging data with the subsurface parameter model through a sparse sampling matrix, to establish a linear observation relationship;

[0139] applying the structure-oriented constraint operator to the subsurface parameter model, to obtain a structure constraint term;

[0140] combining the well logging data and the structure constraint term into a joint observation vector and a joint observation matrix.

[0141] Specifically, the implementation process of the embodiment comprises:

[0142] Thus, a constraint value in the target structure direction is obtained:

[0143] (22) ;

[0144] wherein, is a structure constraint term; is a structure error term; is a horizontal direction structure operator; is a structure error term; is a parameter model (a one-dimensional vector obtained by column-priority flattening a two-dimensional P-wave velocity model);

[0145] combining with the structure constraint term to construct a following joint Gaussian observation system:

[0146] (23);

[0147] wherein, is a structural constraint term; is a well data vector with partial well locations; is a horizontal structural operator; is a sparse sampling matrix, aiming to pick up model data at well locations; is a subsurface parameter model vector; is a structural error term; is a well error term.

[0148] Further, updating the prior probability distribution comprises:

[0149] calculating a posterior mean and a posterior covariance based on the joint observation system via a Bayesian linear Gaussian update formula;

[0150] the posterior mean represents an optimal model estimation that fuses well data and structural constraints;

[0151] the posterior covariance represents a measure of uncertainty of the model estimation.

[0152] Specifically, the implementation process of the embodiment comprises:

[0153] Meanwhile, taking into account well data constraints and structural constraints improves the accuracy of interpolation, and the joint posterior mean and the joint posterior covariance are:

[0154] (24);

[0155] (25);

[0156] wherein, is a joint posterior mean; is a joint posterior covariance; represents a prior mean for the model; represents a prior covariance matrix; is a sparse sampling matrix; T is a matrix transpose; is a well error covariance matrix; is a structural error covariance matrix; is a well data vector with partial well locations; is a horizontal structural operator; is a structural constraint term;

[0157] define the structural error covariance matrix Its diagonal element is Var, and the impact of this parameter on the accuracy of the interpolation results will be discussed in the subsequent data testing section. The advantage of this embodiment is that it not only ensures that the interpolation results are consistent with the actual structure at the well point, but also makes the spatial distribution of the results more consistent with the direction of the stratigraphic structure.

[0158] Furthermore, the posterior probability distribution is used to generate a prediction model for the distribution of underground parameters and to quantify the uncertainty, which is represented by confidence intervals.

[0159] Specifically, the implementation process of this embodiment includes:

[0160] In this embodiment, section 450 is selected as a representative for single-section comparative analysis. Figure 7 (a) shows the comparison curves of different interpolation results and the real model on this profile, respectively. Figure 4 (a) The real model and Figure 5 The 450th model data in (b), (d), and (f). Figure 7 (b) shows Figure 4 The figure shows the 95% confidence interval and posterior mean of well number 450 in (b). In the figure, the black, blue, orange, and green lines and the light blue shading represent the true model, the interpolation results corresponding to the seven well locations, the interpolation results corresponding to the five well locations, the interpolation results corresponding to the three well locations, and the 95% confidence interval, respectively. The results show that... Figure 7 In (a), the interpolation curves of the seven wells almost coincide with the real model, indicating that the interpolation accuracy is extremely high; the interpolation curves of the five wells deviate slightly from the real curves at local structures, but the overall trend is consistent; the interpolation curves of the three wells are slightly worse. Figure 7 In (b), the light blue shading of the confidence interval contains the true model, indicating that the interpolation method can reasonably capture the uncertainty of the model and effectively show the range of uncertainty of the interpolation results.

[0161] In the parameter sensitivity analysis, this embodiment studied the impact of parameter Var on the accuracy of interpolation results. In this embodiment, Var was selected as 10. -3 10 -4 10 -5 10 -6 10 -7 10 -8 Numerical comparison tests were conducted, using interpolation with seven dummy wells under both noise-free and noisy conditions. The interpolation results under noise-free conditions are as follows: Figure 8 As shown in (a)-(f), the corresponding Var values ​​are 10 respectively. -3 10 -4 10 -5 10 -6 10 -7 and 10-8 It is found that the change of the interpolation result is not significant, indicating that the size of the parameter Var between the values has little effect on the interpolation result.

[0162] By calculating the absolute error between the interpolation result and the true model, it is found that as the value of Var decreases from 10 -3 to 10 -8 , the error of the interpolation result gradually increases, indicating that too small Var value leads to the reduction of the interpolation accuracy. It can be seen that in the case of no noise, smaller Var value will lead to poor interpolation effect, which may be caused by numerical instability or model overfitting. In the case of no noise, Var = 10 -3 achieves good interpolation performance and has good model resolution under the experimental setting.

[0163] Next, the influence of the value of Var on the accuracy of the interpolation result in the presence of noise is tested. Random noise with a signal-to-noise ratio of 8dB is added to the original data, and the interpolation results are shown in (a)-(f) of Figure 9 , and the corresponding Var values are 10 -3 , 10 -4 , 10 -5 , 10 -6 , 10 -7 and 10 -8 . It is found that when Var = 10 -3 , 10 -4 , the interpolation image depth area structure is discontinuous, and the accuracy is low, especially when Var = 10 -3 , the interpolation quality is the worst. When Var is reduced to 10 -5 , 10 -6 , 10 -7 and 10 -8 , the interpolation result is significantly improved, the interpolated model structure is more continuous and real, and the shape is closer to the true model, indicating that in the presence of noise, appropriately reducing the value of Var helps to suppress the interference of noise on the interpolation accuracy and improves the interpolation stability.

[0164] In the actual situation of adding noise, selecting appropriate Var value is crucial to improve the accuracy of the interpolation result. Overall, in the presence of noise, Var = 10 -5 , 10 -6 and 10 -7 achieve good interpolation performance and have good model resolution under the experimental setting.

[0165] An application example of the present application is:

[0166] This seismic data profile passes through four real wells, and the well positions are as shown in Figure 10(a) shows the locations of channels 227, 309, 350, and 380, indicated by solid red lines. The seismic profile is as follows. Figure 10 (b) shows 600 seismic traces and a time window of 260 milliseconds for each trace.

[0167] To test the effect of different numbers of wells on the interpolation results, the locations of the three wells are as follows: Figure 11 (a) shows lanes 227, 309, and 380, indicated by solid red lines. Figure 11 (b) is its corresponding interpolation result, and the locations of the two wells are as follows: Figure 11 (c) shows lanes 227 and 380. Figure 11 (d) is its corresponding interpolation result. The locations of the two wells are as follows: Figure 11 (e) shows lanes 309 and 380. Figure 11 f is its corresponding interpolation result.

[0168] To further verify the applicability of the proposed method, the observation error variance was set to Var=10. −3 The data from well No. 350 was used as a blind well for testing, such as... Figure 12 As shown, the black line represents the actual well data for well number 350, the green line represents the posterior mean of the two wells located at well numbers 309 and 380, and the light green shading represents their corresponding 95% confidence intervals. From... Figure 12 The results show that the 95% confidence interval almost completely encompasses the true model, and the posterior mean is largely consistent with the true profile in terms of the location and trend of the main layers. This indicates that the method can maintain good stability and prediction accuracy in real seismic data and has strong engineering applicability.

[0169] This invention discloses a seismic structure feature-guided geostatistical modeling method. By integrating seismic structural features with geostatistical priors, it effectively improves the accuracy and practicality of subsurface parameter modeling. Introducing structural tensor constraints accurately characterizes local structural directions in seismic data, allowing the interpolation process to extend along the actual stratigraphic strike, significantly improving the geological rationality of the interpolation results and overcoming the distortion problem of traditional methods in complex structural regions. The geostatistical framework based on Bayesian theory achieves the organic integration of well logging data and seismic structural information, not only improving the accuracy and stability of inter-well attribute prediction but also quantifying the uncertainty of prediction results, providing a reliable basis for reservoir risk assessment. This method reduces its dependence on the quality of training images, enhancing its applicability under real geological conditions. Furthermore, even with noisy real-world data, it maintains good interpolation stability through parameter adjustment and exhibits strong structural identification capabilities even with a limited number of wells, demonstrating significant engineering practical value.

[0170] The above merely provides the preferred embodiments of the present application, and the protection scope of the present application is not limited thereto, and any changes or substitutions within the technical scope disclosed by the present application should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A seismic structural feature directed geostatistical modeling method characterized by, The method comprises: acquiring seismic data and logging data of a subsurface region; establishing a prior probability distribution of a subsurface parameter model based on geostatistical prior information, the prior probability distribution being a multivariate Gaussian distribution; extracting structure-oriented information from the seismic data, the extracting structure-oriented information comprising calculating a gradient of the seismic data, constructing a structure tensor, and performing eigenvalue decomposition on the structure tensor to obtain a local dominant direction; constructing a structure-oriented constraint operator based on the local dominant direction; combining the logging data and the structure-oriented constraint operator to form a joint observation system; updating the prior probability distribution based on Bayes' theorem using the joint observation system to obtain a posterior probability distribution, the posterior probability distribution comprising a posterior mean and a posterior covariance; the extracting structure-oriented information comprises: calculating a horizontal gradient and a vertical gradient of the seismic data; constructing a structure tensor based on the horizontal gradient and the vertical gradient; performing eigenvalue decomposition on the structure tensor to obtain a large eigenvalue and a small eigenvalue, and corresponding eigenvectors; calculating an angle of the local dominant direction based on the eigenvectors, the angle reflecting an orientation of a formation structure in a seismic profile; the constructing a structure-oriented constraint operator comprises: rotating a Cartesian coordinate system to a direction aligned with the eigenvectors based on the angle of the local dominant direction; constructing a first-order difference operator along a formation direction and a first-order difference operator perpendicular to the formation direction; combining the first-order difference operators to form the structure-oriented constraint operator for smoothing the model along the formation direction; the forming a joint observation system comprises: associating the logging data with the subsurface parameter model through a sparse sampling matrix to establish a linear observation relationship; applying the structure-oriented constraint operator to the subsurface parameter model to obtain a structure constraint term; combining the logging data and the structure constraint term into a joint observation vector and a joint observation matrix.

2. The seismic structural feature-directed geostatistical modeling method of claim 1, wherein, The acquiring seismic data and logging data of a subsurface region comprises: the seismic data is a post-stack seismic profile, which is generated by convolution of a subsurface parameter model reflection coefficient and a Ricker wavelet, and random noise defined by a signal-to-noise ratio can be added to simulate actual observation conditions; the logging data is P-wave velocity data measured at multiple well locations in the subsurface region, and the well locations are distributed at different trace positions of the seismic profile.

3. The seismic structural feature-directed geostatistical modeling method of claim 1, wherein, The establishing a prior probability distribution of a subsurface parameter model comprises: the subsurface parameter model is a P-wave velocity model or a wave impedance model, and the subsurface parameter model is flattened into a one-dimensional vector; the mean of the prior probability distribution reflects a prior cognition of a subsurface parameter overall level; the covariance matrix of the prior probability distribution is calculated by a variogram to describe spatial correlation of the model parameters.

4. The seismic structural feature-directed geostatistical modeling method of claim 1, wherein, The updating the prior probability distribution comprises: calculating a posterior mean and a posterior covariance based on the joint observation system through a Bayes linear Gaussian update formula; the posterior mean represents an optimal model estimation fused with the logging data and the structure constraint; the posterior covariance represents an uncertainty measure of the model estimation.

5. The seismic structural feature-directed geostatistical modeling method of claim 1, wherein, The posterior probability distribution is used to generate a predictive model of the subsurface parameter distribution and a quantification of uncertainty, expressed by confidence intervals.

Images