A Multi-Well Collaborative Productivity Prediction Method and System Based on Inter-well Interference Coupling

Abstract

This invention relates to the field of multi-well productivity prediction technology, specifically to a method and system for multi-well collaborative productivity prediction based on inter-well interference coupling. The method includes: acquiring characteristic parameters, measured data, and geological parameters of each well; determining cross-well coupling relationships based on the characteristic parameters of each well; repairing outliers in the measured data of each well to obtain a multi-well dataset; obtaining a multi-well spatiotemporal feature matrix based on the cross-well coupling relationships and the multi-well dataset; slicing the multi-well spatiotemporal feature matrix; fusing high-frequency and low-frequency interference based on the slicing results to obtain multi-scale spatiotemporal features; obtaining an adaptive prior distribution based on the multi-scale spatiotemporal features; iteratively obtaining a posterior distribution by applying variational inference to the geological parameters and the adaptive prior distribution; interactively correcting the posterior distribution and the multi-scale spatiotemporal features; and iteratively updating the interference radius and inter-well connectivity, thereby achieving multi-well productivity prediction. This invention improves the accuracy of multi-well productivity prediction results.

Description

Technical Field

[0001] This invention relates to the field of multi-well productivity prediction technology, specifically to a multi-well collaborative productivity prediction method and system based on inter-well interference coupling. Background Technology

[0002] Currently, existing technologies for inter-well collaborative production prediction typically employ time series analysis and multi-source data fusion to improve the accuracy of prediction models. These methods correlate surface monitoring data, dynamic observations of production wells, and surrounding geological information. Some schemes dynamically update fracture distribution and pressure diffusion based on the spatiotemporal stratification characteristics of geology or microseismic data. Other schemes construct bottom-level sensitivity factors using multi-statistical and geological methods, and achieve accurate prediction of inter-well stress and flow boundaries by comparing the pressure fields of fractured and unfractured zones. While existing shale oil reservoir production prediction methods can quickly predict shale oil reservoir production, in practical applications, they fail to fully consider the differences in formation environment and production conditions of each well, and the dynamic evolution of inter-well interference makes it impossible to derive a unified and accurate inter-well connectivity relationship and coupling law, resulting in low accuracy of the final production prediction results. Summary of the Invention

[0003] To address the issue of low accuracy in predicting multi-well productivity using existing methods, this invention aims to provide a multi-well collaborative productivity prediction method and system based on inter-well interference coupling. The specific technical solution adopted is as follows: In a first aspect, the present invention provides a multi-well collaborative production capacity prediction method based on inter-well interference coupling, the method comprising the following steps: The characteristic parameters, measured data, and geological parameters of each well are obtained. The characteristic parameters include coordinates and fracturing morphology. The measured data include pressure and production. The geological parameters include well connectivity, permeability, fluid viscosity, fluid compressibility coefficient, formation pressure correlation coefficient, and formation compressibility factor. Based on the characteristic parameters of each well, the cross-well coupling relationship is determined, and the outliers in the measured data of each well are repaired to obtain a multi-well dataset; based on the cross-well coupling relationship and the multi-well dataset, a multi-well spatiotemporal feature matrix is ​​obtained. The multi-well spatiotemporal feature matrix is ​​sliced, and high-frequency interference and low-frequency interference are fused according to the slicing results to obtain multi-scale spatiotemporal features; An adaptive prior distribution is obtained based on the multi-scale spatiotemporal features. A posterior distribution is obtained by iterating the geological parameters and the adaptive prior distribution using variational inference. The posterior distribution and the multi-scale spatiotemporal features are interactively corrected to obtain the corrected multi-scale spatiotemporal features, the iteratively updated interference radius, and the inter-well connectivity. Based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity, the collaborative production capacity of multiple wells is predicted.

[0004] Preferably, determining the cross-well coupling relationship based on the characteristic parameters of each well includes: The relationship between the two wells is quantified by using their relative positions, as well as the permeability and fluid viscosity of the block in which they are located. A cross-well coupling matrix is ​​constructed based on the relationships between every pair of wells. The cross-well coupling matrix is ​​used to reflect the cross-well coupling relationships.

[0005] Preferably, the step of repairing outliers in the measured data of each well to obtain a multi-well dataset includes: For any kind of measured data: The cumulative difference of each well is quantified by combining the well-to-well relationships of each well with the first difference between any measured data of the other wells and the corresponding reference estimate. The differential accumulation is used to identify and remove outliers in any type of measured data from each well, and the inter-well relationship is used to repair missing values ​​in the measured data of the wells. Based on the repaired data, the time axis of each well is standardized by using the ratio between the fracture propagation rate of each well and other wells during the fracturing operation period and the fracturing time of other wells, to obtain the multi-well dataset containing cross-well physical priors.

[0006] Preferably, obtaining the multi-well spatiotemporal feature matrix based on the cross-well coupling relationship and the multi-well dataset includes: Based on the feature parameters of each well and the minimum distance between the well and the grid node within a single time segment, the feature parameters of the well are projected onto the corresponding grid node corresponding to the minimum distance using an inverse distance weighting algorithm, thereby obtaining the feature value of each grid node in a single time segment; the feature values ​​of all grid nodes constitute the initial multi-well spatiotemporal feature matrix. The adjacent well interference radius of each well is determined based on the cross-well coupling relationship; the adaptive correction weight of the grid node is determined based on the influence of the flow boundary pressure, the distance between the well and the grid node on the grid node, the relationship between the distance between the well and the grid node and the corresponding adjacent well interference radius. Based on permeability, porosity, fluid compressibility coefficient, and adaptive correction weights, the eigenvalues ​​of the grid nodes in the initial multi-well spatiotemporal feature matrix are iterated to obtain the multi-well spatiotemporal feature matrix after updating the grid node coupling degree.

[0007] Preferably, the step of slicing the multi-well spatiotemporal feature matrix and fusing high-frequency and low-frequency interferences based on the slicing results to obtain multi-scale spatiotemporal features includes: Various slice combinations are constructed based on the radius of the well and the length of the time window; Based on the multi-well spatiotemporal target feature matrix, the overall distance distribution between wells and grid nodes at all acquisition times within a single time segment, and the relationship between the radius of wells in various slice combinations, the aggregate value of various slice combinations is obtained. Using a local interpolation operator, high-frequency interference generated by inter-well mutual interference under various slice combinations is identified based on the aggregated value; The aggregation value and pressure propagation characteristics are fused using a fusion operator to obtain low-frequency interference under various slice combinations. The pressure propagation characteristics are determined based on permeability and porosity. By fusing the high-frequency interference and the low-frequency interference, multi-scale spatiotemporal features under various slice combinations are obtained.

[0008] Preferably, obtaining the adaptive prior distribution based on the multi-scale spatiotemporal features includes: A priori aggregation function is obtained based on the aforementioned multi-scale spatiotemporal features; Based on the reservoir pressure transmission model and fracture propagation mechanism, the prior aggregation function is transformed into a normalized probability form to obtain an adaptive prior distribution.

[0009] Preferably, the step of using variational inference to iterate the geological parameters and the adaptive prior distribution to obtain the posterior distribution includes: The variational inference objective function is derived based on geological parameters and adaptive prior distribution; By maximizing the value of the variational inference objective function, the initial posterior distribution under the constraints of geological parameters and adaptive prior distribution is obtained. The geological parameters in the initial posterior distribution are projected onto the adaptive prior distribution, and the posterior distribution is obtained by a second iteration correction.

[0010] Preferably, the step of interactively correcting the posterior distribution and multi-scale spatiotemporal features to obtain the corrected multi-scale spatiotemporal features, the iteratively updated interference radius, and the inter-well connectivity includes: The connectivity and permeability in the posterior distribution are mapped to the multi-scale spatiotemporal features. An error metric function is used to identify regions where there are errors in the connectivity and permeability between wells. Based on the connectivity and permeability in the posterior distribution, the corrected interference weight of each grid node in the region where the error is greater than the threshold is determined. The modified interference weights are used to weight the eigenvalues ​​of the grid nodes in the multi-well spatiotemporal feature matrix to obtain the modified multi-well spatiotemporal feature matrix; the modified multi-well spatiotemporal feature matrix is ​​then re-sliced ​​and aggregated to obtain the modified multi-scale spatiotemporal features. The interference radius and connectivity are iteratively updated based on the posterior distribution to obtain the corrected permeability, porosity and fluid compressibility.

[0011] Preferably, the prediction of multi-well collaborative production capacity based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity includes: Based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, inter-well connectivity, and reservoir pressure transmission model, the inter-well mutual interference coupling equation is obtained, and the production intensity of each well under the combined effect of inter-well connectivity and production is predicted. When extreme operational events occur, a short-term reset function is used to adjust the inter-well interference weights, correct the calculation of production intensity in the inter-well mutual interference coupling equation, and obtain the prediction results; the extreme operational events include shutting in wells and opening new wells; The predicted results are compared with the posterior distribution, and the model output is revised using a dynamic correction objective function until the preset conditions are met, thus obtaining the final predicted result of multi-well collaborative production capacity.

[0012] Secondly, the present invention also provides a multi-well collaborative production capacity prediction system based on inter-well interference coupling, which is used to implement the above-mentioned method. The system includes: The data acquisition module is used to acquire characteristic parameters, measured data and geological parameters of each well. The characteristic parameters include coordinates and fracturing morphology. The measured data includes pressure and production. The geological parameters include inter-well connectivity, permeability, fluid viscosity, fluid compressibility coefficient, formation pressure correlation coefficient and formation compressibility factor. The first feature acquisition module is used to determine the cross-well coupling relationship based on the feature parameters of each well, repair outliers in the measured data of each well to obtain a multi-well dataset, and obtain a multi-well spatiotemporal feature matrix based on the cross-well coupling relationship and the multi-well dataset. The second feature acquisition module is used to slice the multi-well spatiotemporal feature matrix and fuse high-frequency interference and low-frequency interference according to the slicing results to obtain multi-scale spatiotemporal features. The correction module is used to obtain an adaptive prior distribution based on the multi-scale spatiotemporal features, and to obtain a posterior distribution by iteratively applying variational inference to the geological parameters and the adaptive prior distribution; the posterior distribution and the multi-scale spatiotemporal features are interactively corrected to obtain the corrected multi-scale spatiotemporal features, the iteratively updated interference radius, and the inter-well connectivity. The production capacity prediction module is used to predict the collaborative production capacity of multiple wells based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity.

[0013] The present invention has at least the following beneficial effects: This invention first evaluates the cross-well coupling relationship based on the characteristic parameters of multiple wells, and constructs a multi-well dataset by correcting outliers in the actual collected pressure and production data. The construction of the multi-well dataset not only considers the data quality of the measured data itself, but also takes into account the impact of the relationship between the wells on the recovery of outliers, so as to better present the real pressure and production data of various geological blocks. Then, the multi-well spatiotemporal feature matrix is ​​generated by combining the cross-well coupling relationship and the multi-well dataset. Based on the slicing results of the multi-well spatiotemporal feature matrix, high-frequency interference and low-frequency interference are evaluated respectively, capturing the mutual interference signal and pressure propagation effect between the wells, and obtaining multi-scale spatiotemporal features. Furthermore, the prior information of the multi-well geological parameters and cross-well coupling is integrated to iterate the prior distribution to obtain the posterior distribution. On this basis, the multi-scale spatiotemporal features, interference radius and inter-well connectivity are corrected, so that the well parameters can be analyzed and corrected in real time at multiple scales. This ensures that the linkage relationship of the multi-wells can be accurately identified and continuously updated under complex geological conditions and variable well control strategies, thereby realizing the accurate prediction of the multi-well collaborative production capacity and improving the accuracy of the multi-well collaborative production capacity prediction results. Attached Figure Description

[0014] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0015] Figure 1 A flowchart of a multi-well collaborative production capacity prediction method based on inter-well interference coupling provided in an embodiment of the present invention; Figure 2 This is a structural block diagram of a multi-well collaborative production capacity prediction system based on inter-well interference coupling provided in an embodiment of the present invention. Detailed Implementation

[0016] To further illustrate the technical means and effects adopted by the present invention to achieve the intended purpose, the following detailed description of a multi-well collaborative production capacity prediction method and system based on inter-well interference coupling proposed in accordance with the present invention is provided in conjunction with the accompanying drawings and preferred embodiments.

[0017] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0018] The following description, in conjunction with the accompanying drawings, details the specific scheme of a multi-well collaborative production capacity prediction method and system based on inter-well interference coupling provided by the present invention.

[0019] Example of a multi-well collaborative productivity prediction method based on inter-well interference coupling: The specific scenario addressed in this embodiment is as follows: When predicting the production capacity of multiple wells, the production capacity is often predicted based solely on the historical data of each well. However, the accuracy of this method is relatively low. In order to improve the accuracy of the prediction results, this embodiment will combine the relationship between different wells, the formation environment of each well, and the production conditions to conduct a comprehensive analysis and obtain more accurate production capacity prediction results.

[0020] This embodiment proposes a multi-well collaborative production capacity prediction method based on inter-well interference coupling, such as... Figure 1 As shown, the multi-well collaborative production capacity prediction method based on inter-well interference coupling in this embodiment includes the following steps: Step S1: Obtain characteristic parameters, measured data, and geological parameters for each well. The characteristic parameters include coordinates and fracturing morphology. The measured data includes pressure and production rate. The geological parameters include well connectivity, permeability, fluid viscosity, fluid compressibility coefficient, formation pressure correlation coefficient, and formation compressibility factor.

[0021] First, characteristic parameters and geological parameters of each well during historical monitoring are retrieved from the database. Characteristic parameters include well coordinates, fracturing morphology, and the block number to which the well belongs. The size of the block is determined by the implementer based on specific circumstances. Fracturing morphology includes crack propagation direction, crack length, and crack width. It should be noted that since cracks may be irregular and their width may vary at different locations, all locations within the crack area are obtained, and a least-squares method is used to fit a straight line to these locations. The direction of the fitted line is taken as the crack propagation direction. The average of all crack width data is taken as the crack width. Geological parameters include well connectivity, permeability, fluid viscosity, fluid compressibility coefficient, formation pressure correlation coefficient, and formation compressibility factor.

[0022] The fluid compressibility coefficient is determined based on the type of fluid; the fluid compressibility coefficient of water is... The fluid compressibility coefficient of crude oil is The fluid compressibility coefficient of natural gas is The compressibility coefficient of water-based drilling fluid is The compressibility coefficient of oil-based drilling fluid is .

[0023] Simultaneously, the pressure and production of each well are collected in real time. It should also be noted that: the production figures mentioned later in this embodiment refer to the total production of the wells within the time period obtained with each time as the last time point. For example, the production at the current time point is actually the total production of the most recent day.

[0024] Step S2: Determine the cross-well coupling relationship based on the characteristic parameters of each well, repair the outliers in the measured data of each well to obtain a multi-well dataset; obtain the multi-well spatiotemporal feature matrix based on the cross-well coupling relationship and the multi-well dataset.

[0025] A cross-well coupling matrix is ​​constructed using the characteristic parameters of each well. The cross-well coupling relationships are then analyzed using this matrix. The specific cross-well coupling matrix is... It can be represented as:

[0026]

[0027] in, Indicates the first Koujing and the first The relationship between wells. Indicates the first The relationship between the wellhead and its own internal wells. This indicates the relationship between the first and second wells. This indicates the relationship between well number 1 and well number M. This indicates the relationship between the second well and the first well. This indicates the relationship between the second well and the second well. This indicates the relationship between the second well and the Mth well. Indicates the first The relationship between well 1 and well 2. Indicates the first The relationship between the first well and the second well. Indicates the first The relationship between the wellhead and its own internal wells. Indicates the number of wells; This represents the coordinates of the i-th well in the geographic coordinate system. This represents the coordinates of the j-th well in the geographic coordinate system. This represents the influence coefficient of the well-to-well distance. Indicates the first Koujing and the first A quantitative indicator of the length of the wellhead in the overlapping area of ​​the fracturing morphology. Indicates crossing the first Koujing and the first The flow path coupling strength in the overlapping region of fracturing morphology in wellheads. This represents an exponential function with the natural constant as its base.

[0028] The well-to-well distance influence coefficient is used to measure the impact of well-to-well distance. Specifically, it is calculated using the permeability and fluid viscosity of the block containing the well. If well i and well j belong to the same block, the ratio of permeability to fluid viscosity corresponding to that block is used as the well-to-well distance influence coefficient. If well i and well j are not in the same block, the ratios of permeability to fluid viscosity corresponding to the block containing well i and well j are calculated separately, and then the two ratios are weighted and summed. This weighted sum is used as the well-to-well distance influence coefficient. In the weighting process, the weights for the two ratios are the normalized areas of the corresponding blocks. This embodiment uses a maximum-minimum value normalization method for area normalization. However, other existing data normalization methods can also be used in other implementations.

[0029] Indicates the first Koujing and the first The coordinate difference norm between wells measures the degree of difference between two coordinate systems; a larger norm indicates a greater degree of difference between the two wells. As one implementation, the coordinate difference norm between wells is represented by the L2 norm, which is the... The coordinates of the well and the first The L2 norm (Euclidean distance) of the well's coordinates is used as the first... Koujing and the first The coordinate difference norm between wells; as another implementation, the coordinate difference norm between wells can also be represented by calculating the L1 norm, that is, the L1 norm. The coordinates of the well and the first The L1 norm (Manhattan distance) of the well's coordinates is used as the first... Koujing and the first The coordinate difference norm between wells.

[0030] No. Koujing and the first The overlapping area of ​​fracturing morphology in a well is the overlapping area of ​​fractures. The length quantification index is calculated using a fracturing morphology vector. Specifically, a Cartesian coordinate system is constructed with the fracture center as the origin and along the fracture extension direction. The length of the first fracture is plotted in this coordinate system. Koujing and the first The smallest polygon of the fracture in the well is used. Boolean operations on the polygons are used to obtain the overlapping region of the fracture polygons in the two wells. The length of the overlapping region is measured along the line connecting the two wells, and this length is taken as the smallest polygon of the fracture in the well. Koujing and the first The length quantification index of the well in the overlapping area of ​​fracturing morphology. The fracturing morphology vector is composed of the crack extension direction, crack length, and crack width.

[0031] Flow path coupling strength Calculations are performed using the permeability of the overlapping region, the average crack width of the overlapping region, and the fluid viscosity. Specifically, the product of the permeability and the average crack width of the overlapping region is first calculated, and then the ratio of this product to the fluid viscosity of the overlapping region is calculated. This ratio is then multiplied by the initial flow path coupling strength value as the value used to cross the first... Koujing and the first The flow path coupling strength in the overlapping area of ​​the fracturing morphology of the well.

[0032] Meanwhile, the well's coordinate vector, fracturing morphology vector, and block number to which the well belongs are mapped to each element in the cross-well coupling matrix to avoid the problem of well disconnection and distortion caused by differences in geological blocks.

[0033] Furthermore, the deviation is analyzed by examining the differences between single-well data and reference data, as well as the cross-well coupling relationship. The cumulative difference is calculated, and then used to identify and remove outliers and jump points in the historical data (measured data) of each well. Multi-well coupling interpolation is used to repair missing values ​​in the historical data of each well. The historical data of each well includes both production and pressure data. This embodiment will use one type of measured data as an example for explanation; the method provided in this embodiment can be used to process other types of measured data. The cumulative difference can be expressed as:

[0034] in, Indicates the first The cumulative difference of wells is used to summarize the inter-well influences for the first well. Indicators for determining the difference value of wells; Indicates the number of wells; Indicates the first Koujing and the first The relationship between wells; Indicates the first Koujing in time Actual measured data (production or pressure); Indicates the first Koujing in time The reference estimate of the measured data, This indicates the absolute value sign.

[0035] No. Koujing in time The reference estimates of the measured data are obtained based on historical characteristics or conventional trend model predictions. As a specific example, the ARIMA model is used to obtain the trend forecast values ​​of output or pressure as reference estimates.

[0036] Based on the above calculations, the difference at all times is accumulated. The value is used as a preset threshold, where Let the standard deviation of the cumulative difference be , if If the preset threshold is exceeded, it indicates that the first... Koujing in time If any jumps occur that do not conform to the physical correlation across wells, the corresponding data should be deleted directly.

[0037] Furthermore, considering the inter-well coupling strength and formation compressibility to construct interpolation weights, multi-well coupling interpolation is calculated based on these weights. The specific calculation formula is as follows:

[0038]

[0039] in, Indicates the first Koujing in time Interpolated data at the location; Indicates the first Koujing in time Available measured data; This represents the combined weighting of cross-well physical correlation and formation compressibility factor during multi-well coupled interpolation. Indicates the formation compressibility factor; Indicates the first Koujing and the first The relationship between wells.

[0040] It should be noted that the available test data is the data that has not been deleted.

[0041] The above method can obtain the interpolated data of each rejection location. The interpolated data is then used for interpolation. The processed data can meet the requirements of cross-well physical linkage and take into account the formation differences between wells.

[0042] Furthermore, the time axis of each well is standardized by using the fracture propagation rate during the fracturing operation period. Specifically, the fracturing fluid injection rate and fracture extension pressure are obtained based on the fracturing operation log, and the fracture propagation rate is calculated by the PKN model. Taking the fracturing start time of the first well as the benchmark, the fracturing time of other wells is calculated as the standardized time according to the fracture propagation rate. That is, the ratio of the fracture propagation rate between the benchmark well and each of the other wells is first calculated, and the product of this ratio and the fracturing time of each of the other wells is used as the standardized time of each of the other wells.

[0043] Aligning cross-well events involves combining the cross-well coupling matrix with the interpolation and correction results obtained in the previous steps to check the interaction consistency between multiple wells, remove data that does not match the fracture propagation rate, and obtain a high-confidence multi-well dataset containing cross-well physical priors.

[0044] A unified reading of the high-confidence multi-well dataset is performed, and the output includes the production and pressure of each well at different times. The time parameters of each well are denoted as... Production is recorded as Pressure is recorded as .

[0045] Furthermore, based on the geographical coordinates of each well, a unified two-dimensional grid coordinate system is constructed. The density and size of the grid nodes are set according to the formation division results, ensuring that each grid node matches the well location distribution. As a specific implementation method, grid nodes are set in the following way: using the circumscribed rectangle of the well cluster distribution area as the boundary, the number of nodes in the horizontal and vertical directions of the grid are set. Specifically, the minimum circumscribed rectangle of the well cluster location is constructed, and grid nodes are set based on this minimum circumscribed rectangle. When setting grid nodes, the following condition must be met: half the distance between adjacent grid nodes is greater than the distance between each well and its nearest grid node. The grid nodes must meet the above conditions during the setting process. The specific setting method and the number of nodes will be determined by the implementer based on the specific circumstances.

[0046] Furthermore, the production and pressure of each well are divided into time segments and projected onto a two-dimensional grid coordinate system to obtain a preliminary multi-well spatiotemporal feature matrix. Specifically, the standardized time is first divided into multiple time segments according to the production cycle, and the duration of each time segment is set by the implementer based on specific circumstances. Within each time segment, the average production and average pressure of a single well are calculated, and the K-nearest neighbor algorithm is used to obtain the nearest multiple wells to each grid node. Combined with the inverse distance weighted algorithm, the well data is projected onto the grid nodes to obtain the feature values ​​of each grid node in different time segments. The specific formula for calculating the feature values ​​is as follows:

[0047] in, Indicates time segment Time grid nodes eigenvalues ​​at that location This indicates that the K-nearest neighbor algorithm is used to obtain the distance to the grid node. The number of the nearest wells; and These represent time segments. Time The mean production and mean pressure of each well. The normalization coefficient representing the pressure-output ratio. Represents grid nodes The distance between the uth well and the well.

[0048] The formula for calculating the normalization coefficient of pressure-output is as follows: ,in and These represent the maximum production and maximum pressure of a single well, respectively. The normalization coefficient between pressure and production is introduced in the formula for calculating the eigenvalue to avoid the influence of dimensional differences on the projection results.

[0049] The purpose of obtaining the feature values ​​of different grid nodes by the above method is to map output and pressure to spatiotemporal nodes according to the division of time segments. In this embodiment, the inverse distance weighting method is used to map the spatiotemporal nodes.

[0050] Furthermore, based on the above projection results, the initial multi-well spatiotemporal feature matrix is ​​obtained as follows:

[0051] in, This represents the initial multi-well spatiotemporal feature matrix; This indicates the number of grid nodes in the horizontal direction, i.e., the number of columns of grid nodes; This indicates the number of grid nodes in the vertical direction, i.e., the number of rows of grid nodes; Indicates the number of a time segment; Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes eigenvalues ​​at that location Indicates time segment Time grid nodes The eigenvalue at that location.

[0052] In this embodiment, the feature values ​​of the grid nodes are constructed by fusing the distance relationship between wells and the multi-well dataset. The contribution of local wells to the grid nodes is extracted from the integrated production and pressure, thereby obtaining the initial multi-well spatiotemporal feature matrix.

[0053] Next, the interference radius and flow boundary of each well are determined, and the grid nodes are adaptively corrected based on the interference radius and flow boundary of the adjacent wells.

[0054] Specifically, all other wells with a well-to-well relationship value greater than a preset relationship threshold are obtained from the cross-well coupling matrix, and the maximum distance between the current well and the other wells is obtained as the neighboring well interference radius of the current well. Considering that there is a significant interference effect between neighboring wells when the coupling strength is greater than 10%, the preset relationship threshold is set to 0.1.

[0055] Based on the flow boundary of the block where each well is located, the difference between the initial pressure of the grid node and the pressure of its nearest flow boundary is calculated. The ratio of this difference to the average pressure gradient of the block where the grid node is located is used as the characteristic distance between the grid node and the flow boundary. Then, the difference between the average bottomhole pressure of each well and the corresponding fluid boundary pressure is calculated. The ratio of this difference to the average pressure gradient of the corresponding block is used as the flow boundary threshold for each well.

[0056] As a specific example, the average pressure gradient of the block where the grid node is located can be obtained by the following method: taking each well in the block where the grid node is located as the center, calculate the difference between the average bottom hole pressure of each well and the average pressure of each other well in the block. Calculate the ratio of this difference to the straight-line distance between each well and each other well in the block where the grid node is located. Take the average of all these ratios as the average pressure gradient of the block where the grid node is located.

[0057] The initial pressure of a grid node can be obtained in the following way: First, the original formation pressure of each block is extracted from the geological database. If the grid node is not located on the edge of the block, the initial pressure of the grid node is the original formation pressure of the block in which it is located. If the grid node is located on the edge of the block, the initial pressure of the grid node is obtained by interpolation of the original formation pressure of the block that shares the edge point. The specific interpolation method is Kriging interpolation.

[0058] Next, the pressure on the grid nodes will be iteratively updated by calculating the adaptive weights of the grid nodes. The formula for calculating the adaptive weights of the grid nodes is:

[0059] in, Represents grid nodes Adaptive weights; This indicates that the K-nearest neighbor algorithm is used to obtain the distance to the grid node. The number of the nearest wells; Indicates the relationship with the first Wellhead-related coupling factors; Indicates the distance attenuation coefficient; Represents grid nodes With the The distance between wells; Indicates the first The interference radius of adjacent wells of the well; Represents grid nodes Characteristic distance from the flow boundary; Indicates the first The flow boundary threshold of the wellhead; This represents the step function.

[0060] The coupling factor associated with each well is determined by the average production of the block in which each well is located and the inter-well coupling relationship. Specifically, firstly, the average production of each well is calculated, and the maximum value among all wells with inter-well relationships to the current well is obtained. The product of the average production and this maximum value is calculated, and this product is recorded as the first product. The coupling factor with the grid nodes is then calculated. The ratio between the first product corresponding to each nearest well and the sum of the first products corresponding to all wells is used as the coupling factor associated with each well. That is, the coupling factor is the result of normalization based on well production intensity.

[0061] The distance decay coefficient is determined by the heterogeneity of the reservoir. The specific method for obtaining it is as follows: First, obtain the average permeability and average fluid viscosity in the block where the well group is located. Use the ratio of average fluid viscosity to average permeability as the distance decay coefficient. The distance decay coefficient corresponding to high permeability reservoirs is smaller and decays more slowly.

[0062] Furthermore, based on the permeability tensor, porosity, and fluid compressibility coefficient, the grid nodes are iteratively analyzed to obtain an updated node coupling degree, thereby refining the initial multi-well spatiotemporal characteristic matrix. The formula for calculating the node coupling degree is:

[0063] in, Represents grid nodes The degree of node coupling; Indicates at grid nodes At any time The pressure; Indicates porosity; Indicates the fluid compressibility coefficient; Represents grid nodes The number of wells within the interference radius of the nearest well; This represents the value of the source term function when the grid node coincides with the well coordinates. The source term function is the Dirac function, meaning that the value of the source term function is 1 when the grid node is at the well location, and 0 when the grid node is not at the well location. The purpose of calculating the source term function is to determine whether the grid node is at the well location. Indicates the first The flow contribution of the wellhead to the node; and This represents the coordinate position of the nth well.

[0064] The flow contribution is determined by calculating the production of each well at each time step, along with the formation thickness and well radius. The relationship between these factors is as follows: ,in For the first The production of the well at time t. and These represent the well radius and the formation thickness of the block, respectively.

[0065] After obtaining the node coupling degree of the grid nodes based on the above calculations, the elements in the initial multi-well spatiotemporal feature matrix are corrected using the finite difference method. Specifically, the weights of the grid nodes are first... Using the initial conditions, the stress at the nodes is calculated iteratively using the finite difference method. The node coupling degree is updated using the calculation formula, and the updated node coupling degree is used to correct the elements in the initial multi-well spatiotemporal feature matrix. This is achieved through calculation... right Make corrections until the difference in coupling between adjacent iteration nodes is less than 1% and convergence is achieved.

[0066] When using the finite difference method for iteration, the input data includes static basic data, geological and fluid parameter data, production dynamic data, and initial pressure of the grid nodes. The static basic data includes grid node coordinates, grid spacing, and grid node weights. The geological and fluid parameter data refers to permeability, porosity, fluid viscosity, fluid compressibility, and formation pressure correlation coefficients obtained by querying the block number. The production dynamic data includes measured production rates and bottom hole pressures for single wells. The specific calculation steps include: (1) Constructing seepage control equations: Based on Darcy's law and the law of conservation of mass, establish pressure control equations including node weights to quantify the relationship between pressure, permeability, and yield; (2) Finite difference discretization: Discretize the continuous reservoir pressure field into algebraic equations for each grid node, approximate the pressure gradient with the pressure difference between adjacent grid nodes, and substitute it into the governing equations; (3) Substitute input parameters: Substitute static basic data, geological and fluid parameter data, and production dynamic data into the discrete equations to form a linear equation system; (4) Iterative solution: The conjugate gradient method is used to solve the system of equations to obtain the grid node pressure in this iteration; (5) Convergence judgment: Calculate the maximum relative error between the pressure of the current round and the pressure of the previous round. If the error is ≤1%, the convergence is achieved and the final pressure of the grid node is output. If the convergence is not achieved, the pressure of the current round is used as the initial value of the next round, and the above calculation is repeated. (6) Results output: The final pressure of each grid node is obtained after convergence.

[0067] The finite difference method is an existing method, and will not be discussed in detail here.

[0068] Using the above method, the grid node pressure was iteratively updated with the help of node coupling, which made the initial multi-well spatiotemporal feature matrix more complete and obtained the multi-well spatiotemporal feature matrix.

[0069] Step S3: Slice the multi-well spatiotemporal feature matrix, and fuse high-frequency interference and low-frequency interference according to the slicing results to obtain multi-scale spatiotemporal features.

[0070] The grid nodes are sliced ​​based on different time segments and well radii. The length of each time segment and the radius of each well are paired to obtain multiple combinations. Each such combination is denoted as a slice combination. Each slice combination contains one time segment length and one well radius. For example: ,in, and These represent the radii of two different well sizes. and These represent the durations of two different time segments. In this embodiment, the durations of the time segments are divided into two types: 30 days and 180 days. In specific applications, the implementer will set the duration according to the specific circumstances.

[0071] The aggregation of individual slice combinations is evaluated based on their distribution. The calculation formulas for slice aggregation under different slice combinations are as follows:

[0072] in, Indicates the first The aggregation value of the slice combination; K represents the first slice combination. The number of time segments corresponding to a certain combination of slices; Indicates the time segment number; Represents a grid node; This represents the element at the corresponding position in the multi-well spatiotemporal feature matrix; Indicates the well distance index; This represents the stratification scale related to well spacing, i.e., the th The radius in a slice combination; This indicates the stratification scale associated with the time axis, i.e., the first... Duration in slice combinations; Indicates the first The set of grid nodes corresponding to a certain combination of slices; This represents a step function used to perform a hard cut within a radius. This represents the weighting function for selecting data within a specific range on the time axis. In this embodiment, the weighting function is a Gaussian function. This represents the center moment of the k-th time segment.

[0073] Well Distance Index The method for obtaining it is: to extract the nodes in the kth time segment. The average distance between all its nearest wells is used as the well distance index. .

[0074] Next, a local interpolation operator is used to identify high-frequency interference caused by inter-well disturbances. High-frequency interference can be represented as:

[0075] in, Indicates the first High-frequency interference corresponding to the combination of slices This represents a local interpolation operator based on a short radius, used to amplify local gradient responses; Indicates the first The aggregation value of the slice combination; Indicates the relationship with the first Mesh nodes under a combination of slices The relevant weighted function values.

[0076] The weighted function value associated with the grid node is used to emphasize highly coupled nodes, that is, nodes with high coupling have a higher weight in the local interpolation calculation. It can be obtained by taking the proportion of the node coupling of each grid node in the sum of the node coupling of all grid nodes in the slice combination as its associated weighted function value.

[0077] By using the weighted function value to the first The weighted sum of all grid nodes under the slice combination is performed, and the result is obtained for the first slice. The high-frequency interference of the slice combination was quantified.

[0078] Furthermore, wide-domain aggregation is performed on the selected slice combination to obtain low-frequency interference. A fusion operator is then used to perform hierarchical weighting of the slowly varying effects and overall pressure propagation over a wide domain. The relationship of the low-frequency interference can be expressed as:

[0079] in, Indicates the first Low-frequency interference corresponding to the combination of slices, and These represent the start and end times of the selected time segment, respectively. Indicates the first Mesh nodes under a combination of slices Thresholds for center and radius The corresponding spatial domain; Indicates the first The aggregation value of the slice combination; This represents an artificial kernel function for slowly varying effects.

[0080] An artificial kernel function for slowly varying effects is used to reflect the overall pressure propagation, and its calculation formula is as follows: ,in, The duration from the selected time segment to the current time. This represents the maximum production time. and These represent permeability and porosity, respectively. The longer the time and the better the permeability of the reservoir, the larger the value of the corresponding kernel function, highlighting the pressure propagation characteristics under long periods.

[0081] Furthermore, unreasonable anomalies are removed by utilizing geological constraint information and known fluid flow patterns in multi-well blocks. High-frequency and low-frequency interferences are weighted and synthesized into multi-scale spatiotemporal features. The weighted synthesis formula is as follows: ,in, Indicates the first Multi-scale spatiotemporal characteristics of slice combinations and Let represent the weights for high-frequency interference weighting and low-frequency interference weighting, respectively, satisfying . In this embodiment and Both are 0.5.

[0082] When synthesizing multi-scale spatiotemporal features, the cross-well coupling matrix and reservoir parameter data are fully utilized. By using geological logic to determine whether there are outliers that deviate from the formation properties at certain locations or time periods, local spikes or sudden drops that do not conform to the actual physical mechanism are eliminated, and multi-scale spatiotemporal features are obtained, forming a multi-scale spatiotemporal feature dataset that has both high-frequency and low-frequency information in the time and spatial domains.

[0083] Step S4: Based on the multi-scale spatiotemporal features, an adaptive prior distribution is obtained. Variational inference is used to iterate the geological parameters and the adaptive prior distribution to obtain a posterior distribution. The posterior distribution and the multi-scale spatiotemporal features are interactively corrected to obtain the corrected multi-scale spatiotemporal features, the iteratively updated interference radius, and the inter-well connectivity.

[0084] Next, an adaptive prior distribution is obtained based on multi-scale spatiotemporal features, and a posterior distribution is obtained by iteratively applying variational inference to the measured data and the adaptive prior distribution.

[0085] Specifically, a priori aggregation function is obtained based on multi-scale spatiotemporal characteristics. By constraining the prior distribution of geological parameters, high connectivity leads to enhanced inter-well interference and a significant increase in the aggregation value of the divided slice combinations, while low permeability reduces long-cycle pressure propagation. In other words, the priori aggregation function weights the multi-scale spatiotemporal characteristics through the weight of the slice combinations in the aggregation process and the sensitivity of geological parameters in the overall distribution, transforming the physical influence of geological parameters into a constraint at the feature level. The aim is to limit the prior distribution of geological parameters, i.e., higher multi-scale spatiotemporal characteristics correspond to a high probability interval of connectivity in the geological parameters. In subsequent calculations, when the priori aggregation function is transformed into a normalized probability distribution using the reservoir pressure transmission model and fracture propagation mechanism, a mathematical mapping between multi-scale spatiotemporal characteristics and geological parameters is established, meaning the prior distribution is directly related to the range of geological parameter values. The priori aggregation function can be expressed as:

[0086] in, Let A represent the prior aggregation function, and let A represent the number of slice combinations. Indicates the first The weight of each slice combination in the aggregation process. Indicates the first Multi-scale spatiotemporal characteristics of slice combinations Indicates the first Multi-scale spatiotemporal characteristics of slice combinations The normalized quantization result is obtained by using the min-max normalization method. This indicates the sensitivity of geological parameters within the overall distribution.

[0087] The sensitivity of geological parameters in the overall distribution typically ranges from 0.1 to 10. In practice, regional analysis can be performed based on historical geological exploration and monitoring data; in this embodiment, the value is 1.2. The method for obtaining the weight of the slice combination in the aggregation process is as follows: [The text abruptly ends here, likely due to an incomplete sentence or a formatting error.] The ratio between the multi-scale spatiotemporal features of the first slice combination and the sum of the multi-scale spatiotemporal features of all slice combinations is used as the ratio of the first slice combination to the sum of the multi-scale spatiotemporal features of all slice combinations. The weight of each slice combination in the aggregation process.

[0088] Furthermore, based on the reservoir pressure transmission model and fracture propagation mechanism, the prior aggregation function is transformed into a normalized probabilistic form to obtain an adaptive prior distribution; the specific transformation relationship is as follows: ,in, Indicates an adaptive prior distribution. Let be the likelihood function. Let F represent the geological parameters, and let F represent the dataset of multi-scale spatiotemporal features. The likelihood function is a probability function that quantifies the rationality of the geological parameter values ​​given the multi-scale spatiotemporal features. Its core function is to transform the reservoir dynamic information contained in the multi-scale spatiotemporal features into constraints on the parameters, ensuring that the prior distribution conforms to both the physical mechanism and the actual reservoir observation characteristics. Its purpose is to convert the aggregation function into a probability distribution for the analysis of adaptive prior distribution.

[0089] Variational inference is derived based on the acquired parameter data and adaptive prior distribution. The measured data includes yield and pressure. The objective function for variational inference is:

[0090] in, This represents the function value of the objective function in variational inference. Indicates the well's first Various geological parameters For the expected term, Represents the log-likelihood function; express Divergence is used to measure the difference between a variational distribution and a prior distribution. This indicates the number of geological parameter types, ensuring that each set of data contributes independently to the update process; This represents the variational distribution, which is also the Gaussian distribution of the parameter data; This represents an adaptive prior distribution.

[0091] Through the By performing maximization iterations, the optimal variational distribution under the current parameter data and prior constraints is obtained. Iterative variational inference is then performed to obtain the initial posterior distribution based on the parameter data and the adaptive prior distribution.

[0092] Furthermore, the geological parameters in the initial posterior distribution are projected back to the adaptive prior distribution to form a secondary correction and iterate. The new adaptive prior function can be expressed as:

[0093] in, This represents the new adaptive prior function. Represents the prior correction mapping function; Represents an adaptive prior distribution; This represents the posterior distribution obtained by the convergence of variational inference; A dataset representing multi-scale spatiotemporal features; This represents a decoupling function for multi-scale spatiotemporal features. This decouples multi-scale features into parameter correction factors; This indicates a normalization operation, ensuring that the output is the normalized prior probability distribution.

[0094] This embodiment corrects the prior by using the initial posterior distribution to obtain the posterior distribution, which keeps the inference process synchronized with the spatiotemporal evolution of actual reservoir interference, strengthens the physical consistency and adaptive characteristics of the iterative process, and improves the accuracy of characterizing the dynamic relationship between wells.

[0095] After obtaining the posterior distribution, the connectivity and permeability in the posterior distribution are mapped to multi-scale spatiotemporal features. An error metric function is used to locate spatiotemporal regions that do not conform to connectivity relationships and pressure transmission patterns. The value of the error metric function is a globally accumulated discrete error index, which can be expressed as:

[0096] in, This represents the error measurement function. This indicates that the first characteristic matrix in the multi-scale feature matrix represents the first characteristic matrix. line, number Feature values ​​of column positions Indicates the 1st... under the posterior distribution line, number The inter-well connectivity corresponding to the column; Indicates the 1st... under the posterior distribution line, number The corresponding penetration rate; This represents a fusion function that maps connectivity and permeability to spatiotemporal features. Let J represent the number of rows in the posterior distribution, and J represent the number of columns in the posterior distribution. and These correspond to N2 and N1 in the grid nodes, respectively.

[0097] The multi-scale feature matrix is ​​composed of the multi-scale spatiotemporal features corresponding to the grid nodes. In other words, the multi-scale spatiotemporal features corresponding to the grid nodes are mapped to the positions of the grid nodes in the matrix. The element values ​​in the multi-scale feature matrix are all multi-scale spatiotemporal features.

[0098] Furthermore, based on the connectivity and permeability in the posterior distribution, the interference weights of each spatiotemporal node of the multi-scale spatiotemporal features are corrected for regions with errors exceeding a threshold. The threshold is set by the implementer according to specific circumstances, which will not be elaborated further here; the corrected interference weights can be expressed as:

[0099] in, Represented as grid nodes The corrected interference weights, Represents grid nodes Adaptive weights; Represents grid nodes The connectivity projected onto the posterior distribution, i.e. the true value after posterior distribution correction; Represents grid nodes The connectivity obtained in interference fusion This represents the correction factor.

[0100] The correction coefficient is used to control the compensation level for the impact of errors. The value ranges from 0.1 to 0.5. In this embodiment, the value is 0.3. In specific applications, the value can be selected according to the accuracy of the compensation for the impact of errors, so as to avoid setting too large a value that may cause oscillation.

[0101] like If the actual connectivity of grid node g is higher than initially estimated, it means that the actual intensity of the inter-well interference on this node is stronger than predicted. In this case, the interference weight of this node needs to be increased to make it more important in subsequent feature calculations. If the actual connectivity of grid node g is lower than initially estimated, it means that the actual intensity of the interference between wells is weaker than predicted. In this case, the interference weight of the node needs to be reduced to avoid it having an excessive impact on feature calculation.

[0102] Grid nodes Connectivity obtained in interference fusion It can be represented as: ,in, This indicates that the K-nearest neighbor algorithm is used to obtain the set of wells closest to grid node g. This represents the set of wells where grid nodes g have significant mutual interference. The connectivity of the u-th well in the middle, Represents a set of wells The weight of each well in the initial connectivity estimate of the node.

[0103] The distance from the well to the node in the set is determined, i.e. ,in, This represents the distance from the u-th well to the grid node g. express The maximum value among all distances from wells to grid node g is the connection between the wells and the grid node. The larger the distance, the smaller the impact of the well's connectivity on the grid node connectivity estimation.

[0104] The above correction process operates on each grid node to eliminate the overestimation and underestimation caused by the single scale and lack of geological coupling in the early stage, so that the entire feature map is more consistent with the posterior results in terms of local mutual interference and connectivity gradient.

[0105] The eigenvalues ​​of grid nodes in the multi-well spatiotemporal feature matrix are weighted using the corrected interference weights. Specifically, the product of the corrected interference weights and the corresponding eigenvalues ​​of the grid nodes in the multi-well spatiotemporal feature matrix is ​​used as the corrected eigenvalue. This corrects the eigenvalues ​​with large errors in the multi-well spatiotemporal feature matrix, resulting in a corrected multi-well spatiotemporal feature matrix. Based on the corrected multi-well spatiotemporal feature matrix, the multi-scale aggregation method in step S3 is used to re-slice and aggregate the matrix. This allows for secondary interpolation and fusion of locally dense well areas and regions with large permeability gradients, yielding the corrected multi-scale spatiotemporal features.

[0106] The disturbance radius and connectivity are iteratively updated based on the posterior distribution obtained from variational inference. The iterative update formulas for the disturbance radius and inter-well connectivity are:

[0107] in, Indicates the first In the next iteration, the grid nodes The interference radius, Indicates the first In the next iteration, the grid nodes The interference radius; Indicates the first The grid nodes obtained from the posterior distribution calculation in the next iteration Inter-well connectivity; Indicates the first The grid nodes obtained from the posterior distribution calculation in the next iteration Inter-well connectivity; represents the amplification factor, used to control the response magnitude to differences in connectivity, and exp represents the exponential function with the natural constant as the base.

[0108] The magnification factor ranges from [0.05, 0.2], and in this embodiment, the magnification factor is 0.12.

[0109] Through multiple iterations, the interference radius and connectivity are continuously brought closer to the true values, generating a more physically consistent spatiotemporal feature distribution and posterior parameters for the final model. The local update information obtained from the latest iteration is fed back to the input of prior experience, enabling the inferred connectivity, permeability distribution, and other reservoir parameters to be continuously adjusted synchronously with the real mutual interference information, ultimately obtaining the iteratively updated interference radius and inter-well connectivity.

[0110] Thus, through the above steps, the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity were obtained.

[0111] Step S5: Based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity, the collaborative production capacity of multiple wells is predicted.

[0112] The disturbance radius, inter-well connectivity, permeability, porosity, and fluid compressibility coefficient, which were iteratively updated in step S4, were uniformly read and combined with the improved variational inference to form a multi-well collaborative prediction framework. The framework includes a detailed description of the degree of inter-well interference during initialization, providing accurate basic input for the pressure and production evolution simulation in subsequent steps.

[0113] Specifically, based on the iteratively updated interference radius, inter-well connectivity, and reservoir pressure transmission model, the inter-well interference coupling equation is obtained to predict the pressure and production changes of each well under different operating conditions. The role of the inter-well interference coupling equation is to quantify the dynamic response relationship between the reservoir pressure field and the production of multiple wells, ultimately solving for the production and reservoir pressure distribution of each well under different operating conditions. The inter-well interference coupling equation is as follows:

[0114] in, Represents grid nodes The degree of node coupling, Indicates the number of wells. Represents grid nodes At the moment The pressure; This indicates the corrected penetration rate; This indicates the corrected porosity; This represents the corrected fluid compressibility coefficient; Represents grid nodes Location of well point , The processing result of the source term indicator function; Indicates the first Inter-well connectivity With output Output intensity under combined effect.

[0115] The specific output intensity is as follows: ;in, No. The radius of the well. For the first The thickness of the formation in the block where the well is located. It represents pi (π).

[0116] By solving the above formula and quantifying the output The response was obtained to show the changes in pressure and production of each well over time under different operating conditions.

[0117] When extreme operational events occur, a short-time reset function is used to adjust the inter-well disturbance weights in real time to assess the impact of transient disturbances. Extreme operational events include sudden well shut-in and opening of a new well. The short-time reset function can be expressed as:

[0118] in, Indicates the first Koujing at all times Interference weights; Indicates the first Initial disturbance weights of wells under normal operating conditions; Indicates the correction range; Indicates at time Indicators of extreme operational events; The first digit obtained from the posterior distribution of variational inference is... The connectivity of the wellhead; Indicates the relationship with the first Pressure gradient associated with wellhead, Indicates the starting point of an extreme operational event.

[0119] The value range is [0.01, 0.1]. In this embodiment, The value is 0.05.

[0120] No. Initial disturbance weight of wellhead under normal operating conditions The interference weights are obtained through the correction in step S5. The specific acquisition process is as follows: for each well, according to the number determined in step S2... The set of nodes within the interference radius of adjacent wells of a given well is weighted according to the distance between the node and the well. That is, the distance from each node in the above set to the well is weighted. The negative correlation mapping result of the distance to the wellhead is related to the distance from all nodes in the set to the first wellhead. The ratio between the sums of the negative correlation mapping results of the distances to wells is used as the weight of each node in the set. The weights of all nodes in the set calculated at this point are then used to perform a weighted summation of the corresponding corrected interference weights. This summation result is the weight of the first node. Initial disturbance weight of wellhead under normal operating conditions As a concrete example, the following method retrieves each node down to the specified position. The negative correlation mapping result of the distance between wells: This is the mapping from each node to the first well, with the natural constant as the base and the negative value as the root. The distance to the well is the value of an exponential function of the exponent, which is the distance from each node to the first well. The negative correlation mapping results of the distance between wells.

[0121] At any moment The indicator value for extreme operational events is either 1 or 0. If a sudden well shutdown or opening of a new well occurs, the value is 1; otherwise, it is 0.

[0122] Regarding the correction magnitude: at the first moment of well shut-in, This means reducing the weight when suddenly shutting down a well to minimize interference from neighboring wells; the first moment of opening a well... This means increasing the weight of a new well when it is opened, thereby increasing interference from adjacent wells.

[0123] Interference weight This is used to correct the calculation of production intensity in the inter-well interference coupling equation. Extreme operations can disrupt the original pressure balance of the reservoir, causing sudden changes in the instantaneous pressure gradient (such as rapid pressure recovery after well shut-in or a sudden drop in surrounding pressure when a new well is put into production). Therefore, a short-time reset function dynamically determines the interference weight to adjust the production intensity. The adjustment relationship is as follows: ,in, This indicates the adjusted output intensity.

[0124] It should be noted that after an extreme operational event occurs, the intensity of the generated event is adjusted based on the above method to reduce the impact of the extreme operational event on the inter-well interference coupling analysis.

[0125] In conjunction with the aforementioned short-term reset function, whenever extreme operations occur, the disturbance weights will be corrected based on the real-time transmitted pressure gradient and connectivity changes, thereby accurately assessing the impact of transient disturbances on the pressure and production curves of adjacent wells.

[0126] Furthermore, the predicted results are compared with the posterior distribution. During multi-well prediction, a dynamic correction objective function is used to continuously revise the model output until the predicted results converge to match the posterior distribution. The dynamic correction objective function can be expressed as:

[0127] in, This represents the dynamic correction objective function, where M represents the number of wells. and These represent the first and second parts of the prediction results. The production and pressure of the well; and They represent the first and second halves of the equation, respectively, obtained from the posterior distribution estimation. The production and pressure of the well; This represents the balance coefficient between output and pressure error; This represents the overall correction amount, used to guide the dynamic correction of the model output.

[0128] As a specific example, the balance coefficient between output and pressure error can be obtained as follows: obtain the mean absolute error of output and normalize it, obtain the mean absolute error of pressure and normalize it, and use the ratio between the normalized value of the mean absolute error of output and the normalized value of the mean absolute error of pressure as the balance coefficient between output and pressure error.

[0129] By minimizing The original prediction results are adjusted in real time so that the model gradually converges to the output and pressure consistent with the posterior distribution during the iteration process.

[0130] Furthermore, adaptive management of boundary uncertainties is achieved based on the dynamic correction objective function and the posterior distribution, enabling the posterior distribution to automatically adapt to the seepage and pressure extension effects.

[0131] The method provided in this embodiment processes the input parameters, trains the prediction model based on the processed parameters, and obtains the final prediction data. The inputs to the prediction model include: (1) Static basic data: well coordinates, geological unit block number, permeability tensor, porosity, fluid compressibility coefficient, formation thickness and other inherent properties of reservoir and fluid.

[0132] (2) Dynamic correction parameters: the corrected multi-scale spatiotemporal feature matrix, the interference radius after iterative update, the inter-well connectivity and permeability.

[0133] (3) Production conditions: contingency plans for extreme operating events (such as potential well shut-in / new well commissioning plans).

[0134] (4) Historical reference data: historical production and pressure data of multiple wells after interpolation in step S2.

[0135] Based on the above inputs, the prediction model uses the inter-well interference coupling equation and dynamic correction objective function to finally output the production prediction data sequence for each well.

[0136] After obtaining the results of multi-well collaborative prediction, the results can be visualized and connected to the production and operation system.

[0137] The obtained multi-well prediction results and parameter value ranges are collected uniformly. The parameter value ranges include the upper and lower limits of production for each well in multiple time segments, pressure distribution curves, and inter-well connectivity classification information. A data queue is established. This includes production With pressure The predicted mean of the core quantity and its high and low confidence intervals are calculated. An index matrix is ​​constructed to correlate the spatiotemporal coordinates of different wells with the output prediction intervals. This step outputs... This provides complete data support for subsequent visualization modules, avoiding omissions of multi-well interactions or time-series evolution results.

[0138] A multidimensional parameter mapping operator is used to visualize production rate, pressure, parameter value ranges, and inter-well disturbance weights in a two-dimensional coordinate system. The visualization vector obtained based on the multidimensional parameter mapping operator is as follows:

[0139] in, Indicates that for the first Koujing at all times The generated visual vector; Represents a multidimensional parameter mapping operator; and They represent the first Koujing at all times Production and pressure; and They represent the first The lower and upper limits of the high and low confidence intervals for wellheads; Indicates the inter-well interference coefficient; Indicates the influence factor of boundary flow; Represents the first of multiple mapping channel functions Path output; This represents the fusion weights output by the mapping channel function.

[0140] As a concrete example, the first The method for obtaining the lower and upper limits of the high and low confidence intervals for well n is as follows: Based on the posterior distribution obtained by variational inference, the posterior mean and standard deviation of the production and pressure of well n are extracted. If the posterior distribution is Gaussian, a 95% confidence level (corresponding to 1.96 times the standard deviation) is taken, and the lower and upper limits of production and pressure are calculated respectively using "mean ± 1.96 × standard deviation", which are the lower and upper limits of the high and low confidence intervals. If the posterior distribution is non-Gaussian, the quantile method is used, taking the lower 2.5% quantile as the lower limit and the upper 97.5% quantile as the upper limit to ensure that 95% of the prediction probability is covered and accurately reflect the reliable range of the prediction results.

[0141] As a specific example, the inter-well interference coefficient can be obtained by the following method: by obtaining the interconnectivity between wells after iterative updates, the mean of the connectivity between the nth well and all wells within the interference radius is calculated, and the mean is normalized, that is, mapped to the interval [0, 1]. The mapping result is the inter-well interference coefficient. The closer the value of the inter-well interference coefficient is to 1, the higher the intensity of inter-well interference.

[0142] As a specific example, the boundary flow influence factor can be obtained by: obtaining the flow boundary threshold and the average pressure gradient of the block containing the grid node closest to the well in step S2, and calculating the ratio of the actual distance from the well to the boundary to the flow boundary threshold. If the distance from the well to the boundary is less than the flow boundary threshold, it indicates that the well is significantly affected by the boundary. The closer the distance, the larger the factor value; if the distance from the well to the boundary is greater than or equal to the flow boundary threshold, the well is considered to be unaffected by the boundary, and the boundary flow influence factor is [not specified]. The value is 0.

[0143] As a concrete example, the k-th output of multiple mapping channel functions can be determined as follows: output / pressure ( , The channel uses a linear mapping function to normalize the output and pressure values ​​to the [0, 1] interval, intuitively reflecting the parameter magnitudes; the confidence interval channel... and Using an interval width mapping function, the narrower the interval (the higher the prediction confidence), the closer the output value is to 1; interference coefficient channel A nonlinear S-shaped function is employed to enhance the discrimination within the 0.3–0.7 range, highlighting the difference between moderate and strong / weak interference; boundary influence channels. Using a step-linear hybrid function, when When the value is greater than 0.5, it increases linearly. The value increases slowly when it is ≤0.5, emphasizing wells that are significantly affected by the boundary. The k-th output is the result of the corresponding channel function's calculation of the input parameters, and each output is mapped to the interval [0, 1].

[0144] As a concrete example, the fusion weights output by the mapping channel function The following methods can be used to determine the fusion weight: when production / pressure is used as the core parameter, the fusion weight is 0.4; the fusion weight of the inter-well interference coefficient is 0.3; and the fusion weight of the confidence interval and the boundary flow influence factor is 0.15.

[0145] By using the multi-mapping fusion method described above, the predicted curves and uncertainty boundaries of each well are rendered in a layered manner in a visual plane or three-dimensional space according to the interference coefficient and boundary influence. This makes the interaction between wells more intuitive and helps observers grasp the multi-well interference characteristics and boundary complexity. An interactive panel is obtained based on the visualization vectors, allowing operators to switch between different time windows and geographic blocks. By setting time scroll bars and geographic filters, the data queue can be dynamically linked. The predicted data and interference coefficients in the data allow for a visual switching between local magnification and overall cruise mode.

[0146] By using bridging operators, the visualization results of the interactive panel are connected to the production and operation system to achieve direct linkage between visualization and on-site production and drainage strategies. By establishing a two-way mapping between well group switching schedules, injection and production control curves and various fracturing strategy parameters, the operation instructions in the visualization interface are sent to the operation platform in real time, and the latest production and pressure observation values ​​are obtained from the operation platform.

[0147] Furthermore, the profit function is used to evaluate the increase in production and costs under different control strategies, and the actual output observations are fed back into the model to complete the closed-loop optimization of prediction and practice. The profit function can be expressed as:

[0148] in, Indicates regulatory strategy The overall benefits; Indicates the first Koujing in Strategy The increased production benefits obtained; Indicates the first The cost of fracturing or adjusting the drainage of a well; Indicates the cost conversion factor; , This represents the set of feasible strategies that satisfy the geological and equipment boundaries.

[0149] This revenue function guides users to quickly compare the visual changes and economic benefits of various fracturing or extraction schemes, thereby enabling them to select the optimal scheduling strategy in real-time operations.

[0150] Furthermore, the control results and actual output observations are imported back into the data queue. By comparing the results and fine-tuning the parameters of the aforementioned prediction models, a closed-loop optimization of prediction and practice is ultimately formed. Operators can view the production increase or disturbance changes after the strategy is implemented in real time on the visualization panel, and input the difference information into the multi-dimensional parameter mapping operator and the benefit function to further correct the disturbance coefficient and cost factor, improve the control over reservoir evolution and inter-well dynamics, and ensure that the prediction model can continuously evolve with the feedback from the production site.

[0151] This embodiment first evaluates the cross-well coupling relationship based on the characteristic parameters of multiple wells, corrects outliers in the actual collected pressure and production data, and constructs a multi-well dataset. The construction of the multi-well dataset not only considers the data quality of the measured data itself, but also takes into account the impact of the relationship between the wells on the recovery of outliers, so as to better present the real pressure and production data of various geological blocks. Then, the multi-well spatiotemporal feature matrix is ​​generated by combining the cross-well coupling relationship and the multi-well dataset. Based on the slicing results of the multi-well spatiotemporal feature matrix, high-frequency interference and low-frequency interference are evaluated respectively, capturing the mutual interference signal and pressure propagation effect between the wells, and obtaining multi-scale spatiotemporal features. Furthermore, the prior information of the multi-well geological parameters and cross-well coupling is integrated to iterate the prior distribution to obtain the posterior distribution. On this basis, the multi-scale spatiotemporal features, interference radius and inter-well connectivity are corrected, so that the geological parameters can be analyzed and corrected in real time at multiple scales. This ensures that the linkage relationship of the multi-wells can be accurately identified and continuously updated under complex geological conditions and variable well control strategies, thereby realizing the accurate prediction of the multi-well collaborative production capacity and improving the accuracy of the multi-well collaborative production capacity prediction results.

[0152] Example of a multi-well collaborative productivity prediction system based on inter-well interference coupling: See Figure 2 The diagram illustrates a structural block diagram of a multi-well collaborative production capacity prediction system based on inter-well interference coupling provided by an embodiment of the present invention. The system may include a data acquisition module, a first feature acquisition module, a second feature acquisition module, a correction module, and a production capacity prediction module.

[0153] The data acquisition module is used to acquire characteristic parameters, measured data and geological parameters of each well. The characteristic parameters include coordinates and fracturing morphology. The measured data includes pressure and production. The geological parameters include inter-well connectivity, permeability, fluid viscosity, fluid compressibility coefficient, formation pressure correlation coefficient and formation compressibility factor. The first feature acquisition module is used to determine the cross-well coupling relationship based on the feature parameters of each well, repair outliers in the measured data of each well to obtain a multi-well dataset, and obtain a multi-well spatiotemporal feature matrix based on the cross-well coupling relationship and the multi-well dataset. The second feature acquisition module is used to slice the multi-well spatiotemporal feature matrix and fuse high-frequency interference and low-frequency interference according to the slicing results to obtain multi-scale spatiotemporal features. The correction module is used to obtain an adaptive prior distribution based on the multi-scale spatiotemporal features, and to obtain a posterior distribution by iteratively applying variational inference to the geological parameters and the adaptive prior distribution; the posterior distribution and the multi-scale spatiotemporal features are interactively corrected to obtain the corrected multi-scale spatiotemporal features, the iteratively updated interference radius, and the inter-well connectivity. The production capacity prediction module is used to predict the collaborative production capacity of multiple wells based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity.

[0154] It should be understood that Figure 2 The structural block diagram and modules of the multi-well collaborative production capacity prediction system based on inter-well interference coupling shown can be implemented in various ways. For example, in some embodiments, the system and its modules can be implemented by hardware, software, or a combination of software and hardware. The hardware portion can be implemented using dedicated logic; the software portion can be stored in memory and executed by an appropriate instruction execution system, such as a microprocessor or dedicated hardware. Those skilled in the art will understand that the above-described methods and systems can be implemented using computer-executable instructions and / or included in processor control code, for example, on a carrier medium such as a disk, CD, or DVD-ROM, a programmable memory such as read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The systems and modules of this specification can be implemented not only by hardware circuits such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field-programmable gate arrays, programmable logic devices, etc., but also by software executed by various types of processors, or by a combination of the above-described hardware circuits and software (e.g., firmware).

[0155] For more details about the above modules, please refer to other parts of this manual; they will not be repeated here.

[0156] The system provided is used to execute the corresponding methods described above. Therefore, the beneficial effects it can achieve can be referred to in the beneficial effects described in the corresponding methods described above, and will not be repeated here.

[0157] It should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A multi-well collaborative production capacity prediction method based on inter-well interference coupling, characterized in that, The method includes the following steps: The characteristic parameters, measured data, and geological parameters of each well are obtained. The characteristic parameters include coordinates and fracturing morphology. The measured data include pressure and production. The geological parameters include well connectivity, permeability, fluid viscosity, fluid compressibility coefficient, formation pressure correlation coefficient, and formation compressibility factor. Based on the characteristic parameters of each well, the cross-well coupling relationship is determined, and the outliers in the measured data of each well are repaired to obtain a multi-well dataset; based on the cross-well coupling relationship and the multi-well dataset, a multi-well spatiotemporal feature matrix is ​​obtained. The multi-well spatiotemporal feature matrix is ​​sliced, and high-frequency interference and low-frequency interference are fused according to the slicing results to obtain multi-scale spatiotemporal features; An adaptive prior distribution is obtained based on the multi-scale spatiotemporal features. A posterior distribution is obtained by iteratively applying variational inference to the geological parameters and the adaptive prior distribution. The posterior distribution and the multi-scale spatiotemporal features are interactively corrected to obtain the corrected multi-scale spatiotemporal features, the iteratively updated interference radius, and the inter-well connectivity. This includes: mapping the connectivity and permeability in the posterior distribution to the multi-scale spatiotemporal features; identifying regions with errors in inter-well connectivity and permeability using an error metric function; determining the corrected interference weight for each grid node in regions with errors greater than a threshold based on the connectivity and permeability in the posterior distribution; weighting the eigenvalues ​​of grid nodes in the multi-well spatiotemporal feature matrix using the corrected interference weights to obtain the corrected multi-well spatiotemporal feature matrix; re-slicing and aggregating the multi-well spatiotemporal feature matrix to obtain the corrected multi-scale spatiotemporal features; and iteratively updating the interference radius and connectivity based on the posterior distribution to obtain the corrected permeability, porosity, and fluid compressibility coefficient. Based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity, the collaborative production capacity of multiple wells is predicted.

2. The multi-well collaborative production capacity prediction method based on inter-well interference coupling according to claim 1, characterized in that, The determination of cross-well coupling relationships based on the characteristic parameters of each well includes: The relationship between the two wells is quantified by using their relative positions, as well as the permeability and fluid viscosity of the block in which they are located. A cross-well coupling matrix is ​​constructed based on the well-to-well relationships of every pair of wells. The cross-well coupling matrix is ​​used to reflect the cross-well coupling relationships.

3. The multi-well collaborative production capacity prediction method based on inter-well interference coupling according to claim 2, characterized in that, The process of repairing outliers in the measured data of each well to obtain a multi-well dataset includes: For any kind of measured data: The cumulative difference of each well is quantified by combining the well-to-well relationships of each well with the first difference between any measured data of the other wells and the corresponding reference estimate. The differential accumulation is used to identify and remove outliers in any type of measured data from each well, and the inter-well relationship is used to repair missing values ​​in the measured data of the well. Based on the repaired data, the time axis of each well is standardized by using the ratio between the fracture propagation rate of each well and other wells during the fracturing operation period and the fracturing time of other wells, to obtain the multi-well dataset containing cross-well physical priors.

4. The multi-well collaborative production capacity prediction method based on inter-well interference coupling according to claim 1, characterized in that, The process of obtaining the multi-well spatiotemporal feature matrix based on the cross-well coupling relationship and the multi-well dataset includes: Based on the feature parameters of each well and the minimum distance between the well and the grid node within a single time segment, the feature parameters of the well are projected onto the corresponding grid node corresponding to the minimum distance using an inverse distance weighting algorithm, thereby obtaining the feature value of each grid node in a single time segment; the feature values ​​of all grid nodes constitute the initial multi-well spatiotemporal feature matrix. The adjacent well interference radius of each well is determined based on the cross-well coupling relationship; the adaptive correction weight of the grid node is determined based on the influence of the flow boundary pressure, the distance between the well and the grid node on the grid node, the relationship between the distance between the well and the grid node and the corresponding adjacent well interference radius. Based on permeability, porosity, fluid compressibility coefficient, and adaptive correction weights, the eigenvalues ​​of the grid nodes in the initial multi-well spatiotemporal feature matrix are iterated to obtain the multi-well spatiotemporal feature matrix after updating the grid node coupling degree.

5. The multi-well collaborative production capacity prediction method based on inter-well interference coupling according to claim 1, characterized in that, The step of slicing the multi-well spatiotemporal feature matrix and fusing high-frequency and low-frequency interference based on the slicing results to obtain multi-scale spatiotemporal features includes: Various slice combinations are constructed based on the radius of the well and the length of the time window; Based on the spatiotemporal feature matrix of the multi-well, the overall distance distribution between wells and grid nodes at all acquisition times within a single time segment, and the relationship between the radius of wells in various slice combinations, the aggregate value of various slice combinations is obtained. Using a local interpolation operator, high-frequency interference generated by inter-well mutual interference under various slice combinations is identified based on the aggregated value; The aggregation value and pressure propagation characteristics are fused using a fusion operator to obtain low-frequency interference under various slice combinations. The pressure propagation characteristics are determined based on permeability and porosity. By fusing the high-frequency interference and the low-frequency interference, multi-scale spatiotemporal features under various slice combinations are obtained.

6. The multi-well collaborative production capacity prediction method based on inter-well interference coupling according to claim 1, characterized in that, The adaptive prior distribution obtained based on the multi-scale spatiotemporal features includes: A priori aggregation function is obtained based on the aforementioned multi-scale spatiotemporal features; Based on the reservoir pressure transmission model and fracture propagation mechanism, the prior aggregation function is transformed into a normalized probability form to obtain an adaptive prior distribution.

7. The multi-well collaborative production capacity prediction method based on inter-well interference coupling according to claim 1, characterized in that, The process of iterating the geological parameters and the adaptive prior distribution using variational inference to obtain the posterior distribution includes: The variational inference objective function is derived based on geological parameters and adaptive prior distribution; By maximizing the value of the variational inference objective function, the initial posterior distribution under the constraints of geological parameters and adaptive prior distribution is obtained. The geological parameters in the initial posterior distribution are projected onto the adaptive prior distribution, and the posterior distribution is obtained by a second iteration correction.

8. The multi-well collaborative production capacity prediction method based on inter-well interference coupling according to claim 7, characterized in that, The prediction of multi-well collaborative production capacity based on the corrected multi-scale spatiotemporal characteristics, iteratively updated interference radius, and inter-well connectivity includes: Based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, inter-well connectivity, and reservoir pressure transmission model, the inter-well mutual interference coupling equation is obtained, and the production intensity of each well under the combined effect of inter-well connectivity and production is predicted. When extreme operational events occur, a short-term reset function is used to adjust the inter-well interference weights, correct the calculation of production intensity in the inter-well mutual interference coupling equation, and obtain the prediction results; the extreme operational events include shutting in wells and opening new wells; The predicted results are compared with the posterior distribution, and the model output is revised using a dynamic correction objective function until the preset conditions are met, thus obtaining the final predicted result of multi-well collaborative production capacity.

9. A multi-well collaborative production capacity prediction system based on inter-well interference coupling, the system being used to implement the method of claim 1, characterized in that, The system includes: The data acquisition module is used to acquire characteristic parameters, measured data and geological parameters of each well. The characteristic parameters include coordinates and fracturing morphology. The measured data includes pressure and production. The geological parameters include inter-well connectivity, permeability, fluid viscosity, fluid compressibility coefficient, formation pressure correlation coefficient and formation compressibility factor. The first feature acquisition module is used to determine the cross-well coupling relationship based on the feature parameters of each well, repair outliers in the measured data of each well to obtain a multi-well dataset, and obtain a multi-well spatiotemporal feature matrix based on the cross-well coupling relationship and the multi-well dataset. The second feature acquisition module is used to slice the multi-well spatiotemporal feature matrix and fuse high-frequency interference and low-frequency interference according to the slicing results to obtain multi-scale spatiotemporal features. The correction module is used to obtain an adaptive prior distribution based on the multi-scale spatiotemporal features, and to obtain a posterior distribution by iteratively applying variational inference to the geological parameters and the adaptive prior distribution; the posterior distribution and the multi-scale spatiotemporal features are interactively corrected to obtain the corrected multi-scale spatiotemporal features, the iteratively updated interference radius, and the inter-well connectivity. The production capacity prediction module is used to predict the collaborative production capacity of multiple wells based on the corrected multi-scale spatiotemporal characteristics, the iteratively updated interference radius, and the inter-well connectivity.

Images