Method and system for identifying the whole picture of leakage channel entrance-channel-escape of dam engineering

Abstract

The application discloses a dam engineering leakage channel entrance-channel-escape full profile identification method and system. The method comprises the following steps: acquiring multi-source monitoring data of a dam engineering area; constructing a background seepage flow field reflecting a macroscopic hydraulic gradient and an anisotropic diffusion tensor field reflecting local permeability characteristics based on the multi-source monitoring data; using a non-steady-state convection-diffusion equation set, taking an entrance geometric feature and an escape temperature anomaly as a source respectively, performing bidirectional probability field evolution calculation to obtain an entrance accessibility field and an escape accessibility field which are continuously distributed in a dam body space; and determining a three-dimensional full profile probability distribution of a leakage channel based on the spatial interaction and overlap characteristics of the two fields. Through the construction of an anisotropic coupling feedback and tensor guiding mechanism, the application realizes physical self-consistent reconstruction from discrete monitoring data to continuous channel entities, and effectively solves the topological ambiguity problem of leakage path identification in a complex medium.

Description

Technical Field

[0001] This invention relates to a method for identifying the overall view of seepage channels in dams, and more particularly to a method and system for identifying the overall view of seepage channels in dam engineering from inlet to outlet. Background Technology

[0002] Accurate identification of seepage channels in dams is a core technical challenge for the reinforcement and safe operation and maintenance of water conservancy projects. Accurately understanding the location of underwater seepage inlets, the complex channel paths within the dam body, and the three-dimensional spatial topology of downstream outlets is of crucial engineering and technical value for precisely guiding anti-seepage curtain grouting, blocking the development of piping and seepage damage, and assessing the overall stability of the dam structure.

[0003] Existing seepage detection techniques typically employ high-density electrical resistivity inversion to derive the resistivity distribution within the dam body to infer low-resistivity anomalies, or utilize artificial tracer experiments to determine seepage connectivity. With the development of computational technology, some studies have begun to attempt to model seepage path identification as a shortest path search problem on a discrete grid. This involves constructing a resistance model using resistivity data and employing Dijkstra's algorithm, A* algorithm, or ant colony algorithm to find the minimum-cost path from the hypothetical inlet to the outlet, thereby simulating potential seepage streamlines.

[0004] However, existing discrete path search methods struggle to accurately characterize the physical constraints of medium anisotropy on the seepage direction and lack a full-space dynamic coupling mechanism for multi-source boundary information. Specifically, traditional resistivity inversion only provides a static scalar field, which is insufficient to describe the tensor characteristics of permeability in porous media, leading to scalar-based search paths often violating hydraulic flow direction laws. Furthermore, graph-based unidirectional discrete search algorithms are susceptible to local noise interference and may get trapped in local optima, and they struggle to handle complex topological situations with multiple inlets / outlets or channel bifurcation and convergence, making it difficult to establish a physically self-consistent three-dimensional connectivity profile between inlet geometry and outflow temperature anomalies. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for identifying the overall view of the inlet-channel-outlet of seepage channels in dam engineering, in order to solve at least one of the aforementioned problems in the prior art.

[0006] Technical solution: A method for identifying the overall view of seepage channels in dam engineering, including:

[0007] Acquire multi-source monitoring data of the dam project area, including point cloud data characterizing the geometric features of the underwater inlet, resistivity data characterizing the medium properties of the dam body, and thermal infrared images characterizing downstream outflow temperature anomalies.

[0008] Based on multi-source monitoring data, a background seepage flow field reflecting the macroscopic hydraulic gradient of the dam body and an anisotropic diffusion tensor field reflecting the local seepage characteristics of the dam body are constructed.

[0009] Based on the background seepage flow field and the anisotropic diffusion tensor field, bidirectional probability field evolution calculations are performed using the inlet geometric features and the outflow temperature anomaly region as sources, respectively, to obtain the inlet accessibility field and the outflow accessibility field that are continuously distributed in the dam body space.

[0010] Based on the spatial interaction and overlap characteristics of the inlet accessibility field and the outlet accessibility field, the probability distribution of the three-dimensional overall appearance of the leakage channel is determined.

[0011] Optionally, the construction process of the anisotropic diffusion tensor field includes:

[0012] The resistivity data is mapped to a three-dimensional resistivity model aligned with the spatial coordinates of the dam body.

[0013] Identify low-resistivity anomalous regions in the three-dimensional resistivity model and construct a spatially varying diffusion coefficient tensor;

[0014] Among them, the eigenvalues ​​of the diffusion coefficient tensor are negatively correlated with the resistivity data, so as to enhance the diffusion trend of the probability flow in the low-resistivity anomaly region in the two-way probability field evolution calculation.

[0015] Optionally, the process of constructing the background seepage flow field includes:

[0016] Based on the geometric boundary of the dam and the preset upstream and downstream water level conditions, the steady-state seepage control equation is solved to obtain the global head distribution.

[0017] The negative gradient of the global head distribution is calculated to obtain the three-dimensional hydraulic gradient field, which constitutes the background seepage flow field.

[0018] The background seepage flow field is used to provide a convection velocity vector characterizing the flow direction guidance effect in two-way probability field evolution calculations.

[0019] Optionally, the two-way probability field evolution calculation is based on the unsteady convection-diffusion equations, which include:

[0020] The diffusion term, which characterizes the guiding effect of the medium, uses an anisotropic diffusion tensor field to control the random walk intensity of the probability density in a non-uniform medium, so that the probability flow preferentially diffuses along the low-resistivity region.

[0021] The convection term, which characterizes the guiding effect of the flow field, is used to construct the convection velocity using the hydraulic gradient vector of the background seepage flow field, thereby controlling the migration trend of the probability density along the macroscopic seepage direction.

[0022] Optionally, the unsteady convection-diffusion equations also include coupled reaction terms characterizing the two-way field interaction;

[0023] When calculating the evolution and update of the inlet accessibility field, the coupled reaction term uses the value of the outflow accessibility field at the current location as a gain factor to enhance the propagation tendency of the inlet probability flow to the outflow region.

[0024] When calculating the evolution update of the escape accessibility field, the coupled reaction term uses the value of the inlet accessibility field at the current location as a gain factor to enhance the tendency of the escape probability flow to the inlet region.

[0025] The coupling reaction terms are weighted and controlled by preset coupling reaction coefficients.

[0026] Optionally, the unsteady convection-diffusion equations also include a self-decay term characterizing information dissipation, which is controlled by a preset self-decay coefficient.

[0027] The coupling reaction coefficient is configured to be less than the self-decay coefficient in order to maintain numerical stability during the iterative process of bidirectional probability field evolution calculation and prevent the probability field from diverging.

[0028] Optionally, the bidirectional probability field evolution calculation is implemented using an independent evolution fusion approach, including:

[0029] Solve the forward evolution equation for the entrance reachability field and the reverse evolution equation for the exit reachability field independently until both reach a steady-state distribution.

[0030] After reaching a steady-state distribution, the inlet accessibility field and the outlet accessibility field are multiplied point by point, and the result is normalized to serve as the probability distribution of the three-dimensional overall appearance of the leakage channel.

[0031] Optionally, after determining the probability distribution of the three-dimensional overall appearance of the leakage channel based on the spatial interaction and overlap characteristics of the inlet accessibility field and the outlet accessibility field, the method further includes:

[0032] Three-dimensional threshold segmentation is performed on the probability distribution of the three-dimensional overall view of the leakage channel to extract high-probability connected regions;

[0033] The centerline of high-probability connected regions is extracted using a 3D skeletonization algorithm and used as the 3D geometric skeleton of the leakage channel.

[0034] Based on connectivity analysis of a three-dimensional geometric skeleton, the branching and confluence points of leakage channels are identified, and the topological correspondence between the inlet and outlet regions is established.

[0035] Optionally, the bidirectional probability field evolution calculation is based on iterative numerical solutions, and its termination conditions include:

[0036] Calculate the relative error norm of the probability field values ​​at the current time step and the previous time step;

[0037] When the relative error norm is less than the preset convergence threshold, the evolution calculation is determined to have reached a steady state, the iteration is stopped, and the final entry reachability field and exit reachability field are output.

[0038] A comprehensive identification system for seepage channels in dam engineering, including:

[0039] The data acquisition module is used to acquire multi-source monitoring data of the dam project area, including point cloud data characterizing the geometric features of the underwater inlet, resistivity data characterizing the medium properties of the dam body, and thermal infrared images characterizing the downstream outflow temperature anomaly.

[0040] The background field construction module is used to construct a background seepage flow field that reflects the macroscopic hydraulic gradient of the dam body, and an anisotropic diffusion tensor field that reflects the local seepage characteristics of the dam body, based on multi-source monitoring data.

[0041] The evolution calculation module is used to perform bidirectional probability field evolution calculations based on the background seepage flow field and the anisotropic diffusion tensor field, with the inlet geometric features and the outflow temperature anomaly region as sources, respectively, to obtain the inlet accessibility field and the outflow accessibility field that are continuously distributed in the dam space.

[0042] The overall view recognition module is used to determine the probability distribution of the three-dimensional overall view of the leakage channel based on the spatial interaction and overlap features of the inlet accessibility field and the outlet accessibility field.

[0043] Beneficial effects: By constructing a heterogeneous field coupling feedback and tensor guidance mechanism, this invention achieves a physically self-consistent reconstruction from discrete monitoring data to continuous channel entities, effectively solving the topological ambiguity problem in identifying leakage paths in complex media. Attached Figure Description

[0044] Figure 1 This is a schematic diagram of the overall process of a method for identifying the inlet-channel-outlet of a seepage channel in a dam project, provided in an embodiment of this application.

[0045] Figure 2 This is a schematic diagram of the construction background seepage flow field provided in the embodiments of this application.

[0046] Figure 3 This is a schematic diagram of the process for extracting the skeleton of the leakage channel and performing topological analysis provided in the embodiments of this application.

[0047] Figure 4 This is a schematic diagram of the process for constructing an anisotropic diffusion tensor field provided in the embodiments of this application. Detailed Implementation

[0048] Example 1 provides an overall framework for a method to identify the inlet-channel-outlet full view of seepage channels in dam engineering, such as... Figure 1As shown, this paper elaborates on the overall logic of identifying the overall appearance of seepage channels in dams based on the bidirectional probability field evolution theory. It overcomes the limitations of traditional discrete path search algorithms in complex media, which are prone to getting stuck in local optima or have difficulty handling multi-channel bifurcation. It finds the globally optimal probabilistic connection from the inlet to the outlet region by constructing a continuous physical evolution field.

[0049] Step 101: Obtain multi-source monitoring data for the dam project area. The multi-source monitoring data includes point cloud data characterizing the geometric features of the underwater inlet, resistivity data characterizing the medium properties of the dam body, and thermal infrared images characterizing downstream outflow temperature anomalies.

[0050] In this embodiment, multi-source monitoring data forms the basis for subsequent calculations. Specifically, point cloud data can be acquired using multibeam sonar or underwater 3D scanners, primarily used to capture the geometry of the underwater slope of the dam, particularly reflecting potential seepage inlets, depressions, or structural defects. Resistivity data is typically collected using a high-density resistivity transilluminator, and after inversion calculations, a three-dimensional resistivity distribution within the dam is obtained. This data directly reflects the porosity and water saturation of the medium inside the dam, serving as a crucial physical basis for identifying potential seepage channels. Thermal infrared images are acquired on the downstream slope using an infrared thermal imager to locate thermal anomaly areas caused by seepage flow, i.e., seepage points.

[0051] In practical applications, to ensure spatial consistency of data from different sources, it is necessary to map the aforementioned data to the same three-dimensional Cartesian coordinate system and perform necessary denoising and normalization processing. For example, resistivity data can be interpolated onto regular three-dimensional grid nodes to form a voxelized data field.

[0052] Step 102: Based on multi-source monitoring data, construct a background seepage flow field that reflects the macroscopic hydraulic gradient of the dam body, and an anisotropic diffusion tensor field that reflects the local seepage characteristics of the dam body.

[0053] In this embodiment, this step provides the physical environment for the evolution of the probability field. The background seepage flow field describes the overall flow trend of water within the dam body driven by macroscopic head differences. It is a steady-state flow field calculated based on Darcy's law and the overall geometric boundaries of the dam, independent of specific seepage channel details. This provides the convection velocity vector for subsequent probability evolution, simulating the downstream migration trend of seepage material. The anisotropic diffusion tensor field transforms resistivity data into a measure of diffusion capacity. Specifically, in low resistivity regions, where the medium has high water content or large porosity, the diffusion coefficient of the probability flow is set to a larger value; while in high resistivity regions, the diffusion coefficient is suppressed. By constructing the tensor field, the geophysical data's detection results of medium properties can be mathematically transformed into spatial operators guiding the probability evolution.

[0054] Step 103: Based on the background seepage flow field and the anisotropic diffusion tensor field, a two-way probability field evolution calculation is performed, taking the inlet geometric features and the outflow temperature anomaly region as sources, respectively, to obtain the inlet accessibility field and the outflow accessibility field that are continuously distributed in the dam body space.

[0055] In this embodiment, the bidirectional probability field evolution calculation no longer involves finding a specific geometric line, but rather calculating the evolution process of two scalar fields over time or pseudotime. The inlet accessibility field is a probability distribution that originates from the inlet's geometric features and grows forward into the dam body under convection and diffusion, representing the probability that water flow from the inlet will reach a certain point in space. The outlet accessibility field is a probability distribution that originates from the outlet temperature anomaly region and flows backward into the dam body against the direction of water flow, representing the probability that water flow at a certain point in space will eventually reach that outlet point. This process can be achieved by solving an unsteady convection-diffusion partial differential equation system.

[0056] It should be noted that, depending on the implementation strategy, the evolution of the two fields can be independent of each other or they can be coupled with feedback. This embodiment serves as a generalization covering both scenarios.

[0057] Step 104: Based on the spatial interaction and overlap characteristics of the inlet accessibility field and the outlet accessibility field, determine the probability distribution of the three-dimensional overall appearance of the leakage channel.

[0058] In this embodiment, the final leakage channel is not determined by a single field, but by the intersection of two fields. Specifically, a point is considered part of the leakage channel only if it is easily accessible from both the inlet and the outlet. The spatial interaction and overlap characteristics can be represented mathematically through product operations or the steady-state solution of coupling terms. The resulting three-dimensional overall probability distribution of the leakage channel is a three-dimensional scalar field. Regions with higher values ​​in the field constitute the three-dimensional solid shape of the leakage channel, including possible main channels, branch channels, and fine seepage networks.

[0059] Example 2, as a refinement of steps 102 and 103 in Example 1, describes how to construct the background flow field using Darcy's law, how to construct the anisotropic diffusion tensor using resistivity data, and how to initialize the probability field using a Gaussian function, providing accurate physical boundaries and initial conditions for bidirectional evolution calculations.

[0060] Step 201: Based on the preset geometric boundaries of the dam and the preset upstream and downstream water level conditions, solve the steady-state seepage control equations to obtain the global head distribution; calculate the negative gradient of the global head distribution to obtain the three-dimensional hydraulic gradient field, which constitutes the background seepage flow field; the background seepage flow field is used to provide the convection velocity vector characterizing the flow direction guidance effect in the two-way probability field evolution calculation, such as... Figure 2As shown.

[0061] In this embodiment, constructing the background seepage flow field is the foundation for subsequent calculations of the convection term. A three-dimensional finite element or finite difference mesh model is established based on the dam's design data or measured topography. For boundary conditions, the submerged portion of the upstream slope is set as a first-type boundary condition, i.e., the head value equals the upstream water level; the portion below the downstream outlet point is set as a constant head boundary corresponding to the downstream water level. For the interior of the dam, it is assumed to satisfy the macroscopic Darcy's law of seepage, i.e., solving the Laplace equation:

[0062] div(K sat *▽(H))=0;

[0063] Where div is the divergence operator, ▽ is the gradient operator, H is the global head distribution, and K... sat K is the saturated permeability coefficient of the dam medium. In the simplified treatment of this embodiment, K... sat An isotropic scalar permeability coefficient is used; in a more general case, K... sat This can be expressed as the anisotropic permeability coefficient tensor K. In this case, the steady-state seepage control equation is correspondingly expanded to div(K·▽(H))=0, and the velocity formula is expanded to v. in =(K·J) / n e v in Let n be the seepage velocity, J be the hydraulic gradient, and n be the flow velocity. e To achieve effective porosity, appropriate treatment methods can be selected based on the anisotropic characteristics of the dam material.

[0064] After solving for the global head distribution H(x,y,z) using numerical methods such as the finite element method, its negative gradient is further calculated to obtain the three-dimensional hydraulic gradient field J:

[0065] J = -▽(H);

[0066] Where J is the three-dimensional hydraulic gradient vector, and its direction represents the main flow direction driven by macroscopic potential energy. The three-dimensional hydraulic gradient field will be directly used to construct the convection velocity vector in subsequent steps.

[0067] Step 202: Map the resistivity data to a three-dimensional resistivity model aligned with the spatial coordinates of the dam body; identify low-resistivity anomaly regions in the three-dimensional resistivity model and construct a spatially varying diffusion coefficient tensor; wherein, the eigenvalues ​​of the diffusion coefficient tensor are negatively correlated with the resistivity data, in order to enhance the diffusion trend of the probability flow in the low-resistivity anomaly regions in the two-way probability field evolution calculation, such as... Figure 4 As shown.

[0068] In this embodiment, the construction of the anisotropic diffusion tensor field is crucial for achieving multi-source data fusion. The inverted three-dimensional resistivity data ρ(x,y,z) is mapped onto the same grid nodes as the seepage calculation. To reflect the physical assumption that low-resistivity regions are more likely to be leakage channels, a negative correlation mapping relationship between the diffusion coefficient and resistivity is constructed. A typical construction formula is as follows:

[0069] D base (x)=D0*(1+α*(ρ max -ρ(x)) / (ρ(x)+ε));

[0070] Among them, D base (x) is the reference diffusion scalar at spatial point x, D0 is the preset background diffusion coefficient, ρ(x) is the resistivity value at that point, and ρ max α is the maximum resistivity across the entire field, α is a positive coefficient for adjusting sensitivity, and ε is a small amount to prevent the denominator from being zero.

[0071] In a typical application scenario, the background diffusion coefficient D0 can range from 0.1 to 1.0 m. 2 For values ​​on the order of / s, the sensitivity adjustment coefficient α can be taken from 0.5 to 2.0, and the small amount ε to prevent the denominator from being zero can be taken from 1.0 to 10.0 Ω·m. The specific values ​​of the above parameters can be adjusted according to the dam material properties and the distribution range of resistivity data.

[0072] When constructing the diffusion tensor D, it is necessary to base it on D. base (x) Determine the two diffusion components along the flow direction and perpendicular to the flow direction, respectively. In this embodiment, the following construction method is adopted:

[0073] D par (x)=D base (x);

[0074] D perp (x)=r aniso ·D base (x);

[0075] Among them, D par (x) represents the diffusion component along the hydraulic gradient direction at the current spatial point, i.e., the preferred leakage direction, and is taken as equal to the reference diffusion scalar D. base (x); D perp (x) represents the diffusion component perpendicular to the hydraulic gradient direction; r aniso The preset anisotropy scaling factor satisfies 0 <r aniso <1, its physical meaning is that diffusion in the vertical direction is suppressed, r aniso The smaller the value, the stronger the anisotropy.

[0076] In one alternative implementation, raniso The value ranges from 0.1 to 0.5, and the specific value can be adjusted adaptively according to the degree of anisotropy of the dam material structure. When r aniso When = 1, the tensor degenerates into the isotropic form D = D base (x)·I.

[0077] Based on this, the isotropic diffusion coefficient D iso Defined as:

[0078] D iso (x)=(D par (x)+D perp (x)) / 2=((1+r aniso ) / 2)·D base (x);

[0079] This definition is used to handle the degradation case where the hydraulic gradient approaches zero, ensuring that the tensor field is numerically stable over the entire domain.

[0080] Furthermore, in order to handle anisotropy and enhance connectivity along the flow direction, a diffusion tensor D can be constructed.

[0081] To address the degradation scenario where the hydraulic gradient approaches zero, i.e., when |J| approaches 0, making it difficult to define the direction vector, this embodiment adopts the following preferred hybrid construction strategy:

[0082] When |J|≥J th Time: D=D par *(e J *e J T )+D perp *(Ie J *e J T );

[0083] When |J| <J th Time: D=D iso *I;

[0084] Among them, J th D is the preset hydraulic gradient threshold. par D represents the diffusion component along the flow direction. perp e represents the diffusion component in the vertical flow direction. J Let I be the unit vector along the hydraulic gradient direction, and let e be the three-dimensional identity matrix. J T For e J The transpose of D iso is the isotropic diffusion coefficient.

[0085] Hydraulic gradient threshold J thUsed to handle the neighborhood of numerical zeros in hydraulic gradient fields, preventing unit vector e from being used. J The calculation is abnormal because the denominator approaches zero. In an alternative implementation, J th The absolute value of the hydraulic gradient in the background seepage flow field can be taken as 1% to 5% of the maximum value of the entire domain, and can be adaptively adjusted according to the specific numerical grid accuracy.

[0086] The above segmented processing ensures the numerical stability and physical rationality of the tensor field across the entire field.

[0087] Regarding the construction of the anisotropic diffusion tensor D, geophysical resistivity data is used to characterize the non-uniformity of medium porosity development, and the hydraulic gradient direction is combined to impart spatial guidance. In this embodiment, the diffusion tensor D is not only a scalar coefficient, but also a second-order symmetric tensor that can guide the probability density to preferentially migrate in a predetermined direction. Its construction formula is as follows:

[0088] D=D par *(e J *e J T )+D perp *(Ie J *e J T );

[0089] Where D is the anisotropic diffusion tensor, D par e is the parallel diffusion coefficient along the hydraulic gradient direction. J Let J be the unit vector of the direction of the three-dimensional hydraulic gradient field, e J T D is the transpose matrix of the direction unit vector. perp I is the vertical diffusion coefficient perpendicular to the hydraulic gradient direction, i.e., the diffusion component in the vertical flow direction, and I is the three-dimensional identity matrix.

[0090] The physical meaning of this tensor is that when the resistivity is low and a significant hydraulic gradient exists at a certain point, D par The value of is significantly amplified, causing the probability flow at that point to primarily follow e. J The leakage propagates in the directional direction but is suppressed in the vertical direction. This mechanism ensures that the leakage paths generated by subsequent evolutionary calculations strictly follow the laws of hydraulics, solving the problem of paths violating physical principles in traditional discrete search algorithms.

[0091] Step 203: Using a three-dimensional Gaussian probability distribution function, the discrete entry points in the point cloud data are mapped to the initial state of a continuous entry accessibility field; using a three-dimensional Gaussian probability distribution function, the escape anomaly regions in the thermal infrared image are mapped to the initial state of a continuous escape accessibility field; wherein, the variance parameter of the three-dimensional Gaussian probability distribution function is set based on the spatial resolution or uncertainty range of the monitoring data.

[0092] In this embodiment, to initiate the unsteady evolution calculation, the discrete observation data must be transformed into the initial condition Φ(x,0) of the continuous field. For the entrance reachability field, it is assumed that the coordinates of the entrance center identified from the point cloud data are x. in The initial field is then defined as:

[0093] Φ in (x,0)=exp(-||xx in || 2 / (2*σ in 2 ));

[0094] Where, Φ in (x,0) represents the probability density of the entrance reachability field at the initial time, ||xx in || represents the Euclidean distance from spatial point x to the center of the entrance, and σ represents the distance from x to the center of the entrance. in The uncertainty parameter for the entry location depends on the spatial resolution of the point cloud data or the geometric scale of the concave features.

[0095] Similarly, for an escape reachability field, assuming the center of the escape region identified by the thermal infrared image is x out The initial field is defined as:

[0096] Φ out (x,0)=exp(-||xx out || 2 / (2*σ out 2 ));

[0097] Where, Φ out (x,0) is the probability density of the escape reachability field at the initial time, σ out This represents the spatial range parameter of the escape region. Gaussian smoothing mapping solves the singularity problem of single-point sources in numerical calculations and objectively reflects the positioning error range inherent in the monitoring data itself.

[0098] Example 3 describes a probabilistic field evolution model based on heterogeneous field coupling feedback. As a preferred embodiment of the present invention, it elaborates in detail how to achieve mutual perception and dynamic guidance of the inlet and outlet reachability fields during the spatiotemporal evolution process by constructing a convection-diffusion equation system containing nonlinear coupling terms. Unlike traditional independent evolution or linear superposition methods, the heterogeneous field coupling feedback mechanism introduced in this embodiment can simulate chemotaxis behavior in biology, allowing the probability flow to automatically focus on the optimal path connecting the inlet and outlet, suppressing background noise and enhancing the ability to identify weak leakage channels.

[0099] Step 301, the bidirectional probability field evolution calculation is based on the unsteady convection-diffusion equation set, which includes: a diffusion term characterizing the guiding effect of the medium, which uses an anisotropic diffusion tensor field to control the random walk intensity of the probability density in the non-uniform medium, so that the probability flow preferentially diffuses along the low-resistivity region; and a convection term characterizing the guiding effect of the flow field, which uses the hydraulic gradient vector of the background seepage flow field to construct the convection velocity and control the migration trend of the probability density along the macroscopic seepage direction.

[0100] In this embodiment, the fundamental physical framework for evolutionary computation is the classical Advance-Diffusion Equation (ADE). This equation physically describes two basic mechanisms of matter or energy transport in a medium: convection and diffusion. Mathematically, for any time t and spatial location x, this set of fundamental equations takes the following form:

[0101] dΦ in / dt=div(D*▽(Φ in ))-div(v in *Φ in )-λ*Φ in ;

[0102] dΦ out / dt=div(D*▽(Φ out ))-div(v out *Φ out )-λ*Φ out ;

[0103] Where, Φ in Let Φ be the probability density of the accessibility field at the entrance. out Let v be the probability density of the escape reachability field, D be the anisotropic diffusion tensor, and v be the probability density of the escape reachability field. in v is the forward convection velocity vector. out λ is the reverse convection velocity vector, and λ is the self-attenuation coefficient, used to simulate the natural dissipation of probability during propagation and prevent the infinite accumulation of values.

[0104] For the diffusion term div(D*▽(Φ)), the spatial differences in resistivity data are used to guide the random walk of the probabilistic flow. In low resistivity regions (corresponding to high D values), the diffusion term is larger, indicating that the probability density tends to penetrate into this region; while in high resistivity regions, diffusion is suppressed. This mechanism ensures that the evolution process can accurately capture the inhomogeneities of the dam medium. For the convection term div(v*Φ), macroscopic head potential energy is used to provide directional guidance, ensuring that the main movement direction of the probabilistic flow conforms to the basic laws of hydraulics.

[0105] Step 302: The unsteady convection-diffusion equations also include coupled reaction terms characterizing the two-way field interaction. When calculating the evolution and update of the inlet accessibility field, the coupled reaction terms use the value of the outflow accessibility field at the current location as a gain factor to enhance the propagation tendency of the inlet probability flow to the outflow region. When calculating the evolution and update of the outflow accessibility field, the coupled reaction terms use the value of the inlet accessibility field at the current location as a gain factor to enhance the source tendency of the outflow probability flow to the inlet region. The coupled reaction terms are weighted and controlled by preset coupled reaction coefficients.

[0106] In this embodiment, the evolution of the bidirectional probability field is not independent; dynamic mutual feedback is achieved through a system of partial differential equations containing nonlinear coupling terms. This system of equations defines the entrance reachability field Φ. in With the reachability field Φ out The coherent enhancement process in the spatiotemporal domain. The complete set of evolution equations is as follows:

[0107] dΦ in / dt=div(D*▽(Φ in ))-div(v in *Φ in )-λ*Φ in +β*Φ out *Φ in ;

[0108] Where, dΦ in / dt is the rate of change of the accessibility field over time, div is the divergence operator, D is the anisotropic diffusion tensor, ▽ is the gradient operator, and v in λ is the forward convection velocity vector, β is the self-attenuation coefficient, and Φ is the coupling reaction coefficient. out This represents the escape reachability field value at the current location.

[0109] dΦ out / dt=div(D*▽(Φ out ))-div(v out *Φ out )-λ*Φ out +β*Φ in *Φ out;

[0110] Where, dΦ out / dt is the rate of change of the escape reachability field with time, v out Φ is the reverse convection velocity vector. in This represents the accessibility field value for the current location.

[0111] The dimensions of λ and β are determined by the dimensional consistency constraints of the convection-diffusion equations. When the time variable t is in seconds (s), the units of λ and β are seconds (s). -1 Within the dimensionless time step framework of the algorithm, λ and β can be understood as the dimensionless decay rate and gain rate within each time step. Appropriate parameter values ​​can be determined by performing corresponding conversions based on the actual physical scale of the dam and the algorithm's time step.

[0112] The last term in the formula (β*Φ) out *Φ in This physically simulates the signal coherence enhancement effect: when the probability signal Φ evolves from the inlet... in Move to the area where it is easy to reach the escape point, i.e., Φ out In larger areas, Φ in The growth rate will be nonlinearly accelerated. This mechanism allows computational resources and probabilistic energy to be rapidly concentrated on the real physical connectivity channel, while the probability value of the noise region that deviates from the channel, lacking the support of the field value from the other side, rapidly approaches zero under the action of the self-decay term λ.

[0113] The mechanism of mutual source terms forms a positive feedback loop: once two fields overlap on a potential channel, their values ​​will mutually stimulate each other and grow rapidly on that channel, while the values ​​in regions away from the channel will relatively decay due to the lack of coherent enhancement.

[0114] In some alternative implementations, the coupling response coefficient β can be designed as a spatial function β(x) rather than a global constant. For example, β(x) can be set to be positively correlated with the local resistivity gradient, increasing the coupling strength at interfaces with significant resistivity abrupt changes, thereby further enhancing the sensitivity to detecting leakage channels at formation boundaries.

[0115] Step 303: The unsteady convection-diffusion equation set also includes a self-decay term characterizing information dissipation. The self-decay term is controlled by a preset self-decay coefficient. The coupling reaction coefficient is configured to be less than the self-decay coefficient to maintain numerical stability during the iterative process of bidirectional probability field evolution calculation and prevent the probability field from diverging.

[0116] In this embodiment, the coupling term +β*Φ out *Φ inIt is a nonlinear positive feedback source, and if left uncontrolled, it can easily lead to an exponential explosion of numerical solutions in local regions, compromising computational stability. Therefore, a competitive suppression mechanism must be introduced.

[0117] The self-decay term -λ*Φ physically represents the natural forgetting of information or energy dissipation. To ensure the boundedness of the system, in the preferred configuration of this embodiment, the parameters satisfy the stability constraint condition:

[0118] β < λ. Where β is the coupling reaction coefficient and λ is the self-decay coefficient.

[0119] In other implementations, the boundedness of the system can also be achieved by dynamically adjusting the β value or by employing other numerical stabilization techniques.

[0120] In a typical application scenario, the self-attenuation coefficient λ can take values ​​from 0.01 to 0.1, and the coupling reaction coefficient β can take values ​​from 0.001 to 0.05, while satisfying the constraint β < λ. For example, when λ = 0.05, β can take values ​​from 0.01 to 0.03.

[0121] The physical implication of this condition is that the energy dissipation rate within the system must be greater than the energy multiplication rate brought about by positive feedback. Therefore, even in the extreme case where the two fields completely overlap, the total probability density can be controlled within a finite range, rather than diverging infinitely.

[0122] Furthermore, to address the potential local hotspot problem in three-dimensional non-uniform media—that is, even if β < λ, transient overshoot may still occur in certain extremely high diffusion regions—this embodiment introduces a dynamic truncation mechanism for the field value as a dual safeguard. Specifically, after the iterative calculation at each time step is completed, the following truncation operation is performed:

[0123] Φ(x,t+1)=min(Φ calculated (x,t+1),Φ max );

[0124] Where Φ(x,t+1) is the final field value after truncation, Φ calculated (x,t+1) represents the calculated original field value, Φ max This is a preset upper limit for the probability density, typically set to 1.0. This engineering approach effectively addresses the risk of local numerical overflow in non-uniform fields.

[0125] Step 304: In the unsteady convection-diffusion equations, the direction of the convection velocity in the convection term is configured as follows: for the evolution of the inlet accessibility field, a positive velocity vector consistent with the flow direction of the background seepage field is used to simulate the probability propagation of the downstream migration of the leaked material; for the evolution of the outflow accessibility field, a reverse velocity vector opposite to the flow direction of the background seepage field is used to simulate the probabilistic backtracking process from the outflow point to the leakage source.

[0126] To achieve the co-evolution of inlet-to-outlet and outlet-to-source flow, the configuration of the convection velocity *v* must possess physical self-consistency. For the inlet field, which simulates downstream material migration, the velocity vector direction is consistent with the background seepage field. For the outlet field, which detects the source of the water flow, counter-current evolution is required. The velocity vector is defined as follows:

[0127] v in =(K sat / n e )*J;

[0128] Among them, v in K represents the positive convection velocity at the inlet field. sat n is the saturated permeability coefficient. e J represents the effective porosity and the three-dimensional hydraulic gradient vector.

[0129] v out =-(K sat / n e )*J;

[0130] Among them, v out The velocity is the reverse convection velocity in the outflow field, and the symbol - indicates that the direction is opposite to the three-dimensional hydraulic gradient vector.

[0131] By configuring the velocity field in a mirrored manner, the information flow can be offset and converged within the same physical coordinate system, thereby improving the accuracy of capturing long-distance seepage channels inside complex dams.

[0132] Step 305: In the two-way probability field evolution calculation, the following boundary conditions are adopted: For the non-seepage boundary of the dam model, a zero flux boundary condition is set to constrain the evolution of the probability flow within the geometry of the dam body; For the permeable boundary except for the inlet and outlet regions, an absorption boundary condition or a Dirichlet zero-value boundary condition is set to simulate the dissipation process of the probability flow after it flows out of the dam body.

[0133] In this embodiment, the boundary conditions directly determine the spatial morphology of the probability field evolution. For physically impermeable boundaries such as the bedrock interface and the cutoff wall interface of the dam, the probability flow should not pass through these boundaries or flow in from the outside. Neumann zero-flux boundary conditions are adopted:

[0134] n·(D*▽(Φ)-v*Φ)=0;

[0135] Where D is the anisotropic diffusion tensor, ▽ is the gradient operator, Φ is the probability field value, v is the convection velocity vector, and n is the normal unit vector of the boundary surface. This condition forces the component of the probability flow vector in the normal direction to be zero, as if a particle colliding with a wall is bounced back into the dam, ensuring that the probability mass is conserved within the solution domain.

[0136] This boundary condition physically forces the probabilistic flow to be difficult to penetrate the impermeable interface, simulating the physical boundary of real seepage.

[0137] The upstream slope (non-inlet region) and downstream slope (non-outlet region) of the dam are physically permeable boundaries. If probabilistic flow reaches these regions, it indicates that the flow has exited the inlet-outlet system and should be considered an invalid path and removed. Dirichlet conditions are used to set the absorbing boundary conditions:

[0138] Φ(x boundary ,t)=0;

[0139] Where, x boundary These are the other permeable boundary points besides the predetermined source area. This condition immediately reduces the probability density of reaching these boundaries to zero, simulating information escape and dissipation, preventing reflection interference from invalid paths, and improving the signal-to-noise ratio of channel identification.

[0140] Example 4 describes a fusion recognition method based on independent bidirectional evolution, providing an alternative scheme based on the principle of linear superposition, parallel to Example 3. The coupling feedback mechanism of Example 3 has advantages in signal-to-noise ratio and focusing ability. In scenarios where computing resources are limited or preliminary results need to be obtained quickly, the independent evolution + post-fusion strategy algorithm of this example has a simple structure, is easy to parallelize, and can also effectively solve the path discontinuity problem of traditional discrete search algorithms.

[0141] Step 401, as another way to implement the unsteady convection-diffusion equations, the bidirectional probability field evolution calculation is implemented by independent evolution fusion, including: solving the forward evolution equation for the inlet reachability field and the reverse evolution equation for the outflow reachability field independently until both reach a steady-state distribution.

[0142] In this embodiment, unlike Embodiment 3, the bidirectional evolution process is completely decoupled. That is, the evolution of the entry reachability field does not depend on the value of the exit reachability field, and vice versa. This independence allows the computation of the two fields to be distributed across different processor cores for parallel execution, improving computational efficiency.

[0143] Specifically, the forward evolution equation for the accessibility field degenerates into a standard linear convection-diffusion equation:

[0144] dΦ in / dt=div(D*▽(Φ in ))-div(v*Φ in )-λ*Φ in ;

[0145] Where, Φ in Let v be the inlet accessibility field, v be the positive velocity vector (i.e., background seepage velocity), D be the diffusion tensor, and λ be the attenuation coefficient.

[0146] The inverse evolution equation for the outflow reachability field also adopts a linear form, but with the convection term in the opposite direction:

[0147] dΦ out / dt=div(D*▽(Φ out ))-div((-v)*Φ out )-λ*Φ out ;

[0148] Where (-v) represents the reverse velocity vector, used to simulate the reverse tracing process.

[0149] In this embodiment, neither equation contains nonlinear coupling terms, i.e., β=0. Evolutionary calculations continue until the rate of change of the two fields over time is less than a preset threshold, reaching a so-called steady-state distribution. At this point, Φ in (x) represents the steady-state probability density of reaching point x in space from the entrance, while Φ out (x) represents the steady-state probability density of the flow from point x in space to the exit point.

[0150] Step 402: Perform point-by-point product operation on the inlet accessibility field and the outlet accessibility field after reaching steady-state distribution, and normalize the operation result to use as the probability distribution of the three-dimensional overall view of the leakage channel.

[0151] In this embodiment, the concept of the law of total probability from probability theory is used to determine the overall picture of the leakage path. If we consider the flow of water through a point x as a random event E... x The probability of this event occurring depends on the joint probability of two independent events: one is that the water flows from the inlet to x (event A), and the other is that the water continues to flow from x to the outlet (event B).

[0152] The probability distribution Ψ(x) of the three-dimensional overall view of the leakage channel can be obtained by pointwise product of two steady-state fields:

[0153] Ψ(x)=Φ in (x)*Φ out (x);

[0154] Here, * represents a scalar multiplication operation corresponding to a spatial location.

[0155] This multiplicative fusion strategy has a significant denoising effect. Specifically:

[0156] In the area near the entrance but far from the exit, Φ in It's very big, but Φ out It is very small, and the product Ψ is close to 0;

[0157] In areas near the exit but far from the entrance, Φ out It's very big, but Φ in It is very small, and the product Ψ is also close to 0;

[0158] Only in areas that can be reached both downstream and upstream, i.e., true connecting channels, Φ in and Φ out Only when all values ​​are relatively large will the product Ψ show a significant high peak.

[0159] To facilitate subsequent visualization and skeleton extraction, Ψ(x) usually needs to be normalized so that its numerical range is mapped to the [0,1] interval:

[0160] Ψ norm (x)=(Ψ(x)-Ψ min ) / (Ψ max -Ψ min );

[0161] Among them, Ψ max and Ψ min These represent the maximum and minimum values ​​of the total product, respectively. The normalized field Ψ norm (x) represents the probability distribution of the three-dimensional overall view of the leakage channel in the final output.

[0162] In some alternative implementations, in addition to product fusion, other fusion operators such as weighted geometric mean or harmonic mean can be used to adjust the sensitivity to channel connectivity. For example, a weighted geometric mean can be used:

[0163] Ψ(x)=(Φ in (x)) w1 *(Φ out (x)) w2 ;

[0164] Here, w1 and w2 are preset weighting coefficients, typically set to w1=w2=0.5. This variant allows assigning a higher weight to one field based on the reliability of the monitoring data, such as the ingress point cloud being more accurate than the exit infrared, thus preserving more features of the high-confidence data source in the fusion result.

[0165] Example 5 describes the specific process of transforming the aforementioned continuous partial differential equations into a computer-executable discrete numerical algorithm. Both the coupled equations in Example 3 and the independent equations in Example 4 need to be solved using the numerical method of this example. This example focuses on solving the numerical stability problem of the discrete convection term and the automatic termination criterion for iterative calculation.

[0166] Step 501: Discretize the continuous physical space in which the bidirectional probability field evolution calculation is performed into a three-dimensional grid system, and use the finite difference method or finite element method to discretize the convection-diffusion equations.

[0167] In this embodiment, the computational grid needs to be constructed first. A regular hexahedral voxel grid is typically used to cover the entire dam area. Let the grid step sizes in the x, y, and z directions be Δx, Δy, and Δz, respectively. For the time dimension, a discrete time step Δt is used for iteration.

[0168] For the diffusion term div(D*▽(Φ)) in the equation, a second-order central difference scheme is usually used for discretization to ensure second-order accuracy. For example, the diffusion component in the x-direction can be approximated as:

[0169] (D*dΦ / dx)| (i+1 / 2) ≈D (i+1 / 2) *(Φ i+1 -Φ i ) / Δx;

[0170] Where the subscript i represents the grid node index.

[0171] For the convection term div(v*Φ) in the equation, since spurious oscillations are prone to occur when convection is dominant, this embodiment preferably adopts a first-order upwind scheme to ensure the monotonicity and stability of the numerical value. Specifically, the forward or backward difference is selected according to the positive or negative direction of the flow velocity v:

[0172] When v x When >0: d(v) x *Φ) / dx≈(v x_i *Φ i -v x_(i-1) *Φ i-1 ) / Δx;

[0173] When v x When <0: d(v) x *Φ) / dx≈(v x_(i+1) *Φ i+1 -v x_i *Φ i ) / Δx;

[0174] This processing method ensures that information always travels in the direction of the wind (i.e., the flow rate), avoiding non-physical numerical back-diffusion.

[0175] In some alternative implementations, to further improve computational accuracy and reduce numerical diffusion, higher-order Total Variation Diminishing (TVD) or Weighted Essentially Non-Oscillatory (WENO) schemes can be used to discretize the convection terms, but this will increase computational costs accordingly.

[0176] Step 502: The bidirectional probability field evolution calculation is based on iterative numerical solution. Its termination conditions include: calculating the relative error norm of the probability field values ​​of the current time step and the previous time step; when the relative error norm is less than the preset convergence threshold, it is determined that the evolution calculation has reached a steady state, the iteration is stopped and the final entry reachability field and exit reachability field are output.

[0177] In this embodiment, the evolutionary computation is a time-stepping process. At each step t->t+Δt, the system updates the Φ of the entire field. in and Φ out Numerical analysis is needed to determine whether a system has completed its evolution, i.e., reached a steady state or a quasi-steady state. A convergence criterion must be introduced.

[0178] Specifically, the relative error norm E(t) is defined as follows:

[0179] E(t)=sqrt(Σ(Φ(x,t)-Φ(x,t-1)) 2 ) / sqrt(Σ(Φ(x,t)) 2 );

[0180] Where Σ represents the summation over all grid nodes x in the entire field, and sqrt is the square root operation. This formula calculates the relative rate of change (L2 norm) of the entire field between two adjacent time steps.

[0181] Set the preset convergence threshold ε c For example, 1.0e-5. During the iteration process, E(t) is monitored in real time.

[0182] If E(t) > ε c This indicates that the field is still evolving rapidly, and the calculation for the next time step should continue.

[0183] If E(t) ≤ ε c If the system has reached dynamic equilibrium and the probability distribution has stabilized, then iteration can be stopped and the current Φ can be output. in (x,t) and Φout (x,t) is the final result.

[0184] This automatic convergence control mechanism avoids the blindness of manually specifying the number of iterations, ensuring both the full evolution of the calculation results and preventing unnecessary waste of computational resources. For the coupled model in Example 3, due to the presence of positive feedback, the convergence process may be slow. In this case, the maximum number of iterations, such as N, can also be used. max =10000, as a fallback termination condition.

[0185] Those skilled in the art will understand that after the bidirectional probability field evolution calculation converges, the normalized value of the product field Ψ(x) along the actual seepage path will be significantly higher than that of the background region far from the path. This characteristic allows the algorithm to identify channel regions in the probability concentration from the continuous probability field. The specific normalization contrast depends on the combined effects of factors such as path length, diffusion coefficient, self-attenuation coefficient, and the non-homogeneity of the dam medium.

[0186] Example 6 describes how to extract the engineering geometrically significant leakage channel entities and their topological structures from a fuzzy probability distribution cloud map. Subsequent processing transforms the mathematical results of the probability field into drawings and data that engineers can directly use, such as... Figure 3 As shown.

[0187] Step 601: Perform three-dimensional threshold segmentation on the probability distribution of the three-dimensional overall view of the leakage channel and extract high-probability connected regions.

[0188] In this embodiment, the input data is the normalized global probability field Ψ. norm (x). To distinguish between channels and background, a probability segmentation threshold T needs to be set. prob Typically, this threshold can be adaptively determined using Otsu's method, or set based on engineering experience, such as T. prob =0.6.

[0189] Perform the following splitting operation:

[0190] When Ψ norm (x)≥T prob At that time, the marker point x is a channel candidate point and is assigned a value of 1;

[0191] When Ψ norm (x) <T prob When the marker point x is a background point, its value is assigned to 0;

[0192] After this step, the continuous probability field is transformed into binary volumetric data. A connected component analysis algorithm is applied to filter out excessively small, isolated noise clumps, retaining only the largest volume or connected components linking known inlet / outlet regions as high-probability connected regions. This region visually demonstrates the three-dimensional spatial extent of the seepage channels within the dam body.

[0193] Step 602: Use a 3D skeletonization algorithm to extract the centerline of the high-probability connected region as the 3D geometric skeleton of the leakage channel.

[0194] In this embodiment, the high-probability connected regions exhibit the thickness and shape of the channels. To perform path length calculations and flow network analysis, they need to be simplified into linear skeletons. A 3D skeletonization algorithm is used to extract the central axis of the connected region while preserving its original topology (such as holes and loops).

[0195] Specifically, iterative erosion or a centerline extraction algorithm based on distance transformation can be used. Taking iterative erosion as an example, this algorithm peels away the outer boundary points of binary voxels layer by layer until only lines the width of a single voxel remain. During the peeling process, the topology preservation principle is strictly adhered to to prevent the continuous lines from being broken. The resulting three-dimensional geometric skeleton is a sequence of interconnected three-dimensional coordinate points that accurately depicts the critical path of the seepage flow.

[0196] Step 603: Based on the connectivity analysis of the three-dimensional geometric skeleton, identify the branching and confluence points of the leakage channels and establish the topological correspondence between the inlet and outlet areas.

[0197] In this embodiment, the extracted skeleton is transformed into a graph structure in graph theory, where voxels on the skeleton are nodes and adjacency relationships are edges. By calculating the degree of each skeleton node, i.e., the number of connected neighbor nodes, key topological feature points can be identified:

[0198] A node with Degree=1 is defined as an endpoint, which typically corresponds to the inlet or outlet of a leak.

[0199] Nodes with Degree=2 are defined as PathNodes, forming the main body of the channel;

[0200] Nodes with a Degree > 2 are defined as BranchPoints or Junctions, indicating that the leakage channel branches or merges at this point.

[0201] Based on these feature points, the system can automatically construct the topology structure tree of the channels. For example, it can identify a complex network relationship where the main entrance A branches into channel 1 and channel 2, with channel 1 leading to exit B and channel 2 leading to exit C. This topological correspondence provides an accurate navigation map for engineering remediation, enabling targeted grouting and sealing projects to focus on sealing key branch points or main trunks.

[0202] After obtaining a stable probability distribution Ψ of the three-dimensional overall view of the leakage channel, this embodiment employs a mature image processing method from the prior art for post-processing. Specifically, the Otsu method is used to determine the globally optimal threshold, transforming the continuous probability field into a binary entity; and the channel's central axis is extracted using a morphological thinning algorithm (Skeletonization). Since the aforementioned evolutionary calculations have already removed noise through physical constraints, in this step, simple topological connectivity analysis can accurately identify the channel's bifurcation and convergence points, eliminating the need for complex heuristic searches.

[0203] Example 7 describes the system architecture and hardware module composition for implementing the above method. This system solidifies the abstract algorithm flow into specific physical devices or software functional modules, facilitating deployment and application in practical engineering.

[0204] A comprehensive identification system for seepage channels in dam engineering, including:

[0205] The data acquisition module is used to acquire multi-source monitoring data of the dam project area. The multi-source monitoring data includes point cloud data characterizing the geometric features of the underwater inlet, resistivity data characterizing the medium properties of the dam body, and thermal infrared images characterizing downstream outflow temperature anomalies.

[0206] In this embodiment, the data acquisition module can be composed of integrated hardware interfaces and preprocessing software. Specifically, it includes a serial interface or network interface for connecting to a multibeam sonar system to read point cloud coordinate files; a data transmission interface for connecting to a high-density resistivity transducer to import apparent resistivity or inverse resistivity data; and an image acquisition card or USB interface for connecting to an infrared thermal imager to acquire thermal infrared streaming media or image sequences. The module also has a built-in coordinate transformation unit for converting the local coordinate systems of different devices into the dam engineering coordinate system, such as the dam axis coordinate system.

[0207] The background field construction module is used to construct a background seepage flow field that reflects the macroscopic hydraulic gradient of the dam body, as well as an anisotropic diffusion tensor field that reflects the local seepage characteristics of the dam body, based on multi-source monitoring data.

[0208] In this embodiment, the background field construction module is typically deployed on a high-performance computing workstation or cloud server. It integrates a finite element mesh generation engine and a seepage field solver, such as a kernel based on open-source or commercial seepage calculation software, to perform the flow field calculations in Embodiment 2. It also includes a tensor mapping unit responsible for performing resistivity-diffusion coefficient conversion formula calculations to generate a voxelized tensor data field.

[0209] The evolution calculation module is used to perform bidirectional probability field evolution calculations based on the background seepage flow field and the anisotropic diffusion tensor field, using the inlet geometric features and the outflow temperature anomaly region as sources, respectively, to obtain the inlet accessibility field and the outflow accessibility field that are continuously distributed in the dam body space.

[0210] In this embodiment, the evolutionary computation module is the core of the system's computation. It is configured to run the convection-diffusion equation solver of Embodiment 3 or Embodiment 4. To improve computational speed, this module preferably uses a graphics processing unit (GPU) for parallel acceleration, utilizing GPU parallel computing technologies such as CUDA or OpenCL to distribute the iterative update task of the 3D mesh across thousands of stream processors for parallel execution. This module also has convergence monitoring capabilities, calculating and displaying the error norm curve in real time.

[0211] The overall view recognition module is used to determine the probability distribution of the three-dimensional overall view of the leakage channel based on the spatial interaction and overlap features of the inlet accessibility field and the outlet accessibility field.

[0212] In this embodiment, the overall image recognition module is mainly responsible for data post-processing and visualization. It executes the field fusion, threshold segmentation, and skeleton extraction algorithms described in Embodiment 6. This module is typically equipped with a professional 3D visualization engine, such as commercial VTK or Unity3D, which can overlay and display the calculated probability cloud map, the extracted channel skeleton, and the original dam model. It supports interactive operations such as rotation, scaling, and slicing for viewing, and ultimately outputs an engineering report containing the precise coordinates and topological relationships of the seepage channels.

[0213] Example 8 further illustrates the adaptive adjustment strategy of the present invention in the face of complex and diverse engineering realities, elaborates in detail how the algorithm handles complex topology situations with multiple inlets and multiple outlets, and the degradation processing scheme in the case of missing monitoring data.

[0214] Step 801: For complex leakage networks with multiple inlets and outlets, perform the same bidirectional probability field evolution calculation as in the single-channel case, without modifying the algorithm structure.

[0215] In this embodiment, the field evolution method employed in this invention can perform natural parallel search. When point cloud data detects M potential entry regions and thermal infrared images detect N escape anomaly regions, the system constructs an initial field during the initialization phase using the principle of linear superposition.

[0216] Specifically, for the multi-entry scenario, the initial entry reachability field is constructed using the following formula:

[0217] Φ in_total (x,0)=Σ k=1 M [w k *exp(-||xx k || 2 / (2*σ k 2 ))];

[0218] Where, Φ in_total (x,0) represents the initial accessibility field after multi-source superposition, M is the number of detected potential entrances, k is the entrance index, Σ represents the summation over all entrances, and w k The confidence weight for the k-th entry point is typically set based on the point cloud concavity depth, x. k Let σ be the center coordinate vector of the k-th entrance. k Let ||xx be the spatial extent parameter of the k-th entrance. k || represents the Euclidean distance from spatial point x to the center of the k-th entrance.

[0219] Similarly, for the multi-exit scenario, the initial escape reachability field is constructed using the following formula:

[0220] Φ out_total (x,0)=Σ j=1 N [v j *exp(-||xx j || 2 / (2*σ j 2 ))];

[0221] Where N is the number of detected escape regions, j is the escape region index, and v j The weight of the j-th escape region is usually set based on the infrared temperature difference amplitude, x j Let σ be the center coordinate vector of the j-th escape region. j Let be the spatial range parameter of the j-th escape region.

[0222] During the evolution process, the probability flow automatically splits and converges. If inlet A is connected to outlet B, and inlet C is connected to outlet D, the evolved product field Ψ will automatically exhibit high values ​​on paths AB and CD, and low values ​​between AD and CB. This automatic pairing allows the present invention to efficiently identify complex tree-like or mesh-like leakage structures.

[0223] Step 802: In the case of a downgrade where some monitoring data is missing, a one-way evolution mode is used for approximate identification.

[0224] In practical engineering projects, it may be difficult to obtain effective downstream thermal infrared images, such as when the water level is too high and the outlet is submerged, or when it is difficult to obtain high-precision point clouds. In such cases, the system automatically determines the completeness of the input data and switches to the corresponding degraded working mode.

[0225] Scenario 1: Missing out data, i.e., Φ cannot be constructed. out source.

[0226] At this point, the system automatically switches to a unidirectional forward evolution mode. The equation system degenerates into solving only Φ. in The convection-diffusion equation. Although it is difficult to obtain the exact path of bidirectional focusing, the evolved Φ in The steady-state field can still reflect the predicted destination of the water flow starting from the known inlet, providing range guidance for downstream investigations.

[0227] Scenario 2: Ingress data is missing, meaning Φ cannot be constructed. in source.

[0228] At this point, the system switches to a one-way reverse tracing mode. Only Φ is solved. out The inverse equation. The evolved Φ out The steady-state field reflects which upstream regions the water flowing out of a known outflow point most likely originates from, guiding upstream exploration efforts.

[0229] This flexible degradation mechanism ensures that the system can still play a supporting role in decision-making even when data is incomplete.

[0230] The following scenario illustrates the identification of seepage channels in a homogeneous earth dam. Assume the dam is 200m long, has a maximum height of 30m, an upstream water level of 25m, and a downstream water level of 5m. Multibeam sonar scanning reveals a depression approximately 0.5m in diameter at coordinates (30,22) on the upstream slope, preliminarily identified as a potential seepage inlet. Thermal infrared imaging reveals a temperature anomaly at coordinates (170,8) on the downstream slope, approximately 2.5℃ lower than the surrounding temperature, identified as a seepage outlet. High-density electrical resistivity inversion (EDI) analysis shows a low-resistivity anomaly zone in the middle of the dam, with a resistivity approximately one-third that of the surrounding area, extending roughly from upstream to downstream.

[0231] After processing using the method of this invention, the inlet accessibility field extends downstream from coordinates (30,22), while the outlet accessibility field traces upstream from coordinates (170,8). Through coupled evolution calculations, the two fields significantly overlap in the low-resistivity anomaly zone region. The final identified three-dimensional panorama of the leakage channel presents as a curved path extending from (30,22) to (170,8), which highly matches the low-resistivity zone derived from resistivity inversion.

[0232] In contrast, if the traditional Dijkstra shortest path algorithm is used, it can only provide a single-pixel-width polyline path because it is based on discrete grid search. This makes it difficult to reflect the true width and probability distribution characteristics of the channel, and it is easily affected by resistivity data noise, resulting in path deviation.

[0233] This application abandons the simple scalar drag model and reconstructs resistivity data into an anisotropic diffusion tensor field. By introducing this tensor operator into the convection-diffusion equation, the evolution of probabilistic flow strictly follows the physical laws of the medium's pore structure and hydraulic gradient, achieving a high-fidelity simulation of the seepage physical process. It solves the problem in existing technologies where discrete path search cannot accurately characterize the constraint of medium anisotropy on the seepage direction.

[0234] This application establishes a bidirectional probability field evolution model incorporating nonlinear coupling feedback terms. Through the mutual sensing of the forward migration of the entry field and the reverse tracing of the exit field, and utilizing the gain effect of the coupling terms, the two fields mutually focus and amplify the signal along the actual connected path, effectively suppressing background noise and ensuring that the result converges to the globally optimal physical connected path. This solves the problems of unidirectional search easily getting trapped in local optima and lacking a multi-source information coupling mechanism.

[0235] This application utilizes the natural bifurcation and convergence characteristics of continuous fields, combined with multi-source Gaussian initialization and 3D skeleton extraction techniques, to automatically adapt and reconstruct a tree-like branch network containing multiple inlets, multiple outlets, and complex internal structures without pre-setting the number of paths or pairing relationships. This achieves accurate identification of the three-dimensional full view of leakage channels, solving the problem of difficulty in parsing complex topologies.

[0236] The preferred embodiments of the present invention have been described in detail above. It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

Claims

1. A method for identifying the overall view of seepage channels in dam engineering, from inlet to outlet, characterized in that... include: Acquire multi-source monitoring data of the dam project area, including point cloud data characterizing the geometric features of the underwater inlet, resistivity data characterizing the medium properties of the dam body, and thermal infrared images characterizing downstream outflow temperature anomalies. Based on multi-source monitoring data, a background seepage flow field reflecting the macroscopic hydraulic gradient of the dam body and an anisotropic diffusion tensor field reflecting the local seepage characteristics of the dam body are constructed. Based on the background seepage flow field and the anisotropic diffusion tensor field, bidirectional probability field evolution calculations are performed using the inlet geometric features and the outflow temperature anomaly region as sources, respectively, to obtain the inlet accessibility field and the outflow accessibility field that are continuously distributed in the dam body space. Based on the spatial interaction and overlap characteristics of the inlet accessibility field and the outlet accessibility field, the probability distribution of the three-dimensional overall appearance of the leakage channel is determined. The construction process of an anisotropic diffusion tensor field includes: The resistivity data is mapped to a three-dimensional resistivity model aligned with the spatial coordinates of the dam body. Identify low-resistivity anomalous regions in the three-dimensional resistivity model and construct a spatially varying diffusion coefficient tensor; Among them, the eigenvalues ​​of the diffusion coefficient tensor are negatively correlated with the resistivity data, so as to enhance the diffusion trend of the probability flow in the low-resistivity anomaly region in the two-way probability field evolution calculation. The process of constructing the background seepage flow field includes: Based on the geometric boundary of the dam and the preset upstream and downstream water level conditions, the steady-state seepage control equation is solved to obtain the global head distribution. The negative gradient of the global head distribution is calculated to obtain the three-dimensional hydraulic gradient field, which constitutes the background seepage flow field. The background seepage flow field is used to provide a convection velocity vector characterizing the flow direction guidance effect in two-way probability field evolution calculations; The two-way probability field evolution calculation is based on the unsteady convection-diffusion equations, which include: The diffusion term, which characterizes the guiding effect of the medium, uses an anisotropic diffusion tensor field to control the random walk intensity of the probability density in a non-uniform medium, so that the probability flow preferentially diffuses along the low-resistivity region. The convection term, which characterizes the guiding effect of the flow field, is used to construct the convection velocity using the hydraulic gradient vector of the background seepage flow field, thereby controlling the migration trend of the probability density along the macroscopic seepage direction. After determining the probability distribution of the three-dimensional overall appearance of the leakage channel based on the spatial interaction and overlap characteristics of the inlet and outlet accessibility fields, the following steps are also included: Three-dimensional threshold segmentation is performed on the probability distribution of the three-dimensional overall view of the leakage channel to extract high-probability connected regions; The centerline of high-probability connected regions is extracted using a 3D skeletonization algorithm and used as the 3D geometric skeleton of the leakage channel. Based on connectivity analysis of a three-dimensional geometric skeleton, the branching and confluence points of leakage channels are identified, and the topological correspondence between the inlet and outlet regions is established.

2. The method according to claim 1, characterized in that, The unsteady convection-diffusion equations also include coupled reaction terms that characterize the two-way field interaction. When calculating the evolution and update of the inlet accessibility field, the coupled reaction term uses the value of the outflow accessibility field at the current location as a gain factor to enhance the propagation tendency of the inlet probability flow to the outflow region. When calculating the evolution update of the escape accessibility field, the coupled reaction term uses the value of the inlet accessibility field at the current location as a gain factor to enhance the tendency of the escape probability flow to the inlet region. The coupling reaction terms are weighted and controlled by preset coupling reaction coefficients.

3. The method according to claim 2, characterized in that, The unsteady convection-diffusion equations also include a self-attenuation term characterizing information dissipation, which is controlled by a preset self-attenuation coefficient. The coupling reaction coefficient is configured to be less than the self-decay coefficient in order to maintain numerical stability during the iterative process of bidirectional probability field evolution calculation and prevent the probability field from diverging.

4. The method according to claim 1, characterized in that, The bidirectional probability field evolution calculation is implemented using an independent evolution fusion approach, including: Solve the forward evolution equation for the entrance reachability field and the reverse evolution equation for the exit reachability field independently until both reach a steady-state distribution. After reaching a steady-state distribution, the inlet accessibility field and the outlet accessibility field are multiplied point by point, and the result is normalized to serve as the probability distribution of the three-dimensional overall appearance of the leakage channel.

5. The method according to claim 1, characterized in that, The evolution calculation of the two-way probability field is based on iterative numerical solutions, and its termination conditions include: Calculate the relative error norm of the probability field values ​​at the current time step and the previous time step; When the relative error norm is less than the preset convergence threshold, the evolution calculation is determined to have reached a steady state, the iteration is stopped, and the final entry reachability field and exit reachability field are output.

6. A comprehensive identification system for the inlet-channel-outlet of seepage channels in dam engineering, characterized in that, include: The data acquisition module is used to acquire multi-source monitoring data of the dam project area, including point cloud data characterizing the geometric features of the underwater inlet, resistivity data characterizing the medium properties of the dam body, and thermal infrared images characterizing the downstream outflow temperature anomaly. The background field construction module is used to construct a background seepage flow field that reflects the macroscopic hydraulic gradient of the dam body, and an anisotropic diffusion tensor field that reflects the local seepage characteristics of the dam body, based on multi-source monitoring data. The evolution calculation module is used to perform bidirectional probability field evolution calculations based on the background seepage flow field and the anisotropic diffusion tensor field, with the inlet geometric features and the outflow temperature anomaly region as sources, respectively, to obtain the inlet accessibility field and the outflow accessibility field that are continuously distributed in the dam space. The overall view recognition module is used to determine the probability distribution of the three-dimensional overall view of the leakage channel based on the spatial interaction and overlap features of the inlet accessibility field and the outlet accessibility field. The system implements a method for identifying the overall view of seepage channels in dam engineering as described in any one of claims 1-5.

Images