A method for locating seepage channels in dam engineering based on intelligent sensing

Abstract

This invention discloses a method for locating seepage channels in dam engineering based on intelligent sensing, comprising: acquiring underwater acoustic scanning data, internal resistivity detection data, and backwater thermal infrared imaging data of the dam area; extracting geometric anomalies at the dam inlet, low-resistivity anomalies inside the dam body, and temperature anomalies at the dam outlet based on multi-source data, and constructing a unified three-dimensional anomaly feature spatial distribution model; introducing a pre-constructed dam seepage physical field model, and searching and constructing a set of three-dimensional seepage candidate paths within the three-dimensional search space reconstructed by the physical field, starting from the inlet features and ending with the outlet features, based on low-resistivity feature guidance and seepage physical field constraints; and calculating the physical consistency index of the paths and selecting the optimal location result based on fluid mechanics and thermodynamic conservation principles. This invention effectively solves the pseudo-connectivity problem in multi-source fusion and improves the physical authenticity and accuracy of the location results.

Description

Technical Field

[0001] This invention relates to the field of safety monitoring of water conservancy projects, and in particular to a method for locating seepage channels in dam projects based on intelligent sensing. Background Technology

[0002] Seepage in dams is a major cause of piping, soil erosion, and dam structural instability. Precise location of internal seepage channels is a prerequisite for targeted grouting reinforcement. Accurately understanding the spatial three-dimensional trajectory of seepage channels, the correspondence between their inlet and outlet points, is of crucial technical significance for assessing dam seepage stability and developing emergency response plans. However, seepage channels are usually deeply buried within the dam structure, characterized by high concealment, complex morphology, and randomness, making it difficult for traditional single-point monitoring methods to reveal their continuous spatial morphology.

[0003] Currently, dam seepage detection mainly employs high-density electrical resistivity tomography (EDT) to detect internal low-resistivity anomalies, underwater side-scan sonar to identify collapse pits in front of the dam, and thermal infrared imaging to capture temperature differences in seepage water behind the dam. The mainstream integrated detection approach typically projects the anomaly features obtained through these different methods (such as inlet point cloud depressions, internal low-resistivity bodies, and outlet thermal anomalies) directly onto the same three-dimensional geometric model. Technicians primarily infer potential seepage pathways based on the geometrical proximity or linear arrangement of various anomaly regions, using manual experience or simple geometric connectivity algorithms.

[0004] Existing technologies for multi-source data fusion mainly face two major technical problems: spatial pseudo-connectivity interference and lack of physical constraints. Specifically, existing geometric superposition methods ignore the inherent constraints of the seepage physical field, relying solely on spatial proximity for correlation. This easily leads to misclassifying spatially close but hydraulically disconnected isolated anomalies as the same seepage channel, resulting in a high false-channel false alarm rate. Furthermore, existing path search algorithms often perform geometric optimization within Euclidean space, failing to consider the fluid dynamics principle of decreasing head potential energy. The generated paths frequently contain segments that violate physical laws, such as countercurrents or crosscurrents. In addition, the lack of end-to-end quantitative verification mechanisms based on mass or energy conservation makes it difficult to perform closed-loop confidence checks on the inferred seepage channels. Summary of the Invention

[0005] The purpose of this invention is to provide a method for locating seepage channels in dam engineering based on intelligent sensing, in order to solve at least one of the aforementioned problems in the existing technology.

[0006] According to one aspect of this application, a method for locating seepage channels in dam engineering based on intelligent sensing includes:

[0007] Acquire underwater acoustic scanning data of the dam area, resistivity detection data inside the dam, and thermal infrared imaging data of the backwater side;

[0008] Based on multi-source sensing data, geometric anomaly features at the dam entrance, low resistivity anomaly features inside the dam body, and temperature anomaly features at the dam exit are extracted respectively, and a three-dimensional spatial distribution model of anomaly features under a unified spatial coordinate is constructed.

[0009] A pre-constructed physical field model of dam seepage is introduced. Based on the three-dimensional anomaly feature spatial distribution model, in the three-dimensional search space, starting from the inlet geometric anomaly feature and ending from the outflow temperature anomaly feature, a set of three-dimensional leakage candidate paths is searched and constructed based on the low resistance anomaly feature and seepage physical field constraints.

[0010] Based on the conservation principles of fluid mechanics and thermodynamics, the physical consistency index of each path in the candidate path set is calculated, and the optimal path is selected as the result of leakage channel location.

[0011] Beneficial effects: This invention effectively solves the pseudo-connectivity problem in multi-source fusion, and improves the physical authenticity and accuracy of the positioning results. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the overall process of the method for locating seepage channels in dam engineering based on intelligent sensing, provided in the embodiments of this application.

[0013] Figure 2 This is a schematic diagram of the process of searching and constructing a set of three-dimensional leakage candidate paths in a hierarchical equipotential surface space based on the reconstruction of a three-dimensional hydraulic head scalar field, provided in an embodiment of this application.

[0014] Figure 3 This is a schematic diagram of the flow balance verification process based on the application quality conservation principle provided in the embodiments of this application.

[0015] Figure 4 This is a schematic diagram of the process for constructing a physical field model of seepage in a dam, provided in an embodiment of this application. Detailed Implementation

[0016] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0017] Example 1 details the overall technical process of the intelligent sensing-based method for locating seepage channels in dam engineering according to the present invention, such as... Figure 1As shown, by introducing a seepage field physical model and an energy / mass conservation criterion, the technical problem of difficulty in identifying pseudo-anomalies that are spatially close but physically disconnected when multi-source data are fused based solely on spatial proximity in existing technologies is solved.

[0018] Step 101: Acquire underwater acoustic scanning data of the dam area, resistivity detection data inside the dam, and thermal infrared imaging data of the backwater side.

[0019] In this embodiment, the acquisition operation encompasses the comprehensive, three-dimensional perception data collection of the dam. Specifically, underwater acoustic scanning data can be acoustic echo signals or 3D point cloud data collected by a multibeam sonar system mounted on an underwater robot (ROV) or unmanned surface vessel, used to finely characterize the micro-topographic features of the underwater slope in front of the dam, especially potential subsidence pits or eddy current scour marks. Dam internal resistivity detection data can be apparent resistivity data collected using a high-density resistivity transducer through a preset electrode array, used to invert the dielectric electrical distribution within the dam body and reveal potential water-bearing channels. Backwater side thermal infrared imaging data can be temperature field images of the slope behind the dam collected by an infrared thermal imager, used to capture local temperature anomalies caused by seepage. In practical engineering applications, the acquisition of these three types of data is typically triggered by a clock synchronization controller to ensure that all data have a unified reference on the timeline. For example, acoustic emission, electrical excitation, and infrared imaging can be triggered simultaneously at certain distances (e.g., every 2 meters) or at certain time intervals (e.g., every 1 second), providing a solid spatiotemporal foundation for subsequent multi-source data registration.

[0020] Step 102: Based on multi-source sensing data, extract geometric anomaly features at the dam inlet, low resistivity anomaly features inside the dam body, and temperature anomaly features at the dam outlet, and construct a three-dimensional spatial distribution model of anomaly features under a unified spatial coordinate system.

[0021] In this step, the data processing task is to extract features directly related to seepage from the raw signal and unify them into a common spatial framework. Geometric anomalies at the dam inlet refer to anomalous depressions, pits, or suction-like structures identified by analyzing the local geometry of underwater topographic point clouds; these are typically the starting points of seepage channels. Low-resistivity anomalies within the dam body refer to three-dimensional spatial regions with resistivity values ​​significantly lower than the surrounding background soil, obtained through resistivity inversion; these usually indicate the main body of seepage channels with high water content. Temperature anomalies at the dam outlet refer to connected regions with temperatures significantly different from the surrounding environment (background temperature), identified through thermal infrared image processing; these regions correspond to the outlet points of seepage water.

[0022] To construct a three-dimensional spatial distribution model of anomaly features, this embodiment employs a pre-defined dam engineering coordinate system, such as the dam axis as the X-axis and the vertical direction as the Z-axis. Multi-source data registration technology is used to map the spatial locations of the three types of features mentioned above into this coordinate system. For example, the coordinates of the acoustic point cloud, the grid coordinates of the resistivity inversion model, and the projected coordinates of the infrared image are transformed and aligned to form a comprehensive three-dimensional dataset containing multiple physical attributes (elevation, resistivity, and temperature). This model not only displays the geometric location of each anomaly feature but also preserves its physical attribute values, providing a data foundation for subsequent physical fusion.

[0023] Step 103: Introduce a pre-constructed dam seepage physical field model. Based on the three-dimensional anomaly feature spatial distribution model, in the three-dimensional search space, starting from the inlet geometric anomaly feature and ending with the outflow temperature anomaly feature, and based on the low resistance anomaly feature and seepage physical field constraints, search and construct a set of three-dimensional seepage candidate paths.

[0024] This step addresses the ambiguity of path identification through physical constraints. The dam seepage physical field model is a model obtained in advance through numerical simulation or analytical calculation, describing the theoretical seepage state of the dam under the current water level conditions. Specifically, it includes a three-dimensional head scalar field and a three-dimensional hydraulic gradient vector field. The seepage physical field constraints refer to the fluid dynamic laws that the seepage channel must follow. The most basic constraints include that the water flow must flow from high head to low head (energy decrease), and the flow direction should roughly follow the direction of the hydraulic gradient. The process of searching and constructing a set of three-dimensional seepage candidate paths is no longer a blind geometric connection in three-dimensional space, but a directed search guided by the physical field model. For example, the algorithm starts from the geometric center where the inlet feature is located, preferentially following the local hydraulic gradient direction, while tending to traverse low-resistivity anomaly regions (i.e., regions with high permeability), gradually extending towards the outflow feature direction, and finally forming one or more three-dimensional paths connecting the inlet and outflow points. By introducing the above physical constraints, pseudo-paths that are spatially close but hydraulically impossible to connect, such as counter-current paths or paths crossing impermeable layers, can be effectively eliminated.

[0025] Step 104: Based on the conservation principles of fluid mechanics and thermodynamics, calculate the physical consistency index of each path in the candidate path set, and select the optimal path as the leakage channel location result based on the physical consistency index.

[0026] In this step, we further introduce conservation laws across physical fields to verify the authenticity of the path. The fluid dynamics and thermodynamics conservation criteria specifically refer to mass conservation (the amount of water flowing in at the inlet should be approximately equal to the amount of water flowing out at the outlet) and energy conservation (the energy loss along the flow path should conform to physical laws). The physical consistency index is a quantitative score used to measure the degree to which the candidate path conforms to the above physical laws. Specifically, in the calculation, we can estimate the Darcy flow rate at the inlet (based on the inlet area and channel permeability) and the thermal convection flow rate at the outlet (based on the intensity of thermal infrared anomalies), and calculate the balance ratio between the two; simultaneously, we check whether the head along the path strictly decreases monotonically.

[0027] Based on the above physical indicators, each path in the candidate path set is scored or ranked, and the path with the highest score or that best conforms to the conservation criterion is identified as the true leakage channel. This method utilizes the inherent connection between Darcy flow and thermal convection, two independent physical processes, to provide strong physical evidence for the identification of leakage channels and improve the credibility of the location results.

[0028] In some optional implementations, the present invention also includes a mechanism for handling abnormal situations. When multi-source sensing data is missing or of abnormal quality, the system first performs a data validity check: if the missing area of ​​acoustic scan data exceeds 20% of the total detection area, supplementary detection is triggered or data repair is performed using an interpolation algorithm; if the resistivity detection data shows obvious poor electrode contact signals, such as negative values ​​or abnormal jumps in apparent resistivity, the data of that measurement point is removed and set as a free boundary in the inversion; if thermal infrared imaging is affected by cloud cover or rainfall, it is delayed until weather conditions improve and then re-acquired. When the path search algorithm fails to find an effective path connecting the inlet and outlet within the preset maximum number of iterations, the system automatically relaxes the search constraints (such as increasing the interlayer search radius or lowering the cost threshold) for a second search; if it is still difficult to obtain an effective path, the system outputs a prompt message "Difficult to confirm the connectivity of the leakage channel" and suggests supplementary detection or manual verification.

[0029] Example 2 details the specific process of achieving high-precision data acquisition and feature extraction through hardware synchronization and joint data processing of acoustic and electrical resistivity devices, solves the inversion artifact problem caused by water depth fluctuations and topographic undulations in traditional electrical resistivity inversion, and provides specific extraction algorithms for inlet geometric anomalies and outlet temperature anomalies.

[0030] Step 201: Use a preset clock synchronization controller to trigger the underwater acoustic scanning device and the resistivity detection device inside the dam to obtain acoustic echo data and potential difference data with a unified time reference.

[0031] In this embodiment, the data acquisition system addresses the temporal and spatial alignment issues of different detection methods. The clock synchronization controller can be a hardware trigger box based on GPS timing or a high-precision crystal oscillator, simultaneously sending trigger pulses to acoustic devices (such as multibeam sonar) and electrical resistivity devices (such as parallel electrical resistivity meters) via physical cables or wireless signals. For example, when the probe vessel travels along the survey line, the synchronization controller sends a trigger signal every 1 meter traveled; the sonar immediately emits sound waves and records the echo, and the electrical resistivity meter immediately performs a power supply and measurement. The recorded acoustic data D... sonar (t) and electrical data D elec (t) All data have a uniform UTC timestamp and differential GPS coordinates. Hardware-level synchronization ensures that the underwater topography measured by acoustic methods and the apparent resistivity profile measured by electrical methods are strictly corresponding in spatial location, eliminating spatial misalignment errors caused by changes in ship speed or equipment response delays.

[0032] In one specific embodiment, the underwater acoustic scanning device employs a multibeam sonar with an operating frequency of 200 to 400 kHz and a beam opening angle of 120 to 150 degrees. A single scan can acquire more than 512 echo beams, with a horizontal resolution better than 0.1 meters and a vertical resolution better than 0.05 meters. The resistivity detection device inside the dam uses a high-density electrical resistivity transducer with an electrode spacing of 1 to 5 meters, a power supply current of 100 to 500 mA, and a measurement accuracy better than 1%, employing conventional device types such as Winner or Schlumberger transducers. The backwater-side thermal infrared imaging device uses an uncooled infrared thermal imager with an operating wavelength of 8 to 14 micrometers, a thermal sensitivity (NETD) better than 50 milliklvin, and a spatial resolution of not less than 640 × 480 pixels. These parameters can be adjusted according to the specific project scale and detection accuracy requirements.

[0033] The clock synchronization triggering mechanism ensures instantaneous consistency of multiphysics data through an external time signal. Since dam seepage is a dynamic process, the surface geometry reflected by acoustic methods and the internal electrical properties reflected by electrical methods must correspond under the same hydraulic conditions. By unifying the time reference, the systemic bias caused by water level fluctuations in data fusion is eliminated.

[0034] Step 202: Based on acoustic echo data, reconstruct the underwater topographic point cloud and generate a water-soil interface elevation model.

[0035] In this step, acoustic data processing is used to obtain accurate underwater geometric boundaries. First, beamforming, sound velocity correction, and attitude compensation are performed on the acoustic echo data to obtain a series of discrete three-dimensional coordinate points, i.e., the underwater terrain point cloud P. bed ={(x i ,y i ,z iUsing triangular network of irregularities (TIN) interpolation or Kriging interpolation algorithms, discrete points are constructed into a continuous digital elevation model (DEM), namely a water-soil interface elevation model. This model accurately describes the interface morphology between the reservoir water and the dam slope, and is a key geometric constraint condition in subsequent electrical resistivity tomography (EDT) inversion.

[0036] Step 203: Extract geometric anomaly features of the dam entrance, specifically including: calculating the local average curvature of the underwater topographic point cloud and identifying areas with curvature lower than a preset concavity threshold as geometric anomaly features of the entrance.

[0037] This step describes the inlet identification algorithm in detail. For each point p in the underwater topographic point cloud, calculate the local mean curvature K in its neighborhood. H Specifically, the local average curvature K of each sampling point in the underwater topographic point cloud is calculated. H Identify its geometric concavity features. The calculation formula is as follows:

[0038] K H =(k1+k2) / 2;

[0039] Among them, K H Let k1 be the local average curvature, k2 be the first principal curvature, and k2 be the second principal curvature.

[0040] For leveling slope protection, K H K approaches 0; for collapsed or scourd pits, K H It is a negative value. Set the indentation threshold K. th When K H <K th It was identified as a suspected entry point.

[0041] A normally flat slope has a curvature close to zero; however, when there are pits or cavities, the local surface will show a concave shape, and the corresponding average curvature will be significantly less than zero (negative value). A concavity threshold K is set. th For example, -0.5 (1 / m), generally can be taken as -0.3 to -1.0 (1 / m). The specific value can be adjusted according to the flatness of the dam slope and the minimum size of the detection target. The flatter the slope or the smaller the detection target, the larger the absolute value of the threshold should be. H <K th Points marked as outliers are then clustered into independent regions using the Euclidean clustering algorithm. For each clustered region, its geometric center coordinates (X, Y, X) are calculated. in ,Y in Z in ) and projected area A inlet The above parameters will serve as the starting point for subsequent path search and the input for traffic estimation.

[0042] Furthermore, as a preferred embodiment, the spatial matching degree M between the center of the inlet recess and the low-resistivity anomaly directly below can also be calculated. spatial The calculation formula is:

[0043] M spatial =exp(-d center 2 / σ 2 )*(A overlap / A inlet );

[0044] Where, d center The horizontal Euclidean distance between the center of the depression and the center of the low-resistivity body is represented by σ; σ is the distance attenuation factor, used to control the weight of the distance's influence on the matching degree; A overlap Let A be the area of ​​overlap between the two on the horizontal plane. inlet This represents the projected area of ​​the indented region at the entrance. If M... spatial If the value is below a preset threshold (e.g., 0.3), the depression is considered to be a surface defect rather than a deep leakage entrance and is therefore rejected to improve the accuracy of entrance identification.

[0045] Step 204: Extract the outflow temperature anomaly features behind the dam. Specifically, this includes: performing atmospheric radiation correction and background temperature removal on the thermal infrared imaging data on the backwater side, and using a region growing algorithm to segment connected regions with temperatures higher than the background temperature as outflow temperature anomaly features.

[0046] This step details the escape point identification algorithm. Atmospheric radiation correction is applied to the original thermal infrared image to eliminate the influence of air humidity and transmission distance on temperature measurement accuracy, resulting in a true surface radiation temperature map. By statistically analyzing the temperature histogram of the image, the average temperature of most non-leaking areas is determined as the background temperature T. bg And set an abnormal threshold ΔT th For example, a threshold of 1.5 degrees Celsius can generally be set between 1.0 and 3.0 degrees Celsius. The specific value can be adjusted according to the season and diurnal temperature range. In environments with large temperature differences, a higher threshold should be used to reduce false alarms. (The image should contain values ​​where T>T.) bg +ΔT th Pixels with high confidence are marked as suspected escape points. A region growing algorithm is used to merge connected anomalous pixels into independent anomalous regions Ω, starting from high-confidence seed points. out For each region, calculate its centroid coordinates (X... out ,Y out Z out ), average temperature difference ΔT and area A out The above characteristics accurately depict the location and intensity of seepage on the back slope.

[0047] Step 205: Construct an inversion grid that includes the water body region and the dam body region, embed the water-soil interface elevation model as the known geometric boundary into the inversion grid, and set the resistivity of the water body region to the fixed value of the predicted water body resistivity.

[0048] When establishing the finite element or finite difference grid for electrical inversion, the horizontal surface is no longer assumed; instead, the water-soil interface elevation model generated in step 202 is directly used to cut the grid. All grid nodes above this interface are labeled as water regions, and all nodes below this interface are labeled as dam regions. Due to the resistivity ρ of the reservoir water... water The resistivity can be directly measured on-site using a portable conductivity meter and is relatively uniformly distributed. During the inversion process, the resistivity values ​​of all grids in the water body area are directly locked to ρ. water It does not participate in iterative updates. This approach is equivalent to adding strong prior constraints to the inversion equation, eliminating the shallow inversion blind zone of the dam body caused by the low resistivity shielding effect of the water body in traditional inversion, and allowing the inversion algorithm to focus on the changes in the internal medium of the dam body below the interface.

[0049] Step 206: Perform joint constraint inversion based on potential difference data and inversion mesh to obtain a three-dimensional resistivity distribution model, and identify low-resistivity anomaly features from the three-dimensional resistivity distribution model.

[0050] In this step, the collected potential difference data d is used obs Inversion calculations are performed using the constructed constrained mesh. The inversion aims to find an optimal dam resistivity distribution model that minimizes the residual between the potential difference calculated in the forward model and the measured data. Due to the introduction of water body constraints, the ambiguity of the inversion problem is significantly reduced. The resulting three-dimensional resistivity distribution model clearly reflects the internal electrical structure of the dam. Based on this model, a low-resistivity threshold (e.g., 0.6 times the background resistivity) is set to extract low-resistivity anomalies, which are the main potential seepage channels.

[0051] Example 3 details the mathematical implementation and processing flow of introducing tubular geometric prior constraints during the resistivity inversion process inside the dam. Addressing the technical problems of traditional smooth constraint inversion leading to volume expansion of low-resistivity anomalies, blurred boundaries, and a tendency towards spherical or blocky shapes, an objective function incorporating a tubular morphology penalty term is constructed. By embedding an anisotropic smoothing operator into the inversion algorithm, the inversion results can be forced to exhibit a slender tubular morphology in space that conforms to the characteristics of real seepage channels, improving the accuracy and spatial convergence of characterizing deep, hidden seepage channels.

[0052] Step 301: Perform joint constraint inversion based on potential difference data and inversion grid, specifically by using an objective function that includes a tubular morphology penalty term for inversion.

[0053] Unlike resistivity inversion using conventional smoothing constraints, this embodiment constructs the inversion problem as a nonlinear optimization problem to recover the true resistivity distribution of the subsurface medium from observational data. Unlike traditional methods that only include data fitting terms and conventional smoothing terms, this embodiment constructs a joint inversion objective function coupled with multiple constraints. Specifically, the objective function Φ(m) can be expressed by the following linear formula:

[0054] Φ(m)=Φ data (m)+λ b *Φ boundary (m)+λ tube *R tube (m)+λ s *Φ smooth (m);

[0055] in:

[0056] Φ(m) represents the total objective function value, and inversion is to find the model parameter m that minimizes Φ(m);

[0057] m represents the parameter vector of the three-dimensional resistivity model to be inverted;

[0058] Φ data (m) represents the data fitting term, used to measure the degree of fit between the model response and the measured data, and is usually defined as ||W d *(d obs -F(m))|| 2 , where d obs Let W be the measured potential difference data vector, F(m) be the potential difference response calculated using forward modeling, and W be the potential difference response. d Weight the data into a matrix;

[0059] Φ boundary (m) represents the water interface constraint term, which uses the geometric information of the water-soil interface to force the resistivity of the water area to remain at a fixed value;

[0060] R tube (m) represents the tubular shape penalty term in this embodiment, which is used to quantify the degree to which the low-resistivity anomaly deviates from the ideal tubular shape in the current model;

[0061] Φ smooth (m) represents the model smoothing constraint term, which is used to ensure the spatial continuity of the model parameters;

[0062] λ b , λ tube , λ s These are the regularization weight coefficients for each constraint term, used to balance the influence weight of each constraint term on the inversion result.

[0063] Furthermore, the tubular morphology penalty term Rtube (m) is constructed based on the geometric features of the leakage channel's cross-section. In each iteration, the algorithm extracts several vertical cross-sections along the main extension direction of the channel and calculates the roundness deviation of the low-resistivity region on each cross-section. Specifically, R... tube The calculation of (m) can be achieved using the following formula:

[0064] R tube (m)=Σ s [1-(4*π*A s ) / (P s 2 )] 2 ;

[0065] Among them, R tube (m) represents the tubular morphology penalty term;

[0066] Σ s This represents the summation of all sampled cross sections s extracted along the channel path;

[0067] A s This represents the area of ​​the low-resistivity abnormal region on the s-th cross-section, which can be obtained by summing the areas of the grid cells on this cross-section whose resistivity values ​​are lower than the preset low-resistivity threshold.

[0068] P s This represents the perimeter of the low-resistivity anomalous region on the s-th cross-section. Specifically, it can be obtained by extracting the edge contour length of the low-resistivity region on that cross-section. The divergence expansion of non-circular cross-sections is penalized by roundness deviation measurement.

[0069] π represents the mathematical constant pi.

[0070] Term (4*π*A) s ) / (P s 2 ) is the roundness formula in plane geometry. For an ideal circle, the value is 1; for flat or irregular shapes, the value is less than 1.

[0071] As can be seen from the above formula, the closer the inverted low-resistivity body cross-section is to a circle (i.e., exhibiting tubular characteristics), the higher R... tube When the value of (m) approaches 0, the penalty to the objective function is smaller; conversely, if the low-resistivity body exhibits a plate-like, sheet-like, or diffuse distribution, R... tube The value of (m) will increase significantly, driving the optimization algorithm to adjust the model parameter m, causing it to converge in the direction of generating tubular structures.

[0072] In some alternative implementations, to improve computational efficiency, R tubeThe calculation of (m) can employ a simplified ellipse fitting method. For example, the moment of inertia tensor of the low-resistivity region can be calculated, and the roundness deviation can be approximated by the ratio of the major axis to the minor axis of the moment of inertia. If the lengths of the major and minor axes are close, it is considered to have good tubular characteristics; if the major axis is much larger than the minor axis, a larger penalty weight is applied. Furthermore, the weighting coefficient λ... tube An adaptive adjustment strategy can be adopted, setting a small value at the beginning of the inversion to allow for a wide range of model searches, and gradually increasing λ as the number of iterations increases. tube The value of is used to strengthen geometric constraints in the later stages of inversion and achieve fine shaping of the anomalous body.

[0073] Step 302, wherein the tubular morphology penalty term includes an anisotropic smoothing operator perpendicular to the principal direction of channel extension determined by principal component analysis of the low resistivity anomaly region in the current inversion iteration result. The anisotropic smoothing operator is configured to suppress the diffusion of the resistivity anomaly region in the cross-sectional direction during the iterative inversion process, so as to enhance the spatial convergence of the low resistivity anomaly in the three-dimensional resistivity distribution model.

[0074] In this embodiment, to directly implement tubular constraints at the gradient update level, an anisotropic smoothing mechanism based on a local coordinate system is designed. First, the extension direction of the leakage channel needs to be identified in real time. Specifically, after each iteration update, the set of low-resistivity anomaly voxels in the current model m is extracted, and the first principal component direction of this point cloud dataset is calculated using Principal Component Analysis (PCA). This first principal component direction is defined as the extension direction of the channel, denoted as the unit vector n. axial A local orthogonal coordinate system (ξ, η, ζ) is established based on this principal direction, where the ξ axis is parallel to the n axis. axial The η-ζ plane is perpendicular to n. axial .

[0075] The anisotropic smoothing operator strongly suppresses diffusion perpendicular to the channel direction through local coordinate transformation. Its gradient constraint formula is as follows:

[0076] ;

[0077] Where, Φ smooth (m) represents the model smoothing constraint term, w axial The vertical smoothing weighting coefficients are... w represents the partial derivatives of the model parameters along the channel extension direction. cross The horizontal smoothing weighting coefficients are used. and These are the partial derivatives of the model parameters along the two orthogonal axes of the channel cross-section, respectively. Let w... cross For w axial More than 10 times that of the anomalous body, achieving spatial convergence.

[0078] By setting Wcross The norm is significantly greater than W. axial norm, for example W cross The element value is W axial This algorithm, which is 10 to 20 times more efficient than traditional methods, suppresses drastic changes or dispersion of model parameters in the cross-sectional direction during iteration, forcing resistivity contour lines to converge rapidly in the direction perpendicular to the channel, forming steep boundaries. Simultaneously, it allows model parameters to maintain relatively gentle and continuous changes along the channel direction. This mechanism tightens the low-resistivity anomaly that might otherwise diverge due to inversion ambiguity, reconstructing a clear, continuous leakage channel model with well-defined tubular boundaries in three-dimensional space.

[0079] In some implementations, considering that the actual seepage channels of the dam may be curved, the main direction n is extended. axial The orientation is not constant across the entire region. This embodiment can employ a block-based or voxel-based local orientation estimation strategy. Specifically, for each grid node in the model, the local characteristic orientation is calculated using the structure tensor in its neighborhood, and used as n at that location. axial Correspondingly, the anisotropic smoothing operator also rotates dynamically with the spatial position, always acting on the cross section perpendicular to the local streamlines. Through adaptive directional guidance, it can accurately adapt to the inversion of leakage channels with complex shapes such as bends and twists, ensuring that good tubular convergence effect can still be maintained at the channel bends.

[0080] For example, suppose a dam area undergoes resistivity inversion using conventional smoothing constraints. The initially obtained low-resistivity anomaly is a flattened ellipsoid, with its principal axis along the longitudinal direction of the dam body. Its major axis is approximately 8 meters, and its minor axis is approximately 2 meters, with a major-to-minor axis ratio of 4:1. After introducing the tubular morphology penalty term in this embodiment, a lateral smoothing weight W is set. cross The element value is the vertical weight W axial The cross-sectional area of ​​the low-resistivity anomaly obtained after 20 iterations of optimization was 15 times larger than that of a single anomaly. After these iterations, the major axis of the cross-section shrank to 3.5 meters, while the minor axis expanded to 3 meters, resulting in a major-to-minor axis ratio of 1.17:1, approaching a circular shape. Simultaneously, the spatial continuity along the channel's extension direction was significantly improved, with the anomaly converging from multiple discrete clumps into a continuous tubular structure approximately 45 meters long with an average equivalent diameter of about 3.2 meters. This result highly matches the actual leakage channel morphology revealed by subsequent borehole verification, validating the effectiveness of the tubular prior constraint.

[0081] Example 4 details the implementation process of the Penetrating Path Search (LPS) algorithm based on equipotential surface layering of the present invention. By reconstructing the traditional Cartesian coordinate search space into an equipotential surface layering space that conforms to the laws of fluid mechanics, and defining a specific interlayer penetration cost function, the technical problem that traditional path search algorithms are prone to generating reverse flow paths or unclear physical meanings is solved, thereby realizing three-dimensional intelligent tracking of leakage channels.

[0082] Step 401: The seepage physical field model of the dam is pre-constructed through the following steps: establishing a three-dimensional geometric model of the dam, including the dam body outline, dam foundation range, and upstream and downstream water level boundary conditions; under the steady-state seepage assumption, applying constant head boundary conditions and impermeable boundary conditions to the three-dimensional geometric model of the dam, and solving the seepage continuity equation to obtain nodal head distribution data covering the dam analysis area; calculating the spatial gradient based on the nodal head distribution data to generate a three-dimensional hydraulic gradient vector field, such as... Figure 4 As shown.

[0083] In this embodiment, the construction of the physical field model is the foundation for subsequent searches and belongs to the offline preparation stage. Specifically, firstly, based on the dam's design drawings or the latest measured topographic data, the three-dimensional geometric entity of the dam is constructed using computer-aided design software or finite element preprocessing software. When setting boundary conditions, the upstream phreatic surface is set as the first type of boundary condition (constant head boundary), with a head value H. up It equals the current upstream water level elevation; sets the downstream outflow surface as the boundary where seepage may occur; sets the dam base surface and both bank slopes as the second type of boundary condition (waterproof boundary, i.e., normal flow is 0).

[0084] Under the steady-state assumption, the water head distribution inside the dam follows the seepage continuity equation. The formula is as follows:

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

[0086] Where div is the divergence operator, K is the permeability coefficient tensor, ▽ is the gradient operator, and H is the total head scalar value to be determined. The three-dimensional head field generated based on the solution of this equation is the physical basis for subsequent equipotential surface hierarchical search.

[0087] Numerical solutions to the above equations using the finite element method (FEM) or the finite difference method (FDM) can yield the head value H for all grid nodes within the computational domain. i Based on this, by calculating the spatial negative gradient of the head scalar field H(x,y,z), the three-dimensional hydraulic gradient vector field J(x,y,z) can be obtained. At each point in space, this vector field J points in the direction of the fastest head drop, representing the direction of the seepage driving force under ideal conditions.

[0088] Step 402: The physical field model of dam seepage includes a three-dimensional head scalar field and a three-dimensional hydraulic gradient vector field covering the dam analysis area; a set of three-dimensional seepage candidate paths is searched and constructed in the three-dimensional search space, specifically in the hierarchical equipotential surface space reconstructed based on the three-dimensional head scalar field.

[0089] Traditional search algorithms typically operate directly on regular 3D meshes (Voxels) or point clouds. This isotropic spatial structure makes it difficult to embed the strong constraint that water flow can only flow downstream. This embodiment utilizes the 3D head scalar field H calculated in step 401 to discretize the continuous space into a series of equipotential surfaces. Specifically, a head interval ΔH is set, for example, 0.5 meters, starting from the upstream water level H... up Initially, an equal-head surface is extracted every ΔH, forming an ordered family of surfaces. This family of surfaces constitutes a layered equipotential surface space, whose topological structure is similar to a layered cake, with each layer representing a specific potential energy level.

[0090] Step 403: Search and construct a set of three-dimensional leakage candidate paths in the hierarchical equipotential surface space reconstructed based on the three-dimensional head scalar field, such as... Figure 2 As shown, this specifically includes: based on a three-dimensional head scalar field, extracting a set of discrete equipotential surfaces with strictly monotonically decreasing head values, and defining the set of discrete equipotential surfaces as a hierarchical node structure in the three-dimensional search space.

[0091] In this embodiment, to enable the computer to process continuous curved surfaces, they need to be discretized. Specifically, the Moving Cubes algorithm can be used to extract the k-th equipotential surface S from the three-dimensional head scalar field. k The corresponding head value is H k The equipotential surface S k It consists of a series of triangular facets or discrete points. To construct the node structure for graph search, at each equipotential surface S... k Uniform sampling is performed on the top to generate a node set N. k ={n k,1 ,n k,2 ,...n k,m' (m' is the total number). Therefore, the entire search space is reconstructed into a hierarchical structure: Layer1 (highest head), Layer2, ..., Layer... N (Minimum head). Each node n k,i Each carries its three-dimensional coordinates (x, y, z), its level index k, and the normal vector vec of that point. n (Pointing in the direction of decreasing water head).

[0092] Step 404: Map the inlet geometric anomaly features to the high-head initiation layer in the discrete equipotential surface set, and map the outflow temperature anomaly features to the low-head termination layer in the discrete equipotential surface set.

[0093] In this step, the start and end points of the search task are established. For each identified entry point geometric anomaly feature, such as the crater center P... in Calculate the water head H at its location. in And find the water head value closest to H in the hierarchical node structure. in The equipotential surface layer is denoted as the initial layer. start Find the distance P in this layer. in The nearest node is used as the starting point for the search, node n. start Similarly, for anomaly characteristics of the escape temperature, such as the geometric center P of the escape region... out Find the corresponding termination layer. end and the endpoint node n end This process anchors anomalous features in geometric space to a specific potential energy surface in physical field space.

[0094] Step 405: Perform a layer-by-layer penetration search in the hierarchical node structure. Construct candidate paths by connecting nodes between adjacent equipotential surface layers. The layer-by-layer penetration search is configured to allow only unidirectional expansion from the high head layer to the adjacent low head layer and prohibits reverse hierarchical jumps.

[0095] This step details the execution logic of the LPS algorithm. Unlike A* or Dijkstra's algorithms, which allow expansion in any direction, Layer-by-Layer Penetration Search (LPS) enforces a one-way constraint. The algorithm maintains a queue of active nodes, initially containing only the starting node n. start At each step of the search, assume the current node n curr Located at layer k, the algorithm is only allowed to probe nodes at layer k+1, and is strictly prohibited from probing nodes at layer k or k-1. This hard constraint prohibiting reverse hierarchical jumps directly solidifies the physical law of water flowing from high to low in the algorithm structure, avoiding the possibility of generating reverse flow paths (i.e., water flowing from low to high), and improving the physical interpretability of the results.

[0096] Step 406, the layer-by-layer penetration search is configured to construct candidate paths based on the criterion of minimizing the cumulative interlayer crossing cost; the interlayer crossing cost consists of multiple weighted summation sub-costs, which include at least: crossing angle cost, which is calculated based on the angle between the interlayer connection path vector and the normal vector of the equipotential surface node, used to penalize interlayer jumps that deviate from the seepage normal; and resistivity integral cost, which is calculated based on the resistivity value of the spatial region traversed by the interlayer connection path vector in the low-resistivity anomaly feature, used to reward interlayer jumps that traverse the low-resistivity region.

[0097] In this embodiment, starting from node n at layer k... i To the (k+1)th level node n j Whether the connection is reasonable depends on the cost of inter-layer traversal, C. step This is used for measurement. The design of this cost function incorporates features from both fluid dynamics and geophysics.

[0098] Specifically, the cost of crossing angle C angle The calculation formula can be expressed as:

[0099] C angle =1-(vec v *vec n ) / (|vec v |*|vec n |);

[0100] Among them, vec v =r j -r i For node n i Pointing to n j displacement vector; vec n For node n i The normal vector of the equipotential surface (i.e., the direction of the local streamline); |...| represents the magnitude of the vector; * represents the dot product.

[0101] The physical meaning of this formula is that when crossing the direction vec v vec with ideal seepage direction n When they are completely identical, the included angle is 0 degrees, the cosine value is 1, and the cost is C. angle The value is 0 (optimal); as the angle between the two increases, the cost increases rapidly.

[0102] resistivity integral cost C resist The calculation formula can be expressed as:

[0103] C resist =(ρ avg -ρ min ) / (ρ max -ρ min );

[0104] Where, ρ avg For the connection path n i to n j The average resistivity along the line segment is obtained by sampling along the line in a three-dimensional resistivity model; ρ min and ρ max These are the minimum and maximum resistivity values ​​across the entire field.

[0105] The meaning of this formula is that if the path passes through a low-resistivity area (a suspected leakage area), ρavg Smaller, cost C resist If the value is relatively small, the algorithm tends to choose that path; otherwise, it will penalize it.

[0106] Cost of inter-layer travel C step This determines the evolution trend of the path in the physical field. The formula is as follows:

[0107] C step =w1*(1-(v*n) / (|v|*|n|))+w2*((ρ avg -ρ min ) / (ρ max -ρ min ));

[0108] Among them, C step The cost of single-step inter-layer traversal is given by w1, where w1 is the directional constraint weight, v is the connection vector between adjacent layer nodes, n is the normal vector of the current layer node, |v| and |n| are the magnitudes of the corresponding vectors, w2 is the electrical constraint weight, and ρ is the directional constraint weight. avg ρ is the average resistivity of the region covered by the connecting vector. min and ρ max These are the minimum and maximum resistivity values ​​across the entire field, respectively.

[0109] For example, suppose there is a node A at layer k (head 50.0m), and two candidate nodes B and C at layer k+1 (head 49.5m). The normal vector at node A is vertically downward. Vector AB makes a 10-degree angle with the normal, and the average resistivity of the region it passes through is 20 Ω·m; vector AC makes a 5-degree angle with the normal, but the average resistivity of the region it passes through is 100 Ω·m. The minimum resistivity across the entire field is 10 Ω·m, and the maximum is 1000 Ω·m.

[0110] For path AB:

[0111] C angle (AB)=1-cos(10°)≈1-0.985=0.015C resist (AB)=(20-10) / (1000-10)≈0.010;

[0112] Total Cost AB =0.3*0.015+0.7*0.010=0.0115;

[0113] For path AC:

[0114] C angle (AC)=1-cos(5°)≈1-0.996=0.004C resist (AC)=(100-10) / (1000-10)≈0.090;

[0115] Total Cost AC =0.3*0.004+0.7*0.090=0.0642.

[0116] It can be seen that although path AC is straighter in direction (lower angular cost), the overall cost is higher because path AB passes through a more obvious low-resistivity region (lower resistivity cost). AB Much less than Cost AC The algorithm intelligently selects path AB, which is consistent with the physical fact that seepage water tends to flow through areas with high permeability (low resistance).

[0117] Example 5 details an alternative implementation for searching and constructing a set of three-dimensional seepage candidate paths. As an alternative to the LPS algorithm, it is suitable for scenarios where it is difficult to construct continuous equipotential surfaces or where rapid coarse path planning is required. It adopts a heuristic search based on regular grids (such as the A* algorithm) and introduces a physical field direction factor into the cost function to achieve approximate constraints on the physical laws of seepage within a general graph search framework.

[0118] Step 501: The physical field model of dam seepage includes a three-dimensional hydraulic gradient vector field. The search and construction of a three-dimensional seepage candidate path set in the three-dimensional search space is carried out in a three-dimensional regular grid space constructed by discretizing the dam analysis area using voxels.

[0119] In this embodiment, the search space is no longer reconstructed as a complex curved hierarchical structure, but instead directly uses a regular three-dimensional voxel grid. Specifically, the dam analysis region is discretized into uniformly sized cubic units, such as voxels of 0.5 m × 0.5 m × 0.5 m. Each voxel unit V ijk Each stores the three-dimensional coordinates (x, y, z) of its center point and the local resistivity value ρ obtained by mapping from the resistivity inversion model. ijk and the local hydraulic gradient vector J obtained from numerical simulation ijk This regular grid structure has the advantages of simple data structure and easy indexing. At this point, the three-dimensional leakage candidate path is defined as a sequence of adjacent voxel units V1, V2, ..., V... m .

[0120] Step 502: Use a heuristic search algorithm to search for the optimal node sequence that connects the inlet geometric anomaly features and the outlet temperature anomaly features in the three-dimensional regular grid space.

[0121] This step employs the classic A* (A-Star) heuristic search algorithm as its core engine. This algorithm manages the search process by maintaining an open list and a closed list. For each node n to be expanded, its evaluation function F(n) is defined as:

[0122] F(n) = G(n) + H(n);

[0123] Where: G(n) represents the actual cumulative cost of reaching node n from the starting node (corresponding to the inlet geometric anomaly feature) along the current path; H(n) represents the estimated cost of reaching the ending node (corresponding to the outflow temperature anomaly feature) from node n, which is usually estimated by multiplying the Euclidean distance or Manhattan distance by the minimum unit movement cost.

[0124] The algorithm selects the node with the smallest F(n) value from the open list each time for expansion until the endpoint is reached. Unlike the traditional A* algorithm, this embodiment introduces fluid dynamics constraints into the calculation of G(n), meaning that the transition cost between nodes is no longer just geometric distance.

[0125] Step 503, wherein the heuristic search algorithm employs a node transfer cost function that includes a hydraulic gradient direction factor. The hydraulic gradient direction factor is configured to dynamically adjust the cost value based on the consistency between the search expansion direction and the direction of the local hydraulic gradient vector, in order to penalize path expansion in the reverse hydraulic gradient direction.

[0126] In this embodiment, the single-step transfer cost Cost(n,m) for moving from the current node n to the adjacent node m is defined as a comprehensive function of geometric distance, medium properties, and physical field direction. Its specific linear calculation formula is as follows:

[0127] Cost(n,m)=D(n,m)*[w ρ *C ρ (m)+w j *C j [(n,m)];

[0128] in:

[0129] D(n,m) represents the Euclidean distance between node n and node m (usually 1 or sqrt(2) for adjacent grids).

[0130] C ρ (m) represents the resistivity cost term, used to guide the path towards the low-resistivity region. Its calculation formula can be the normalized resistivity: (ρ) m -ρ min ) / (ρ max -ρ min )(ρ m(where m is the resistivity of node m). This term approaches 0 when node m is located in the low-resistivity anomaly region;

[0131] C j (n,m) represents the hydraulic gradient direction factor, and its calculation formula is as follows:

[0132] C j (n,m)=1-cos(θ)=1-(vec v *vec j ) / (|vec v |*|vec j |);

[0133] Here, vec v vec is the direction vector of movement from n to m. j It is the local hydraulic gradient vector at node n, and θ is the angle between the two.

[0134] w ρ and w j These are the resistivity weight and the hydraulic gradient weight, respectively, which can be set to 0.5 and 0.5, for example, respectively.

[0135] To more clearly illustrate the mechanism of this hydraulic gradient direction factor, a specific numerical example is given below:

[0136] Suppose that at a certain node n, the hydraulic gradient vector vec j Pointing due east (i.e., the water flows eastward).

[0137] Scenario A (Downstream Search): The algorithm attempts to reach the neighbor node m in the due east direction. east Expand. The direction of movement at this point is vec. v With vec j The included angle θ is 0 degrees, cos(θ) = 1, C j =1-1=0. This means that the cost of moving downstream in the physical field direction is 0 (no penalty).

[0138] Scenario B (Crossflow Search): The algorithm attempts to reach neighbor node m in the due north direction. north Extended. At this point, the included angle θ is 90 degrees, cos(θ) = 0, C j =1-0=1. This indicates that movement laterally across a streamline will be subject to a moderate penalty.

[0139] Scenario C (Reverse Search): The algorithm attempts to find the neighboring node m to the due west direction. west Extended. At this point, the included angle θ is 180 degrees, cos(θ) = -1, C j =1-(-1)=2. This means that moving against the current will incur the greatest penalty (the cost is twice or more of the base value).

[0140] By dynamically adjusting the cost value, the search algorithm can perceive the resistance of the water flow even in a regular grid space. When there is a path flowing downstream through a low-resistance zone, its cumulative cost G(n) will be much smaller than that of a path flowing upstream or meandering. Although this method does not forcibly prohibit upstream flow through topology like the LPS algorithm, it can mathematically screen out leakage channels that conform to physical laws with a relatively high probability by imposing a high cost (soft constraint) on upstream flow, serving as a computationally simple and adaptable alternative.

[0141] Example 6 details how, after generating a set of three-dimensional seepage candidate paths, the authenticity of the paths is verified through multi-physics consistency scoring and flow balance verification, and how discrete paths are smoothed and deduplicated to finally output a high-precision three-dimensional seepage channel model. By introducing a Darcy flow-thermal convection cross-physics coupling verification mechanism, the problem of strong multiple solutions and high false alarm rate of single geophysical methods is solved, providing quantitative confidence basis for engineering decisions.

[0142] Step 601: The physical consistency index includes at least a hydraulic gradient consistency score and a head monotonically decreasing score. The hydraulic gradient consistency score is calculated based on the cosine of the angle between the tangent direction vector at each point along the candidate path and the local hydraulic gradient vector at that point in the dam seepage physical field model. It is used to characterize the degree of agreement between the path direction and the direction of the seepage driving force. The head monotonically decreasing score is calculated based on the sequence of head values ​​at each point along the candidate path in the dam seepage physical field model. It is used to penalize paths with reverse head increase phenomena.

[0143] In this embodiment, the evaluation of candidate paths no longer focuses solely on whether they traverse low-resistivity regions, but incorporates dynamic indices from fluid dynamics. For each candidate path P containing N nodes, its hydraulic gradient consistency score S is first calculated. j This score is used to quantify whether the pathway follows the natural flow of groundwater.

[0144] Specifically, for the i-th node on the path, calculate its path tangent vector t. i The hydraulic gradient vector J of this point and its neighborhood i The included angle θ i S j The calculation formula can be expressed as:

[0145] S j =(1 / N)*Σ i [cos(θ i Alternatively, in a more preferred embodiment, a Gaussian weighted form is used to strengthen the penalty for large angular deviations:

[0146] Sj =(1 / N)*Σ i [exp(-(θ i / θ ref ) 2 )];

[0147] Where N is the number of evaluation points; θ ref The preset angle tolerance parameter (e.g., 30 degrees or π / 6 radians).

[0148] If S j If the value is close to 1, it indicates that the entire path flows downstream; if S j A lower value indicates that the path contains a large number of cross-flow or even counter-flow segments, raising doubts about its physical authenticity.

[0149] Simultaneously, calculate the monotonically decreasing head score S. H Although the LPS algorithm enforces hierarchical monotonicity during construction, head anomalies can still occur in local micropaths (especially after interpolation or smoothing). The sequence of values ​​H1, H2, ..., H at each point along the path in the head scalar field is extracted. n The statistics show that H appears in the data. i+1 >H i The number of events N (i.e., reverse flow) reverse S H It can be defined as:

[0150] S H =1-(N reverse / (N-1)).

[0151] This indicator is a hard veto indicator; for genuine natural seepage channels, S H It should be strictly equal to 1 or very close to 1.

[0152] Step 602: Calculate the physical consistency index based on the conservation principles of fluid mechanics and thermodynamics, such as... Figure 3 As shown, this specifically includes applying the mass conservation principle to verify flow balance:

[0153] Based on the effective flow area of ​​the geometric anomaly characteristics at the dam inlet and the equivalent permeability coefficient mapped by the average resistivity of the candidate path in the low-resistivity anomaly characteristics inside the dam body, the estimated value of the inlet seepage flow is calculated using Darcy's law.

[0154] Based on the heat dissipation area of ​​the outflow temperature anomaly characteristics behind the dam and the intensity of the temperature anomaly relative to the ambient temperature, the estimated value of the outflow rate is calculated using the convective heat transfer formula.

[0155] Calculate the flow balance ratio between the estimated inlet inflow and the estimated outflow, and use the degree to which the flow balance ratio is close to 1 as the evaluation criterion for selecting the optimal path.

[0156] This step and subsequent steps constitute the cross-field verification stage of this invention. A bridge needs to be established between resistivity and hydraulic parameters. Although resistivity is not directly equivalent to permeability coefficient, there is a significant negative correlation between the two in water-rich dam media. This embodiment uses a preset empirical power-law formula to calculate the equivalent permeability coefficient K. eq :

[0157] K eq =α*(ρ avg ) (-β) ;

[0158] Among them, K eq ρ is the equivalent permeability coefficient, α is the base value of the permeability conversion coefficient, and ρ is the base value of the permeability conversion coefficient. avg β is the average resistivity of the channel region, and β is the permeability conversion index, which is usually taken as 1.0 to 2.0.

[0159] For different types of dam media, typical parameter ranges are as follows: for clay core dams, α ranges from 1e-5 to 1e-3, and β ranges from 1.2 to 1.8; for homogeneous earth dams, α ranges from 1e-4 to 1e-2, and β ranges from 1.0 to 1.5; for gravel dams, α ranges from 1e-3 to 1e-1, and β ranges from 0.8 to 1.2. For example, for clay core dams, α might be 1e-4, and β might be 1.5.

[0160] The estimated inlet seepage flow rate Q was calculated based on Darcy's Law. in :

[0161] Q in =K eq *A in *i avg ;

[0162] Among them, A in i represents the effective overcurrent area (in square meters) obtained from the inlet geometric anomaly features based on point cloud computing; avg K represents the hydraulic gradient modulus along the path at the inlet, which is directly extracted from the hydraulic gradient vector field generated in step 401. eq The equivalent permeability coefficient is obtained by converting the channel average resistivity through a preset empirical relationship between resistivity and permeability.

[0163] For example, suppose a candidate path passes through a region with an average resistivity ρ avg =50Ω·m, empirical parameters α=0.01, β=1.0, then K eq =0.01*50 (-1) =0.0002m / s. If the area of ​​the crater measured from the entrance point cloud is A in =2.0 square meters, local hydraulic gradient i at the inlet avg=0.1. Therefore, the theoretical seepage flow rate Q at the inlet is... in The estimate is:

[0164] Q in =0.0002 * 2.0 * 0.1 = 0.00004m 3 / s=40ml / s.

[0165] Step 603, calculating the estimated egress flow rate using the convection heat transfer formula, is based on the following formula:

[0166] Q out =(h c *A out *δT) / (ρ w *c w *ΔT w );

[0167] Among them, Q out h is the estimated outflow rate. c A is the surface convective heat transfer coefficient. out ρ represents the area of ​​the escape temperature anomaly region, δT represents the temperature difference intensity of the escape region relative to the background environment, and ρ represents the area of ​​the escape temperature anomaly region. w c is the density of water. w Let ΔT be the specific heat capacity of water. w This refers to the temperature difference between the water temperature at the leakage source and the ambient air temperature.

[0168] In this step, the flow rate is calculated using thermodynamic principles. When seepage water flows out and reaches the dam surface, it changes the local surface temperature through convective heat transfer. This process follows Newton's law of cooling and the law of conservation of energy. The specific meanings of the parameters in the formula are as follows:

[0169] h c Surface convective heat transfer coefficient (W / (m²)) 2 ·K), this value is greatly affected by wind speed, and is usually selected based on empirical values ​​of on-site wind speed observation data (e.g., 5-10 under windless conditions, and 15-30 under windy conditions).

[0170] A out Area (m²) of the escape temperature anomaly 2 (), derived from infrared image segmentation results.

[0171] δT: The difference (K or °C) between the average temperature of the abnormal area and the background temperature on the infrared image.

[0172] ρ w The density of water is usually taken as 1000 kg / m³. 3 .

[0173] c w The specific heat capacity of water is usually taken as 4200 J / (kg·K).

[0174] ΔT w The difference between the temperature of the leak source water (usually approximating the reservoir water temperature) and the ambient air temperature.

[0175] This formula represents the sensible heat flux (Q) carried out by the outflow. out *ρ w *c w *ΔT w It should be equal to the heat flux lost through convection at the Earth's surface (h). c *A out *δT).

[0176] Continuing from the previous example, suppose a thermal anomaly is found on the back slope, with an area of ​​A. out =5.0 square meters, abnormal intensity δT=2.0℃. On-site environmental parameters are: heat transfer coefficient h c =10W / (m 2 Given that the reservoir water temperature is 15℃ and the air temperature is 5℃, the temperature difference ΔT is... w =10℃. Substitute into the formula to calculate:

[0177] Q out =(10*5.0*2.0) / (1000*4200*10)=100 / 42,000,000≈0.00000238m 3 / s≈2.38ml / s.

[0178] Note: The values ​​in this example are for demonstration purposes only. In actual operating conditions, Q... in With Q out They should be on the same order of magnitude.

[0179] Step 604: Calculate the flow balance ratio between the estimated inlet inflow and the estimated outflow, and use the degree to which the flow balance ratio is close to 1 as the evaluation criterion for selecting the optimal path.

[0180] Finally, the traffic balancing ratio is used to determine the physical consistency of the path. The formula is as follows:

[0181] R Q =|ln(Q in / Q out )|;

[0182] Among them, R Q Q is the flow balance consistency index, ln is the natural logarithm function, and Q is the flow balance consistency index. in Q is the estimated inlet infiltration flow rate. out This is an estimated outflow rate. R Q The closer it is to 0, the higher the degree of mass conservation and the more realistic the path.

[0183] Step 605: Use cubic B-spline curves to perform three-dimensional spatial fitting and smoothing on the discrete node sequence contained in the candidate path to generate a centerline model of the leakage channel with continuous curvature; and construct a tubular three-dimensional channel entity along the centerline model based on the equivalent radius of the low-resistivity anomaly characteristics of the area through which the candidate path passes.

[0184] This step involves visualizing and refining the model. The original path output by the LPS search or grid search is a series of discrete points P1, P2, ..., P connected by a polyline. n In engineering demonstrations, this approach appears rigid and does not conform to the streamlined characteristics of water flow. A smooth curve C(u) with continuous second derivatives is generated using a cubic B-spline interpolation algorithm, with discrete points as control points.

[0185] Furthermore, for each point on the curve, its corresponding low-resistivity anomaly range in the resistivity inversion model is queried, and an equivalent radius r(u) is estimated, for example, by dividing the cross-sectional area of ​​the low-resistivity region by π and taking the square root. A pipe surface with radius r(u) is generated by scanning along the centerline C(u), ultimately forming a three-dimensional tubular channel solid model with realistic spatial morphology and varying thickness, which is directly output to engineers for grouting and sealing design.

[0186] Step 606, constructing a set of three-dimensional leakage candidate paths also includes deduplication and merging of multiple paths: calculating the Hausdorff distance between any two candidate paths; when the Hausdorff distance is less than a preset channel diameter threshold, the two paths are determined to be overlapping paths, and the one with the higher physical consistency index is retained as the valid path.

[0187] This step addresses the redundancy issue of the algorithm potentially generating multiple similar paths. Hausdorff distance is a measure of similarity between two sets of points, defined as the maximum of the nearest distances from any point in set A to set B. For two paths P... a and P b If their Hausdorff distance H(P) a ,P b If the diameter of the two paths is less than a preset channel diameter threshold (e.g., 1.0 meter), then physically, these two paths are considered to actually indicate the same leakage channel (high spatial overlap). In this case, the overall score for physical consistency between the two paths is compared:

[0188] Score=w j *S j +w H *S H +w Q *(1 / |1-R Q |);

[0189] Among them, w j w is the hydraulic gradient consistency weight. H For the decreasing consistency weight of the water head, w Q Consistency weights are used to balance traffic.

[0190] The system retains the highest-scoring entry and removes the lowest-scoring one. Through this deduplication mechanism, the final output will be a limited list of independent leakage channels with the highest physical confidence.

[0191] Hausdorff distance is used to measure the distance between two three-dimensional paths P. a With P b Spatial overlap. The formula is as follows:

[0192] H(P a ,P b )=max(d(P a ,P b ),d(P b ,P a ));

[0193] Among them, H(P) a ,P b ) represents the Hausdorff distance between the two paths, max is the function to maximize the distance, and d(P) a ,P b ) represents set P a All points to set P b The maximum value of the shortest distance. If this value is less than the preset channel radius, path merging is performed.

[0194] According to one aspect of this application, the path probability assessment based on the physical consistency of seepage can also be as follows:

[0195] For each candidate path, the degree of consistency with the physical laws of seepage is calculated, and the probability value is used to characterize the likelihood that the path is a real leakage channel.

[0196] Calculate the hydraulic gradient consistency score.

[0197] For candidate path L k Calculate the consistency score of the hydraulic gradient along the path. Discretize the path into N evaluation points {P1, P2, ..., P...}. n} Calculate the path tangent direction t at each evaluation point. i With respect to the hydraulic gradient direction J at that point i The included angle θ i :

[0198] θ i =arccos[(t i ·Ji ) / (|t i |×|J i |)];

[0199] Where: t i J is the unit vector of the tangent line to the path at the i-th evaluation point; i Let be the hydraulic gradient vector at the i-th evaluation point; · represents the vector dot product; || represents the vector magnitude.

[0200] Hydraulic gradient consistency score S j Defined as:

[0201] S j =(1 / N)×Σ i=1 N exp[-(θ i / θ ref ) 2 ];

[0202] Where N is the number of evaluation points; θ ref The reference deflection angle is π / 6 (30°); exp is the natural exponential function. When the path direction is consistent with the hydraulic gradient direction (θ... i When the path direction is perpendicular to or opposite to the hydraulic gradient direction, the score approaches 0.

[0203] Calculate the head decrease consistency score (head monotonically decreasing score).

[0204] In a true leakage path, the head must decrease monotonically along the seepage direction. For N evaluation points of the candidate path, extract the head value {H1, H2, ..., H...} at each point. n} Calculate the head decrease consistency score S H :

[0205] S H =N dec / (N-1);

[0206] Where, N dec The number of times the water head decreases between adjacent assessment points, i.e., satisfying H i >H i+1 The number of point pairs; N-1 is the total number of adjacent point pairs. If the head decreases strictly monotonically along the flow path, then S H =1; if there is a reverse increase in water head, then S H <1.

[0207] Calculate the consistency score for permeability anomalies.

[0208] The permeability coefficient of the seepage channel area should be significantly higher than that of the surrounding normal soil. Based on the resistivity inversion results, the relative permeability coefficient at each point along the path is estimated using the empirical relationship between soil resistivity and permeability. Let the permeability coefficient of the normal embankment soil be K0, and the estimated permeability coefficient at the i-th evaluation point along the path be K. i Then the consistency score of permeability anomaly S k Defined as:

[0209] S k =(1 / N)×Σ i=1 N [1-exp(-K i / K0)];

[0210] Where K0 is the baseline value of the permeability coefficient of normal soil, in m / s; K i Let be the estimated permeability coefficient for the i-th evaluation point, in m / s. When the permeability coefficient of the path area is significantly higher than the normal value, the score approaches 1; when the permeability coefficient is comparable to the normal value, the score approaches 0.63.

[0211] The candidate path P is calculated by weighting and combining the above three scores. k The probability value of the actual leakage path:

[0212] Prob(P k )=w j ×S j +w H ×S H +w k ×S k ;

[0213] Among them, Prob(P k ) is the candidate path P k The overall probability value, ranging from [0,1]; w j For hydraulic gradient consistency weights; w H For decreasing head consistency weights; w k The weights represent the consistency of permeability anomalies; the three weights satisfy w j +w H +w k =1.

[0214] The recommended weight value is w. j =0.4, w H =0.3, w k =0.3. The basis for this weight allocation is that the consistency of the hydraulic gradient directly reflects the physical rationality of the seepage direction, and has the highest weight; the decrease in hydraulic head and the anomaly of permeability are auxiliary verification conditions, and have equal weight.

[0215] This application employs a hardware-level clock synchronization triggering mechanism and utilizes the water-soil interface extracted by high-precision sonar as the hard boundary for electrical resistivity tomography (EDT) inversion, achieving mechanism-level registration and fusion of data from different dimensions in the same spatial coordinate system. This solves the problem of unifying multi-source heterogeneous data.

[0216] The proposed solution abandons the traditional blind geometric connection method and proposes a penetration-based search algorithm (LPS) based on equipotential surfaces. This approach reconstructs the search space into a hierarchical structure that conforms to the energy gradient distribution. By imposing a hard physical constraint that only allows unidirectional crossings from high-head to low-head water, it eliminates potential backflow or crossflow artifacts in the path. This solves the problem of the lack of physical rationality in the leakage path.

[0217] The scheme introduces a Darcy flow-thermal convection flow balance verification mechanism across physical fields. By calculating the consistency between inlet seepage flow and outlet heat loss, the path is checked in a closed-loop manner from the perspective of mass conservation. This solves the problems of high false alarm rate and insufficient confidence.

[0218] The shift from pure geometric superposition to physical constraints makes the positioning results not only spatially connected but also physically consistent, improving the accuracy and reliability of detecting hidden channels in complex seepage environments.

[0219] 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 locating seepage channels in dam engineering based on intelligent sensing, characterized in that, include: Acquire underwater acoustic scanning data of the dam area, resistivity detection data inside the dam, and thermal infrared imaging data of the backwater side; Based on multi-source sensing data, geometric anomaly features at the dam entrance, low resistivity anomaly features inside the dam body, and temperature anomaly features at the dam exit are extracted respectively, and a three-dimensional spatial distribution model of anomaly features under a unified spatial coordinate is constructed. A pre-constructed physical field model of dam seepage is introduced. Based on the three-dimensional anomaly feature spatial distribution model, in the three-dimensional search space, starting from the inlet geometric anomaly feature and ending from the outflow temperature anomaly feature, a set of three-dimensional leakage candidate paths is searched and constructed based on the low resistance anomaly feature and seepage physical field constraints. Based on the conservation principles of fluid mechanics and thermodynamics, the physical consistency index of each path in the candidate path set is calculated, and the optimal path is selected as the result of leakage channel location. The physical field model of dam seepage includes a three-dimensional scalar field of hydraulic head and a three-dimensional hydraulic gradient vector field covering the analysis area of ​​the dam; The search and construction of a set of three-dimensional leakage candidate paths is carried out in a three-dimensional search space, specifically in a hierarchical equipotential surface space reconstructed based on a three-dimensional hydraulic head scalar field. A search is conducted in the hierarchical equipotential surface space reconstructed based on a 3D head scalar field to construct a set of 3D leakage candidate paths, including: Based on the three-dimensional head scalar field, a set of discrete equipotential surfaces with monotonically decreasing head values ​​is extracted and defined as a hierarchical node structure in the three-dimensional search space. The inlet geometric anomaly features are mapped to the high-head initiation layer in the discrete equipotential surface set, and the outflow temperature anomaly features are mapped to the low-head termination layer in the discrete equipotential surface set. In the hierarchical node structure, a layer-by-layer traversal search is performed, and candidate paths are constructed by connecting nodes between adjacent equipotential surface layers. Among them, the layer-by-layer penetration search is configured to only allow unidirectional expansion from the high head layer to the adjacent low head layer, and prohibits reverse hierarchical jumps; Calculation of physical consistency indices based on fluid mechanics and thermodynamics conservation principles, including flow balance verification using mass conservation principles: Based on the effective flow area of ​​the geometric anomaly characteristics at the dam inlet and the equivalent permeability coefficient mapped by the average resistivity of the candidate path in the low-resistivity anomaly characteristics inside the dam body, the estimated value of the inlet seepage flow is calculated using Darcy's law. Based on the heat dissipation area of ​​the outflow temperature anomaly characteristics behind the dam and the intensity of the temperature anomaly relative to the ambient temperature, the estimated value of the outflow rate is calculated using the convective heat transfer formula. Calculate the flow balance ratio between the estimated inlet inflow and the estimated outflow, and use the degree to which it approaches 1 as the evaluation criterion for selecting the optimal path.

2. The method according to claim 1, characterized in that, Acquire underwater acoustic scanning data of the dam area, resistivity detection data of the dam interior, and thermal infrared imaging data of the backwater side, including: The underwater acoustic scanning device and the resistivity detection device inside the dam are triggered by a clock synchronization controller to obtain acoustic echo data and potential difference data with a unified time reference. Geometric anomaly features at the dam inlet, low resistivity anomaly features inside the dam body, and temperature anomaly features at the dam outlet were extracted, including: Based on acoustic echo data, underwater topographic point cloud is reconstructed to generate a water-soil interface elevation model. An inversion grid containing the water body region and the dam body region is constructed. The water-soil interface elevation model is embedded into the inversion grid as a known geometric boundary. The resistivity of the water body region is set to the fixed value of the predicted water body resistivity. A three-dimensional resistivity distribution model is obtained by performing joint constraint inversion based on potential difference data and inversion mesh, and low resistivity anomaly characteristics are identified from it.

3. The method according to claim 2, characterized in that, Joint constraint inversion is performed based on potential difference data and inversion grid, specifically by using an objective function that includes a tubular morphology penalty term for inversion. The tubular morphology penalty term includes an anisotropic smoothing operator perpendicular to the principal direction of channel extension determined by principal component analysis of the low-resistivity anomaly region in the current inversion iteration results. It is configured to suppress the diffusion of the resistivity anomaly region in the cross-sectional direction during the iterative inversion process, so as to enhance the spatial convergence of the low-resistivity anomaly in the three-dimensional resistivity distribution model.

4. The method according to claim 1, characterized in that, Layer-by-layer traversal search is configured to construct candidate paths based on the criterion of minimizing the cumulative inter-layer traversal cost; The cost of inter-layer traversal consists of multiple cost components calculated using a weighted summation method. These cost components include at least the following: Crossing angle cost, which is calculated based on the angle between the interlayer connection path vector and the normal vector of the equipotential surface node, is used to penalize interlayer jumps that deviate from the seepage normal. The resistivity integral cost, calculated based on the resistivity value of the spatial region traversed by the interlayer connection path vector in the low-resistivity anomaly feature, is used to reward interlayer transitions that traverse the low-resistivity region.

5. The method according to claim 1, characterized in that, The physical field model of dam seepage includes a three-dimensional hydraulic gradient vector field. A set of three-dimensional seepage candidate paths is searched and constructed within a three-dimensional search space, which is performed in a three-dimensional regular grid space; including: A heuristic search algorithm is used to search for the optimal node sequence connecting the inlet geometric anomaly feature and the outlet temperature anomaly feature in a three-dimensional regular grid space; Among them, the heuristic search algorithm adopts a node transfer cost function that includes a hydraulic gradient direction factor. The hydraulic gradient direction factor is configured to dynamically adjust the cost value according to the consistency between the search expansion direction and the direction of the local hydraulic gradient vector, so as to penalize path expansion in the opposite hydraulic gradient direction.

6. The method according to claim 1, characterized in that, Physical consistency indicators should include at least the hydraulic gradient consistency score and the head monotonically decreasing score; The hydraulic gradient consistency score is calculated based on the cosine of the angle between the tangent direction vector at each point along the candidate path and the local hydraulic gradient vector at that point in the dam seepage physical field model. It is used to characterize the degree of agreement between the path direction and the direction of the seepage driving force. The monotonically decreasing head score is calculated based on the sequence of head values ​​at each point along the candidate path in the seepage physical field model of the dam, and is used to penalize paths with reverse head increase.

7. The method according to claim 1, characterized in that, Inlet seepage flow rate estimate Q in =K eq *A in *i avg ; Among them, A in i represents the effective inlet flow area calculated based on point cloud geometry. avg Let K be the hydraulic gradient modulus along the path at the inlet. eq The equivalent permeability coefficient is obtained by converting the channel average resistivity through a preset empirical relationship between resistivity and permeability. The estimated escape flow rate is calculated using the convective heat transfer formula, based on the formula Q. out =(h c ×A out ×δT) / (ρ w ×c w ×ΔT w The process was carried out by ) Among them, Q out h is the estimated outflow rate. c A is the surface convective heat transfer coefficient. out ρ represents the area of ​​the escape temperature anomaly region, δT represents the temperature difference intensity of the escape region relative to the background environment, and ρ represents the area of ​​the escape temperature anomaly region. w c is the density of water. w Let ΔT be the specific heat capacity of water. w This refers to the temperature difference between the water temperature at the leakage source and the ambient air temperature.

Images