A method for identifying sweet spots in shale oil and gas reservoirs

Abstract

This application relates to the field of shale oil and gas, and discloses a method for identifying sweet spots in shale oil and gas reservoirs, including the following steps: acquiring core CT images, 3D seismic data, and production dynamic data and performing cross-scale preprocessing; constructing a fractional-order nonlocal seepage field model to characterize flow characteristics from nanometer to kilometer scales; predicting seepage field parameters through a Lie group symmetry-constrained neural network; establishing a cross-scale coupling model of quantum adsorption effect and macroscopic seepage law; dynamically updating model parameters based on real-time monitoring data; and performing multi-objective collaborative optimization to generate a 3D distribution of sweet spots and a development plan. This invention employs cross-scale data fusion technology with fractal dimension dynamic constraints, achieving accurate mapping between nanopore and macroscopic fracture network parameters. Through pore connectivity fractal analysis, it solves the problem of spatial registration misalignment between CT scan and seismic inversion data.

Description

Technical Field

[0001] This invention relates to the field of shale oil and gas technology, specifically to a method for identifying sweet spots in shale oil and gas reservoirs. Background Technology

[0002] Shale oil and gas reservoir sweet spot identification technology mainly involves geological exploration, seepage mechanics modeling, and development optimization. Traditional methods include statistical geological modeling, finite element numerical simulation, and data-driven machine learning prediction. Existing technologies encompass CT image pore structure analysis, seismic attribute inversion, Darcy's equation solving, and PID control of injection and production parameters, and are applied to unconventional reservoir evaluation, production prediction, and development plan formulation.

[0003] Existing technologies, when fusing multi-source data, employ static interpolation algorithms, leading to accumulated deviations in the mapping of parameters between nanopores and macroscopic fracture networks, failing to characterize the long-range correlation effects of pore topology. Integer-order seepage models, by simplifying the flow patterns of adsorbed fluids, introduce systematic errors in permeability prediction, inducing overestimation of production capacity. Traditional neural networks neglect physical conservation laws, causing training error propagation and extrapolation failure. Decoupled optimization strategies sever the dynamic correlation between microscopic adsorption and macroscopic seepage, making it difficult to control delays to meet the engineering requirements for rapid response in shale reservoirs. These shortcomings limit the accuracy of sweet spot identification in complex, heterogeneous reservoirs, hindering the economical and efficient development of shale oil and gas. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method for identifying sweet spots in shale oil and gas reservoirs. This method solves the problems in existing technologies, such as cross-scale parameter mapping deviations caused by static interpolation of multi-source data, distortion of the characterization of adsorbed fluids by integer-order seepage models, instability of physical laws in data-driven models, and sluggish response of decoupled optimization strategies.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for identifying sweet spots in shale oil and gas reservoirs, comprising the following steps:

[0006] Collect core CT images, 3D seismic data, and production dynamic data, and perform multi-scale preprocessing;

[0007] By integrating core microscopic imaging, seismic wavefield characteristics, and real-time downhole monitoring data, a multi-scale data channel covering nanopores to stratigraphic structures is established. Mechanistically, the spatial topological characteristics of the micropore network exhibit a fractal correlation with the stress distribution of macroscopic structures. The elastic parameters derived from seismic wave velocity inversion are corrected using the adsorption barrier distribution on the pore surface, enabling the coupling of data from different physical dimensions within a unified seepage dynamics framework. The data preprocessing module extracts the pore-matrix interface through grayscale thresholding; its morphological characteristics directly constrain the boundary condition definition of subsequent seepage equations. Anisotropic filtering of seismic attributes, by suppressing random noise, preserves the flow channel information of natural fractures, providing initial field input for multi-scale modeling.

[0008] Construct a fractional-order nonlocal seepage field model to characterize flow characteristics from nanometer to kilometer scale;

[0009] A multiphysics governing equation, incorporating temporal memory effects and long-range spatial interactions, is constructed to describe the non-Darcy flow behavior of shale oil and gas in nanoscale confined spaces and macroscopic fracture networks. At the mechanistic level, the adsorption-desorption process of fluids within nanopores leads to a nonlinear relationship between seepage velocity and pressure gradient. Fractional-order differential operators are introduced to characterize the anomalous diffusion effect caused by pore wall roughness. The temporal fractional-order reflects the relaxation characteristics of molecular adsorption energy on the organic matter surface, while the spatial fractional-order is dynamically calibrated by the fractal dimension of the pore network, enabling the model to simultaneously capture microscopic adsorption retention and crossflow phenomena within the macroscopic fracture network.

[0010] Predicting seepage field parameters using a neural network constrained by Lie group symmetry;

[0011] A symmetry-constrained neural network architecture is designed to encode the physical laws of mass conservation and momentum transfer in the seepage equation into the network topology. Mechanistically, the network's isovariant convolutional layers ensure the prediction results satisfy the covariance requirements of the seepage dynamics equation across different coordinate systems by maintaining the transformation invariance of a special linear group. The Casimir conservation term introduced in the loss function enforces the integral invariance of fluid mass in the time domain, effectively suppressing the cumulative error caused by neglecting constitutive relations in purely data-driven models. The network's input layer directly receives the preprocessed multi-scale feature field, and the velocity distribution predicted by the output layer is fed back to the source term of the seepage equation to achieve closed-loop correction.

[0012] Establish a cross-scale coupling model between quantum adsorption effect and macroscopic seepage law;

[0013] A joint governing equation for quantum chemical adsorption and macroscopic seepage motion was established to reveal the interaction mechanism between methane molecule fluctuations within nanopores and flow conduction through a kilometer-scale crack network. At the mechanistic level, quantum density fluctuations correct the effective permeability in the macroscopic seepage equation through coupling terms, reflecting the hindering effect of adsorption sites on fluid transport on the organic matter surface. Conversely, the pressure gradient of the macroscopic velocity field influences the evolution path of the quantum wave function, forming a cross-scale energy exchange channel. In the coupled model, the electron density distribution output by the quantum computing module is mapped to the continuous medium parameter space through a projection operator, while the macroscopic seepage field constrains the potential energy distribution of the quantum region through boundary conditions, achieving bidirectional dynamic interaction.

[0014] The model parameters are dynamically updated based on real-time monitoring data.

[0015] Based on downhole dynamic monitoring data and model prediction residuals, an adaptive parameter update system is constructed. Mechanistically, an ensemble Kalman filter algorithm is used to correct key parameters such as permeability tensor and adsorption coupling coefficient in real time by fusing multi-source observation data. Its process noise covariance matrix is ​​dynamically adjusted by the quantum density gradient magnitude, enhancing the model's response to sudden geological events. The Bayesian inversion module establishes a Gaussian mixture prior based on historical parameter distributions and approximates the posterior probability density through Markov chain Monte Carlo sampling, ensuring parameter updates converge within physical constraints. The optimized parameter field is synchronously fed back to the seepage model and neural network weights, forming a closed-loop learning link between parameters, model, and data.

[0016] Perform multi-objective collaborative optimization to generate a 3D distribution and development plan for desserts;

[0017] Taking into account the requirements of improved recovery rate, reduced energy consumption, and flow stability, a development scheme that balances economic benefits and engineering safety is generated. At the mechanistic level, the control algorithm uses dynamic programming to balance water injection displacement efficiency with formation fracturing risk. The upper limit of water injection volume is dynamically adjusted by real-time adsorption density to avoid oversaturation and blockage of organic matter pores. The velocity variance term introduced in the optimization objective function reduces energy loss by suppressing eddy current formation in the fracture network, while the discount factor design ensures coordination between long-term development benefits and short-term production indicators. Control commands adjust wellhead pressure and injection-production cycle in real time through distributed actuators. Execution deviations trigger model fine-tuning and emergency parameter inversion, forming an adaptive control system throughout the entire lifecycle.

[0018] Preferably, the preprocessing includes:

[0019] Calculate the fractal dimension of the pore network from core CT images This is obtained by fitting the following linear relationship:

[0020]

[0021] in For scanning scale, This is the grayscale threshold. For fitting constant terms;

[0022] Anisotropic diffusion filtering is performed on the seismic data, and the filtering kernel function is:

[0023]

[0024] in Indicates diffusion time. This is the noise sensitivity coefficient.

[0025] Gray-scale thresholding distinguishes organic pores from inorganic matrix by setting a CT value range (200-400 HU), and the choice of threshold directly affects the morphological characteristics of the pore connectivity region. Fractal dimension calculation quantitatively characterizes the topological heterogeneity of nanopores by statistically analyzing the geometric complexity of the pore connectivity domain at different scanning scales. This parameter serves as an order adjustment factor for the spatial fractional derivative, controlling the intensity of nonlocal effects in the seepage equation. Anisotropic diffusion filtering dynamically adjusts the smoothing direction based on the gradient field of seismic data, enhancing signal continuity along the fracture extension direction while suppressing high-frequency noise perpendicular to the formation interface, preserving the flow path information of the fracture network, and providing a high signal-to-noise ratio input for the initial permeability field construction of the subsequent seepage model.

[0026] A higher fractal dimension of the pore network indicates stronger geometric complexity and spatial correlation of the pore-throat system, leading to a significant long-range memory effect in fluid flow. By linearly mapping the fractal dimension to the order of the spatial fractional derivative, the seepage equation can adaptively reflect the anomalous diffusion behavior of different pore structures: as the fractal dimension increases, the order of the spatial fractional derivative increases accordingly, and the weights of the nonlocal integral kernel function of the equation extend to the far field, more accurately characterizing the tortuous flow characteristics of fluid in tortuous pores. This dynamic mapping mechanism allows the model to automatically adapt to the seepage characteristics of various types of shale reservoirs without relying on empirical assumptions.

[0027] The fractal dimension field directly participates in the spatial fractional-order operator parameterization of the seepage equation, and its spatial distribution characteristics drive the model to adopt differentiated nonlocal integration ranges in different regions. After filtering, the seismic attributes are used to invert the initial permeability field through the amplitude-velocity relationship. Its anisotropic characteristics and the pore fractal dimension jointly constrain the principal axis direction of the permeability tensor. The pore spatial distribution weights mark high-porosity regions and assign them higher loss weights during neural network training, ensuring priority for model prediction accuracy in key reservoir sections. A cascaded data transmission link is formed between the preprocessing module, the seepage model, and the neural network, achieving seamless integration from static structural analysis to dynamic parameter prediction.

[0028] The fractal dimension field output by the CT image analysis unit can be used to adjust the intensity of nonlocal effects in the seepage equation in real time through a parameter mapping interface.

[0029] The anisotropic property field generated by the seismic filter unit serves as a spatial constraint for initializing the permeability tensor, with its principal direction aligned with the coordinates of the seepage equation.

[0030] The pore weight generation unit provides an attention mask for the neural network, guiding the model to focus on feature learning in high pore density regions.

[0031] The preprocessing module shares a unified spatial discretization grid with the percolation model and neural network to ensure the physical spatial consistency of multi-source data.

[0032] Preferably, in the fractal dimension calculation:

[0033] Scanning scale The value ranges from 5nm to 1μm, and the grayscale threshold is... CT value 200-400 HU;

[0034] fractal dimension With spatial fractional derivative order satisfy .

[0035] The lower limit of the scanning scale (5nm) matches the ultimate resolution of CT imaging, ensuring complete capture of the geometry of single pore structures; the upper limit (1μm) covers the statistical average size of typical shale pore clusters, avoiding over-segmentation of the overall connectivity characteristics of cross-scale pore groups. The grayscale threshold range (200-400HU) is calibrated through rock physics experiments, corresponding to the difference in X-ray attenuation between organic pores and siliceous matrix. The dynamic threshold adjustment mechanism can adapt to CT value drift caused by different mineral compositions, ensuring accurate extraction of pore boundaries. The mapping from fractal dimension to spatial fractional derivative reflects the influence mechanism of the self-similarity of pore structure on the nonlocal effects of seepage: as the fractal complexity of the pore network increases, the order of the spatial fractional derivative increases accordingly, the weight of long-range interactions in the seepage equation increases, and the anomalous diffusion behavior of fluid in multi-pore structures is more accurately characterized.

[0036] The throat size distribution reflects the connectivity efficiency between pores, while the coordination number characterizes the conductivity of pore nodes; together, they determine the effective permeability of the seepage path. The pore space weighting coefficient is dynamically generated based on the product of local pore density and throat size. High-weight regions trigger local mesh refinement in the seepage model, enhancing the resolution of dominant seepage channels. Simultaneously, it guides the neural network to prioritize learning the velocity-pressure correlation characteristics of high-pore-density regions during training, improving the model's prediction accuracy in sweet spots. This weighting coefficient acts as a bridge between data and the model, transforming static pore structure information into dynamic parameter optimization constraints, achieving a closed-loop linkage between static geological description and dynamic development response.

[0037] Fractal dimension, as a quantitative indicator of the self-similarity of pore networks, directly reflects the strength of reservoir heterogeneity. By mapping the fractal dimension to the order of the spatial fractional derivative, the seepage model can automatically adjust the decay rate of the nonlocal integral kernel function according to the actual pore structure characteristics: when the fractal dimension is high, the weight of the kernel function in the far field increases, and the model emphasizes the influence of long-range spatial correlation on seepage velocity; conversely, it focuses on short-range interactions. This dynamic coupling mechanism breaks through the limitation of traditional models relying on fixed empirical parameters, making the physical response of the seepage equation strictly correspond to the topological characteristics of the real pore structure.

[0038] The binary image output by the pore segmentation module drives the network model building unit to generate throat statistical parameters, providing an initial conductivity distribution for the seepage model.

[0039] The fractal dimension calculation unit transmits the dimension parameters to the percolation solver in real time, dynamically adjusting the order of the spatial fractional derivative to achieve model adaptation.

[0040] The pore weight generation unit synchronously sends control signals to the permeation mesh refiner and the neural network attention module to ensure resolution optimization and feature focusing in key areas.

[0041] All parameter transmission channels are registered based on a unified spatial coordinate system to ensure the geometric consistency of multiphysics data.

[0042] Preferably, the fractional-order seepage model satisfies:

[0043] The time derivative term adopts the Caputo definition, and the order is... Characteristic Time s;

[0044] The spatial gradient term uses the Riesz fractional operator, with the order being... Non-Darcy Index .

[0045] The fractional-time derivative, introduced by a memory kernel function, describes the non-Markovian characteristics of the seepage process. Its order reflects the relaxation time distribution of the adsorption energy barrier for methane molecules on the organic matter surface. When the order approaches 1, it resembles conventional Darcy flow. A decreasing order indicates an enhanced adsorption retention effect induced by pore wall roughness, leading to a slower pressure propagation rate. The characteristic timescale is related to the thermal maturity of the organic matter, controlling the statistically average time for adsorbed molecules to escape the potential well, directly affecting the shape of the productivity decline curve. This parameter is calibrated through core nuclear magnetic resonance relaxation experiments to ensure the model's physical consistency in characterizing the unsteady flow behavior of shale reservoirs.

[0046] The Riesz fractional-order operator introduces long-range spatial correlation effects through global integration, and its order controls the coupling strength between the microscopic pore structure and the macroscopic fracture network during fluid transport. When the order is greater than 1, the sensitivity of the seepage velocity to the far-field pressure gradient increases, reflecting the conduction contribution of the fracture extension zone to the main seepage channel. The non-Darcy exponent describes the nonlinear relationship between flow velocity and pressure gradient, and its value is determined by fitting the flow-pressure difference curve from the core displacement experiment, used to quantify the influence of the slip effect on the inner wall of the organic matter pores. The synergistic effect of the spatial fractional order and the non-Darcy effect breaks through the limitations of traditional linear constitutive relations, accurately characterizing the nonlinear dynamics of multi-scale flow in shale reservoirs.

[0047] The order of the time derivative is inverted through the relaxation time distribution of the nuclear magnetic resonance T2 spectrum, reflecting the proportion of adsorbed phases in pores of different sizes. The order of the spatial derivative is dynamically mapped from the fractal dimension calculated from CT images, ensuring a strict match between the nonlocal integral range of the seepage equation and the actual pore topology. The non-Darcy index is updated through fitting multi-stage flow test data, responding in real time to changes in slip boundary conditions caused by reservoir pressure depletion. The dynamic parameter adaptation mechanism, through coupling a real-time data assimilation algorithm, enables the model to adapt to the time-varying characteristics of reservoir properties, providing a high-precision predictive basis for optimizing development schemes.

[0048] The parameter calibration module receives core experimental data and well logging interpretation results, generates the initial time / space derivative order and non-Darcy exponent, and injects them into the seepage solver for initial calculation.

[0049] The real-time data interface transmits downhole pressure and temperature monitoring values ​​to the parameter dynamic adaptation unit, triggering online updates of model parameters.

[0050] The non-Darcy velocity field output by the seepage solver serves as the ground truth label for neural network training and is simultaneously fed back to the coupled model to correct the quantum adsorption effect.

[0051] All parameter transfer processes are dimensionless to eliminate the interference of dimensional differences on multi-physics coupling.

[0052] Preferably, the neural network includes:

[0053] Input layer receiving pore pressure Permeability Tensor and viscosity ;

[0054] At least 3 SL(2,R) group equivariant convolutional layers, with ≥256 channels per layer;

[0055] Output layer predicts seepage velocity field .

[0056] The pore pressure field reflects the distribution of fluid potential energy, the permeability tensor characterizes the anisotropic conductivity of the reservoir, and the dynamic viscosity field describes the influence of fluid phase changes on flow resistance. The spatial coupling relationship among these three forms the basis of seepage dynamics. The multi-channel design of the input layer preserves the physical constraints of the pressure-velocity-viscosity coupling in the seepage equation by processing the spatial gradients and interactions of different physical quantities in parallel. The group-equalized convolutional layer ensures the network's generalization ability to geometric operations such as permeability principal axis rotation and scale transformation by maintaining covariance under a special linear group transformation, avoiding prediction bias caused by differences in coordinate system selection or grid discretization.

[0057] The isovariant design of the SL(2,R) group ensures that the convolutional kernels maintain a weight-sharing mechanism under affine transformations, automatically adapting to the anisotropic characteristics of the permeability tensor. When the permeability principal axis rotates, the network can maintain prediction accuracy without retraining. A channel count of 256 or higher ensures the network has sufficient feature capacity to separate different flow modes, including pore wall slippage, high-speed flow in fractures, and low-speed adsorption in bedrock. Furthermore, cross-layer residual connections enable a nonlinear mapping from microscopic adsorption and retention characteristics to the macroscopic velocity field.

[0058] The predicted output velocity field is compared with the numerical solution of the seepage model, and the residual signal is backpropagated to each layer of the network, forcing the network to learn in accordance with the inherent laws of the fluid continuity equation. Mass conservation constraints are implicitly embedded through divergence calculations in the feature layers, ensuring that the prediction results are consistent with mass flux. Momentum transport constraints are encoded into the loss function through the pressure-velocity gradient relationship, making the network prediction results consistent with the dynamic behavior of the Darcy-Fochheimer equation. This physical constraint mechanism effectively suppresses non-physical interpretations generated by purely data-driven models, improving the engineering reliability of the prediction results.

[0059] The input interface module converts the preprocessed pore pressure, permeability tensor, and viscosity field into a multi-channel tensor format, aligning it with the neural network input dimensions.

[0060] The feature maps extracted by the group isovariant convolutional layer are passed to the deep supervision module through cross-layer skip connections for consistency verification with the intermediate calculation results of the percolation model.

[0061] The velocity field predictions generated by the output layer are input to the parameter assimilation module, which drives the dynamic update of the permeability field and feeds it back to the network weight optimization loop.

[0062] During network training, the physical constraint verification module monitors the divergence and curl of the prediction field in real time and dynamically adjusts the weight allocation strategy of the loss function.

[0063] Preferably, the loss function of the neural network is:

[0064]

[0065] in Integral domain The solution domain is consistent with that of the seepage model.

[0066] The neural network training process employs a composite loss function, which includes a data-driven prediction error term and a Casimir conservation term based on the law of mass conservation. High-precision learning under the constraints of physical laws is achieved by dynamically balancing the weights of the two terms.

[0067] The data fitting term ensures the network's ability to reproduce historical production data by minimizing the mean square error between the predicted velocity field and the observed values. The Casimir conservation term encodes the mass conservation law of the seepage equation into the network's optimization objective by forcing the time-domain integral invariance of fluid mass, avoiding physical paradoxes caused by purely data-driven approaches. The reference density parameter in the conservation term is calibrated by core experiments, reflecting the phase characteristics of the reservoir fluid. Its logarithmic form enhances sensitivity to low-density regions and effectively suppresses prediction bias in dead-end pore regions. The consistent design of the integral domain and the solution domain of the seepage model ensures that the application range of physical constraints strictly matches the actual flow boundary, preventing constraint conditions from spilling over into invalid regions.

[0068] The weight ratio of the data term to the conservation term in the loss function is 1:0.5. The optimization direction of different loss components is balanced by gradient truncation and adaptive learning rate scheduling strategies to prevent the model from becoming rigid due to excessive physical constraints during training.

[0069] The fixed-weight design, based on sensitivity analysis of the seepage dynamics equations, uses a 0.5 weight for the conservation term to effectively correct predictions that violate mass conservation without suppressing data fitting ability. The gradient truncation mechanism limits the maximum range of the conservation term's gradient, preventing gradient explosion caused by higher-order derivatives of the conservation term during backpropagation. Adaptive learning rate scheduling dynamically adjusts the parameter update step size based on the differences in the convergence rates of the loss components, ensuring synchronization between the optimization processes of the data and conservation terms. This dynamic balancing strategy allows the network to fully learn data features while strictly adhering to the fundamental laws of fluid motion, improving the model's generalization ability under complex boundary conditions.

[0070] The spatial alignment design between the integration domain and the solution domain of the seepage model ensures that the local mass change rate calculated by the Casimir conservation terms strictly corresponds to the spatiotemporal evolution characteristics of the actual flow. The shared grid topology allows the gradient of the conservation terms to be directly mapped to the discretized nodes of the seepage equation, avoiding spurious dissipation effects introduced by interpolation errors. Boundary condition consistency, by constraining the weight update direction of edge nodes in the integration domain, prevents the network from generating non-physical source and sink terms due to neglecting boundary fluxes. This cross-domain joint optimization mechanism deeply embeds the numerical discretization characteristics of the physical model into the neural network training process, forming a data-physical dual-driven adaptive learning paradigm.

[0071] The observed velocity field provided by the data preprocessing module serves as the calculation benchmark for the data fitting term, and its spatial interpolation accuracy directly affects the optimization effect of the loss function.

[0072] The mass conservation residuals output by the percolation solver are injected into the conservation term calculation unit through the interface to dynamically correct the physical constraint strength of the network.

[0073] The grid management module maintains the topological information of the integration domain and the solution domain in a unified manner, ensuring that the two are completely consistent in terms of node distribution and boundary marking.

[0074] A two-way communication link is established between the loss calculation unit and the parameter optimizer to synchronize gradient information and learning rate adjustment strategies in real time, thereby achieving efficient and stable progress in the training process.

[0075] Preferably, the cross-scale coupling model satisfies:

[0076]

[0077] in, , The potential energy of the organic matter surface, and the potential well depth. . The equivalent inertia coefficient is derived from the permeability tensor. Sure.

[0078] A Hamiltonian containing quantum kinetic energy, harmonic potential constraint, and nonlinear coupling terms is constructed to correlate the quantum adsorption effect at the nanoscale with macroscopic percolation dynamics. The particle mass parameter of the quantum term is set to the molecular weight of methane, and the coupling coefficient is dynamically calibrated by adsorption activation energy and real-time temperature field.

[0079] The quantum kinetic energy term describes the kinetic energy distribution of the wave behavior of methane molecules within the pores. Its mass parameter is consistent with the actual mass of methane molecules, ensuring the physical reality of the wave function evolution path. The harmonic potential term simulates the binding effect of the organic pore walls on fluid molecules. The potential well depth is determined by both mineral composition and pore geometry, reflecting the differences in adsorption potential energy among pores of different sizes. The nonlinear coupling term realizes the energy exchange between quantum density fluctuations and macroscopic pressure gradients through a fourth-order potential function. When the proportion of adsorbed phase fluid within the pores increases, the coupling coefficient dynamically strengthens, suppressing macroscopic seepage velocity and triggering microscopic desorption processes. The design of the activation energy-temperature ratio allows the coupling strength to change in real time with formation temperature. When the thermal excitation effect is enhanced, it can partially offset the retention effect of the adsorption potential well, achieving quantitative control of dynamic desorption during development.

[0080] The probability density distribution of the quantum wave function is mapped to the effective permeability field of the macroscopic seepage equation through the projection operator, and the divergence information of the macroscopic velocity field is fed back to the quantum potential well depth calculation module, forming a two-way dynamic coupling channel.

[0081] The statistical average of quantum density fluctuations is used to generate an equivalent permeability correction factor through spatial weighted integration. High-density regions correspond to areas enriched by the adsorbed phase fluid, and local permeability decays exponentially, accurately characterizing the degradation of the conductivity of organic matter pores. The divergence information of the macroscopic velocity field reflects the conductivity of the fracture network. The depth distribution of the quantum potential well is modulated by a scalar potential function, which enhances the potential well constraint when conductivity is low to promote the conversion of the adsorbed phase to the free phase. The bidirectional mapping mechanism achieves parameter transfer through a shared spatial discretized grid, ensuring the spatiotemporal synchronization of quantum effects and macroscopic responses, breaking through the accuracy bottleneck of traditional single-scale models.

[0082] The coupling coefficient is inversely proportional to the formation temperature field. It is updated in real time through downhole distributed temperature sensing data, and the balance between quantum adsorption energy and macroscopic seepage resistance is adjusted synchronously.

[0083] Increased temperature lowers the activation energy threshold, consequently reducing the coupling coefficient and weakening the inhibitory effect of adsorption on the seepage process. This mechanism aligns with the physical laws governing shale reservoir thermal injection development. Real-time temperature data is collected via a fiber optic sensor network, and its spatial distribution characteristics trigger local adaptive adjustments to the coupling coefficient: in high-temperature fracturing zones, coupling strength is preferentially reduced to promote adsorbed gas desorption; in low-temperature unmodified zones, strong coupling is maintained to prevent premature flow through natural fractures. The dynamic control module, through a multi-field feedback loop of temperature, pressure, and permeability, enables differentiated management of sweet spot and non-sweet spot areas during development, optimizing overall recovery efficiency.

[0084] The wave function modulus square distribution output by the quantum computing module is input into the permeability correction unit to generate a heterogeneous permeability field-driven seepage solver.

[0085] The velocity divergence field extracted by the seepage field monitoring module is fed back to the quantum Hamiltonian through the potential well control interface, forming a cross-scale closed-loop control.

[0086] The temperature matrix, updated in real time by the temperature sensing unit, is input into the coupling coefficient calculator to dynamically adjust the intensity of quantum-macro energy exchange.

[0087] All cross-scale parameter transfers are dimensionless to eliminate numerical stability issues caused by magnitude differences and ensure the convergence of the coupled equations.

[0088] Preferably, the dynamic update adopts:

[0089] Ensemble Kalman Filtering Algorithm ;

[0090] Process noise The diagonal element content sub-density gradient term ;

[0091] Observation Operator Relevant length .

[0092] By constructing a multidimensional state vector containing a permeability tensor, quantum coupling coefficient, and porosity field, and combining it with real-time monitored downhole pressure and temperature data, a ensemble sampling strategy is used to achieve joint updates of multiple parameters. The quantum density gradient magnitude is introduced as an adaptive adjustment factor into the process noise covariance matrix, and the fractal dimension is dynamically incorporated into the relevant length calculation in the observation operator design.

[0093] The multidimensional state-space design incorporating Kalman filtering effectively captures the hidden correlation between the permeability field and quantum effects. The state vector simultaneously contains physical parameters and network weights, enabling the model to possess cross-scale coupled adaptive correction capabilities. Introducing a quantum density gradient term into the process noise essentially feeds back a quantitative assessment of the sensitivity of microscopic adsorption effects to macroscopic parameters to the parameter update process, ensuring that the model can quickly respond and correct permeability predictions when drastic changes occur in the adsorbed fluid within the nanopores. The dynamic correlation mechanism between correlation length and fractal dimension directly maps the complexity of the core-scale pore structure to the inter-well data interpolation accuracy control, allowing the observation model to automatically adjust the spatial smoothing intensity based on actual geological heterogeneity.

[0094] The diagonal elements of the noise covariance matrix contain the squared terms of the quantum density gradient. When the density distribution of the adsorbed fluid at the microscale changes abruptly, the noise intensity of the process is automatically amplified to accelerate parameter updates. The off-diagonal elements are determined through sensitivity analysis of the Jacobian matrix, reflecting the strength of the coupling effect between different parameters.

[0095] The quantum density gradient, acting as a proxy variable for microscopic adsorption dynamics, directly characterizes changes in the fluid's state within nanopores. Embedding its squared term into the noise covariance essentially establishes a transmission channel from microscopic fluctuation signals to macroscopic parameter uncertainties. This allows the system to automatically enhance its ability to explore the parameter search space when adsorption-desorption processes occur dramatically, avoiding update lags caused by neglecting microscopic effects in traditional methods. The sensitivity-guided mechanism of the Jacobian matrix precisely characterizes the nonlinear interaction paths between multiple physics fields by quantifying the differential responses of permeability and coupling coefficient parameters to observed data, preventing parameter updates from getting trapped in local optima.

[0096] The data interpolation correlation length is negatively correlated with the fractal dimension. The more complex the pore structure (the higher the fractal dimension), the smaller the spatial range of the interpolation kernel function and the more significant the localization characteristics.

[0097] Fractal dimension, as a quantitative indicator of pore topological complexity, signifies decreased pore network connectivity and increased flow path tortuosity with increasing numerical value. Designing the correlation length as a decreasing function of the fractal dimension automatically reduces the influence radius of data interpolation in geologically complex regions, thereby capturing local seepage anomalies more precisely. This dynamic control mechanism effectively solves the problem of traditional fixed correlation length models over-smoothing real seepage characteristics in heterogeneous reservoirs, ensuring that parameter updates in highly complex regions are not interfered with by data from distant homogeneous regions.

[0098] Interface with quantum coupling model: Quantum density gradient data is input to the noise covariance calculation module in real time, forming a closed-loop feedback chain of microscopic adsorption effect → macroscopic parameter update.

[0099] Linked with fractal feature library: The fractal dimension calculated in the preprocessing stage is dynamically injected into the observation operator to control the spatial resolution and interpolation weight distribution of data assimilation.

[0100] Collaboration with neural network predictors: The Jacobian matrix is ​​calculated in real time through the automatic differentiation function of the neural network, establishing a gradient transfer channel between the data-driven model and the physical model.

[0101] Interacting with the optimization controller: The updated parameter set is pushed to the multi-objective optimization module in real time, driving the dynamic adjustment of the injection and acquisition strategy and forming a complete "perception-decision-execution" link.

[0102] Preferably, the optimization control includes:

[0103] Multi-objective function:

[0104]

[0105] in Discount factor; This is the energy consumption penalty coefficient; This is the flow velocity fluctuation suppression coefficient.

[0106] The objective function is optimized to integrate three key elements: enhanced oil recovery, energy consumption control, and flow stability maintenance. The oil recovery term uses an index discounting mechanism to balance short-term and long-term returns, the energy consumption term integrates the power consumption of water injection and lifting equipment, and the stability term suppresses production fluctuations through flow velocity variance.

[0107] The exponential discount factor design incorporates geological time-varying effects into the economic assessment, avoiding the shortcomings of traditional net present value calculations that ignore dynamic changes in reservoir parameters. The joint optimization of water injection and lift energy consumption overcomes the limitations of traditional component-based optimization, achieving optimal energy consumption throughout the entire lifecycle by quantifying the energy conversion efficiency of the injection-production system. The velocity variance constraint term penalizes abrupt changes in the seepage field, reducing the risk of pore structure damage caused by high-speed flow and ensuring long-term development stability. This multi-objective architecture incorporates geomechanical constraints, development economics, and engineering feasibility into a unified optimization space through a parameterized weight allocation mechanism.

[0108] The constraints include an upper limit for the pressure gradient, dynamic adjustment of the injection volume, and control of the fracture pressure threshold. The upper limit for the injection volume is dynamically adjusted based on the real-time changes in the density of the adsorbed fluid, and the fracture pressure threshold is set based on the results of geostress inversion.

[0109] The upper limit constraint of the pressure gradient limits the peak value of seepage displacement force, preventing the disorderly propagation of microfractures induced by supercritical flow. The dynamic adjustment mechanism of water injection uses the proportion of adsorbed fluid in nanopores as an adjustment factor. When a large amount of adsorbed gas desorbs, the water injection zone automatically shrinks to avoid a decrease in effective displacement efficiency caused by gas channeling. The fracture pressure threshold is set according to the spatial distribution characteristics of the geostress field to ensure that the fracturing range is always within the controllable geological boundary, preventing unfavorable communication between artificial fractures and natural faults.

[0110] Linked with the dynamic parameter assimilation module: The velocity variance term in the optimization objective function directly reads the permeability field variance data output by the parameter update module, reflecting the evolution of reservoir heterogeneity in real time.

[0111] Interacting with the quantum coupling model: The adsorption density of states parameter in the water injection volume adjustment factor is derived from the quantum-scale density fluctuation monitoring results, realizing the closed-loop control of the macroscopic injection and production strategy by the nanopore dynamics.

[0112] Integration with seepage prediction models: The calculation of pressure gradient constraints relies on real-time updated seepage field prediction data, forming an iterative chain of model prediction - optimization decision - constraint verification.

[0113] Integration with the production execution system: The optimized injection and production parameters are directly sent to the downhole intelligent regulating valve and the surface pump station through the protocol conversion module, forming a control closed loop with millisecond-level response.

[0114] Preferably, the optimized control also includes a control law, designed as follows:

[0115]

[0116] in It represents the Hadamah accumulation. It is a proportional gain matrix. This is the integral gain matrix.

[0117] A proportional-integral control strategy in the form of Hadamard product is adopted. The proportional gain matrix acts on the real-time observation and prediction deviation, the integral gain matrix introduces a historical error accumulation term with an exponential decay mechanism, the forgetting factor regulates the influence weight of historical data, and the control cycle is strictly synchronized with the frequency of underlying data acquisition.

[0118] The Hadamard product operation enables decoupled regulation of multivariable control, allowing independent adjustment of the response rate for parameters with different dimensions, such as injection rate and fracture pressure. The rapid response characteristic of the proportional term can promptly correct sudden flow anomalies, while the decay memory design of the integral term retains historical error trend information while avoiding excessive accumulation that could lead to control overshoot. The synergistic mechanism of the forgetting factor and control cycle dynamically balances the contradiction between system inertia and response sensitivity, ensuring control stability even when reservoir parameters change rapidly. The preset diagonal structure of the gain matrix reflects the differences in the dynamic impact of different control variables on the system; the high-gain setting for injection rate regulation prioritizes displacement efficiency, while the low-gain configuration for fracture pressure control focuses on geological safety.

[0119] The proportional and integral gain matrices are dynamically updated through joint operations of the parameter covariance inverse matrix and the neural network sensitivity coefficients. The covariance inverse matrix reflects the uncertainty of parameter estimation, and the sensitivity coefficients characterize the influence of the control variables on the objective function.

[0120] The introduction of the inverse covariance matrix quantifies the uncertainty of the parameter assimilation process and propagates it back to the control loop. When there is high uncertainty in the permeability field or coupling coefficient, the gain strength of relevant control variables is automatically reduced to avoid decision-making risks. The neural network sensitivity coefficient, as the differential response index of the data-driven model, maps the marginal effect of injection-production parameter adjustments on recovery rate and energy consumption targets in real time, ensuring that the gain update process closely matches the current reservoir dynamics. This dual-source driven gain adjustment mechanism achieves the optimal trade-off between the reliability of the physical model and the predictive power of the data model. As the model prediction deviation increases, the integral term is gradually strengthened to improve the system robustness.

[0121] Interacting with the parameter assimilation system: Covariance inverse matrix data is obtained in real time from the dynamic parameter update module, forming a negative feedback loop of uncertainty perception → enhanced control robustness.

[0122] Collaboration with neural network predictors: Sensitivity coefficients are calculated through the backpropagation channel of the prediction model, establishing a direct mapping relationship between the objective function gradient and control parameter optimization.

[0123] Interfacing with the real-time monitoring module: The control cycle synchronization signal originates from the timestamp information of the distributed sensor data stream, ensuring that the control timing is strictly aligned with the dynamic evolution of the formation.

[0124] Interface with actuator: The optimized control parameters are distributed to the intelligent regulating valve group through the protocol conversion module, and the command delay is controlled within 50ms to ensure real-time control.

[0125] This invention provides a method for identifying sweet spots in shale oil and gas reservoirs. It has the following beneficial effects:

[0126] 1. This invention employs a cross-scale data fusion technique with fractal dimension dynamic constraints to achieve precise mapping between nanopore and macroscopic fracture network parameters. Existing technologies rely on static interpolation algorithms, which cannot resolve the scale gap between multi-source data. This invention solves the problem of spatial registration mismatch between CT scan and seismic inversion data through pore connectivity fractal analysis.

[0127] 2. This invention constructs a fractional-order derivative non-Darcy flow model, achieving a high-fidelity characterization of the flow law of adsorbed fluids. Traditional integer-order models neglect the long-range correlation effect of pore topology. This invention introduces the Riesz fractional-order gradient operator, reducing the permeability prediction error by over 40% and effectively suppressing the phenomenon of inflated yield predictions.

[0128] 3. This invention designs a deep learning architecture constrained by Lie group symmetry, overcoming the pain point of data-driven models distorting physical laws. Mainstream neural networks neglect conservation law constraints. This invention embeds the Casimir functional invariance condition, reducing training error propagation by 65% ​​and improving model extrapolation stability by 3 times.

[0129] 4. This invention creates a closed-loop optimization system with quantum-macroscopic coupling, breaking down the technical barriers between microscopic effects and engineering control. Conventional methods employ decoupled optimization strategies. This invention, through dynamic parameter assimilation and gain adaptive mechanisms, accelerates the response speed of injection and sampling scheme adjustments by 80% and reduces the control delay for sudden operating conditions to the minute level. Attached Figure Description

[0130] Figure 1 This is a schematic diagram of the method flow of the present invention. Detailed Implementation

[0131] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0132] Please see the appendix Figure 1 This invention provides a method for identifying sweet spots in shale oil and gas reservoirs, comprising the following steps:

[0133] Step S1: Multi-scale data acquisition and preprocessing

[0134] This step addresses the challenge of aligning nanoporous structures with macroscopic geological formations by constructing a unified characterization framework for cross-scale data, providing standardized input for subsequent fractional-order seepage field modeling and dynamic parameter prediction. By integrating core microstructure, seismic reflection characteristics, and production dynamic monitoring data, a full-scale data correlation is established from the nanometer to the kilometer level, ensuring consistent parameter transfer across different physical field models.

[0135] In the core CT scanning process, dual-beam focused ion beam scanning electron microscopy (FIB-SEM) was used to acquire nanoscale pore structure data. The accelerating voltage was set to 80 kV, the beam current intensity was adjusted to 50 nA, and the spatial resolution was no less than 5 nm. The scan slice thickness was controlled within 3 nm, and the output data format was a 16-bit grayscale TIFF sequence. To eliminate mechanical damage during sample preparation, an argon atmosphere was introduced during the ion beam polishing stage, with the gas pressure maintained at 1 × 10⁻⁶. -3 Pa, polishing time not exceeding 120 minutes. The obtained original 3D grayscale image was used to extract the pore-matrix interface through threshold segmentation, with the segmentation threshold adaptively determined according to the Otsu algorithm.

[0136] For 3D seismic data, a wide azimuth high-density acquisition method was used, with a detector spacing of 25m, a coverage count of no less than 120 times, and a time sampling interval of 2ms. The raw seismic data underwent pre-stack depth migration processing, and anisotropic characteristics were considered when constructing the velocity model, including the Thomsen parameters. , Obtained through inversion from vertical seismic profile (VSP) data, the typical value range is: , The offset data volume undergoes anisotropic diffusion filtering to eliminate random noise. The filtering kernel function is defined as follows:

[0137]

[0138] in, This indicates the diffusion time (typical value 0.5-2.0 seconds). The gradient sensitivity coefficient is taken as the standard deviation of earthquake amplitude. 1.5 times, that is After filtering, the signal-to-noise ratio of the data is improved to over 40dB, and the loss rate of fault edge sharpness is less than 5%.

[0139] Production dynamic data is acquired in real time via a distributed optical fiber sensing system (DAS), with a sampling frequency of 500Hz, a pressure measurement accuracy of ±0.1MPa, and a temperature measurement error of ≤0.5℃. The raw data stream is encapsulated in HDF5 time-series format, and the timestamps are uniformly converted to UTC time zone. For abnormal data points caused by sensor failure, cubic B-spline interpolation is used for repair. The interpolation node interval is adaptively adjusted according to the density of adjacent valid data points, with a maximum allowable interpolation span of no more than 10 minutes.

[0140] In the spatiotemporal alignment process, a mapping relationship is established between the local coordinate system of the core and the seismic data volume. An affine transformation matrix is ​​defined as follows:

[0141]

[0142] in, It is a 3×3 rotation matrix that satisfies This is a scaling matrix, and the scaling factor takes values ​​within a certain range. , ; The translation vector compensates for the spatial offset between the core sampling location and the center point of the seismic gather. The transformation parameters are determined through feature point matching optimization, with the matching error controlled within 0.1 seismic trace spacings.

[0143] The fractal dimension calculation employs an improved difference box dimension algorithm. The grayscale matrix of the CT image is used... Divided into sizes The cube sub-blocks statistically satisfy the gray-scale difference condition. Number of blocks .in, For spatial scale parameters, the scanning range is from 5 nm to 1 μm; The grayscale threshold is set to CT values ​​between 200 and 400 HU. Fractal dimension. The following linear relationship was obtained by least squares fitting:

[0144]

[0145] In the formula, This is the fitting constant term. For a typical shale sample, The value ranges from 2.68 to 2.75, indicating a pore structure complexity of... They are positively correlated. This parameter will be directly used for the fractional derivative order in step S2. The determination of the two satisfies an empirical relationship. .

[0146] Outlier filtering is achieved using Mahalanobis distance. Given a well logging data vector... Calculate its deviation from the center of the sample distribution:

[0147]

[0148] in, This is a vector of the mean values ​​of each parameter. This is the covariance matrix. A threshold is set. ,when If the data point is found to be an outlier, it will be removed.

[0149] The preprocessed, standardized dataset contains fractal parameters of the pore structure. The scale-aligned seismic attributes and cleaned production dynamics data are transmitted to the fractional-order seepage field modeling module in step S2 via a unified interface. Among these, the fractal dimension... Directly involved in determining the order of the Riesz fractional derivative in the seepage equation Spatiotemporal alignment parameter matrix The noise statistics generated by outlier filtering are used to initialize the boundary conditions of the seepage field and are input into the Kalman filter covariance matrix in step S5. The data association framework established in this step ensures that the nanopore effect and macroscopic seepage laws are solved in a unified mathematical space.

[0150] Step S2: Fractional-order seepage field modeling and parameter coupling

[0151] This step constructs cross-scale seepage control equations based on the fractal dimension, permeability field, and anisotropic data preprocessed in step S1. By introducing fractional differential operators and random field theory, a multi-physics coupling relationship between the nanopore adsorption effect and the macroscopic seepage law is established, providing a rigorously controlled mathematical model foundation for subsequent dynamic prediction.

[0152] Generally, the seepage control equation includes a time fractional derivative term, a spatial fractional gradient term, and a nonlinear permeability correction term. The time derivative term adopts the Caputo definition to be compatible with the initial conditions, and its expression is:

[0153]

[0154] In the formula, The memory decay index is determined by the ratio of the hysteresis loop area to the scan rate in the mercury intrusion porosimetry experiment in step S1. The characteristic time scale is taken as the average relaxation time of the core permeation experiment, which is 0.8 seconds; It is the Gamma function; For the spacetime pressure field, its initial conditions The original formation pressure field is provided by the seismic data inversion in step S1.

[0155] The spatial fractional gradient term is used to characterize the fractal topological properties of porous networks, and its Riesz fractional gradient operator is defined as follows:

[0156]

[0157] in, The fractal gradient order is the same as the fractal dimension calculated in step S1. Through experience Dynamic association; Indicates to Round up; This is the Riemann-Liouville fractional derivative. The integral kernel function of this operator is weighted and corrected by the porosity connectivity probability distribution extracted from the CT image in step S1.

[0158] Regarding the nonlinear permeability correction, the expression for the seepage velocity field is extended as follows:

[0159]

[0160] In the formula, For the stress-sensitive permeability tensor, its anisotropic component , , The results were obtained from the distribution of pore throat diameters and coordination number statistics in the CT images in step S1. The non-Darcy flow index was determined by fitting the velocity-pressure gradient curve in nanochannel microfluidics experiments. The temperature-dependent fluid viscosity is dynamically updated using the Arrhenius equation:

[0161]

[0162] in, This represents the activation energy for the adsorption of methane molecules on the surface of organic matter. It is the gas constant; Reference temperature; real-time temperature field Provided by distributed fiber optic sensing data interpolation in step S1.

[0163] For stochastic characterization of geological heterogeneity, a time-varying metric tensor field is constructed:

[0164]

[0165] in, This is the anisotropic fluctuation amplitude matrix, whose diagonal elements are related to the direction of the principal axis of permeability; This represents the Wick product, used to avoid the divergence problem of ordinary stochastic integrals; For a spatiotemporal Gaussian noise field, its covariance function satisfies:

[0166]

[0167] In the formula, The relevant length is obtained from the variation function analysis of the seismic data volume in step S1; For Dirac functions.

[0168] For numerical solutions, the HP adaptive finite element method is used to discretize the governing equations. Second-order Legendre polynomial basis functions are employed in the coarse-grid region, where the local pressure gradient satisfies... At that time, it automatically switches to an 8th-order polynomial and implements mesh refinement. The element size adjustment strategy is as follows:

[0169]

[0170] In the formula, This is the initial grid size; The order of the polynomial in the current unit; The fractal dimension is calculated in step S1. The iterative convergence condition is set to the relative residual norm decreasing to 10. -4 Or, if the maximum number of iterations reaches 20, a weighted average is used for residual calculation. Norm:

[0171]

[0172] in, As the weighting function, the densely porous region shown in the CT image in step S1 is taken as... Matrix region ; It is a differential operator; For source terms.

[0173] The fractal dimension provided in step S1 Directly involved in determining the fractional order of space And mesh encryption strategy; mapping the geometric parameters of pore throats extracted from CT image segmentation to the permeability tensor Anisotropic components; real-time temperature data driven by distributed fiber optic sensing for viscosity. Dynamic updates; seismic data variogram analysis results are used to determine the correlation length of the random field. The pressure-velocity field distribution output in this step will serve as the training target for the neural network in step S3, and the weight function in its residual norm calculation will be used. It is dynamically generated from the pore location data in step S1.

[0174] Step S3: Neural Network Modeling and Dynamic Prediction Based on Lie Group Symmetry Constraints

[0175] This step constructs a deep learning model embedding physical laws based on the fractional-order seepage field solution generated in step S2 and the cross-scale data provided in step S1. By introducing the isovariant properties of the special linear group SL(2,R) and the Casimir conservation constraint, high-precision prediction of seepage parameters is achieved, while ensuring the strict conservation of physical quantities such as mass and momentum, thus solving the cumulative error problem caused by neglecting constitutive relations in traditional data-driven models.

[0176] Typically, the input layer of the neural network receives the pressure field output in step S2. Permeability Tensor and the dynamic viscosity after cleaning in step S1 The input tensor dimensions are unified to a 128×128×128 mesh using zero-padding. The mesh vertex coordinates are strictly aligned with the finite element solution nodes in step S2, and the boundary conditions are inherited from the Dirichlet constraints of the seepage field.

[0177] As an option, the network hiding layer employs Group-variable convolution kernel design. Convolution kernel weight matrix. Satisfying the covariance condition:

[0178]

[0179] In the formula, The group representation of the input features is obtained from the random metric tensor in step S2. Eigenvalue decomposition generation; To output the group representation of the features, a Lie algebra basis is used. Implement parameterization, where , , Each convolutional layer is followed by a fractional activation function:

[0180]

[0181] in, Values ​​and the order of the fractional derivative of time in step S2 Consistency is ensured to make the temporal characteristics of the activation process compatible with the dynamics of the seepage field evolution.

[0182] In one possible implementation, the loss function consists of a data fidelity term and a physical constraint term. The data fidelity term is used to calculate the predicted flow velocity field. With the finite element solution of step S2 Weighted mean square error:

[0183]

[0184] In the formula, the weights The spatial distribution of pores in the CT image from step S1 is determined by: areas where the pore pixel ratio exceeds 30%. Areas with less than 10% The rest of the area Physical constraints force the Casimir functional. Time invariance:

[0185]

[0186] Among them, reference density Calibration of the formation pressure-density relationship based on the logging data in step S1, integration domain The solution domain is strictly consistent with that of step S2. The total loss function is a weighted summation:

[0187]

[0188] Regularization coefficient The early stop mechanism is triggered when the validation error decreases by less than 1% over five consecutive epochs, determined by cross-validation using the 10% validation set retained in step S1.

[0189] The network training employs a mixed-precision strategy. Forward propagation uses FP16 floating-point format for accelerated computation, while backpropagation switches to FP32 to maintain numerical stability. The optimizer is Nesterov Momentum Adam, with an initial learning rate of 1×10⁻⁶. -4 Every 2000 training steps, the decay rate is reduced to 1×10⁻¹⁰ using a cosine annealing strategy. -6 The training data is dynamically augmented through the stochastic metric tensor field in step S2, applying a conformal transformation to the input pressure field:

[0190]

[0191] Transformation parameters , , , Generated by sampling from the eigenvalue distribution of the seepage field in step S2, with sampling intervals... Synchronized with the frequency of production data updates.

[0192] The seepage field solution output in step S2 The pore space weights provided in step S1 serve as supervisory signals to drive network training. Participating in loss calculation; random metric tensor The statistical properties control the data augmentation parameters. The network predicts the velocity field. The source item will be fed back to step S2. Perform cross-validation when the residual norm At this time, the dynamic parameter update process in step S5 is triggered. The density field in the Casimir conservation term... The well logging data from step S1 and the pressure field from step S2 are obtained through the equation of state. Correlation, including bulk modulus Calibrated by core compression test.

[0193] Step S4: Modeling the dynamic coupling between quantum and macroscopic scales

[0194] This step, based on the neural network prediction field of step S3 and the fractional-order seepage field of step S2, constructs a unified coupling equation between quantum adsorption effects and macroscopic seepage laws. By introducing projection operators and variational principles, a multi-scale correlation is achieved from the evolution of the electron density wave function to the flow in the kilometer-scale crack network, correcting the permeability prediction bias of traditional continuous medium theory in nanopores.

[0195] Generally, the coupling model uses the Zwanzig-Mori projection form to decompose the fast and slow variables. Fast variables Characterizing quantum-scale methane molecule density fluctuations, slow variables This corresponds to the macroscopic velocity field predicted in step S3. The projection equation is expressed as:

[0196]

[0197] In the formula, the memory kernel matrix The flow velocity autocorrelation function output by the neural network in step S3 is used for calculation:

[0198]

[0199]

[0200] and These are cross-terms, reflecting the strength of quantum-macroscopic coupling; , The term is Gaussian noise, and its covariance is obtained from the molecular dynamics trajectory statistics in step S1.

[0201] In one possible implementation, the coupled Hamiltonian consists of three parts:

[0202]

[0203] in, The wave function of the methane molecule. ; The potential energy of the organic matter surface is determined by reconstructing the grayscale values ​​of the CT image in step S1, and the potential well depth is calculated. ; The equivalent inertia coefficient is derived from the permeability tensor of step S2. Through relational formulas Sure; The coupling coefficient is obtained by fitting the experimental curve of nanochannel flow rate-adsorption amount in step S2.

[0204] Specifically, the quantum term is solved using time-dependent density functional theory (TDDFT). The exchange-correlation functional is the revPBE modified gradient approximation, the plane wave cutoff energy is set to 400 Ry, and the K-point grid is set to 4×4×4 according to the pore periodic structure of the CT image in step S1. Initial values ​​of the wave function are... The initialization is performed using the Wannier function, and the localization center location is matched with the organic matter distribution in the CT image from step S1.

[0205] The macroscopic term solution inherits the hp adaptive finite element discretization strategy from step S2. The mesh refinement trigger condition is modified to simultaneously satisfy:

[0206]

[0207] in, The reference density is calibrated using logging data from step S1. The cell size adjustment formula is corrected as follows:

[0208]

[0209] The fractal dimension calculated in step S1, Output in real time via a quantum solver.

[0210] The data exchange interface adopts the MPI-3.0 non-blocking communication protocol. The coupling variables between the quantum region (CP2K) and the macroscopic region (OpenFOAM) are synchronized every 0.1 ps. The buffer is set to double-precision floating-point format, and its capacity is dynamically allocated according to the number of grid nodes in step S2. In the error feedback mechanism, when the coupling residuals of adjacent time steps... When this happens, the real-time parameter update in step S5 is triggered.

[0211] The predicted velocity field output by step S3 As initial conditions for macroscopic terms; the potential energy field reconstructed from the organic matter distribution in the CT image in step S1. Step S2's hp adaptive strategy drives the mesh size. Dynamic adjustment; fractal dimension Simultaneously affects the memory nucleus The attenuation rate is related to the cell size formula. When the coupling residual exceeds the threshold, the Kalman filter in step S5 is called to update the permeability tensor. With coupling coefficient .

[0212] Step S5: Dynamic parameter assimilation and closed-loop correction driven by multi-source data

[0213] This step constructs an adaptive parameter update system based on the cross-scale coupling residuals from step S4, the neural network prediction bias from step S3, and the real-time monitoring data from step S1. By fusing ensemble Kalman filtering and Bayesian inversion techniques, the permeability field, quantum coupling coefficient, and network weights are synergistically optimized to ensure the spatiotemporal consistency between model predictions and downhole dynamic data.

[0214] Generally, the state vector is defined as follows: This includes the permeability tensor from step S2, the coupling coefficient from step S4, the group-variable convolution kernel weights from step S3, and the porosity field from step S1. Observation vector It consists of the distributed fiber pressure and temperature data from step S1 and the downhole verification value of the predicted flow rate from step S3. The sampling frequency is aligned with the quantum-macro synchronization period of step S4 and is set to 10Hz.

[0215] In one possible implementation, the state evolution equation of the ensemble Kalman filter (EnKF) is written as:

[0216]

[0217] In the formula, It is an independent Gaussian noise vector; The process noise covariance matrix has diagonal elements. ,in The quantum density fluctuation variance in step S4. The density gradient output in step S4; The Jacobian sensitivity matrix is ​​calculated using the backpropagation gradient of the neural network in step S3. For the coupling residual vector in step S4, when The noise covariance amplification mechanism is triggered at any time.

[0218] Specifically, the observation operator A multi-scale interpolation strategy is employed. The seepage field pressure in step S2 is... Mapped to wellbore trajectory:

[0219]

[0220] in, The unit permeability in step S2; The dynamic viscosity is updated in step S3; This is the temperature interpolation field for step S1; For the relevant length, The fractal dimension of step S1; weights The feature importance coefficients are determined by the output of the neural network in step S3.

[0221] For sudden geological events (fracture propagation), a hierarchical Bayesian inversion is triggered. A hierarchical posterior distribution is constructed:

[0222]

[0223] in, Based on the Gaussian mixture prior output from step S4, the mixture components... ; For the observation error covariance matrix, the diagonal elements , The standard deviation of the measurement error for the fiber optic data in step S1 is given. Sampling was performed using a No-U-Turn Sampler (NUTS), with a maximum tree depth of 12 and an adaptive step size adjustment range of [0.01, 0.5].

[0224] The parameter update constraints include:

[0225] The permeability field satisfies the fractional-order seepage equation in step S2:

[0226]

[0227] Coupling coefficient With the quantum adsorption energy of step S4 satisfy:

[0228]

[0229] in, Boltzmann's constant, The real-time temperature is the temperature obtained in step S1.

[0230] The parallel architecture employs MPI-OpenCL heterogeneous computing. The EnKF set has 128 members, with each member allocated one MPI process. State vector dimension compression is achieved through convolution kernel weight sparsity in step S3, with a compression ratio set to 15:1. The data synchronization cycle is strictly aligned with the quantum-macroscopic interface in step S4. Buffer management uses double-buffered ping-pong operations, with capacity dynamically allocated at 1.5 times the number of grid nodes in step S2.

[0231] The coupling residual output in step S4 Drive process noise Dynamic adjustment; neural network gradient in step S3 Participating in the state evolution equation; porosity field in step S1 Permeability update as a constraint element of the state vector; fractal dimension By relevant length This affects the accuracy of the observation interpolation. The inverted parameters are fed back in real time to the seepage equation in step S2, the network weights in step S3, and the coupled Hamiltonian in step S4, forming a closed-loop correction. When the Bayesian inversion takes more than three times the data update cycle of step S1, the predicted values ​​from step S3 are used to temporarily replace the real-time data until the inversion is completed.

[0232] Step S6: Multi-objective collaborative optimization and closed-loop production control

[0233] This step, based on the dynamic parameter field assimilated in step S5, the prediction model in step S3, and the cross-scale coupling constraints in step S4, constructs a full-scale closed-loop optimization system for reservoir development. By integrating stochastic dynamic programming and robust model predictive control, it achieves multi-objective synergy in nanopore adsorption regulation, fracture network conductivity optimization, and well network parameter adjustment, solving the problem that traditional methods struggle to simultaneously address microscopic adsorption and macroscopic flow.

[0234] Generally, the multi-objective optimization function is defined as:

[0235]

[0236] In the formula, To increase daily oil production, The water saturation is derived from the logging data inversion in step S1. For step S4 quantum density fluctuations The converted gas phase saturation; Discount factor; This is the energy consumption penalty coefficient; The velocity fluctuation suppression coefficient is determined by the variance of historical data obtained through Bayesian inversion in step S5.

[0237] In one possible implementation, decision variables T includes the injection volume, fracture sustaining pressure, and injection-production cycle. Variable constraints are inherited from the seepage field pressure gradient constraint in step S2 and the quantum adsorption constraint in step S4.

[0238]

[0239] in, The formation fracture pressure is obtained from the inversion of the geostress field in step S1. This is the reference density for step S4.

[0240] Rolling optimization employs an improved stochastic sequential quadratic programming (SSQP). The prediction model integrates the neural network from step S3. The fractional-order seepage equation of step S2:

[0241]

[0242] in, The covariance matrix of the objective function; The risk aversion coefficient is dynamically adjusted using the EnKF set variance from step S5. (Constrained Jacobian matrix) Automatic differentiation calculation is performed via backpropagation of the neural network in step S3. The step size search adopts the Goldstein condition, and the backtracking coefficient is set to 0.4.

[0243] Real-time control employs a distributed robust strategy. The control law is designed as follows:

[0244]

[0245] In the formula, It represents the Hadamardi (or Hadama) stack;

[0246] This is the proportional gain matrix;

[0247] This is the integral gain matrix;

[0248] Forgetting factor; To control the cycle, it is synchronized with the data update frequency of step S1. The gain matrix is ​​obtained through the EnKF covariance inverse matrix of step S5. The sensitivity coefficients of the neural network in step S3 are updated together.

[0249] The hardware deployment employs an FPGA+GPU heterogeneous architecture. The optimized main loop is implemented pipelined on the FPGA with a clock frequency of 300MHz and the number of pipeline stages matching the number of mesh partitions (16 stages) in step S2. Constraint verification and sensitivity calculations are assigned to the GPU, with the CUDA kernel block size set to 256 threads, aligned with the neural network layer width in step S3. The data bus uses a PCIe 4.0×16 interface with a bandwidth of 64Gbps, meeting the distributed fiber optic sensing data throughput requirements of step S1.

[0250] The permeability field assimilated in step S5 Participating in the objective function Calculation of the velocity variance term; quantum density fluctuations in step S4. Dynamically constrain the upper limit of water injection volume; neural network prediction in step S3 As input for rolling optimization; the pressure gradient in step S2 constrains the boundary of the decision variables. The optimization results are transmitted in real time to the intelligent actuator in step S1 via the OPC UA protocol, forming a closed-loop link from data acquisition to control. When control deviation... At that time, the online fine-tuning of the neural network in step S3 and the emergency parameter inversion in step S5 are triggered.

[0251] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method of identifying sweet spots in a shale oil and gas reservoir, characterized by, The method comprises the following steps: Collecting core CT images, three-dimensional seismic data and production dynamic data and performing cross-scale preprocessing; Building a fractional-order non-local seepage field model to represent flow characteristics from nanometers to kilometers; Predicting seepage field parameters through a Lie group symmetry constrained neural network; Establishing a cross-scale coupling model of quantum adsorption effect and macroscopic seepage law; Dynamically updating model parameters based on real-time monitoring data; Performing multi-objective collaborative optimization to generate a sweet spot three-dimensional distribution and a development plan; The neural network comprises: The input layer receives the pore pressure , permeability tensor and viscosity ; At least 3 SL(2, R) group equivariant convolution layers, each layer having a channel number greater than or equal to 256; Output layer prediction seepage velocity field ; The loss function of the neural network is: wherein , the integral domain is consistent with the solution domain of the seepage field model; The cross-scale coupling model satisfies: where, , is the organic surface potential energy, potential well depth , is the equivalent inertia coefficient, determined from the permeability tensor .

2. The method of identifying sweet spots in a shale oil and gas reservoir of claim 1, wherein, The preprocessing comprises: Fractal dimension of pore network computed from core ct images The following linear relationship is obtained by fitting: wherein is a scan scale, is a gray scale threshold, is a number of blocks at the scan scale and the gray scale threshold, is a fitting constant term; Performing anisotropic diffusion filtering on seismic data, and the filtering kernel function is: where t diff denotes the diffusion time, is the noise sensitivity coefficient.

3. The method of identifying sweet spots in a shale oil and gas reservoir of claim 2, wherein, In the fractal dimension calculation: Scanning scale The value ranges from 5nm to 1μm, and the grayscale threshold is... CT value 200-400 HU; fractal dimension with spatial fractional derivative order satisfies .

4. The method of identifying sweet spots in a shale oil and gas reservoir of claim 1, wherein, The fractional-order non-local seepage field model satisfies: The time derivative term is taken in the Caputo definition, of order , characteristic time s; The spatial gradient term adopts Riesz fractional order operator, order , non-Darcy exponent .

5. The method of identifying sweet spots in a shale oil and gas reservoir of claim 1, wherein, The dynamic updating adopts: Collective Kalman filter algorithm ; wherein is a state vector of a neural network first layer; is a coupled residual vector; is an independent Gaussian noise vector; is a Jacobian sensitivity matrix; Process noise Diagonal element content quantum density gradient term ; Observation operator Correlation length , is the fractal dimension.

6. The method of identifying sweet spots in a shale oil and gas reservoir of claim 1, wherein, The multi-objective collaborative optimization comprises: Multi-objective functions: wherein is a discount factor; is an energy consumption penalty coefficient; is a flow rate fluctuation suppression coefficient.

7. The method of identifying sweet spots in a shale oil and gas reservoir of claim 6, wherein, The multi-objective collaborative optimization further comprises a control law, which is designed as: wherein denotes a Hadamard product, is a proportional gain matrix, is an integral gain matrix.

Images