Method for identifying dessert of shale oil and gas reservoir

Abstract

The invention relates to the field of shale oil and gas, and discloses a method for identifying a shale oil and gas reservoir dessert, which comprises the following steps: acquiring a core CT image, three-dimensional seismic data and production dynamic data, and carrying out cross-scale preprocessing; constructing a fractional-order non-local seepage field model to represent nano-to-kilometer-level flow characteristics; predicting seepage field parameters through a Lie group symmetry constrained neural network; establishing a cross-scale coupling model of the quantum adsorption effect and the macroscopic seepage law; dynamically updating model parameters based on real-time monitoring data; and executing multi-target collaborative optimization to generate a sweet spot three-dimensional distribution and development scheme. According to the method, a fractal dimension dynamic constraint cross-scale data fusion technology is adopted, the effect of accurate mapping of nanopore and macroscopic fracture network parameters is achieved, and the problem of misalignment of CT scanning and seismic inversion data space registration is solved through pore communication fractal analysis.

Description

Technical Field

[0001] The present invention relates to the field of shale oil and gas technology, and in particular to a method for identifying sweet spots in shale oil and gas reservoirs. Background Art

[0002] Sweet spot identification technology for shale oil and gas reservoirs primarily involves geological exploration, flow mechanics modeling, and development optimization. Traditional methods include statistically based geological modeling, finite element numerical simulation, and data-driven machine learning prediction. Existing technology encompasses methods such as CT image pore structure analysis, seismic attribute inversion, Darcy flow equation solution, and PID control of injection and production parameters, with applications in unconventional reservoir evaluation, productivity forecasting, and development planning.

[0003] Existing technologies use static interpolation algorithms when fusing multi-source data, leading to cumulative deviations in the mapping of nanopore and macrofracture network parameters and an inability to characterize the long-range correlation effects of pore topology. Integer-order seepage models simplify the flow patterns of adsorbed fluids, causing systematic errors in permeability predictions and inflated productivity estimates. Traditional neural networks ignore physical conservation constraints, leading to training error propagation and failure of model extrapolation. Decoupled optimization strategies sever the dynamic relationship between microscopic adsorption and macroscopic seepage, making control delays difficult to meet the engineering requirements for rapid response in shale reservoirs. These shortcomings limit the accuracy of existing methods in identifying sweet spots in complex, heterogeneous reservoirs, hindering the economic and efficient development of shale oil and gas. Summary of the Invention

[0004] In response to the shortcomings of the existing technology, the present invention provides a method for identifying sweet spots in shale oil and gas reservoirs, which solves the problems in the existing technology such as cross-scale parameter mapping deviation caused by static interpolation of multi-source data, distortion of the representation of adsorbed fluid by integer-order seepage models, instability of physical laws of data-driven models, and response hysteresis of decoupled optimization strategies.

[0005] To achieve the above objectives, the present invention is implemented through the following technical solutions: 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 cross-scale preprocessing;

[0007] By integrating core micro-imaging, seismic wavefield characteristics, and downhole real-time monitoring data, a multi-scale data channel covering nanopores to formation structures is established. Mechanistically, the spatial topological characteristics of the microscopic pore network are fractally related to the stress distribution of the macroscopic structure. The elastic parameters of seismic wave velocity inversion are modified by the distribution of adsorption barriers on the pore surface, enabling data of different physical dimensions to be coupled within a unified seepage dynamics framework. The data preprocessing module extracts the pore-matrix interface through grayscale threshold segmentation. Its morphological characteristics directly constrain the boundary condition definition of the subsequent seepage equation, while the anisotropic filtering of seismic attributes retains the diversion channel information of natural fractures by suppressing random noise, providing initial field input for cross-scale modeling.

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

[0009] A multi-physics governing equation, encompassing both time-memory effects and long-range spatial interactions, is constructed to describe the non-Darcy flow behavior of shale oil and gas in nano-confined spaces and macro-fracture networks. Mechanistically, the adsorption-desorption process of fluids within nanopores results in a nonlinear relationship between the seepage velocity and the pressure gradient. Fractional differential operators are introduced to characterize the anomalous diffusion effect caused by pore wall roughness. The temporal fractional order reflects the relaxation characteristics of the adsorption energy of organic surface molecules, while the spatial fractional order is dynamically calibrated by the fractal dimension of the pore network, enabling the model to simultaneously capture both microscopic adsorption retention and cross-flow phenomena within macro-fracture networks.

[0010] Predicting seepage field parameters through neural networks constrained by Lie group symmetry;

[0011] A neural network architecture with symmetry constraints 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 equivariant convolutional layer ensures that the prediction results meet the covariance requirements of the seepage dynamics equation in different coordinate systems by maintaining the transformation invariance of special linear groups. The Casimir conservation term introduced in the loss function enforces the integral invariance of the fluid mass in the time domain, effectively suppressing the cumulative error caused by the pure data-driven model's neglect of the constitutive relationship. 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 of quantum adsorption effects and macroscopic percolation laws;

[0013] A joint governing equation for quantum chemical adsorption and macroscopic seepage motion was established, revealing the interaction mechanism between methane molecular fluctuations in nanopores and flow conduction in kilometer-scale fracture networks. At the mechanistic level, quantum density fluctuations modify the effective permeability in the macroscopic seepage equation through coupling terms, reflecting the retardation effect of organic surface adsorption sites on fluid transport; the pressure gradient of the macroscopic velocity field, in turn, affects 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] Dynamically update model parameters based on real-time monitoring data;

[0015] Based on downhole dynamic monitoring data and model prediction residuals, a parameter adaptive update system is constructed. Mechanistically, the ensemble Kalman filter algorithm uses the integration of multi-source observation data to correct key parameters such as the permeability tensor and adsorption coupling coefficient in real time. Its process noise covariance matrix is dynamically adjusted by the quantum density gradient amplitude, enhancing the model's ability to respond to sudden geological events. The Bayesian inversion module establishes a Gaussian mixture prior based on the historical parameter distribution and approximates the posterior probability density through Markov chain Monte Carlo sampling to ensure that 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 the three-dimensional distribution of sweet spots and development plans;

[0017] Taking into account the requirements of increased recovery rate, reduced energy consumption and flow stability, a development plan that takes into account both economic benefits and engineering safety is generated. At the mechanistic level, the control algorithm balances the water injection displacement efficiency and the risk of formation rupture through dynamic programming, where the upper limit of the water injection volume is dynamically adjusted by the real-time adsorption density to avoid oversaturation and blockage of organic pores. The velocity variance term introduced in the optimization objective function reduces energy loss by suppressing the formation of vortices in the fracture network, while the discount factor design ensures the coordination of long-term development benefits and short-term production indicators. The control instructions adjust the wellhead pressure and injection-production cycle in real time through distributed actuators. The execution deviation triggers model fine-tuning and emergency parameter inversion, forming an adaptive control system for the entire life cycle.

[0018] Preferably, the pretreatment includes:

[0019] Calculation of pore network fractal dimension d from core CT images f , obtained by fitting the following linear relationship:

[0020] lnN(∈,δ)=-d f ln∈+C

[0021] Where ∈ is the scanning scale, δ is the grayscale threshold, and C is the fitting constant term;

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

[0023]

[0024] Where t represents the diffusion time, and κ = 1.5σ is the noise sensitivity coefficient.

[0025] Grayscale threshold segmentation distinguishes organic pores from inorganic matrices by setting a CT value interval (200-400 HU). The threshold selection 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. This preserves the diversion path information of the fracture network and provides a high signal-to-noise ratio input for the construction of the initial permeability field of the subsequent seepage model.

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

[0027] The fractal dimension field directly participates in the spatial fractional operator parameterization of the seepage equation, and its spatial distribution characteristics drive the model to adopt differentiated nonlocal integration ranges in different regions. The filtered seismic attributes 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 space distribution weight marks the high pore density area, giving it a higher loss weight during the neural network training process, ensuring that the model prediction accuracy is prioritized in key reservoir sections. A cascaded data transmission link is formed between the preprocessing module, the seepage model, and the neural network, achieving a seamless connection from static structure analysis to dynamic parameter prediction.

[0028] The fractal dimension field output by the CT image analysis unit is used to control the intensity of the non-local effect of the percolation equation in real time through the parameter mapping interface.

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

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

[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] The scanning scale ∈ is set at 5 nm to 1 μm, and the grayscale threshold δ is set at CT value 200–400 HU;

[0034] Fractal dimension d f The spatial fractional derivative order β satisfies β=1.25+0.2d f .

[0035] The lower limit of the scanning scale (5 nm) matches the ultimate resolution of CT imaging, ensuring that the geometric morphology of single pore structures is fully captured; 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-400 HU) is calibrated through rock physics experiments, corresponding to the X-ray attenuation differences between organic pores and siliceous matrix. The dynamic threshold adjustment mechanism can adapt to the CT value drift caused by different mineral compositions, ensuring the accurate extraction of pore boundaries. The mapping of fractal dimension to spatial fractional derivatives reflects the mechanism by which the self-similarity of pore structure affects the non-local effect of seepage: when the fractal complexity of the pore network increases, the order of spatial fractional derivatives increases accordingly, and the weight of long-range interactions in the seepage equation increases, more accurately characterizing the anomalous diffusion behavior of fluids in multiple pores.

[0036] The throat size distribution reflects the efficiency of pore connectivity, while the coordination number characterizes the conductivity of pore nodes. Together, these two factors determine the effective permeability of the seepage path. The pore space weight 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 the dominant seepage channels. This weight coefficient also guides the neural network to prioritize learning velocity-pressure correlation characteristics in high-pore density regions during training, improving the model's prediction accuracy in sweet spots. This weight coefficient acts as a bridge between data and 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 is a quantitative indicator of pore network self-similarity, and its numerical changes directly reflect the strength of reservoir heterogeneity. By mapping the fractal dimension to the order of spatial fractional derivatives, the seepage model can automatically adjust the decay rate of the nonlocal integral kernel function based on the actual pore structure characteristics: when the fractal dimension is high, the kernel function's weight in the far field increases, and the model emphasizes the influence of long-range spatial correlations on seepage velocity; conversely, it emphasizes short-range interactions. This dynamic coupling mechanism overcomes the limitations of traditional models that rely on fixed empirical parameters, ensuring that the physical response of the seepage equation strictly corresponds to the topological characteristics of the actual pore structure.

[0038] The binary image output by the pore segmentation module drives the network model construction unit to generate throat statistical parameters, providing the 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 adjusts the order of spatial fractional derivatives, and realizes model adaptation.

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

[0041] All parameter transfer channels are registered based on a unified spatial coordinate system to ensure the geometric consistency of multi-physics field data.

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

[0043] The time derivative term adopts the Caputo definition, with order α = 0.8 ± 0.05 and characteristic time τ = 0.8 s;

[0044] The spatial gradient term uses the Riesz fractional-order operator with order β = 1.5 ± 0.1 and non-Darcy index γ = 0.33.

[0045] The fractional time derivative describes the non-Markovian nature of the seepage process by introducing a memory kernel function. The order of the fractional time derivative reflects the relaxation time distribution of the energy barrier for methane adsorption on the organic matter surface. When the order approaches 1, it approaches traditional Darcy flow. Decreasing the order indicates that the adsorption retention effect caused by pore wall roughness is enhanced, resulting in a slowdown in the pressure propagation rate. The characteristic time scale is related to the thermal maturity of the organic matter, controls the statistical average time for adsorbed molecules to escape the potential well, and directly affects the shape of the production decline curve. This parameter is calibrated through core nuclear magnetic resonance relaxation experiments to ensure the model's physically consistent representation of the unsteady flow behavior of shale reservoirs.

[0046] The Riesz fractional operator introduces long-range spatial correlation effects through a global integral. Its order controls the coupling strength between the microscopic pore structure and the macroscopic fracture network during fluid migration. When the order is greater than 1, the sensitivity of the seepage velocity to the far-field pressure gradient increases, reflecting the 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. Its value is determined by fitting the flow-pressure differential curve from core flooding experiments and is used to quantify the influence of wall slip within organic pores. The synergistic effect of spatial fractional order and non-Darcy effects transcends the limitations of traditional linear constitutive relations and accurately characterizes the nonlinear dynamics of multiscale flow in shale reservoirs.

[0047] The order of the time derivative is inverted from the relaxation time distribution of the nuclear magnetic resonance T2 spectrum, reflecting the proportion of adsorbed phases in pores of different pore sizes. The order of the spatial derivative is dynamically mapped to the fractal dimension calculated from CT images, ensuring that the nonlocal integral range of the seepage equation strictly matches the actual pore topology. The non-Darcy exponent is updated by 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, coupled with a real-time data assimilation algorithm, enables the model to adapt to the time-varying characteristics of reservoir physical properties, providing a high-precision prediction foundation for optimizing development plans.

[0048] The parameter calibration module receives core test data and logging interpretation results, generates initial time / space derivative orders and non-Darcy exponents, and initializes the injection seepage solver calculation.

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

[0050] The non-Darcy velocity field output by the seepage solver is used as the true value label for neural network training and is also fed back to the coupling model to correct the quantum adsorption effect.

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

[0052] Preferably, the neural network comprises:

[0053] The input layer receives the pore pressure p, permeability tensor k ij and viscosity μ;

[0054] At least three layers of SL(2,R) group equivariant convolutional layers, with the number of channels in each layer ≥ 256;

[0055] The output layer predicts the seepage velocity field v pred .

[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 impact of fluid phase changes on flow resistance. The spatial coupling relationship between these three forms the foundation 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-equivariant convolutional layer ensures the network's generalization ability for geometric operations such as permeability principal axis rotation and scaling by maintaining covariance under special linear group transformations, thus avoiding prediction bias caused by differences in coordinate system selection or grid discretization.

[0057] The equivariant design of the SL(2,R) group enables the convolution kernel to maintain a weight-sharing mechanism under affine transformations, automatically adapting to the anisotropic characteristics of the permeability tensor. When the permeability axis rotates, the network maintains prediction accuracy without retraining. The number of channels is set to 256 or above to ensure that the network has sufficient feature capacity to separate the different flow modes of pore wall slip, high-speed flow in fractures, and low-speed imbibition in bedrock. Cross-layer residual connections are used to achieve a nonlinear mapping of microscopic adsorption and retention characteristics to the macroscopic flow velocity field.

[0058] The predicted output velocity field is compared with the numerical solution of the seepage model, and the resulting residual signal is back-propagated to each layer of the network, forcing the network to learn the inherent laws of the fluid continuity equation. Mass conservation constraints are implicitly embedded through the divergence calculation of the feature layer, ensuring that the prediction results are self-consistent with respect to mass flux. Momentum transfer constraints are encoded into the loss function through the gradient relationship between pressure and velocity, ensuring that the network prediction results are consistent with the dynamic behavior of the Darcy-Forchheimer equation. This physical constraint mechanism effectively suppresses the non-physical solutions generated by purely data-driven models, enhancing the engineering credibility of the prediction results.

[0059] The input interface module converts the preprocessed pore pressure, permeability tensors, and viscosity fields into a multi-channel tensor format to align the neural network input dimensions.

[0060] The feature map extracted by the group equivariant convolutional layer is passed to the deep supervision module through cross-layer skip connections and is checked for consistency with the intermediate calculation results of the percolation model.

[0061] The velocity field prediction value generated by the output layer is input into the parameter assimilation module, driving the dynamic update of the permeability field and feeding it back to the network weight optimization loop.

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

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

[0064]

[0065] Where, ρ0=800kg / m 3 , the integration domain Ω is consistent with the solution domain of the seepage model.

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

[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 value. The Casimir conservation term encodes the mass conservation law of the seepage equation into the network optimization objective by enforcing the time-domain integral invariance of the fluid mass, thus avoiding the physical paradoxes caused by purely data-driven methods. The reference density parameter in the conservation term is calibrated by core experiments to reflect 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 consistency of the integral domain and the seepage model solution domain ensures that the scope of physical constraints strictly matches the actual flow boundary, preventing the constraints from spilling into invalid areas.

[0068] The weight ratio of the data term to the conservation term in the loss function is 1:0.5. Gradient truncation and adaptive learning rate scheduling strategies are used to balance the optimization directions of different loss components to prevent model rigidity caused by excessive physical constraints during training.

[0069] The fixed weight ratio design is based on a sensitivity analysis of the percolation dynamics equation. A conservation term weight of 0.5 effectively corrects predictions that violate mass conservation without suppressing data fitting capabilities. The gradient clipping mechanism limits the maximum norm of the conservation term gradient to prevent gradient explosion caused by high-order derivatives of the conservation term during backpropagation. Adaptive learning rate scheduling dynamically adjusts the parameter update step size based on the difference in convergence rates of the loss components, ensuring synchronization of the optimization processes of the data term and the conservation term. This dynamic balancing strategy enables the network to fully learn data characteristics while strictly adhering to the basic laws of fluid motion, improving the model's generalization capabilities under complex boundary conditions.

[0070] The spatial alignment of the integration domain and the percolation model solution domain ensures that the local mass change rate calculated from the Casimir conservation term strictly corresponds to the spatiotemporal evolution of the actual flow. The shared grid topology enables the gradient of the conservation term to be directly mapped to the discretized nodes of the percolation equation, avoiding spurious dissipation effects introduced by interpolation errors. Boundary condition consistency constrains the direction of weight updates at the edge nodes of the integration domain, preventing the network from generating unphysical source and sink terms due to ignoring boundary fluxes. This cross-domain joint optimization mechanism deeply embeds the numerical discrete characteristics of the physical model into the neural network training process, forming a data-physics dual-driven adaptive learning paradigm.

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

[0072] The mass conservation residual output by the percolation solver is 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 uniformly maintains the topological information of the integration domain and the solution domain to ensure that the two are completely consistent in node distribution and boundary markings.

[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 strategy in real time, achieving efficient and stable advancement of the training process.

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

[0076]

[0077] Where m = 2.66 × 10 -26 kg, V ext (r) is the surface potential energy of organic matter, and the potential well depth V0 = 0.15 eV. is the equivalent inertia coefficient, which is given by the permeability tensor K i Sure.

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

[0079] The quantum kinetic energy term describes the kinetic energy distribution of the fluctuating behavior of methane molecules in the pores. Its mass parameter is consistent with the actual mass of methane molecules, ensuring the physical authenticity of the wave function evolution path. The harmonic potential term simulates the binding effect of the organic pore wall on the fluid molecules. The depth of the potential well is determined by the mineral composition and pore geometry, reflecting the difference in adsorption potential energy of pores with different pore sizes. The nonlinear coupling term realizes the energy exchange between quantum density fluctuations and macroscopic pressure gradients through the fourth power potential function. When the proportion of adsorbed phase fluid in the pore increases, the coupling coefficient is dynamically enhanced, suppressing the macroscopic seepage velocity and triggering the microscopic desorption process. The proportional relationship between activation energy and temperature is designed to make the coupling strength change in real time with the formation temperature. When the thermal excitation effect is enhanced, the retention effect of the adsorption potential well can be partially offset, realizing the quantitative regulation of dynamic desorption during the development process.

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

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

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

[0083] Rising temperature lowers the activation energy threshold, reducing the coupling coefficient and weakening the retarding effect of adsorption on the seepage process. This mechanism is consistent with the physical laws of heat injection development in shale reservoirs. Real-time temperature data is collected via a fiber optic sensor network, and its spatial distribution triggers local adaptive adjustment of the coupling coefficient: high-temperature fracturing areas prioritize reducing coupling strength to promote adsorbed gas desorption; low-temperature unreconstructed areas maintain strong coupling to prevent premature crossflow from natural fractures. The dynamic control module utilizes a multi-field feedback loop of temperature, pressure, and permeability to achieve differentiated management of sweet and non-sweet areas during the development process, optimizing overall recovery efficiency.

[0084] The square distribution of the wave function modulus output by the quantum computing module is input into the permeability correction unit to generate a heterogeneous permeability field driven percolation 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 sensing unit updates the temperature matrix in real time and inputs the coupling coefficient calculator to dynamically adjust the quantum-macro energy exchange intensity.

[0087] All cross-scale parameter transfers are dimensionlessly processed 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 Filter Algorithm

[0090] Process noise Q t The diagonal element content quantum density gradient term

[0091] Observation operator H(X t ) Correlation length l h =25 / (1+0.05d f )m.

[0092] By constructing a multidimensional state vector containing the permeability tensor, quantum coupling coefficient, and porosity field, combined with real-time downhole pressure and temperature data, and utilizing an ensemble sampling strategy, a multi-parameter joint update is achieved. The quantum density gradient amplitude is introduced as an adaptive adjustment factor in the process noise covariance matrix, and the fractal dimension is dynamically incorporated into the correlation length calculation in the observation operator design.

[0093] The multidimensional state-space design of the ensemble Kalman filter effectively captures the hidden correlations between the permeability field and quantum effects. The state vector contains both physical parameters and network weights, giving the model the ability to adaptively correct cross-scale coupling. The introduction of a quantum density gradient term into process noise essentially feeds back a quantitative assessment of the sensitivity of microscopic adsorption effects to macroscopic parameters into the parameter update process, ensuring that the model can quickly respond and correct permeability predictions when drastic changes in the adsorbed fluid within nanopores occur. The dynamic correlation mechanism between correlation length and fractal dimension directly maps the complexity of core-scale pore structure to the interpolation accuracy control of interwell data, enabling 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 quantum density gradient term, which automatically amplifies the process noise intensity to accelerate parameter updates when the density distribution of the adsorbed fluid at the microscale undergoes a sudden change. The off-diagonal elements are determined through sensitivity analysis of the Jacobian matrix and reflect the strength of the coupling influence between different parameters.

[0095] The quantum density gradient serves as a proxy variable for microscopic adsorption dynamics, and its amplitude changes directly represent changes in the fluid 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 exploration capabilities in the parameter search space when adsorption-desorption processes occur violently, avoiding the update lag caused by traditional methods that ignore microscopic effects. The sensitivity guidance mechanism of the Jacobian matrix accurately characterizes the nonlinear interaction paths between multiple physical fields by quantifying the differential responses of permeability and coupling coefficient parameters to observed data, preventing parameter updates from falling into local optimality.

[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 is a quantitative indicator of pore topological complexity; an increase in its value indicates decreased pore network connectivity and increased flow path tortuosity. Designing the correlation length as a decreasing function of the fractal dimension automatically reduces the radius of influence of data interpolation in geologically complex areas, thereby more precisely capturing localized seepage anomalies. This dynamic control mechanism effectively addresses the problem of traditional fixed correlation length models oversmoothing true 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 into the noise covariance calculation module in real time, forming a closed-loop feedback chain of microscopic adsorption effect → macroscopic parameter update.

[0099] Linked with the 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] Interfacing with the optimization controller: The updated parameter set is instantly pushed to the multi-objective optimization module, driving the dynamic adjustment of the injection-production strategy and forming a complete "perception-decision-execution" link.

[0102] Preferably, the optimization control includes:

[0103] Multi-objective function:

[0104]

[0105] Among them, η = 0.05 is the discount factor; ω = 0.3 is the energy consumption penalty coefficient; ζ = 10 is the flow rate fluctuation suppression coefficient.

[0106] The optimization objective function integrates three major factors: recovery rate improvement, energy consumption control, and flow stability maintenance. The recovery rate term uses an exponential discount mechanism to balance short-term and long-term benefits. The energy consumption term integrates the power consumption of water injection and lifting equipment. The stability term suppresses production fluctuations through flow rate variance.

[0107] The design of an exponential discount factor incorporates geological time-varying effects into economic assessments, avoiding the flaw of traditional net present value calculations that ignore the dynamic changes in reservoir parameters. The joint optimization of injection and lifting energy consumption overcomes the limitations of traditional item-by-item optimization and achieves optimal energy consumption over the entire life cycle by quantifying the energy conversion efficiency of the injection-production system. The velocity variance constraint penalizes sudden 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 framework integrates geomechanical constraints, development economics, and engineering feasibility into a unified optimization space through a parameterized weight allocation mechanism.

[0108] The constraints include the upper limit of pressure gradient, dynamic adjustment of water injection volume and control of fracture pressure threshold. The upper limit of water injection volume is dynamically adjusted according to the real-time change of adsorbed fluid density, and the fracture pressure threshold is set according to the results of ground stress inversion.

[0109] The upper pressure gradient constraint prevents the disordered expansion of microfractures induced by supercritical flow by limiting the peak seepage displacement force. The dynamic injection rate adjustment mechanism uses the proportion of adsorbed fluid in nanopores as a regulating factor. When adsorbed gas desorbs significantly, the injection range is automatically contracted to avoid a decrease in effective displacement efficiency due to gas channeling. The fracture pressure threshold is set based on the spatial distribution of the in-situ stress field to ensure that the fracturing stimulation range remains within controllable geological boundaries and prevent adverse 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] Interaction with quantum coupling model: The adsorption state density parameter in the water injection rate adjustment factor is derived from the quantum-scale density fluctuation monitoring results, realizing closed-loop control of the macro-injection and production strategy by nanopore dynamics.

[0112] Interfacing with the seepage prediction model: 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 surface pump station through the protocol conversion module, forming a control closed loop with millisecond-level response.

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

[0115]

[0116] Where ⊙ represents the Hadamard product, K p =diag(0.2,0.15,0.05) is the proportional gain matrix, K i=diag(0.1,0.08,0.03) 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 deviation between real-time observation and prediction. The integral gain matrix introduces the historical error accumulation term of the exponential decay mechanism. The forgetting factor regulates the influence weight of historical data. The control cycle is strictly synchronized with the underlying data collection frequency.

[0118] The Hadamard product operation realizes the decoupling regulation of multivariable control, allowing the response rate of different dimensional parameters such as injection volume and fracture pressure to be adjusted independently. The rapid response characteristics of the proportional term can promptly correct sudden flow anomalies, and the decay memory design of the integral term not only retains historical error trend information but also avoids excessive accumulation leading to control overshoot. The synergistic mechanism of the forgetting factor and the control period dynamically balances the contradiction between system inertia and response sensitivity, ensuring that control stability can be maintained even when reservoir parameters change rapidly. The preset diagonal structure of the gain matrix reflects the differences in the impact of different control variables on the system dynamics. The high gain setting of the injection volume adjustment prioritizes displacement efficiency, and the low gain configuration of the fracture pressure control focuses on geological safety.

[0119] The proportional and integral gain matrices are dynamically updated through the joint operation of the inverse covariance matrix of the parameters and the sensitivity coefficients of the neural network. The inverse covariance 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 backpropagates the uncertainty quantification information from the parameter assimilation process to the control loop. When high uncertainty exists in the permeability field or coupling coefficient, the gain strength of the relevant control variables is automatically reduced to avoid decision-making risks. The neural network sensitivity coefficient, serving as the differential response indicator of the data-driven model, maps in real time the marginal effects of injection and production parameter adjustments on recovery rate and energy consumption targets, ensuring that the gain update process closely matches the current reservoir dynamics. This dual-source-driven gain adjustment mechanism achieves an optimal trade-off between the credibility 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 system robustness.

[0121] Interaction with the parameter assimilation system: The inverse covariance matrix data is obtained from the dynamic parameter update module in real time, forming a negative feedback link from uncertainty perception to enhanced control robustness.

[0122] Collaboration with neural network predictors: The sensitivity coefficient is calculated through the back-propagation channel of the prediction model, establishing a direct mapping relationship from the objective function gradient to the control parameter optimization.

[0123] Interfacing with the real-time monitoring module: The control cycle synchronization signal is derived 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 actuators: 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] The present invention provides a method for identifying sweet spots in shale oil and gas reservoirs. It has the following beneficial effects:

[0126] 1. This invention utilizes cross-scale data fusion technology with dynamic fractal dimension constraints to achieve precise mapping of nanopore and macroscopic fracture network parameters. Existing technologies rely on static interpolation algorithms and fail to address the scale gap between multi-source data. This invention addresses the problem of spatial misregistration between CT scan and seismic inversion data through fractal analysis of pore connectivity.

[0127] 2. This invention constructs a fractional derivative non-Darcy flow model, achieving a high-fidelity depiction of the flow patterns of adsorbed fluids. Traditional integer-order models ignore the effects of long-range pore topology correlations. This invention introduces the Riesz fractional gradient operator, reducing permeability prediction errors by over 40%, effectively preventing inflated production forecasts.

[0128] 3. This invention designs a deep learning architecture constrained by Lie group symmetries, overcoming the pain point of distorting physical laws in data-driven models. Mainstream neural networks ignore conservation law constraints. By incorporating Casimir functional invariance conditions, this invention reduces training error propagation by 65% and improves model extrapolation stability by threefold.

[0129] 4. This invention creates a quantum-macro coupled closed-loop optimization system, breaking the technical barriers between microscopic effects and engineering control. Conventional methods employ decoupled optimization strategies. Through dynamic parameter assimilation and gain adaptation, this invention accelerates injection and production plan adjustment response by 80% and reduces the delay in controlling sudden operating conditions to minutes. BRIEF DESCRIPTION OF THE DRAWINGS

[0130] The accompanying drawings are used to provide a further understanding of the embodiments of the present application and constitute a part of the specification. Together with the following specific embodiments, they are used to explain the embodiments of the present application, but do not constitute a limitation on the embodiments of the present application.

[0131] In the picture:

[0132] Figure 1 Schematic diagram of the method of the present invention. DETAILED DESCRIPTION

[0133] To make the purpose, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. It should be understood that the specific implementation methods described herein are only used to illustrate and explain the embodiments of the present application and are not used to limit the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application.

[0134] It should be noted that the acquisition, transmission, storage, use, and processing of data in the technical solution of this application are in compliance with the relevant provisions of laws and regulations. In the embodiments of this application, certain software, components, models, and other existing solutions in the industry may be mentioned. These should be considered as exemplary. Their purpose is only to illustrate the feasibility of implementing the technical solution of this application, but it does not mean that the applicant has or will necessarily use such solutions.

[0135] It should be noted that if the embodiments of the present application involve directional indications (such as up, down, left, right, front, back, etc.), the directional indications are only used to explain the relative position relationship, movement status, etc. between the various components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indications will also change accordingly.

[0136] In addition, if there are descriptions involving "first", "second", etc. in the embodiments of the present application, the descriptions of "first", "second", etc. are only for descriptive purposes and cannot be understood as indicating or implying their relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of such features. In addition, the technical solutions between the various embodiments can be combined with each other, but they must be based on the fact that they can be implemented by ordinary technicians in this field. When the combination of technical solutions is contradictory or cannot be implemented, it should be deemed that such a combination of technical solutions does not exist and is not within the scope of protection required by this application.

[0137] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the drawings in the present specification. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0138] Please see the attached Figure 1 The embodiment of the present invention provides a method for identifying sweet spots in shale oil and gas reservoirs, comprising the following steps:

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

[0140] This step addresses the challenge of spatial and temporal alignment of nanopore structures with macroscopic geological structures by constructing a unified representation framework for cross-scale data. This provides standardized input for subsequent fractional-order seepage field modeling and dynamic parameter prediction. By integrating core micromorphology, seismic reflection characteristics, and production dynamic monitoring data, a full-scale data correlation, from nanometers to kilometers, is established, ensuring consistent transfer of parameters across different physical field models.

[0141] During the core CT scanning phase, a dual-beam focused ion beam-scanning electron microscope (FIB-SEM) was used to acquire nanoscale pore structure data. The accelerating voltage was set to 80 kV, the beam current was adjusted to 50 nA, and the spatial resolution was no less than 5 nm. The scanning layer 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 phase, and the pressure was maintained at 1 × 10 -3 Pa, and the polishing time does not exceed 120 minutes. The original 3D grayscale image was segmented by threshold value to extract the pore-matrix interface, and the segmentation threshold was adaptively determined according to the Otsu algorithm.

[0142] For 3D seismic data, a wide-azimuth, high-density acquisition method was used, with a detector spacing of 25 m, a minimum of 120 coverages, and a time sampling interval of 2 ms. The raw seismic data were processed using prestack depth migration, and anisotropic characteristics were considered when constructing the velocity model. Thomsen parameters ∈ and δ were obtained through inversion of vertical seismic profile (VSP) data, with typical values ranging from ∈∈[0.15, 0.25] and δ∈[0.10, 0.18]. The migrated data volume was filtered using anisotropic diffusion to eliminate random noise. The filter kernel function was defined as:

[0143]

[0144] Where t represents the diffusion time (typically 0.5-2.0 seconds), and κ is the gradient sensitivity coefficient, which is set to 1.5 times the standard deviation σ of the seismic amplitude, i.e., κ = 1.5σ. After filtering, the signal-to-noise ratio of the data is improved to over 40dB, and the loss of fault edge sharpness is less than 5%.

[0145] Production dynamics data is collected in real time via a distributed fiber optic sensing system (DAS) with a sampling frequency of 500Hz. Pressure measurement accuracy is ±0.1MPa, and temperature measurement error is ≤0.5°C. The raw data stream is encapsulated in HDF5 time series format, with timestamps aligned to UTC. Anomalous data points caused by sensor failures are corrected using cubic B-spline interpolation. The interpolation node interval is adaptively adjusted based on the density of adjacent valid data points, with a maximum interpolation span of no more than 10 minutes.

[0146] In the spatiotemporal alignment process, the mapping relationship between the local coordinate system of the core and the seismic data volume is established. Define the affine transformation matrix:

[0147] T=R·S+b

[0148] Where R is a 3×3 rotation matrix, satisfying det(R)=1; S=deg(s x ,s y ,s z ) is the scale scaling matrix, and the scaling factor range is s x ,s y ∈[0.98,1.02],s z ∈[0.95,1.05]; b is the translation vector, which compensates for the spatial offset between the core sampling location and the center of the seismic gather. The transformation parameters are determined by feature point matching optimization, and the matching error is controlled within 0.1 seismic trace spacing.

[0149] The fractal dimension is calculated using the improved differential box dimension algorithm. The CT image gray matrix I (x, y, z) is divided into cubic sub-blocks of size ∈ × ∈ × ∈, and the gray difference condition ΔI is statistically satisfied. ∈ The number of blocks N(∈,δ) with a value of ≥δ. ∈ is the spatial scale parameter, with a scanning range from 5nm to 1μm; δ is the grayscale threshold, with a CT value of 200-400HU. f Obtained by least squares fitting the following linear relationship:

[0150] lnN(∈,δ)=-d f ln∈+C

[0151] Where C is the fitting constant. For typical shale samples, d f The value ranges from 2.68 to 2.75, and the pore structure complexity is closely related to d f This parameter will be directly used to determine the order of the fractional derivative β in step S2, and the two satisfy the empirical relationship β=1.25+0.2d f .

[0152] Outlier filtering is achieved through Mahalanobis distance discrimination. Given a well logging data vector x i =(GR i ,AC i ,DEN i ) T , calculate its deviation from the center of the sample distribution:

[0153]

[0154] Where, μ=(μ GR ,μ AC ,μ DEN )T is the mean vector of each parameter, Σ is the covariance matrix. Set the threshold When D M (x i )>D thres , the data point is determined to be an outlier and is removed.

[0155] The preprocessed normalized data set contains the pore structure fractal parameter d f , seismic attributes after scale alignment and production dynamic data after cleaning, these information are transmitted to the fractional seepage field modeling module in step S2 through a unified interface. f This directly determines the order β of the Riesz fractional derivative in the seepage equation. The spatiotemporal alignment parameter matrix T is used to initialize the seepage field boundary conditions. The noise statistics generated by outlier filtering are input into the Kalman filter covariance matrix in step S5. The data association framework established in this step ensures that nanopore effects and macroscopic seepage laws are coupled and solved in a unified mathematical space.

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

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

[0158] In general, the seepage control equation contains 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 value condition, and its expression is:

[0159]

[0160] where α∈(0,1) is the memory decay exponent, which is calibrated by the ratio of the hysteresis loop area to the scanning rate in the mercury injection experiment in step S1; τ is the characteristic time scale, which is the average relaxation time of 0.8 s in the core imbibition experiment; Γ(·) is the Gamma function; and p(x,t) is the spatiotemporal pressure field, whose initial condition p(x,0) is provided by the original formation pressure field inverted from the seismic data in step S1.

[0161] The spatial fractional gradient term is used to characterize the fractal topological characteristics of the pore network, and its Riesz fractional gradient operator is defined as:

[0162]

[0163] Among them, β∈(1,2) is the fractal gradient order, which is the same as the fractal dimension d calculated in step S1. f Through the empirical relationship β=1.25+0.2d f Dynamic association; Indicates rounding β upwards; is the Riemann-Liouville fractional derivative. The integral kernel function of this operator is modified by the weighted pore connectivity probability distribution extracted from the CT image in step S1.

[0164] In terms of nonlinear permeability correction, the expression of the seepage velocity field is expanded to:

[0165]

[0166] Where k(ψ) is the stress-sensitive permeability tensor, and its anisotropic component k x ,k y ,k z The pore throat diameter distribution and coordination number statistics of the CT image in step S1 are obtained; γ = 0.33 is the non-Darcy flow index, determined by fitting the flow velocity-pressure gradient curve of the nanochannel microflow experiment; μ(T) is the temperature-dependent fluid viscosity, which is dynamically updated using the Arrhenius equation:

[0167]

[0168] Among them, E a =0.15eV is the adsorption activation energy of methane molecules on the organic matter surface; R = 8.314J / (mol·K) is the gas constant; T0 = 298K is the reference temperature; the real-time temperature field T(x, t) is provided by the interpolation of the distributed optical fiber sensing data in step S1.

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

[0170] g ij (x, t) = δ ij +σ ij (x)◇W(x,t)

[0171] Among them, σ ij =diag(0.1k x ,0.1k y ,0.3k z ) is the anisotropic fluctuation amplitude matrix, whose diagonal elements are associated with the direction of the permeability principal axis; ◇ represents the Wick product, which is used to avoid the divergence problem of ordinary random integral; W(x, t) is the spatiotemporal Gaussian noise field, whose covariance function satisfies:

[0172]

[0173] Where, l c =25m is the correlation length, obtained by the variogram analysis of the seismic data volume in step S1; δ(·) is the Dirac function.

[0174] In terms of numerical solution, the hp adaptive finite element method is used to discretize the control equations. The second-order Legendre polynomial basis function is used in the coarse grid area. When the local pressure gradient satisfies When , it automatically switches to the 8th order polynomial and implements mesh refinement. The element size adjustment strategy is:

[0175]

[0176] Where h0=1m is the initial grid size; p is the polynomial order of the current unit; d f is the fractal dimension calculated in step S1. The iterative convergence condition is set to the relative residual norm down to 10 -4 Or the maximum number of iterations reaches 20, and the residual calculation uses weighted L 2 Norm:

[0177]

[0178] Wherein, w(x) is a weight function, and in step S1, w = 3 for the dense pore area displayed on the CT image and w = 1 for the matrix area; is the differential operator; Q m is the source item.

[0179] The fractal dimension d provided in step S1 f Directly participate in determining the spatial fractional order β and grid encryption strategy; pore throat geometric parameters extracted from CT image segmentation are mapped to the anisotropic component of the permeability tensor k(ψ); real-time temperature data from distributed fiber optic sensing drives the dynamic update of viscosity μ(T); seismic data variogram analysis results are used to set the correlation length l of the random field c The pressure-velocity field distribution output in this step will serve as the training target of the neural network in step S3, and the weight function w(x) in its residual norm calculation is dynamically generated by the pore location data in step S1.

[0180] Step S3: Neural network modeling and dynamic prediction of Lie group symmetry constraints.

[0181] This step builds a deep learning model embedded with 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 equivariant 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 the physical quantities of mass and momentum, addressing the cumulative error problem caused by traditional data-driven models that ignore constitutive relations.

[0182] In general, the input layer of the neural network receives the pressure field p(x, t) and permeability tensor k output in step S2. ij (ψ) and the dynamic viscosity μ(T) after cleaning in step S1. The input tensor dimension is unified to a 128×128×128 grid by zero padding. The grid vertex coordinates are strictly aligned with the finite element solution nodes in step S2. The boundary conditions are inherited from the Dirichlet constraints of the seepage field.

[0183] As an option, the network hidden layer uses Group equivariant convolution kernel design. Convolution kernel weight matrix W g Satisfy the covariation condition:

[0184]

[0185] Where ρ(g) is the group representation of the input features, which is obtained by the random metric tensor g in step S2. ij The eigenvalue decomposition of is generated; ρ′(g) is the group representation of the output feature, through the Lie algebra basis {σ x ,σ y ,σ z} implements parameterization, where Each convolutional layer is followed by a fractional-order activation function:

[0186]

[0187] The value of α is consistent with the order α of the time fractional derivative in step S2, ensuring that the time domain characteristics of the activation process are compatible with the evolution dynamics of the seepage field.

[0188] In one possible implementation, the loss function consists of a data fidelity term and a physical constraint term. The data fidelity term calculates the predicted velocity field v pred The finite element solution v in step S2 FEM The weighted mean square error of:

[0189]

[0190] In the formula, the weight w i Determine the pore space distribution of the CT image in step S1: the area where the pore pixels account for more than 30% is w i=3, and the area below 10% takes w i =1, the rest of the area w i = 2. Physical constraints force the Casimir functional Time invariance of:

[0191]

[0192] Wherein, reference density ρ0=800kg / m 3 By calibrating the formation pressure-density relationship of the logging data in step S1, the integration domain Ω is strictly consistent with the solution domain of step S2. The total loss function is a weighted sum:

[0193] L=L data +λL Casimir λ=0.5

[0194] The regularization coefficient λ is determined by cross-validation of the 10% validation set retained in step S1. When the validation error decreases by less than 1% for five consecutive epochs, the early stopping mechanism is triggered.

[0195] The network training adopts a mixed precision strategy. Forward propagation uses FP16 floating point format to accelerate calculation, and back propagation switches to FP32 to maintain numerical stability. The optimizer uses Nesterov momentum Adam, and the initial learning rate is set to 1×10 -4 , every 2000 training steps, the cosine annealing strategy is used to decay to 1×10 -6 The training data is dynamically augmented by the random metric tensor field in step S2, which applies a conformal transformation to the input pressure field:

[0196]

[0197] The transformation parameters a(t), b(t), c(t), and d(t) are sampled and generated from the eigenvalue distribution of the seepage field in step S2, and the sampling interval Δt=1h is synchronized with the production data update frequency.

[0198] The seepage field solution v output from step S2 FEM As a supervisory signal to drive network training; the pore space weight w provided in step S1 i Participate in loss calculation; random metric tensor g ij The statistical characteristics of the data augmentation parameters are controlled by the network. The velocity field v predicted by the network pred Feedback to the source term Q in step S2 m Perform cross validation, when the residual norm ||v pred -v FEMWhen || > 5%, the dynamic parameter update process in step S5 is triggered. The density field ρ(x, t) in the Casimir conservation term is related to the well logging data from step S1 and the pressure field in step S2 via the equation of state ρ = ρ0exp(p / K), where the bulk modulus K = 2.3 GPa is calibrated by core compression experiments.

[0199] Step S4: Quantum-macro cross-scale dynamic coupling modeling.

[0200] This step builds a unified coupling equation for the quantum adsorption effect and macroscopic seepage laws based on the neural network prediction field from step S3 and the fractional-order seepage field from step S2. By introducing a projection operator and the variational principle, a multiscale correlation is achieved, from the evolution of the electron density wave function to the flow in kilometer-scale fracture networks. This corrects the permeability prediction bias of traditional continuum theory in nanopores.

[0201] In general, the coupling model uses the Zwanzig-Mori projection form to decompose the fast and slow variables. q (r,t) represents the quantum scale fluctuation of methane molecular density, the slow variable v m (x, t) corresponds to the macroscopic velocity field predicted in step S3. The projection equation is expressed as:

[0202]

[0203] Where, the memory kernel matrix K ij (t) Calculation of the velocity autocorrelation function output by the neural network in step S3:

[0204]

[0205] K 12 (t) and K 21 (t) is the cross term, reflecting the quantum-macroscopic coupling strength; ξ q (t),ξ m (t) is a Gaussian noise term, and its covariance is obtained from the molecular dynamics trajectory statistics of step S1.

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

[0207]

[0208] Where ψ(r,t) is the wave function of methane molecules, m = 2.66×10 -26 kg; V ext (r) is the surface potential energy of organic matter, which is reconstructed by the grayscale value of the CT image in step S1, and the potential well depth V0 = 0.15 eV; μ is the equivalent inertia coefficient, which is obtained from the permeability tensor K in step S2. i Through the relationship Determine; λ = 0.32 (eV / s / m) 2 is the coupling coefficient, which is obtained by fitting the nanochannel flow rate-adsorption amount experimental curve in step S2.

[0209] Specifically, the quantum term was solved using time-dependent density functional theory (TDDFT). The exchange-correlation functional used the revPBE modified gradient approximation, with a plane wave cutoff of 400 Ry. The K-point grid was set to a 4×4×4 grid based on the periodic pore structure of the CT image in step S1. The initial wave function value ψ(r,0) was initialized using the Wannier function expansion, and the localization center position matched the organic matter distribution in the CT image in step S1.

[0210] The macro-term solution inherits the hp adaptive finite element discretization strategy of step S2. The mesh refinement trigger condition is modified to simultaneously meet the following conditions:

[0211] And δρ q >0.1ρ0

[0212] Where, ρ0=800kg / m 3 The reference density is calibrated by the well logging data in step S1. The cell size adjustment formula is modified to:

[0213]

[0214] d f is the fractal dimension calculated in step S1, Real-time output via quantum solver.

[0215] The data exchange interface uses the non-blocking communication protocol of MPI-3.0. The coupling variables of the quantum domain (CP2K) and the macro domain (OpenFOAM) are synchronized every 0.1ps. The buffer is set to double-precision floating point format, and the 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 are , triggers the real-time parameter update of step S5.

[0216] The predicted velocity field v output in step S3 m As the initial condition of the macro term; the organic matter distribution of the CT image in step S1 reconstructs the potential energy field V ext ; The hp adaptive strategy in step S2 drives the grid size h p Dynamic adjustment; fractal dimension d f Also affects the memory core K ij The decay rate of (t) 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 k ij and the coupling coefficient λ.

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

[0218] This step builds 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 integrating ensemble Kalman filtering with Bayesian inversion techniques, a coordinated optimization of the permeability field, quantum coupling coefficient, and network weights is achieved, ensuring spatiotemporal consistency between the model predictions and downhole dynamic data.

[0219] In general, the state vector is defined as X t =(k ij ,λ,W g ,φ) T Contains the permeability tensor of step S2, the coupling coefficient of step S4, the group equivariant convolution kernel weight of step S3 and the porosity field of step S1. t It consists of the distributed optical fiber pressure and temperature data of step S1 and the downhole verification value of the predicted flow rate of step S3. The sampling frequency is aligned with the quantum-macro synchronization period of step S4 and is set to 10 Hz.

[0220] In one possible implementation, the state evolution equation of the Ensemble Kalman Filter (EnKF) is written as:

[0221]

[0222] Where, is an independent Gaussian noise vector; Q t is the process noise covariance matrix, whose diagonal elements are in is the quantum density fluctuation variance in step S4, is the density gradient output in step S4; is the Jacobian sensitivity matrix, which is calculated by back-propagation gradient of the neural network in step S3; ΔR t is the coupling residual vector of step S4, when ||ΔR t ||>10 -3 The noise covariance amplification mechanism is triggered when

[0223] Specifically, the observation operator H(·) adopts a multi-scale interpolation strategy to map the seepage field pressure p(x,t) in step S2 to the wellbore trajectory:

[0224]

[0225] Among them, k m is the unit permeability of step S2; μ m T is the dynamic viscosity updated in step S3; mis the temperature interpolation field of step S1; h =25 / (1+0.05d f )m is the correlation length, d f is the fractal dimension of step S1; weight w m The feature importance coefficient output by the neural network in step S3 is determined.

[0226] For sudden geological events (crack expansion), hierarchical Bayesian inversion is triggered. Construct hierarchical posterior distribution:

[0227]

[0228] in, is a Gaussian mixture prior based on the historical output of step S4, with k=5 mixed components; R s is the observation error covariance matrix, the diagonal elements σ DAS is the standard deviation of the measurement error of the fiber data in step S1. The sampling adopts the No-U-Turn Sampler (NUTS), the maximum tree depth is set to 12, and the step size is adaptively adjusted in the range of [0.01, 0.5].

[0229] Parameter update constraints include:

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

[0231]

[0232] The coupling coefficient λ and the quantum adsorption energy E in step S4 a satisfy:

[0233] λ=0.15E a / (k B T m )

[0234] Among them, k B is the Boltzmann constant, T m is the real-time temperature of step S1.

[0235] The parallel architecture utilizes heterogeneous MPI-OpenCL computing. The number of EnKF ensemble members is set to 128, with one MPI process assigned to each member. State vector dimensionality compression is achieved through convolution kernel weight sparsification in step S3, with a compression ratio of 15:1. The data synchronization cycle is strictly aligned with the quantum-macro interface in step S4. Buffer management utilizes a double-buffered ping-pong operation, with capacity dynamically allocated based on 1.5 times the number of grid nodes in step S2.

[0236] The coupling residual ΔR output in step S4 t Driving process noise Qt Dynamic adjustment of the neural network gradient J in step S3 t Participate in the state evolution equation; the porosity field φ in step S1 is used as a state vector element to constrain the permeability update; the fractal dimension d f By the correlation length l h This affects the accuracy of observation interpolation. The inverted parameters are fed back in real time to the percolation equation in step S2, the network weights in step S3, and the coupled Hamiltonian in step S4, forming a closed-loop correction. If the Bayesian inversion takes longer than three times the data update period in step S1, the predicted values in step S3 are used to temporarily replace the real-time data until the inversion is complete.

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

[0238] This step builds a full-scale closed-loop optimization system for reservoir development based on the dynamic parameter field assimilated in step S5, the predictive model from step S3, and the cross-scale coupling constraints from step S4. By integrating stochastic dynamic programming with robust model predictive control, it achieves multi-objective coordination of nanopore adsorption control, fracture network conductivity optimization, and well pattern parameter adjustment, addressing the difficulty of traditional methods in balancing microscopic adsorption and macroscopic flow.

[0239] In general, the multi-objective optimization function is defined as:

[0240]

[0241] Where ΔQ oil (t) = q prod (t)·(1-S w (t)-S g (t)) is the daily oil production increment, S w is the water saturation inverted from the logging data in step S1, S g is the quantum density fluctuation δρ in step S4 q The converted gas phase saturation; η = 0.05 is the discount factor; ω = 0.3 is the energy consumption penalty coefficient; ζ = 10 is the flow rate fluctuation suppression coefficient, which is determined by the Bayesian inversion historical data variance in step S5.

[0242] In one possible implementation, the decision variable u t =(q inj ,P fracture ,Δt cycle ) T T includes the injection volume, fracture maintenance pressure, and injection-production cycle. The variable constraints are inherited from the seepage field pressure gradient constraint in step S2 and the quantum adsorption constraint in step S4:

[0243]

[0244] Among them, σ geo =45MPa is the formation fracture pressure obtained by inversion of the ground stress field in step S1; ρ0 = 8001g / m 3 is the reference density of step S4.

[0245] The rolling optimization adopts the improved stochastic sequential quadratic programming (SSQP). The neural network v of the prediction model integration step S3 pred =f NN (X t ) and the fractional-order seepage equation in step S2:

[0246]

[0247] Among them, Σ J is the covariance matrix of the objective function; λ = 0.7 is the risk aversion coefficient, which is dynamically adjusted by the EnKF set variance in step S5. The constrained Jacobian matrix J t The back propagation automatic differentiation calculation of the neural network in step S3 is performed, the step size search adopts the Goldstein condition, and the backtracking coefficient is set to 0.4.

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

[0249]

[0250] Where, ⊙ represents the Hadamard product;

[0251] K p =diag(0.2,0.15,0.05) is the proportional gain matrix;

[0252] K i =diag(0.1,0.08,0.03) is the integral gain matrix;

[0253] ξ=0.1h -1 is the forgetting factor; T c = 24h is the control period, which is synchronized with the data update frequency of step S1. The gain matrix is obtained by the EnKF covariance inverse matrix of step S5. Updated jointly with the sensitivity coefficient of the neural network in step S3.

[0254] The hardware deployment utilizes a heterogeneous FPGA+GPU architecture. The optimization main loop is pipelined on the FPGA, with a clock frequency set to 300 MHz and the number of pipeline stages matching the number of grid partitions in step S2 (16). Constraint checking and sensitivity calculation are offloaded to the GPU, with a CUDA kernel block size set to 256 threads, aligning with the neural network layer width in step S3. The data bus utilizes a PCIe 4.0×16 interface with a bandwidth of 64 Gbps, meeting the data throughput requirements of the distributed fiber optic sensing in step S1.

[0255] The permeability field k assimilated in step S5 ij Participate in the calculation of the velocity variance term of the objective function J; the quantum density fluctuation δρ in step S4 q Dynamically constrain the upper limit of water injection volume; the neural network prediction v in step S3 pred As the rolling optimization input; the pressure gradient in step S2 limits the decision variable boundary. The optimization results are sent to the intelligent actuator in step S1 in real time through the OPC UA protocol, forming a closed-loop link from data acquisition to control. When , the neural network online fine-tuning in step S3 and the parameter emergency inversion in step S5 are triggered.

[0256] 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. A method for identifying sweet spots in shale oil and gas reservoirs, characterized in that: The following steps are involved: Collect core CT images, 3D seismic data, and production dynamic data and perform cross-scale preprocessing; Construct a fractional-order nonlocal seepage field model to characterize flow characteristics at the nanometer to kilometer scale; Predicting seepage field parameters through neural networks constrained by Lie group symmetry; Establish a cross-scale coupling model of quantum adsorption effects and macroscopic percolation laws; Dynamically update model parameters based on real-time monitoring data; Perform multi-objective collaborative optimization to generate the three-dimensional distribution of sweet spots and development plans.

2. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 1, wherein: The pretreatment includes: Calculation of pore network fractal dimension d from core CT images f , obtained by fitting the following linear relationship: lnN(∈,δ)=-d f ln∈+C Among them, ∈ is the scanning scale, δ is the grayscale threshold, and C is the fitting constant term; Anisotropic diffusion filtering is performed on seismic data, and the filter kernel function is: Where t represents the diffusion time, and κ = 1.5σ is the noise sensitivity coefficient.

3. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 2, characterized in that: In the fractal dimension calculation: The scanning scale ∈ is set at 5 nm to 1 μm, and the grayscale threshold δ is set at CT value 200–400 HU; Fractal dimension d f The spatial fractional derivative order β satisfies β=1.25+0.2d f .

4. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 1, wherein: The fractional order seepage model satisfies: The time derivative term adopts the Caputo definition, with order α = 0.8 ± 0.05 and characteristic time τ = 0.8 s; The spatial gradient term uses the Riesz fractional-order operator with order β = 1.5 ± 0.1 and non-Darcy index γ = 0.

33.

5. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 1, characterized in that: The neural network comprises: The input layer receives the pore pressure p, permeability tensor k ij and viscosity μ; At least three layers of SL(2,R) group equivariant convolutional layers, with the number of channels in each layer ≥ 256; The output layer predicts the seepage velocity field v pred .

6. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 5, characterized in that: The loss function of the neural network is: Where, ρ0=800kg / m 3 , the integration domain Ω is consistent with the solution domain of the seepage model.

7. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 1, characterized in that: The cross-scale coupling model satisfies: Where m = 2.66 × 10 -26 kg, V ext (r) is the surface potential energy of organic matter, the potential well depth V0 = 0.15 eV, is the equivalent inertia coefficient, which is given by the permeability tensor K i Sure.

8. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 1, characterized in that: The dynamic update adopts: Ensemble Kalman Filter Algorithm Process noise Q t The diagonal element content quantum density gradient term Observation operator H(X t ) Correlation length l h =25 / (1+0.05d f )m.

9. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 1, characterized in that: The optimization control includes: Multi-objective function: Among them, η = 0.05 is the discount factor; ω = 0.3 is the energy consumption penalty coefficient; ζ = 10 is the flow rate fluctuation suppression coefficient.

10. The method for identifying sweet spots in shale oil and gas reservoirs according to claim 9, characterized in that: The optimization control also includes a control law designed as follows: Where ⊙ represents the Hadamard product, K p =diag(0.2,0.15,0.05) is the proportional gain matrix, K i =diag(0.1,0.08,0.03) is the integral gain matrix.

Images