Methods for estimating termite density and inferring termite tracts in earth-rockfill dams

Abstract

The present invention discloses a method for estimating termite density and inferring termite paths in earth-rockfill dams, belonging to the technical field of earth-rockfill dam safety monitoring. The method comprises: processing acquired multi-source monitoring data to generate a standardized multi-source data set; performing phase-preserving parallel frequency-domain time mapping on dielectric constant time series data to generate time-synchronized data that retains transient characteristics; utilizing the time-synchronized data and resistivity spatial distribution data to extract instantaneous density response characteristics that characterize the intensity of termite activity and termite path branching probability characteristics that indicate the tendency of termite path branching; constructing and solving the objective function of density-path mutual feedback based on the two types of characteristics to obtain the instantaneous termite density distribution and termite path branching prediction results. The present invention realizes the coordinated optimization of multi-physical field information through the time mapping and density-path mutual feedback inversion framework, significantly improving the accuracy of density estimation and the reliability of termite path inference.

Description

Technical Field

[0001] The present invention relates to a safety monitoring technology for earth-rock dams, in particular to a method for estimating termite density and inferring ant paths in earth-rock dams. Background Art

[0002] Termites are among the oldest living organisms on Earth. They are known for their elusive nature, deep nests, difficulty in detection, and significant safety hazards to earth-rockfill dams. Therefore, termite control has long been a focus and a challenge in engineering research. Termites nest densely within dams and reproduce rapidly, creating extensive tunnels that damage the dam's internal structure and can easily lead to leakage, piping, sinkholes, and even dam collapse.

[0003] Currently, nondestructive detection of termites in earth-rockfill dams relies primarily on the integrated application of multiple geophysical exploration methods. For example, the high-density resistivity method measures the resistivity distribution within the dam to identify low-resistance anomalies caused by ant tunnels and ant nests. This method provides relatively clear spatial structural information and is suitable for macroscopic delineation of ant tunnel networks. Furthermore, the geological radar (GPR) method utilizes the propagation characteristics of high-frequency electromagnetic waves to detect underground targets by analyzing differences in dielectric constant. It is sensitive to changes in soil moisture and porosity caused by termite activity and can indicate signs of termite activity. At the data processing level, when combining multiple detection methods, methods such as time window interpolation or sliding average are often used to align monitoring data with different sampling rates for comprehensive analysis. In actual engineering applications, experienced technicians typically perform manual interpretation and comprehensive analysis based on the anomaly patterns of different physical fields to roughly determine the area of ​​termite infestation and its development trend.

[0004] However, existing technologies still face many challenges in achieving high-precision and quantitative assessment of termite activity, including:

[0005] Existing time synchronization methods (such as interpolation or averaging) are inherently lossy when aligning monitoring data at different frequencies. They irreversibly smooth or eliminate key amplitude and phase information in the high-frequency dielectric pulse signals that characterize transient termite activity. This loss of key dynamic features leads to significant inaccuracy in quantitative assessments of termite activity intensity, making it impossible to accurately capture the true extent of the instantaneous damage.

[0006] Existing analytical methods lack inherent synergistic constraints and typically rely on a simple overlay and comparison of individual inversion results, leading to a disconnect between density estimation and path inference. This inability to achieve a closed-loop optimization where density drives the path and the path constrains the density leads to physically inconsistent results and often produces spurious anomalies (e.g., false density peaks or illogical paths) that contradict the actual geology and activity. Summary of the Invention

[0007] The purpose of the invention is to provide a method for estimating termite density and inferring ant tunnels in earth-rock dams in order to solve the technical problems existing in the prior art.

[0008] The technical solution is a method for estimating termite density and inferring termite tracts in earth-rock dams, including:

[0009] Process acquired multi-source monitoring data to generate standardized multi-source data sets, including dielectric constant time series data, resistivity spatial distribution data, alarm records, and geological parameters;

[0010] Perform phase-preserving parallel frequency-domain time mapping on the dielectric constant time series data to generate time-synchronized data;

[0011] Using time-synchronized data and resistivity spatial distribution data, we extracted the instantaneous density response characteristics that characterize the intensity of termite activity and the ant trail branching probability characteristics that indicate the tendency of ant trail branching.

[0012] Based on the instantaneous density response characteristics and the ant trail branching probability characteristics, the objective function of density-path feedback was constructed and solved to obtain the instantaneous termite density distribution and ant trail branching prediction results.

[0013] Preferably, performing phase-preserving parallel frequency-domain time mapping comprises:

[0014] Parse the dielectric constant time series data into multiple independent frequency bands to obtain frequency band dielectric data;

[0015] In the process of parallel time mapping of the frequency-band dielectric data, the instantaneous phase difference between each frequency band is detected and quantified in real time;

[0016] When the instantaneous phase difference deviates from the preset coupling relationship, the phase of each frequency band is dynamically compensated and corrected to generate time-synchronized multi-band coupling data and phase-related characteristic parameters; and the two (multi-band coupling data and phase-related characteristic parameters) are used to extract the instantaneous density response characteristics.

[0017] Preferably, the steps of constructing and solving the objective function of density-path mutual feedback include:

[0018] Based on the instantaneous density response characteristics and the ant path branching probability characteristics, a joint objective function is constructed that integrates the dielectric inversion consistency, the density gradient-resistivity curl coupling, and the path topology rationality.

[0019] An iterative optimization algorithm is used to solve the joint objective function, and each iteration includes interaction and bidirectional constraint mechanisms;

[0020] Based on the current density distribution estimation results, spatial weights are applied to the ant path branch prediction to guide the generation of paths;

[0021] Based on the updated ant path branch prediction results, topological constraints are imposed on the density distribution estimation in reverse order to correct the spatial distribution of density.

[0022] Continue executing until convergence to obtain the instantaneous termite density distribution and ant path branch prediction results.

[0023] Preferably, based on the instantaneous density response characteristics and the ant path branching probability characteristics, an objective function of density-path mutual feedback is constructed, including:

[0024] The instantaneous density response characteristics and ant path branching probability characteristics are spatially aligned in the same coordinate system, so that each grid point in the coordinate system carries both density and probability information.

[0025] The density-gradient field coupling operator is called to convolute the spatially aligned instantaneous density response features with the ant path branching probability features through a preset spatiotemporal correlation function. A density-path coupling feature set is generated that integrates the dielectric field and resistivity field information. This (density-path coupling feature set) is used to construct a joint objective function.

[0026] Preferably, the process of parallel time mapping includes:

[0027] The dielectric loss tangent value of the soil layer corresponding to each frequency band is extracted from the standardized multi-source data set;

[0028] An independent nonlinear time mapping function is constructed for each frequency band, and the compression degree of the time mapping is modulated by the extracted dielectric loss tangent value;

[0029] Dynamically calculate adaptive weights for each frequency band based on soil depth;

[0030] Adaptive weights are used to perform weighted combination of the independent mapping results of each frequency band to generate a time mapping result that adapts to the heterogeneity of the soil layer and is used for real-time detection of instantaneous phase differences.

[0031] Preferably, the step of extracting the ant trail branching probability feature indicating the ant trail branching tendency includes:

[0032] Based on the resistivity spatial distribution data, the three-dimensional gradient vector field representing the spatial change rate of the physical field is calculated;

[0033] Perform a curl operation on the three-dimensional gradient vector field to obtain a curl field that describes the local rotation strength of the gradient field;

[0034] By establishing a quantitative correlation model between the intensity of the curl field and the branching tendency of the ant path, the ant path branching probability characteristics are generated.

[0035] Preferably, the joint objective function includes the following constraints:

[0036] Dielectric inversion consistency constraints are used to ensure that the final density estimation results are consistent with the dielectric response analysis results;

[0037] Cross-physics coupling constraints to mathematically relate the density space gradient to the resistivity gradient field curl;

[0038] Topological rationality regularization term, which is used to impose structural rationality penalties on the ant path network to avoid generating illogical paths;

[0039] Time consistency constraints are used to maintain the smooth continuity of the solution in the time series.

[0040] Preferably, parallel time mapping of the frequency-band dielectric data is achieved by a parallel computing architecture, which includes:

[0041] An independent parallel processing unit set for each frequency band, responsible for performing the time mapping calculation of the corresponding frequency band;

[0042] A shared memory mechanism for data exchange between multiple parallel processing units;

[0043] A synchronization barrier is used to ensure that all parallel processing units complete the calculation before entering the next step, thereby eliminating the cumulative delay of serial calculation.

[0044] Preferably, dynamic compensation correction is performed on the phase of each frequency band, including:

[0045] Comparing the instantaneous phase difference detected in real time with a preset phase difference threshold;

[0046] When the absolute value of the instantaneous phase difference exceeds the threshold, phase compensation is started;

[0047] In addition, phase compensation is achieved through a correction algorithm based on the hyperbolic tangent function.

[0048] Preferably, the weight coefficients of the constraints in the joint objective function are determined by an adaptive method, specifically:

[0049] The L-curve method is used to find the optimal balance between the smoothness and fitting error of the solution under different weight coefficient combinations.

[0050] The weight coefficients corresponding to the optimal balance point are combined as the optimal weight in the joint objective function.

[0051] Preferably, the solution is obtained by iterative optimization, including:

[0052] In each iteration, the density distribution and path probability are optimized and updated alternately;

[0053] Each optimization update is performed by calculating the current gradient and combining it with the previous search direction weighted by the conjugate coefficient to determine the current optimal search direction.

[0054] Along the optimal search direction, the optimal step size is determined by line search to complete the update of the density distribution or path probability.

[0055] Preferably, the method further includes a dynamic update and error correction mechanism, specifically:

[0056] When new monitoring data is received, the prediction error between the prediction results output by the current model and the actual measured conditions reflected by the new data is quantified;

[0057] When the prediction error exceeds a dynamic threshold, adaptive correction is initiated:

[0058] The key parameters in the joint objective function are updated, and the joint inversion solution is re-executed based on the updated parameters to generate error-corrected termite instantaneous density distribution and ant path branching prediction results.

[0059] Preferably, dynamic compensation correction is performed on the phase of each frequency band, and predictive compensation is also included, specifically:

[0060] A phase change trend prediction model is established. This model predicts the phase change trend at the next moment based on the compensated and corrected historical phase sequence at multiple moments in the past.

[0061] According to the predicted phase change trend, the coupling parameters used in the compensation correction are adjusted in advance to achieve forward-looking phase compensation.

[0062] Beneficial effect: Through the innovative time mapping and density-path mutual feedback inversion framework, the collaborative optimization of multi-physical field information is achieved, which significantly improves the accuracy of density estimation and the reliability of ant path inference. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 It is a flow chart of the present invention.

[0064] Figure 2 4 is a flow chart of performing phase-preserving parallel frequency-domain time mapping according to the present invention.

[0065] Figure 3 It is a flow chart of the objective function of constructing and solving density-path mutual feedback in the present invention.

[0066] Figure 4 It is a flow chart of the objective function of constructing density-path mutual feedback in the present invention. DETAILED DESCRIPTION

[0067] In order to make the purpose, technical solutions and advantages of the present invention clearer, the following Figures 1 to 4 The present invention is further described in detail. First, the applicant has conducted an in-depth analysis of the problems existing in the prior art. The analysis process is not found in the existing literature and is as follows:

[0068] Research has found that the dielectric constant response, which characterizes termite activity such as instantaneous gnawing and movement, is a high-frequency (kHz-level) pulse signal rich in transient characteristics; whereas the resistivity distribution, which characterizes the macroscopic spatial structure of termite tunnels, is a low-frequency (minute-level) slowly varying signal. Existing synchronization methods, such as time window averaging or linear interpolation, employ lossy compression. In the process of aligning high-frequency data with low-frequency data, they irreversibly smooth or even eliminate transient pulse amplitude and phase information, which is crucial for quantitatively assessing activity intensity. This leads to a significant underestimation of termite activity intensity, making the so-called synchronization effectively at the expense of the most critical dynamic characteristics and unable to provide effective data support for subsequent precise quantitative inversion.

[0069] Furthermore, existing joint analyses often rely on qualitative superposition and comparison of results, lacking a mathematical framework capable of intrinsically coupling information from different physical fields. For example, even if a dielectric anomaly signal is detected within a low-resistivity anomaly, this density information cannot be used to reverse-correct or constrain the inference of the ant path. Similarly, the inferred ant path topology cannot constrain the rationality of the spatial distribution of density estimates. This one-way, loose coupling approach results in the final interpretation results being physically inconsistent, often leading to problems such as false density peaks or isolated paths unrelated to actual activity, severely limiting the accuracy and reliability of the detection results.

[0070] To this end, the following solutions are provided:

[0071] First, the main characters used in the relevant embodiments are as follows:

[0072] ε' is the real part of the dielectric constant; ε'' is the imaginary part of the dielectric constant; E is the resistivity; ρ is the instantaneous density of termites; P is the probability of ant tunnel branching; t raw is the original sampling time; t sync is the mapping time after synchronization; f is the signal frequency; tan(δ) is the dielectric loss tangent; α i is the time mapping compression coefficient of the i-th frequency band; w i (z) is the adaptive weight of the i-th frequency band at depth z; φ is the signal phase; Δφ ij is the instantaneous phase difference between frequency band i and frequency band j; β is the phase compensation intensity coefficient; θ thresholdis the phase difference compensation threshold; grad is the Hamiltonian operator (used to calculate gradient and curl); J is the joint objective function; λ, μ, γ are the weight coefficients of each constraint in the objective function; Ω (topology) is the topological rationality regularization term; R (consistency) is the temporal consistency constraint term; M density is the density space weight mask matrix; M path is the path topology mask matrix.

[0073] Example 1: A method for estimating termite density and inferring termite tracts in an earth-rock dam, the method comprising the following steps:

[0074] Step S100 , processing the acquired multi-source monitoring data to generate a standardized multi-source data set, including dielectric constant time series data, resistivity spatial distribution data, alarm records and geological parameters.

[0075] Standardized multi-source data sets refer to converting raw monitoring data from different sources, sampling frequencies, and time and space references into a data set with a unified timestamp, unified spatial coordinate system, and standard data format through a series of preprocessing. In this embodiment, the specific data processing process includes:

[0076] Raw dielectric constant time series data (including the real and imaginary parts ε'', with a sampling frequency of 1-10kHz) were read from the earth-rock dam monitoring system. Data integrity was checked and outliers were removed (for example, using the 3σ criterion to identify and replace outliers with cubic spline interpolation). At the same time, resistivity spatial distribution data (spatial resolution 0.5m) and monitoring point alarm records were obtained and formatted in a standardized manner.

[0077] The soil geological exploration data was read to extract parameters such as the dielectric loss tangent tan(δ), soil moisture content, and density of each depth layer. These parameters were calibrated through field calibration experiments in known termite activity areas to establish a calibration geological parameter library.

[0078] Using GPS timestamp alignment algorithm and spatial coordinate transformation, all data are mapped to a unified spatiotemporal reference system and the final standardized multi-source dataset is output.

[0079] This step aims to address the data fusion challenge caused by the heterogeneity of multi-source monitoring data. Without this step, the temporal and spatial mismatches between different physical quantities will make subsequent coupled analysis and joint inversion impossible or result in significant errors.

[0080] Step S200 , performing phase-preserving parallel frequency-domain time mapping on the dielectric constant time series data to generate time-synchronized data.

[0081] Unlike traditional interpolation or averaging downsampling, the core of phase-preserving parallel frequency-domain time mapping is to decompose the high-frequency signal into multiple frequency bands for parallel processing. While performing nonlinear time-scale compression, it also uses a dynamic compensation mechanism to maximize the preservation of the relative phase relationship between the frequency bands. The specific implementation method will be detailed in Example 2.

[0082] This step aims to address the three- to four-order-of-magnitude timescale discrepancy between the dielectric constant data (kHz-level) and the resistivity data (minute-level). Traditional interpolation methods lose the phase information and transient characteristics of the pulse signals, which are crucial for characterizing transient termite activity. This step, through a phase-preserving mechanism, achieves time synchronization while preserving the temporal accuracy of density fluctuations, enabling accurate estimation of instantaneous density.

[0083] Step S300 : Using the time synchronization data and the resistivity spatial distribution data, extract the instantaneous density response characteristics representing the intensity of termite activity and the ant trail branching probability characteristics indicating the tendency of ant trail branching.

[0084] The instantaneous density response characteristic is a characteristic quantity calculated by analyzing the time-synchronized dielectric response data, which quantitatively represents the intensity of termite activity at a specific point in time and space. The ant trail branching probability characteristic is a probability value calculated by analyzing the spatial gradient field and curl field of the resistivity data, which represents the probability of the occurrence of an ant trail branch at that spatial point.

[0085] In this embodiment, it mainly includes two processes: extracting instantaneous density response features and ant path branch probability features (which will be described in detail in the fourth embodiment).

[0086] Extracting the instantaneous density response characteristics mainly uses time-synchronized multi-band coupling data through a calibrated relationship model (for example, ρ(t, ω)=Σ[A i ×ε' i (t, ω)×cos(ω×t+φ corrected-i )]), calculate the instantaneous density response value. Among them, ρ(t, ω) represents the instantaneous density response characteristic, which means the instantaneous termite activity density estimated by the dielectric response model at time t and angular frequency ω. i The response amplitude coefficient of the i-th frequency band is obtained by measuring in a calibration area with known termite density and performing multiple linear regression calibration. It reflects the correlation strength between the dielectric response in a specific frequency band and the termite density.

[0087] ε' i (t, ω) represents the real part of the dielectric constant of the i-th frequency band at time t and angular frequency ω, which is one of the core input data obtained from the monitoring equipment and pre-processed. corrected-irepresents the phase of the i-th frequency band after dynamic compensation correction. The specific calculation method is described in detail in sub-step S213 of the second embodiment, which is calculated by φ corrected-i =φ i +β×tanh(Δφ ij / θ threshold ) is calculated. t is time and ω is angular frequency.

[0088] This step aims to extract key information from different physical fields that can be correlated and verified. The dielectric response provides high-temporal-resolution information on activity intensity, while the resistivity field provides high-spatial-resolution information on structural orientation. Extracting these two types of features is a prerequisite for the subsequent density-path feedback (density-path bidirectional constraint) inversion.

[0089] Step S400 : Based on the instantaneous density response characteristics and the ant trail branch probability characteristics, an objective function of density-path interaction is constructed and solved to obtain the instantaneous termite density distribution and ant trail branch prediction results.

[0090] The objective function of the density-path feedback loop is a mathematical optimization model. It not only incorporates data fitting and regularization terms for multiple physical fields but, more importantly, establishes a closed-loop feedback optimization mechanism where density drives the path and the path constrains the density. The specific construction and solution process will be detailed in Example 3.

[0091] To address the issues of insufficient accuracy and poor robustness of single physical field detection, a mutual feedback mechanism is used. Density information provides real-time spatial weights for path prediction, while the path topology, in turn, provides reasonable spatial constraints for density estimation. This achieves a 1+1>2 technical fusion effect, resulting in more accurate and reliable density distribution and ant path prediction results than any single method.

[0092] Embodiment 2: This embodiment is a detailed explanation of step S200 in embodiment 1, and its purpose is to solve the problem of huge time scale differences between different monitoring signals, while avoiding the loss of key transient feature information in traditional methods.

[0093] In step S210, the dielectric constant time series data is parsed into multiple independent frequency bands to obtain frequency-band dielectric data. During the parallel time mapping of the frequency-band dielectric data, the instantaneous phase difference between the frequency bands is detected and quantified in real time. When the instantaneous phase difference deviates from the preset coupling relationship, the phase of each frequency band is dynamically compensated and corrected to generate time-synchronized multi-band coupling data and phase-related characteristic parameters. These two (multi-band coupling data and phase-related characteristic parameters) are used to extract the instantaneous density response characteristics.

[0094] Phase preservation means that when performing time-scale transformation on a signal, not only the timestamp of the signal is adjusted, but also the relative phase relationship between the frequency components is actively monitored and maintained without distortion, which is crucial for preserving the transient waveform characteristics of the signal.

[0095] Termite gnawing and movement, among other activities, generate non-periodic, short-lived pulses in the dielectric constant signal. Traditional time synchronization methods, such as linear interpolation, average these pulses, resulting in loss of phase information and making it impossible to identify the specific moment and intensity of the activity. This step, through phase-preserving nonlinear time mapping in the frequency domain, precisely compresses the high-frequency signal's time axis to align with the low-frequency signal, while preserving these critical transient pulse characteristics.

[0096] Furthermore, in order to implement the above steps, the following steps are further performed:

[0097] Step S211, the parallel time mapping process, includes: extracting the dielectric loss tangent value of the soil layer corresponding to each frequency band from a standardized multi-source data set; constructing an independent nonlinear time mapping function for each frequency band, wherein the compression degree of the time mapping is modulated by the extracted dielectric loss tangent value; dynamically calculating an adaptive weight for each frequency band according to the soil layer depth; using the adaptive weight to perform a weighted combination of the independent mapping results of each frequency band to generate a time mapping result that adapts to the heterogeneity of the soil layer and is used for real-time detection of instantaneous phase difference.

[0098] Specifically, it includes steps such as extracting frequency band decomposition and parameters, constructing nonlinear time mapping functions, and combining adaptive weights, such as:

[0099] The dielectric constant time series data is subjected to fast Fourier transform (FFT) and decomposed into n characteristic frequency bands (for example, n=8) at logarithmic intervals. At the same time, the dielectric loss tangent value tan(δ i (z)).

[0100] Construct an independent nonlinear time mapping function for each frequency band i, the specific form is: sync-i =t raw ×[1+α i ×(log(ε i ′ (f i )+1)) / (tan(δ i )+ε reg )];

[0101] Where: t sync-i Indicates the synchronization time of the i-th frequency band after mapping. raw Indicates the original sampling time. α iis the time mapping compression coefficient of the i-th frequency band.

[0102] This coefficient is preferably obtained by the empirical formula α i =α base ×(σ i / μ i )×[1+0.1×log(f i / f ref )] is determined, where α base is the base coefficient (e.g. 0.15), σ i and μ i are the standard deviation and mean of the dielectric strength in this frequency band, f ref is the reference frequency (e.g. 1kHz). ε' i (f i ) is the real part of the dielectric constant of the i-th frequency band. tan(δ i ) is the dielectric loss tangent value corresponding to the i-th frequency band. ε reg is a small regularization parameter (e.g. 0.001) to prevent the denominator from being zero.

[0103] Calculate the adaptive weight w of each frequency band at depth z i (z), for example, w i (z)=tan(δ(z)) -0.5 ×exp(-z / λ i ), where λ i is the frequency band attenuation coefficient. The final time mapping result t _sync is the mapping result of each frequency band t _sync-i Weighted combination of: t sync =Σ[w i (z)×t sync-i ].

[0104] In some optional implementations, the nonlinear time mapping function may also take other forms, such as based on an exponential function or a power function, as long as it can reflect the nonlinear relationship between the dielectric response intensity and the degree of time compression.

[0105] Step S212, performing parallel time mapping on the frequency-band dielectric data, is implemented through a parallel computing architecture, which includes: an independent parallel processing unit set for each frequency band, responsible for performing the corresponding frequency band time mapping calculation; a shared memory mechanism for exchanging data between multiple parallel processing units; and a synchronization barrier for ensuring that all parallel processing units complete the calculation before entering the subsequent steps, so as to eliminate the cumulative delay of serial calculation.

[0106] For example, on a GPU supporting CUDA or OpenCL, a computing core (or thread block) can be launched for each of the eight frequency bands. Each core reads the raw dielectric data and geological parameters from shared memory, independently performs the calculations in sub-step S211, and writes the results back to the shared memory. After all cores complete the calculations, a synchronization barrier (e.g., the _syncthreads() instruction) is used to ensure data consistency before subsequent phase compensation and weighted combination are performed.

[0107] This parallel architecture aims to address the computational latency issues associated with multi-band processing. If serial computation is used, the total latency is the sum of the latency of each band, which would not meet real-time requirements. This parallel architecture reduces the total computation time to the maximum processing time of a single band, significantly improving algorithm execution efficiency.

[0108] Step S213, dynamically compensating and correcting the phase of each frequency band, including: comparing the instantaneous phase difference detected in real time with a preset phase difference threshold; when the absolute value of the instantaneous phase difference exceeds the threshold, starting phase compensation; wherein, phase compensation is achieved through a correction algorithm based on the hyperbolic tangent function.

[0109] In some embodiments, the process includes real-time monitoring of phase difference, triggering compensation, and phase correction. The details are as follows:

[0110] Calculate the instantaneous phase difference Δφ between any two frequency bands i and j ij (t)=φ i (t)-φ j (t).

[0111] |Δφ ij (t)|with a phase difference threshold θ threshold The threshold is preferably adaptive, for example, θ threshold (t)=θ base +k×std(Δφ ij ), where θ base is the base threshold (e.g., π / 6), and k is the standard deviation coefficient (e.g., 0.5).

[0112] When |Δφ ij (t)|>θ threshold When , the phase of one of the frequency bands is corrected, the specific form is:

[0113] φ corrected-i =φ i +β×tanh(Δφ ij / θ threshold ); where: φ correctedi :Indicates the corrected phase. φ iRepresents the original phase. β is the phase compensation strength coefficient, an empirical parameter with a value range of [0.1, 0.5], for example, 0.3. The principle of choosing this value is to ensure rapid phase convergence while avoiding excessive oscillation. tanh(): This is the hyperbolic tangent function, which provides a bounded, nonlinear correction. When the phase difference is much larger than the threshold, the correction tends to saturate, preventing overmodulation.

[0114] The specific calculation process can be found in step S2 of Example 6 below.

[0115] Step S214, dynamically compensates and corrects the phase of each frequency band, and also includes predictive compensation, specifically: establishing a phase change trend prediction model, which is based on the historical phase sequence after compensation and correction at multiple moments in the past to predict the phase change trend at the next moment; and according to the predicted phase change trend, adjusting the coupling parameters used in the compensation correction in advance to achieve forward-looking phase compensation.

[0116] In some embodiments, the prediction model may be a weighted autoregressive model, such as:

[0117] Δφ pred (t+1)=Δφ(t)+Σ[w k ×(Δφ(t-k+1)-Δφ(tk))], where Δφ(t) is the phase difference at time t, w k is the time weight that decays exponentially with time k.

[0118] According to the predicted Δφ pred (t+1), the compensation intensity coefficient β or threshold θ in sub-step S213 can be adjusted in advance threshold , achieving active and proactive phase drift correction.

[0119] This step aims to address the lag issue in phase compensation. When a termite colony moves rapidly, phase changes occur rapidly, and compensation based solely on real-time monitoring can introduce delays. Predictive compensation anticipates future phase changes and enables proactive correction, maintaining extremely high synchronization accuracy even in high-speed, dynamic scenarios.

[0120] Through the above steps, this embodiment describes in detail how to achieve time synchronization under phase preservation. The time-synchronized data finally output is not only temporally aligned with the low-frequency resistivity data, but also retains the transient characteristics of the high-frequency dielectric signal that are critical for termite activity analysis.

[0121] Example 3: This example is a detailed description of step S400 in Example 1. Its purpose is to integrate multi-physics field information and obtain more accurate and reliable termite density and ant path results than any single physics field analysis through an optimization framework with physical constraints.

[0122] Step S410 constructs the objective function of density-path mutual feedback based on the instantaneous density response characteristics and the ant path branching probability characteristics. This includes: aligning the spatial coordinates of the instantaneous density response characteristics and the ant path branching probability characteristics; calling the density-gradient field coupling operator to perform convolution coupling on the spatially aligned features to generate a density-path coupling feature set that integrates the dielectric field and resistivity field information. The specific process can be as follows:

[0123] The nearest neighbor interpolation method is used to unify the instantaneous density response characteristics and ant path branching probability characteristics of different spatial resolutions into the same grid (e.g., 0.25m grid), ensuring that each grid point has both density and probability information.

[0124] Design density-gradient field coupling operator Ψ, Ψ(ρ, grad E) = [ρ(t, ω)⊙P branch (x, y, z)] ○ C(x, y, z, t); where: ρ(t, ω) is the instantaneous density response characteristic. P branch (x, y, z) is the ant path branch probability feature. ⊙ represents the element-by-element product. ○ represents the convolution operation. C(x, y, z, t) is the spatiotemporal correlation correction function, which describes the strength of the feature association in space and time.

[0125] Preferably, it is of the form C(x, y, z, t) = exp(-ds / L spatial )×exp(-dt / T temporal ),in,

[0126] ds and dt are spatial distance and time interval respectively, L spatial and T temporal (e.g. 2m and 300s) are the spatial and temporal correlation lengths determined by historical data statistics.

[0127] This step does not simply juxtapose the two features, but rather mathematically establishes a deep correlation between them through a physically meaningful coupling operator. This allows the subsequent objective function to be optimized based on this fused feature, rather than processing two isolated information sources.

[0128] Step S420 , constructing and solving the objective function of density-path mutual feedback, which integrates the dielectric inversion consistency, density gradient-resistivity curl coupling, and path topology rationality, and is solved through iterative optimization.

[0129] Dielectric inversion consistency constraints are used to ensure that the final density estimation results are consistent with the dielectric response analysis results;

[0130] Cross-physics coupling constraints to mathematically relate the density space gradient to the resistivity gradient field curl;

[0131] Topological rationality regularization term, which is used to impose structural rationality penalties on the ant path network to avoid generating illogical paths;

[0132] Time consistency constraints are used to maintain the smooth continuity of the solution in the time series.

[0133] Feedback refers to the iterative solution process in which: 1) the current density distribution estimate serves as spatial prior information to guide and constrain the prediction of the ant path; 2) conversely, the updated ant path topology serves as structural prior information to correct and constrain the density distribution estimate. This forms a closed-loop, interactive, and mutually reinforcing optimization process.

[0134] Specifically, in this application, the following process can be adopted:

[0135] Constructing a joint objective function J: Based on the coupling feature set generated in step S410, construct the following joint objective function: J = || |ρ est -ρ diel ∣∣ 2 +λ||grad ρest -curl(grad E)|| 2 +μΩ(topology)+γR(consistency);

[0136] The above items are dielectric inversion consistency, cross-physics coupling, topological rationality and time consistency;

[0137] Where: est is the instantaneous density distribution of termites to be solved. diel is the instantaneous density response characteristic obtained from the dielectric response analysis. grad ρ est is the spatial gradient of the density distribution to be determined. curl (grad E) is the gradient curl calculated from the resistivity data. λ, μ, and γ are the weight coefficients for each constraint. Ω (topology) is the topological rationality regularization term.

[0138] In this embodiment, its specific form is Ω(topology)=w1Σκ 2 +w2Σ(1 / L seg ), where κ is the curvature of the predicted path, L segis the length of the isolated path segment, which is used to penalize sharp turns that do not conform to biological habits and meaningless short path segments. w1 and w2 are internal weights. R (consistency) is the time consistency constraint. In this embodiment, the specific form is R (consistency) = || ρ est (t)-ρ est (t-1)|| 2 , which is used to ensure the smoothness of the density solution in time.

[0139] Preferably, the L-curve method is used to find the inflection point of the L-curve by calculating the smoothness and fitting error of the solution under different weight combinations. The weight combination corresponding to the inflection point is the optimal weight λ, μ, γ.

[0140] In other words, the L-curve method is used to find the optimal balance point between the smoothness and fitting error of the solution under different weight coefficient combinations; and the weight coefficient combination corresponding to the optimal balance point is used as the optimal weight in the joint objective function.

[0141] Preferably, a modified two-variable conjugate gradient method is used for the solution. In each iteration, the density distribution ρ and the path probability P are alternately optimized and updated. Each update is performed by calculating the current gradient and combining it with the previous search direction weighted by the conjugate coefficient calculated using the Polak-Ribière formula to determine the current optimal search direction. An Armijo line search is then performed along this direction to determine the optimal step size to complete the update. Specifically, the line search is used to determine the optimal step size and update the density distribution or path probability.

[0142] In the iterative process, mutual feedback is achieved through the mask matrix. For example, using the current density distribution ρ k Generate mask M density (The area with density greater than the threshold is 1, and the rest is 0), then the path is updated to P k '=P k ×M density Then use the updated path P k 'Generate topology mask M path (path connectivity area is 1), reverse constraint density update ρ k+1 '=ρ k+1 ×M path .

[0143] For the specific objective function construction and iteration ideas, please refer to step S4 in Example 6.

[0144] In traditional methods, density estimation and path prediction are two independent processes, prone to false density peaks or lagging paths. This embodiment couples these two problems through a mutually fed inversion framework, leveraging path structural information to eliminate unreasonable density estimates and real-time density information to guide path generation. This achieves information complementarity and significantly improves the accuracy and robustness of the solution.

[0145] Embodiment 4: This embodiment provides a detailed description of a part of step S300 in embodiment 1, namely, the extraction of ant path branch probability features.

[0146] Step S310, the step of extracting the ant trail branching probability feature indicating the ant trail branching tendency includes: calculating a three-dimensional gradient vector field representing the spatial change rate of the physical field based on the resistivity spatial distribution data; performing a curl operation on the three-dimensional gradient vector field to obtain a curl field describing the local rotation intensity of the gradient field; and generating the ant trail branching probability feature by establishing a quantitative correlation model between the intensity of the curl field and the ant trail branching tendency.

[0147] It mainly includes the calculation of three-dimensional gradient field, curl field and branch probability, which can be specifically as follows:

[0148] The standardized resistivity spatial distribution data E(x, y, z) is read, and its partial derivatives in three directions are calculated using the central difference format to obtain the three-dimensional gradient vector field grad E=(dE / dx, dE / dy, dE / dz).

[0149] Based on the gradient vector field grad E, the curl field is calculated by the curl operator grad ×grad E, and the curl vector curl(grad E)=(curl x , curl y , curl z ) and calculate its modulus |curl(grad E)|.

[0150] A model is established to map the curl modulus to the branch probability. In this embodiment, the Sigmoid function is preferably used, which has the form: branch =1 / (1+exp(-(a×∣curl(grad E)∣+b×div factor )) ). Among them: P branch is the final generated ant trail branch probability feature. a, b are the model coefficients obtained by parameter fitting in the calibration field with known ant trail branch positions. factor is an optional scatter factor used to enhance the discrimination of the model.

[0151] The specific calculation process can be found in step S3 of Example 6. The gradient of resistivity indicates how fast it changes, and the curl of the gradient indicates the degree of distortion or rotation of this change. When termites build branching ant tunnels, they often cause drastic and complex disturbances in local geological electrical parameters, which are reflected in the physical field as high-curl areas. This step creatively uses the curl, a topological feature of the physical field, to predict the branching behavior of biological paths, providing key spatial structure prior information for subsequent joint inversion.

[0152] Embodiment 5: This embodiment provides a detailed description of a system-level enhancement function of the entire method.

[0153] Step S500, the method also includes a dynamic update and error correction mechanism, specifically: when new monitoring data is received, quantify the prediction error between the prediction result output by the current model and the actual measured situation reflected by the new data; when the prediction error exceeds the dynamic threshold, start adaptive correction; update the key parameters in the joint objective function, and re-execute the joint inversion solution based on the updated parameters to generate error-corrected termite instantaneous density distribution and ant path branch prediction results.

[0154] In this embodiment, the following process may be used:

[0155] When receiving a new standardized multi-source dataset, calculate the multi-dimensional prediction error. For example, the density error Eρ = ||ρ predicted -ρ measured ||2, path error E P (can be calculated based on new alarm records) and form a comprehensive error index E total =w ρ ×E ρ +w P ×E P .

[0156] E total With a dynamic threshold E threshold In comparison, the threshold can be set based on the mean and standard deviation of historical errors, such as E threshold =E baseline +k×σ historical .

[0157] When E total >E threshold When θ is θ, the adaptive parameter update is started. Preferably, the parameter that contributes most to the error (such as the weight coefficient λ or the calibration coefficient a) is determined through sensitivity analysis and updated using the gradient descent method: new =θ old -η×(dE total / dθ), where η is the update step size. Then, based on the new parameter θ new, re-execute the iterative inversion optimization solution process of Example 3.

[0158] Because earth-rock dams are dynamic systems and termite activity is time-varying, any fixed model is subject to deviation over time. The dynamic update and error correction mechanisms provided in this embodiment enable the entire method to be self-learning and self-correcting. It continuously utilizes newly acquired data to calibrate the model, ensuring the accuracy and reliability of long-term monitoring and achieving a technological leap from static diagnosis to dynamic tracking and early warning.

[0159] Example 6: In order to make the technical solution of the present invention clearer, a simplified numerical calculation example running through the core steps is provided below.

[0160] Assume that in a monitoring area of ​​an earth-rock dam, we focus on a 2 2, the grid point coordinates are (1, 1), (1, 2), (2, 1), (2, 2).

[0161] Step S1: Acquire and process data, including:

[0162] Resistivity data E: The resistivity spatial distribution data (unit: Ω·m) measured at a certain moment is:

[0163] E=(120, 80; 150, 160);

[0164] Dielectric constant data: At the grid point (1, 2), a high-frequency dielectric constant time series data is collected, where when t raw =10s, the real parts of the dielectric constants of the two main frequency bands f1 and f2 are ε'1(10)=15 and ε'2(10)=18, and the phases are φ1(10)=1.2rad and φ2(10)=1.8rad, respectively.

[0165] Geological parameters: It is known that the tan(δ) of this area is tan(δ1)=0.1 and tan(δ2)=0.12 in two frequency bands.

[0166] Alarm data, text data recording the occurrence of termites.

[0167] Specifically, the alarm text data is processed and converted into a binary time series corresponding to a spatial grid. Alarms occurring at specific time points and monitoring point locations are marked as 1, and otherwise as 0. This alarm event sequence can be used as labels for model training or ground truth for validation in subsequent steps, thereby converting unstructured event information into quantitative data suitable for mathematical modeling.

[0168] Step S2: performing phase-preserving parallel frequency-domain time mapping, including:

[0169] Calculate the nonlinear time map: Assume that the compression coefficients α1=0.15, α2=0.12 obtained by calibration, and the regularization parameter ε reg =0.001.

[0170] t sync_1 =10×[1+0.15×log(15+1) / (0.1+0.001)]≈51.1s;

[0171] t sync_2 =10×[1+0.12× / (0.12+0.001)]≈39.5s;

[0172] Perform phase compensation: Calculate the instantaneous phase difference Δφ 12 =1.8-1.2=0.6rad. Assuming the phase difference threshold θ threshold =0.5rad, due to |Δφ 12 |>θ threshold , start compensation. Assuming the compensation intensity coefficient β = 0.3, the phase of frequency band 1 is corrected:

[0173] φ corrected_1= φ1+β×tanh(Δφ 12 / θ threshold )=1.2+0.3×tanh(0.6 / 0.5)≈1.45rad.

[0174] The time-synchronized and phase-corrected data are obtained for subsequent calculations.

[0175] Step S3: Extract instantaneous density and branch probability features

[0176] Extract the ant path branch probability feature P:

[0177] Calculate the gradient field grad E: for 2 The resistivity data E of 2 can be used to estimate the gradient of each point using the central difference method (for simplicity, only the rate of change of the z-direction gradient in the x and y directions of point (1, 1) is given schematically here).

[0178] Calculate the curl field curl(grad E): Calculates the curl based on the gradient field. Assume that the curl modulus |curl(grad E)| at the grid point (1, 2) is 0.8.

[0179] Calculate branch probability: Assume that the calibrated model is P branch =1 / (1+exp(-(2.0×|curl(grad E)|-1.0))).

[0180] P(1,2)=1 / (1+exp(-(2.0×0.8-1.0)))≈0.65.

[0181] Extract the instantaneous density response feature ρ:

[0182] Using the result of step S2, it is assumed that the frequency band response amplitude coefficients A1=0.05, A2=0.04, and the depth weights w1(z)=1, w2(z)=0.9 obtained by multivariate linear regression calibration.

[0183] ρ(t sync ,ω)=A1×ε'1×cos(ω1t sync +φ corrected_1 )×w1+

[0184] A2×ε'2×cos(ω2t sync +φ corrected_2 )×w2;

[0185] Substituting the numerical value, the specific instantaneous density response value can be calculated, for example, ρ≈0.82.

[0186] Step S4: Construct and solve the joint objective function

[0187] Construct the objective function: At this point, we have ρ dielectric ≈0.82, curl(grad E)≈0.8, and ρ to be solved estimated and path topology. The specific form of the joint objective function J is: J= ρ estimated -0.82|| 2 +λ||gradρ estimated -0.8|| 2 +μ·Ω(topology)+γ·R(consistency), where Ω(topology) calculates a penalty value based on the path prediction result of the current iteration, and R(consistency) calculates a penalty value based on the density result of the previous time step.

[0188] Feedback iteration:

[0189] First iteration (density → path): Currently there is a higher density estimate of 0.82 at point (1, 2), so the mask matrix M is generated density , increasing the probability of generating a path at point (1, 2).

[0190] Second iteration (path → density): Generate a topology mask M based on the updated path prediction result that is more likely to pass through point (1, 2) path, the mask will reversely constrain the density estimation, making the density distribution in the path connected area smoother and suppressing isolated false density peaks.

[0191] After multiple iterations until convergence, the converged density distribution map and ant path branch prediction map are finally output.

[0192] If the traditional time window averaging method is used, the instantaneous peak of the dielectric constant in step S2 will be smoothed out, resulting in the instantaneous density ρ calculated in step S3 being likely to be only around 0.3, seriously underestimating the intensity of termite activity. Furthermore, due to the lack of mutual feedback, the path of ant tunnels predicted solely by resistivity anomalies may be inaccurate. Through the above steps, the present invention achieves high-fidelity capture of instantaneous activity and collaborative inversion of multi-physics field information, significantly improving estimation accuracy (e.g., RMSE <15%) compared to traditional methods (e.g., RMSE >40%).

[0193] According to one aspect of the present application, the alarm record is regarded as a sparse ground truth label in the steps of the present application, such as Example 3. Although it cannot provide a complete ant tunnel network, it can confirm the fact that termite damage exists at a specific time and space point.

[0194] After processing the alarm records into binary ground truth label maps (alarm locations are marked as 1 and the rest are 0), its application methods include:

[0195] After the iterative solution is completed, the final predicted ant path branch probability map P can be compared with the ground truth label map to calculate the model's performance indicators such as prediction accuracy and recall rate, thereby conducting an independent and objective quantitative evaluation of the reliability of the inversion results.

[0196] In some optional implementations, the alarm record can also be directly introduced as a strong constraint into the joint objective function J. For example, a penalty term δ·||M can be added. alarm × (1-P)|| 2 Among them, M alarm is a mask matrix that is set to 1 only at alarm locations, and δ is its weight coefficient. The significance of this penalty term is that if the model's predicted path probability P at an alarm location approaches 0, the objective function value will increase sharply, driving the optimization process to converge in a direction that better matches the alarm facts.

[0197] In this way, scattered, non-continuous alarm event information is effectively integrated into the continuous physical field inversion framework. This effectively uses known, deterministic point information to constrain and verify inferences about unknown surface information. This not only provides an objective criterion for evaluating model performance, but also, by imposing hard constraints, forces the model solution to evolve in a direction consistent with known facts. This effectively avoids the unpredictable, erroneous predictions that can arise when the model lacks sufficient constraints, further improving the accuracy and real-world consistency of the overall solution.

[0198] In summary, this invention solves the challenge of synchronizing multi-time-scale monitoring data by constructing a phase-preserving parallel frequency-domain time mapping mechanism, while preserving the transient characteristics of termite activity with high fidelity. This provides highly reliable data input for subsequent precise quantitative estimation. This effect is achieved through the synergistic effect of the following technical features: First, by decomposing the dielectric signal into multiple frequency bands and performing soil-layer-adaptive nonlinear time mapping, the alignment of time scales is ensured to accommodate the heterogeneity of the earth-rockfill dam medium. Second, and most critically, by monitoring the instantaneous phase difference between each frequency band in real time and correcting it using a dynamic compensation algorithm based on the hyperbolic tangent function, the destruction of signal phase information by traditional interpolation or averaging methods is fundamentally avoided. In the context of termite detection in earth-rockfill dams, the phase and amplitude information of the weak and short-lived pulse signals generated by the instantaneous biting activity of termites directly serves as the basis for determining the intensity of the activity. The mechanism of the present invention can ensure that these fingerprint-level characteristic signals remain clearly distinguishable after synchronization across time scales of several orders of magnitude, avoiding the underestimate of activity intensity and enabling the estimation results of instantaneous density to truly reflect the intensity of termite activity. The temporal resolution and accuracy are improved by several orders of magnitude compared to traditional methods.

[0199] Furthermore, the present invention establishes an inversion framework with density-path mutual feedback, coupling the previously independent processes of density estimation and path inference into a collaboratively optimized whole. This improves the physical consistency and spatial rationality of the final solution and significantly reduces the error rate. First, a joint objective function is constructed that incorporates multiple constraints, including dielectric inversion consistency, cross-physics coupling constraints, and topological rationality, to ensure that the solution process is guided by strict physical rules and prior knowledge. Second, a mutual feedback mechanism is introduced during the iterative solution process. Specifically, in the earth-rock dam scenario, high-density activity areas initially estimated from dielectric data provide a high-priority search space for ant path prediction through methods such as masking matrices, avoiding blind path generation. Conversely, connected paths with reasonable topological structures predicted from resistivity data also constrain the density distribution in reverse, effectively eliminating spurious density peaks located in physically unreasonable locations (such as isolated points). This closed-loop feedback optimization mechanism allows the density and path information to verify and correct each other. The final density distribution map and ant path map are not only highly accurate in their own right, but also perfectly match each other. The reliability of the results far exceeds any single physical field analysis or simple result superposition.

[0200] Finally, by proposing a feature that quantifies the branching probability of ant trails based on the curl of the resistivity gradient field, this approach provides more physically insightful prior information for path prediction, improving the ability to predict the complex morphology of ant trail networks, particularly their development trends. This approach creatively applies a fundamental physical field vector analysis tool (curl operations) to biological behavior prediction. In the specific context of earth-rockfill dams, the resistivity gradient field generated by homogeneous soil or gentle moisture variations is typically irrotational or low-rotation. However, when termites construct branching ant trails or nests, they cause dramatic and highly complex local disturbances in the surrounding soil structure, moisture content, and porosity. This complex three-dimensional disturbance manifests itself in the physical field as a high-rotation phenomenon that causes severe distortion and rotation of the gradient vector. By calculating the curl modulus and quantitatively correlating it with branching probability, this paper successfully identifies high-rotation regions, a calculable physical feature, as a biological indicator of the high probability of ant trail branching. Compared with the traditional method of only using high and low resistivity values ​​to judge anomalies, this is a leap from scalar analysis to vector field topology analysis. It provides the joint inversion model with valuable predictive information about where the path will develop, which is of great significance for achieving the technical goal from diagnosis to early warning.

[0201] The preferred embodiments of the present invention are described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within 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 scope of protection of the present invention.

Claims

1. A method for estimating termite density and inferring termite tracts in earth-rock dams, characterized by: include: Process acquired multi-source monitoring data to generate standardized multi-source data sets, including dielectric constant time series data, resistivity spatial distribution data, alarm records, and geological parameters; Perform phase-preserving parallel frequency-domain time mapping on the dielectric constant time series data to generate time-synchronized data; Using time-synchronized data and resistivity spatial distribution data, we extracted the instantaneous density response characteristics that characterize the intensity of termite activity and the ant trail branching probability characteristics that indicate the tendency of ant trail branching. Based on the instantaneous density response characteristics and the ant trail branching probability characteristics, the objective function of density-path feedback was constructed and solved to obtain the instantaneous termite density distribution and ant trail branching prediction results.

2. The method according to claim 1, characterized in that Performs phase-preserving parallel frequency-domain time mapping, including: Parse the dielectric constant time series data into multiple independent frequency bands to obtain frequency band dielectric data; In the process of parallel time mapping of the frequency-band dielectric data, the instantaneous phase difference between each frequency band is detected and quantified in real time; When the instantaneous phase difference deviates from the preset coupling relationship, the phase of each frequency band is dynamically compensated and corrected to generate time-synchronized multi-band coupling data and phase-related characteristic parameters, which are then used to extract the instantaneous density response characteristics.

3. The method according to claim 1, characterized in that The steps to construct and solve the objective function of density-path interaction include: Based on the instantaneous density response characteristics and the ant path branching probability characteristics, a joint objective function is constructed that integrates the dielectric inversion consistency, the density gradient-resistivity curl coupling, and the path topology rationality. Solve the joint objective function through an iterative optimization algorithm, and each iteration covers the interaction and bidirectional constraints; Based on the current density distribution estimation results, spatial weights are applied to the ant path branch prediction to guide the generation of paths; Based on the updated ant path branch prediction results, topological constraints are imposed on the density distribution estimation in reverse order to correct the spatial distribution of density. Continue executing until convergence to obtain the instantaneous termite density distribution and ant path branch prediction results.

4. The method according to claim 1, wherein Based on the instantaneous density response characteristics and the ant path branching probability characteristics, the objective function of density-path mutual feedback is constructed, including: The instantaneous density response characteristics and ant path branching probability characteristics are spatially aligned in the same coordinate system to ensure that each grid point carries both density and probability information. The density-gradient field coupling operator is called, and the spatially aligned instantaneous density response characteristics are convolved with the ant path branching probability characteristics through the preset spatiotemporal correlation function. A density-path coupling feature set that integrates the dielectric field and resistivity field information is constructed, and applied to the construction of the joint objective function.

5. The method according to claim 2, characterized in that The process of parallel time mapping of frequency-band dielectric data includes: The dielectric loss tangent value of the soil layer corresponding to each frequency band is extracted from the standardized multi-source data set; An independent nonlinear time mapping function is constructed for each frequency band, where the compression degree of the time mapping is modulated by the extracted dielectric loss tangent value; Dynamically calculate adaptive weights for each frequency band based on soil depth; Adaptive weights are used to perform weighted combination of the independent mapping results of each frequency band to generate a time mapping result that adapts to the heterogeneity of the soil layer and is used for real-time detection of instantaneous phase differences.

6. The method according to claim 1, characterized in that The steps of extracting the ant trail branching probability feature indicating the ant trail branching tendency include: Based on the resistivity spatial distribution data, the three-dimensional gradient vector field representing the spatial change rate of the physical field is calculated; Perform a curl operation on the three-dimensional gradient vector field to obtain a curl field that describes the local rotation strength of the gradient field; By establishing a quantitative correlation model between the intensity of the curl field and the branching tendency of the ant path, the ant path branching probability characteristics are generated.

7. The method according to claim 3, characterized in that The joint objective function includes the following constraints: Dielectric inversion consistency constraints are used to ensure that the final density estimation results are consistent with the dielectric response analysis results; Cross-physics coupling constraints to mathematically relate the density space gradient to the resistivity gradient field curl; Topological rationality regularization term, which is used to impose structural rationality penalties on the ant path network to avoid generating illogical paths; Time consistency constraints are used to maintain the smooth continuity of the solution in the time series.

8. The method according to claim 2, characterized in that Parallel time mapping of the frequency-band dielectric data is achieved through a parallel computing architecture that includes: An independent parallel processing unit set for each frequency band, responsible for performing the time mapping calculation of the corresponding frequency band; A shared memory mechanism for data exchange between multiple parallel processing units; A synchronization barrier is used to ensure that all parallel processing units complete the calculation before entering the next step, thereby eliminating the cumulative delay of serial calculation.

9. The method according to claim 2, characterized in that Dynamic compensation correction is performed on the phase of each frequency band, including: Comparing the instantaneous phase difference detected in real time with a preset phase difference threshold; When the absolute value of the instantaneous phase difference exceeds the threshold, phase compensation is started; Among them, phase compensation is achieved through a correction algorithm based on the hyperbolic tangent function.

10. The method according to claim 7, characterized in that The weight coefficients of each constraint in the joint objective function are determined by an adaptive method, specifically: The L-curve method is used to find and determine the optimal balance between the smoothness and fitting error of the solution under different weight coefficient combinations; The weight coefficients corresponding to the optimal balance point are combined as the optimal weight in the joint objective function.

Images