Oily sandstone prediction method based on CNN-BILSTM mixed model

Abstract

The invention relates to the technical field of reservoir prediction, and discloses an oil-containing sandstone prediction method based on a CNN-BILSTM hybrid model, and the method comprises the following steps: data preparation: carrying out the standardization processing of seismic data; constructing a three-dimensional sample library, and realizing sample equalization through virtual sample generation and dynamic boundary adjustment; hybrid model construction: constructing a deep learning model comprising a multi-scale convolutional layer, a bidirectional LSTM layer and a cross-modal attention fusion layer; model training and optimization: carrying out model training by adopting a dynamic batch strategy and a mixed loss function; and performing whole-region prediction, inputting seismic data into the trained model, outputting probability distribution of the oil-containing sandstone, and performing post-processing to generate a final prediction result. According to the method, a three-dimensional sample library with strong geological representativeness is constructed through a virtual sample generation and dynamic equilibrium strategy of geological-geophysical cooperative constraint, the problems of actual sample sparsity and category imbalance are solved, and the sensitivity and discrimination capability of the model to hidden oil-bearing sandstone under a complex geological background are remarkably enhanced.

Description

Technical Field

[0001] The present invention relates to the technical field of reservoir prediction, and in particular to an oil-bearing sandstone prediction method based on a CNN-BILSTM hybrid model. Background Art

[0002] Oil-bearing sandstone prediction is a core component of reservoir characterization in oil and gas exploration, directly determining the optimization of drilling targets and the evaluation of development benefits. As exploration targets extend to concealed reservoirs with complex lithologies and strong heterogeneity, traditional seismic interpretation methods, relying on manual experience and single-attribute analysis, struggle to quantitatively characterize the spatial distribution patterns of targets such as thin interbedded sand bodies and fault-blocked reservoirs. The industry urgently needs to develop intelligent prediction technologies to overcome multiple bottlenecks, including low manual interpretation efficiency, poor model generalization, and insufficient geological adaptability, to achieve high-precision identification and three-dimensional visualization of oil-bearing sandstones under complex geological conditions.

[0003] The current technology system mainly includes two types of methods: classification models based on seismic attribute optimization and statistical learning, and end-to-end prediction methods based on deep learning. The former constructs a feature space by manually screening seismic attributes and combines classifiers such as support vector machines and random forests to achieve oil-bearing discrimination, but its prediction accuracy is limited by attribute redundancy and weakened geophysical correlation; the latter uses convolutional neural networks to directly extract features from seismic data. Although it achieves high accuracy in local work areas, it is highly sensitive to the number and distribution balance of samples and lacks explicit modeling of the dynamic evolution of geological sedimentation and structure. It is worth noting that existing methods mostly rely on near-well seismic traces to construct training samples, which makes it difficult to solve the problems of sample sparsity in areas with few wells and generalization across work areas. It is also unable to effectively integrate geological prior knowledge such as faults and sedimentary phases to optimize prediction results.

[0004] The core flaw of existing technologies lies in the disconnect between geological laws and data-driven models, which is specifically manifested in the insufficient embedding of geological prior knowledge in the intelligent prediction process. Since virtual sample generation lacks sedimentary phase control and physical property constraints, model training is easily affected by sample bias, resulting in distortion of sand body morphology during cross-work area predictions. The model architecture design does not take into account both spatial sedimentary feature extraction and geological time series evolution modeling, making it difficult to capture the spatial distribution patterns of reservoirs under tectonic transformation effects such as fault cutting and stratigraphic erosion. In addition, the post-processing link of the prediction results lacks a geological credibility calibration mechanism, making it impossible to suppress non-geological anomalies through operations such as fault shielding and sedimentary phase probability correction, which restricts the geological interpretability and practical application value of the results. This defect makes it difficult for existing technologies to meet the comprehensive requirements of algorithm accuracy, generalization ability, and geological adaptability for complex and concealed reservoir exploration, severely limiting the large-scale application of intelligent oil and gas exploration technologies. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for predicting oil-bearing sandstone based on a CNN-BILSTM hybrid model, which solves the problems of low prediction accuracy, poor cross-area generalization, and insufficient geological interpretability in existing oil-bearing sandstone prediction technologies caused by the separation of geological laws and data-driven models.

[0006] To achieve the above objectives, the present invention is implemented through the following technical solutions: The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model includes the following steps: Data preparation: Load seismic data and well data, establish depth-time mapping through well-seismic time-depth calibration, and standardize the seismic data; 3D sample library construction: Delineate oil-bearing areas based on geological knowledge, label oil-bearing / non-oil-bearing samples based on well interpretation, and achieve sample balance through virtual sample generation and dynamic boundary adjustment; Hybrid model construction: Build a deep learning model that includes multi-scale convolutional layers, bidirectional LSTM layers, and cross-modal attention fusion layers; Model training and optimization: Dynamic batch strategy and hybrid loss function are used for model training; Full-area prediction: seismic data is fed into the trained model, which outputs the probability distribution of oil-bearing sandstones and performs post-processing to generate the final prediction results.

[0007] Preferably, the well seismic time-depth calibration specifically includes: Generate synthetic seismic records of the well bypass channel, use the least squares matching algorithm to perform waveform matching between the synthetic seismic records and the actual seismic trace, and establish a conversion relationship between the depth domain and the time domain; The standardization process uses global mean variance normalization, and the formula is: Where μ and σ are the mean and standard deviation of the seismic data volume, respectively.

[0008] Preferably, the well data includes: Wellhead data: well name, plane coordinates; Well trajectory data: depth, well inclination, azimuth; Time-depth curve: the corresponding relationship between depth and time; Interpretation conclusion: start depth, end depth and label code of the oil-bearing layer.

[0009] Preferably, the virtual sample generation includes: Simulating the physical properties of the virtual well based on the sedimentary facies model, including porosity, permeability, and water saturation; and generating seismic response data of the virtual well through forward modeling of the acoustic wave equation. Automatically label the oil content of virtual samples according to the threshold value of physical property parameters.

[0010] Preferably, the dynamic boundary adjustment is achieved by iteratively shrinking the training area boundary, so that the ratio of oil-containing samples to non-oil-containing samples is optimized from an initial 1:24 to 1:2.5-1:3.

[0011] Preferably, the multi-scale convolutional layer includes: Three asymmetric convolution kernels (3×3, 5×5, and 7×7) are used in parallel to extract local-to-global features. The multi-scale feature maps are spliced ​​along the channel dimension to form comprehensive spatial features.

[0012] Preferably, the cross-modal attention fusion layer is implemented in the following manner: Perform tensor outer product operations on convolution features and temporal features to generate a three-dimensional attention weight matrix; The spatial features and temporal features are dynamically weighted and fused based on the weight matrix.

[0013] Preferably, the hybrid loss function is composed of focal loss and contrast loss weighted by a preset weight coefficient, and its expression is: L=αL focal +βL contrast ; in: α and β are weighted coefficients, satisfying α+β=1, and α∈[0.5,0.9], β∈[0.1,0.5], L focal =-∑y(1-p) γ logp, γ is the adjustment factor, the value range is γ∈[1,3], which is used to control the weight difference of difficult and easy samples; Represents the square loss of feature distance between similar samples.

[0014] Preferably, the dynamic batch strategy adjusts the batch size according to the following rules: The initial batch size is 64, which increases by 8 every 100 iterations, with an upper limit of 128; The learning rate adopts the cosine annealing strategy, with an initial value of 3×10 -4 , the minimum value is 1×10 -5 .

[0015] Preferably, the post-processing includes For the oil-containing probability body P oil ∈[0,1] performs three-dimensional median filtering with a filter window size of 3×3×3; Set the probability threshold P th ∈[0.6,0.8], the effective oil-bearing area is divided into {(x,y,t)|P oil (x,y,t)≥P th}; The confidence level of the prediction result C∈[0,1] is calculated based on information entropy. The formula is: C=1-[-P oil logP oil -(1-P oil )log(1-P oil )]; Among them, the area with C ≥ 0.8 is marked as a high-confidence oil-bearing area.

[0016] In summary, the present invention includes at least one of the following beneficial technical effects: 1. This invention utilizes a virtual sample generation and dynamic balancing strategy based on geological and geophysical collaborative constraints to construct a three-dimensional sample library with strong geological representation, effectively addressing the sparsity and class imbalance issues of actual samples. Combining a hybrid model architecture of multi-scale CNN and bidirectional LSTM, this model achieves collaborative modeling of spatial sedimentary characteristics and temporal evolution patterns, significantly enhancing the model's sensitivity and ability to identify hidden oil-bearing sandstones in complex geological settings.

[0017] 2. This method employs a joint loss function and dynamic gradient optimization strategy to simultaneously optimize classification accuracy, virtual-to-real sample distribution alignment, and model complexity constraints. Through adversarial training and DropPath regularization, the model's robustness to interference factors such as seismic noise and missing data is enhanced, ensuring stable prediction performance across work areas and geological settings.

[0018] 3. The present invention uses a post-processing algorithm based on temperature-scaling confidence calibration and fault-sedimentary phase constraints to seamlessly integrate geological prior knowledge into the prediction process, suppressing anomalous results caused by non-geological factors, ensuring the consistency of the spatial distribution of oil-bearing probability with sedimentary patterns and tectonic evolution, and improving the geological credibility and interpretability of the prediction results. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 This is the flow chart of the CNN-BILSTM hybrid deep learning oil-bearing sandstone prediction process of the present invention; Figure 2 This is a diagram showing changes in the training area range of the present invention; Figure 3 This is a structural diagram of the CNN-BILSTM hybrid deep learning network model of the present invention. DETAILED DESCRIPTION

[0020] The following is combined with Figure 1 -Attached Figure 3 , the present invention is described in further detail.

[0021] The present invention provides a method for predicting oil-bearing sandstone based on a CNN-BILSTM hybrid model, which may include the following steps: Step S1: Data preparation and preprocessing In this embodiment, the data preparation and preprocessing stage is used to achieve high-precision spatiotemporal alignment and normalization of well data and seismic data. The technical details are fully disclosed below: 1. Seismic wavelet extraction and reflection coefficient calculation Seismic wavelet extraction: The time-varying seismic wavelet is extracted from the seismic trace near the well, specifically using the statistical wavelet estimation method: The seismic data within a 10×10 plane around the well seismic trace were selected as the training set; The eigenvectors of the covariance matrix are calculated by principal component analysis (PCA), and the waveforms corresponding to the principal components are extracted as wavelet basis functions. The wavelet basis functions are subjected to time-varying correction to generate a time-varying wavelet W(T) that matches the absorption and attenuation characteristics of the formation. The expression is: W(T) = W0·e -αT ·cos(2πf c T+φ); Where W0 is the initial amplitude, α is the attenuation coefficient, and f c is the main frequency, and φ is the phase angle.

[0022] Reflection coefficient calculation: Based on the acoustic transit time (AC) and density (DEN) logging curves, the reflection coefficient sequence is calculated layer by layer: The well logging curves were resampled in the depth domain with a sampling interval of 0.125 m; The exact reflection coefficient is calculated using the Zoeppritz equation: Where ρ is the density, v is the longitudinal wave velocity, μ is the shear modulus, and Δμ is the difference in shear modulus between layers; The reflection coefficient sequence is band-pass filtered (5Hz-80Hz) to remove high-frequency noise and DC offset.

[0023] 2. Synthetic seismic record generation and waveform matching Synthetic record generation: Perform one-dimensional time-domain convolution of the reflection coefficient sequence with the seismic wavelet to generate a synthetic seismic record: Among them, N w is the wavelet length, and Δt is the time sampling interval (2ms).

[0024] Waveform matching optimization: Levenberg-Marquardt algorithm is used to minimize the residual energy between synthetic records and actual seismic traces: The objective function is defined as mean square error (MSE): Where θ={T align ,α,f c ,φ} is the parameter to be optimized; Calculate the partial derivatives of the objective function with respect to the parameters and construct the Jacobian matrix J; Iteratively update parameters: θ (n+1) =θ (n) -(J T J+λI) -1 J T ∈; Among them, λ is the damping factor, ∈ is the residual vector; Termination condition: the parameter change is less than 10-6 or the number of iterations exceeds 100.

[0025] Data standardization: 1. Global mean and standard deviation calculation Three-dimensional data volume traversal calculation: seismic data volume Perform point-by-point traversal: Mean calculation: Standard deviation calculation: 2. Channel-by-channel normalization and truncation processing: For each seismic trace D raw (x,y,:) is normalized: Truncate the normalized data: Well data loading and structured processing: 1. Wellhead data analysis and coordinate conversion: Parse the well header file (ASCII format) and extract the fields: wellname(string, unique identifier); Plane coordinates (x, y) (floating point, unit: meter); The coordinate system is converted to UTM projection (such as EPSG:32648). The conversion formula is: x UTM =x WGS84 +k0·cosφ·(λ-λ0); Where k0 is the scale factor and λ0 is the central meridian.

[0026] 2. Spatial interpolation of well trajectory data: Perform cubic spline interpolation on depth (MD), well inclination (θ), and azimuth (φ): Input the original measurement point data Sampling interval ≤ 0.5 m; Construct spline basis functions: S(MD)=a i +b i (MD-MD i )+c i (MD-MD i ) 2 +d i (MD-MD i ) 3 ; Solve the coefficient a through boundary conditions (natural spline) i ,b i ,c i ,d i .

[0027] Non-uniform sampling processing of time-depth curve: Time-depth curve Perform resampling: Linear interpolation was used to generate a time-depth relationship table with uniform depth intervals (0.125 m); Akima interpolation was performed on missing segments to ensure continuity.

[0028] Seismic data quality control: 1. Weighted interpolation algorithm for missing road repair For the missing trace D(x,y,:), calculate the weighted average of the eight adjacent valid traces: Among them, the weight (Avoid division by zero).

[0029] 2. Adaptive median filtering for high-frequency noise suppression Use a 3×3×3 three-dimensional median filter: For each sampling point (x, y, t), extract 27 data points within the neighborhood window; After sorting, take the median to replace the current point: D filtered (x, y, t) = median({D(x+i, y+j, t+k) | i, j, k∈{-1, 0, 1}}); 3. Data consistency check Verify the integrity of the seismic data volume: Verify the number of tracks N in the SGY file header traces Is it consistent with the actual data dimension? Check whether the time sampling interval Δt is consistent with the header file declaration (error < 0.01ms); Detect amplitude overflow (|amplitude|>32767) and record the abnormal location.

[0030] Step S2: Construction of 3D sample library In this embodiment, the three-dimensional sample library is constructed through sample annotation based on geological and geophysical collaborative constraints, sedimentary facies-controlled virtual sample generation, and dynamic boundary iterative optimization to achieve high-precision modeling and balanced processing of oil-bearing sandstone samples. The specific technical content is fully disclosed as follows: Multi-dimensional constraints on sample labeling rules: Mathematical modeling of planar position constraints: Reservoir boundary Ω oil It is defined by the spatial intersection of the structural interpretation polygon and the sedimentary facies boundary. The structural interpretation polygon passes through the fault trace line Γ fault and the formation closing line Γ closure The topological stitching generation of , the specific steps include: 1. Perform B-spline curve fitting on the fault trajectory to generate a continuous and smooth fault surface S fault (x,y); 2. Using the stratum thickness contour Γ isopach and structural contour Γ contour 3. Use the Delaunay triangulation algorithm to construct a planar mesh of the closed area, and finally extract the boundary from the outer contour of the mesh: Dynamic calculation of longitudinal time windows: The time window length Δt is determined based on the acoustic time difference characteristics of the target sand body and the difference in surrounding rock velocity. Specifically, the time window is calculated based on the sand body thickness hsand and average porosity φ interpreted by well logging, combined with the surrounding rock velocity vshale: Among them, v shale is the sandstone velocity, according to the porosity-velocity relationship Calculate, v matrix is the skeleton velocity, v fluid is the fluid velocity.

[0031] Fault-tolerant mechanism for label assignment: For seismic traces near the well, if there is a local contradiction between the interpretation conclusion and the seismic response (for example, weak amplitude corresponding to the oil-bearing layer), the Bayesian probability model is used to correct the label: Among them, A seis is the amplitude characteristic of the seismic trace, the prior probability P(y=1) is determined by the well interpretation statistics, and the likelihood function P(Aseis |y) is modeled by kernel density estimation.

[0032] Physical property-seismic response coupling modeling of virtual sample generation: In this embodiment, the virtual sample is generated by modeling the physical parameters controlled by sedimentary phases and forward modeling of high-precision wave equations to achieve the physical mechanism coupling between the physical parameters and the seismic response. The specific technical content is fully disclosed as follows: Mathematical Implementation of Joint Distribution Modeling of Physical Property Parameters Sedimentary facies-controlled mixed distribution model of porosity: The probability density function of porosity φ is defined by the sedimentary facies-dependent Gaussian mixture model (GMM): in: θ depo are sedimentary facies parameters (e.g., sand-to-ground ratio, sedimentary energy index); π k is the mixing coefficient, which is calculated by the sand-to-ground ratio SGR through the Softmax function: μ k is the mean of the kth Gaussian component, and is related to the deposition energy index E depo Linear correlation: μ k =α k ·E depo +β k ; is the variance, which is determined by the sediment heterogeneity index H heter adjust: Physical property coupling equation of permeability: The permeability k is modeled by the nonlinear relationship between the porosity φ and the particle size parameter: lnk=a0+a1φ+a2ln(d 50 )+a3φln(d 50 )+∈ noise ; in: d 50 is the median particle size, which obeys the log-normal distribution: is the core experiment error term; The coefficients a0, a1, a2, and a3 are calibrated by least squares regression.

[0033] Capillary pressure constraint model for water saturation: Water saturation S w By capillary pressure P cCombined with the fluid potential energy balance equation to determine: in: P entry is the mercury injection pressure, and the pore throat radius r t calculate: λ is the pore structure index, which is related to porosity and sorting coefficient: λ=c1φ+c2·Sorting+c3; 2. Details of numerical solution of wave equation forward modeling Discretization of three-dimensional elastic wave equation: The elastic wave equation is discretized using the Virieux finite difference scheme with velocity-stress staggered grid: in: v i is the particle velocity component; σ ij is the stress tensor component; λ and μ are Lame constants, which are given by the velocity v p 、v s And density ρ calculation: Discrete format for space and time: Spatial discretization: Use the 4th order central difference format, for example Approximation: Time discretization: Use the second-order leapfrog scheme to update the velocity and stress fields: Among them, D j Represents the spatial difference operator. The mathematical expression of the absorbing boundary condition (PML) is: A PML layer is set at the boundary of the calculation area, and its complex coordinate stretching function is: in: d0 is the maximum attenuation coefficient; L PML is the thickness of the PML layer; ω is the angular frequency.

[0034] The modified form of the wave equation in the PML region is: Noise injection model of virtual seismic traces Frequency domain modeling of Gaussian white noise: Generate zero-mean Gaussian white noise in the frequency domain, and its power spectral density (PSD) is constant: Among them, A signal is the mean amplitude of the seismic signal, and SNR is the signal-to-noise ratio.

[0035] Synthesis method of coherent noise: simulate the interference of multiple waves and surface waves and generate them through time-varying filtering: Among them, the Morlet wavelet basis function is: Parameter τ k is a random delay time, which obeys the uniform distribution τ k ~U(0,T max ).

[0036] Time-varying gain model of amplitude distortion: random time-varying gain perturbations are applied to the seismic trace to simulate surface consistency anomalies: G(t) = 1 + γ·sin(2πf g t+φ g )·e -βt ; Among them, f g is the disturbance frequency, φ g is the random phase, and β is the attenuation coefficient.

[0037] Gradient-driven iterative algorithm with dynamic boundary adjustment: Quantitative evaluation of sample imbalance ratio Definition of sample imbalance ratio R imb =N neg / N pos , where N pos With N neg are the number of positive and negative samples respectively.

[0038] The gradient analytical calculation of the boundary parameters will be the training area boundary B train Parameterized as an implicit level set function φ(x,y), the derivative of the sample size with respect to the boundary is calculated using shape derivative theory: Where n is the boundary normal vector and δ(φ) is the Dirac function. The imbalance ratio gradient is: The adaptive step size control of boundary update uses Armijo line search to determine the optimal shrinkage step size η, ensuring that the objective function R imb Monotonically decreasing: 1. Initialize step size η0 = 0.1, decay factor τ = 0.5; 2. Iterative search satisfies The minimum k such that η=η0τ k ; 3. Update boundaries:

[0039] 4. Gradual Modeling of Transition Zone Labels At boundary B train Set the transition zone Ω on the outside buffer , whose label value decays nonlinearly according to the geodesic distance: Among them, d geod (x,y) is the shortest path distance from the sample point to the boundary (considering the fault blocking effect), σ b is the attenuation bandwidth parameter.

[0040] Step S3: CNN-BiLSTM Hybrid Model Construction In this embodiment, the CNN-BiLSTM hybrid model achieves collaborative characterization of seismic data spatial and temporal features and high-precision prediction of oil-bearing sandstones through multi-scale spatial feature extraction, bidirectional temporal dependency modeling, and cross-modal dynamic feature fusion. The technical implementation details of the model architecture are fully disclosed below: The spatial feature extraction mechanism of the multi-scale convolution layer is composed of three convolution kernels of different sizes deployed in parallel, aiming to capture multi-level geological features in seismic data through differentiated receptive fields. The specific implementation includes: Local texture feature extraction (3×3 convolution kernel): A small-scale convolution kernel focuses on local details of the seismic waveform (such as peak-trough combinations and amplitude jumps). The output feature map reflects subtle differences in lithologic interfaces. Such features are critical for identifying thin interbeds between sandstone and mudstone.

[0041] Medium-scale structural feature extraction (5×5 convolution kernel): A medium receptive field covers geological structures such as sandstone boundaries and fault occlusions. This increases sensitivity to discontinuous geological interfaces by expanding the spatial context perception range. This type of feature is effective in characterizing sedimentary phenomena such as river channel incision and sandstone pinchout.

[0042] Macro-sedimentary pattern feature extraction (7×7 convolution kernel): The large-scale convolution kernel captures macro-sedimentary patterns such as sheet sand distribution and river channel migration. Its output feature map is strongly correlated with the spatial distribution of sedimentary facies, providing background constraints for regional oil content prediction.

[0043] The output feature maps of each convolution kernel are spliced ​​together in the channel dimension to achieve multi-scale fusion and form a comprehensive spatial feature tensor: F cnn =Concat(Conv 3×3 (X; W1), Conv 5×5 (X; W2), Conv 7×7 (X; W3))∈R H×W×3C Where W1, W2, and W3 are the weight matrices of each convolution kernel, H and W are the spatial dimensions of the feature map, and C is the number of output channels of a single convolution kernel. Preferably, each convolution layer is followed by batch normalization (BatchNorm) and a ReLU activation function to accelerate training convergence and enhance nonlinear expression capabilities.

[0044] The bidirectional LSTM layer simulates the forward evolution of the stratigraphic sedimentary sequence and the reverse impact of geological transformation (such as erosion and faulting) through forward and reverse time propagation mechanisms, respectively, to fully depict the dynamics of geological time series: - Forward LSTM unit: Propagates forward along the seismic time axis, simulating the vertical inheritance of stratigraphic thickness and lithologic combinations in the sedimentary sequence. Its hidden state update formula is: f t =σ(W f ·[H t-1 ,X t ]+b f ); i t =σ(W i ·[H t-1 ,X t ]+b i ); o t =σ(W o ·[H t-1 ,X t ]+b o ); Among them, f t ,i t ,o t are the activation values ​​of the forget gate, input gate and output gate respectively, C t is the cell state, and ⊙ represents element-by-element multiplication.

[0045] Reverse LSTM unit: propagates in the reverse direction along the earthquake time axis, simulating the transformation effect of tectonic movement (such as uplift, erosion, and fault activity) on the original sedimentary sequence. Its hidden state update formula is consistent with the forward unit structure, but the input sequence order is reversed: The forward and reverse outputs are superimposed through residual connections to generate comprehensive time series features: in, It represents element-by-element addition, preserving the complementarity and consistency of bidirectional temporal information, T is the number of time steps, and D is the hidden layer dimension.

[0046] Dynamic feature interaction mechanism of the cross-modal attention fusion layer The cross-modal attention fusion layer achieves collaborative optimization of CNN and BiLSTM output features through interactive modeling of spatial-temporal features and dynamic weight allocation: 1. Feature interaction modeling: transform the spatial feature F cnn ∈R H×W×3C Expanded into a vector sequence along the spatial dimension Timing characteristics H stack ∈R T×D Expand along the time dimension to Construct the space-time correlation matrix through tensor outer product operation: 2. Attention weight generation: Perform three-dimensional convolution and nonlinear transformation on the correlation matrix to generate normalized attention weight: W att =σ(Conv3D(M;K))∈[0,1] H×W×T ; Among them, Conv3D uses a 3×3×3 convolution kernel K to capture local spatial-temporal correlation; σ is the Sigmoid function, which maps the weights to the [0,1] interval.

[0047] 3. Dynamic feature fusion: Spatial and temporal features are weighted and aggregated according to attention weights. The fusion formula is: in, The element-by-element multiplication of the representation vectors enhances the fine-grained correlation of feature interactions. Preferably, the fusion result is mapped to the target dimension through a fully connected layer, and the oil content probability prediction is output.

[0048] Step S4: Model training and optimization In this embodiment, the model training and optimization phase achieves efficient training and improved generalization capabilities of the CNN-BiLSTM hybrid model by combining loss function design, dynamic weight allocation strategy, and regularization constraints. The specific technical content is fully disclosed as follows: Joint loss function: The joint loss function consists of cross entropy loss, KL divergence constraint and L2 regularization term, aiming to simultaneously optimize classification accuracy and model robustness: Weighted cross entropy loss (main loss term): To address the problem of sample category imbalance, the category weight coefficient α is introduced pos With α neg , dynamically adjust the loss contribution of positive and negative samples: in, To balance the impact of sample size differences on gradient updates.

[0049] KL divergence constraint (distribution alignment term): By minimizing the KL divergence between the predicted distribution of virtual samples and actual samples, the model's generalization ability for virtual data is enhanced: Among them, p vir With p real are the predicted probability distributions of virtual samples and actual samples respectively, and the KL divergence is calculated as: Preferably, gradient locality sensitive hashing (LSH) is used to approximate the KL divergence to reduce the computational complexity.

[0050] L2 regularization term (model complexity constraint): suppresses model overfitting and imposes L2 norm penalty on all weight parameters W: Where λ is the regularization strength coefficient, which is preferably dynamically adjusted through Bayesian optimization. The joint loss function is a linear combination of the above three terms: L total =L CE +βL KL +L reg Among them, β is the balance coefficient of the KL term, the initial value is set to 0.1, and it decays exponentially with the number of training rounds.

[0051] Gradient optimization mechanism of dynamic weight allocation strategy: To solve the gradient conflict problem between different modules of the model (CNN and BiLSTM), a dynamic weight allocation strategy is adopted to coordinate the optimization direction of feature extraction and time series modeling: Gradient similarity metric: Calculate the gradient cosine similarity between CNN and BiLSTM to evaluate the optimization consistency between modules: Dynamic weight adjustment: Adjust the learning rate ratio of CNN and BiLSTM according to similarity: Wherein, γ is the attenuation coefficient, and preferably, the adjustment process is smoothed by the sliding window mean.

[0052] In addition to L2 regularization, the collaborative optimization method of the regularization strategy uses the following regularization techniques to suppress overfitting: DropPath regularization: randomly discard part of the gradient flow on the residual connection path to enhance the model's robustness to feature redundancy.

[0053] Specifically, the residual branch of the lth layer is set to zero with probability p: H (l+1) =H (l) +Mask(p)·F(H (l) ); Among them, Mask(p) is the mask matrix that obeys the Bernoulli distribution, and the drop probability p increases linearly with the number of training rounds.

[0054] Adversarial Training: Apply a small perturbation δ to the input seismic data to generate adversarial samples and participate in training to improve the model's robustness to noise: Wherein, ∈ is the perturbation intensity coefficient. Preferably, the projected gradient descent (PGD) algorithm is used to iteratively generate multi-step adversarial samples.

[0055] Adaptive learning rate scheduling and early stopping mechanism: Cosine annealing learning rate scheduling: The learning rate is adjusted according to the cosine function period to balance the convergence speed and stability: Among them, T cycle is the cycle length, preferably, set to 1 / 4 of the total training rounds.

[0056] Early Stopping: Monitor the validation set loss and terminate training if it does not decrease after K consecutive rounds to avoid overfitting. Preferably, the model parameters with the lowest validation loss are retained as the final model.

[0057] Step S5: Oil-bearing sandstone prediction and 3D visualization In this embodiment, the oil-bearing sandstone prediction and 3D visualization stage achieves the generation of a highly reliable oil-bearing probability volume and the visualization of its 3D spatial distribution through model reasoning, confidence calibration, and geologically constrained post-processing algorithms. The specific technical content is fully disclosed as follows: Sliding window prediction method for three-dimensional oil-bearing probability volume: Based on the trained CNN-BiLSTM hybrid model, a sliding window inference with full spatial coverage is performed on the seismic data volume of the target work area. The sliding window size is consistent with the training sample size, and the window moves along the inline, crossline, and time axis with a preset step size to ensure spatial continuity. For overlapping window areas, a weighted average fusion strategy is used to eliminate edge effects: in, is the local prediction probability of the kth overlapping window, weight w kInversely proportional to the Euclidean distance from the window center to the current voxel: Preferably, ∈ is a very small value (such as 1e-6) to avoid division by zero errors.

[0058] Bayesian probability correction algorithm for confidence calibration In order to improve the geological interpretability of the predicted probability, the temperature scaling method based on confidence calibration is used to correct the original output probability: Where z is the Logit value of the model output, T is the calibration temperature coefficient, and it is optimized by maximizing the negative log-likelihood function of the validation set probability: Preferably, the L-BFGS algorithm is used to solve the optimal T to ensure the statistical consistency of the calibrated probability and the true oiliness label.

[0059] Constraints on geological interpretation results: 1. Fault shielding: Based on the fault polygon Ωfault from structural interpretation, an attenuation coefficient is applied to the predicted probability within the fault zone: P masked (x,y,t)=P calibrated (x,y,t)·exp(-λ·d fault (x,y,t)); Among them, d fault is the spatial distance from the voxel to the nearest fault, and λ is the attenuation intensity coefficient, which is preferably dynamically adjusted with the fault activity level.

[0060] 2. Sedimentary facies-controlled probability enhancement: The sedimentary facies model S(x, y, t) is used as prior knowledge and the oil-bearing probability is modified using the Bayesian theorem: Where P(S|Oil) is the conditional probability of the sedimentary facies under oil-bearing conditions, which is obtained from well data statistics; P(S) is the prior probability of the sedimentary facies.

[0061] 3. Morphological post-processing: Use 3D morphological opening (erosion followed by dilation) to remove isolated noise points. The structural element is a 3×3×3 sphere: in, represents the corrosion operation, represents the expansion operation, and B is the structural element.

[0062] 3D visualization and interactive interpretation tools: 1. Volume Rendering and Transparency Transfer Function In this embodiment, the ray casting algorithm maps the oil probability value to the visual attribute through the color and transparency transfer function. The specific mathematical formula is defined as follows: The mapping of the oil probability value P∈[0,1] of the color component mapping function to the RGB color space adopts piecewise linear interpolation: Low probability area (P < 0.3): the green component increases linearly with the probability, and the red component is zero, simulating an oil-free background; Medium probability region (0.3≤P<0.7): The red component increases linearly, and the green component decreases, achieving a gradient from yellow (R=0.5, G=0.5) to orange-red (R=1, G=0); High probability area (P ≥ 0.7): The red component is the maximum value and the green component is zero, highlighting the high-confidence oil-bearing area.

[0063] The transparency mapping function transparency α(P)∈[0,1] is segmented controlled by probability value: Low probability area: completely transparent to avoid occlusion of irrelevant areas; Medium and high probability areas: Transparency increases linearly with probability, ensuring a hierarchical expression of the spatial distribution of oil-bearing sandstone bodies.

[0064] Light accumulation formula: Accumulate color and transparency along the ray path r(s) to calculate the final pixel color C final : in: C(P i )=(R(P i ),G(P i ),B(P i )) is the color of the i-th sampling point; Indicates the cumulative forward transmittance; The sampling step size Δs is set according to the data volume resolution. Preferably, an adaptive step size optimization algorithm is used to balance rendering quality and computational efficiency.

[0065] 2. Isosurface extraction and geological attribute overlay: Extract the isosurface of the specified probability threshold (such as 0.5) and fuse it with the seismic amplitude attribute (such as RMS amplitude) to display: A fused (x,y,t)=A RMS (x,y,t)·exp(-α·(1-P clean (x,y,t))); Among them, α is the fusion coefficient, which controls the superposition strength of seismic attributes and oil-bearing probability.

[0066] 3. Interactive Profile Analysis Tool: It supports the generation of well-connected profiles along any azimuth angle, and synchronously displays seismic traces, oil probability curves, and well interpretation conclusions. The coordinate transformation formula is: x′=xcosθ-ysinθ, y′=xsinθ+ycosθ; Where θ is the profile azimuth, which is adjusted in real time using the interactive rotation control.

[0067] In summary, a three-dimensional sample library was constructed through geological-geophysical collaborative virtual sample generation and dynamic balancing strategy. A multi-scale CNN was used to capture spatial sedimentary features, a bidirectional LSTM was used to model geological temporal evolution, and a cross-modal attention mechanism was combined to achieve feature fusion. Model training was completed based on a joint loss function and dynamic gradient optimization. High-confidence probability volumes were generated through post-processing with sliding window inference, temperature scaling calibration, and fault-sedimentary phase constraints. Finally, a ray casting algorithm and interactive tools were used to visualize the three-dimensional spatial distribution of oil-bearing sandstones. This systematically solves the problems of sample sparsity, model generalization, and geological interpretability faced by the identification of hidden reservoirs in seismic data.

[0068] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model is characterized by: The following steps are involved: Data preparation: Load seismic data and well data, establish depth-time mapping through well-seismic time-depth calibration, and standardize the seismic data; 3D sample library construction: Delineate oil-bearing areas based on geological knowledge, label oil-bearing / non-oil-bearing samples based on well interpretation, and achieve sample balance through virtual sample generation and dynamic boundary adjustment; Hybrid model construction: Build a deep learning model that includes multi-scale convolutional layers, bidirectional LSTM layers, and cross-modal attention fusion layers; Model training and optimization: Dynamic batch strategy and hybrid loss function are used for model training; Full-area prediction: seismic data is fed into the trained model, which outputs the probability distribution of oil-bearing sandstones and performs post-processing to generate the final prediction results.

2. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The well seismic time-depth calibration specifically includes: Generate synthetic seismic records of the well bypass channel, use the least squares matching algorithm to perform waveform matching between the synthetic seismic records and the actual seismic trace, and establish a conversion relationship between the depth domain and the time domain; The standardization process uses global mean variance normalization, and the formula is: Where μ and σ are the mean and standard deviation of the seismic data volume, respectively.

3. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The well data includes: Wellhead data: well name, plane coordinates; Well trajectory data: depth, well inclination, azimuth; Time-depth curve: the corresponding relationship between depth and time; Interpretation conclusion: start depth, end depth and label code of the oil-bearing layer.

4. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The virtual sample generation includes: Simulating the physical properties of the virtual well based on the sedimentary facies model, including porosity, permeability, and water saturation; and generating seismic response data of the virtual well through forward modeling of the acoustic wave equation. Automatically label the oil content of virtual samples according to the threshold value of physical property parameters.

5. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The dynamic boundary adjustment is achieved by iteratively shrinking the training area boundary, so that the ratio of oil-containing samples to non-oil-containing samples is optimized from the initial 1:24 to 1:2.5–1:

3.

6. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The multi-scale convolutional layer includes: Three asymmetric convolution kernels (3×3, 5×5, and 7×7) are used in parallel to extract local-to-global features. The multi-scale feature maps are spliced ​​along the channel dimension to form comprehensive spatial features.

7. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The cross-modal attention fusion layer is implemented in the following way: Perform tensor outer product operations on convolution features and temporal features to generate a three-dimensional attention weight matrix; The spatial features and temporal features are dynamically weighted and fused based on the weight matrix.

8. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The hybrid loss function is composed of focal loss and contrast loss weighted by a preset weight coefficient, and its expression is: L=αL focal +βL contrast ; in: α and β are weighted coefficients, satisfying α+β=1, and α∈[0.5,0.9], β∈[0.1,0.5], L focal =-∑y(1-p) γ logp, γ is the adjustment factor, the value range is γ∈[1,3], which is used to control the weight difference of difficult and easy samples; Represents the square loss of feature distance between similar samples.

9. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The dynamic batch strategy adjusts the batch size according to the following rules: The initial batch size is 64, which increases by 8 every 100 iterations, with an upper limit of 128; The learning rate adopts the cosine annealing strategy, with an initial value of 3×10 -4 , the minimum value is 1×10 -5 .

10. The oil-bearing sandstone prediction method based on the CNN-BILSTM hybrid model according to claim 1 is characterized in that: The post-processing includes For the oil-containing probability body P oil ∈[0,1] performs three-dimensional median filtering with a filter window size of 3×3×3; Set the probability threshold P th ∈[0.6,0.8], the effective oil-bearing area is divided into {(x,y,t)|P oil (x,y,t)≥P th }; The confidence level of the prediction result C∈[0,1] is calculated based on information entropy. The formula is: C=1-[-P oil logP oil -(1-P oil )log(1-P oil )]; Among them, the area with C ≥ 0.8 is marked as a high-confidence oil-bearing area.

Images