Seismic horizon optimization and uncertainty quantitative characterization method and device

Abstract

The invention belongs to the technical field of oil and gas field seismic exploration, and discloses a seismic horizon optimization and uncertainty quantitative characterization method and device, and the method comprises the steps: determining a target seismic horizon initial surface; based on the initial surface of the target seismic horizon, selecting a plurality of one-dimensional seismic samples to carry out network model training; based on the trained network model, obtaining corresponding embedded vectors of all the one-dimensional seismic samples in the submerged space; calculating the average distance between the one-dimensional seismic sample and the seed point in the submerged space according to the embedded vector; calculating seismic horizon probability information corresponding to the one-dimensional seismic sample according to the average distance in the submerged space; and according to the seismic horizon probability information, calculating horizon information entropy corresponding to the one-dimensional seismic sample, and determining the actual position of the target seismic horizon. By calculating the structural tensor in the seismic direction, the pickup quality and stability of the initial surface of the horizon are improved; according to the method, uncertainty description of seismic horizon probability information is realized by introducing an entropy function.

Description

Technical Field

[0001] The present invention belongs to the technical field of oil and gas field seismic exploration, and in particular relates to a seismic horizon optimization and uncertainty quantification characterization method and device. Background Art

[0002] 3D seismic horizon tracking is crucial for improving the accuracy and efficiency of seismic data interpretation. This is especially true given the current trend of rapidly increasing exploration seismic data volumes. Manual interpretation is extremely time-consuming, and in areas with complex seismic reflections, accuracy still falls short of production requirements. Therefore, developing automated seismic horizon tracking technology is crucial for improving the accuracy and efficiency of oil and gas exploration.

[0003] A large number of horizon tracking algorithms have been widely adopted and applied. This, coupled with rapid advances in computing power, has greatly reduced the workload for interpreters. Examples include horizon tracking techniques based on seismic waveform similarity comparison, seismic phase contrast, and seismic structural tensor-based horizon tracking. In particular, the global horizon tracking method based on the seismic structural tensor has significantly improved the accuracy and stability of automatic seismic horizon tracking under steep structural conditions, providing crucial technical support for oil and gas exploration in complex geological conditions.

[0004] However, in practical applications, existing methods still have certain problems that need to be solved, as follows:

[0005] 1. Under complex geological conditions, the accuracy of seismic structure tensor attribute extraction will be greatly reduced, seriously restricting the accuracy of seismic layer tracking.

[0006] 2. The global layer tracking algorithm is based on a global optimization mechanism, which has the natural disadvantage of insufficient local layer optimization.

[0007] 3. The layer tracking algorithm based on the seismic structure tensor completely abandons the guidance of seismic waveform similarity, resulting in certain errors in the seismic phase of the automatically tracked layers, that is, it cannot achieve strict geological isochronism. Summary of the Invention

[0008] To address the above problems, the present invention provides a method and device for seismic horizon optimization and uncertainty quantification, which adopts the following technical solutions:

[0009] A method for seismic layer optimization and uncertainty quantification characterization comprises the following steps: determining an initial surface of a target seismic layer based on three-dimensional seismic data, geological data and a seismic direction structure tensor of a work area; selecting multiple one-dimensional seismic samples within a set time window for network model training based on the initial surface of the target seismic layer; obtaining, based on the trained network model, corresponding embedded vectors of all one-dimensional seismic samples within the set time window in a latent space; calculating, based on the embedded vectors corresponding to the one-dimensional seismic samples in the latent space, an average distance between the one-dimensional seismic samples and a seed point in the latent space; calculating, based on the average distance between the one-dimensional seismic samples and the seed point in the latent space, seismic layer probability information corresponding to the one-dimensional seismic sample; calculating, based on the seismic layer probability information corresponding to the one-dimensional seismic sample, a layer information entropy corresponding to the one-dimensional seismic sample, and determining the actual position of the target seismic layer according to the layer information entropy.

[0010] Furthermore, the initial surface of the target seismic layer is determined based on the three-dimensional seismic data, geological data and seismic direction structure tensor of the work area, including the following steps:

[0011] Based on the 3D seismic data volume in the 3D seismic data, the target seismic layer is determined according to the geophysical information, core information, logging information and geological information of the work area;

[0012] Based on several seed points of the target seismic horizon, the seismic structure tensor matrix is ​​used to calculate the initial surface of the seismic horizon;

[0013] Based on the earthquake direction structure tensor, a new earthquake data grid space is constructed;

[0014] Calculate the local seismic slope in the new seismic line direction and the local seismic slope in the new trace direction according to the new seismic data grid space;

[0015] The initial surface of the target seismic layer is determined according to the local seismic slope in the new seismic line direction, the local seismic slope in the new trace direction and the initial surface of the seismic layer calculated using the seismic structure tensor matrix.

[0016] Furthermore, based on the earthquake directional structure tensor, a new earthquake data grid space is constructed, which includes the following steps:

[0017] The seismic direction structure tensor matrix is ​​constructed, and the eigenvector corresponding to the maximum eigenvalue in the seismic direction structure tensor matrix is ​​extracted. The eigenvector corresponding to the maximum eigenvalue is mapped to the seismic data grid space to obtain a new seismic data grid space.

[0018] Furthermore, based on the initial surface of the target seismic horizon, multiple one-dimensional seismic samples are selected within a set time window for network model training, including the following steps:

[0019] In the three-dimensional seismic data volume, with the initial surface of the target seismic layer as the center, multiple one-dimensional seismic samples are randomly selected within the set time window length and input into the vector variational autoencoder model for network training.

[0020] Furthermore, based on the trained network model, the corresponding embedded vectors of all one-dimensional earthquake samples within the set time window in the latent space are obtained, including the following steps:

[0021] The encoder of the trained vector variational autoencoder model is extracted. In the three-dimensional seismic data volume, all one-dimensional seismic samples are extracted within a set time window with the initial surface of the target seismic layer as the center, and the corresponding embedded vectors are calculated.

[0022] Furthermore, according to the seismic layer probability information corresponding to the one-dimensional seismic sample, the layer information entropy corresponding to the one-dimensional seismic sample is calculated, and the actual position of the target seismic layer is determined according to the layer information entropy, including the following steps:

[0023] According to the normalized seismic layer probability information, the layer information entropy corresponding to the one-dimensional seismic sample is calculated;

[0024] For one-dimensional earthquake samples whose horizon information entropy is greater than the set value, the time window range is narrowed to recalculate the earthquake horizon probability information, and the position with the maximum probability is determined as the actual position of the target earthquake horizon.

[0025] Furthermore, the average distance between the one-dimensional seismic sample and the seed point in the latent space is calculated based on the embedded vector corresponding to the one-dimensional seismic sample in the latent space, as follows:

[0026]

[0027] Where d(γ) represents the average distance between the one-dimensional earthquake sample at coordinate γ and the earthquake labels at all seed points in the latent space, ||.||2 represents the Euclidean distance, N represents the number of seed points, i=1:N, F all (γ) represents the embedded vector corresponding to the one-dimensional earthquake sample at coordinate γ, F i Represents the embedded vector corresponding to the one-dimensional earthquake sample at the seed point.

[0028] Furthermore, the seismic layer probability information corresponding to the one-dimensional seismic sample is calculated based on the average distance between the one-dimensional seismic sample and the seed point in the latent space, as follows:

[0029] p(γ)=e -d(γ)

[0030] Where p(γ) represents the probability information of the earthquake layer corresponding to the one-dimensional earthquake sample at the coordinate γ, and d(γ) represents the average distance between the one-dimensional earthquake sample at the coordinate γ and the earthquake labels at all seed points in the latent space.

[0031] Furthermore, based on the normalized seismic horizon probability information, the horizon information entropy corresponding to the one-dimensional seismic sample is calculated as follows:

[0032]

[0033] Where Epy(γ) represents the horizon information entropy corresponding to the one-dimensional seismic sample at coordinate γ, p j (γ) represents the normalized seismic layer probability information of the one-dimensional seismic sample at the coordinate γ, M represents the number of one-dimensional seismic samples within the vertical time window, j = 1, ..., M.

[0034] The present invention also provides a seismic horizon optimization and uncertainty quantification characterization device, comprising:

[0035] The first calculation module is used to determine the initial surface of the target seismic layer based on the three-dimensional seismic data, geological data and seismic direction structure tensor of the work area;

[0036] The model training module is used to select multiple one-dimensional earthquake samples within a set time window for network model training based on the initial surface of the target earthquake layer;

[0037] The data acquisition module is used to obtain the corresponding embedded vectors of all one-dimensional seismic samples in the latent space within a set time window based on the trained network model;

[0038] The second calculation module is used to calculate the average distance between the one-dimensional seismic sample and the seed point in the latent space according to the embedded vector corresponding to the one-dimensional seismic sample in the latent space;

[0039] The third calculation module is used to calculate the seismic layer probability information corresponding to the one-dimensional seismic sample based on the average distance between the one-dimensional seismic sample and the seed point in the latent space;

[0040] The fourth calculation module is used to calculate the horizon information entropy corresponding to the one-dimensional seismic sample according to the seismic horizon probability information corresponding to the one-dimensional seismic sample, and determine the actual position of the target seismic horizon according to the horizon information entropy.

[0041] Furthermore, the first calculation module is specifically configured to:

[0042] Based on the 3D seismic data volume in the 3D seismic data, the target seismic layer is determined according to the geophysical information, core information, logging information and geological information of the work area;

[0043] Based on several seed points of the target seismic horizon, the seismic structure tensor matrix is ​​used to calculate the initial surface of the seismic horizon;

[0044] Based on the earthquake direction structure tensor, a new earthquake data grid space is constructed;

[0045] Calculate the local seismic slope in the new seismic line direction and the local seismic slope in the new trace direction according to the new seismic data grid space;

[0046] The initial surface of the target seismic layer is determined according to the local seismic slope in the new seismic line direction, the local seismic slope in the new trace direction and the initial surface of the seismic layer calculated using the seismic structure tensor matrix.

[0047] Furthermore, the fourth calculation module is specifically configured to:

[0048] According to the normalized seismic layer probability information, the layer information entropy corresponding to the one-dimensional seismic sample is calculated;

[0049] For one-dimensional earthquake samples whose horizon information entropy is greater than the set value, the time window range is narrowed to recalculate the earthquake horizon probability information, and the position with the maximum probability is determined as the actual position of the target earthquake horizon.

[0050] Beneficial effects of the present invention:

[0051] 1. By calculating the seismic directional structure tensor, the present invention effectively further improves the calculation accuracy of the local seismic slope on the basis of the traditional seismic structure tensor, and greatly improves the picking quality and stability of the initial plane of the horizon;

[0052] 2. The present invention considers the quantitative characterization of the uncertainty of the automatic tracking results of seismic horizons, and realizes the uncertainty description of seismic horizon probability information by introducing the entropy function, thereby achieving further optimization of the tracking results on this basis;

[0053] 3. The present invention extracts deep learning features of one-dimensional seismic samples through a vector variational autoencoder, which can minimize the information redundancy between samples and reduce the uncertainty of the seismic horizon probability obtained.

[0054] 4. This invention considers the impact of one-dimensional seismic sample length on horizon tracing results. By appropriately adjusting the sample length, the number of channels in the deep learning network is increased. Through these adjustments and optimizations, deep learning feature extraction of seismic data is improved, further enhancing the accuracy of horizon tracing.

[0055] 5. The algorithm of the present invention is simple and easy to implement, has high operating efficiency, and is easy to promote and apply.

[0056] Other features and advantages of the present invention will be described in the following description, and in part will become apparent from the description, or will be understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures pointed out in the description and the drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0058] Figure 1 A schematic flow chart of a method for seismic layer optimization and uncertainty quantification characterization according to an embodiment of the present invention is shown;

[0059] Figure 2 A schematic diagram of a deep learning network based on an improved one-dimensional vector variational autoencoder according to an embodiment of the present invention is shown;

[0060] Figure 3 A schematic diagram showing a two-dimensional seismic profile and its corresponding target horizon seed points and initial horizon according to an embodiment of the present invention is shown;

[0061] Figure 4 Shows the use of Figure 2 Schematic diagram of the changes in reconstruction error during deep learning network training;

[0062] Figure 5 A schematic diagram of latent space distance and seismic horizon probability obtained by combining three model samples according to an embodiment of the present invention is shown;

[0063] Figure 6 Shown according to Figure 5 The information entropy values ​​corresponding to the three earthquake probabilities;

[0064] Figure 7 Shown in Figure 3 Schematic diagram of seismic layer tracing results obtained based on the VQ-VAE-C8+L32 combination and the VQ-VAE-C16+L128 combination on the basis of 2D seismic profiles;

[0065] Figure 8 Shown Figure 7 The entropy value corresponding to the seismic layer tracing result;

[0066] Figure 9 FIG2 shows a schematic diagram of a three-dimensional seismic data volume used according to an embodiment of the present invention;

[0067] Figure 10 A schematic diagram of an initial surface of a three-dimensional seismic horizon used in an embodiment of the present invention is shown;

[0068] Figure 11 The VQ-VAE-C8+L32 Figure 10 Schematic diagram of the results after layer tracking optimization of the initial surface of the middle layer;

[0069] Figure 12 The VQ-VAE-C16+L128 pair is shown. Figure 10 Schematic diagram of the results after layer tracking optimization of the initial surface of the middle layer;

[0070] Figure 13 The VQ-VAE-C8+L32 Figure 10 The entropy value corresponding to the result of layer tracking optimization of the initial surface of the middle layer;

[0071] Figure 14 The VQ-VAE-C16+L128 pair is shown. Figure 10 The entropy value corresponding to the result of layer tracking optimization of the initial surface of the middle layer;

[0072] Figure 15 Shown in Figure 14 Under entropy constraints Figure 12 3D view of the results after further optimization of the horizon;

[0073] Figure 16 Shown Figure 10 and Figure 15 Comparison chart of 3D seismic horizon results;

[0074] Figure 17 A structural schematic diagram of a seismic horizon optimization and uncertainty quantification characterization device according to an embodiment of the present invention is shown. DETAILED DESCRIPTION

[0075] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0076] It should be noted that the terms "first," "second," etc., in this application are used to distinguish similar objects, and are not necessarily used to describe a specific order or precedence. It should be understood that the terms used in this manner are interchangeable where appropriate, so as to facilitate the embodiments of the present application described herein.

[0077] To address the above problems, the present invention proposes a method and device for seismic layer optimization and uncertainty quantification, which carry out intelligent optimization of the initial surface of the layer and uncertainty quantification, can greatly highlight the characteristics of the target seismic layer, thereby achieving high-precision tracking of the seismic layer and improving the seismic phase consistency of the layer; in addition, by introducing the concept of information entropy to quantitatively characterize the uncertainty of earthquake probability, the phase consistency and comprehensive accuracy of automatic tracking of seismic layers are comprehensively improved.

[0078] like Figure 1 As shown, a seismic layer optimization and uncertainty quantification characterization method includes the following steps:

[0079] S1. Determine the initial surface of the target seismic horizon based on the 3D seismic data, geological data and seismic direction structure tensor of the work area, as follows:

[0080] S11. Based on the 3D seismic data volume s(t,x,y) in the 3D seismic data, the target seismic layer is determined according to the geophysical information, core information, logging information and geological information of the work area.

[0081] S12, manually extract several seed points of the target seismic horizon, denoted as seed(t i ,x i ,y i )(i=1:N), where (t i ,x i ,y i )(i=1:N) represent the two-way travel time and plane coordinates corresponding to these N seed points in the three-dimensional seismic data volume.

[0082] S13. Based on several seed points of the target seismic horizon, the seismic structure tensor matrix is ​​used to calculate the initial surface of the seismic horizon. Assume that the three-dimensional seismic structure tensor T can be expressed as:

[0083]

[0084] Among them, g1, g2, g3 are the gradients of seismic data along the vertical, line and trace directions respectively; <.> represents filtering processing. By performing eigendecomposition on the structure tensor matrix T, we can obtain:

[0085] T=λ u UU T +λ v VV T +λ w WW T (2)

[0086] Where U, V, and W are the eigenvectors of the structure tensor matrix T in the vertical, line, and track directions, respectively, and the corresponding eigenvalues ​​are λu ,λ v ,λ w And satisfy λ u ≥λ v ≥λ w Assume that the vector U = [u1, u2, u3] is downward (u1>0), then the local seismic slope along the seismic line and the trace direction (respectively denoted as s l and s c ) can be expressed as:

[0087]

[0088] The initial surface of the seismic layer h(x,y) can be obtained by iteratively solving the following formula:

[0089]

[0090] Where h i (x,y) represents the seismic layer plane after the i-th iteration; h 0 (x, y) is the initial value formed by interpolation of the above N seed points; α is a small constant used to constrain the smoothness of the layer surface.

[0091] The above solution also has certain issues: when underground geological conditions are complex, the calculation accuracy of the seismic structure tensor is generally low, which in turn affects the accuracy of the horizon h(x,y). To address this issue, the present invention introduces the seismic directional structure tensor to perform the initial horizon plane calculation. Compared to the traditional structure tensor, the seismic directional structure tensor considers the difference between seismic data gradient and waveform gradient, further improving the calculation accuracy of seismic gradient.

[0092] S14. Constructing a new seismic data grid space based on the seismic direction structure tensor, specifically: constructing a seismic direction structure tensor matrix, extracting the eigenvector corresponding to the maximum eigenvalue in the seismic direction structure tensor matrix, and mapping the eigenvector corresponding to the maximum eigenvalue to the seismic data grid space to obtain a new seismic data grid space, including the following steps:

[0093] To calculate the earthquake directional structure tensor, we first need to construct a UPQ space, where U is the eigenvector in equation (2), and P and Q can be expressed as:

[0094]

[0095] In the formula, the three unit vectors U, P, and Q are perpendicular to each other. On this basis, the new seismic gradient is recalculated in the UPQ space:

[0096]

[0097] Where γ represents the spatial coordinates of 3D seismic data (t, x, y); g u 、g p 、g q They are the seismic gradients in the three vector directions of U, P, and Q. Since U, P, and Q are space-varying, the calculated seismic gradients have higher accuracy and stability in complex geological conditions. npq It can be expressed as:

[0098]

[0099] T npq Carry out eigendecomposition to extract the eigenvector U corresponding to the maximum eigenvalue upq =[u u ,u p ,u q ]; then U upq Mapping from the UPQ space to the traditional seismic data grid space, we obtain a new seismic data grid space, denoted as U':

[0100]

[0101] S15. Calculate the local seismic slope in the new seismic line direction and the local seismic slope in the new trace direction based on the new seismic data grid space, as follows:

[0102] The local seismic slope s′ of the new seismic line direction can be calculated with higher accuracy through U′=[u′1,u′2,u′3] l and the local seismic slope s' of the new trace direction c :

[0103]

[0104] S16. Determine the initial surface of the target seismic horizon based on the local seismic slope in the new seismic line direction, the local seismic slope in the new trace direction, and the initial surface of the seismic horizon calculated using the seismic structure tensor matrix, as follows:

[0105] s′ l and s' c By substituting it into formula (4), we can obtain the initial layer surface with higher accuracy.

[0106] In step S1, a higher-precision initial layer surface is extracted based on the seismic directional structure tensor; compared with the traditional seismic structure tensor, the seismic directional structure tensor takes into account the difference between the seismic data gradient and the waveform gradient, and can more accurately extract seismic slope information in complex seismic reflection areas, thereby improving the accuracy of picking the initial surface of the seismic layer.

[0107] S2. Based on the initial surface of the target seismic horizon, multiple one-dimensional seismic samples are selected within the set time window for network model training, as follows:

[0108] In the three-dimensional seismic data volume, with the initial surface of the target seismic layer as the center, multiple one-dimensional seismic samples are randomly selected within the set time window length and input into the Vector Quantised Variational AutoEncoder (VQ-VAE) model for network training.

[0109] In this step, after obtaining the initial seismic horizon h(x,y) in step S1, multiple one-dimensional seismic samples ψ can be randomly sampled within a certain time window length with h(x,y) as the center in the three-dimensional seismic data volume to carry out vector quantised variational autoencoder (VQ-VAE, vector quantised variational autoencoder, discrete feature learning model) network training, such as Figure 3 As shown, Figure 3 The two-dimensional seismic profile used in the present invention and its corresponding target horizon seed points (circles) and initial horizons (dashed lines) are shown.

[0110] like Figure 2 As shown, the present invention adopts a deep learning network of an improved one-dimensional vector variational autoencoder as a self-supervised network. The VQ-VAE model training mechanism is as follows:

[0111]

[0112] Among them, Encoder VQ-VAE (.) represents the encoder of the VQ-VAE model, Decoder VQ-VAE (.) is the decoder; F represents the embedded vector corresponding to the one-dimensional earthquake sample in the latent space constructed by the VQ-VAE model; ψ' represents the vector based on F and the decoder Decoder VQ-VAE (.) is the approximate value of the reconstructed one-dimensional earthquake label ψ; the smaller the difference between ψ and ψ', the higher the model training accuracy. Generally, with more training iterations, the difference between ψ and ψ' gradually decreases. When the model is trained to a certain level, the error tends to stabilize and is less than a specific threshold, and the model training is considered complete.

[0113] like Figure 4 As shown, Figure 4 It is adopted Figure 2The changes in reconstruction error (Reconstruction loss) during the deep learning network training process in , where VQ-VAE-C8+L32 means that the training adopts an 8-channel VQ-VAE model and a one-dimensional seismic sample with a length of 32 samples; VQ-VAE-C8+L128 means that the training adopts an 8-channel VQ-VAE model and a one-dimensional seismic sample with a length of 128 samples; VQ-VAE-C16+L128 means that the training adopts a 16-channel VQ-VAE model and a one-dimensional seismic sample with a length of 128 samples. Figure 4 It can be seen that the reconstruction error of the model sample combination of VQ-VAE-C16+L128 is significantly reduced compared with other combinations, which has obvious advantages.

[0114] Figure 5 Figures 2 and 3 show the latent space distances (left) and seismic layer probabilities (right) obtained for the three aforementioned combinations. (a, b) correspond to the results of VQ-VAE-C8+L32; (c, d) to VQ-VAE-C8+L128; and (e, f) to VQ-VAE-C16+L128. The figure shows that the VQ-VAE-C16+L128 model sample combination has a clear advantage over the other combinations in that the target layer is more prominent in the latent space, highlighting the target layer events.

[0115] S3. Based on the trained network model, the corresponding embedded vectors of all one-dimensional earthquake samples in the latent space within the set time window are obtained as follows:

[0116] Extract the encoder of the trained vector variational autoencoder model VQ-VAE (.); Then, in the 3D seismic data volume, with the initial surface of the target seismic layer h(x,y) as the center, all 1D seismic samples ψ are extracted within the set time window range. all And calculate the corresponding embedded vector F all , specifically as follows:

[0117] F all =Encoder VQ-VAE (ψ all ) (11)

[0118] After obtaining the embedded vector F corresponding to all one-dimensional seismic samples within the time window, all After that, we can further calculate the relationship between each embedded vector and the seed point seed(t i ,x i ,y i )(i=1:N) is the spatial distance of the embedded vector corresponding to the position.

[0119] The present invention introduces a vector variational autoencoder (VQ-VAE) network model to carry out intelligent optimization and uncertainty quantification characterization of the above-mentioned initial surface of the layer; the vector variational autoencoder can extract the embedded vector features of the one-dimensional seismic label in the latent space. By mapping the original seismic data into the latent space to become an embedded vector, it can greatly highlight the characteristics of the target seismic layer, fully expand the spatial distance between the seismic reflection at the location of the seismic layer and other seismic reflections, thereby achieving high-precision tracking of the seismic layer and improving the seismic phase consistency of the layer.

[0120] S4. Calculate the average distance between the one-dimensional seismic sample and the seed point in the latent space based on the embedded vector corresponding to the one-dimensional seismic sample in the latent space, as follows:

[0121] Let seed point seed(t i ,x i ,y i )(i=1,2,...,N) one-dimensional earthquake samples ψ i The corresponding embedded vector is F i ,Right now:

[0122] F i =Encoder VQ-VAE (ψ i ) (i=1,...,N) (12)

[0123] Then, in the 3D seismic data volume, with h(x,y) as the center, the average distance between any sample and the seed point in the latent space within the set time window length can be expressed as:

[0124]

[0125] Where d(γ) represents the average distance between the one-dimensional earthquake sample at coordinate γ and the earthquake labels at all seed points in the latent space, ||.||2 represents the Euclidean distance, N represents the number of seed points, i=1:N, F all (γ) represents the embedded vector corresponding to the one-dimensional earthquake sample at coordinate γ, F i Represents the embedded vector corresponding to the one-dimensional earthquake sample at the seed point.

[0126] S5. Calculate the seismic layer probability information p(γ) corresponding to the one-dimensional seismic sample based on the average distance between the one-dimensional seismic sample and the seed point in the latent space, as follows:

[0127] p(γ)=e -d(γ) (14)

[0128] Where p(γ) represents the seismic layer probability information corresponding to the one-dimensional earthquake sample at coordinate γ.

[0129] Then, for any seismic data, assume that there are M one-dimensional seismic samples within the vertical time window, and perform normalization on the seismic layer probability information p(γ) corresponding to each data:

[0130] p(γ j )=p(γ j ) / ∑p(γ) (j=1,...,M) (15)

[0131] S6. Calculate the horizon information entropy corresponding to the one-dimensional seismic sample based on the seismic horizon probability information corresponding to the one-dimensional seismic sample, and determine the actual position of the target seismic horizon based on the horizon information entropy, as follows:

[0132] S61, according to the normalized seismic layer probability information, calculate the layer information entropy corresponding to the one-dimensional seismic sample, for any seismic data, after obtaining the normalized seismic probability p(γ j )(j=1,...,M), calculate the horizon information entropy Epy(γ) corresponding to the seismic data, as follows:

[0133]

[0134] Where Epy(γ) represents the horizon information entropy corresponding to the one-dimensional seismic sample at coordinate γ, p j (γ) represents the normalized seismic layer probability information of the one-dimensional seismic sample at the coordinate γ, and M represents the number of one-dimensional seismic samples within the vertical time window.

[0135] S62: For one-dimensional earthquake samples where the horizon information entropy is greater than a set value, the time window range is narrowed and the probability information of the earthquake horizon is recalculated, and the position with the maximum probability is determined as the actual position of the target earthquake horizon, as follows:

[0136] The smaller the entropy value, the greater the certainty of the seismic layer calculated on the seismic trace. For seismic data traces with large entropy values, the time window range can be further narrowed to recalculate the seismic layer probability information; finally, based on this, the position with the maximum probability is selected as the actual position of the target seismic layer:

[0137] h(x,y)=t j ; satisfy p(t j ,x,y)=max(p(t1:t M ,x,y))(17)

[0138] The seismic traces with entropy values ​​greater than a given threshold are recalculated to obtain the final target layer tracking results.

[0139] Figure 6 is based on Figure 5 The information entropy values ​​corresponding to the three earthquake probabilities in . Figure 6 It can be seen that the model sample combination of VQ-VAE-C16+L128 has a smaller information entropy than other combinations, and the uncertainty of the results is lower, which has obvious advantages.

[0140] Figure 7 The present invention Figure 3 Seismic horizon tracing results obtained based on the VQ-VAE-C8+L32 combination and the VQ-VAE-C16+L128 combination on the basis of 2D seismic profiles; Figure 7 On the left are the seed points (dots) and initial values ​​(dashed lines) corresponding to the above layers. Figure 7 The right side shows the tracking results (the corresponding results of VQ-VAE-C8+L32 are dark curves, and the corresponding results of VQ-VAE-C16+L128 are light curves). Figure 7 It can be seen that the model sample combination of VQ-VAE-C16+L128 has higher accuracy and better stability than other combinations, which shows obvious advantages.

[0141] Figure 8 yes Figure 7 The entropy value corresponding to the result. Figure 8 It can be seen that the information entropy calculated by the model sample combination of VQ-VAE-C16+L128 is smaller than that of other combinations, so the uncertainty of the result is lower, which has obvious advantages.

[0142] Figure 9 The three-dimensional seismic data volume used in the application example of the present invention is Figure 9 The target layer is pointed by the arrow in the middle. Figure 9 It can be seen that faults are developed near this layer and the fault throw is relatively large.

[0143] Figure 10 The initial surface of the three-dimensional seismic horizon used in the application example of the present invention is: Figure 10 The middle dot is the seed point corresponding to the horizon; the initial horizon is the automatic tracking result based on the seismic structure tensor and the above three seed points. Figure 10 It can be seen that due to the well-developed faults near this layer and the relatively large fault throw, the above initial layer surface has serious layer penetration and phase errors.

[0144] Figure 11 and Figure 12 To use VQ-VAE-C8+L32(a) and VQ-VAE-C16+L128(b) respectively Figure 10 The results of layer tracking optimization of the initial surface of the middle layer; the entropy values ​​corresponding to the results are respectively Figure 13 and Figure 14 .Depend on Figure 11-14It can be seen that the layer tracking results obtained by VQ-VAE-C16+L128 are more accurate and have a smaller entropy value than those obtained by VQ-VAE-C8+L32, which can effectively overcome the influence of large fault distance on the layer tracking results, verifying the advantages of this technical method.

[0145] Figure 15 For Figure 14 Under entropy constraints Figure 12 The resulting 3D view of the layers after further optimization. The accuracy of layer tracking was further improved by setting a specific threshold (0.86 in this case) and then re-tracking the layers with entropy values ​​greater than the threshold.

[0146] Figure 16 for Figure 10 and Figure 15 From the plane display of the results comparison, it can be seen that after the application of the technology of the present invention, the accuracy of the tracking results from the initial layer plane to the final layer plane is greatly improved.

[0147] The present invention introduces the concept of information entropy to quantitatively characterize the uncertainty of earthquake probability, and realizes secondary optimization by re-adjusting the window range of the area with higher uncertainty, thereby comprehensively improving the phase consistency and comprehensive accuracy of automatic tracking of seismic horizons.

[0148] Based on the above seismic layer optimization and uncertainty quantification characterization methods, such as Figure 17 As shown, the present invention also provides a seismic layer optimization and uncertainty quantification characterization device, including a first calculation module, a model training module, a data acquisition module, a second calculation module, a third calculation module and a fourth calculation module.

[0149] Among them, the first calculation module is used to determine the initial surface of the target seismic layer based on the three-dimensional seismic data, geological data and seismic direction structure tensor of the work area; the model training module is used to select multiple one-dimensional seismic samples within the set time window for network model training based on the initial surface of the target seismic layer; the data acquisition module is used to obtain the embedded vectors corresponding to all one-dimensional seismic samples within the set time window in the latent space based on the trained network model; the second calculation module is used to calculate the average distance between the one-dimensional seismic sample and the seed point in the latent space based on the embedded vector corresponding to the one-dimensional seismic sample in the latent space; the third calculation module is used to calculate the seismic layer probability information corresponding to the one-dimensional seismic sample based on the average distance between the one-dimensional seismic sample and the seed point in the latent space; the fourth calculation module is used to calculate the layer information entropy corresponding to the one-dimensional seismic sample based on the seismic layer probability information corresponding to the one-dimensional seismic sample, and determine the actual position of the target seismic layer based on the layer information entropy.

[0150] This invention significantly improves the accuracy of automatic 3D layer tracking and enables quantitative evaluation of the uncertainty of tracking results, achieving excellent results in multiple real-world data applications. Even with a small number of artificial seed points or large fault distances, the invention can achieve stable tracking results, demonstrating satisfactory results in multiple field tests at home and abroad.

[0151] Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for seismic horizon optimization and uncertainty quantification, characterized in that: The following steps are involved: Determine the initial surface of the target seismic layer based on the three-dimensional seismic data, geological data and seismic direction structure tensor of the work area; Based on the initial surface of the target seismic horizon, multiple one-dimensional seismic samples are selected within the set time window for network model training; Based on the trained network model, the corresponding embedded vectors of all one-dimensional earthquake samples within the set time window are obtained in the latent space; According to the embedded vector corresponding to the one-dimensional earthquake sample in the latent space, the average distance between the one-dimensional earthquake sample and the seed point in the latent space is calculated; According to the average distance between the one-dimensional seismic sample and the seed point in the latent space, the seismic layer probability information corresponding to the one-dimensional seismic sample is calculated; According to the seismic layer probability information corresponding to the one-dimensional seismic sample, the layer information entropy corresponding to the one-dimensional seismic sample is calculated, and the actual position of the target seismic layer is determined according to the layer information entropy.

2. The seismic horizon optimization and uncertainty quantification method according to claim 1, wherein: Determine the initial surface of the target seismic horizon based on the 3D seismic data, geological data and seismic directional structure tensor of the work area, including the following steps: Based on the 3D seismic data volume in the 3D seismic data, the target seismic layer is determined according to the geophysical information, core information, logging information and geological information of the work area; Based on several seed points of the target seismic horizon, the seismic structure tensor matrix is ​​used to calculate the initial surface of the seismic horizon; Based on the earthquake direction structure tensor, a new earthquake data grid space is constructed; Calculate the local seismic slope in the new seismic line direction and the local seismic slope in the new trace direction according to the new seismic data grid space; The initial surface of the target seismic layer is determined according to the local seismic slope in the new seismic line direction, the local seismic slope in the new trace direction and the initial surface of the seismic layer calculated using the seismic structure tensor matrix.

3. The seismic horizon optimization and uncertainty quantification characterization method according to claim 2, characterized in that: Based on the seismic directional structure tensor, a new seismic data grid space is constructed, which includes the following steps: The seismic direction structure tensor matrix is ​​constructed, and the eigenvector corresponding to the maximum eigenvalue in the seismic direction structure tensor matrix is ​​extracted. The eigenvector corresponding to the maximum eigenvalue is mapped to the seismic data grid space to obtain a new seismic data grid space.

4. The seismic horizon optimization and uncertainty quantification characterization method according to claim 1, characterized in that: Based on the initial surface of the target seismic horizon, multiple one-dimensional seismic samples are selected within the set time window for network model training, including the following steps: In the three-dimensional seismic data volume, with the initial surface of the target seismic layer as the center, multiple one-dimensional seismic samples are randomly selected within the set time window length and input into the vector variational autoencoder model for network training.

5. The seismic horizon optimization and uncertainty quantification method according to claim 4, characterized in that: Based on the trained network model, the corresponding embedded vectors of all one-dimensional earthquake samples within the set time window in the latent space are obtained, including the following steps: The encoder of the trained vector variational autoencoder model is extracted. In the three-dimensional seismic data volume, all one-dimensional seismic samples are extracted within a set time window with the initial surface of the target seismic layer as the center, and the corresponding embedded vectors are calculated.

6. The method for seismic horizon optimization and uncertainty quantification according to any one of claims 1 to 5, characterized in that: According to the seismic layer probability information corresponding to the one-dimensional seismic sample, the layer information entropy corresponding to the one-dimensional seismic sample is calculated, and the actual position of the target seismic layer is determined according to the layer information entropy, including the following steps: According to the normalized seismic layer probability information, the layer information entropy corresponding to the one-dimensional seismic sample is calculated; For one-dimensional earthquake samples whose horizon information entropy is greater than the set value, the time window range is narrowed to recalculate the earthquake horizon probability information, and the position with the maximum probability is determined as the actual position of the target earthquake horizon.

7. The seismic horizon optimization and uncertainty quantification method according to claim 1, characterized in that: According to the embedded vector corresponding to the one-dimensional earthquake sample in the latent space, the average distance between the one-dimensional earthquake sample and the seed point in the latent space is calculated as follows: Where d(γ) represents the average distance between the one-dimensional earthquake sample at coordinate γ and the earthquake labels at all seed points in the latent space, ||.||2 represents the Euclidean distance, N represents the number of seed points, i=1:N, F all (γ) represents the embedded vector corresponding to the one-dimensional earthquake sample at coordinate γ, F i Represents the embedded vector corresponding to the one-dimensional earthquake sample at the seed point.

8. The seismic horizon optimization and uncertainty quantification method according to claim 1, characterized in that: According to the average distance between the one-dimensional seismic sample and the seed point in the latent space, the seismic layer probability information corresponding to the one-dimensional seismic sample is calculated as follows: p(γ)=e -d(γ) Where p(γ) represents the probability information of the earthquake layer corresponding to the one-dimensional earthquake sample at the coordinate γ, and d(γ) represents the average distance between the one-dimensional earthquake sample at the coordinate γ and the earthquake labels at all seed points in the latent space.

9. The seismic horizon optimization and uncertainty quantification method according to claim 6, characterized in that: According to the normalized seismic horizon probability information, the horizon information entropy corresponding to the one-dimensional seismic sample is calculated as follows: Where Epy(γ) represents the horizon information entropy corresponding to the one-dimensional seismic sample at coordinate Υ, p j (Y) represents the normalized seismic layer probability information of the one-dimensional seismic sample at the coordinate Y, M represents the number of one-dimensional seismic samples within the vertical time window, j = 1, ..., M.

10. A seismic horizon optimization and uncertainty quantification device, characterized in that: include: The first calculation module is used to determine the initial surface of the target seismic layer based on the three-dimensional seismic data, geological data and seismic direction structure tensor of the work area; The model training module is used to select multiple one-dimensional earthquake samples within a set time window for network model training based on the initial surface of the target earthquake layer; The data acquisition module is used to obtain the corresponding embedded vectors of all one-dimensional seismic samples in the latent space within a set time window based on the trained network model; The second calculation module is used to calculate the average distance between the one-dimensional seismic sample and the seed point in the latent space according to the embedded vector corresponding to the one-dimensional seismic sample in the latent space; The third calculation module is used to calculate the seismic layer probability information corresponding to the one-dimensional seismic sample based on the average distance between the one-dimensional seismic sample and the seed point in the latent space; The fourth calculation module is used to calculate the horizon information entropy corresponding to the one-dimensional seismic sample according to the seismic horizon probability information corresponding to the one-dimensional seismic sample, and determine the actual position of the target seismic horizon according to the horizon information entropy.

11. The seismic horizon optimization and uncertainty quantification characterization device according to claim 10, characterized in that: The first calculation module is specifically used for: Based on the 3D seismic data volume in the 3D seismic data, the target seismic layer is determined according to the geophysical information, core information, logging information and geological information of the work area; Based on several seed points of the target seismic horizon, the seismic structure tensor matrix is ​​used to calculate the initial surface of the seismic horizon; Based on the earthquake direction structure tensor, a new earthquake data grid space is constructed; Calculate the local seismic slope in the new seismic line direction and the local seismic slope in the new trace direction according to the new seismic data grid space; The initial surface of the target seismic layer is determined according to the local seismic slope in the new seismic line direction, the local seismic slope in the new trace direction and the initial surface of the seismic layer calculated using the seismic structure tensor matrix.

12. The seismic horizon optimization and uncertainty quantification device according to claim 10 or 11, characterized in that: The fourth calculation module is specifically used for: According to the normalized seismic layer probability information, the layer information entropy corresponding to the one-dimensional seismic sample is calculated; For one-dimensional earthquake samples whose horizon information entropy is greater than the set value, the time window range is narrowed to recalculate the earthquake horizon probability information, and the position with the maximum probability is determined as the actual position of the target earthquake horizon.

Images