A deep learning assisted accelerated fracturing operation parameter intelligent real-time optimization method

Abstract

The application discloses a deep learning assisted accelerated fracturing construction parameter intelligent real-time optimization method, and belongs to the field of unconventional oil and gas field development, and the steps are as follows: a structured sample database covering geological properties, engineering parameters and fracturing effect responses is constructed, and unified pretreatment is carried out on multi-source data; a multi-modal space-time collaborative prediction model is constructed, combined with multi-expansion time convolution network, three-dimensional residual network and cross attention mechanism; supervised learning training is carried out based on the database, sparse regularization and learning rate scheduling are introduced, and the model precision and generalization ability are improved; a hierarchical collaborative optimization framework of outer layer Bayesian search-inner layer CMA-ES refinement is proposed, reservoir reconstruction volume maximization and construction feasibility are realized under complex constraints; compared with the existing method, the optimization parameters of the application not only include cluster distance and section distance, but also include construction parameters such as displacement, sand concentration and fracturing fluid type, and the application can realize real-time optimization of fracturing construction parameters by using a machine learning agent model.

Description

Technical Field

[0001] This invention relates to the field of unconventional oil and gas field development technology, specifically to a deep learning-assisted, accelerated intelligent real-time optimization method for fracturing operation parameters. Background Technology

[0002] The efficient development of unconventional oil and gas resources is a key research direction in the current energy field. Fracturing, as a crucial technology for enhancing the productivity of unconventional reservoirs, exhibits complex characteristics with high dimensionality and strong nonlinearity due to the coupling effect of multiple factors, including geological parameters and fracturing operation parameters. Existing fracturing parameter optimization methods mainly rely on numerical simulation technology, optimizing fracture morphology parameters by calling post-fracturing production simulations. However, these methods typically cannot achieve real-time and intelligent optimization of operation parameters during fracturing. Such methods require the establishment of detailed geological models, consuming significant computational resources and involving long cycles. Especially when facing complex geological conditions, multi-variable collaboration, and multi-objective constraints in fracturing optimization design, traditional simulation methods are inefficient and have slow response times, making it difficult to support rapid iteration and accurate optimization of fracturing schemes in dynamic reservoir environments. Summary of the Invention

[0003] To address the aforementioned problems in existing technologies, this invention proposes a deep learning-assisted, accelerated intelligent real-time optimization method for fracturing construction parameters. This method is rationally designed, overcomes the shortcomings of existing technologies, and demonstrates excellent results.

[0004] To achieve the above objectives, the present invention adopts the following technical solution:

[0005] A deep learning-assisted, accelerated intelligent real-time optimization method for fracturing operation parameters includes the following steps:

[0006] Step 1: Based on actual geological and engineering data, construct a multi-stage horizontal well fracturing simulation model using numerical simulation; by adjusting reservoir geological property parameters and fracturing construction parameters, conduct multiple fracturing numerical simulations to obtain the corresponding reservoir stimulation volume simulation results;

[0007] Step 2: Perform unified preprocessing on multi-source data, including normalization of reservoir geological property parameters, fracturing operation parameters and reservoir stimulation volume, as well as thresholding, resampling and feature stitching of wellbore location images and natural fracture images.

[0008] Step 3: Construct a multimodal spatiotemporal collaborative prediction model, which integrates a multi-dilated temporal convolutional network, a 3D residual network with spatial attention, and a bidirectional cross-attention mechanism to achieve end-to-end high-precision real-time prediction of reservoir stimulation volume.

[0009] Step 4: Design the loss function and build a training strategy based on AdamW and cosine annealing. Introduce structured sparse regularization constraints to enhance the identification of key parameters, and combine adaptive optimization and learning rate scheduling to comprehensively improve prediction accuracy and generalization ability.

[0010] Step 5: Integrate the outer Bayesian global search, the inner CMA-ES continuous refinement, and the discrete neighborhood switching into a collaborative optimization framework, and call the objective function to quickly evaluate the candidate solutions, so as to achieve the organic unity of global direction, local refinement and category exit, thereby obtaining a real-time optimized fracturing construction scheme that takes into account both reservoir stimulation volume maximization and engineering workability under complex constraints.

[0011] Furthermore, in step 1, the reservoir geological property parameters include minimum horizontal principal stress, maximum horizontal principal stress, vertical stress, porosity, permeability, Young's modulus, and Poisson's ratio;

[0012] The fracturing operation parameters include the pumping rate, fracturing fluid type, proppant type, sand concentration parameters, fluid volume, and cluster spacing configuration within the horizontal well section.

[0013] Furthermore, in step 2, the StandardScaler function is used to process the reservoir geological property parameters and fracturing operation parameters, and the MinMaxScaler function is used to process the reservoir stimulation volume results to eliminate the adverse effects caused by singular sample data.

[0014] Both the wellbore location image and the natural fracture distribution image are converted into grayscale images. A threshold judgment method is used to judge the grayscale values ​​of the wellbore location image and the natural fracture distribution image, and grayscale value classification is performed. The wellbore location image and the natural fracture distribution image are resampled to the same size as the reservoir attribute spatial resolution. The features of the processed wellbore location image and the natural fracture distribution image are stitched with the reservoir attributes to obtain the reservoir fusion feature data, thereby expanding the feature dimension of the reservoir data volume.

[0015] Furthermore, step 3 includes the following sub-steps:

[0016] Step 3.1: Based on the temporal characteristics of fracturing operation parameters, a multi-dilation temporal convolutional network is constructed to extract dynamic features, specifically as follows:

[0017] Step 3.1.1: Extract features from fracturing operation parameters and construct a tensor of shape [B,24,11], where B is the batch size, 24 is the time step, and 11 is the fracturing operation feature dimension;

[0018] Step 3.1.2: Map the 11-dimensional features to a 32-dimensional high-dimensional space through a fully connected layer; then, perform normalization using LayerNorm and apply the ReLU activation function for non-linear activation;

[0019] Step 3.1.3: Design TCNBlock. Each TCNBlock uses Conv1d to perform one-dimensional convolution operations with a kernel size of 3, and applies dilated convolution to expand the receptive field. Each convolutional layer is processed by batch normalization, ReLU activation function, Dropout, and residual connections.

[0020] Step 3.1.4: Stack 3 TCNBlocks, with expansion rates of 1, 2, and 4 respectively, to gradually increase the receptive field in order to capture dependencies at different time scales;

[0021] Step 3.1.5: After processing by 3 TCNBlocks, the time-series feature tensor is output with shape [B,32,24]. Mean pooling and max pooling are used to extract the global features of the time-series feature tensor, and the two pooling results are concatenated into a feature vector. Finally, the fully connected layer is used for further processing to obtain the fracturing construction time-series feature vector with shape [B,32].

[0022] Step 3.2: Based on the spatial structure of reservoir geological attributes, a three-dimensional convolutional residual network is constructed to extract spatial features, specifically as follows:

[0023] Step 3.2.1: Extract features from reservoir fusion feature data and construct a tensor of shape [B,9,5,239,76], where 9 is the feature dimension, 5 is the spatial dimension in the Z direction, and 239 and 76 represent the spatial dimensions in the X and Y directions;

[0024] Step 3.2.2: Introduce a 3D convolutional layer with a kernel size of 3×3×3 and padding=1 to map the input 9-channel features to a 16-channel output, resulting in a feature tensor shape of [B,16,5,239,76].

[0025] Step 3.2.3: Construct a two-layer 3D residual block. Each residual block contains two layers of 3D convolutions with a kernel size of 3×3×3. Batch normalization and GELU activation functions are added after each layer. The residual block adopts a residual connection structure with identity mapping. When the number of input and output channels is inconsistent, channel matching is achieved through 1×1×1 convolution to ensure consistent dimensions. The first residual block expands the number of channels from 16 to 32, and the second residual block further expands it to 64. The final output feature tensor shape is [B,64,5,239,76].

[0026] Step 3.2.4: Introduce a spatial attention mechanism. Calculate the global context of the feature tensor output by the two-layer 3D residual block through local average pooling. Then, use 1×1×1 convolution and the Sigmoid function to generate a spatial attention weight map. Weight the input features point by point to obtain the weighted spatial feature tensor with the shape [B,64,5,239,76].

[0027] Step 3.2.5: Employing a multi-scale feature aggregation strategy, global average pooling is applied to the weighted spatial feature tensor to reduce the spatial dimension to a single global feature vector, with an output shape of [B, 64]. Spatial pooling is then applied to the weighted spatial feature tensor to retain some spatial information, resulting in an output shape of [B, 1152]. The global pooling feature and the spatial pooling feature are then concatenated to obtain a fused feature of [B, 1216]. After processing through a fully connected layer, layer normalization is performed to compress the fused feature dimension to [B, 128], outputting the reservoir spatial feature vector.

[0028] Step 3.3: Introduce a bidirectional cross-attention mechanism to model the nonlinear interaction between fracturing operations and reservoir properties, specifically:

[0029] Step 3.3.1: Map the fracturing operation timing characteristics and reservoir spatial characteristics to a unified 64-dimensional space through linear projection;

[0030] Step 3.3.2: Feature interaction is performed through a multi-head attention mechanism, which is divided into two paths. The first path uses reservoir spatial features as the query and fracturing operation time sequence features as the key and value, to extract response features related to key areas of the reservoir from the fracturing operation dynamics. The second path uses fracturing operation features as the query and reservoir features as the key and value, to capture the influence intensity of different fracturing operation stages on the reservoir spatial area.

[0031] Step 3.3.3: Employ a multi-head attention mechanism to combine different subspaces. For each query vector in Step 3.3.2, multiple attention heads will learn different attention weights and perform a weighted average within each head to obtain two 64-dimensional output features.

[0032] Step 3.3.4: The two 64-dimensional output features are respectively processed by linear transformation, GELU activation function, and Dropout, and finally passed through a linear layer to output two 64-dimensional interaction feature vectors;

[0033] Step 3.3.5: Concatenate the two interactive feature vectors obtained in step 3.3.4 to obtain a 128-dimensional feature vector with shape [B, 128].

[0034] Step 3.4: The reservoir spatial characteristics of [B, 128], the fracturing operation timing characteristics of [B, 128], and the bidirectional attention interaction characteristics of [B, 128] are fused and spliced ​​to obtain the multi-source spliced ​​tensor of [B, 288]. The multi-source spliced ​​tensor is compressed to a low-dimensional representation through a fully connected layer, with the shape [B, 64]. Layer normalization and ReLU activation function processing are then performed. Finally, the prediction head is used to map this representation to the final reservoir stimulation volume prediction value.

[0035] Furthermore, step 4 includes the following sub-steps:

[0036] Step 4.1: Use mean squared error as the main loss function, and add L1 sparse regularization to the feature projection layer of the encoder of construction parameters to ensure that the model can automatically identify and focus on key construction parameters. The expression is:

[0037] (3)

[0038] in, Mean square error, Indicates the first One sample, To modify the predicted volume, To modify the actual volume value, The number of samples;

[0039] (4)

[0040] in, For sparse regularization terms, For the model One parameter, The regularization coefficient is used.

[0041] Step 4.2: Using the AdamW optimizer and cosine annealing learning rate scheduling mechanism, the learning rates of the fracturing operation time encoder and the reservoir space encoder are set differently for model training.

[0042] Furthermore, step 5 includes the following sub-steps:

[0043] Step 5.1: Establish the hard constraint optimization range of fracturing construction parameters, convert the fracturing construction parameters into a hybrid vector according to the stage hard constraints, and generate multiple sets of feasible solutions by random sampling within the feasible range;

[0044] Step 5.2: Use the outer Bayesian global search to find all feasible solutions, score all feasible solutions using the sampling function, and use phantom sampling to select the p candidate solutions with the best scores as the output of the outer global search;

[0045] Step 5.3: Press The top q optimal solutions are selected from high to low, and these q optimal solutions are respectively entered into the inner local optimization. Continuous variables are adaptively shrunk and refined within the trust region using CMA-ES, and discrete variables are flipped in a small range using the neighborhood switching operator.

[0046] Step 5.4: Input the q optimal solutions optimized from the inner layer into the multimodal spatiotemporal collaborative prediction model, and obtain a new global optimal solution by calculating and comparing the objective function values. Let the objective function values ​​of the new global optimal solution and the current global optimal solution be denoted as follows: and If satisfied ,in To minimize the improvement threshold, the global optimal solution is updated, and then the process returns to step 5.2 with the updated solution. Recalculate the acquisition function and repeat the above process; if the conditions are not met, return to step 5.2, add the selected candidate solutions to the forbidden acquisition set, and continue the subsequent process; until the loop reaches the upper limit of iterations, output the current value. The corresponding fracturing construction parameters.

[0047] Furthermore, step 5.1 includes the following sub-steps:

[0048] Step 5.1.1: Based on the fracturing process, the fracturing construction parameters are divided into four key stages: pre-flush fluid stage, initial stage of proppant-carrying fluid stage, late stage of proppant-carrying fluid stage, and displacement fluid stage. Hard constraints are applied to the variables of each stage to optimize the range. The variables include fracturing fluid type, proppant type, pumping rate, proppant concentration, and fluid volume.

[0049] Step 5.1.2: Construct shape functions for pumping rate, sand concentration, and liquid volume for each stage. Use the starting value, ending value, and trend parameters to generate the pumping rate curve, sand concentration change curve, and liquid volume distribution curve for each stage. Then, map the curves to the hard constraint optimization range through interval mapping. Integrate the cluster spacing, stage discrete variables, and shape functions. Randomly sample within the hard constraint optimization range to obtain a batch of mixed variables H, where the stage discrete variables are fracturing fluid type and proppant type.

[0050] Step 5.1.3: Expand the mixture vector H into specific parameters for 24 time steps, and perform interval pruning, sand banning, support banning, adjacent jump amplitude limiting, and consistency fine-tuning on those that do not meet the hard constraint optimization range, and record metadata. This includes the extent of repair, the number of times the limit was exceeded, the proportion of sand transported beyond the boundary, and the number of times the category was corrected.

[0051] Furthermore, in step 5.2, designing the acquisition function includes the following sub-steps:

[0052] Step 5.2.1: Design the objective function and introduce a multidimensional penalty term, specifically as follows:

[0053] The reservoir stimulation volume and engineering constructability are combined into a single objective function. :

[0054] (5)

[0055] in, This is the output of the multimodal spatiotemporal collaborative prediction model. The fracturing parameters are for H after repair. The penalty coefficient is... A soft penalty for the rationality of fracturing construction parameters;

[0056] Four types of penalty functions are designed to suppress and regulate inappropriate behavior in the objective function:

[0057] (6)

[0058] in, The weight for adjacent transition penalties. The adjacent transition penalty function, Weights are used to penalize smoothness. Here is the smoothness penalty function. Penalty weight for sand transport capacity, Let the sand transport capacity penalty function be... Weights are used to penalize complexity. This is a complexity penalty function;

[0059] Design an adjacent jump penalty to limit the rate and sand concentration variation between adjacent time steps:

[0060] (7)

[0061] in, Index for the stage; Representative stage A single time step index within; For the stage The set of all time step indices in the middle; This represents the difference in pumping rate between adjacent time steps. For the stage Threshold for maximum pumping rate variation; This represents the difference in sand concentration between adjacent time steps. For the stage The threshold for the maximum sand concentration change in the middle; Representation phase Does it contain sand? Indicates the addition of sand. This indicates that no sand is added;

[0062] Design a smoothness penalty to limit the rate and frequency of sand concentration changes between adjacent time steps:

[0063] (8)

[0064] in, For the stage The maximum range of change in medium speed; For the stage Maximum variation range of medium sand concentration; For the stage The maximum range of variation in liquid volume; The weighting coefficients are used to adjust the total variation and the second difference weights.

[0065] The design imposes a penalty on the sand conveying capacity, limiting the instantaneous sand conveying rate to no more than the equipment's load capacity.

[0066] (9)

[0067] in, For the stage The upper limit of medium sand conveying capacity indicates the maximum sand conveying rate that the equipment or wellbore can withstand. For sand concentration, This refers to the pumping rate;

[0068] Design complexity penalty to avoid overly complex fracturing construction schemes:

[0069] (10)

[0070] (11)

[0071] (12)

[0072] in, For sand concentration changes Sign judgment, The threshold for sand concentration change, when hour, ,when hour, ,when hour, ; For determining the sign of the change in pumping rate, The threshold for the change in pumping rate, when hour, ,when and hour, ,when and hour, ; For the stage The number of peak concentrations of medium sand; For the stage The number of inflection points at medium speed; For the stage Maximum number of peak values ​​allowed in the middle; For the stage The maximum number of inflection points allowed in the middle;

[0073] Step 5.2.2: Utilize the metadata from Step 5.1.3 The feasibility model is trained using a gradient boosting decision tree, and the probability that the mixed variable H satisfies the hard constraints set in step 5.1.1 is output.

[0074] Step 5.2.3: with To achieve the objective, a mixed-variable Bayesian optimization model is established, and a multimodal spatiotemporal collaborative prediction model is used to calculate the predicted mean from the output. ,variance Provides sampling and decision-making support for mixed-variable Bayesian optimization models;

[0075] Step 5.2.4: Combine with desired improvements The probability output by the feasibility model Design a constraint-sensing acquisition function :

[0076] (13)

[0077] (14)

[0078] in, It is the cumulative distribution function of the standard normal distribution; It is the probability density function of the standard normal distribution; This is the current optimal solution; To explore constants; It is a potential improvement after standardization.

[0079] Furthermore, step 5.3 includes the following sub-steps:

[0080] Step 5.3.1: From the p candidate solutions obtained in the outer optimization, sort them according to... Select the top q optimal solutions from highest to lowest;

[0081] Step 5.3.2: Expand the q optimal solutions in the outer optimization into specific time step parameters, and use the Sigmoid function to automatically repair abnormal parameters;

[0082] Step 5.3.3: Use the CMA-ES algorithm to perform local optimization on the continuous variables of the q optimal solutions, and adjust the search strategy through the adaptive trust region algorithm, where the initial step size is set to 0.3, the trust region radius is set to 20% of the interval width, and the population size is set to 5.

[0083] Step 5.3.4: Use the neighborhood switching operator to flip the local class of the discrete variable.

[0084] The beneficial technical effects of this invention are as follows:

[0085] (i) Achieve high-precision modeling of geological-engineering multimodal characteristics and improve the reliability of reservoir stimulation volume prediction;

[0086] Overcoming the limitations of traditional models in accurately depicting the coupling relationship between geological static attributes and engineering dynamic parameters, a high-precision model for predicting reservoir transformation volume was constructed using multimodal feature extraction and interactive modeling techniques. A 3D convolutional-residual network was employed to deeply mine the spatial distribution characteristics of reservoir attribute volumes, combined with a multi-expansion rate temporal convolutional network to capture the multi-scale dynamic changes in the pumping procedure. Furthermore, a bidirectional cross-attention mechanism was used to fully integrate key information from geological and engineering modes, addressing the shortcomings of traditional data-driven models in insufficient consideration of multi-source parameter coupling and weak generalization.

[0087] (ii) Construct an intelligent optimization architecture for construction parameters to achieve intelligent real-time optimization of construction parameters;

[0088] To address the shortcomings of traditional numerical simulation optimization, such as high computational resource consumption and long cycles, and the inability to perform coarse-grained optimization of global indicators while failing to achieve time-step collaborative optimization of multi-dimensional construction parameters, this invention proposes an intelligent optimization architecture. The outer layer employs constraint-aware Bayesian optimization to achieve global exploration of the mixed variable space, while the inner layer combines CMA-ES for fine-tuning of continuous variables and handles discontinuous variables through discrete neighborhood switching. This method can optimize multi-dimensional parameters such as liquid ratio, pumping rate, and sand concentration in real time at each time step, ensuring the rationality and stability of the construction scheme throughout its dynamic evolution.

[0089] (iii) Achieving synergy between "maximizing the transformation effect and engineering constraint feasibility" and possessing practical value on site;

[0090] This invention focuses on maximizing reservoir stimulation volume. It integrates the rigid constraints of the fracturing process into the optimization engine through multiple constraint penalty functions, ensuring that the output scheme conforms to reservoir characteristics and the operating range of field equipment, thus avoiding the problem of theoretically optimal solutions that are impractical for construction. Furthermore, the scheme design does not rely on complex geological modeling and high-configuration computing resources; it can be quickly generated directly based on field geological and construction data, facilitating understanding and implementation by field engineers and providing technical support for the efficient development of unconventional oil and gas resources. Attached Figure Description

[0091] Figure 1 This is a flowchart of a deep learning-assisted, accelerated intelligent real-time optimization method for fracturing construction parameters in this invention.

[0092] Figure 2 This is a schematic diagram of the multi-dilation rate temporal convolutional network in this invention.

[0093] Figure 3 This is a schematic diagram of the 3D residual network with spatial attention in this invention.

[0094] Figure 4 This is a schematic diagram of the multimodal spatiotemporal collaborative prediction model in this invention.

[0095] Figure 5 This is the multi-scale hybrid variable optimization architecture of outer Bayesian layer, inner CMA-ES layer, and discrete neighborhood switching in this invention.

[0096] Figure 6 Accuracy map of the multimodal spatiotemporal collaborative prediction model for reservoir stimulation volume.

[0097] Figure 7 This is a graph validating the performance of a multi-scale mixed variable intelligent optimization model. Detailed Implementation

[0098] The specific embodiments of the present invention will be further described below with reference to specific examples:

[0099] A deep learning-assisted, accelerated intelligent real-time optimization method for fracturing operation parameters, such as... Figure 1 As shown, it includes the following steps:

[0100] Step 1: Based on actual geological and engineering data, construct a multi-stage horizontal well fracturing simulation model using numerical simulation; conduct multiple fracturing numerical simulations by adjusting reservoir geological property parameters and fracturing construction parameters, obtain the corresponding reservoir stimulation volume simulation results, and establish a structured database covering geological properties, engineering parameters, and reservoir stimulation volume.

[0101] Based on logging data from multiple actual wells, multiple horizontal well fracturing models were constructed using a numerical simulator. Key geological attribute parameters of different reservoir layers were extracted. These reservoir geological attribute parameters include minimum horizontal principal stress, maximum horizontal principal stress, vertical stress, porosity, permeability, Young's modulus, and Poisson's ratio. The data were saved in .npy format to establish a horizontal well reservoir attribute dataset.

[0102] The fracturing operation parameters include the pumping rate, fracturing fluid type, proppant type, sand concentration parameters, fluid volume and cluster spacing configuration in the horizontal well section. 100 fracturing pumping schemes are designed and saved as .xlsx format data to establish a horizontal well fracturing operation design dataset.

[0103] Step 2: Perform unified preprocessing on multi-source data, including normalization of reservoir geological property parameters, fracturing operation parameters and reservoir stimulation volume, as well as thresholding, resampling and feature stitching of wellbore location images and natural fracture images.

[0104] The StandardScaler function is used to process reservoir geological property parameters and fracturing operation parameters to eliminate the adverse effects caused by singular sample data in reservoir parameter properties and fracturing pump design. The expression is as follows:

[0105] (1)

[0106] in, Different parameters are designed for different reservoir properties (or pumping). The normalized value; Different parameters are designed for various reservoir properties (or for pumping). For data The mean; For data Standard deviation;

[0107] The MinMaxScaler function is used to process reservoir stimulation volume results, eliminating the adverse effects caused by outlier data. The expression is as follows:

[0108] (2)

[0109] in, For reservoir stimulation volume The normalized value; For reservoir stimulation volume; For reservoir stimulation volume The minimum value; For reservoir stimulation volume The maximum value;

[0110] Both the wellbore borehole location image and the natural fracture distribution image are converted into grayscale images, and corresponding grayscale images of the wellbore borehole location and the natural fracture distribution are obtained.

[0111] A threshold judgment method was used to determine the grayscale values ​​of the wellbore borehole location grayscale image and the natural fracture distribution grayscale image, and then perform grayscale value classification processing. Specifically:

[0112] In the grayscale image of the wellbore hole location, values ​​less than 35 are converted to 0, values ​​greater than or equal to 35 and less than 55 are converted to 2, and values ​​greater than or equal to 55 are converted to 1. Here, 0 represents the background, 1 represents the wellbore location, and 2 represents the first hole in the segment.

[0113] In the grayscale image of natural crack distribution, values ​​with a grayscale value less than 35 are converted to 0, and values ​​greater than or equal to 35 are converted to 3, where 0 represents the background and 3 represents natural cracks.

[0114] The grayscale images of wellbore locations and natural fracture distributions are resampled to the same size as the spatial resolution of reservoir attributes. The processed images of wellbore locations and natural fracture distributions are then stitched together with the reservoir attributes to obtain reservoir fusion feature data, thereby expanding the feature dimensions of the reservoir data volume. In this embodiment, the feature dimensions of the reservoir data volume are expanded from (7,5,239,76) to (9,5,239,76).

[0115] Step 3: Construct a multimodal spatiotemporal collaborative prediction model, which integrates a multi-dilated temporal convolutional network, a 3D residual network with spatial attention, and a bidirectional cross-attention mechanism to achieve end-to-end high-precision real-time prediction of reservoir stimulation volume.

[0116] Includes the following sub-steps:

[0117] Step 3.1: Based on the temporal characteristics of fracturing operation parameters, construct a multi-expansion temporal convolutional network to extract dynamic features, such as... Figures 2-4 As shown, specifically:

[0118] Step 3.1.1: Extract features from fracturing operation parameters and construct a tensor of shape [B,24,11], where B is the batch size, 24 is the time step, and 11 is the fracturing operation feature dimension. Since the fracturing fluid type and proppant type are one-hot encoded, the fracturing fluid type has 4 dimensions and the proppant type has 3 dimensions, so there are a total of 11 feature dimensions for the 6 operation parameters.

[0119] Step 3.1.2: Map the 11-dimensional features to a 32-dimensional high-dimensional space through a fully connected layer; then, perform normalization using LayerNorm and apply the ReLU activation function for non-linear activation;

[0120] Step 3.1.3: Design TCNBlock. Each TCNBlock uses Conv1d to perform one-dimensional convolution operations with a kernel size of 3, and applies dilated convolution to expand the receptive field. Each convolutional layer is processed by batch normalization, ReLU activation function, Dropout, and residual connections.

[0121] Step 3.1.4: Stack 3 TCNBlocks, with expansion rates of 1, 2, and 4 respectively, to gradually increase the receptive field in order to capture dependencies at different time scales;

[0122] Step 3.1.5: After processing by 3 TCNBlocks, the time-series feature tensor is output with shape [B,32,24]. Mean pooling and max pooling are used to extract the global features of the time-series feature tensor, and the two pooling results are concatenated into a feature vector. Finally, the fully connected layer is used for further processing to obtain the fracturing construction time-series feature vector with shape [B,32].

[0123] Step 3.2: Based on the spatial structure of reservoir geological attributes, a three-dimensional convolutional residual network is constructed to extract spatial features, specifically as follows:

[0124] Step 3.2.1: Extract features from reservoir fusion feature data and construct a tensor of shape [B,9,5,239,76], where 9 is the feature dimension, 5 is the spatial dimension in the Z direction, and 239 and 76 represent the spatial dimensions in the X and Y directions;

[0125] Step 3.2.2: Introduce a 3D convolutional layer with a kernel size of 3×3×3 and padding=1 to map the input 9-channel features to a 16-channel output, resulting in a feature tensor shape of [B,16,5,239,76].

[0126] Step 3.2.3: Construct two-layer 3D residual blocks. Each residual block contains two layers of 3D convolutions with a kernel size of 3×3×3. Batch normalization and GELU activation functions are added after each layer. The residual blocks adopt a residual connection structure with identity mapping. When the number of input and output channels is inconsistent, 1×1×1 convolutions are used to achieve channel matching to ensure dimension consistency, thereby ensuring efficient gradient propagation and alleviating the degradation problem in deep training. During this process, the network gradually expands the number of channels: the first residual block expands the number of channels from 16 to 32, and the second residual block further expands it to 64. The final output feature tensor shape is [B,64,5,239,76].

[0127] Step 3.2.4: Introduce a spatial attention mechanism. Calculate the global context of the feature tensor output by the two-layer 3D residual block using local average pooling. Then, use 1×1×1 convolution and the Sigmoid function to generate a spatial attention weight map. Weight the input features point by point to obtain a weighted spatial feature tensor with the shape [B,64,5,239,76], thus enhancing the feature response of key regions.

[0128] Step 3.2.5: Employing a multi-scale feature aggregation strategy, global average pooling is applied to the weighted spatial feature tensor to reduce the spatial dimension to a single global feature vector, with an output shape of [B, 64]. Spatial pooling is then applied to the weighted spatial feature tensor to retain some spatial information, resulting in an output shape of [B, 1152]. The global pooling feature and the spatial pooling feature are then concatenated to obtain a fused feature of [B, 1216]. After processing through a fully connected layer, layer normalization is performed to compress the fused feature dimension to [B, 128], outputting the reservoir spatial feature vector.

[0129] Step 3.3: Introduce a bidirectional cross-attention mechanism to model the nonlinear interaction between fracturing operations and reservoir properties, specifically:

[0130] Step 3.3.1: Map the fracturing operation time sequence characteristics (32-dimensional) and reservoir spatial characteristics (128-dimensional) to a unified 64-dimensional space through linear projection;

[0131] Step 3.3.2: Feature interaction is performed through a multi-head attention mechanism, which is divided into two paths. The first path uses reservoir spatial features as the query and fracturing operation time sequence features as the key and value, to extract response features related to key areas of the reservoir from the fracturing operation dynamics. The second path uses fracturing operation pumping features as the query and reservoir features as the key and value, to capture the intensity of the impact of different fracturing operation stages on the reservoir spatial area.

[0132] Step 3.3.3: Employ a multi-head attention mechanism to combine different subspaces. For each query vector in Step 3.3.2, multiple attention heads will learn different attention weights and perform a weighted average within each head to obtain two 64-dimensional output features.

[0133] Step 3.3.4: The two 64-dimensional output features are respectively processed by linear transformation, GELU activation function, and Dropout, and finally passed through a linear layer to output two 64-dimensional interaction feature vectors;

[0134] Step 3.3.5: Concatenate the two interactive feature vectors (each with a dimension of 64) obtained in step 3.3.4 to obtain a 128-dimensional feature vector with richer information and a shape of [B, 128].

[0135] Step 3.4: The reservoir spatial characteristics of [B, 128], the fracturing operation timing characteristics of [B, 128], and the bidirectional attention interaction characteristics of [B, 128] are fused and spliced ​​to obtain the multi-source spliced ​​tensor of [B, 288]. The multi-source spliced ​​tensor is compressed to a low-dimensional representation through a fully connected layer, with the shape [B, 64]. Layer normalization and ReLU activation function processing are then performed. Finally, the prediction head is used to map this representation to the final reservoir stimulation volume prediction value.

[0136] Step 4: Design the loss function and build a training strategy based on AdamW and cosine annealing. Introduce structured sparse regularization constraints to enhance the identification of key parameters, and combine adaptive optimization and learning rate scheduling to comprehensively improve prediction accuracy and generalization ability.

[0137] Includes the following sub-steps:

[0138] Step 4.1: Use mean squared error as the main loss function, and add L1 sparse regularization to the feature projection layer of the encoder of construction parameters to ensure that the model can automatically identify and focus on key construction parameters. The expression is:

[0139] (3)

[0140] in, Mean square error, Indicates the first One sample, To modify the predicted volume, To modify the actual volume value, The number of samples;

[0141] (4)

[0142] in, For sparse regularization terms, For the model One parameter, The regularization coefficient is used.

[0143] Step 4.2: Using the AdamW optimizer and cosine annealing learning rate scheduling mechanism, the learning rates of the fracturing operation time encoder and the reservoir space encoder are set differently for model training.

[0144] Includes the following sub-steps:

[0145] Step 4.2.1: Use the AdamW optimizer with the following configuration: initial learning rate 0.0005, weight decay 0.0005, momentum parameters β1 and β2 0.9 and 0.999 respectively, and smoothing parameter 10⁻⁸.

[0146] Step 4.2.2: Use the cosine annealing learning rate scheduling mechanism with the following configuration: initial learning rate 0.0005, minimum learning rate 0.00001, initial period 20, and period multiplication factor 2.

[0147] Step 4.2.3: Differentiate the learning rates of the pumping timing encoder and the reservoir space encoder to 0.0005 and 0.0001, respectively;

[0148] Step 4.2.4: Set the batch size to 16, use 80% of the data as the training set and 20% as the validation set, set the number of training rounds to 300, and add an early stopping mechanism to train the model;

[0149] Step 5: Integrate the outer Bayesian global search, the inner CMA-ES continuous refinement, and the discrete neighborhood switching into a collaborative optimization framework, such as... Figure 5 As shown, the objective function is called to quickly evaluate the candidate solutions, achieving an organic unity of global direction, local refinement and category exit, thereby obtaining a real-time optimized fracturing construction scheme that balances reservoir stimulation volume maximization and engineering workability under complex constraints.

[0150] Includes the following sub-steps:

[0151] Step 5.1: Establish the hard constraint optimization range of fracturing construction parameters, convert the fracturing construction parameters into a hybrid vector according to the stage hard constraints, and generate multiple sets of feasible solutions by random sampling within the feasible range;

[0152] Step 5.1.1: Based on the fracturing process, the fracturing operation parameters are divided into four key stages: pre-flush fluid stage, initial proppant-carrying fluid stage, late proppant-carrying fluid stage, and displacement fluid stage. The distinction between the initial and late proppant-carrying fluid stages is based on the amount of proppant added; the amount added is smaller in the early stage and larger in the later stage. Hard constraints are applied to the optimization range of variables for each stage. These variables include fracturing fluid type, proppant type, pumping rate, proppant concentration, and fluid volume. The hard constraint optimization range is shown in Table 1, specifically:

[0153] During the pre-fracturing stage, the fracturing fluid type is hydrochloric acid, slickwater, or cross-linked adhesive; the proppant type is not specified; and the pumping rate is 2–6 m / s. 3 / min, sand concentration 0, liquid volume 10~30 m 3 ;

[0154] In the initial stage of proppant delivery, the fracturing fluid type is slickwater or guar gum, the proppant type is 70 / 140 mesh quartz sand or 40 / 70 mesh quartz sand, and the pumping rate is 6~18 m. 3 / min, sand concentration 0~300 kgPA, liquid volume 30~60 m³ 3 ;

[0155] In the later stage of proppant-carrying, the fracturing fluid type is slickwater or guar gum, the proppant type is 70 / 140 mesh quartz sand or 40 / 70 mesh quartz sand, and the pumping rate is 8~2 m. 3 / min, sand concentration 0~300 kgPA, liquid volume 50~100 m³ 3 ;

[0156] In the later stage of proppant-carrying, the fracturing fluid type is slickwater or guar gum, the proppant type is not specified, and the pumping rate is 2~16m. 3 / min, sand concentration 0, liquid volume 0~100 m 3 ;

[0157] Table 1. Optimization range for each stage with hard constraints applied;

[0158] ;

[0159] Step 5.1.2: Construct shape functions for pumping rate, sand concentration, and liquid volume for each stage. Use the starting value, ending value, and trend parameters to generate the pumping rate curve, sand concentration change curve, and liquid volume distribution curve for each stage. Then, map the curves to the hard constraint optimization range through interval mapping. Integrate the cluster spacing, stage discrete variables, and shape functions. Randomly sample within the hard constraint optimization range to obtain a batch of mixed variables H, where the stage discrete variables are fracturing fluid type and proppant type.

[0160] Step 5.1.3: Expand the mixture vector H into specific parameters for 24 time steps, and perform interval pruning, sand banning, support banning, adjacent jump amplitude limiting, and consistency fine-tuning on those that do not meet the hard constraint optimization range, and record metadata. This includes the extent of repair, the number of times the limit was exceeded, the proportion of sand transported beyond the boundary, and the number of times the category was corrected.

[0161] Step 5.2: Use the outer Bayesian global search to find all feasible solutions, score all feasible solutions using the sampling function, and use phantom sampling to select the p candidate solutions with the best scores as the output of the outer global search;

[0162] Designing the data acquisition function includes the following sub-steps:

[0163] Step 5.2.1: Design the objective function and introduce a multidimensional penalty term, specifically as follows:

[0164] The reservoir stimulation volume and engineering constructability are combined into a single objective function. :

[0165] (5)

[0166] in, This is the output of the multimodal spatiotemporal collaborative prediction model. The fracturing parameters are for H after repair. The penalty coefficient is... A soft penalty for the rationality of fracturing construction parameters;

[0167] Four types of penalty functions are designed to suppress and regulate inappropriate behavior in the objective function:

[0168] (6)

[0169] in, The weight for adjacent transition penalties. The adjacent transition penalty function, Weights are used to penalize smoothness. Here is the smoothness penalty function. Penalty weight for sand transport capacity, Let the sand transport capacity penalty function be... Weights are used to penalize complexity. This is a complexity penalty function;

[0170] Design an adjacent jump penalty to limit the rate and sand concentration variation between adjacent time steps:

[0171] (7)

[0172] in, Index for the stage; Representative stage A single time step index within; For the stage The set of all time step indices in the middle; This represents the difference in pumping rate between adjacent time steps. For the stage Threshold for maximum pumping rate variation; This represents the difference in sand concentration between adjacent time steps. For the stage The threshold for the maximum sand concentration change in the middle; Representation phase Does it contain sand? Indicates the addition of sand. This indicates that no sand is added;

[0173] Design a smoothness penalty to limit the rate and frequency of sand concentration changes between adjacent time steps:

[0174] (8)

[0175] in, For the stage The maximum range of change in medium speed; For the stage Maximum variation range of medium sand concentration; For the stage The maximum range of variation in liquid volume; The weighting coefficients are used to adjust the total variation and the second difference weights.

[0176] The design imposes a penalty on the sand conveying capacity, limiting the instantaneous sand conveying rate to no more than the equipment's load capacity.

[0177] (9)

[0178] in, For the stage The upper limit of medium sand conveying capacity indicates the maximum sand conveying rate that the equipment or wellbore can withstand. For sand concentration, This refers to the pumping rate;

[0179] Design complexity penalty to avoid overly complex fracturing construction schemes:

[0180] (10)

[0181] (11)

[0182] (12)

[0183] in, For sand concentration changes Sign judgment, The threshold for sand concentration change, when hour, ,when hour, ,when hour, ; For determining the sign of the change in pumping rate, The threshold for the change in pumping rate, when hour, ,when and hour, ,when and hour, ; For the stage The number of peak concentrations of medium sand; For the stage The number of inflection points at medium speed; For the stage Maximum number of peak values ​​allowed in the middle; For the stage The maximum number of inflection points allowed in the middle;

[0184] Step 5.2.2: Utilize the metadata from Step 5.1.3 The feasibility model is trained using a gradient boosting decision tree, and the probability that the mixed variable H satisfies the hard constraints set in step 5.1.1 is output.

[0185] Step 5.2.3: with To achieve the objective, a mixed-variable Bayesian optimization model is established, and a multimodal spatiotemporal collaborative prediction model is used to calculate the predicted mean from the output. ,variance Provides sampling and decision-making support for mixed-variable Bayesian optimization models;

[0186] Step 5.2.4: Combine with desired improvements The probability output by the feasibility model Design a constraint-sensing acquisition function :

[0187] (13)

[0188] (14)

[0189] in, It is the cumulative distribution function of the standard normal distribution; It is the probability density function of the standard normal distribution; This is the current optimal solution; To explore the constant, it is usually taken as 0.01; It is a potential improvement after standardization.

[0190] Step 5.3: Press The top q optimal solutions are selected from high to low, and these q optimal solutions are respectively entered into the inner local optimization. Continuous variables are adaptively shrunk and refined within the trust region using CMA-ES, and discrete variables are flipped in a small range using the neighborhood switching operator.

[0191] Step 5.3.1: From the p candidate solutions obtained in the outer optimization, sort them according to... Select the top q optimal solutions from highest to lowest;

[0192] Step 5.3.2: Expand the q optimal solutions in the outer optimization into specific time step parameters, and use the Sigmoid function to automatically repair abnormal parameters;

[0193] Step 5.3.3: Use the CMA-ES algorithm to perform local optimization on the continuous variables of the q optimal solutions, and adjust the search strategy through the adaptive trust region algorithm, where the initial step size is set to 0.3, the trust region radius is set to 20% of the interval width, and the population size is set to 5;

[0194] Step 5.3.4: Use the neighborhood switching operator to flip the local class of the discrete variable.

[0195] Step 5.4: Input the q optimal solutions optimized from the inner layer into the multimodal spatiotemporal collaborative prediction model, and obtain a new global optimal solution by calculating and comparing the objective function values. Let the objective function values ​​of the new global optimal solution and the current global optimal solution be denoted as follows: and If satisfied ,in To minimize the improvement threshold, update the global optimal solution. Then return to step 5.2 with the updated solution. Recalculate the acquisition function and repeat the above process. If the conditions are not met, return to step 5.2, add the selected candidate solutions to the forbidden acquisition set, and continue the subsequent process; until the loop reaches the upper limit of iterations, output the current value. The corresponding fracturing construction parameters.

[0196] By selecting real geological conditions, the prediction accuracy of the multimodal spatiotemporal collaborative prediction model and the numerical simulator on the optimization results is compared, and the practicality of the model of this invention is verified.

[0197] Taking a typical horizontal well in a shale development block of Jianghan Oilfield as the research object, and selecting the core geological parameters of the reservoir as the benchmark conditions, 20 pumping schemes with systematic differences were designed. Numerical simulators and intelligent prediction models were used to calculate the reservoir stimulation volume.

[0198] Using the reservoir stimulation volume output by the numerical simulator as the "true value," the average relative error and maximum relative error are used as quantitative evaluation indicators, as shown in the following formula:

[0199] Mean relative error:

[0200] (15)

[0201] Maximum relative error:

[0202] (16)

[0203] in, For the number of schemes, For the surrogate model's predicted values, This is the true value for the numerical simulator.

[0204] The overall average relative error was 1.45%, the maximum relative error was 2.8%, and the error distribution was uniform with no obvious systematic bias. Figure 6 As shown, the prediction accuracy is high, and the prediction time for a single scheme of the prediction model (0.1s) is about 18,000 times faster than that of the numerical simulator (30min), achieving the dual advantages of "high accuracy and high efficiency". This verifies that it can be used as an effective alternative to the numerical simulator and supports the reliability assessment of subsequent real-time optimization schemes.

[0205] To comprehensively evaluate the superiority of the optimization method of this invention, three sets of comparative schemes were constructed:

[0206] Historical plans: Select the fracturing construction plans that have been implemented for this well. The plan parameters are derived from the field construction records, as shown in Table 2.

[0207] Conventional optimization scheme: The industry-standard fracturing scheme optimization method is adopted—with the goal of maximizing the reservoir stimulation volume, single-parameter traversal optimization is performed through numerical simulator to generate conventional optimization schemes, as shown in Table 3;

[0208] Intelligent optimization scheme: Based on the intelligent optimization framework constructed in step five, the optimal fracturing scheme is generated by combining the multimodal proxy model, as shown in Table 4.

[0209] Table 2 Historical Schemes;

[0210] ;

[0211] Table 3. Conventional optimization schemes;

[0212] ;

[0213] Table 4 Intelligent Optimization Scheme;

[0214] ;

[0215] To eliminate the interference of differences in computing platforms on the evaluation results, the above three schemes were uniformly imported into the multi-stage horizontal well fracturing simulation model in step 1 to calculate the fracturing effect of each scheme. The total design time of the scheme of this invention is reduced by 87.5% compared with historical schemes and by more than 95% compared with conventional optimized schemes. Figure 7 As shown, this method completely breaks through the bottleneck of time-consuming traditional methods, enabling real-time optimization of construction parameters. Furthermore, the reservoir stimulation volume reaches 430 m³, an increase of 8.86% compared to historical schemes and 3.86% compared to conventional optimization schemes, achieving improved reservoir stimulation results while shortening the design cycle.

[0216] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for intelligent real-time optimization of fracturing operation parameters assisted by deep learning, characterized in that, Includes the following steps: Step 1: Based on actual geological and engineering data, construct a multi-stage horizontal well fracturing simulation model using numerical simulation; by adjusting reservoir geological property parameters and fracturing construction parameters, conduct multiple fracturing numerical simulations to obtain the corresponding reservoir stimulation volume simulation results; Step 2: Perform unified preprocessing on multi-source data, including normalization of reservoir geological property parameters, fracturing operation parameters and reservoir stimulation volume, as well as thresholding, resampling and feature stitching of wellbore location images and natural fracture images. Step 3: Construct a multimodal spatiotemporal collaborative prediction model, which integrates a multi-dilated temporal convolutional network, a 3D residual network with spatial attention, and a bidirectional cross-attention mechanism to achieve end-to-end high-precision real-time prediction of reservoir stimulation volume. Step 4: Design the loss function and build a training strategy based on AdamW and cosine annealing. Introduce structured sparse regularization constraints to enhance the identification of key parameters, and combine adaptive optimization and learning rate scheduling to comprehensively improve prediction accuracy and generalization ability. Step 5: Integrate the outer Bayesian global search, the inner CMA-ES continuous refinement, and the discrete neighborhood switching into a collaborative optimization framework, and call the objective function to quickly evaluate the candidate solutions, so as to achieve the organic unity of global direction, local refinement and category exit, thereby obtaining a real-time optimized fracturing construction scheme that balances reservoir stimulation volume maximization and engineering workability under complex constraints. Step 3 includes the following sub-steps: Step 3.1: Based on the temporal characteristics of fracturing operation parameters, construct a multi-expansion temporal convolutional network to extract the temporal features of fracturing operation; Step 3.2: Based on the spatial structure of reservoir geological properties, construct a three-dimensional convolutional residual network to extract reservoir spatial features; Step 3.3: Introduce a bidirectional cross-attention mechanism to model the nonlinear interaction between fracturing operations and reservoir properties, and obtain bidirectional attention interaction characteristics; Step 3.4: The reservoir spatial characteristics, fracturing operation time sequence characteristics and bidirectional attention interaction characteristics are fused and stitched together to obtain a multi-source stitched tensor. The multi-source stitched tensor is compressed to a low-dimensional representation through a fully connected layer, and layer normalization and ReLU activation function processing are performed. Finally, the prediction head is used to map the representation into the final reservoir stimulation volume prediction value. Step 5 includes the following sub-steps: Step 5.1: Establish the hard constraint optimization range of fracturing construction parameters, convert the fracturing construction parameters into a hybrid vector according to the stage hard constraints, and generate multiple sets of feasible solutions by random sampling within the feasible range; Step 5.1 includes the following sub-steps: Step 5.1.1: Based on the fracturing process, the fracturing construction parameters are divided into four key stages: pre-flush fluid stage, initial stage of proppant-carrying fluid stage, late stage of proppant-carrying fluid stage, and displacement fluid stage. Hard constraints are applied to the variables of each stage to optimize the range. The variables include fracturing fluid type, proppant type, pumping rate, proppant concentration, and fluid volume. Step 5.1.2: Construct shape functions for pumping rate, sand concentration, and liquid volume for each stage. Use the starting value, ending value, and trend parameters to generate the pumping rate curve, sand concentration change curve, and liquid volume distribution curve for each stage. Then, map the curves to the hard constraint optimization range through interval mapping. Integrate the cluster spacing, stage discrete variables, and shape functions. Randomly sample within the hard constraint optimization range to obtain a batch of mixed variables H, where the stage discrete variables are fracturing fluid type and proppant type. Step 5.1.3: Expand the mixture vector H into specific parameters for 24 time steps, and perform interval pruning, sand banning, support banning, adjacent jump amplitude limiting, and consistency fine-tuning on those that do not meet the hard constraint optimization range, and record metadata. This includes the extent of repair, the number of times the limit was exceeded, the proportion of sand transported beyond the boundary, and the number of times the category was corrected; Step 5.2: Use the outer Bayesian global search to search for all feasible solutions, score all feasible solutions using the sampling function, and use sampling to select the p candidate solutions with the best scores as the output of the outer global search; In step 5.2, designing the acquisition function includes the following sub-steps: Step 5.2.1: Design the objective function and introduce a multidimensional penalty term, specifically as follows: The reservoir stimulation volume and engineering constructability are combined into a single objective function. : ；(5) in, This is the output of the multimodal spatiotemporal collaborative prediction model. The fracturing parameters are for H after repair. The penalty coefficient is... A soft penalty for the rationality of fracturing construction parameters; Four types of penalty functions are designed to suppress and regulate inappropriate behavior in the objective function: ；(6) in, The weight for adjacent transition penalties. The adjacent transition penalty function, Weights are used to penalize smoothness. Here is the smoothness penalty function. Penalty weight for sand transport capacity, Let the sand transport capacity penalty function be... Weights are used to penalize complexity. This is a complexity penalty function; Step 5.2.2: Utilize the metadata from Step 5.1.3 The feasibility model is trained using a gradient boosting decision tree, and the probability that the mixed variable H satisfies the hard constraint optimization range set in step 5.1.1 is output. Step 5.2.3: with To achieve the objective, a mixed-variable Bayesian optimization model is established, and a multimodal spatiotemporal collaborative prediction model is used to calculate the predicted mean from the output. ,variance Provides sampling and decision-making support for mixed-variable Bayesian optimization models; Step 5.2.4: Combine with desired improvements The probability output by the feasibility model Design a constraint-sensing acquisition function : ；(13) ；(14) in, It is the cumulative distribution function of the standard normal distribution; It is the probability density function of the standard normal distribution; This is the current optimal solution; To explore constants; It is a potential improvement after standardization; Step 5.3: Press The top q optimal solutions are selected from high to low, and these q optimal solutions are respectively entered into the inner local optimization. Continuous variables are adaptively shrunk and refined within the trust region using CMA-ES, and discrete variables are flipped in a small range using the neighborhood switching operator. Step 5.4: Input the q optimal solutions optimized from the inner layer into the multimodal spatiotemporal collaborative prediction model, and obtain a new global optimal solution by calculating and comparing the objective function values. Let the objective function values ​​of the new global optimal solution and the current global optimal solution be denoted as follows: and If satisfied ,in To minimize the improvement threshold, the global optimal solution is updated, and then the process returns to step 5.2 with the updated solution. Recalculate the acquisition function and repeat steps 5.2 to 5.4; if the conditions are not met, return to step 5.2, add the selected candidate solutions to the forbidden acquisition set, and continue executing step 5.3; until the loop reaches the upper limit of iterations, output the current value. The corresponding fracturing construction parameters.

2. The intelligent real-time optimization method for fracturing construction parameters assisted by deep learning as described in claim 1, characterized in that, In step 1, the reservoir geological property parameters include minimum horizontal principal stress, maximum horizontal principal stress, vertical stress, porosity, permeability, Young's modulus, and Poisson's ratio. The fracturing operation parameters include the pumping rate, fracturing fluid type, proppant type, sand concentration parameters, fluid volume, and cluster spacing configuration within the horizontal well section.

3. The intelligent real-time optimization method for fracturing construction parameters assisted by deep learning as described in claim 1, characterized in that, In step 2, the StandardScaler function is used to process reservoir geological property parameters and fracturing operation parameters, and the MinMaxScaler function is used to process reservoir stimulation volume results to eliminate the adverse effects caused by singular sample data. Both the wellbore location image and the natural fracture distribution image are converted into grayscale images. A threshold judgment method is used to judge the grayscale values ​​of the wellbore location image and the natural fracture distribution image, and grayscale value classification is performed. The wellbore location image and the natural fracture distribution image are resampled to the same size as the reservoir attribute spatial resolution. The features of the processed wellbore location image and the natural fracture distribution image are stitched with the reservoir attributes to obtain the reservoir fusion feature data, thereby expanding the feature dimension of the reservoir data volume.

4. The intelligent real-time optimization method for fracturing construction parameters assisted by deep learning as described in claim 1, characterized in that, Step 3.1 specifically involves: Step 3.1.1: Extract features from fracturing operation parameters and construct a tensor of shape [B,24,11], where B is the batch size, 24 is the time step, and 11 is the fracturing operation feature dimension; Step 3.1.2: Map the 11-dimensional features to a 32-dimensional high-dimensional space through a fully connected layer; then, perform normalization using LayerNorm and apply the ReLU activation function for non-linear activation; Step 3.1.3: Design TCNBlock. Each TCNBlock uses Conv1d to perform one-dimensional convolution operations with a kernel size of 3, and applies dilated convolution to expand the receptive field. Each convolutional layer is processed by batch normalization, ReLU activation function, Dropout, and residual connections. Step 3.1.4: Stack 3 TCNBlocks, with expansion rates of 1, 2, and 4 respectively, to gradually increase the receptive field in order to capture dependencies at different time scales; Step 3.1.5: After processing by 3 TCNBlocks, the time-series feature tensor is output with shape [B,32,24]. Mean pooling and max pooling are used to extract the global features of the time-series feature tensor, and the two pooling results are concatenated into a feature vector. Finally, the fully connected layer is used for further processing to obtain the fracturing construction time-series feature vector with shape [B,32]. Step 3.2 specifically involves: Step 3.2.1: Extract features from reservoir fusion feature data and construct a tensor of shape [B,9,5,239,76], where 9 is the feature dimension, 5 is the spatial dimension in the Z direction, and 239 and 76 represent the spatial dimensions in the X and Y directions; Step 3.2.2: Introduce a 3D convolutional layer with a kernel size of 3×3×3 and padding=1 to map the input 9-channel features to a 16-channel output, resulting in a feature tensor shape of [B,16,5,239,76]. Step 3.2.3: Construct a two-layer 3D residual block. Each residual block contains two layers of 3D convolutions with a kernel size of 3×3×3. Batch normalization and GELU activation functions are added after each layer. The residual block adopts a residual connection structure with identity mapping. When the number of input and output channels is inconsistent, channel matching is achieved through 1×1×1 convolution to ensure consistent dimensions. The first residual block expands the number of channels from 16 to 32, and the second residual block further expands it to 64. The final output feature tensor shape is [B,64,5,239,76]. Step 3.2.4: Introduce a spatial attention mechanism. Calculate the global context of the feature tensor output by the two-layer 3D residual block through local average pooling. Then, use 1×1×1 convolution and the Sigmoid function to generate a spatial attention weight map. Weight the input features point by point to obtain the weighted spatial feature tensor with the shape [B,64,5,239,76]. Step 3.2.5: Employing a multi-scale feature aggregation strategy, global average pooling is applied to the weighted spatial feature tensor to reduce the spatial dimension to a single global feature vector, with an output shape of [B, 64]. Spatial pooling is then applied to the weighted spatial feature tensor to retain some spatial information, resulting in an output shape of [B, 1152]. The global pooling feature and the spatial pooling feature are then concatenated to obtain a fused feature of [B, 1216]. After processing through a fully connected layer, layer normalization is performed to compress the fused feature dimension to [B, 128], outputting the reservoir spatial feature vector. Step 3.3 specifically involves: Step 3.3.1: Map the fracturing operation timing characteristics and reservoir spatial characteristics to a unified 64-dimensional space through linear projection; Step 3.3.2: Feature interaction is performed through a multi-head attention mechanism, which is divided into two paths. The first path uses reservoir spatial features as the query and fracturing operation time sequence features as the key and value, which is used to extract response features related to key areas of the reservoir from the fracturing operation dynamics. The second approach uses fracturing operation characteristics as the query and reservoir characteristics as the key and value to capture the intensity of the impact of different fracturing operation stages on the reservoir space. Step 3.3.3: Employ a multi-head attention mechanism to combine different subspaces. For each query vector in Step 3.3.2, multiple attention heads will learn different attention weights and perform a weighted average within each head to obtain two 64-dimensional output features. Step 3.3.4: The two 64-dimensional output features are respectively processed by linear transformation, GELU activation function, and Dropout, and finally passed through a linear layer to output two 64-dimensional interaction feature vectors; Step 3.3.5: Concatenate the two interactive feature vectors obtained in step 3.3.4 to obtain a 128-dimensional feature vector with shape [B, 128].

5. The intelligent real-time optimization method for fracturing construction parameters assisted by deep learning as described in claim 1, characterized in that, Step 4 includes the following sub-steps: Step 4.1: Use mean squared error as the main loss function, and add L1 sparse regularization to the feature projection layer of the encoder of construction parameters to ensure that the model can automatically identify and focus on key construction parameters. The expression is: ；(3) in, Mean square error, Indicates the first One sample, To modify the predicted volume, To modify the actual volume value, The number of samples; ；(4) in, For sparse regularization terms, For the model One parameter, The regularization coefficient is used. Step 4.2: Using the AdamW optimizer and cosine annealing learning rate scheduling mechanism, the learning rates of the fracturing operation time encoder and the reservoir space encoder are set differently for model training.

6. The intelligent real-time optimization method for fracturing construction parameters assisted and accelerated by deep learning according to claim 1, characterized in that, In step 5.2.1, an adjacent jump penalty is designed to limit the rate and sand concentration variation between adjacent time steps: ；(7) in, Index for the stage; Representative stage A single time step index within; For the stage The set of all time step indices in the middle; This represents the difference in pumping rate between adjacent time steps. For the stage Threshold for maximum pumping rate variation; This represents the difference in sand concentration between adjacent time steps. For the stage The threshold for the maximum sand concentration change in the middle; Representation phase Does it contain sand? Indicates the addition of sand. This indicates that no sand is added; Design a smoothness penalty to limit the rate and frequency of sand concentration changes between adjacent time steps: ；(8) in, For the stage The maximum range of change in medium speed; For the stage Maximum variation range of medium sand concentration; For the stage The maximum range of variation in liquid volume; The weighting coefficients are used to adjust the total variation and the second difference weights. The design imposes a penalty on the sand conveying capacity, limiting the instantaneous sand conveying rate to no more than the equipment's load capacity. ；(9) in, For the stage The upper limit of medium sand conveying capacity indicates the maximum sand conveying rate that the equipment or wellbore can withstand. For sand concentration, This refers to the pumping rate; Design complexity penalty to avoid overly complex fracturing construction schemes: ；(10) ；(11) ；(12) in, For sand concentration changes Sign judgment, The threshold for sand concentration change, when hour, ,when hour, ,when hour, ; For determining the sign of the change in pumping rate, The threshold for the change in pumping rate, when hour, ,when and hour, ,when and hour, ; For the stage The number of peak concentrations of medium sand; For the stage The number of inflection points at medium speed; For the stage Maximum number of peak values ​​allowed in the middle; For the stage The maximum number of inflection points allowed in the middle.

7. The intelligent real-time optimization method for fracturing construction parameters assisted by deep learning as described in claim 1, characterized in that, Step 5.3 includes the following sub-steps: Step 5.3.1: From the p candidate solutions obtained in the outer optimization, sort them according to... Select the top q optimal solutions from highest to lowest; Step 5.3.2: Expand the q optimal solutions in the outer optimization into specific time step parameters, and use the Sigmoid function to automatically repair abnormal parameters; Step 5.3.3: Use the CMA-ES algorithm to perform local optimization on the continuous variables of the q optimal solutions, and adjust the search strategy through the adaptive trust region algorithm, where the initial step size is set to 0.3, the trust region radius is set to 20% of the interval width, and the population size is set to 5. Step 5.3.4: Use the neighborhood switching operator to flip the local class of the discrete variable.

Images