Method and system for identifying linear high-permeability channels in an aquifer based on optimization of a penalty term

Abstract

The application relates to the technical field of computer numerical simulation and data optimization, and discloses a method and system for identifying linear high-permeability channels in aquifers based on a penalty term optimization, which comprises the following steps: performing numerical water flow simulation based on a two-dimensional numerical model of a preset stratum region to obtain a simulated water head, calculating the difference between the simulated water head and pre-acquired observation data to obtain a fitting term, constructing a penalty term according to the grid data of the two-dimensional numerical model, performing weighted fusion on the fitting term and the penalty term to obtain a target function value, taking minimizing the target function value as an objective to iteratively optimize a preset conductivity field to obtain an optimized conductivity field, generating a two-dimensional distribution map of the preset stratum region according to the optimized conductivity field, identifying linear high-permeability channels in the preset stratum region according to the two-dimensional distribution map, and obtaining a final identification result. The application can improve the accuracy of identifying linear high-permeability channels in aquifers in strata and rock masses.

Description

Technical Field

[0001] The present invention relates to the technical field of computer numerical simulation and data optimization, and in particular to a method and system for identifying linear high-permeability channels in an aquifer based on penalty term optimization. Background Art

[0002] In the field of geological data processing and algorithm optimization, accurately identifying highly permeable channels is a significant technical challenge with significant scientific and engineering value. These channels, often present as fissures or conduits, occupy only approximately 10% of the formation volume but control over 90% of water flow, significantly impacting groundwater flow, solute transport, and formation stability. Their precise identification is crucial in areas such as groundwater resource development (e.g., fissures / karst formations), tunnel construction safety, reservoir leakage assessment, and pollution spread prediction.

[0003] Hydraulic tomography uses a simulation-optimization framework to infer formation hydraulic properties. Its core approach is to calibrate numerical models based on measured hydraulic response data. This technology has proven effective in theoretical simulations, laboratory validation, and practical projects (such as investigations of fractured aquifers). This technology has been used to identify linear high-permeability pathways, and both numerical simulations and field experiments (such as in a Canadian alluvial aquifer) have demonstrated its potential.

[0004] However, existing methods often use models and optimization methods that tend to generate smooth conductivity distributions, making it difficult to accurately reflect the geometric characteristics of linear channels, resulting in ambiguous inversion results. Their recognition accuracy is highly dependent on the number and distribution of observation wells, and the effect decreases significantly when the number of wells is insufficient. In actual applications, drilling costs limit the dense arrangement of well networks. There is also a lack of a dedicated constraint mechanism for linear characteristics, which makes the characterization results less accurate in complex heterogeneous formations. Summary of the Invention

[0005] The present invention provides a method and system for identifying linear high-permeability channels in aquifers based on penalty term optimization, which realizes accurate identification of high-permeability channels in strata and rock masses through numerical simulation and iterative optimization. Its main purpose is to solve the problem of low accuracy in the existing identification of linear high-permeability channels in aquifers in strata and rock masses.

[0006] To achieve the above-mentioned object, the present invention provides a method for identifying linear high-permeability channels in aquifers based on penalty term optimization, comprising:

[0007] Construct a two-dimensional numerical model of the pre-set stratigraphic area;

[0008] Performing numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head;

[0009] Calculating the difference between the simulated hydraulic head and the previously acquired observation data to obtain a fitting term;

[0010] Constructing a penalty term according to the grid data of the two-dimensional numerical model, performing weighted fusion on the fitting term and the penalty term to obtain an objective function value;

[0011] Iteratively optimizing a preset conductivity field with the goal of minimizing the objective function value to obtain an optimized conductivity field;

[0012] A two-dimensional distribution map of the preset formation area is generated according to the optimized conductivity field, and linear high-permeability channels are identified in the preset formation area according to the two-dimensional distribution map to obtain a final identification result.

[0013] Optionally, constructing a two-dimensional numerical model of a predetermined stratum area includes:

[0014] Performing grid division on the preset stratum area to obtain an initial grid model;

[0015] Based on a preset number of channel grids, a preset number of linear high-permeability channels are randomly generated in the initial grid model to obtain a random channel model;

[0016] Setting a first simulation well, a second simulation well, and a third simulation well for the random channel model according to preset coordinates, and setting a fourth simulation well based on the position of the linear high permeability channel in the random channel model to obtain an initial two-dimensional model;

[0017] Setting a buffer zone of a preset area around the initial two-dimensional model to obtain a complete two-dimensional model;

[0018] A first logarithmic conductivity is set for the linear high permeability channel in the complete two-dimensional model, a second logarithmic conductivity is set for the matrix other than the linear high permeability channel in the complete two-dimensional model, and a third logarithmic conductivity is set for the buffer zone in the complete two-dimensional model to obtain a two-dimensional numerical model.

[0019] Optionally, performing numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head includes:

[0020] Setting a finite difference equation for each grid in the two-dimensional numerical model by combining Darcy's law and the continuity equation;

[0021] Calculate the coefficients of the difference equation for each grid based on the conductivity field of its neighboring grids;

[0022] Construct a sparse matrix based on the coefficients of the difference equation for each grid;

[0023] Solving the steady-state hydraulic head of each grid based on the sparse matrix pair to obtain a steady-state hydraulic head set;

[0024] Based on the coordinates of the first simulation well, the second simulation well, the third simulation well, and the fourth simulation well, the steady-state hydraulic heads of corresponding coordinates in the steady-state hydraulic head set are extracted to obtain simulated hydraulic heads.

[0025] Optionally, constructing a penalty term according to the grid data of the two-dimensional numerical model includes:

[0026] Obtaining the total number of grids of the two-dimensional numerical model;

[0027] Obtaining coordinates of all grids in the two-dimensional numerical model to obtain a grid coordinate set;

[0028] constructing a distance matrix based on the grid coordinate set;

[0029] Constructing a covariance matrix according to the distance matrix and a preset correlation length;

[0030] Calculate a penalty term based on the total number of grids and the covariance matrix;

[0031] Optionally, the penalty term is calculated as follows:

[0032] ;

[0033] ;

[0034] in, Represents the grid in the two-dimensional numerical model and grid The spatial correlation of is the base of natural logarithms, Represents the grid in the two-dimensional numerical model and grid The Euclidean distance of is the preset correlation length, is the penalty term, represents the n-dimensional conductivity parameter vector of the two-dimensional numerical model, n is the total number of grids, is an n-dimensional all-1 vector, is the preset conductivity field average value, is the covariance matrix, represents the inverse matrix of the covariance matrix.

[0035] Optionally, the objective function value is expressed using the following formula:

[0036] ;

[0037] in, is the objective function value, is the fitting term, is the penalty term, is the preset weight coefficient, represents the n-dimensional conductivity parameter vector of the two-dimensional numerical model, where n is the total number of grids.

[0038] Optionally, the iterative optimization of the preset conductivity field with the goal of minimizing the objective function value includes:

[0039] Updating the conductivity parameter vector according to a preset parameter update amount;

[0040] Recalculate the penalty term based on the updated conductivity parameter vector to obtain an iterative penalty term;

[0041] Recalculate the objective function value according to the iterative penalty term and the fitting term to obtain a post-iteration function value;

[0042] Calculating the gradient of the objective function value and the iterated function value;

[0043] Determine whether the gradient value is less than or equal to a preset gradient threshold or whether the current iteration number reaches a preset iteration number;

[0044] If the gradient value is greater than the gradient threshold and the current iteration number does not reach a preset iteration number, updating the conductivity field according to a preset step size, recalculating the fitting term according to the updated conductivity field, and returning to the step of updating the conductivity parameter vector according to a preset parameter update amount;

[0045] If the gradient value is less than or equal to the gradient threshold or the current iteration number reaches a preset iteration number, it is confirmed that the iteration is completed and the conductivity field is confirmed to be the optimized conductivity field.

[0046] Optionally, generating a two-dimensional distribution map of the preset formation area according to the optimized conductivity field includes:

[0047] Setting geographic coordinates for each grid of the two-dimensional numerical model;

[0048] Draw a two-dimensional map based on the geographic coordinates of each grid to obtain an initial two-dimensional map;

[0049] confirming the conductivity of each grid of the two-dimensional numerical model according to the optimized conductivity field;

[0050] identifying high conductivity grids in the two-dimensional numerical model based on a preset conductivity threshold and the conductivity of each grid;

[0051] According to the identified high conductivity grids, a highlighted area is drawn in the initial two-dimensional map to obtain a two-dimensional distribution map.

[0052] Optionally, the identifying of linear high permeability channels in the preset formation area according to the two-dimensional distribution map to obtain a final identification result includes:

[0053] Identifying the highlighted area of ​​the two-dimensional distribution map to obtain a high permeability channel identification result;

[0054] Performing coherence identification and boundary clarity identification on the high permeability channel identification result to obtain a quality identification result;

[0055] Determining whether the two-dimensional distribution map meets the requirements according to the quality identification result;

[0056] If the two-dimensional distribution graph does not meet the requirements, determining that the accuracy of the two-dimensional distribution graph is low;

[0057] If the two-dimensional distribution map meets the requirements, then calculating the similarity between the high permeability channel identification result and the pre-acquired field survey data;

[0058] Determining whether the similarity meets preset similarity requirements;

[0059] If the similarity does not meet the requirement, it is determined that the accuracy of the high permeability channel identification result is low;

[0060] If the similarity meets the requirements, the high permeability channel identification result is confirmed as the final identification result.

[0061] In order to solve the above problems, the present invention also provides a linear high-permeability channel identification system for aquifers based on penalty term optimization, the system comprising a numerical modeling module, a simulation calculation module, a function value calculation module, an inversion iteration module and a channel identification module, wherein:

[0062] The numerical modeling module is used to construct a two-dimensional numerical model of a preset stratum area;

[0063] The simulation calculation module is used to perform numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head, calculate the difference between the simulated water head and pre-acquired observation data, and obtain a fitting term;

[0064] The function value calculation module is used to construct a penalty term according to the grid data of the two-dimensional numerical model, and perform weighted fusion on the fitting term and the penalty term to obtain an objective function value;

[0065] The inversion iteration module is used to iteratively optimize the preset conductivity field with the goal of minimizing the objective function value to obtain an optimized conductivity field;

[0066] The channel identification module is used to generate a two-dimensional distribution map of the preset formation area according to the optimized conductivity field, and perform linear high permeability channel identification on the preset formation area according to the two-dimensional distribution map to obtain a final identification result.

[0067] The embodiment of the present invention constructs a two-dimensional numerical model of a preset formation area, combines Darcy's law with the continuity equation to perform numerical water flow simulation, effectively obtains simulated water head and calculates fitting terms; constructs a penalty term through grid data, utilizes the covariance matrix to constrain the spatial consistency of the conductivity field, and weightedly fuses the fitting term and the penalty term to form an objective function value; it iteratively optimizes the conductivity field with the goal of minimizing the objective function value, generates a two-dimensional distribution map, and realizes the identification of linear high permeability channels through consistency, boundary clarity and similarity verification; by introducing a special constraint mechanism for linear characteristics, it breaks through the limitations of traditional methods such as smooth conductivity distribution, dependence on the number of observation wells and inaccurate characterization results, eliminates the need for dense well network layout, avoids high cost and complex operations, and significantly improves the accuracy and reliability of the identification of linear high permeability channels in aquifers based on penalty term optimization through full-process model optimization and multiple verifications. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 A schematic flow chart of a method for identifying linear high-permeability channels in aquifers based on penalty term optimization according to an embodiment of the present invention;

[0069] Figure 2 A plan layout diagram of a predetermined stratum area provided in one embodiment of the present invention;

[0070] Figure 3 A conductivity field comparison diagram of an inversion result provided by one embodiment of the present invention;

[0071] Figure 4 This is a functional module diagram of a system for identifying linear high-permeability channels in aquifers based on penalty term optimization provided by one embodiment of the present invention.

[0072] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

[0073] It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0074] The embodiment of the present application provides a method for identifying linear high-permeability channels in aquifers based on penalty item optimization. The execution subject of the method for identifying linear high-permeability channels in aquifers based on penalty item optimization includes but is not limited to at least one of the electronic devices such as a server and a terminal that can be configured to execute the method provided in the embodiment of the present application. In other words, the method for identifying linear high-permeability channels in aquifers based on penalty item optimization can be executed by software or hardware installed on a terminal device or a server device, and the software can be a blockchain platform. The server includes but is not limited to: a single server, a server cluster, a cloud server or a cloud server cluster, etc. The server can be an independent server, or it can be a cloud server that provides basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communications, middleware services, domain name services, security services, content delivery networks (CDNs), and big data and artificial intelligence platforms.

[0075] Reference Figure 1 FIG. 1 is a flow chart of a method for identifying linear high permeability channels in aquifers based on penalty term optimization according to an embodiment of the present invention. In this embodiment, the method for identifying linear high permeability channels in aquifers based on penalty term optimization includes:

[0076] S1. Construct a two-dimensional numerical model of the preset stratigraphic area.

[0077] In the embodiment of the present invention, Figure 2 As shown, the predetermined formation area may be an area with a size of 100×100 square meters, which is used to represent a test site with a size of 200×200 square meters and contains linear high permeability channels (such as cracks and pipes).

[0078] In an embodiment of the present invention, constructing a two-dimensional numerical model of a predetermined stratum region includes:

[0079] Performing grid division on the preset stratum area to obtain an initial grid model;

[0080] Based on a preset number of channel grids, a preset number of linear high-permeability channels are randomly generated in the initial grid model to obtain a random channel model;

[0081] Setting a first simulation well, a second simulation well, and a third simulation well for the random channel model according to preset coordinates, and setting a fourth simulation well based on the position of the linear high permeability channel in the random channel model to obtain an initial two-dimensional model;

[0082] Setting a buffer zone of a preset area around the initial two-dimensional model to obtain a complete two-dimensional model;

[0083] A first logarithmic conductivity is set for the linear high permeability channel in the complete two-dimensional model, a second logarithmic conductivity is set for the matrix other than the linear high permeability channel in the complete two-dimensional model, and a third logarithmic conductivity is set for the buffer zone in the complete two-dimensional model to obtain a two-dimensional numerical model.

[0084] In detail, the coordinates of the first simulation well may be (-45, -45); the coordinates of the second simulation well may be (-45, 0); the coordinates of the third simulation well may be (-45, 45); and the coordinates of the fourth simulation well may be (0, 45).

[0085] In detail, see Figure 2 Figure 1 shows a plan layout of a predetermined formation area according to an embodiment of the present invention. The coordinates of four simulated wells are indicated in the figure, where P1 represents the first simulated well, P2 represents the second simulated well, P3 represents the third simulated well, and P4 represents the fourth simulated well. K represents the permeability coefficient, a physical parameter that describes the conductivity of water flow in a porous medium (such as soil, rock, or groundwater aquifer). A higher K value indicates higher permeability, allowing water to flow more easily. A lower K value indicates greater resistance to water flow. The rock matrix refers to the relatively homogeneous main part of an underground aquifer, typically composed of dense rock or soil. A hyperpermeable channel refers to an area in an aquifer with significantly higher permeability than the matrix, typically consisting of fractures, karst conduits, or highly permeable sand layers.

[0086] In the embodiment of the present invention, Figure 2 As shown, gridding the preset stratum area refers to dividing the preset stratum area into 50*50 square grids (a total of 2500 grid cells, each cell area is 4m²).

[0087] In the embodiment of the present invention, conductivity is a physical quantity that characterizes the permeability of a formation, and its unit is m / s. A larger value indicates that it is easier for water to flow through the region (for example, the conductivity of a fracture / pipeline is much higher than that of the matrix).

[0088] In the embodiment of the present invention, Figure 2 As shown, the number of preset channel grids may be 250, and the remaining 2250 are matrix grids. The preset number of linear high permeability channels may be 5 linear high permeability channels. The preset formation area contains five randomly distributed linear high permeability channels, which are randomly distributed in position, length and direction to form a network with a branch and ring structure, representing the typical linear high permeability channel characteristics in natural formations. Figure 2 The direction of the fracture zone inferred from geological surveys is also superimposed.

[0089] In an embodiment of the present invention, setting the fourth simulation well based on the position of the linear high permeability channel in the random channel model means setting the fourth simulation well above the linear high permeability channel, which can amplify the contribution of the channel to the head change and make subsequent simulation data more significantly reflect the channel characteristics.

[0090] In detail, the first logarithmic conductivity may be -3 m / s, the second logarithmic conductivity may be -5 m / s, and the third logarithmic conductivity may be -4 m / s. That is, the logarithmic value of the conductivity of the buffer zone is -4 m / s, the logarithmic value of the conductivity of the channel is -3 m / s, and the logarithmic value of the conductivity of the matrix is ​​-5 m / s. At the same time, it can be deduced that the channel conductivity is 100 times that of the matrix conductivity.

[0091] In the embodiment of the present invention, the logarithmic conductivity is a logarithmic form of conductivity, which is convenient for numerical calculation.

[0092] Furthermore, the buffer zone in the complete two-dimensional model is set to 2000×2000 square meters.

[0093] S2. Perform numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head.

[0094] In the embodiment of the present invention, the simulated water head refers to the water head value when the first simulation well, the second simulation well, the third simulation well, and the fourth simulation well reach a steady state during the numerical water flow simulation.

[0095] In an embodiment of the present invention, the numerical water flow simulation based on the two-dimensional numerical model can be a numerical water flow simulation based on a pumping speed of 6 liters / minute, or a numerical water flow simulation based on a pumping speed of 10 liters / minute, wherein the pumping speed is adjusted according to the water volume on the site.

[0096] In detail, the hydraulic head is the mechanical energy possessed by unit weight of groundwater, which comprehensively reflects the contribution of the position, pressure and flow rate of groundwater to energy, and the unit is meter (m) (equivalent to the water level).

[0097] In an embodiment of the present invention, the numerical water flow simulation based on the two-dimensional numerical model to obtain the simulated water head includes:

[0098] Setting a finite difference equation for each grid in the two-dimensional numerical model by combining Darcy's law and the continuity equation;

[0099] Calculate the coefficients of the difference equation for each grid based on the conductivity field of its neighboring grids;

[0100] Construct a sparse matrix based on the coefficients of the difference equation for each grid;

[0101] Solving the steady-state hydraulic head of each grid based on the sparse matrix pair to obtain a steady-state hydraulic head set;

[0102] Based on the coordinates of the first simulation well, the second simulation well, the third simulation well, and the fourth simulation well, the steady-state hydraulic heads of corresponding coordinates in the steady-state hydraulic head set are extracted to obtain simulated hydraulic heads.

[0103] In the embodiment of the present invention, the finite difference equation is set for each grid in the two-dimensional numerical model in combination with Darcy's law and the continuity equation in order to transform the continuous groundwater flow problem into a computable discrete grid model.

[0104] In the embodiment of the present invention, the calculation of the differential equation coefficient of each grid based on the conductivity field of the adjacent grids of each grid is to quantify the influence of the adjacent grids on the water flow of the central grid.

[0105] In the embodiment of the present invention, the sparse matrix is ​​constructed according to the differential equation coefficients of each grid, and the non-zero elements and their positions are stored in a sparse matrix format to improve storage efficiency.

[0106] In the embodiment of the present invention, solving the steady-state hydraulic head of each grid based on the sparse matrix refers to obtaining the hydraulic head values ​​of all grids after the pumping test is stabilized.

[0107] In detail, Darcy's law is a basic theory in the field of seepage mechanics, which describes the seepage law of fluid when it moves in laminar flow in porous media (such as soil, rock, etc.).

[0108] S3. Calculate the difference between the simulated water head and the pre-obtained observation data to obtain a fitting term.

[0109] In the embodiment of the present invention, the pre-acquired observation data may refer to a groundwater head value obtained by performing field measurements on the preset stratum area.

[0110] In detail, the observation data is obtained by setting a monitoring well in the preset formation area and recording the steady-state head value of the preset formation area based on the monitoring well.

[0111] In the embodiment of the present invention, calculating the difference between the simulated water head and the pre-acquired observation data to obtain the fitting term refers to calculating the root mean square error between the simulated water head and the pre-acquired observation data.

[0112] In detail, the root mean square error is used to quantify the overall deviation between the simulation and measured data, and has the characteristic of being insensitive to outliers.

[0113] S4. Construct a penalty term according to the grid data of the two-dimensional numerical model, and perform weighted fusion on the fitting term and the penalty term to obtain an objective function value.

[0114] In this embodiment of the present invention, the penalty term is an optimization technique based on regularization theory, which is used to impose additional constraints on model parameters to prevent the results from deviating from physical reality or overfitting. In the identification of high-permeability channels in aquifers, the penalty term forces the spatial distribution of the conductivity field to conform to the geometric characteristics of linear channels (such as coherence and clear boundaries), thereby improving the rationality and accuracy of the inversion results.

[0115] In the embodiment of the present invention, by introducing a penalty term, the problems of insufficient linear feature characterization, image fuzziness, and strong dependence on the number of observation wells in the existing hydraulic tomography technology when identifying linear high permeability channels can be overcome.

[0116] In an embodiment of the present invention, constructing a penalty term based on the grid data of the two-dimensional numerical model includes:

[0117] Obtaining the total number of grids of the two-dimensional numerical model;

[0118] Obtaining coordinates of all grids in the two-dimensional numerical model to obtain a grid coordinate set;

[0119] constructing a distance matrix based on the grid coordinate set;

[0120] Constructing a covariance matrix according to the distance matrix and a preset correlation length;

[0121] Calculate a penalty term based on the total number of grids and the covariance matrix;

[0122] In detail, the calculation formula of the penalty term is as follows:

[0123] ;

[0124] ;

[0125] in, Represents the grid in the two-dimensional numerical model and grid The spatial correlation of is the base of natural logarithms, Represents the grid in the two-dimensional numerical model and grid The Euclidean distance of is the preset correlation length, is the penalty term, represents the n-dimensional conductivity parameter vector of the two-dimensional numerical model, where n is the total number of grids, is an n-dimensional all-1 vector, is the preset conductivity field average value, is the covariance matrix, represents the inverse matrix of the covariance matrix.

[0126] In detail, the preset conductivity field average value may be -4 m / s.

[0127] In detail, in the calculation formula of the penalty term, the penalty term is based on regularization theory and originates from the multivariate Gaussian distribution assumption. The spatial consistency of the conductivity field is constrained by the covariance matrix. Its form measures the spatial irrationality of the conductivity deviation from the mean: < hour, If the conductivity changes suddenly (such as nonlinear distribution), the penalty term will increase significantly. If it meets the linear coherence, it will be kept at a low value to limit nonlinear mutations within a short distance.

[0128] Specifically, the preset correlation length can be 2 meters. This differs from traditional methods, where the correlation length is typically larger (10-50 meters), resulting in global smoothing and image blur. The present invention matches the grid size to the typical widths of fractures and conduits, ensuring the constraints are focused on the channel scale.

[0129] In detail, since the typical width of a crack or pipe is 1m to 3m, setting the correlation length to 2m can ensure that only the conductivity mutations of adjacent grids are constrained, avoiding the global smoothing caused by the larger correlation length in traditional methods and preserving the local high heterogeneity characteristics of the channel.

[0130] Specifically, setting a penalty term can enhance channel coherence and boundary clarity. In the two-dimensional numerical model, conventional methods produce a gradual transition at the interface between the channel and the matrix, resulting in blurred images. However, the present invention penalizes sudden changes and adjusts the conductivity distribution, concentrating high-value areas in the channel and creating sharper boundaries.

[0131] In the embodiment of the present invention, the conductivity parameter vector is a one-dimensional vector used to describe the conductivity distribution of each grid cell in the study area. The two-dimensional conductivity field of the study area is flattened into a one-dimensional vector by rows or columns. Each element corresponds to the logarithmic value of the conductivity of a grid cell. The n-dimensional conductivity parameter vector can represent the conductivity parameter vector corresponding to a 50×50 grid cell in the two-dimensional numerical model.

[0132] In the embodiment of the present invention, the objective function value is expressed by the following formula:

[0133] ;

[0134] in, is the objective function value, is the fitting term, is the penalty term, is the preset weight coefficient, represents the n-dimensional conductivity parameter vector of the two-dimensional numerical model, where n is the total number of grids.

[0135] S5. Iteratively optimize the preset conductivity field with the goal of minimizing the objective function value to obtain an optimized conductivity field.

[0136] In an embodiment of the present invention, the preset conductivity field may be -4 m / s.

[0137] In the embodiment of the present invention, the conductivity field is the distribution of conductivity (or logarithmic conductivity) in space, and is usually expressed as the conductivity value of each unit in a two-dimensional or three-dimensional grid.

[0138] In an embodiment of the present invention, the iterative optimization of the preset conductivity field with the goal of minimizing the objective function value includes:

[0139] Updating the conductivity parameter vector according to a preset parameter update amount;

[0140] Recalculate the penalty term based on the updated conductivity parameter vector to obtain an iterative penalty term;

[0141] Recalculate the objective function value according to the iterative penalty term and the fitting term to obtain a post-iteration function value;

[0142] Calculating the gradient of the objective function value and the iterated function value;

[0143] Determine whether the gradient value is less than or equal to a preset gradient threshold or whether the current iteration number reaches a preset iteration number;

[0144] If the gradient value is greater than the gradient threshold and the current iteration number does not reach a preset iteration number, updating the conductivity field according to a preset step size, recalculating the fitting term according to the updated conductivity field, and returning to the step of updating the conductivity parameter vector according to a preset parameter update amount;

[0145] If the gradient value is less than or equal to the gradient threshold or the current iteration number reaches a preset iteration number, it is confirmed that the iteration is completed and the conductivity field is confirmed to be the optimized conductivity field.

[0146] In detail, the gradient threshold may be 0.01.

[0147] In detail, the preset number of iterations may be 500 times.

[0148] S6. Generate a two-dimensional distribution map of the preset formation area according to the optimized conductivity field, and identify linear high-permeability channels in the preset formation area according to the two-dimensional distribution map to obtain a final identification result.

[0149] In an embodiment of the present invention, generating a two-dimensional distribution map of the preset formation area according to the optimized conductivity field includes:

[0150] Setting geographic coordinates for each grid of the two-dimensional numerical model;

[0151] Draw a two-dimensional map based on the geographic coordinates of each grid to obtain an initial two-dimensional map;

[0152] confirming the conductivity of each grid of the two-dimensional numerical model according to the optimized conductivity field;

[0153] identifying high conductivity grids in the two-dimensional numerical model based on a preset conductivity threshold and the conductivity of each grid;

[0154] According to the identified high conductivity grids, a highlighted area is drawn in the initial two-dimensional map to obtain a two-dimensional distribution map.

[0155] In an embodiment of the present invention, the identifying of linear high permeability channels in the preset formation area according to the two-dimensional distribution map to obtain a final identification result includes:

[0156] Identifying the highlighted area of ​​the two-dimensional distribution map to obtain a high permeability channel identification result;

[0157] Performing coherence identification and boundary clarity identification on the high permeability channel identification result to obtain a quality identification result;

[0158] Determining whether the two-dimensional distribution map meets the requirements according to the quality identification result;

[0159] If the two-dimensional distribution graph does not meet the requirements, determining that the accuracy of the two-dimensional distribution graph is low;

[0160] If the two-dimensional distribution map meets the requirements, then calculating the similarity between the high permeability channel identification result and the pre-acquired field survey data;

[0161] Determining whether the similarity meets preset similarity requirements;

[0162] If the similarity does not meet the requirement, it is determined that the accuracy of the high permeability channel identification result is low;

[0163] If the similarity meets the requirements, the high permeability channel identification result is confirmed as the final identification result.

[0164] In detail, see Figure 3As shown in the figure, a conductivity field comparison diagram of the inversion result is provided in an embodiment of the present invention, showing schematic diagrams of different conductivity fields obtained by inversion under two environments with and without penalty terms. When the inversion is performed under the penalty term environment, two northeast-oriented linear high permeability channels are identified in the preset formation area, with a width of about 2-3 meters and a conductivity concentrated at 10 -3 m / s to 10 -4 m / s range, the matrix conductivity is close to 10 -5 m / s, where the channel path is Figure 2 The direction of the fracture zone is consistent with the inferred direction.

[0165] In an embodiment of the present invention, when it is determined that the accuracy of the two-dimensional distribution map is low, the preset correlation length can be reduced and the above-mentioned step of generating the two-dimensional distribution map can be re-executed to improve the consistency and boundary clarity of the high-permeability channel identification results.

[0166] In an embodiment of the present invention, when it is determined that the accuracy of the high permeability channel identification result is low, the weight coefficient can be increased to control the proportion of the penalty term in the objective function, thereby enhancing the linear feature and re-executing the above-mentioned step of generating a two-dimensional distribution map.

[0167] In detail, the continuity identification and boundary clarity identification of the high permeability channel identification results refer to determining whether each highlighted area in the high permeability channel identification results has at least one adjacent highlighted area. If there is a highlighted area without an adjacent highlighted area, it is confirmed to be incoherent; and calculating the conductivity gradient of the boundary grid and determining the clarity based on the gradient threshold.

[0168] In detail, the similarity between the high permeability channel identification result and the pre-acquired field survey data can be calculated by using the Jaccard index to measure the overlap between the identification result and the field survey data.

[0169] like Figure 4 , which is a functional module diagram of a system for identifying linear high-permeability channels in aquifers based on penalty term optimization provided by one embodiment of the present invention.

[0170] The penalty-optimization-based linear high-permeability channel identification system 100 described in the present invention can be installed in an electronic device. Depending on the functionality implemented, the penalty-optimization-based linear high-permeability channel identification system 100 can include a numerical modeling module 101, a simulation calculation module 102, a function value calculation module 103, an inversion iteration module 104, and a channel identification module 105. A module, also referred to as a unit, is a series of computer program segments that can be executed by a processor in an electronic device and perform a fixed function. These modules are stored in the memory of the electronic device.

[0171] In this embodiment, the functions of each module / unit are as follows:

[0172] The numerical modeling module 101 is used to construct a two-dimensional numerical model of a preset stratum area;

[0173] The simulation calculation module 102 is used to perform numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head, calculate the difference between the simulated water head and pre-acquired observation data, and obtain a fitting term;

[0174] The function value calculation module 103 is used to construct a penalty term according to the grid data of the two-dimensional numerical model, and perform weighted fusion on the fitting term and the penalty term to obtain an objective function value;

[0175] The inversion iteration module 104 is configured to iteratively optimize a preset conductivity field with the goal of minimizing the objective function value to obtain an optimized conductivity field;

[0176] The channel identification module 105 is configured to generate a two-dimensional distribution map of the preset formation area according to the optimized conductivity field, and identify linear high-permeability channels in the preset formation area according to the two-dimensional distribution map to obtain a final identification result.

[0177] In detail, each module in the penalty-optimized aquifer linear high permeability channel identification system 100 according to the embodiment of the present invention is used in the same manner as described above. Figure 1 The same technical means as the method for identifying linear high permeability channels in aquifers based on penalty term optimization described in , and can produce the same technical effects, will not be repeated here.

[0178] In the embodiments provided herein, it should be understood that the disclosed devices, systems, and methods may be implemented in other ways. For example, the system embodiments described above are merely illustrative. For example, the module division is merely a logical functional division, and actual implementation may employ other division methods.

[0179] The modules described as separate components may or may not be physically separate, and the components shown as modules may or may not be physical units, that is, they may be located in one place or distributed across multiple network elements. Some or all of the modules may be selected to achieve the purpose of the solution of this embodiment according to actual needs.

[0180] In addition, the functional modules in various embodiments of the present invention may be integrated into a single processing unit, each unit may exist physically separately, or two or more units may be integrated into a single unit. The aforementioned integrated units may be implemented in the form of hardware or hardware plus software functional modules.

[0181] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention.

[0182] Therefore, the embodiments should be considered in all respects as illustrative and non-restrictive, and the scope of the invention is defined by the appended claims rather than the foregoing description, and all changes that come within the meaning and range of equivalents of the claims are intended to be embraced therein. Any reference to a figure in a claim should not be construed as limiting the claim to which it relates.

[0183] The embodiments of the present application can acquire and process relevant data based on artificial intelligence technology. Artificial intelligence (AI) refers to the theories, methods, techniques, and application systems that use digital computers or machines controlled by digital computers to simulate, extend, and expand human intelligence, perceive the environment, acquire knowledge, and use that knowledge to achieve optimal results.

[0184] Furthermore, it is clear that the word "comprising" does not exclude other units or steps, and the singular does not exclude the plural. Multiple units or systems recited in a system claim may also be implemented by a single unit or system through software or hardware. Terms such as "first" and "second" are used to indicate names and do not imply any particular order.

[0185] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for identifying linear high permeability channels in aquifers based on penalty term optimization, characterized in that: The method comprises: Gridding a preset formation area to obtain an initial grid model, randomly generating a preset number of linear high-permeability channels in the initial grid model based on a preset number of channel grids to obtain a random channel model, setting a first simulation well, a second simulation well, and a third simulation well for the random channel model according to preset coordinates, setting a fourth simulation well based on the position of the linear high-permeability channel in the random channel model to obtain an initial two-dimensional model, setting a buffer zone of a preset area around the initial two-dimensional model to obtain a complete two-dimensional model, setting a first logarithmic conductivity for the linear high-permeability channels in the complete two-dimensional model, setting a second logarithmic conductivity for the matrix excluding the linear high-permeability channels in the complete two-dimensional model, and setting a third logarithmic conductivity for the buffer zone in the complete two-dimensional model to obtain a two-dimensional numerical model; Performing numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head; Calculating the difference between the simulated hydraulic head and the previously acquired observation data to obtain a fitting term; Constructing a penalty term according to the grid data of the two-dimensional numerical model, performing weighted fusion on the fitting term and the penalty term to obtain an objective function value; Iteratively optimizing a preset conductivity field with the goal of minimizing the objective function value to obtain an optimized conductivity field; A two-dimensional distribution map of the preset formation area is generated according to the optimized conductivity field, and linear high-permeability channels are identified in the preset formation area according to the two-dimensional distribution map to obtain a final identification result.

2. The method for identifying linear high permeability channels in aquifers based on penalty term optimization according to claim 1, wherein: The performing of numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head includes: Setting a finite difference equation for each grid in the two-dimensional numerical model by combining Darcy's law and the continuity equation; Calculate the coefficients of the difference equation for each grid based on the conductivity field of its neighboring grids; Construct a sparse matrix based on the coefficients of the difference equation for each grid; Solving the steady-state hydraulic head of each grid based on the sparse matrix pair to obtain a steady-state hydraulic head set; Based on the coordinates of the first simulation well, the second simulation well, the third simulation well, and the fourth simulation well, the steady-state hydraulic heads of corresponding coordinates in the steady-state hydraulic head set are extracted to obtain simulated hydraulic heads.

3. The method for identifying linear high permeability channels in aquifers based on penalty term optimization according to claim 1, wherein: The constructing of a penalty term according to the grid data of the two-dimensional numerical model comprises: Obtaining the total number of grids of the two-dimensional numerical model; Obtaining coordinates of all grids in the two-dimensional numerical model to obtain a grid coordinate set; constructing a distance matrix based on the grid coordinate set; Constructing a covariance matrix according to the distance matrix and a preset correlation length; A penalty term is calculated according to the total number of grids and the covariance matrix.

4. The method for identifying linear high permeability channels in aquifers based on penalty term optimization according to claim 3, wherein: The penalty term is calculated as follows: ； ； in, Represents the grid in the two-dimensional numerical model and grid The spatial correlation of is the base of natural logarithms, Represents the grid in the two-dimensional numerical model and grid The Euclidean distance of is the preset correlation length, is the penalty term, represents the n-dimensional conductivity parameter vector of the two-dimensional numerical model, n is the total number of grids, is an n-dimensional all-1 vector, is the preset conductivity field average value, is the covariance matrix, represents the inverse matrix of the covariance matrix.

5. The method for identifying linear high permeability channels in aquifers based on penalty term optimization according to claim 4, wherein: The objective function value is expressed using the following formula: ； in, is the objective function value, is the fitting term, is the penalty term, is the preset weight coefficient, represents the n-dimensional conductivity parameter vector of the two-dimensional numerical model, where n is the total number of grids.

6. The method for identifying linear high permeability channels in aquifers based on penalty term optimization according to claim 5, wherein: The iterative optimization of the preset conductivity field with the goal of minimizing the objective function value includes: Updating the conductivity parameter vector according to a preset parameter update amount; Recalculate the penalty term based on the updated conductivity parameter vector to obtain an iterative penalty term; Recalculate the objective function value according to the iterative penalty term and the fitting term to obtain a post-iteration function value; Calculating the gradient of the objective function value and the iterated function value; Determine whether the gradient value is less than or equal to a preset gradient threshold or whether the current iteration number reaches a preset iteration number; If the gradient value is greater than the gradient threshold and the current iteration number does not reach a preset iteration number, updating the conductivity field according to a preset step size, recalculating the fitting term according to the updated conductivity field, and returning to the step of updating the conductivity parameter vector according to a preset parameter update amount; If the gradient value is less than or equal to the gradient threshold or the current iteration number reaches a preset iteration number, it is confirmed that the iteration is completed and the conductivity field is confirmed to be the optimized conductivity field.

7. The method for identifying linear high permeability channels in aquifers based on penalty term optimization according to claim 1, wherein: Generating a two-dimensional distribution map of the preset formation area according to the optimized conductivity field includes: Setting geographic coordinates for each grid of the two-dimensional numerical model; Draw a two-dimensional map based on the geographic coordinates of each grid to obtain an initial two-dimensional map; confirming the conductivity of each grid of the two-dimensional numerical model according to the optimized conductivity field; identifying high conductivity grids in the two-dimensional numerical model based on a preset conductivity threshold and the conductivity of each grid; According to the identified high conductivity grids, a highlighted area is drawn in the initial two-dimensional map to obtain a two-dimensional distribution map.

8. The method for identifying linear high permeability channels in aquifers based on penalty term optimization according to claim 1, wherein: The identifying of linear high permeability channels in the preset formation area according to the two-dimensional distribution map to obtain a final identification result includes: Identifying the highlighted area of ​​the two-dimensional distribution map to obtain a high permeability channel identification result; Performing coherence identification and boundary clarity identification on the high permeability channel identification result to obtain a quality identification result; Determining whether the two-dimensional distribution map meets the requirements according to the quality identification result; If the two-dimensional distribution graph does not meet the requirements, determining that the accuracy of the two-dimensional distribution graph is low; If the two-dimensional distribution map meets the requirements, then calculating the similarity between the high permeability channel identification result and the pre-acquired field survey data; Determining whether the similarity meets preset similarity requirements; If the similarity does not meet the requirement, it is determined that the accuracy of the high permeability channel identification result is low; If the similarity meets the requirements, the high permeability channel identification result is confirmed as the final identification result.

9. A linear high permeability channel identification system for aquifers based on penalty term optimization, characterized in that: The system includes a numerical modeling module, a simulation calculation module, a function value calculation module, an inversion iteration module and a channel identification module, wherein: The numerical modeling module is configured to grid a preset formation area to obtain an initial grid model, randomly generate a preset number of linear high-permeability channels in the initial grid model based on a preset number of channel grids to obtain a random channel model, set a first simulation well, a second simulation well, and a third simulation well for the random channel model according to preset coordinates, set a fourth simulation well based on the position of the linear high-permeability channel in the random channel model to obtain an initial two-dimensional model, set a buffer zone of a preset area around the initial two-dimensional model to obtain a complete two-dimensional model, set a first logarithmic conductivity for the linear high-permeability channels in the complete two-dimensional model, set a second logarithmic conductivity for the matrix in the complete two-dimensional model excluding the linear high-permeability channels, and set a third logarithmic conductivity for the buffer zone in the complete two-dimensional model to obtain a two-dimensional numerical model; The simulation calculation module is used to perform numerical water flow simulation based on the two-dimensional numerical model to obtain a simulated water head, calculate the difference between the simulated water head and pre-acquired observation data, and obtain a fitting term; The function value calculation module is used to construct a penalty term according to the grid data of the two-dimensional numerical model, and perform weighted fusion on the fitting term and the penalty term to obtain an objective function value; The inversion iteration module is used to iteratively optimize the preset conductivity field with the goal of minimizing the objective function value to obtain an optimized conductivity field; The channel identification module is used to generate a two-dimensional distribution map of the preset formation area according to the optimized conductivity field, and perform linear high permeability channel identification on the preset formation area according to the two-dimensional distribution map to obtain a final identification result.

Images