Five-dimensional seismic data prediction method and device

Abstract

The invention relates to the technical field of seismic exploration, and provides a five-dimensional seismic data prediction method and device. The method comprises the following steps: carrying out inclination angle imaging enhancement processing on five-dimensional seismic data; carrying out azimuth time window analysis based on the inclination angle imaging data, and calculating coherent attributes; carrying out gradient mode attributes on the basis of coherent attributes, and carrying out orientation attribute fusion; five-dimensional phase control step-by-step inversion is carried out under fusion attribute constraint; and predicting an oil and gas reservoir of the fractured volcanic rock based on an inversion result. According to the method, the fractured oil and gas reservoir is predicted based on the five-dimensional seismic data, interrelation between adjacent seismic traces is considered, and underground medium parameters are accurately obtained in combination with seismic sensitive attributes and a multi-trace inversion method, so that the problems that single attribute prediction is high in multiplicity of solutions and space description of a fracture system is not comprehensive enough are solved, and the prediction accuracy of the fractured oil and gas reservoir is improved. And meanwhile, the seismic inversion stability and continuity are improved, and high-precision fractured reservoir prediction is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration, and particularly relates to a five-dimensional seismic data prediction method, a five-dimensional seismic data prediction device, and a computer device. Background Art

[0002] Complex fractured hydrocarbon reservoirs such as carbonate rocks, tight sandstones, shales, and volcanic rocks are key areas for oil and gas exploration and development. The characteristics of multi-field, multi-phase tectonic movements, and multi-basin types result in a more complex microstructure of fractured reservoirs and diverse fracture development characteristics, urgently requiring innovation in geophysical prediction technologies for fractured hydrocarbon reservoirs.

[0003] Currently, the prediction technologies for fractured reservoirs in oilfield exploration and development mainly include two categories. One is the attribute analysis technology based on seismic data, which mainly uses attributes such as coherence and curvature to qualitatively predict medium and large-scale fractures, and uses tensor fields and maximum likelihood attributes to indirectly reflect the fracture development situation macroscopically. However, the post-stack geometric attributes lack consideration of azimuth and offset information, and the spatial characterization of the fault system is not comprehensive enough, resulting in a low success rate of oil and gas drilling. The other is the seismic inversion technology for quantitatively predicting reservoirs. Currently, in production, it is mainly seismic inversion based on the convolution model. Due to the limitation of the wavelet, the resolution can only be equivalent to that of the seismic data, and it is affected by the model in inclined strata and complex structures, prone to problems such as poor continuity and blurred boundaries, resulting in the uncertainty of the inversion results and unable to meet the requirements of the exploration and development of fractured hydrocarbon reservoirs for inversion resolution. Summary of the Invention

[0004] To solve the defects of the existing technology, the present invention provides a five-dimensional seismic data prediction method to improve the prediction accuracy of fractured reservoirs.

[0005] On the one hand, the present invention provides a five-dimensional seismic data prediction method, including:

[0006] Performing dip imaging enhancement processing on five-dimensional seismic data;

[0007] Based on the dip imaging data, performing azimuth time window analysis and calculating the coherence attribute;

[0008] Carrying out gradient modulus attribute on the basis of the coherence attribute and performing azimuth attribute fusion;

[0009] Performing five-dimensional phase-controlled step-by-step inversion under the constraint of the fusion attribute;

[0010] Predicting the hydrocarbon reservoir of fractured volcanic rocks based on the inversion result.

[0011] In the embodiment of the present invention, the dip imaging enhancement processing of the five-dimensional seismic data includes: performing multi-scale multi-azimuth dip imaging enhancement on the original seismic data, and calculating the similarity of adjacent traces by using the changes in dip and azimuth angles.

[0012] In the embodiment of the present invention, the performing multi-scale multi-azimuth dip imaging enhancement on the original seismic data and calculating the similarity of adjacent traces by using the changes in dip and azimuth angles includes:

[0013] Defining the azimuth- and scale-dependent characteristics of seismic data for the three-dimensional seismic data volume;

[0014] Based on the seismic data with different azimuths and scales, defining a tensor matrix according to the correlation of the decomposed data;

[0015] Representing the tensor matrix as eigenvalues and eigenvectors;

[0016] Obtaining dip attributes based on eigenvalues and eigenvectors;

[0017] Calculating the similarity of adjacent traces by using the changes in dip and azimuth angles.

[0018] In the embodiment of the present invention, the performing azimuth time-window analysis based on dip imaging data and calculating coherence attributes includes: calculating coherence values along the direction consistent with the observation azimuth for seismic data with different azimuths to highlight the fracture characteristics developed perpendicular to the observation azimuth.

[0019] In the embodiment of the present invention, the performing azimuth time-window analysis based on dip imaging data and calculating coherence attributes further includes: assigning different weights to seismic traces at different positions during the calculation of coherence attributes, increasing the contribution of seismic traces closer to the central trace and reducing the corresponding weights of seismic traces farther away; calculating the weighted correlation of seismic traces in the same direction as the coherence value of the central trace in that direction.

[0020] In the embodiment of the present invention, the performing gradient modulus attribute based on the coherence attribute and conducting azimuth attribute fusion includes: using the kernel principal component analysis method for azimuth attribute fusion, calculating the kernel matrix according to the kernel function, performing eigen-decomposition on the kernel matrix, obtaining the eigenvector corresponding to the largest eigenvalue and the corresponding weight vector, and performing operations on the eigenvector and azimuth attributes to obtain the azimuth fusion attribute.

[0021] In the embodiment of the present invention, the performing five-dimensional phased step-by-step inversion under the constraint of the fusion attribute includes:

[0022] Evaluating the quality and rationality of seismic data by using the gradient modulus attribute;

[0023] Constructing an optimization operator and a forward model under the constraint of the fusion attribute;

[0024] Perform five-dimensional phased step-by-step inversion based on the optimization operator and the forward model.

[0025] In the embodiments of the present invention, an optimization operator and a forward model are constructed under the fusion attribute constraint, including:

[0026] Construct the PP-wave reflection coefficient equation characterized by the P-wave to S-wave velocity ratio in the HTI medium;

[0027] Linearize the PP-wave reflection coefficient equation and express it as a matrix-vector form of the reflectivity containing the P-wave to S-wave velocity ratio, P-wave velocity, and density;

[0028] Represent the P-wave to S-wave velocity ratio, P-wave velocity, and density reflectivity through the first-order difference operator;

[0029] Construct a forward model containing multiple inversion parameters according to the matrix-vector form.

[0030] In the embodiments of the present invention, performing five-dimensional phased step-by-step inversion based on the optimization operator and the forward model includes: constructing the objective functional for multi-parameter inversion in the HTI medium and the objective functional for anisotropic parameter update, obtaining the minimum solution of the objective functional, and inversely obtaining the P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters; using the inversely obtained P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters as the initial forward model and performing iterative loops on the initial forward model until the algorithm converges.

[0031] In the embodiments of the present invention, the objective functional for multi-parameter inversion in the HTI medium is expressed as:

[0032]

[0033] where m lfm is the phased initial model obtained by fusing azimuth attributes and consists of the P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters;

[0034] ξ1>0 is the regularization parameter to improve the noise resistance of the inversion algorithm;

[0035] C1 and C2 are coefficient matrices and are expressed as:

[0036] C1 = kron(ν, I n×n ), C2 = kron(υ, I n×n );

[0037] where I n×n is the n-order identity matrix, n represents the number of sampling points per trace of the seismic data, and kron() is the Kronecker product operator;

[0038] ν = [1, 1, 1, 0, 0, 0] T , υ = [0, 0, 0, 1, 1, 1]T ;

[0039] The specific forms of the coefficient matrices C1 and C2 are as follows:

[0040]

[0041] In the embodiment of the present invention, the PP-wave reflection coefficient equation is:

[0042]

[0043] Wherein,

[0044]

[0045] The matrix-vector form of the PP-wave reflection coefficient equation is:

[0046]

[0047] Wherein, are respectively the reflectivities of the P-to-S wave velocity ratio, the P-wave velocity, and the density;

[0048] The reflectivities of the P-to-S wave velocity ratio, the P-wave velocity, and the density represented by the first-order difference operator are respectively:

[0049]

[0050] Wherein, D is the first-order difference operator.

[0051] On the other hand, the present invention provides a five-dimensional seismic data prediction device, including:

[0052] An dip imaging processing module, configured to perform dip imaging enhancement processing on five-dimensional seismic data;

[0053] An azimuth time window analysis module, configured to perform azimuth time window analysis based on dip imaging data and calculate coherence attributes;

[0054] An attribute fusion module, configured to develop gradient modulus attributes based on coherence attributes and perform azimuth attribute fusion;

[0055] An inversion module, configured to perform five-dimensional phased step-by-step inversion under the constraint of fused attributes;

[0056] A prediction module, configured to predict the oil and gas reservoir of fractured volcanic rocks based on the inversion results.

[0057] The present invention also provides a computer device, including: a memory, a processor, and a computer program, where the computer program is stored in the memory and is configured to be executed by the processor to implement the above five-dimensional seismic data prediction method.

[0058] The present invention predicts fractured hydrocarbon reservoir formations based on five-dimensional seismic data, takes into account the mutual connection between adjacent seismic traces, combines seismic sensitive attributes and multi-trace inversion methods to accurately obtain underground medium parameters, so as to avoid the problems of strong multi-solution in single-attribute prediction and insufficient comprehensive spatial characterization of the fracture system, and at the same time improve the stability and continuity of seismic inversion, and achieve high-precision prediction of fractured reservoir formations.

[0059] Other features and advantages of the technical solution of the present invention will be described in detail in the following specific implementation section. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] The drawings described herein are used to provide a further understanding of the present invention, and constitute a part of the present invention. The schematic embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0061] Figure 1 is a flowchart of a method for predicting fractured hydrocarbon reservoir formations based on five-dimensional seismic data provided by an embodiment of the present invention;

[0062] Figure 2 is a predicted plan view of different azimuth gradient modulus attributes in a specific example of the present invention;

[0063] Figure 3 is a multi-trace inversion plan view under attribute constraints in a specific example of the present invention;

[0064] Figure 4 is a block diagram of a five-dimensional seismic data prediction device provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0065] In order to make the technical solutions and advantages in the embodiments of the present invention clearer, the following further describes the exemplary embodiments of the present invention in detail with reference to the drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than an exhaustive list of all embodiments. It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other.

[0066] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by terms such as "center", "longitudinal", "transverse", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the drawings. It is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be construed as a limitation to the present invention. In addition, the terms "first" and "second" are only used for descriptive purposes and cannot be understood as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include one or more of such features. In the description of the present invention, the meaning of "a plurality" is at least two, such as two, three, etc., unless otherwise specifically defined.

[0067] As introduced in the background art, the existing attribute analysis technology based on seismic data mainly uses attributes such as coherence and curvature to qualitatively predict medium and large-scale fractures, and uses tensor fields and maximum likelihood attributes to indirectly reflect the fracture development situation macroscopically. However, the post-stack geometric attributes lack consideration of azimuth and offset information, and the spatial characterization of the fracture system is not comprehensive enough, resulting in a low success rate of oil and gas drilling. The existing method for quantitatively predicting reservoirs using seismic inversion technology is mainly seismic inversion based on the convolution model. Due to the limitation of the wavelet, the resolution can only be equivalent to that of the seismic data, and it is affected by the model in inclined strata and complex structures, and is prone to problems such as poor continuity and blurred boundaries, resulting in uncertainty in the inversion results and unable to meet the requirements of fracture-type oil and gas reservoir exploration and development for inversion resolution.

[0068] Based on the above problems, the embodiments of the present invention provide a method for predicting fracture-type oil and gas reservoir reservoirs based on five-dimensional seismic data, which considers the mutual connection between adjacent seismic traces, combines seismic sensitive attributes and multi-trace inversion methods to accurately obtain underground medium parameters, so as to avoid the problem of strong multi-solution of single-attribute prediction and insufficient comprehensive spatial characterization of the fracture system, and at the same time improve the stability and continuity of seismic inversion, and achieve high-precision fracture-type reservoir prediction.

[0069] Figure 1 is a flowchart of the method for predicting fracture-type oil and gas reservoir reservoirs based on five-dimensional seismic data provided by the embodiments of the present invention. As Figure 1 shown, the method for predicting fracture-type oil and gas reservoir reservoirs based on five-dimensional seismic data provided by this embodiment includes the following steps:

[0070] S110, perform dip imaging enhancement processing on the five-dimensional seismic data;

[0071] S120. Conduct azimuth time window analysis based on dip imaging data and calculate coherence attributes;

[0072] S130. Develop gradient modulus attributes on the basis of coherence attributes and conduct azimuth attribute fusion;

[0073] S140. Conduct five-dimensional phased step-by-step inversion under the constraint of fused attributes;

[0074] S150. Predict oil and gas reservoir of fractured volcanic rock based on inversion results.

[0075] In the above step S110, preprocess the seismic data, specifically including multi-scale multi-azimuth dip imaging enhancement of the original seismic data, and calculate the similarity of adjacent traces using the changes in dip angle and azimuth angle to improve the lateral signal-to-noise ratio of the seismic data. First, define the azimuth- and scale-dependent characteristics ψ u,v (x) of the seismic data for the 3D seismic data volume, and its expression is:

[0076]

[0077] In Equation (1), i is a complex operator, σ is the filter bandwidth, v corresponds to the scale of the filter, u corresponds to the direction of the filter, and ||*|| represents the modulus;

[0078] k u,v = k v (cosθ, sinθ) T , k v = 2 (2-v) / 2 π, θ = ugπ / K;

[0079] For the given 3D seismic data I(x), for any given spatial point x = (in, cl, t), define the seismic data of different azimuths and scales as D u,v (x):

[0080] D u,v (x) = ψ u,v (x) * I(x) (2)

[0081] In Equation (2), D u,v (x) is the seismic data of different azimuths and scales obtained by filters determined according to different (u, v);

[0082] Based on the seismic data of different azimuths and scales, on the basis of filtering, define the tensor matrix according to the correlation of the decomposed data:

[0083]

[0084]

[0085] The tensor matrix is represented by eigenvalues and eigenvectors as follows:

[0086]

[0087] Based on the above eigenvalues and eigenvectors, the dip attribute is obtained. The similarity between adjacent traces is calculated using the changes in dip and azimuth to improve the vertical resolution of seismic data, retain its discontinuous wave group characteristics, thereby smoothing the seismic data and making discontinuities (such as faults) easier to identify.

[0088] In the above step S120, for seismic data in different azimuths, the coherence value is calculated along the direction consistent with the observation azimuth to highlight the fracture characteristics developed perpendicular to the observation azimuth. Different weights are assigned to seismic traces at different positions during the coherence attribute calculation, increasing the contribution of seismic traces closer to the central trace and decreasing the corresponding weights of seismic traces farther away. The weighted correlation of seismic traces is calculated in the same direction as the coherence value of the central trace in that direction.

[0089] In actual azimuthal seismic processing, the azimuth angle can be divided at sampling intervals of 30° or 22.5°. Generally, 6 - 8 azimuth data volumes are divided. In conventional three-dimensional analysis windows of normal size, only four directions can be taken. To achieve multi-azimuth coherence, the planar size of the window needs to be enlarged, which improves the calculation stability to a certain extent. By setting different planar sizes of the analysis window, the coherence discontinuity is enhanced, making the fracture effect more obvious. To avoid the resolution decrease caused by a large time window analysis window, resulting in some information becoming blurred, anti-distance weighting is added in the azimuth time window analysis, that is, different weights are assigned to seismic traces at different positions during the coherence attribute calculation, increasing the contribution of seismic traces closer to the central trace and decreasing the corresponding weights of seismic traces farther away, effectively reducing the average effect brought by the enlarged analysis window. The weighted correlation of seismic traces is calculated in the same direction as the coherence value of the central trace in that direction.

[0090] Suppose an azimuth analysis window containing X seismic traces (X > 5). Except for the central trace and the two adjacent traces, the remaining seismic traces use the Gaussian kernel function of formula (6) as the weight function to calculate the weight of each trace respectively, and ensure that the sum of the weights on each side is 1 through formula (7). Among them, the calculation of the weight parameter is as follows:

[0091]

[0092]

[0093] In the above formula, w x'$w_x$ is the weight of the $x$-th seismic trace in the azimuth analysis window during calculation. $X$ is the total number of seismic traces in the time window, and $\sigma$ is the weight factor. The smaller $\sigma$ is, the closer the distance from the center trace, and its weight decreases. Conversely, the larger $\sigma$ is, the relatively smoothly the weight decreases with the increase of the distance;

[0094] After obtaining the weights of each seismic trace, weighted summation is performed on both sides respectively:

[0095]

[0096]

[0097] Among them, $d_1$ and $d_2$ represent the weighted average traces on both sides, and $S$ x is the actual seismic trace, and $w$ x ' is the weight value corresponding to the seismic trace. During the large time window analysis, the weighted average of the seismic traces on both sides is used as the new seismic traces at both ends. Considering the contributions of seismic traces at different positions, the number of seismic traces in the azimuth analysis window is kept at 5, that is, 5×5 seismic traces can obtain 4 different directions.

[0098] In the above step S130, coherent attribute calculation is carried out in the dip strengthening process, and then gradient modulus attribute is carried out on the basis of coherence to perform azimuth attribute fusion. The kernel principal component analysis method (KPCA) is used to carry out azimuth attribute fusion, that is, the original data is mapped to a high-dimensional feature space $H$ through a non-linear mapping function $\varphi(c)$, and principal component analysis (PCA) is carried out in the high-dimensional space. When using KPCA to fuse multi-attributes, first calculate the kernel matrix $K$ according to the kernel function, then perform eigenvalue decomposition on the kernel matrix to obtain the eigenvector $a$ corresponding to the largest eigenvalue, and then obtain the corresponding weight vector $P$. The formula is:

[0099]

[0100] The eigenvector $a$ is operated with the azimuth attribute to obtain the azimuth fusion attribute containing most of the information. The process of azimuth coherence attribute fusion is expressed as:

[0101] $C$ multi $=\varphi(C)$ T $\alpha=\varphi(C)$ T $\varphi(C)\alpha = K\alpha=\lambda_1\alpha\ (10)$

[0102] In the formula, $C$ represents the matrix composed of the coherence of different azimuths in the same trace, $C$ multiIt is defined as multi-directional coherence. From the comparison of the two methods, both fusion methods can identify underground geological bodies. However, the KPCA method can better depict local fractures and is more conducive to identifying the planar characteristics of the lithofacies of geological bodies. Using the combined attributes in different directions to depict underground geological bodies, compared with depicting from the attribute plane of the original seismic data, the seismic attributes in different directions can better reflect the detailed changes of geological bodies. From the prediction results of the stacked attributes in different directions, the structural characteristics of geological bodies mainly develop along the fault zones on the plane, showing chaotic reflections and circular structures. This fused attribute greatly improves the prediction accuracy of geological bodies. For the description and identification of concealed geological bodies, the geological morphology is clearer and the position is more accurate.

[0103] In the above step S140, five-dimensional phase-controlled step-by-step inversion is carried out under the constraint of the fused attribute, including: evaluating the quality and rationality of seismic data using the gradient modulus attribute; constructing an optimization operator and a forward model under the constraint of the fused attribute; and performing five-dimensional phase-controlled step-by-step inversion based on the optimization operator and the forward model. Specifically, construct the PP-wave reflection coefficient equation characterized by the P-wave to S-wave velocity ratio for the HTI medium; linearize the PP-wave reflection coefficient equation, which is expressed as a matrix-vector form including the P-wave to S-wave velocity ratio, P-wave velocity, and density reflectivity; represent the P-wave to S-wave velocity ratio, P-wave velocity, and density reflectivity through the first-order difference operator; construct a forward model containing multiple inversion parameters according to the matrix-vector form. Construct the objective functional for multi-parameter inversion of the HTI medium and the objective functional for updating anisotropic parameters, obtain the minimum solution of the objective functional, and invert to obtain the P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters; use the inverted P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters as the initial forward model, and perform iterative loops until the algorithm converges.

[0104] In a specific embodiment, construct the PP-wave reflection coefficient equation characterized by the P-wave to S-wave velocity ratio for the HTI medium, and this equation is expressed as:

[0105]

[0106] where,

[0107]

[0108]

[0109] Since the P-wave to S-wave velocity ratio appears in both the coefficient term and the perturbation term, Equation (12) is non-linear. In conventional prestack seismic inversion, usually a reference value is set for the P-wave to S-wave velocity ratio in the coefficient term to linearize the equation. The matrix-vector form of the PP-wave reflection coefficient equation (11) is:

[0110]

[0111] In the formula, The reflectivities of the P-wave to S-wave velocity ratio, P-wave velocity, and density, respectively;

[0112] The reflectivities of the P-wave to S-wave velocity ratio, P-wave velocity, and density represented by the first-order difference operator are respectively:

[0113]

[0114] where D is the first-order difference operator;

[0115] Construct a forward model containing multiple inversion parameters according to the matrix-vector formula;

[0116] Since Δε = Dε, Δγ = Dγ, Δδ = Dδ, the forward model can be expressed as:

[0117]

[0118] In the formula, S is the seismic data composed of two azimuths and h incident angle gathers, W is the wavelet kernel matrix, A is the forward operator, and m is the inversion parameter.

[0119] Since there are many parameters to be inverted, this embodiment adopts a step-by-step inversion strategy to estimate the P-wave to S-wave velocity ratio. The main steps include: 1) Estimate three parameters that contribute greatly to the reflection coefficient, namely the P-wave to S-wave velocity ratio, P-wave velocity, and density; 2) Estimate the Thomsen anisotropy parameters; 3) Alternately cycle until the algorithm converges.

[0120] This embodiment of the invention designs two coefficient matrices to implement this step-by-step inversion. In step 1, a target functional for multi-parameter inversion in HTI media is constructed, which is expressed as:

[0121]

[0122] where m lfm is the phased initial model obtained by azimuth attribute fusion, which consists of the P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropy parameters;

[0123] ξ1>0 is the regularization parameter to improve the noise resistance of the inversion algorithm;

[0124] C1 and C2 are coefficient matrices used to control the parameters participating in the inversion, and they are expressed as:

[0125] C1 = kron(ν, I n×n ), C2 = kron(υ, I n×n ) (16)

[0126] where I n×n is the n-order identity matrix, n represents the number of sampling points per trace of the seismic data, and kron() is the Kronecker product operator;

[0127] ν = [1, 1, 1, 0, 0, 0] T , υ = [0, 0, 0, 1, 1, 1] T ;

[0128] The specific forms of the coefficient matrices C1 and C2 are as follows:

[0129]

[0130] In step 2, the objective functional for updating the anisotropic parameters is:

[0131]

[0132] Both equations (16) and (18) are convex optimization problems, and the minimum solutions of the objective functional can be obtained by taking partial derivatives;

[0133] In step 3, the P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters inverted in step 2 are used as the initial model, i.e., m lfm = m2, and steps 1 and 2 are executed, and the loop iteration ensures the convergence of the algorithm.

[0134] In a specific application example, the method of the present invention is applied to fracture prediction in volcanic rock oil and gas reservoirs in the Junggar Basin, and the specific process is as follows:

[0135] (1) First, preprocess the seismic data, specifically including multi-scale multi-azimuth dip imaging enhancement of the original seismic data, and calculating the similarity of adjacent traces using the changes in dip and azimuth to improve the lateral signal-to-noise ratio of the seismic data.

[0136] (2) For seismic data in different azimuths, calculate the coherence value along the direction consistent with the observation azimuth to highlight the fracture characteristics developed perpendicular to the observation azimuth. As Figure 2 shown, in actual azimuthal seismic processing, the azimuth angle is divided at a 30° sampling interval to divide 6 azimuth data volumes. Then, inverse distance weighting is added in the azimuth time window analysis to effectively reduce the averaging effect caused by the expansion of the analysis window. Calculate the weighted correlation of seismic traces in the same direction as the coherence value of the central trace in this direction.

[0137] (3) Calculate the coherence attribute in the dip enhancement process, then calculate the gradient modulus attribute on the basis of coherence, and finally perform azimuth attribute fusion using KPCA.

[0138] (4) Use the gradient modulus attribute calculated above to evaluate the quality and rationality of the seismic data, then construct an optimization operator and a low-frequency model, and then perform five-dimensional phase-controlled step-by-step inversion under attribute constraints.

[0139] The technical method of the present invention can be applied to the reservoir prediction of complex fractured volcanic rock oil and gas reservoirs, and can solve the problems of unclear lithofacies boundaries and low coincidence rate of reservoir quantitative prediction when using a single amplitude attribute to characterize the lithofacies of fractured volcanic rock reservoirs. The present invention makes full use of the characteristics of rich drilling, logging data and five-dimensional seismic data information in the study area, applies the five-dimensional seismic data prediction method based on fractured oil and gas reservoirs, and adopts the multi-attribute fusion prediction. The macroscopic characteristics of volcanic lithofacies are clear and the boundaries are clearer. Under the constraint of volcanic rock sensitive attributes, the coincidence rate of multi-channel inversion prediction of high-quality fractured volcanic rock reservoirs reaches 90% (as Figure 3 shown).

[0140] As Figure 4 shown, an embodiment of the present invention provides a five-dimensional seismic data prediction device, including: a dip imaging processing module, an azimuth time window analysis module, an attribute fusion module, an inversion module, and a prediction module. The dip imaging processing module is used to perform dip imaging enhancement processing on the five-dimensional seismic data. The azimuth time window analysis module is used to perform azimuth time window analysis based on the dip imaging data and calculate the coherence attribute. The attribute fusion module is used to carry out gradient modulus attributes on the basis of the coherence attribute and perform azimuth attribute fusion. The inversion module is used to perform five-dimensional phase-controlled step-by-step inversion under the constraint of the fusion attribute. The prediction module is used to predict the oil and gas reservoir of the fractured volcanic rock based on the inversion result.

[0141] In a specific embodiment, the dip imaging processing module performs multi-scale multi-azimuth dip imaging enhancement on the original seismic data, and calculates the similarity of adjacent traces by using the changes of dip angle and azimuth angle to improve the lateral signal-to-noise ratio of the seismic data. Specifically, it includes: defining the azimuth and scale characteristics of seismic data for the seismic three-dimensional data volume; based on the seismic data of different azimuths and scales, defining a tensor matrix according to the correlation of the decomposed data; representing the tensor matrix as eigenvalues and eigenvectors; obtaining the dip attribute based on the eigenvalues and eigenvectors; calculating the similarity of adjacent traces by using the changes of dip angle and azimuth angle. The dip imaging processing module calculates the similarity of adjacent traces by using the changes of dip angle and azimuth angle, which can improve the vertical resolution of the seismic data, retain its discontinuous wave group characteristics, so that the seismic data is smoothed and the discontinuity (such as faults) is easier to identify.

[0142] In a specific embodiment, the azimuth time window analysis module calculates the coherence value along the direction consistent with the observation azimuth for seismic data in different azimuths to highlight the fracture features developed perpendicular to the observation azimuth. When calculating the coherence attribute, different weights are assigned to seismic traces at different positions, increasing the contribution of seismic traces closer to the central trace and reducing the corresponding weights of seismic traces farther away. The weighted correlation of seismic traces is calculated in the same direction as the coherence value of the central trace in this direction. In actual azimuthal seismic processing, the azimuth angle can be divided at a sampling interval of 30° or 22.5°. Generally, 6-8 azimuth data volumes are divided. However, for a conventional three-dimensional analysis time window, only four directions can be taken. To achieve multi-azimuth coherence, the planar size of the window needs to be enlarged, which improves the calculation stability to a certain extent. By setting different planar sizes of the analysis window, the coherence discontinuity is enhanced, making the fracture effect more obvious. To avoid the resolution decrease caused by the large time window analysis window, which makes some information blurred, inverse distance weighting is added in the azimuth time window analysis. That is, different weights are assigned to seismic traces at different positions when calculating the coherence attribute, increasing the contribution of seismic traces closer to the central trace and reducing the corresponding weights of seismic traces farther away, effectively reducing the averaging effect brought by the enlarged analysis window. The weighted correlation of seismic traces is calculated in the same direction as the coherence value of the central trace in this direction.

[0143] In a specific embodiment, the attribute fusion module calculates the coherence attribute in the dip strengthening process and then calculates the gradient modulus attribute on the basis of coherence for azimuth attribute fusion. Specifically, the kernel principal component analysis method (KPCA) can be used to carry out azimuth attribute fusion, that is, the original data is mapped to a high-dimensional feature space through a non-linear mapping function, and principal component analysis (PCA) is carried out in the high-dimensional space. When using KPCA to fuse multi-attributes, the kernel matrix is first calculated according to the kernel function, and then the kernel matrix is eigen-decomposed to obtain the eigenvector a corresponding to the largest eigenvalue, and then the corresponding weight vector P is obtained. The eigenvector a is operated with the azimuth attribute to obtain the azimuth fusion attribute containing most of the information. Using the azimuthal attribute combination to depict the subsurface geological body, compared with depicting from the attribute plane of the original seismic data, the seismic attributes in different azimuths can better reflect the detailed changes of the geological body. From the results predicted by the azimuthal stacked attributes, the structural features of the geological body mainly develop along the fracture zone on the plane, showing chaotic reflections and circular structures. This fusion attribute greatly improves the prediction accuracy of the geological body. For the description and identification of hidden geological bodies, the geological morphology is clearer and the position is more accurate.

[0144] In a specific embodiment, the inversion module evaluates the quality and rationality of seismic data using the gradient modulus attribute, constructs an optimization operator and a forward model under the constraint of fusion attributes, and performs five-dimensional phased step-by-step inversion based on the optimization operator and the forward model. Specifically, it includes: constructing an equation for the PP-wave reflection coefficient characterized by the P-wave to S-wave velocity ratio in an HTI medium; linearizing the PP-wave reflection coefficient equation and expressing it as a matrix-vector form of the reflectivity containing the P-wave to S-wave velocity ratio, P-wave velocity, and density; representing the P-wave to S-wave velocity ratio, P-wave velocity, and density reflectivity through a first-order difference operator; constructing a forward model containing multiple inversion parameters according to the matrix-vector form. Constructing an objective functional for multi-parameter inversion in an HTI medium and an objective functional for anisotropic parameter update, obtaining the minimum solution of the objective functional, and inversely obtaining the P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters; using the inversely obtained P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropic parameters as the initial forward model and performing iterative loops until the algorithm converges.

[0145] The embodiment of the present invention also provides a computer device, including: a memory, a processor, and a computer program, where the computer program is stored in the memory and is configured to be executed by the processor to implement the above-mentioned method for suppressing free-surface multiples of seafloor node data.

[0146] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk memories, CD-ROMs, optical memories, etc.) containing computer-usable program codes. The solutions in the embodiments of the present invention can be implemented using various computer languages. For example, object-oriented programming languages such as Java and interpreted scripting languages such as JavaScript.

[0147] The present invention is described with reference to the flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to the embodiments of the present invention. It should be understood that each process and / or block in the flowcharts and / or block diagrams, and the combination of processes and / or blocks in the flowcharts and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, such that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the functions specified in Figure 1 one process or multiple processes and / or blocks Figure 1 one block or multiple blocks.

[0148] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to work in a particular manner, such that the instructions stored in the computer-readable memory produce a manufacture including an instruction device that implements the functions specified in the process Figure 1 one process or a plurality of processes and / or blocks Figure 1 specified in one block or a plurality of blocks.

[0149] These computer program instructions can also be loaded onto a computer or other programmable data processing apparatus, so that a series of operation steps are executed on the computer or other programmable apparatus to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable apparatus provide steps for implementing the functions specified in the process Figure 1 one process or a plurality of processes and / or blocks Figure 1 specified in one block or a plurality of blocks.

[0150] Although the preferred embodiments of the present invention have been described, additional changes and modifications can be made by those skilled in the art once they learn of the basic inventive concept. Therefore, the appended claims are intended to be construed to include the preferred embodiments as well as all changes and modifications falling within the scope of the present invention. Obviously, those skilled in the art can make various changes and variations to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalent technologies, the present invention is also intended to include these modifications and variations.

Claims

1. A five-dimensional seismic data prediction method, characterized in that, Including: Performing dip imaging enhancement processing on five-dimensional seismic data; Conducting azimuth time window analysis based on dip imaging data and calculating coherence attributes; Developing gradient modulus attributes on the basis of coherence attributes and performing azimuth attribute fusion; Performing five-dimensional phase-controlled step-by-step inversion under the constraint of fused attributes; Predicting hydrocarbon reservoir of fractured volcanic rock based on inversion results.

2. The five-dimensional seismic data prediction method according to claim 1, characterized in that The performing dip imaging enhancement processing on five-dimensional seismic data includes: Performing multi-scale multi-azimuth dip imaging enhancement on original seismic data and calculating the similarity of adjacent traces by using the changes of dip angle and azimuth angle.

3. The five-dimensional seismic data prediction method according to claim 2, characterized in that The performing multi-scale multi-azimuth dip imaging enhancement on original seismic data and calculating the similarity of adjacent traces by using the changes of dip angle and azimuth angle includes: Defining the azimuth and scale characteristics of seismic data for the three-dimensional seismic data volume; Defining a tensor matrix based on the seismic data of different azimuths and scales according to the correlation of decomposed data; Representing the tensor matrix as eigenvalues and eigenvectors; Obtaining dip attributes based on eigenvalues and eigenvectors; Calculating the similarity of adjacent traces by using the changes of dip angle and azimuth angle.

4. The five-dimensional seismic data prediction method according to claim 1, wherein The conducting azimuth time window analysis based on dip imaging data and calculating coherence attributes includes: Calculating coherence values along the direction consistent with the observation azimuth for seismic data of different azimuths to highlight the fracture characteristics developed perpendicular to the observation azimuth.

5. The five-dimensional seismic data prediction method according to claim 4, wherein The conducting azimuth time window analysis based on dip imaging data and calculating coherence attributes also includes: Assigning different weights to seismic traces at different positions during coherence attribute calculation, increasing the contribution of seismic traces closer to the central trace and decreasing the corresponding weights of seismic traces farther away; Calculating the weighted correlation of seismic traces in the same direction as the coherence value of the central trace in this direction.

6. The five-dimensional seismic data prediction method according to claim 1, characterized in that The developing gradient modulus attributes on the basis of coherence attributes and performing azimuth attribute fusion includes: Adopting the kernel principal component analysis method for azimuth attribute fusion, calculating the kernel matrix according to the kernel function, performing eigenvalue decomposition on the kernel matrix, and obtaining the eigenvector corresponding to the maximum eigenvalue and the corresponding weight vector; Performing operations on the eigenvector and azimuth attributes to obtain azimuth fusion attributes.

7. The five-dimensional seismic data prediction method according to claim 1, wherein The performing five-dimensional phase-controlled step-by-step inversion under the constraint of fused attributes includes: Evaluating the quality and rationality of seismic data by using gradient modulus attributes; Constructing an optimization operator and a forward model under the constraint of fused attributes; Performing five-dimensional phase-controlled step-by-step inversion based on the optimization operator and the forward model.

8. The five-dimensional seismic data prediction method according to claim 7, wherein, The constructing an optimization operator and a forward model under the constraint of fused attributes includes: Constructing the PP-wave reflection coefficient equation of HTI medium characterized by the P-to-S wave velocity ratio; Linearizing the PP-wave reflection coefficient equation and expressing it as a matrix-vector form containing the P-to-S wave velocity ratio, P-wave velocity, and density reflectivity; Representing the P-to-S wave velocity ratio, P-wave velocity, and density reflectivity by a first-order difference operator; Constructing a forward model containing multiple inversion parameters according to the matrix-vector form.

9. The five-dimensional seismic data prediction method according to claim 8, wherein The performing five-dimensional phase-controlled step-by-step inversion based on the optimization operator and the forward model includes: Constructing the objective functional of multi-parameter inversion of HTI medium and the objective functional of anisotropic parameter update, obtaining the minimum solution of the objective functional, and inversely obtaining the P-to-S wave velocity ratio, P-wave velocity, density, and anisotropic parameters. Taking the P-wave to S-wave velocity ratio, P-wave velocity, density, and anisotropy parameters obtained by inversion as the initial forward modeling model, the initial forward modeling model is iteratively looped until the algorithm converges.

10. The five-dimensional seismic data prediction method according to claim 9, characterized in that The objective function of the multi-parameter inversion in the HTI medium is expressed as: where m lfm is the phased initial model obtained by fusing azimuth attributes, and is composed of the P-S wave velocity ratio, P wave velocity, density, and anisotropy parameters; ξ1>0 is the regularization parameter to improve the noise resistance of the inversion algorithm; C1 and C2 are coefficient matrices, expressed as: C1 = kron(ν, I n×n ), C2 = kron(υ, I n×n ); where, I n×n is an n - order identity matrix, n represents the number of sampling points per trace of seismic data, and kron() is the Kronecker product operator; ν = [1, 1, 1, 0, 0, 0] T , υ = [0, 0, 0, 1, 1, 1] T ; The specific forms of the coefficient matrices C1 and C2 are:

11. The five-dimensional seismic data prediction method according to claim 8, wherein The PP-wave reflection coefficient equation is: where, The matrix-vector form of the PP-wave reflection coefficient equation is: wherein, are the reflectivities of the P-wave to S-wave velocity ratio, P-wave velocity, and density, respectively; The reflectivities of the P-wave to S-wave velocity ratio, P-wave velocity, and density represented by the first-order difference operator are: where D is the first-order difference operator.

12. A five-dimensional seismic data prediction device, characterized in that Including: An inclination imaging processing module for performing inclination imaging enhancement processing on five-dimensional seismic data; An azimuth time window analysis module for performing azimuth time window analysis based on the inclination imaging data to calculate the coherence attribute; An attribute fusion module for developing the gradient modulus attribute based on the coherence attribute and performing azimuth attribute fusion; An inversion module for performing five-dimensional phase-controlled step-by-step inversion under the constraint of the fused attribute; A prediction module for predicting the oil and gas reservoir of fractured volcanic rocks based on the inversion result.

13. A computer device, characterized in that, Including: A memory; A processor; and A computer program; wherein, the computer program is stored in the memory and is configured to be executed by the processor to implement the five-dimensional seismic data prediction method according to any one of claims 1-11.

Images