A Termite Nest Irrigation Planning Method Based on Ground Penetrating Radar and Physical Flow Characteristics

Abstract

This invention constructs a termite nest infill planning method based on ground-penetrating radar and physical flow characteristics. The method includes: acquiring ground-penetrating radar detection data of the target area, extracting normalized multi-attribute features, and classifying anomalies based on a feature fingerprint database to obtain an anomaly spatial label dataset; combining preset biological morphological prior constraints, performing topological association search and self-consistency verification and correction in the dataset to reconstruct the three-dimensional structural connectivity graph of the nest; determining the directed reachability of paths in the connectivity graph based on physical flow characteristics to construct a directed reachable flow graph; selecting an infill set according to the maximum coverage principle, calculating the allocation amount, and arranging the execution sequence to obtain a multi-point infill planning scheme. This invention solves the problems of low accuracy in complex nest reconstruction and lack of physical basis for infill schemes, improves the three-dimensional reconstruction accuracy of underground nests and the success rate of infill filling, and can be applied to scenarios such as liquid metal infill for nest specimen preparation and grouting sealing.

Description

Technical Field

[0001] This invention relates to the field of underground biological nest detection and non-destructive reconstruction of three-dimensional structures, and in particular to a termite nest irrigation planning method based on ground-penetrating radar and physical flow characteristics. Background Technology

[0002] Termites and other underground nesting organisms pose a serious safety threat to water conservancy projects, dams, forestry resources, and underground cables. Because their nest networks are deeply buried underground and possess highly complex spatial topologies, accurately identifying the three-dimensional connectivity of the nests and scientifically planning injection schemes (such as liquid metal injection to obtain nest samples or grouting for sealing) is a prerequisite for termite control and ensuring project safety. The combination of non-destructive testing and precise injection not only avoids structural damage caused by large-scale excavation but also provides a scientific basis for nest morphology research and engineering restoration.

[0003] Currently, ground-penetrating radar (GPR) is widely used for detecting shallow subsurface anomalies due to its high resolution and non-destructive characteristics. Existing nest detection technologies mainly rely on manual interpretation of two-dimensional radar profiles or conventional image connectivity component analysis (VPC) methods to extract subsurface cavities. After obtaining the location of the subsurface cavity, a single borehole is typically used, relying on gravity-driven natural permeation for the injection of reagents or grout. When faced with irregularly shaped and tortuous micro-biological channels, existing three-dimensional reconstruction methods are often susceptible to interference from soil heterogeneity and GPR signals, resulting in numerous breaks or false connections in the reconstructed subsurface network that do not conform to biological common sense.

[0004] Furthermore, most existing grouting strategies only consider the geometric volume and gravitational connectivity of underground channels, neglecting the loss of fluidity due to changes in physical properties as the fluid flows through a network of microchannels. For example, in high-temperature liquid metal grouting, rapid heat dissipation of the fluid within the microchannels can cause solidification and narrowing, eventually blocking the channels. Similarly, in cement-based grouting, the initial setting of the grout during hydration also leads to a decrease in fluidity over time. In practical operations, problems frequently arise such as premature loss of fluidity leading to channel blockage, airlock caused by simultaneous grouting of multiple channels, or uneven filling due to differences in flow resistance.

[0005] Therefore, it is necessary to study a spatiotemporal planning method that can improve the accuracy of reconstructing complex underground biological channels and scientifically guide multi-point injection operations. Summary of the Invention

[0006] Purpose of the invention: To provide a termite nest irrigation planning method based on ground-penetrating radar and physical flow characteristics, in order to solve the above-mentioned technical problems existing in the prior art.

[0007] Technical solution: A termite nest inundation planning method based on ground-penetrating radar and physical flow characteristics, comprising:

[0008] The ground-penetrating radar detection data of the target area is acquired, normalized multi-attribute features are extracted, and anomalies are classified based on a pre-built feature fingerprint database to obtain an anomaly spatial label dataset.

[0009] Combining pre-set biological morphological prior constraints, topological association search and self-consistency verification correction are performed in the anomaly spatial label dataset to reconstruct the three-dimensional structure connectivity graph of the nest.

[0010] Based on the physical flow characteristics, the directed reachability of paths in the three-dimensional structure connectivity graph of the nest is determined, and a directed reachable flow graph is constructed.

[0011] Based on the directed reachable flow graph, the set of injection ports is selected according to the maximum coverage principle, the allocation volume is calculated and the execution sequence is arranged to obtain a multi-point injection planning scheme.

[0012] In one exemplary embodiment, extracting normalized multi-attribute features includes:

[0013] The anomalous body region is identified from the ground-penetrating radar detection data, and an exclusion buffer zone is delineated outward from the anomalous body region as a reference.

[0014] An annular candidate window is defined in the outer region of the exclusion buffer zone, and background signal samples are extracted within the annular candidate window;

[0015] Calculate the attribute variation coefficient of the background signal sample, and perform a homogeneity test on the annular candidate window in combination with a preset homogeneity threshold;

[0016] If the homogeneity test fails, the range of the annular candidate window is expanded outward, and sampling and testing are repeated until the homogeneity test is passed.

[0017] Background reference values ​​are calculated based on background signal samples within the annular candidate window that have passed the homogeneity test.

[0018] By using background reference values ​​to perform benchmark conversion on the original attribute values ​​of the anomaly region, normalized multi-attribute features are obtained.

[0019] In one exemplary embodiment, performing a topological association search specifically includes:

[0020] Determine the primary nest starting node from the anomaly spatial label dataset;

[0021] Based on the starting node of the main nest, and combined with the path angle deflection constraints and cumulative length constraints specified in the biological morphology prior constraints, nodes are searched and associated in the anomaly spatial label dataset to construct the initial topological skeleton.

[0022] In one exemplary embodiment, searching and associating nodes specifically employs a multi-path tracing search mechanism, including:

[0023] Construct a path seed with the main nest starting node and assign it an initial direction of movement, then place it into a priority queue sorted by confidence.

[0024] Extract active paths from the priority queue and retrieve candidate nodes within the neighborhood of their end nodes;

[0025] For candidate nodes, perform local direction deflection verification and cumulative length verification. Add candidate nodes that meet the path angle deflection limit and cumulative length limit as new branches to the current path, update the local forward direction, and put them back into the priority queue.

[0026] Complete paths are extracted based on preset path termination conditions and node affiliation conflict resolution rules to construct the initial topology skeleton.

[0027] In an exemplary embodiment, the directed reachability of paths in the three-dimensional connectivity graph of the nest is determined based on physical flow characteristics, and a directed reachable flow graph is constructed, including:

[0028] Extract the path segments from the 3D connected graph of the nest structure;

[0029] Based on the static properties, the gravitational potential difference and capillary resistance of each path segment are evaluated to determine the candidate segments that meet the static reachability conditions.

[0030] Based on solidification kinetics, time-varying gating constraint verification is performed on candidate segments to verify whether the fluid can traverse the entire length of the segment before it completely solidifies and blocks it.

[0031] The path segments that simultaneously satisfy both static reachability conditions and time-varying gating constraint verification are used as reachable directions to construct a directed reachable flow graph.

[0032] In an exemplary embodiment, performing time-varying gating constraint verification on candidate segments based on solidification kinetics characteristics includes:

[0033] Based on the thermophysical properties of the wall material corresponding to the candidate segments, the solidification rate coefficient characterizing the fluid's solidification migration from the pipe wall to the center is obtained;

[0034] Based on the solidification rate coefficient, the equivalent diameter of the cross section of the candidate segment, and the minimum channel threshold required for fluid flow, the flow duration corresponding to the candidate segment is calculated.

[0035] By combining the inclination angle of the candidate segments with the time-varying reduction characteristics of the fluid velocity, the thermally permeable distance that the fluid can reach during the flow duration is estimated.

[0036] Compare the actual segment length of the candidate segment with the thermally passable distance. If the actual segment length does not exceed the thermally passable distance, the candidate segment is determined to have passed the time-varying gating constraint check.

[0037] In one exemplary embodiment, calculating the allocation amount specifically includes:

[0038] For exclusive nodes that can be reached independently from a single injection port in a directed reachable flow graph, their entire cavity volume is included in the total allocation of the corresponding injection port.

[0039] For a shared node that is simultaneously reached by multiple injection ports through different paths, each path is treated as a series pipe and the flow resistance corresponding to each path is calculated separately.

[0040] The volume allocation coefficient is determined based on the reciprocal proportion of the flow resistance corresponding to each path, and the cavity volume of the shared node is allocated to the corresponding injection port according to the volume allocation coefficient, so as to obtain the allocation amount of each injection port.

[0041] In one exemplary embodiment, the execution timing is arranged as follows:

[0042] Identify adjacent injection ports with shared path segments based on directed reachable flow graphs;

[0043] Obtain the longest permissible interval between the fluid injected from the preceding inlet solidifying into a blockage in the shared path segment;

[0044] Obtain the shortest permissible interval time required for the fluid injected from the preceding infusion port to flow to the far end of the shared path segment;

[0045] By combining the longest and shortest allowed interval times, execution timing with time interval window constraints is arranged for adjacent injection ports.

[0046] In an exemplary embodiment, the pre-built feature fingerprint library is obtained through the following steps:

[0047] Obtain measured dielectric data of material samples under different moisture content gradients;

[0048] A multi-component mixed dielectric model is constructed to characterize the arrangement of subphases within the material. The equivalent properties of the multi-component mixed dielectric model are jointly determined by the volume fraction and shape factor of each subphase.

[0049] Using the measured dielectric data of material samples as the inversion target, the error minimization solution is performed on the multi-component mixed dielectric model to invert and calibrate the optimal value of the shape factor;

[0050] Forward modeling of various preset nest configurations is performed using a calibrated multi-component mixed dielectric model, and response features are extracted to compile a pre-built feature fingerprint library.

[0051] According to another aspect of the embodiments of this application, a computer-readable storage medium is also provided, wherein a computer program is stored in the computer program, and the computer program is configured to execute the above-described termite nest infill planning method based on ground-penetrating radar and physical flow characteristics when it is run.

[0052] Beneficial effects: This invention solves the problems of low accuracy in reconstructing complex nests and lack of physical basis in the injection scheme, improves the three-dimensional reconstruction of underground nests and the success rate of injection and filling, and can be applied to scenarios such as liquid metal injection to obtain nest specimens and grouting and sealing. Attached Figure Description

[0053] Figure 1 This is a flowchart illustrating the steps of the termite nest irrigation planning method based on ground-penetrating radar and physical flow characteristics in the embodiments of this application.

[0054] Figure 2 This is a flowchart illustrating the steps for extracting normalized multi-attribute features in an embodiment of this application.

[0055] Figure 3 This is a flowchart illustrating the steps of performing a topological association search in an embodiment of this application.

[0056] Figure 4 This is a flowchart illustrating the steps involved in constructing a directed reachable flow graph in an embodiment of this application.

[0057] Figure 5 This is a flowchart illustrating the steps of performing time-varying gating constraint verification on candidate segments based on solidification kinetics in an embodiment of this application. Detailed Implementation

[0058] Example 1: This example provides a termite nest irrigation planning method based on ground-penetrating radar and physical flow characteristics, such as... Figure 1 As shown, the method includes the following steps:

[0059] Step 101: Obtain ground-penetrating radar detection data of the target area, extract normalized multi-attribute features, and classify anomalies based on a pre-built feature fingerprint database to obtain an anomaly spatial label dataset; wherein, the feature fingerprint database can be established based on the ground-penetrating radar response features of excavated verification samples, grouted samples, or artificially constructed simulation samples.

[0060] Specifically, ground-penetrating radar (GPR) data refers to the reflected signals from the subsurface medium obtained by non-destructive electromagnetic wave scanning of the surface of a target area. A two-stage adaptive scanning strategy can be employed to acquire this data. In the first stage, a sparse grid of GPR survey lines is used for initial scanning, and sections with reflected energy higher than the background average are marked using a sliding window energy detection method. In the second stage, denser survey lines are used for directional, intensified scanning of the marked sections to acquire high spatial resolution raw detection data.

[0061] Accordingly, the process of extracting normalized multi-attribute features is used to eliminate the absolute reflection amplitude deviation caused by differences in soil base moisture content and target burial depth under different detection scenarios, so that the extracted signal features have physical comparability across scenarios.

[0062] Furthermore, a pre-constructed feature fingerprint database stores the standard feature distributions of known categories of anomalies. This feature fingerprint database can be pre-established based on the ground-penetrating radar response features of excavated verification samples, cast-in-place specimens, or artificially constructed simulation samples; its specific construction method will be detailed in subsequent embodiments. By measuring the distance between the extracted normalized multi-attribute features and this feature fingerprint database, a classification label can be assigned to each suspected anomaly in the detection data. The specific categories of anomaly classification can include main nests, secondary nests (including fungal gardens), and tree root interference bodies. Since the cross-sectional scale of ant trails is usually much smaller than the spatial resolution capability of ground-penetrating radar at the corresponding burial depth, ant trails are not considered as an independent identification category; the connection relationship between ant trails between cavities will be inferred and constructed based on the topological rules of nests in the subsequent topological association search step. The classification confidence is determined by the reciprocal of the distance measurement. Finally, all anomaly information with three-dimensional spatial coordinates, classification labels, and classification confidence is summarized to form an anomaly spatial label dataset.

[0063] Step 102: Combining the preset biological morphological prior constraints, perform topological association search and self-consistency verification and correction in the anomaly spatial label dataset to reconstruct the three-dimensional structure connectivity graph of the nest.

[0064] In this embodiment, the preset biological morphological prior constraints refer to transforming the common knowledge of termite nest architecture into computer-readable numerical boundary conditions. These constraints include the allowable range of the equivalent diameter of the main nest cavity, the allowable range of the equivalent diameter of the secondary nest cavity, the maximum allowable distance between adjacent cavities, the maximum allowable single-segment bending angle of the connecting path between cavities, and the default reference value range for the ant tunnel cross-section based on species experience. The specific values ​​of the above parameters are configured specifically according to the dominant termite species, soil type, and historical excavation verification data of the target area, and the values ​​may differ under different engineering scenarios.

[0065] As an optional implementation, the main nest anchor point is used as the starting point when performing topological association search. The location logic of the main nest anchor point is to select the centroid coordinates of the anomaly with the highest classification confidence and burial depth falling within the typical biological range from the anomaly spatial label dataset. The preset typical burial depth range of the main nest is between 0.3 and 1.5 m. Using this coordinate as the center, a spatial search algorithm is used to perform topological association between the identified cavity nodes. For cavity nodes such as the main nest and secondary nest that can be identified by ground penetrating radar, a node set is established based on their spatial coordinates and classification labels. For the ant trail connections between cavity nodes, since their cross-sectional scale is lower than the resolution capability of ground penetrating radar, a topological inference method is used for construction: based on prior rules such as the spatial distance between adjacent cavities, the depth hierarchy relationship, and the central radiation characteristics of the main nest, under the condition of satisfying the maximum allowable spacing between adjacent cavities and the path bending angle constraints, the inferred ant trail connection edges between cavities are automatically generated, and default cross-sectional reference values ​​based on species experience are assigned to the inferred ant trails. This forms an initial topological skeleton composed of nest chamber nodes and inferred ant trail edges.

[0066] Optionally, due to the potential for misjudgment caused by the influence of complex underground media on radar signals, the initial topological framework may exhibit local structures that violate biological common sense. A self-consistency verification and correction are then performed. The deviation of the geometric parameters of each cavity node in the current framework from its allowable boundary is calculated, and unreasonable labels are iteratively corrected based on this deviation. For example, nodes with equivalent diameters exceeding the upper limit of the allowable range for secondary nests are reclassified, or presumed connection paths with significantly excessive crossing distances are removed. For secondary nests that still have isolated cavity nodes after correction (i.e., those not forming any connection path with the main network), presumed ant trail edges are automatically added based on spatial proximity to connect them to the main network, and the newly added presumed connections are assigned a low-confidence label. After iterative correction, a three-dimensional connectivity graph of the nest structure with clear three-dimensional coordinates, topological connections, and geometric dimensions is output.

[0067] Step 103: Based on the physical flow characteristics, determine the directed reachability of the paths in the three-dimensional structure connectivity graph of the nest, and construct a directed reachable flow graph.

[0068] As an alternative implementation, the three-dimensional connectivity diagram of the nest only represents the geometric connection state. However, during actual grouting operations, the flow of fluid in the underground network is subject to a combination of limitations, including gravity, capillary resistance, and the fluidity decay caused by changes in the fluid's own physical properties over time. For high-temperature liquid metal grouting, this fluidity decay mainly stems from the solidification and shrinkage effect caused by temperature loss; for grouting and sealing, it mainly stems from the increased consistency effect caused by the initial setting of the grout hydration. Therefore, it is necessary to introduce physical flow characteristics for a secondary determination at the fluid dynamics level.

[0069] Accordingly, the physical flow characteristics include static and solidification kinetic characteristics. Regarding static characteristics, the elevation difference generated by the upward-sloping ant tunnel segment is evaluated to determine whether it exceeds the balance point between capillary resistance and hydrostatic head that liquid aluminum can overcome at the current cross-sectional diameter. Regarding flow attenuation characteristics, the effective channel cross-section reduction effect caused by heat dissipation solidification or initial hydration of the injected fluid is considered, and the limiting distance the fluid can travel before losing flowability is calculated. Only paths that simultaneously satisfy the gravity climb condition and can traverse the entire length of the segment before loss of flowability are marked as effective connectivity directions. By traversing all paths in the network, the original undirected geometrically connected graph is transformed into a directed reachable flow graph with directional properties that reflects the actual flow capacity of the fluid.

[0070] Step 104: Based on the cavity volume and path connectivity of each node in the directed reachable flow graph, select the set of injection ports according to the maximum coverage principle, calculate the allocation amount and arrange the execution sequence to obtain the multi-point injection planning scheme.

[0071] In this embodiment, based on a directed reachable flow graph, the irrigation network is abstracted as a directed graph structure with candidate ground nodes as sources and underground cells as sinks. The process of selecting the set of irrigation ports according to the maximum coverage principle maximizes coverage of all underground cell nodes using as few ground boreholes as possible. Specifically, a greedy incremental strategy can be adopted, successively selecting the surface projection point that can cover the most uncovered cells as the actual irrigation port.

[0072] Furthermore, for each selected injection port, the corresponding fluid distribution volume is calculated. For a shared nest chamber accessible by multiple injection ports, a distribution model based on fluid flow resistance is adopted. The cavity volume is allocated according to the inverse ratio of the flow resistance of each converging path, and the volumes of each ant tunnel cavity filled by that injection port are added together to form the single-port distribution volume.

[0073] Accordingly, the execution sequence is arranged to avoid incomplete filling or solidification blockage of the shared channel by the pre-cooled fluid when multiple ports are injected simultaneously or in a disordered manner. Based on the physical length of the directed path and the cooling time window, the order and time interval constraints of the operation of each injection port are planned. The final output is a multi-point injection planning scheme that includes the ground coordinates of the injection ports, the fluid distribution volume, and the execution sequence with time constraints.

[0074] Example 2: This example provides a detailed explanation of the process for extracting normalized multi-attribute features, such as... Figure 2 As shown, it specifically includes:

[0075] Step 201: Identify the anomalous body region in the ground penetrating radar detection data, and delineate an exclusion buffer zone outward from the anomalous body region as a reference.

[0076] As an alternative implementation, in the profile images presented by ground-penetrating radar (GPR) data, anomalies typically appear as a combination of waveforms with energy convergence or phase abrupt changes. For example, in a two-dimensional profile corresponding to a single survey line, it may appear as a hyperbolic diffraction arc or a group of strongly reflected waves. The following uses a two-dimensional profile of a single survey line as the operating object to illustrate the specific process of normalized feature extraction. The fusion of results from multiple survey lines into a three-dimensional spatial dataset will be described later.

[0077] Specifically, the operation to determine the anomalous region involves establishing a bounding box or equivalent boundary along the outermost envelope of the diffraction arc or wave group. The space enclosed by this equivalent boundary is defined as the anomalous region, and its horizontal span is defined as the equivalent horizontal scale of the anomalous body.

[0078] Furthermore, when extracting the background signal for normalization, it is necessary to avoid the waveform tailing of the anomaly itself and near-field scattering interference. Therefore, based on the equivalent boundary of the determined anomaly region, an exclusion buffer zone is delineated by extending a certain distance horizontally outward. This extension distance can be set as a preset multiple of the equivalent horizontal scale of the anomaly; the specific multiple can be determined by those skilled in the art based on the soil homogeneity conditions of the target area. Detection data located within the exclusion buffer zone are not considered as background signal sampling objects.

[0079] Step 202: Define an annular candidate window in the outer region of the exclusion buffer zone, and extract background signal samples within the annular candidate window.

[0080] Specifically, the annular candidate window is a transitional spatial interval outside the exclusion buffer zone, used for collecting background signals. Its inner boundary is adjacent to the outer edge of the exclusion buffer zone, and its outer boundary is determined by the maximum reference distance. This maximum reference distance can also be set as a preset multiple of the equivalent horizontal scale of the anomaly. Within the area defined by this annular candidate window, radar echo signals from multiple locations along the detection line are extracted as background signal samples. To ensure statistical validity, the number of extracted background signal samples should meet a preset lower limit, not less than a preset minimum number of channels, which can be determined by those skilled in the art based on statistical significance requirements.

[0081] Step 203: Calculate the attribute variation coefficient of the background signal sample, and perform a homogeneity test on the annular candidate window in combination with the preset homogeneity threshold.

[0082] Accordingly, after obtaining the background signal samples, it is necessary to verify whether they represent the response of a homogeneous substrate medium. For each extracted background signal sample, the value of its specific physical property is calculated. Specific physical properties may include instantaneous amplitude and dominant frequency. The standard deviation and mean of all samples on the specific physical property are calculated, and the quotient of the standard deviation and mean is used to obtain the coefficient of variation of that property; the calculation formula is as follows:

[0083] CV attr =SD attr / Mean attr ;

[0084] Among them, CV attr SD is the attribute variation coefficient. attr Mean is the standard deviation of all samples on this attribute. attr This is the mean of all samples on this attribute.

[0085] Specifically, the calculated attribute variation coefficient is compared numerically with a preset homogeneity threshold. The preset homogeneity threshold reflects the maximum allowable fluctuation range of the background medium and can be determined by those skilled in the art through verification experiments in representative test areas. If the attribute variation coefficient is not greater than the preset homogeneity threshold, the current annular candidate window is deemed to have passed the homogeneity test, indicating that the selected area has a uniform medium and no obvious interference; if it is greater than the threshold, it is deemed to have failed.

[0086] Step 204: If the homogeneity test is not passed, the range of the annular candidate window is expanded outward, and sampling and testing are repeated until the homogeneity test is passed.

[0087] As an optional implementation, the number of times the outward expansion range is expanded is limited to a preset upper limit. If the homogeneity test is still not passed after reaching the maximum number of executions, the expansion operation is terminated, and a low reliability label is set for the corresponding anomaly region. This causes the normalized multi-attribute features calculated for the anomaly region to carry the low reliability label, thereby reducing its confidence weight in subsequent anomaly classification steps.

[0088] Specifically, when the signal samples within the annular candidate window exhibit strong heterogeneity, the area may contain soil interfaces or unidentified minor disturbances. In this case, it is necessary to expand the sampling range to dilute the influence of local anomalies. The specific operation involves shifting the outer boundary of the annular candidate window horizontally outward, for example, doubling the original maximum reference distance. Background signal samples are then re-extracted within the expanded window, and the coefficient of variation is calculated and compared again. This expansion and testing process is repeated until the recalculated attribute coefficient of variation meets the homogeneity threshold condition.

[0089] In some optional implementations, the expansion range of the annular candidate window is limited by a maximum number of executions. If the continuous expansion operation reaches the maximum number of executions, and the recalculated attribute variation coefficient is still greater than a preset homogeneity threshold after the preset maximum number of executions, the expansion operation is terminated. In this case, the normalized multi-attribute features extracted from the anomaly region will be assigned a low-reliability label. When subsequent anomaly classification is performed, the classification confidence score of features with low-reliability labels will be multiplied by a decay coefficient of less than 1 to reduce the interference of this poor-quality sample on topology reconstruction.

[0090] Step 205: Calculate the background reference value based on the background signal samples within the annular candidate window that have passed the homogeneity test.

[0091] In this embodiment, after the annular candidate window passes the homogeneity test, the attribute values ​​of all background signal samples contained within the window are arithmetically averaged, and this average value is used as the formal background reference value for the specific depth layer in the detection scenario. For the extracted multiple attributes (such as amplitude and dominant frequency), the corresponding formal background reference values ​​are calculated separately.

[0092] Correspondingly, if the annular candidate window fails the homogeneity test after reaching the maximum number of executions, it indicates that there is heterogeneous interference in the medium environment surrounding the anomaly, and the acquired background signal samples cannot represent the true response of the uniform substrate. In this case, the mean of the samples within the window is not used as the formal background reference value, but is only recorded as a reference coarse estimate. This coarse estimate is only used to complete the formal benchmark conversion operation in the normalization calculation below, but its output normalized multi-attribute features will carry a low reliability label according to the provisions of step 204, so as to reduce the contribution of the poor sample to the topology reconstruction by attenuating the confidence weight in the anomaly classification stage.

[0093] Step 206: Use background reference values ​​to perform benchmark conversion on the original attribute values ​​of the anomaly region to obtain normalized multi-attribute features.

[0094] In one alternative implementation, the original attribute values ​​of the anomalous body region are extracted from the ground-penetrating radar detection data, and the original attribute values ​​are benchmarked and converted using background reference values ​​to obtain normalized multi-attribute features.

[0095] Specifically, for anomaly regions, the original attribute values ​​of their own radar echo signals are extracted. Since the numerical characteristics and physical meanings of different physical attributes differ fundamentally, using a uniform division operation for normalization may distort the physical meaning of some attributes or cause numerical overflow when the background reference value approaches zero. Therefore, a matching normalization method is selected for benchmark conversion based on the physical type of each attribute. For amplitude and gradient attributes (such as instantaneous amplitude peak value and signal attenuation gradient), ratio normalization is used, that is, dividing the attribute value of the anomaly region by the corresponding background reference value; when the absolute value of the background reference value is lower than a preset lower limit threshold for numerical stability, this lower limit threshold is used instead for division to prevent numerical overflow. For frequency shift attributes (such as dominant frequency offset), difference normalization is used, that is, calculating the difference between the measured dominant frequency value of the anomaly region and the dominant frequency value of the background reference. For polarity-related attributes (such as waveform polarity patterns), since they contain directional information and their values ​​can be positive or negative, numerical division or difference operations are not applicable. Instead, an independent symbolic encoding method is used for calculation, the specific encoding method of which will be explained later in this step. For scale-related attributes (such as the width of a hyperbolic diffraction aperture), relative deviation normalization is used, i.e., the difference between the attribute value of the anomaly region and the background reference value is calculated and then divided by the background reference value. After being converted according to attribute type, the normalized results of each dimension together constitute a normalized multi-attribute feature that eliminates the influence of the environmental background base. This feature can truly reflect the physical property differences of the anomaly itself relative to the surrounding medium from multiple physical dimensions.

[0096] Furthermore, the extraction of normalized multi-attribute features can also include a feature encoding process. For example, encoding waveform polarity patterns. Since the polarity of a ground-penetrating radar wavelet reverses when it encounters a dielectric constant transition interface, this characteristic is scalarized and encoded; the calculation formula is as follows:

[0097] x polar =s*r;

[0098] Where, x polar The polarity mode scalar is defined as follows: 's' represents the polarity sign of the first half-cycle, and 'r' represents the absolute ratio of the amplitude of the second half-cycle to that of the first half-cycle. The polarity sign of 's' in the first half-cycle is positive when the dielectric constant increases and negative when it decreases. This scalar encoding can simultaneously characterize the dielectric gradient direction and asymmetry of the interface, serving as a dimension for normalized multi-attribute features in subsequent classification. The hyperbolic diffraction aperture width is defined as the horizontal distance measured in the profile with the attenuation threshold of the peak amplitude (e.g., -6 dB) as the cutoff boundary.

[0099] In practical 3D nest reconstruction applications, multiple parallel or intersecting ground-penetrating radar (GPR) survey lines are laid out in a grid within the target area, and each survey line generates an independent 2D profile. The above processing procedure is performed on each 2D profile of all survey lines to obtain the anomalous body regions and their corresponding normalized multi-attribute features on each profile.

[0100] Furthermore, after independently processing all profiles, the two-dimensional results need to be fused into three-dimensional spatial data. Specifically, based on the known ground spatial coordinates of each survey line and the depth-horizontal distance coordinate system within the profile, the identified anomalous body regions on each profile are uniformly projected into the three-dimensional spatial coordinate system of the target area. For anomalous body regions that are spatially adjacent and have similar normalized multi-attribute features on different survey line profiles, they are determined to be the cross-sectional manifestations of the same subsurface anomalous body on different profiles, and cross-profile association is performed. The cross-sectional contours of the same anomalous body on each profile after association are fused using three-dimensional interpolation to estimate its three-dimensional spatial envelope and centroid coordinates. Finally, the fused three-dimensional spatial coordinates, classification labels, classification confidence scores, and normalized multi-attribute features of each anomalous body are summarized to form an anomalous body spatial label dataset with three-dimensional spatial coordinates, which is used for topological association search and self-consistency verification and correction.

[0101] Example 3: This example further details the steps involved in performing the topological association search, such as... Figure 3 As shown, it specifically includes:

[0102] Step 301: Determine the starting node of the main nest from the anomaly spatial label dataset.

[0103] Specifically, when performing topological reconstruction in 3D space, a highly reliable reference point is established. All label points initially classified as primary nest candidates are obtained from the anomaly spatial label dataset. The classification confidence score and corresponding 3D coordinates of each primary nest candidate are extracted. The depth component of these 3D coordinates is extracted, and it is determined whether this depth component lies within a preset typical burial depth range for primary nests. This preset typical burial depth range can be set to 0.3~1.5m. The anomaly with the highest classification confidence score and meeting the depth component requirements is selected, and its 3D coordinates are determined as the starting node for constructing the entire connected graph's primary nest.

[0104] Step 302: Based on the starting node of the main nest, and combined with the path angle deflection constraint and cumulative length constraint specified in the biological morphology prior constraint, search and associate nodes in the anomaly spatial label dataset to construct an initial topological skeleton, and input the initial topological skeleton into the self-consistency verification and correction step to reconstruct the three-dimensional structure connectivity graph of the nest.

[0105] Accordingly, after determining the starting node of the main nest, the discretely distributed anomaly tags in the surrounding space need to be organized into a connected graph according to the actual physical connectivity. This connected graph presents a radial topology with the main nest as the core, but it is not limited to a tree structure. Lateral connecting edges are allowed between chambers at the same level to accommodate the mixed topological forms that may occur in real nests, which are between tree and network structures.

[0106] In one optional basic implementation, the horizontal plane can be divided into multiple fixed search sectors at fixed angular intervals, and nodes meeting the conditions can be searched radially outward within each sector. However, real underground material channels often exhibit high curvature bends, and the fixed sector search method is prone to search path breaks when the channel crosses the sector boundary, leading to missed detections of connected structures. To overcome the limitations of fixed sectors, this embodiment preferably adopts an adaptive multi-path tracking search mechanism, transforming the global angular division constraint into a local segment-by-segment angular deflection verification, allowing the search path to naturally adapt to the direction of the curved channel and complete the construction of the initial topology skeleton.

[0107] Optionally, nodes are searched and associated, specifically using a multi-path tracing search mechanism, including:

[0108] Step 303: Construct a path seed with the main nest starting node and assign it an initial direction of movement, and place it into a priority queue sorted by confidence.

[0109] Optionally, the adaptive multi-path tracing search uses path seeds as the smallest computational unit for expansion. With the main nest starting node as the core, all anomaly tag points belonging to the corresponding category are retrieved within a preset initial search radius. The main nest starting node is paired with each retrieved tag point to generate multiple independent path seeds. An initial forward direction is calculated for each path seed, which is represented as a three-dimensional spatial vector pointing from the main nest starting node to the corresponding tag point. A priority queue data structure is established, and all generated path seeds are placed into this priority queue. The priority queue is sorted in descending order based on the classification confidence of the current end node of the path seed, so that path seeds with higher confidence are processed first.

[0110] Step 304: Extract active paths from the priority queue and retrieve candidate nodes in the neighborhood of their end nodes.

[0111] Specifically, in each round of expansion iteration, the highest priority path seed is extracted from the head of the priority queue and marked as an active path. The 3D coordinates of the terminal node of the active path are extracted, and a spherical neighborhood is constructed with these 3D coordinates as the center and a preset maximum allowable single segment length as the radius. In the anomaly spatial label dataset, the anomaly label points located within this spherical neighborhood that have not yet been fixed are extracted as candidate nodes for the active path to continue extending forward.

[0112] Step 305: Perform local direction deflection verification and cumulative length verification on candidate nodes. Add candidate nodes that meet the path angle deflection limit and cumulative length limit as new branches to the current path, update the local forward direction, and put them back into the priority queue. If all candidate nodes do not meet the verification, mark the current active path as a broken path. Repeat the above process of extracting active paths from the priority queue, retrieving candidate nodes, and adding branches until the priority queue is empty.

[0113] Further, for each retrieved candidate node, the local forward direction vector of the current active path and the target direction vector from the end node of the active path to the candidate node are obtained. The spatial angle between the local forward direction vector and the target direction vector is calculated. This spatial angle is numerically compared with the path angle deflection limit in the biological morphology prior constraints. If the spatial angle is not greater than the limit value, the local direction deflection check is passed. The cumulative total length of the entire path after including the candidate node is calculated, and it is verified whether the cumulative total length is not greater than the cumulative length limit. Only candidate nodes that pass both of the above checks are confirmed as legitimate extension points. An independent copy branch of the current active path is created for each legitimate extension point, and the legitimate extension point is appended to the end of the sequence of the copy branch. The target direction vector is recalculated using the three-dimensional coordinates of the newly added node and the preceding node, and the local forward direction of the branch is updated accordingly. After the update is completed, all newly generated branches are placed back into the priority queue to participate in subsequent iterations.

[0114] Step 306: Extract complete paths based on preset path termination conditions and node affiliation conflict resolution rules to construct the initial topology skeleton.

[0115] Specifically, as active paths expand, conditions need to be set to determine the path's completion status. The preset path termination conditions include three states. The first state is when a node classified as a secondary nest or other non-primary nest cell is detected in the terminal neighborhood of the active path. In this case, the node is included at the end of the path, and the path is marked as a successfully extracted complete pathway. The second state is when the cumulative length of the path has reached the length limit but has not yet encountered the terminal chamber node. In this case, it is determined to be a blind-end path and temporarily retained. After all active paths have been processed, during the self-consistency verification and correction phase, the blind-end path is used as a structure to be verified in the calculation of the total penalty value for constraint violations. If the geometric parameters of its isolated terminal node meet the allowable range of the nest cell category constraint, it is updated to a secondary nest node; otherwise, the blind-end path is removed from the initial topology skeleton. The third state is when there are no candidate nodes that meet the verification in the terminal neighborhood. In this case, the path is determined to be broken, and the branch is terminated.

[0116] Accordingly, when the priority queue is emptied, all complete pathways and blind paths are collected, and their topological connectivity states are merged to form a radial initial topological skeleton containing node coordinates and connecting edges. This radial skeleton reflects the radial expansion result with the main nest as the core, and the possible lateral connectivity relationships between chambers at the same level will be supplemented in the subsequent lateral connectivity discovery steps.

[0117] Optionally, the complete path is extracted based on preset path termination conditions and node affiliation conflict resolution rules, specifically including:

[0118] Optionally, when the same candidate node in the anomaly spatial label dataset is simultaneously included in the competition by multiple active paths, the path fit of each active path to the candidate node is evaluated separately.

[0119] Furthermore, since the search process is multi-branched and parallel, in densely populated spatial regions, the same candidate node often appears in the search neighborhood of multiple active paths, all satisfying the aforementioned deflection and length checks. This leads to node ownership conflicts. A node ownership conflict resolution rule is triggered, and for each active path competing for the candidate node, its technical suitability index for receiving the node, i.e., path suitability, is calculated independently.

[0120] Optionally, the path fit is determined based on the classification confidence of the candidate node and the local orientation deflection angle caused by including the candidate node.

[0121] In one alternative implementation, for each active path, the path fitness value corresponding to the active path is calculated by dividing the classification confidence of the candidate node by the local directional deflection angle caused by including the candidate node.

[0122] Specifically, path fit reflects the overall structural morphology's combined benefits in terms of classification certainty and geometric smoothness after including a node in the path. The classification confidence parameter of the candidate node itself is obtained, along with the specific value of the local directional deflection angle resulting from including the candidate node in the current active path. The following two parameters are used for calculation:

[0123] F path =C node / (Θ deflect +δ);

[0124] Among them, F path For path adaptability, C node The classification confidence of candidate nodes, Θ deflect To include the local directional deflection angle caused by the candidate node, δ is a small positive number to prevent the denominator from being zero, for example, δ = 0.001 radians.

[0125] Optionally, the final path assignment of the candidate node is determined based on the highest path fit, and the remaining active paths that were not selected are truncated at the corresponding positions.

[0126] Furthermore, for the candidate nodes currently experiencing a conflict, the path fitness values ​​calculated from each competing path are compared. The active path that produces the highest path fitness value is determined as the unique home path for that candidate node, allowing the node to be added to the end of this optimal path for subsequent expansion.

[0127] Correspondingly, for the remaining active paths that participated in the competition but failed to achieve the highest fitness value, their extension in the direction of this candidate node is deemed invalid. Subsequent expansion of these unselected paths is stopped at the corresponding location to ensure that the backbone topology generated during the radial search phase does not exhibit redundant overlap due to search competition. It should be noted that the single-attribute rule here only applies to node competition resolution during the radial search process. To ensure the stability of backbone path extraction, lateral connectivity relationships that may exist between chambers at the same level are not covered by the radial search and will be handled separately in subsequent steps through a dedicated lateral connectivity discovery mechanism.

[0128] According to one aspect of this application, it also includes: performing lateral connectivity discovery based on the radial initial topological skeleton to supplement any non-radial connectivity relationships that may exist between chambers at the same level.

[0129] The true topology of termite nests is not a strictly tree-like structure; lateral tunnels sometimes connect secondary nests or fungal gardens at the same depth level. The radial search performed in the above steps, centered on the main nest and extending in a radial direction, ensures the complete extraction of the main pathways, but it cannot cover the lateral paths between nodes at the same level. Therefore, a dedicated lateral connection discovery is needed after the radial search is completed.

[0130] Specifically, the process iterates through all confirmed chamber nodes in the initial radial topology framework, extracting the 3D coordinates and depth level of each node. For any two chamber nodes whose depth difference is no greater than a preset depth tolerance at the same level, their 3D spatial distance is calculated. If this spatial distance is no greater than the maximum allowable spacing between adjacent chambers specified in the biological morphology prior constraints, and the angle between the direction of the line connecting the two nodes and the horizontal plane is no greater than a preset lateral connection inclination angle threshold (e.g., 30°), then it is determined that the node pair has the possibility of lateral connectivity. A lateral connection edge is established between the node pair, and this lateral connection edge is assigned the attribute of a presumed ant trail and an initial confidence level lower than that of the radial trunk edge. The addition of the lateral connection edge expands the initial topology framework from a purely radial radial structure to a hybrid topology network that allows for the existence of local loops, which is closer to the biological reality of termite nests. The newly added lateral connection edge will be included in the subsequent self-consistency verification and correction process along with the radial connection edge. If its geometric parameters violate the biological morphology prior constraints, it will be automatically removed during the verification and correction stage.

[0131] Example 4: This example further details the self-consistency verification and correction process, specifically including:

[0132] As an optional implementation, the geometric parameters of each node obtained through topological association search are extracted, and the dimensionless deviation of each geometric parameter from the corresponding allowable range in the prior constraints of biological morphology is calculated.

[0133] Specifically, the extracted geometric parameters are verified using the equivalent diameter of the nest chamber nodes obtainable by ground-penetrating radar. For ant trail connecting edges generated based on topological patterns, since their cross-sectional parameters are default values ​​assigned by species experience rather than measured values, their equivalent cross-sectional diameters are not included in the penalty calculation; only the path length of the inferred ant trail edge is used as an auxiliary verification parameter in the constraint evaluation.

[0134] Optionally, the prior constraints on biological morphology include allowable lower and upper bounds for each of the aforementioned geometric parameters. Since different geometric parameters have different dimensions and magnitudes, calculating the absolute difference between the estimated value and the allowable boundary can cause the error of larger numerical parameters to mask the error of smaller numerical parameters. Therefore, the deviations of each parameter need to be dimensionless. The estimated value for each geometric parameter is obtained, and the absolute violation amount by which this estimated value exceeds the allowable lower or upper bound is calculated. This absolute violation amount is then divided by the allowable range width of the corresponding geometric parameter. The allowable range width is the difference between the allowable upper and lower bounds.

[0135] As an optional implementation, a nonlinear penalty is applied to the degree of dimensionless deviation, and the total penalty value for constraint violation is obtained by summing the results.

[0136] Furthermore, to improve the sensitivity of the self-consistency verification model to anomalous structures, a nonlinear function is used to map the aforementioned dimensionless deviation. In a preferred embodiment, a nonlinear penalty is applied to the dimensionless deviation, including: for the deviation portion exceeding the corresponding allowable range boundary, a squared calculation is used to enhance the penalty effect, so that the geometric parameter that deviates further from the allowable range occupies a higher sensitivity weight in the total penalty value for constraint violation. The deviations of all parameters to be verified after the squared penalty are summed; the calculation formula is as follows:

[0137] P total =Σ i ((max(0,a iL -v i ) / (a iU -a iL )+max(0,v i -a iU ) / (a iU -a iL )) 2 );

[0138] Among them, P total To constrain the total penalty value for violations, Σ i This represents the summation of all geometric parameters to be verified, v i Let a be the estimated value of the i-th geometric parameter. iL Let a be the allowable lower bound of the i-th geometric parameter in the biological morphology prior constraints. iU Let be the upper bound allowed for the i-th geometric parameter in the prior constraints of biological morphology, and max() be the maximum value function. When the estimated value falls between the lower and upper bounds allowed, the output of the maximum value function is 0, and this parameter does not contribute to the total penalty value; when the estimated value exceeds the bounds, the cross-bounds distance is dimensionless and squared, and included in the sum as a positive penalty term.

[0139] As an optional implementation, iterative correction is performed based on the total penalty value for constraint violation until the total penalty value for constraint violation is not greater than a preset tolerance threshold. If the correction operation causes a change in the node category in the anomaly spatial label dataset, then, with the changed node as the center, a local topological association re-search is triggered within a local spatial range centered on the changed node and with the maximum allowed total length of a single ant trail as the radius, and the connection graph of the 3D structure of the nest is updated. If the correction operation does not cause a change in the node category, then the current topological connection relationship is retained, the total penalty value for constraint violation is recalculated based on the numerically adjusted geometric parameters, and the next iteration begins.

[0140] Accordingly, a tolerance threshold is pre-defined to characterize fault tolerance. The calculated total penalty value for constraint violation is compared with this tolerance threshold. If the total penalty value for constraint violation is greater than the tolerance threshold, an iterative correction process is initiated. The geometric parameter that contributes the most to the total penalty value in the current topology skeleton is located, and a category correction is performed based on the violation type. For example, if the estimated equivalent diameter of the cross section of an antway edge is greater than its corresponding upper allowable bound, the anomalous body feature at that location is determined to be more consistent with the secondary nest morphology, and the node category of the anomalous body in the spatial label dataset is changed from antway to secondary nest.

[0141] As an optional implementation, to ensure the algorithm terminates within a finite time, the iterative correction process has a preset maximum iteration limit. If the total penalty for constraint violations still exceeds the tolerance threshold after reaching the maximum iteration limit, the iteration is forcibly terminated. The topological skeleton obtained in the current iteration is taken as the final output, and nodes in the structure that still have constraint violations are assigned a low reliability label to reduce their reference weight in subsequent injection planning steps. The maximum iteration limit can be adaptively set by those skilled in the art based on the number of anomalies and computational resource conditions in the actual operating scenario.

[0142] Optionally, since changes in node categories in the underlying dataset disrupt the connection logic of the original topological paths, a backtracking trigger mechanism needs to be introduced to maintain structural consistency. The 3D coordinates of the node whose category has changed are obtained, and the affected local spatial range is extracted using these coordinates as the center. The radius of this local spatial range can be set to the maximum allowed total length of a single ant trail. The original topological connection state within this local spatial range is cleared, and the aforementioned multi-path tracing search is re-initiated from the nearest adjacent known nest node to construct a new local topological relationship. After the re-search is completed, the new local connections are spliced ​​back into the global network, and the total penalty value for constraint violations is recalculated until the total penalty value is no greater than a preset tolerance threshold. Finally, a closed-loop 3D nest structure connectivity graph in terms of geometric attributes and category logic is output.

[0143] In some optional implementations, if the iterative correction operation in the current round only performs numerical adjustments to the geometric parameters, such as forcibly reducing the equivalent volume of a node to the allowable upper bound, without causing a change in the node category in the anomaly spatial label dataset, then the underlying search label data is determined to be unaffected. In this case, the backtracking operation that triggers a re-search of local topological associations is skipped, the existing topological connections are preserved, and the calculation of the total penalty value for constraint violation in the next round is carried out based on the numerically adjusted geometric parameters.

[0144] Example 5: This example further details the process of determining the directed reachability of paths in the three-dimensional connectivity graph of the nest based on physical flow characteristics and constructing a directed reachable flow graph.

[0145] In this embodiment, the preparation of nest specimens by high-temperature liquid aluminum injection is used as a typical application scenario for detailed explanation. Due to its high-temperature melting characteristics, liquid aluminum faces significant physical constraints such as solidification and diameter contraction when flowing in underground micro-channels, making it an injection scenario with high requirements for modeling the time-varying decay of fluidity. When the injection fluid is a room-temperature cement-based slurry used for sealing purposes, the solidification kinetics process in the physical flow model below is replaced by the initial setting process of slurry hydration, and the corresponding thermophysical parameters (such as melting point, latent heat, and wall thermal conductivity) are replaced by rheological parameters such as the initial setting time and consistency time-varying function of the slurry. The solidification rate coefficient is replaced by the consistency growth rate coefficient, and the overall judgment logic and step framework remain consistent.

[0146] In one optional implementation, the directed reachability of paths in the three-dimensional connectivity graph of the nest is determined based on physical flow characteristics. The cavity volume and path geometry of each node in the three-dimensional connectivity graph of the nest are preserved to construct a directed reachable flow graph, such as... Figure 4 As shown, it includes the following steps:

[0147] Step 501: Extract the path segments from the three-dimensional connected graph of the nest structure.

[0148] Specifically, the three-dimensional structural connectivity graph of the nest consists of nodes representing nest chambers and edges representing ant trails. Extracting each path segment involves traversing and extracting each edge in the connectivity graph and reading the physical and geometric properties assigned to that edge. These physical and geometric properties include the equivalent diameter of the cross-section, the actual segment length, the three-dimensional spatial orientation vector, and the corresponding wall material type. The extracted path segments constitute the basic computational units for fluid dynamics analysis.

[0149] Step 502: Based on preset fluid physical property parameters (including fluid density, surface tension and contact angle) and static characteristics, evaluate the gravitational potential difference and capillary resistance of each path segment, and determine the candidate segments that meet the static reachability conditions.

[0150] As an alternative implementation, when fluid flows through underground micro-channels, it is driven not only by gravity but also by capillary forces caused by interfacial tension. In this embodiment, high-temperature liquid aluminum is used as the injection fluid. Since liquid aluminum and wall materials such as soil typically form a non-wetting system, capillary forces manifest as resistance to fluid movement. The static characteristics are evaluated, and the elevation difference between the two ends of the path segment is calculated based on the three-dimensional spatial direction vector. For downward-sloping path segments, gravity does positive work, and the path is determined to be statically reachable.

[0151] Accordingly, for upward-sloping path segments, it is necessary to assess whether the hydrostatic head at the injection port can overcome the capillary resistance within that segment. Specifically, based on the fluid's surface tension, the contact angle between the fluid and the wall material, and the equivalent diameter of the path segment's cross-section, the equivalent head of capillary resistance is calculated using the Young-Laplace equation. This head is then subtracted from the hydrostatic head at the injection port to obtain the actual maximum rise height of the fluid within that segment. The smaller the equivalent diameter, the greater the capillary resistance and the lower the maximum rise height.

[0152] In some scenarios, for upward-sloping path segments, it is necessary to calculate the actual maximum rise height after deducting the equivalent head due to capillary resistance from the hydrostatic head at the injection port; the calculation can be simplified as follows:

[0153] h max =H0-4*γ Al *abs(cos(θ c )) / (ρ Al *g*d);

[0154] Among them, h max The maximum climbing height is given by H0, where H0 is the hydrostatic head of the fluid column at the injection port, and γ is the maximum climbing height. Al Let θ be the surface tension of the fluid, abs() be the absolute value function, and θ be the absolute value function. c ρ is the contact angle between the fluid and the wall material. Al Let be the fluid density, g be the gravitational acceleration, and d be the equivalent diameter of the cross-section of the path segment. The hydrostatic head H0 of the fluid column at the injection port is equal to the actual height difference between the liquid surface of the liquid storage container of the injection equipment and the injection port. It can be calculated based on the equipment layout plan after on-site survey, or the corresponding equivalent liquid column height can be obtained by pre-setting the target pressure.

[0155] Furthermore, the calculated maximum climb height is compared with the actual upward elevation difference of the path segment. If the actual upward elevation difference is not greater than the maximum climb height, the fluid is deemed capable of overcoming gravity and capillary resistance. Based on the above evaluation, path segments with the possibility of hydrostatic flow are retained as candidate segments that meet the hydrostatic reachability condition.

[0156] Step 503: Perform time-varying gating constraint verification on the candidate segment based on the wall material thermophysical parameters and solidification kinetics characteristics of the candidate segment to verify whether the fluid can pass through the entire length of the segment before it completely solidifies and blocks the channel.

[0157] Specifically, static property assessment alone cannot reflect the heat loss process of a fluid in the surrounding low-temperature medium. Once a high-temperature fluid is injected into an underground network at ambient temperature, intense heat exchange occurs at the pipe walls, causing the fluid to gradually solidify from the outside in. The effective flow cross-section continuously shrinks over time until the channel is closed. Therefore, solidification kinetics must be introduced to perform gating constraint verification that includes a time dimension, in order to intercept statically feasible paths that are kinetically inaccessible due to premature solidification.

[0158] Optionally, time-varying gating constraint verification is performed on the candidate segments based on solidification kinetics, such as... Figure 5 As shown, it includes:

[0159] Step 503a: Based on the thermophysical parameters of the wall material corresponding to the candidate segment, obtain the solidification rate coefficient that characterizes the solidification and diffusion of the fluid from the pipe wall to the center.

[0160] As an optional implementation, a planar solidification approximation model is used to calculate the propagation velocity of the solidification front to compute the solidification process. This velocity is primarily determined by the thermal conductivity of the wall material. The wall material type obtained during the classification and identification phase of the candidate segment is read. If the wall material is undisturbed soil, the equivalent thermal conductivity of undisturbed soil is used; if the wall material is a compacted nest wall, due to its organic binder content, the corresponding nest wall thermal conductivity needs to be used.

[0161] Furthermore, by obtaining the fluid's melting point, latent heat of fusion, and other thermophysical parameters, and based on the approximate solution to the classical Stefan solidification problem, and by combining the thermal conductivity of the wall material, the temperature difference between the initial wall temperature and the fluid's melting point, as well as the fluid's density and latent heat of fusion, the solidification rate coefficient, which characterizes the rate at which the solidification front moves from the pipe wall to the center, is calculated. The higher the thermal conductivity of the wall material or the greater the temperature difference, the larger the solidification rate coefficient, resulting in faster solidification and a shorter usable flow time window for the fluid.

[0162] In some scenarios, to obtain parameters such as the melting point and latent heat of a fluid and calculate the solidification rate coefficient, the following formula can be used:

[0163] β=sqrt(2*k wall *(T m -T wall ) / (ρ Al *L f ));

[0164] Where β is the solidification rate coefficient, k wall T represents the thermal conductivity of the corresponding wall material. m T is the melting point of the fluid. wall ρ is the initial temperature of the wall material. Al For fluid density, L fLet be the latent heat of fusion of the fluid, and sqrt() be the square root function.

[0165] Step 503b: Calculate the flow duration corresponding to the candidate segment based on the solidification rate coefficient, the equivalent diameter of the cross section of the candidate segment, and the minimum channel threshold required for fluid flow.

[0166] Correspondingly, the solidified shell grows inward from the pipe wall, causing the effective inner diameter for fluid flow to continuously decrease. When this effective inner diameter decays to the minimum channel threshold required for fluid flow, the capillary inlet pressure tends to infinity, and fluid flow ceases. The preset minimum channel threshold d min It can be determined by calculation based on the minimum capillary inlet pressure of the fluid.

[0167] Based on the solidification rate coefficient calculated above, and combined with the difference between the initial cross-sectional equivalent diameter of the candidate segment and the minimum channel threshold, the time it takes for the solidified shell to grow from the pipe wall inward until the fluid flow cross-section shrinks to the minimum channel threshold is estimated, i.e., the flow duration. The larger the initial cross-sectional equivalent diameter or the smaller the solidification rate coefficient, the longer the flow duration, and the more abundant the flow window for the fluid.

[0168] In some scenarios, based on the solidification rate coefficient calculated above, the time it takes for the equivalent diameter of the initial cross-section to shrink to the minimum channel threshold, i.e., the flow duration, can be calculated using the following formula:

[0169] T flow =((dd min ) / (2*β)) 2 ;

[0170] Among them, T flow Let d be the flow duration, and d be the initial cross-sectional equivalent diameter of the candidate segment. min β is the minimum channel threshold, and β is the solidification rate coefficient.

[0171] Step 503c: Combine the inclination angle of the candidate segments with the time-varying reduction characteristics of the fluid velocity to estimate the thermally permeable distance that the fluid can reach during the flow duration.

[0172] Specifically, during the flow duration, the instantaneous velocity of the fluid is not constant but decreases nonlinearly as the effective pipe diameter decreases. Treating the fluid motion within the path segments as laminar flow, a time-varying velocity function is constructed based on the dynamic viscosity of the fluid, the angular component of the path segments, and the time-varying relationship of the effective pipe diameter due to solidification. Integrating this time-varying velocity function over the entire flow duration, the maximum distance the fluid can travel before the channel becomes blocked is derived, i.e., the thermally permeable distance. This distance is influenced by a combination of factors: the path angle, the initial pipe diameter, the fluid viscosity, and the solidification rate coefficient. A larger angle results in a stronger gravity-driven component, and a larger initial pipe diameter results in a longer effective flow time, both of which contribute to increasing the thermally permeable distance. Conversely, higher fluid viscosity or a faster solidification rate coefficient results in a shorter thermally permeable distance.

[0173] Step 503d: Compare the actual segment length of the candidate segment with the thermally passable distance. If the actual segment length does not exceed the thermally passable distance, the candidate segment is determined to have passed the time-varying gating constraint verification.

[0174] Optionally, the actual segment length parameter of the candidate segment is obtained and compared with the calculated thermally permeable distance. If the actual segment length is not greater than the thermally permeable distance, it indicates that the flow front has reached the far end of the path segment before the fluid solidifies inside the segment. In this case, the candidate segment is determined to be kinetically unobstructed and formally passes the time-varying gating constraint verification. Conversely, if the actual segment length is greater than the thermally permeable distance, the path will be marked as blocked and eliminated.

[0175] Step 504: The path segments that simultaneously satisfy the static reachability condition and the time-varying gating constraint verification are taken as reachable directions, and a directed reachability flow graph is constructed.

[0176] Optionally, after the aforementioned dual filtering mechanism, all path segments that can overcome static mechanical resistance and complete the entire journey within the solidification time window are selected and retained. The 3D connecting nodes and corresponding unidirectional fluid flow direction vectors of these valid path segments are extracted. The nodes in the connected graph are used as vertices, and the retained unidirectional flow direction vectors are used as directed edges, updating and outputting a directed reachable flow graph. Since each edge in this graph has been verified by calculations of fluid dynamics and heat transfer parameters, its reachability is highly consistent with the actual physical flow, providing a data foundation for subsequent injection volume allocation.

[0177] Example 6: This example further details the steps of calculating the allocation quantity and arranging the execution sequence, specifically including:

[0178] Step 601: For exclusive nodes that can be reached independently from a single injection port in the directed reachable flow graph, the cavity volume recorded in the three-dimensional connectivity graph of the nest is entirely allocated to the total amount of the corresponding injection port.

[0179] Accordingly, the directed reachability graph determines the topological reachability between the underground cell network and the candidate surface injection points. During the actual allocation of injection fluid, some edge cell nodes have only one physical path to the surface; these nodes are defined as exclusive nodes. For such exclusive nodes, the cavity volume estimated during the topology reconstruction phase is obtained, and this cavity volume is treated as an indivisible whole and fully added to the initial allocation total for the corresponding unique injection point.

[0180] Step 602: For a shared node that is simultaneously reached by multiple injection ports through different paths, each path is regarded as a series pipe and the flow resistance corresponding to each path is calculated separately.

[0181] As an alternative implementation, underground main nests or large fungal nurseries are often located at transportation hubs, with multiple ant tunnels connecting to different surface infusion points; such nodes are defined as shared nodes. In a basic implementation, the cavity volume of the node can be roughly allocated according to the proportion of the terminal cross-sectional area of ​​each converging path entering the shared node. However, the cross-sectional area allocation method ignores the changes in flow resistance caused by differences in path length, which can easily lead to uneven distribution. This embodiment preferably adopts a parallel flow distribution model based on flow resistance weighting.

[0182] Specifically, the different physical paths from each injection port to the same shared node are abstracted into parallel macroscopic branches. For each macroscopic branch, the multiple ant trail segments connected in series within it are regarded as serial pipes.

[0183] Accordingly, under the laminar flow assumption, and based on Poiseuille's law, the flow resistance of a single ant-path segment is calculated using the fluid's dynamic viscosity, the actual length of each ant-path segment, and the equivalent diameter of the cross-section. Since the flow resistance is inversely proportional to the fourth power of the pipe diameter and directly proportional to the segment length, segments with smaller cross-sections or larger lengths will contribute higher flow resistance. The total flow resistance of the macro-branch is obtained by linearly superimposing the flow resistances of all ant-path segments on the macro-branch.

[0184] In some scenarios, calculation examples are as follows:

[0185] Z path =Σ e ((128*μ Al *L e ) / (π*d e 4 ));

[0186] Among them, Zpath Σ represents the total flow resistance corresponding to a single convergence path. e This represents the summation of all ant trail segments e along the path, μ Al L is the dynamic viscosity of the fluid. e Let d be the actual segment length of segment e in the antway. e Let be the equivalent diameter of the cross section of segment e of the ant tunnel.

[0187] Step 603: Determine the volume allocation coefficient based on the inverse ratio of the flow resistance corresponding to each path, and allocate the cavity volume of the shared nodes recorded in the three-dimensional structure connectivity diagram of the nest to the corresponding injection port according to the volume allocation coefficient, so as to summarize the allocation amount of each injection port.

[0188] Specifically, fluid always preferentially flows to channels with lower resistance; therefore, the flow rate of a parallel branch is inversely proportional to its flow resistance. Obtain the total flow resistance for all paths leading to the shared node, and calculate the reciprocal of each total flow resistance, i.e., the conductance. Divide the conductance of a single path by the sum of the conductances of all converging paths to calculate the volume distribution factor specific to that path; the calculation formula is as follows:

[0189] λ i =(1 / Z path_i ) / Σ j (1 / Z path_j );

[0190] Where, λ i Z is the volume distribution coefficient for the i-th convergence path. path_i Let Σ be the total flow resistance of the i-th convergence path. j This represents the summation of the path flows to all paths leading to the shared node.

[0191] Further, the overall cavity volume of the shared node is multiplied by the calculated volume allocation coefficients to obtain several sub-volumes after the node is divided. Each sub-volume is then assigned to its corresponding injection port allocation sequence. All exclusive and shared nodes in the network are traversed, and the node sub-volumes borne by each injection port and the pipe volumes of the ant trails they pass through are globally accumulated to finally obtain the precise total fluid allocation for each injection port. Optionally, the execution timing is arranged, specifically including:

[0192] Step 604: Identify adjacent injection ports that share a path segment based on the directed reachable flow graph.

[0193] Optionally, when multiple ports are injected simultaneously, filling failure will occur if the fluid from different ports converges or is blocked within the same physical channel. Therefore, the execution timing needs to be arranged. A directed reachable flow graph is obtained, and the subgraph network that each injection port can cover is extracted. By comparing the subgraph networks of different injection ports, subgraph pairs with overlapping edges are extracted. Two injection ports with overlapping edges, i.e., sharing a path segment, are marked as adjacent injection ports with a temporal dependency. For sets of independent injection ports that do not share a path segment, it is determined that there is no temporal interference, and parallel execution can be arranged.

[0194] Step 605: Based on the equivalent diameter of the cross section of the shared path segment and the thermophysical properties of its wall material, obtain the longest permissible interval time for the fluid injected into the preceding injection port to solidify into blockage in the shared path segment.

[0195] Optionally, for adjacent injection ports with a time-dependent relationship, the order of execution (preceding injection port and subsequent injection port) is determined based on the topological relationship. Fluid from the preceding injection port flows first through the shared path segment and begins to solidify at the pipe wall. The subsequent injection port must complete fluid injection before the shared path segment is blocked by the solidified shell of the preceding port. The flow duration formula calculated based on the solidification rate coefficient and minimum channel threshold is applied, and the equivalent diameter of the shared path segment is substituted to obtain the maximum allowable interval time. If the start time of the subsequent injection port is later than the sum of the start time of the preceding injection port and this maximum allowable interval time, the shared path segment will lose its flow capacity.

[0196] Step 606: Based on the cumulative path length from the preceding injection port to the far end of the shared path segment and the average fluid velocity in the directed reachable flow graph, obtain the shortest allowable interval time required for the fluid injected from the preceding injection port to flow to the far end of the shared path segment.

[0197] Optionally, in addition to preventing solidification and blockage, it is also necessary to prevent the two fluid streams from colliding within the shared path segment, as this collision can trigger an airlock effect, leading to voids. Therefore, the initiation of a subsequent injection port must wait for the fluid from the preceding injection port to pass through the shared path segment. The cumulative spatial length from the preceding injection port to the far end of the shared path segment is obtained, and combined with the average fluid velocity, the time required for the fluid to traverse this distance is calculated; this time is defined as the minimum permissible interval.

[0198] Step 607: Combining the longest and shortest allowable interval times, arrange the execution timing with time interval window constraints for adjacent injection ports.

[0199] Optionally, the longest allowed interval provides the upper limit boundary for the timing arrangement, and the shortest allowed interval provides the lower limit boundary; together, they constitute an effective time interval window constraint. This constraint is attached to the control logic of the corresponding adjacent injection ports. When generating a multi-point injection planning scheme, not only is the start-up order of each injection port output, but the delay waiting time interval between adjacent actions is also clearly defined. For example, the delay waiting time interval between subsequent actions is set to T after the execution of the preceding action. min To T max Initiating within the designated range avoids the risks of counter-current gas lock and solidification blockage at the physical level.

[0200] For example, the execution process of the present invention is illustrated in a specific application scenario: Three anomalies are identified in the ground-penetrating radar data of the target area, and are sequentially labeled as node A (candidate for main nest, classification confidence 0.92, burial depth 0.8m), node B (candidate for ant trail, confidence 0.85), and node C (candidate for secondary nest, confidence 0.78). After data acquisition and classification, an anomaly spatial label dataset containing the three nodes and their three-dimensional coordinates and confidence labels is formed. During topology reconstruction, node A is used as the starting node for the main nest, and multi-path tracking search is used to establish an initial topological skeleton along the ABC direction. Self-consistency verification calculates the dimensionless deviation of each geometric parameter. The estimated equivalent diameter of node B's cross-section does not exceed the allowable range of the ant trail, and the total penalty for constraint violation is within the tolerance threshold. The skeleton passes verification, resulting in a three-dimensional connectivity graph of the nests on both sides of the three nodes. In the directed reachability determination, segment AB is a downward-sloping segment, statically reachable; the solidification rate coefficient β is calculated based on the wall thermal properties, and the flow duration T... flow The time required for fluid to traverse the segment is greater than the time required for fluid to pass through it, and the time-varying gating verification passes; the segment is marked as a reachable direction and added to the directed reachable flow graph. In the infusion scheme planning, node C is reached by only one path and is an exclusive node; its entire cavity volume is allocated to the corresponding infusion port. The final output includes the infusion port coordinates, allocation volume, and execution sequence of the multi-point infusion planning scheme. Since ABC is a single-chain path, the shortest and longest allowable interval time windows between the two infusion ports are calculated according to the above steps and arranged in the order of injection from end A first, waiting for the interval window, and then injection from end C.

[0201] Example 7: This example further details the process of pre-constructing the feature fingerprint database, specifically including:

[0202] Step 701: Obtain measured dielectric data of material samples under different moisture content gradients.

[0203] As an optional implementation method, the construction of the feature fingerprint database relies on the precise electromagnetic physical properties of the subsurface medium. In a laboratory environment, standard specimens for five types of subsurface media were prepared, specifically including scrapings from the walls of airborne ant tunnels, fungal garden fillers, compacted nest wall blocks, undisturbed soil surrounding nests, and woody root samples invading the nest area. For each type of sample, multiple uniform moisture content gradients were configured within a physical range from naturally air-dried moisture content to capillary saturation moisture content. Using a vector network analyzer with an open-ended coaxial probe, frequency sweep detection was performed on each of the configured specimens within the frequency band corresponding to the actual transmission frequency of the ground-penetrating radar, obtaining the real and imaginary parts of the complex permittivity. The obtained relative permittivity and loss factor were categorized and summarized according to medium type and moisture content gradient to obtain the measured dielectric data of the material samples.

[0204] Step 702: Based on the dielectric type corresponding to the measured dielectric data of the material sample and its internal subphase composition, construct a multi-component mixed dielectric model to characterize the arrangement of internal subphases of the material. The equivalent properties of the multi-component mixed dielectric model are jointly determined by the volume fraction and shape factor of each subphase.

[0205] Optionally, the dielectric constant of the undisturbed soil can be obtained by third-order fitting using a conventional empirical polynomial, while the internal dielectric constant of the air tunnel is directly assigned a near-vacuum constant. For the fungal garden filling and compacted nest walls unique to termite nests, due to their strong internal heterogeneity, conventional empirical models cannot accurately predict their properties. Therefore, this complex medium is considered a mixture composed of multiple subphases. Specifically, the fungal garden filling is divided into three subphases: mycelial organic matrix, mineral particles, and pore water; the compacted nest walls are divided into two subphases: compacted mineral matrix and pore water. A multi-component mixed dielectric model is constructed based on the volume fraction weighted complex refractive index mixing theory; its calculation formula is as follows:

[0206] ε eff α =Σ j (f j *ε j α );

[0207] Where, ε eff Σ is the equivalent relative permittivity, α is the shape factor describing the geometric arrangement of each subphase, and Σ is the Σ-phase. j f represents the summation over all subphases within a component. j Let ε be the volume fraction of the j-th subphase. j Let be the intrinsic relative permittivity of the j-th subphase. In this model, the change in the overall moisture content of the medium is converted into an increase or decrease in the volume fraction of the pore water subphase.

[0208] Step 703: Using the measured dielectric data of the material sample as the inversion target, perform error minimization solution on the multi-component mixed dielectric model to invert and calibrate the optimal value of the shape factor.

[0209] Correspondingly, traditional hybrid models typically use fixed empirical values ​​for the shape factor, which can lead to model truncation errors when dealing with complex biological nest materials. This embodiment uses the shape factor as a free parameter to be optimized, extracting the measured relative permittivity for the fungal garden and nest wall materials under multiple moisture content gradients. A least squares inversion mechanism is used to minimize the sum of squared deviations between the theoretically calculated equivalent relative permittivity and the measured relative permittivity. The calculation formula is as follows:

[0210] α opt =argmin α (Σ k ((ε eff (θ vk ,α)-ε meas_k ) 2 ));

[0211] Where, α opt argmin is the optimal shape factor. α Let Σ represent the variable that minimizes the objective function. k ε represents the summation of the total number of samples over the moisture content gradient. eff θ is the equivalent relative permittivity calculated based on the theoretical model. vk Let α be the volumetric water content parameter of the k-th sample, α be the shape factor to be optimized, and ε be the volumetric water content parameter. meas_k Let be the measured relative permittivity of the k-th sample.

[0212] Accordingly, the optimal parameters obtained will be embedded into the multi-component mixed dielectric model. Furthermore, an optional error compensation mechanism is provided. Measured dielectric data from reserved material samples are extracted for cross-validation, and the root mean square relative error is calculated. If the root mean square relative error exceeds a preset accuracy tolerance limit, a second-order moisture content correction term is superimposed on the original multi-component mixed dielectric model output; the calculation formula is as follows:

[0213] △ ε =c1*θ v +c2*θ v 2 ;

[0214] Among them, △ ε Here, c1 is the dielectric constant compensation factor, c2 is the first-order correction factor, and θ is the second-order correction factor. v This represents the volumetric water content.

[0215] Step 704: Perform forward modeling on a variety of preset nest configurations using the calibrated multi-component mixed dielectric model, extract response features to compile and generate a pre-built feature fingerprint library.

[0216] Specifically, a typical nest configuration with multiple variable combinations is pre-defined. These configuration variables include the burial depth of the main nest, the equivalent diameter of the main nest cavity, the inclination angle of the ant tunnels, the location of tree root disturbances, and the soil background moisture content. These variables are orthogonally combined to generate multiple independent test scenarios. Corresponding dielectric distribution areas are arranged in a two-dimensional computational grid according to the test scenarios, and the calibrated multi-component hybrid dielectric model is used to assign precise conductivity and dielectric constant to the computational grid.

[0217] Accordingly, the excitation waveform and antenna center frequency are set, and the Maxwell's equations in the computational grid are solved using the finite-difference time-domain method to simulate the propagation, reflection, and attenuation processes of electromagnetic waves. The simulated synthetic received signals are then stitched together to form a two-dimensional radar spectrum profile. For various anomalous regions in this two-dimensional radar spectrum profile, the aforementioned background adaptive normalization process is executed to extract multi-attribute features such as instantaneous peak amplitude ratio, dominant frequency offset, and polarity mode scalar. The arithmetic mean of the feature vectors of the same type of anomalous body under different configuration scenarios is calculated. The various mean feature vectors, category labels, and discrimination weight vectors determined by Fisher discriminant analysis are stored together to generate a pre-built feature fingerprint library for online classification.

[0218] Example 8: Based on the above examples, this example provides a method for field execution and non-destructive verification of a multi-point injection planning scheme. In one possible implementation, after obtaining the multi-point injection planning scheme, the following steps are further included:

[0219] As an optional implementation method, multi-point high-temperature liquid aluminum injection operations are performed in the target area according to the multi-point injection planning scheme.

[0220] Specifically, the ground spatial coordinates of each injection port in the multi-point injection plan are extracted. Working holes are drilled to the target depth at the corresponding ground spatial coordinates in the target area to establish a physical connection between the surface injection equipment and the underground medium channel. The execution sequence with time interval window constraints arranged in the multi-point injection plan is obtained. According to the order specified in the execution sequence, a predetermined amount of liquid aluminum is injected into each injection port within the corresponding time window. After the injection operation is completed, the physical state of the site is kept stable, and the liquid aluminum is allowed to dissipate heat and solidify in the underground medium.

[0221] As an optional implementation, ground-penetrating radar resurvey data is collected at the same survey line location as the one used to acquire ground-penetrating radar detection data.

[0222] Further, the spatial layout parameters and hardware configuration parameters of the survey line recorded during the initial acquisition of ground-penetrating radar (GPR) data are obtained. After the liquid aluminum solidifies, the GPR equipment is redeployed in the target area according to the spatial layout parameters of the survey line. Using the same center frequency and sampling interval settings as the initial detection, continuous electromagnetic wave scanning is performed along the original survey line trajectory to acquire the original re-measurement signal sequence. The acquired original re-measurement signal sequence undergoes the same standardized preprocessing operations as the initial data. The preprocessing operations specifically include zero-time correction, background clutter removal, and gain compensation. After preprocessing, spatially aligned GPR re-measurement data is output.

[0223] As an optional implementation method, the numerical deviation between the ground-penetrating radar re-measurement data and the ground-penetrating radar detection data is calculated to generate a difference profile, and the filling integrity evaluation result is output by combining the three-dimensional structure connectivity diagram of the nest.

[0224] Optionally, aluminum, as a high-conductivity medium, will strongly reflect electromagnetic waves after being filled into the underground cavity, and its echo characteristics will have an absolute physical difference from the original cavity medium response before filling. Ground-penetrating radar (GPR) re-measurement data with strictly aligned spatial positions and initial GPR detection data are extracted, and subtraction is performed one by one according to the survey line number and sampling time point to generate a difference profile characterizing the change in electromagnetic response. In the generated difference profile, the area successfully filled with liquid aluminum will appear as a high-amplitude concentrated anomaly zone. The three-dimensional spatial distribution boundary of this high-amplitude concentrated anomaly zone is extracted and spatially superimposed and compared with the three-dimensional spatial positions of each chamber and ant tunnel recorded in the three-dimensional structural connectivity map of the nest.

[0225] Based on the spatial overlay comparison results, the actual volume of areas exhibiting high-amplitude anomalies was statistically analyzed, and the filling rate index was calculated; the calculation formula is as follows:

[0226] R fill =V filled / V total ;

[0227] Among them, R fill V is the filling rate indicator. filled V represents the estimated filling volume mapped by high-amplitude anomaly regions in the difference profile to the 3D connectivity graph of the nest structure. total The overall network cavity volume is recorded in the three-dimensional connectivity graph of the nest structure.

[0228] Furthermore, in an optional implementation, areas marked as usable cavities in the three-dimensional connectivity diagram of the nest but where no high-amplitude anomaly response is observed at the corresponding spatial location in the differential profile are identified as unfilled dead zones. The three-dimensional spatial coordinates of these unfilled dead zones are extracted and packaged together with the calculated filling rate index, and then output. This filling integrity evaluation result, through a differential calculation mechanism using preceding and following radar data, enables the objective calculation of the fluid filling quality of the underground three-dimensional pipe network using non-destructive methods without large-area destructive excavation.

[0229] This invention introduces an adaptive multi-path tracking search mechanism, which transforms global constraints into local segment-by-segment deflection verification, enabling high-fidelity topology extraction of complex curved meshes and solving the problem of path breakage caused by high curvature channels that is difficult to track by conventional detection algorithms.

[0230] Furthermore, for structural distortions caused by noise interference, a dimensionless squared penalty function is used to verify the self-consistency of geometric features, and a local backtracking mechanism is used to eliminate false connectivity, reducing topological logic fallacies in the reconstructed network. To address the issue that traditional static planning neglects the time-varying decay of fluid flow, heat transfer and fluid mechanics are integrated, a time-varying flow decay model is introduced to establish time-varying reachability determination, and a flow resistance-weighted allocation model and time window constraints are used to arrange the timing sequence.

[0231] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A termite nest irrigation planning method based on ground-penetrating radar and physical flow characteristics, characterized in that, include: The ground-penetrating radar detection data of the target area is acquired, normalized multi-attribute features are extracted, and anomalies are classified based on a pre-built feature fingerprint database to obtain an anomaly spatial label dataset. Combining pre-set biological morphological prior constraints, topological association search and self-consistency verification correction are performed in the anomaly spatial label dataset to reconstruct the three-dimensional structure connectivity graph of the nest. Based on the physical flow characteristics, the directed reachability of paths in the three-dimensional structure connectivity graph of the nest is determined, and a directed reachable flow graph is constructed. Based on the directed reachable flow graph, the set of infusion ports is selected according to the maximum coverage principle, the allocation volume is calculated and the execution sequence is arranged to obtain a multi-point infusion planning scheme. The construction of a directed reachable flow graph includes: Extract the path segments from the 3D connected graph of the nest structure; Based on the static properties, the gravitational potential difference and capillary resistance of each path segment are evaluated to determine the candidate segments that meet the static reachability conditions. Based on solidification kinetics, time-varying gating constraint verification is performed on candidate segments to verify whether the fluid can traverse the entire length of the segment before it completely solidifies and blocks it. The path segments that simultaneously satisfy both static reachability conditions and time-varying gating constraint verification are used as reachable directions to construct a directed reachable flow graph; Among them, the time-varying gating constraint verification of candidate segments based on solidification kinetics includes: Based on the thermophysical properties of the wall material corresponding to the candidate segments, the solidification rate coefficient characterizing the fluid's solidification migration from the pipe wall to the center is obtained; Based on the solidification rate coefficient, the equivalent diameter of the cross section of the candidate segment, and the minimum channel threshold required for fluid flow, the flow duration corresponding to the candidate segment is calculated. By combining the inclination angle of the candidate segments with the time-varying reduction characteristics of the fluid velocity, the thermally permeable distance that the fluid can reach during the flow duration is estimated. Compare the actual segment length of the candidate segment with the thermally passable distance. If the actual segment length does not exceed the thermally passable distance, the candidate segment is determined to have passed the time-varying gating constraint check. A two-stage adaptive scanning strategy was used to acquire this data. In the first stage, a sparse grid of ground-penetrating radar lines was used for initial scanning, and sections with reflected energy higher than the background mean were marked by the sliding window energy detection method. In the second phase, a more directional and encrypted scan is conducted on the marked sections using denser survey lines to obtain raw detection data with high spatial resolution.

2. The method according to claim 1, characterized in that, Extract normalized multi-attribute features, including: The anomalous body region is identified from the ground-penetrating radar detection data, and an exclusion buffer zone is delineated outward from the anomalous body region as a reference. An annular candidate window is defined in the outer region of the exclusion buffer zone, and background signal samples are extracted within the annular candidate window; Calculate the attribute variation coefficient of the background signal sample, and perform a homogeneity test on the annular candidate window in combination with a preset homogeneity threshold; If the homogeneity test fails, the range of the annular candidate window is expanded outward, and sampling and testing are repeated until the homogeneity test is passed. Background reference values ​​are calculated based on background signal samples within the annular candidate window that have passed the homogeneity test. By using background reference values ​​to perform benchmark conversion on the original attribute values ​​of the anomaly region, normalized multi-attribute features are obtained.

3. The method according to claim 1, characterized in that, Performing a topological association search specifically includes: Determine the primary nest starting node from the anomaly spatial label dataset; Based on the starting node of the main nest, and combined with the path angle deflection constraints and cumulative length constraints specified in the biological morphology prior constraints, nodes are searched and associated in the anomaly spatial label dataset to construct the initial topological skeleton.

4. The method according to claim 3, characterized in that, The search and association of nodes employs a multi-path tracing search mechanism, including: Construct a path seed with the main nest starting node and assign it an initial direction of movement, then place it into a priority queue sorted by confidence. Extract active paths from the priority queue and retrieve candidate nodes within the neighborhood of their end nodes; For candidate nodes, perform local direction deflection verification and cumulative length verification. Add candidate nodes that meet the path angle deflection limit and cumulative length limit as new branches to the current path, update the local forward direction, and put them back into the priority queue. Complete paths are extracted based on preset path termination conditions and node affiliation conflict resolution rules to construct the initial topology skeleton.

5. The method according to claim 1, characterized in that, The calculation of the allocation amount specifically includes: For exclusive nodes that can be reached independently from a single injection port in a directed reachable flow graph, their entire cavity volume is included in the total allocation of the corresponding injection port. For a shared node that is simultaneously reached by multiple injection ports through different paths, each path is treated as a series pipe and the flow resistance corresponding to each path is calculated separately. The volume allocation coefficient is determined based on the reciprocal proportion of the flow resistance corresponding to each path, and the cavity volume of the shared node is allocated to the corresponding injection port according to the volume allocation coefficient, so as to obtain the allocation amount of each injection port.

6. The method according to claim 1, characterized in that, Arrange the execution sequence, including: Identify adjacent injection ports with shared path segments based on directed reachable flow graphs; Obtain the longest permissible interval between the fluid injected from the preceding inlet solidifying into a blockage in the shared path segment; Obtain the shortest permissible interval time required for the fluid injected from the preceding infusion port to flow to the far end of the shared path segment; By combining the longest and shortest allowed interval times, execution timing with time interval window constraints is arranged for adjacent injection ports.

7. The method according to claim 1, characterized in that, The pre-built feature fingerprint library is obtained through the following steps: Obtain measured dielectric data of material samples under different moisture content gradients; A multi-component mixed dielectric model is constructed to characterize the arrangement of subphases within the material. The equivalent properties of the multi-component mixed dielectric model are jointly determined by the volume fraction and shape factor of each subphase. Using the measured dielectric data of material samples as the inversion target, the error minimization solution is performed on the multi-component mixed dielectric model to invert and calibrate the optimal value of the shape factor; Forward modeling of various preset nest configurations is performed using a calibrated multi-component mixed dielectric model, and response features are extracted to compile a pre-built feature fingerprint library.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored program, wherein the program, when executed, performs the method of any one of claims 1 to 7.

Images