A soil and groundwater collaborative monitoring and hidden trouble elimination method and system

Abstract

The application discloses a soil and groundwater cooperative monitoring and hidden trouble investigation method and system, and belongs to the field of environmental monitoring and pollution prevention and control. The method comprises the following steps: acquiring a first data set based on a pre-constructed three-dimensional space monitoring network; processing the first data set by executing a deterministic multiphase flow coupling and a deep gradient evaluation logic to generate a second data set, wherein the second data set comprises a three-dimensional heterogeneous space-time grid and a dynamic sampling control instruction; extracting hydrogeological parameters from the first data set and executing a reverse hydrogeochemical convection and diffusion tracking logic based on the hydrogeological parameters to generate three-dimensional coordinates of a leakage source; and outputting the three-dimensional coordinates of the leakage source and the second data set. The application takes into account the high time resolution advantage of sensor data and the high precision advantage of test data.

Description

Technical Field

[0001] This invention relates to the field of environmental monitoring and pollution prevention, and in particular to a method and system for the coordinated monitoring and hazard investigation of soil and groundwater. Background Technology

[0002] Soil and groundwater pollution is characterized by its high degree of concealment, complex migration pathways, long duration of pollution, and high remediation costs. Pollutants typically do not exist in isolation within a single medium; rather, they may enter the surface soil layer from a surface pollution source, migrate downwards through the vadose zone, and ultimately enter the groundwater aquifer, spreading laterally along the groundwater flow direction. Therefore, in actual pollution hazard investigations, monitoring only surface soil or groundwater is often insufficient to promptly identify the entire process of pollutant migration from shallow to deep layers, and it is also difficult to accurately determine the starting location and migration path of the pollution source.

[0003] In existing technologies, soil monitoring and groundwater monitoring are typically implemented separately. Soil monitoring relies heavily on surface sampling, borehole sampling, and laboratory analysis; groundwater monitoring relies heavily on well water sampling, water level observation, and analysis of groundwater chemical indicators. While these methods can obtain high-precision data for local locations, they lack spatial continuity and have low temporal resolution, making it difficult to form a unified three-dimensional monitoring data structure covering the surface soil layer, vadose zone, and groundwater aquifer. When pollutants are in the transitional stage of migration from the surface to groundwater, traditional separate monitoring methods often fail to accurately capture the precursory signals of migration, such as shallow anomalies followed by deep responses. Summary of the Invention

[0004] One of the objectives of this invention is to provide a method for the coordinated monitoring and hazard investigation of soil and groundwater, in order to solve the problems of cross-boundary interpolation errors caused by complex geological conditions and insufficient positioning accuracy due to the lack of deterministic reverse tracing methods based on physical constraints in pollution source tracing in the existing technology.

[0005] This invention is achieved through the following technical solution: a method for collaborative monitoring and hazard investigation of soil and groundwater, comprising the following steps: Based on a pre-constructed three-dimensional spatial monitoring network covering the surface soil layer, vadose zone, and groundwater aquifer, a first dataset is acquired. The first dataset includes in-situ sensing sequence data, non-contact geophysical exploration boundary data, and borehole sampling final value data. The first dataset is processed by executing deterministic multiphase flow coupling and depth gradient evaluation logic to generate a second dataset. The second dataset includes a three-dimensional heterogeneous spatiotemporal grid and dynamic sampling control instructions. When the bottom layer data corresponding to the groundwater aquifer in the three-dimensional heterogeneous spatiotemporal grid meets the deep anomaly threshold, hydrogeological parameters are extracted from the first dataset, and based on the hydrogeological parameters, reverse hydrogeochemical convection-diffusion tracing logic is executed to generate the three-dimensional coordinates of the leakage source. The three-dimensional coordinates of the leakage source and the second dataset are output.

[0006] Furthermore, when the dynamic sampling control command increases the sampling frequency of the deep node to exceed a preset high-frequency threshold, the corresponding edge computing node executes a data flow control mechanism based on differential incremental hashing, including: dividing the continuously sampled data stream into data blocks of fixed byte size; calculating hash digest values ​​for data blocks in adjacent sampling periods, performing a bitwise XOR operation on the two hash digest values, and calculating the Hamming distance to obtain a data change metric between adjacent data blocks; when the data change metric exceeds a preset data change tolerance threshold, writing the data block into a transmission buffer pool and uploading it to the cloud; when the data change metric does not exceed the data change tolerance threshold, overwriting the data block locally.

[0007] Furthermore, the hydrogeological parameters include groundwater flow direction vector, groundwater flow velocity scalar, soil porosity scalar, and permeability coefficient scalar.

[0008] Furthermore, the three-dimensional heterogeneous spatiotemporal grid is generated through the following steps: using the non-contact geophysical exploration boundary data as three-dimensional spatial geometric constraints, the three-dimensional spatial monitoring network is divided into multiple three-dimensional voxel units; based on the hydraulic distance interpolation algorithm of geological permeation topology, the in-situ sensing sequence data is mapped to the corresponding three-dimensional voxel units to generate an initial voxel state vector; the borehole sampling final value data is extracted as a hard calibration factor, and a weighted residual correction calculation is performed on the initial voxel state vector to generate the three-dimensional heterogeneous spatiotemporal grid.

[0009] Furthermore, the hydraulic distance interpolation algorithm based on geological permeability topology includes: for each three-dimensional voxel unit, calculating the hydraulic distance from each in-situ sensor to the three-dimensional voxel unit, wherein the hydraulic distance is defined as a scalar value obtained by taking the reciprocal of the permeability at each point along the spatial path from the in-situ sensor position to the target three-dimensional voxel unit and then performing line integration, wherein the permeability is obtained by inversion processing of the non-contact geophysical exploration boundary data; using the negative power of the hydraulic distance as the contribution weight of each in-situ sensor to the three-dimensional voxel unit, and performing a weighted average of the actual measurement values ​​of each in-situ sensor at the current moment to obtain the initial voxel state vector of the three-dimensional voxel unit.

[0010] Further, the weighted residual correction calculation includes: calculating the deviation between the final value data of the borehole sampling at the time of laboratory sampling and the interpolated state value of the corresponding three-dimensional voxel unit at the same time, as the calibration residual; weighting the calibration residual spatially using a spatial reliability tensor matrix, wherein the spatial reliability tensor matrix independently adjusts the calibration intensity in the horizontal and vertical directions; attenuating the spatially weighted calibration residual temporally using an exponential time decay factor, wherein the exponential time decay factor approaches zero as the time interval between the current time and the laboratory sampling time increases; and superimposing the spatially weighted and time-decayed calibration residual onto the interpolated state value at the current time to obtain the calibrated environmental state value and writing it into the three-dimensional heterogeneous spatiotemporal grid.

[0011] Furthermore, the time decay constant in the exponential time decay factor is adjusted according to the severity of changes in the monitoring site environment: in sites where environmental conditions change rapidly, the time decay constant is taken to make the calibration effectiveness decay rapidly; in sites where environmental conditions are relatively stable, the time decay constant is taken to make the calibration effectiveness last longer.

[0012] Furthermore, the deterministic multiphase flow coupling includes: establishing a set of multiphase flow equations on the three-dimensional heterogeneous spatiotemporal grid and performing spatiotemporal evolution dynamics calculations, wherein: for the vadose zone region, the Richards equation is used to describe the unsaturated water flow flux; for the groundwater aquifer region, Darcy's law is used to calculate the saturated water flux; at the physical interface between the vadose zone and the groundwater aquifer, a mass conservation boundary condition is set, which forces the normal mass flux on both sides of the interface to be strictly equal.

[0013] Furthermore, the mass conservation boundary condition specifically states that: on the upper side of the physical interface, the unsaturated vertical flux calculated by the Richards equation includes the combined contribution of capillary driving force and gravitational components; on the lower side of the physical interface, the saturated vertical flux calculated by Darcy's law is driven by the head gradient; the mass conservation boundary condition forces the unsaturated vertical flux on the upper side to be equal to the saturated vertical flux on the lower side, and uses the unsaturated mass flux output by the Richards equation as the inflow flux input to the Darcy's law calculation grid, thereby achieving seamless coupling of pollutants across phase interfaces.

[0014] Furthermore, the depth gradient assessment logic includes: extracting a first depth data vector corresponding to the depth range of the surface soil layer and a second depth data vector corresponding to the depth range from the lower part of the vadose zone to the groundwater aquifer from the vertical direction of the three-dimensional heterogeneous spatiotemporal grid; calculating the rate of change of the first time derivative of the concentration of a specific pollutant in the first depth data vector; comparing the rate of change of the first time derivative with a dynamic correction threshold; and based on the comparison result and the current state of the deep monitoring nodes in the second depth data vector, executing a cascade early warning judgment and generating the dynamic sampling control command.

[0015] Furthermore, the dynamic correction threshold is determined as follows: Meteorological precipitation sequence data is acquired synchronously, and based on the baseline warning threshold, an equivalent concentration dilution bias caused by rainfall infiltration is superimposed to obtain the dynamic correction threshold. The equivalent concentration dilution bias is calculated as follows: the rainfall amount with time delay is multiplied by the soil infiltration coefficient to obtain the actual infiltrated water volume, then divided by the shallow soil porosity to obtain the concentration dilution amplitude caused by rainwater injection, and then multiplied by the rainfall dilution response coefficient for scaling. The time delay reflects the time required for rainfall to infiltrate from the surface to the shallow sensor installation depth, and its value is determined by the infiltration path length and surface soil permeability.

[0016] Furthermore, the tiered early warning determination adopts a shallow-deep joint criterion: when the absolute value of the rate of change of the first-order time derivative exceeds the dynamic correction threshold, and the current state value of the deep monitoring node in the second depth data vector is lower than the preset deep safety threshold, the first-level tiered early warning is determined to be triggered; when both shallow and deep anomalies occur simultaneously, the first-level tiered early warning is not triggered.

[0017] Furthermore, after the first-level tiered early warning is triggered, the sampling frequency of each deep node in the dynamic sampling control command is adjusted between the minimum reference sampling frequency and the maximum allowable sampling frequency through the product of the spatial distance factor and the dynamic change factor. Specifically: the spatial distance factor uses a Gaussian kernel distance decay function, with the square of the three-dimensional Euclidean distance between the shallow anomaly point that triggered the early warning and each deep node as the independent variable. Deep nodes that are closer in distance receive a greater frequency boost, and the spatial decay scale parameter of the Gaussian kernel controls the effective radiation radius of the frequency boost effect. The dynamic change factor uses a hyperbolic tangent function, with the product of the current concentration gradient field magnitude and the concentration gradient response coefficient as the independent variable. This factor approaches zero when the concentration gradient is zero and approaches one when the concentration gradient increases.

[0018] Furthermore, the execution of the reverse hydrogeochemical convection-diffusion tracing logic includes: using hydrogeological parameters as tracing boundary parameters, and using three-dimensional voxel units in the three-dimensional heterogeneous spatiotemporal grid that satisfy the deep anomaly threshold as initial source and sink terms to construct a reverse convection-diffusion equation; solving the reverse convection-diffusion equation on the three-dimensional heterogeneous spatiotemporal grid to obtain a spatial probability density field; and extracting the peak cluster center coordinates of the spatial probability density field to determine the three-dimensional coordinates of the leakage source.

[0019] Furthermore, the reverse convection-diffusion equation is constructed by reversing the time variable into a reverse time variable, which increases in the past direction from the current moment when the anomaly is triggered. The reverse convection-diffusion equation includes a dispersion term, a convection term, and a source-sink term, wherein: the dispersion term describes the diffusion and distribution effect of pollutants caused by medium inhomogeneity and molecular diffusion using a hydrodynamic dispersion coefficient tensor; the convection term describes the propagation of the spatial probability density field in the opposite direction of the groundwater flow field using an effective reverse velocity vector; the source-sink term includes a spatial pruning operator used to exclude physically impermeable regions from the solution domain; the initial condition of the spatial probability density field is to take a high probability density value at the three-dimensional voxel unit corresponding to the initial source-sink term, and zero or near-zero values ​​at other locations.

[0020] Furthermore, the effective reverse velocity vector is corrected using a geochemical blocking factor, including: based on real-time pH sensing data in the first dataset, using a quadratic polynomial fitting of the nonlinear relationship between the soil adsorption partition coefficient and pH value, combined with soil dry bulk density and porosity, to calculate the pH-dependent geochemical blocking factor; dividing the Darcy velocity vector by the product of porosity and the geochemical blocking factor to obtain the effective reverse velocity vector; wherein, the fitting coefficient of the quadratic polynomial is calibrated using laboratory batch adsorption experimental data.

[0021] Furthermore, the spatial pruning operator is implemented in the following way: reading the non-contact geophysical exploration boundary data from the three-dimensional heterogeneous spatiotemporal grid, identifying three-dimensional voxel units with permeability lower than a preset minimum permeability threshold as absolute aquitard voxels; in each iteration cycle of solving the reverse convection diffusion equation, applying a maximum negative penalty term to the source-sink term of the absolute aquitard voxel, so that the probability density of the absolute aquitard voxel decays to zero during the solution process.

[0022] Furthermore, the reverse convection diffusion equation is solved using the finite volume method. By integrating the reverse convection diffusion equation on each three-dimensional voxel unit, a discrete set of equations is established. The solution is then gradually advanced with a preset reverse time step to ensure that the mass flux balance of each three-dimensional voxel unit is maintained.

[0023] Furthermore, extracting the peak cluster center coordinates of the spatial probability density field includes: when the spatial probability density field has a unimodal distribution, performing a weighted spatial integral on the spatial probability density field in the solution domain to calculate its centroid position, which is used as the three-dimensional coordinates of the leakage source; when the spatial probability density field has a multimodal distribution, performing density-based clustering analysis on the spatial probability density field to extract the cluster center coordinates of each peak as the three-dimensional coordinates of multiple candidate leakage sources.

[0024] Another aspect of the present invention provides a soil and groundwater co-monitoring and hazard investigation system, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements any of the soil and groundwater co-monitoring and hazard investigation methods described above.

[0025] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0026] 1. This invention uses hydraulic distance based on geological permeability topology to replace traditional Euclidean distance for spatial interpolation. The line integral of the inverse of permeability along the path is used as the basis for weight calculation. Geological boundary information obtained from geophysical exploration is used as hard constraints for interpolation. When there is a low-permeability aquitard in the path, the contribution weight of the corresponding sensor automatically tends to zero. This eliminates the cross-boundary interpolation error caused by ignoring geological barriers under complex strata conditions and improves the reconstruction accuracy of the three-dimensional state field.

[0027] 2. This invention designs a weighted residual correction mechanism with an exponential time decay factor and a spatial reliability tensor. It uses borehole test data as a high-precision benchmark to correct the systematic bias of sensor interpolation results. At the same time, it uses time decay to reasonably reduce the calibration effectiveness over time, taking into account both the high temporal resolution advantage of sensor data and the high precision advantage of test data. Furthermore, by establishing mass conservation forced boundary conditions for the unsaturated and saturated zone control equations at the groundwater surface, it achieves a deterministic physical characterization of the cross-phase migration process of pollutants. Simultaneously, it introduces a dynamic correction threshold based on rainfall infiltration compensation and a shallow-deep dual-condition joint criterion, which effectively eliminates false alarms caused by meteorological interference and improves the accuracy and reliability of cascade early warning.

[0028] 3. This invention constructs a reverse convection-diffusion equation and introduces a pH-dependent geochemical blocking factor to dynamically correct the reverse flow velocity. Combined with a spatial pruning mechanism, it eliminates false diffusion in impermeable areas, achieving precise physical reverse tracing from deep anomaly detection points to the surface leakage source, thus improving the accuracy and reliability of source tracing and location. At the same time, it designs an adaptive sampling frequency adjustment mechanism based on spatial distance attenuation and concentration gradient response, as well as a data flow control mechanism based on differential incremental hashing. This effectively controls the resource consumption of edge nodes while ensuring that abnormal peak data is not lost, achieving a balance between monitoring accuracy and computing resources. Attached Figure Description

[0029] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

[0030] Figure 1 This is a flowchart of the method provided in Embodiment 1 of the present invention.

[0031] Figure 2 This is a schematic diagram of the spatial layering of the three-dimensional collaborative monitoring network provided in Embodiment 1 of the present invention.

[0032] Figure 3 This is a spatial error comparison diagram provided for Embodiment 1 of the present invention.

[0033] Figure 4 The graph shows the attenuation curve of the distance attenuation power parameter as a function of hydraulic resistance, as provided in Embodiment 1 of the present invention.

[0034] Figure 5 This is a voxel pollutant concentration-time curve provided in Embodiment 1 of the present invention.

[0035] Figure 6 The attenuation effect curve provided in Embodiment 1 of the present invention.

[0036] Figure 7This is a schematic diagram illustrating the verification of mass flux continuity at the coupling interface provided in Embodiment 1 of the present invention.

[0037] Figure 8 This is a comparison chart of false alarms in a rainfall event provided in Embodiment 1 of the present invention.

[0038] Figure 9 This is a comparison chart of edge node cache usage provided in Embodiment 1 of the present invention.

[0039] Figure 10 This is the differential incremental hash and dynamic threshold verification upload determination diagram provided in Embodiment 1 of the present invention.

[0040] Figure 11 This is a comparison chart of the convergence speed of the inverse inversion residuals before and after spatial pruning, as provided in Embodiment 1 of the present invention.

[0041] Figure 12 This is a two-dimensional cross-sectional view of the spatial probability distribution and leakage source location provided in Embodiment 1 of the present invention. Detailed Implementation

[0042] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0043] Example 1

[0044] This embodiment discloses a method for collaborative monitoring and hazard mitigation of soil and groundwater. The basic implementation involves constructing a deterministic mathematical model system with physical equation constraints at its core, establishing a three-dimensional spatial monitoring network covering the surface soil layer, vadose zone, and groundwater aquifer. In-situ sensor sequence data, non-contact geophysical exploration boundary data, and borehole sampling final value data are acquired through multimodal environmental data acquisition as unified inputs. The core data processing flow is divided into three progressive modules: The first module, based on hydraulic distance interpolation constrained by geological permeability topology and combined with the exponential decay calibration mechanism of laboratory data, fuses heterogeneous multi-source data into a physically self-consistent three-dimensional heterogeneous spatiotemporal grid. On this basis, the Richards equation and Darcy's law are used to establish a mass-conserving coupling boundary between the unsaturated and saturated zones, deterministically characterizing the cross-phase migration trajectory of pollutants. The second module detects abnormal signals through a tiered early warning mechanism based on a shallow-deep joint criterion. It uses dynamic correction thresholds to eliminate false alarms caused by external hydrological turbulence such as rainfall, and dynamically adjusts the sampling frequency of deep nodes through a dual adaptive mechanism based on Gaussian kernel distance decay and concentration gradient response. Simultaneously, it employs a differential incremental hash flow control mechanism to prevent memory overflow at edge nodes due to high-frequency sampling. The third module uses the reverse convection-diffusion partial differential equation as the core of source tracing. It dynamically corrects the geochemical hindrance factor using real-time pH sensor data to adjust the reverse flow velocity, and combines this with a space pruning machine to prune impermeable areas to accelerate convergence and improve positioning accuracy. Finally, it achieves deterministic output of the three-dimensional coordinates of the leakage source through peak clustering of the reverse probability density field. Figure 1 The overall method flowchart in this embodiment is shown, combined with Figure 1 As can be seen, this embodiment includes the following steps:

[0045] Step 1: Based on the pre-constructed three-dimensional spatial monitoring network covering the surface soil layer, vadose zone and groundwater aquifer, obtain a multimodal environmental dataset as the first dataset.

[0046] Among them, the three-dimensional spatial monitoring network refers to a three-dimensional monitoring network composed of multi-layered and multi-type sensing and detection devices deployed in the horizontal and vertical directions according to a preset spatial density within the target monitoring area. This three-dimensional spatial monitoring network covers the complete stratigraphic profile from the surface to the depths in the vertical direction, specifically including the surface soil layer, the vadose zone (i.e., the unsaturated zone), and the groundwater aquifer (i.e., the saturated zone).

[0047] The surface soil layer refers to the shallow soil medium located below the surface and down to the bottom of the plant root zone. It is usually the first contact layer for pollutant infiltration, and its thickness varies from tens of centimeters to one to two meters depending on site conditions.

[0048] The vadose zone refers to the unsaturated stratum region located between the bottom of the surface soil layer and the groundwater level. Within this region, both gas and liquid phases exist simultaneously in the soil pores, and water movement is controlled by both capillary action and gravity. The vadose zone is the essential pathway for the vertical migration of pollutants from the surface to the groundwater aquifer, and its seepage characteristics are crucial for determining the infiltration rate and time it takes for pollutants to reach groundwater.

[0049] A groundwater aquifer is a saturated geological region located below the groundwater level, where the pores of the soil or rock are completely filled with water. Water flow within this region is primarily driven by the hydraulic gradient and follows Darcy's law. Groundwater aquifers are the carriers of groundwater resources and also the main sites for the lateral diffusion and long-distance migration of pollutants once they arrive.

[0050] A multimodal environment dataset refers to a collection of raw observation data covering a variety of physical quantities and chemical parameters, which are acquired synchronously or asynchronously in a three-dimensional spatial monitoring network through different types of monitoring methods. This dataset serves as a unified input source for all subsequent calculations and analyses.

[0051] The first dataset contains in-situ sensing sequence data, non-contact geophysical exploration boundary data, and borehole sampling final value data.

[0052] In-situ sensing sequence data refers to time-series monitoring data continuously collected at fixed or adjustable time intervals by in-situ sensors buried at different depths underground. In-situ sensors include, but are not limited to, soil temperature sensors, soil moisture sensors, conductivity sensors, pH sensors, dissolved oxygen sensors, and sensors for specific pollutant concentrations. In-situ sensing sequence data is characterized by high temporal resolution but limited spatial coverage; that is, each sensor can only represent the environmental conditions within a limited spatial range near its installation location.

[0053] Non-contact geophysical boundary data refers to data reflecting the spatial distribution of physical properties of subsurface media, acquired non-invasively using geophysical exploration devices deployed on or above the Earth's surface. This type of data includes, but is not limited to, apparent resistivity profiles obtained using the high-density resistivity method, dielectric constant distributions obtained using ground-penetrating radar, and conductive strata obtained using transient electromagnetic methods. The core value of non-contact geophysical boundary data lies in providing three-dimensional spatial geometric constraints on subsurface geological structures, including the location of stratigraphic interfaces, the distribution and thickness of aquitards, and the orientation of fault zones. This information constitutes the hard boundary conditions for subsequent three-dimensional spatial interpolation and physical modeling.

[0054] Borehole sampling final value data refers to high-precision chemical analysis results obtained by drilling boreholes and extracting soil or groundwater samples within the monitoring area, followed by laboratory testing and analysis. Borehole sampling final value data includes, but is not limited to, absolute concentrations of specific pollutants, soil organic matter content, and heavy metal content. This type of data is characterized by high accuracy but low sampling frequency and time lag, as it typically takes several hours to several days from sampling to the issuance of a laboratory report. In this method, borehole sampling final value data is used as a high-precision calibration benchmark to correct for systematic drift or cumulative errors that may exist in in-situ sensor data.

[0055] In this embodiment, in-situ sensing sequence data can be transmitted in real time to an edge computing gateway or cloud server via wired or wireless means for aggregation; non-contact geophysical exploration boundary data can be collected periodically or on demand and then imported offline into the system; borehole sampling final value data can be imported into the system manually or through an automated interface after laboratory testing.

[0056] In this embodiment, the deployment density and depth range of the in-situ sensors are determined based on the geological conditions, pollution risk level, and regulatory requirements of the monitoring site to ensure that the three-dimensional spatial monitoring network can effectively cover all potential pollution migration paths. Figure 2 This diagram illustrates the spatial stratification of the three-dimensional collaborative monitoring network of soil-vadose zone-groundwater in this embodiment. Figure 2 The middle layer shows the surface soil layer, vadose zone, and groundwater aquifer, and marks the spatial distribution of shallow sensors, vadose zone sensors, groundwater sensors, and borehole sampling points.

[0057] Step 2: Process the first dataset by executing the deterministic multiphase flow coupling and deep gradient evaluation logic to generate a second dataset containing a three-dimensional heterogeneous spatiotemporal grid and dynamic sampling control instructions.

[0058] Deterministic multiphase flow coupling refers to the physical calculation process that describes the movement of pollutants in the unsaturated vadose zone and saturated aquifers and their cross-phase coupling relationships in the form of deterministic partial differential equations, based on the classical hydrogeological equation system (including Richards' equations and Darcy's law). This process differs from empirical models based on statistical correlations, as each step of the calculation has a clear physical meaning and is interpretable.

[0059] The depth gradient assessment logic refers to a set of deterministic judgment rules that perform gradient analysis and joint judgment on shallow and deep environmental state parameters along the vertical direction of the monitoring grid. It is used to capture early signals of pollutants migrating from shallow to deep layers in the time dimension and trigger corresponding early warning and sampling scheduling actions accordingly.

[0060] The second dataset refers to the comprehensive output data set generated after the above-mentioned deterministic multiphase flow coupling calculation and deep gradient evaluation logic processing. This dataset contains two core components: a three-dimensional heterogeneous spatiotemporal grid and dynamic sampling control instructions.

[0061] A three-dimensional heterogeneous spatiotemporal grid refers to a discretized spatiotemporal data structure within the geometric framework of a three-dimensional spatial monitoring network. This structure discretizes the monitoring area into multiple three-dimensional voxel units, assigning environmental state parameter values, calculated through multi-source data fusion and physical equation constraints, to each voxel unit. This grid is called "heterogeneous" because the data sources constituting the grid's state values ​​come from various sensing devices with different accuracy characteristics, temporal resolution, and spatial representativeness. Furthermore, the physical equations controlling material transport also have different mathematical forms in different depth regions (unsaturated and saturated zones).

[0062] Dynamic sampling control commands refer to real-time control commands that the system automatically sends to the edge nodes of the corresponding physical sensors based on the judgment results of the deep gradient evaluation logic, in order to adjust the sensor sampling frequency and data transmission strategy.

[0063] The following provides a detailed explanation of the sub-steps included in this step:

[0064] Sub-step 1: Using the boundary data of non-contact geophysical exploration as the three-dimensional spatial geometric constraint, the pre-constructed three-dimensional spatial monitoring network is divided into multiple three-dimensional voxel units of fixed size.

[0065] Among them, the three-dimensional spatial geometric constraints refer to the hard boundary information reflecting the spatial distribution of underground geological structures obtained by non-contact geophysical exploration. This information is used to constrain the spatial division range and topological relationships of the three-dimensional grid. Specifically, the location of stratigraphic interfaces, the spatial distribution range of aquitards, and the orientation of fault zones identified in the geophysical exploration data determine the spatial arrangement and neighborhood connectivity of the three-dimensional voxel units.

[0066] A three-dimensional voxel unit refers to the smallest three-dimensional spatial unit formed by uniformly discretizing the underground space area covered by a three-dimensional spatial monitoring network along the horizontal and vertical directions according to a preset spatial resolution. Each voxel unit has fixed geometric dimensions (such as length, width, and height), and its position in space is indexed by three-dimensional coordinates. The only certainty.

[0067] Understandably, the selection of voxel size needs to strike a balance between spatial resolution and computational efficiency: the smaller the size, the higher the spatial resolution, but the computational load increases accordingly; the larger the size, the higher the computational efficiency, but local stratigraphic details may be missed.

[0068] In this embodiment, the vertical dimension of the three-dimensional voxel unit can be divided into non-equidistant parts according to the fineness of the stratigraphic structure. Smaller vertical dimensions are used in areas with drastic stratigraphic changes, while larger vertical dimensions are used in areas with uniform stratigraphic structures.

[0069] Sub-step 2: Based on the hydraulic distance interpolation algorithm of geological permeability topology, the in-situ sensing sequence data is mapped to the corresponding three-dimensional voxel unit to generate the initial voxel state vector.

[0070] The initial voxel state vector refers to a set of environmental state parameter values ​​assigned to each three-dimensional voxel unit after spatial interpolation. This state vector includes at least physical quantities such as temperature, humidity, electrical conductivity, and concentration of specific pollutants.

[0071] In this embodiment, the initial voxel state vector may also include parameters such as dissolved oxygen concentration, redox potential, and specific heavy metal ion concentration, the specific composition of which is determined by the monitoring target and sensor configuration.

[0072] To resolve the contradiction between the spatial sparsity of in-situ sensors and the spatial continuity of three-dimensional voxel units, this step employs a hydraulic distance interpolation algorithm based on geological permeability topology to replace the traditional inverse Euclidean distance attenuation interpolation. In the traditional inverse distance weighted interpolation method, the contribution weight of each sensor to the target voxel is directly determined by a negative power of the Euclidean distance (i.e., the geometric straight-line distance) between them, which is reasonable under ideal conditions of uniform geological conditions. However, in actual underground environments, there are often geological structures such as impermeable rock layers, dense clay layers, or fault zones between the sensors and the target voxels. The Euclidean distance cannot truly reflect the actual path and resistance encountered by fluids or contaminants migrating from the sensor location to the target voxel.

[0073] Therefore, this step abandons Euclidean distance and instead adopts a hydraulic distance based on geological permeability topology as the basis for weight calculation. The core definition of this hydraulic distance is as follows: along the spatial path from the sensor location to the target voxel, the inverse of the permeability at each point on the path is taken as the line integral, thus obtaining a scalar value characterizing the total resistance encountered by the fluid through the path. For areas with lower permeability, the equivalent hydraulic distance per unit length of the path is larger; when there is a nearly impermeable aquitard in the path, the integral value of that segment of the path tends to be extremely large, making the sensor's contribution weight to the voxel on the other side of the path approach zero. In this way, the unreasonable influence of sensors isolated by geological barriers on the estimation of the target voxel state can be naturally eliminated from the physical mechanism, without the need for manual pre-marking of barrier locations or manual removal of sensors.

[0074] Exemplarily, in this embodiment, three-dimensional voxels At any moment The initial environmental state can be calculated using the following formula:

[0075]

[0076] The weighting function is defined as follows:

[0077]

[0078] in, For target voxels At any moment Interpolated state values, such as pollutant concentration or conductivity; The total number of physical sensors involved in the interpolation calculation; For the first A physical sensor at time The actual measured value; For the first The contribution weight of each sensor to the target voxel; To from the sensor To the target voxel Spatial integration path; For each point on the path The relative permeability at a given point, the permeability field of which is extracted from boundary data obtained by non-contact geophysical exploration after inversion processing; This is a distance decay power parameter used to control the rate at which the weight decreases as the hydraulic distance increases. Figure 3 This diagram shows a comparison of the spatial errors between the hydraulic distance and the traditional Euclidean distance interpolation in this embodiment. Figure 3 Using the actual concentration distribution as a baseline, it can be seen that near the impermeable layer, the Euclidean distance interpolation exhibits severe over-limit penetration errors (the curve deviates from the true value). However, in this embodiment, the hydraulic distance interpolation naturally falls back to a reasonable range due to the weighting suppression of the inverse integral of permeability, demonstrating the hydraulic distance weighting function based on geological permeability topology. Under complex geological conditions, it can effectively eliminate interpolation out-of-bounds errors caused by ignoring the aquitard layer, and the accuracy is significantly improved compared with traditional Euclidean distance interpolation. Figure 4 The graphs showing the attenuation curves of different distance attenuation power parameters as a function of hydraulic resistance in this embodiment are illustrated. Figure 4 This paper presents a family of weighted decay curves under different power-law parameters for distance decay, and marks the typical resistance ranges corresponding to high-permeability paths (sand layer) and low-permeability paths (clay / impermeable layer). It can be seen that the power-law parameter... The effect of regulating the rate of weight decay, and the physical rationale for the surge in hydraulic resistance leading to the weight approaching zero when the path passes through a low-permeability region.

[0079] It is understandable that in the above formula, It is the line integral of the inverse of permeability along the spatial path, and its physical meaning is the fluid flow from the sensor... The installation location was moved to the target voxel. The total hydraulic resistance experienced. When the path passes through a high-permeability sandy aquifer, the integral value for this segment is small, and the corresponding hydraulic distance is short; when the path passes through low-permeability dense clay or rock, the integral value for this segment increases significantly, and the corresponding hydraulic distance is greatly lengthened. The integral value is then raised to a negative power. Subsequently, sensors with greater hydraulic resistance have smaller weights, thus enabling spatial state estimation based on geological permeability topology (rather than simple geometric distance). This approach uses geological boundary information obtained from geophysical exploration as a hard constraint for interpolation calculations, fundamentally eliminating interpolation out-of-bounds errors caused by ignoring geological barriers under complex strata conditions.

[0080] Sub-step 3: Extract the final value data of the borehole sampling as a hard calibration factor, and perform weighted residual correction calculation on the initial voxel state vector to generate a three-dimensional heterogeneous spatiotemporal mesh.

[0081] Among them, the hard calibration factor refers to the data element that directly uses the high-precision measured values ​​obtained from borehole sampling and testing as the benchmark true value of the environmental state in the calibration calculation. It is called "hard" because in the calibration calculation, the accuracy level of the test data is considered higher than that of the sensor interpolation result. When a deviation occurs between the two, the test data is used as the standard for correction.

[0082] Weighted residual correction calculation refers to the process of extracting the deviation between the true value of the test and the interpolation result of the sensor at the same time, using this deviation as a correction amount, and then superimposing it on the interpolation result at the current time after time decay weighting, so as to obtain the final calibrated environmental state value.

[0083] After obtaining the aforementioned spatial interpolation results based on hydraulic distance, further calibration and correction using high-precision borehole sampling and analysis data is required. The motivation is as follows: while in-situ sensor data has the advantage of high temporal resolution, providing continuous monitoring data at intervals of seconds or minutes, sensor elements may develop systematic errors such as zero-point drift or sensitivity decay over time. On the other hand, while borehole sampling and analysis data has extremely high precision (typically comparable to laboratory instrument analysis), its sampling frequency is low (usually once every few days to weeks) and there is a time lag between sampling and report generation. Therefore, a mechanism is needed to combine the advantages of these two data sources: utilizing the high precision of the analysis data to correct for systematic biases in the sensor data, while simultaneously using the high temporal resolution of the sensor data to fill the temporal gaps in the analysis data.

[0084] The core idea of ​​the calibration mechanism designed in this invention is to use laboratory data as the baseline truth value, calculate the deviation (i.e., residual) between the laboratory value at the sampling time and the sensor interpolation value, and then add this residual as a correction amount to the sensor interpolation result at the current time. Simultaneously, considering the time-sensitivity of laboratory data—that is, the further away from the sampling time, the weaker the representativeness of the laboratory data to the current state due to changes in environmental conditions—an exponential time decay factor is introduced to control the decrease of the calibration amount over time, so that the calibration effectiveness gradually declines after the sampling time.

[0085] For example, in this embodiment, the calibrated final environmental state value can be calculated using the following formula:

[0086]

[0087]

[0088] in, voxels At any moment The final environmental state value after calibration is the actual state parameter written into the three-dimensional heterogeneous spatiotemporal grid. The current state value is obtained by the aforementioned hydraulic distance interpolation algorithm. To be at the time of laboratory sampling High-precision measured values ​​obtained through borehole sampling and analysis; To be at the time of laboratory sampling Spatial interpolation; This is a spatial confidence tensor matrix, used to independently adjust the calibration intensity in different spatial dimensions; This is the Hadamard product operator, representing element-wise multiplication. It is a time decay constant that controls the rate at which the calibration effectiveness decays over time; It is the exponential time decay factor. Figure 5 The following is a graph showing the voxel contaminant concentration over time before and after calibration using the calibration mechanism in this embodiment. Figure 5 The study compared the actual trend, the in-situ sensor interpolation curve, the borehole test points, and the calibration curve after residual correction. It can be seen that this scheme corrects the residual of the in-situ sensor interpolation results by using the final value data of borehole sampling, which can reduce sensor drift and interpolation deviation while maintaining real-time performance. Figure 6 The following diagram illustrates the decay effect curves for different time decay constants in this embodiment. Figure 6 This demonstrates the proof of exponential time decay. It can reasonably control the timeliness of the calibration effect of test data, so that the weight of recent test data is greater and the weight of long-term test data naturally declines, thus avoiding unreasonable interference of outdated calibration data on the current state.

[0089] It is understandable that in the above formula, This represents the deviation between the true value of the test and the sensor interpolation result at the time of sample collection. This deviation reflects the systematic error of the sensor data relative to the reference true value at that time. It is an exponential decay factor, when the current time Exactly equal to the time of laboratory sampling At that time, the factor takes a value of 1, indicating that the test calibration is completely effective and the correction amount is fully added; as the current time... Gradually moving away from the sampling time The value of this factor gradually approaches zero from 1, indicating that the calibration power of the test data gradually becomes invalid. The introduction of the spatial reliability tensor matrix allows the system to set different calibration intensity weights in the horizontal and vertical directions. For example, areas with drastic vertical geological changes can be given greater calibration weights, while areas with relatively uniform geological conditions in the horizontal direction can be given smaller calibration weights. By combining time decay with spatial weighting, the high-precision characteristics of the test data are propagated forward along the time axis, ensuring reasonable calibration decay while maintaining real-time performance, thus solving the problems of insufficient accuracy or timeliness of a single data source.

[0090] In this embodiment, the time decay constant Adjustments can be made based on the severity of environmental changes at the monitoring site: in sites where environmental conditions change rapidly (such as industrially active areas). Taking a larger value causes the calibration effectiveness to decay rapidly, prompting the system to rely more on real-time sensor data; in relatively stable environmental conditions (such as remote natural areas). Choose a smaller value to maintain the calibration effectiveness for a longer period of time.

[0091] Sub-step 4: Establish a set of multiphase flow equations and perform spatiotemporal evolution dynamics calculations on a three-dimensional heterogeneous spatiotemporal grid.

[0092] After completing the aforementioned multi-source data fusion and calibration and generating the initial state of the three-dimensional heterogeneous spatiotemporal grid, in order to achieve a deterministic characterization of the spatiotemporal evolution of pollutants in the underground environment, it is also necessary to establish a set of physical equations describing the motion of fluids and pollutants in the unsaturated vadose zone and saturated aquifer, and solve them on the three-dimensional heterogeneous spatiotemporal grid.

[0093] The core feature of the multiphase flow equation set established in this invention is that it adopts appropriate governing equations for two regions with drastically different physical properties, namely the unsaturated vadose zone and the saturated aquifer, and forces the setting of mass conservation boundary conditions at the physical interface between the two regions, thereby achieving seamless coupling calculation across phase states.

[0094] Specifically, in the vadose zone (i.e., the unsaturated zone), both gas and liquid phases exist simultaneously in the soil pores, and water movement is controlled by both capillary suction and gravity. Therefore, the Richards equation is used to describe the unsaturated water flux in this region. In the groundwater aquifer (saturated zone), the soil or rock pores are completely filled with water, and water flow is mainly driven by the hydraulic gradient. Therefore, Darcy's law is used to calculate the saturated water flux in this region.

[0095] For example, in this embodiment, the equation for water movement in the unsaturated vadose zone (Richards equation) can be expressed as follows:

[0096]

[0097] in, Unsaturated hydraulic conductivity is the water content. The function of is physically defined as the amount of water passing through a unit cross-sectional area under a unit hydraulic gradient driven by a given water content. Its value decreases sharply as the water content decreases. The pressure head represents the energy state of soil water under capillary action, and is negative under unsaturated conditions (i.e., suction). The position is the water head, i.e., the vertical coordinate, representing the gravitational potential energy; The total hydraulic gradient is the direction of the resultant force of capillary driving force and gravitational driving force; Specific water capacity is defined as the change in water content caused by a unit change in pressure head, and it characterizes the water storage and release capacity of soil under different humidity conditions.

[0098] Understandably, the physical essence of the Richards equation is a mass conservation description of water movement under unsaturated conditions: the left side of the equation... This represents the net change in water flux flowing into or out of a spatial element per unit time; the right side of the equation... This represents the rate of change of water content within the spatial element over time. When the inflow rate on the left is greater than the outflow rate, the water content within the element increases; conversely, it decreases.

[0099] For example, in this embodiment, the water flow in the saturated aquifer follows Darcy's law, and its flux can be expressed as follows:

[0100]

[0101] in, The Darcy flux vector for the saturation zone represents the amount of water passing through a unit cross-sectional area per unit time and its direction. Saturated hydraulic conductivity is a constant related to the pore structure and viscosity of the aqueous medium. The total head is equal to the sum of the pressure head and the position head. This represents the head gradient, pointing in the direction of the fastest decrease in head. Figure 7 This diagram illustrates the verification of mass flux continuity at the coupling interface of Richards' equation and Darcy's law in this embodiment. Figure 7 The two curves in the figure represent the flux calculated by Richards' equation in the unsaturated zone and the flux calculated by Darcy's law in the saturated zone, respectively. They are precisely connected at the groundwater surface interface (the values ​​are continuous without jumps). The contrasting curve shows the flux distribution without applying coupled boundary conditions, which shows obvious flux discontinuities at the interface. This eliminates the discontinuity of mass flux across phase interfaces and provides a physical basis for the uniqueness and stability of subsequent reverse tracing.

[0102] Understandably, Darcy's law states that the direction of water flow in the saturation zone is determined by the head gradient, and the flow velocity is directly proportional to the head gradient and the permeability of the medium. The negative sign indicates that the water flows in the direction of decreasing head.

[0103] In this invention, considering that there is a physical boundary between the unsaturated vadose zone and the saturated aquifer at the groundwater surface, when pollutants migrate downwards from the surface, they will undergo a sudden change in the hydrodynamic phase when crossing this boundary, that is, from unsaturated flow controlled by capillary force and gravity to saturated flow dominated by hydraulic gradient, and thus experience a nonlinear change in motion characteristics. If strict mass conservation constraints are not applied at this boundary, it will lead to the false generation or disappearance of matter during the calculation process, thereby causing the environmental state values ​​in the three-dimensional heterogeneous spatiotemporal grid to lose physical meaning.

[0104] Therefore, in this embodiment, a mass conservation boundary condition is forcibly set at the physical interface between the vadose zone and the groundwater aquifer, which can be mathematically expressed as follows:

[0105]

[0106] Among them, subscript This indicates an approach from above the groundwater surface, i.e., from the unsaturated zone side; This indicates approaching from below the groundwater level, i.e., from the saturation zone side. The symbol represents the identity. This equation mandates that the normal mass fluxes on both sides of the interface be strictly equal.

[0107] Understandably, the left side of the equation This represents the unsaturated vertical flux calculated above the interface using the Richards equation, with the added term representing the contribution of the gravitational component; the right side of the equation... This represents the saturated vertical flux calculated below the interface using Darcy's law. By forcibly making them equal, this invention uses the unsaturated mass flux output by the Richards equation as the inflow flux input to the Darcy's law calculation grid, thereby deterministically characterizing the cross-phase migration trajectory of pollutants from top to bottom across the phase interface. This forced boundary condition eliminates the mass non-conservation error that may occur when pollutants cross the phase interface due to abrupt changes in hydrodynamic properties, laying a physical foundation for the uniqueness and stability of subsequent reverse tracing.

[0108] Sub-step 5: Extract the first depth data vector corresponding to the surface soil layer and the second depth data vector corresponding to the vadose zone and groundwater aquifer from the three-dimensional heterogeneous spatiotemporal grid.

[0109] The first depth data vector refers to the sequence of environmental state parameters extracted vertically from a three-dimensional heterogeneous spatiotemporal grid, corresponding to all voxel units within the depth range of the surface soil layer. Each element in this vector represents the state value (such as the concentration of a specific pollutant, conductivity, etc.) of a voxel unit located at a certain depth in the shallow layer at the current moment.

[0110] The second depth data vector refers to the sequence of environmental state parameters extracted vertically from a three-dimensional heterogeneous spatiotemporal grid, corresponding to all voxel units from the lower part of the vadose zone to the depth of the groundwater aquifer. This vector covers a complete depth profile from the shallow subsurface to the deep aquifer.

[0111] Understandably, separating the data in the three-dimensional heterogeneous spatiotemporal grid into two independent data vectors according to depth intervals aims to accurately capture the precursor signals of pollutant migration from shallow to deep layers by performing independent time-series analysis and joint judgment on shallow and deep data. When the shallow layer undergoes abnormal changes while the deep layer remains normal, this often indicates an ongoing pollution infiltration event.

[0112] Sub-step 6: Calculate the rate of change of the concentration of a specific pollutant in the first depth data vector on the first time derivative, and compare it with the dynamic correction threshold.

[0113] The rate of change on the first-order time derivative refers to the instantaneous rate at which state parameters such as the concentration of specific pollutants or conductivity at shallow soil nodes change over time. It is used to quantify the dynamic evolution trend and speed of the shallow environmental state. The larger the absolute value of this rate of change, the faster the shallow environmental state is undergoing change.

[0114] For example, in this embodiment, the rate of change of state of shallow soil nodes can be approximately calculated using the following formula:

[0115]

[0116] in, This represents the current state measurement of the shallow node. This is the state measurement value at the previous sampling time; The sampling time interval is given. This formula uses a first-order backward difference approximation of the continuous-time derivative.

[0117] In this embodiment, in order to address the interference of hydrological changes caused by extreme rainfall on shallow sensor data and eliminate the resulting false alarm problem, the present invention introduces a false alarm compensation logic in the process of comparing the rate of change with the threshold.

[0118] The design concept of this false alarm compensation logic is as follows: When rainfall occurs, a large amount of rainwater infiltrates into the shallow soil through the surface, which dilutes the concentration of dissolved pollutants in the shallow soil, or changes the conductivity reading due to the difference in mineralization between the rainwater itself and the in-situ soil pore water. These changes in shallow parameters caused by rainfall are not caused by actual pollution leakage events. However, if the warning threshold remains unchanged, the system will be unable to distinguish between normal hydrological responses caused by rainfall and actual pollutant migration signals, thus generating a large number of false alarms, seriously affecting the reliability and operational efficiency of the warning system.

[0119] To address this, the present invention designs a dynamic correction threshold model. Its core idea is to synchronously acquire meteorological precipitation sequence data, calculate the equivalent concentration dilution caused by rainfall infiltration based on real-time rainfall information, and then superimpose this dilution onto a baseline warning threshold, allowing the warning threshold to dynamically adjust with changes in rainfall intensity. During rainfall, the threshold automatically rises to accommodate normal concentration fluctuations caused by rainfall; under no-rainfall conditions, the threshold reverts to the baseline level to maintain sensitivity to actual anomalies.

[0120] For example, in this embodiment, the dynamic correction threshold can be calculated using the following formula:

[0121]

[0122] in, For a moment Dynamic early warning threshold; The baseline warning threshold is set under no-rainfall conditions. This baseline threshold is determined based on historical monitoring data and statistical analysis of site background values. The rainfall dilution response coefficient is used to scale the magnitude of the effect of rainfall infiltration on shallow layer concentrations. For time delay Rainfall, including time delay This reflects the time required for rainfall to infiltrate from the surface to the actual installation depth of the shallow sensor, a delay determined by the length of the infiltration path and the permeability of the surface soil. The soil infiltration coefficient characterizes the ability of topsoil to infiltrate rainfall. This refers to the porosity of shallow soil. Figure 8 The diagram shows a comparison of false alarms in rainfall events between the dynamic correction threshold warning and the traditional fixed threshold warning in this embodiment. Figure 8 The image shows the time series of concentration change rates, a fixed threshold level, and a dynamic correction threshold curve during a typical rainfall event. The fixed threshold was frequently exceeded during the rainfall (triggering multiple false alarms, marked with a red cross); the dynamic correction threshold automatically increased with the rainfall, successfully avoiding false alarms. In other words, the dynamic warning threshold... By superimposing the rainfall dilution bias, false alarms caused by sudden changes in sensor conductivity / concentration due to rainfall are effectively eliminated.

[0123] It is understandable that in the above formula, The significance is as follows: After the rainfall is converted into actual infiltrated water volume using the infiltration coefficient, it is divided by the pore volume of the shallow soil to obtain the concentration dilution caused by rainwater infiltration. That is, when a certain amount of rainwater infiltrates into shallow soil with a specific porosity, the proportion of pore volume occupied by this infiltrated water determines the degree to which the original pore water pollutant concentration is diluted. By superimposing this dilution bias onto the baseline threshold, the warning threshold is automatically raised during rainfall, thereby avoiding misjudging normal concentration decreases caused by rainfall as pollution leakage events. In this embodiment, the time delay... The delay can be dynamically estimated based on the surface soil texture and sensor burial depth. For sandy soil, the delay is shorter; for clay soil, the delay is longer.

[0124] Sub-step 7: Perform tiered early warning judgment and adaptive sampling frequency adjustment.

[0125] After completing the above-mentioned rate of change calculation and dynamic threshold correction, this invention employs a shallow-deep joint criterion to trigger tiered early warning. The core design idea of ​​this joint criterion is: only when the rate of change of the shallow state exceeds the dynamic correction threshold, and the current value of the deep node is still within the preset safe range, is it determined to be a suspected real leakage event and the first-level tiered early warning is triggered. This dual condition requirement ensures that: if anomalies occur simultaneously in both the shallow and deep layers (e.g., in a large-scale hydrological response caused by a regional rainstorm), the system will not misjudge it as a local leakage event; and the early warning will only be triggered when the "shallow-to-deep" migration pattern, in which local anomalies occur in the shallow layer but the deep layer has not yet been affected, is captured.

[0126] For example, in this embodiment, the triggering condition for the first-level tiered warning can be as shown in the following formula:

[0127]

[0128] in, This is the indicator value for triggering the first-level warning. A value of 1 indicates that the warning has been triggered, and a value of 0 indicates that it has not been triggered. This represents the absolute value of the rate of change of the shallow node state. The aforementioned dynamic correction threshold; This represents the current state measurement value of the deep monitoring node; This is the deep security threshold, which is the upper limit of the state parameters that deep nodes should not exceed under normal conditions.

[0129] After the first-level tiered warning is triggered, in order to track the migration process of potential pollutants from shallow to deep layers, it is necessary to intensify the monitoring of deep nodes below the anomaly area, i.e., increase their sampling frequency. However, simply and indiscriminately increasing the sampling frequency of all deep nodes to the highest level will lead to the generation of a large amount of unnecessary data, exacerbating the power consumption burden and data storage pressure of edge computing nodes. To this end, this invention designs an adaptive sampling frequency adjustment mechanism based on spatial distance attenuation. Its core idea is to dynamically adjust the sampling frequency of each deep node through two dimensions: first, spatial proximity—the deeper the deep node is in three-dimensional space (mainly in the vertical projection direction) from the shallow anomaly point, the greater the increase in its sampling frequency, because this node is most likely to be affected by infiltrating pollutants first; second, the intensity of dynamic change—the larger the magnitude of the current concentration gradient field, the more likely the pollution front is to be advancing rapidly. At this time, all deep nodes in the affected area need a higher sampling frequency to capture its dynamic changes.

[0130] For example, in this embodiment, deep nodes The adaptive sampling frequency can be calculated using the following formula:

[0131]

[0132] in, For deep nodes At any moment The sampling frequency; This is the lowest baseline sampling frequency, i.e., the normal sampling frequency of deep nodes in the absence of warnings; This is the highest allowed sampling frequency, which is the upper limit of the frequency that the sensor hardware and edge node processing capabilities can support; The three-dimensional spatial coordinate vector of the shallow anomaly point that triggers the early warning; For deep nodes 3D spatial coordinate vector; The square of the three-dimensional Euclidean distance between the two; The spatial decay scale parameter of the Gaussian kernel controls the spatial influence range of the sampling frequency enhancement effect. The concentration gradient response coefficient; This represents the modulus of the current concentration gradient field, i.e., the rate of change of concentration in the direction of the fastest change in three-dimensional space.

[0133] It is understandable that in the above formula, It is a Gaussian kernel distance decay function with a range of . When deep nodes When the spatial location coincides with a shallow anomaly point (i.e. The function takes a maximum value of 1, meaning that the deepest node closest to the anomaly receives the largest frequency boost. As the distance between the deep node and the shallow anomaly increases, the function value decays exponentially in a Gaussian manner; deep nodes far from the anomaly region are almost unaffected by the frequency adjustment. Parameters The effective radiation radius of the frequency boosting effect is determined. The larger the value, the wider the range of deep nodes affected.

[0134] It is a hyperbolic tangent response function with a range of . When the concentration gradient When the concentration gradient is zero (i.e., the pollution field is completely uniform and unchanged in space), the function value approaches zero, meaning that even if a node is very close to an anomaly, there is no need to increase the sampling frequency, because there is no dynamic process of pollutant migration at this time; when the concentration gradient is zero... As the value increases, the function approaches 1, meaning that when a pollution front is advancing rapidly and the spatial concentration distribution is changing drastically, the sampling frequency should be increased by the maximum permissible amount. Through the product of these two factors, dual adaptive adjustment of the sampling frequency is achieved in both the spatial dimension (distance attenuation) and the dynamic change dimension (gradient response).

[0135] The system calculates the sampling frequency value for each deep node based on the above formula, generates a corresponding dynamic sampling control command, and sends it to the edge controller of the physical sensor corresponding to each deep node to adjust its wake-up frequency and sampling clock pulse.

[0136] Sub-step 8: Addressing edge node memory overflow scenarios caused by high-concurrency dynamic sampling.

[0137] In this embodiment, when the dynamic sampling control command increases the sampling frequency of a deep node to exceed a preset high-frequency threshold, the edge computing node corresponding to that deep node will face the risk of memory overflow and network bandwidth congestion due to a surge in data volume. To effectively control the resource consumption of edge nodes while ensuring that abnormal peak data is not lost, this invention further designs a data flow control mechanism based on differential incremental hashing.

[0138] The core idea of ​​this flow control mechanism is that, in high-frequency sampling mode, not every sampling cycle contains substantial new information. When the pollution front has not yet reached the deep-seated node, the data from multiple consecutive sampling cycles may be almost identical; only when there is a substantial change in the environmental state (such as the arrival of the pollution front, a sudden increase in concentration, etc.) will the data show meaningful changes. Therefore, the degree of change in data between adjacent sampling cycles can be quantified to determine whether data needs to be uploaded to the cloud for storage: only when the data undergoes substantial changes is the upload triggered; otherwise, the old data is directly overwritten locally, thus achieving a balance between data preservation and resource conservation.

[0139] Specifically, when the dynamic sampling control command increases the sampling frequency of deep nodes to exceed a preset high-frequency threshold, the system issues a data flow control block division command to the corresponding edge nodes. The edge nodes, based on this block division command, divide the continuously sampled data stream into data blocks of fixed byte size. For each data block, the edge nodes perform a customized differential incremental hash calculation.

[0140] For example, in this embodiment, the data change metric between adjacent sampling periods can be calculated using the following formula:

[0141]

[0142] in, For a moment Collected data blocks; MD5 hash digest function maps data blocks of arbitrary length to fixed-length 128-bit hash values; This is the bitwise XOR operator, which performs a bitwise XOR operation on two hash values ​​of the same length. For Hamming distance operation, calculate the number of bits with a value of 1 in the XOR result, that is, the total number of different bits between two hash values; The hash Hamming distance between adjacent data blocks is used to measure the substantial change in data between adjacent sampling periods.

[0143] It is understandable that when the data blocks in two adjacent sampling periods have completely identical content, their hash values ​​are completely identical, and the XOR result is always zero. This corresponds to the Hamming distance. When the content of a data block undergoes minor changes, the hash value will change significantly due to the avalanche effect of the hash function, and the Hamming distance will increase accordingly. By setting an appropriate change tolerance threshold, it is possible to distinguish between normal minor fluctuations in data and substantial mutations.

[0144] Edge nodes maintain a dynamic threshold verification table. Only when the calculated differential incremental hash value exceeds the current tolerance of the verification table will the data block be written to the transmission buffer pool and uploaded to the cloud; otherwise, the data block will be overwritten locally.

[0145] For example, in this embodiment, the flow control strategy based on data change measurement can be as follows:

[0146]

[0147] in, This is the tolerance threshold for data changes, i.e., the current tolerance value in the dynamic threshold verification table.

[0148] Understandably, this flow control mechanism based on differential incremental hashing, in high-frequency sampling mode, only triggers data upload when there is a substantial change in the monitored data (such as a pollution front reaching deep nodes causing a sudden increase in concentration), while local overwriting is performed directly when the data is stable and static. This strategy effectively controls the bandwidth consumption and storage pressure brought by high-frequency sampling while ensuring that key abnormal peak data is not lost, thereby avoiding memory overflow at edge nodes. Figure 9 This example shows a comparison chart of edge node cache usage. Figure 9 The comparison showed that the cache quickly approached the overflow threshold under all upload strategies, while the cache usage remained within a safe range under the flow control strategy of this solution. It can be seen that the edge node data flow control mechanism of this solution can effectively suppress cache expansion and avoid edge node memory overflow in high-frequency sampling scenarios. Figure 10 This embodiment illustrates the differential incremental hashing and dynamic threshold verification upload decision diagram. Figure 10 The demonstration shows that stable data blocks are locally overwritten, while abnormal peak data blocks are uploaded after exceeding a dynamic threshold. It can be seen that this solution does not simply compress data, but determines and retains key abnormal peaks based on substantial data changes, thereby achieving the technical effect of not losing key data and not uploading redundant data.

[0149] Step 3: When the bottom layer data of the three-dimensional heterogeneous spatiotemporal grid in the second dataset meets the preset deep anomaly threshold, extract the hydrogeological parameters in the first dataset and execute the reverse hydrogeochemical convection diffusion tracing logic to generate the three-dimensional coordinates of the leakage source, and output the three-dimensional coordinates of the leakage source and the second dataset.

[0150] The deep anomaly threshold is a pre-set reference value used to determine whether the environmental state parameters of the bottom layer of the three-dimensional heterogeneous spatiotemporal grid (i.e., the groundwater aquifer region) are abnormally high. When the concentration of a specific pollutant or other key indicators of the bottom voxel unit exceeds this threshold, it indicates that the pollutant has migrated from the shallow layer to the deep aquifer region, and source tracing needs to be initiated to locate the source of pollution.

[0151] The reverse hydrogeochemical convection-diffusion tracing logic refers to a deterministic physical source tracing method that uses classical hydrogeological equations as a physical constraint framework and reverses the time dimension by performing reverse extrapolation (i.e., taking negative values ​​for the time variable and reversing the velocity vector of the convection term) to trace the source path and starting location of pollutants backward from deep anomaly detection points towards the surface.

[0152] The three-dimensional coordinates of the leak source refer to the coordinates in three-dimensional space of the starting position of the pollutants when they first entered the underground environment, determined through reverse tracing calculations.

[0153] The following will provide a detailed explanation of the sub-steps included in this step:

[0154] Sub-step 1: Extract hydrogeological parameters from the first dataset as tracking boundary parameters.

[0155] The tracing boundary parameters refer to a set of key hydrogeological parameters used to constrain and drive the solution of the inverse tracing equations. These parameters collectively define the groundwater flow field and the migration conditions of pollutants. The tracing boundary parameters include the groundwater flow direction vector, the groundwater velocity scalar, the soil porosity scalar, and the permeability coefficient scalar.

[0156] The groundwater flow direction vector refers to the direction vector of groundwater flow in three-dimensional space, determined by the hydraulic head gradient. The direction determines the direction in which the water head decreases the fastest.

[0157] Groundwater velocity scalar refers to the actual flow rate of groundwater in porous media, which is calculated by dividing Darcy flux by effective porosity.

[0158] Soil porosity scalar refers to the proportion of pore volume to the total soil volume. This parameter determines the effective space available for the migration of water and pollutants in the soil.

[0159] The permeability coefficient is a scalar parameter characterizing the ability of soil or rock to conduct water, and is related to the aforementioned saturated hydraulic conductivity. Numerically equivalent.

[0160] In this embodiment, the aforementioned tracking boundary parameters can be extracted from in-situ sensor measurement data and borehole sampling data in the first dataset, or they can be obtained through inversion calculation in combination with geophysical exploration data.

[0161] Sub-step 2: Using the deep anomaly threshold trigger point as the initial source and sink term, construct the reverse convection diffusion equation.

[0162] Once the deep anomaly threshold is triggered, the system uses the deep high-concentration voxel that triggered the anomaly as the starting point for reverse tracing (i.e., the initial source-sink term), and performs reverse deduction in the time dimension to construct the reverse convection-diffusion partial differential equation.

[0163] In the forward physics process, pollutants migrate and diffuse downstream from the leak source, with concentration distribution spreading over time. The core idea of ​​reverse tracing is to reverse the time variable, making the inverse time variable... This means that time reversal is mathematically achieved. In the inverse equation, the velocity vector sign of the convection term (i.e., the term describing the directional transport of pollutants with the water flow) is opposite to that of the forward equation, causing the probability density field to propagate in the opposite direction of the groundwater flow field, gradually reversing from the downstream detection point to the upstream source.

[0164] Meanwhile, in the reverse simulation, the reduced effective migration velocity caused by the adsorption of pollutants by soil particles during their migration to deeper layers also needs to be considered. In forward migration, the adsorption of pollutants by the soil makes the actual transport velocity of pollutants lower than the groundwater flow velocity; in the reverse simulation, in order to accurately estimate the starting location and migration distance of pollutants, the reverse flow velocity needs to be increased accordingly to compensate for this hindering effect.

[0165] For example, in this embodiment, the reverse convection diffusion tracking equation can be expressed as follows:

[0166]

[0167] in, The spatial probability density field for reverse tracing represents the pollutants in the reverse time... The probability density of a time that may occur at various locations in space is initially set to a high probability density value at deep anomaly detection points and a zero or near-zero value at other locations. It is a reverse time variable, which increments in the past from the current moment when the exception was triggered; The hydrodynamic dispersion coefficient tensor contains longitudinal and lateral dispersion components, describing the diffusion and distribution effect of pollutants during their migration with water flow due to medium inhomogeneity and molecular diffusion. This is the effective reverse velocity vector after correction by the retardation factor; It is a spatial pruning operator used to exclude physically impenetrable regions from the solution domain.

[0168] Sub-step 3: Calculate the geochemical blocking factor using real-time pH sensor data and soil parameters to correct the reverse flow velocity.

[0169] In this embodiment, in order to improve the distance estimation accuracy of reverse tracing, the geochemical retardation factor is dynamically calculated by using the soil pH value contained in the first dataset, and the reverse flow velocity is corrected by the retardation factor.

[0170] The design motivation lies in the fact that the adsorption capacity of soil particles for pollutants varies significantly under different pH conditions. Under strongly acidic or alkaline conditions, the adsorption capacity of certain metal ions can change by orders of magnitude. Using a fixed resistance factor in reverse tracing would lead to severe distance estimation biases in sites with large pH spatial variability. Therefore, dynamically calculating the resistance factor using real-time pH sensing data allows reverse velocity correction to adapt to the actual geochemical conditions at different locations within the site.

[0171] For example, in this embodiment, the pH-dependent blocking factor and the effective reverse flow rate can be calculated using the following formula:

[0172]

[0173]

[0174] in, It is a pH-dependent lag factor, representing the lag factor of the actual migration rate of pollutants relative to the groundwater flow rate; Soil dry bulk density, which is the mass of dry matter per unit volume of soil. Porosity; , , These are empirical fitting coefficients, calibrated using batch adsorption experimental data from the laboratory. This is a real-time pH sensor measurement; This is the Darcy velocity vector of groundwater calculated using the aforementioned Darcy's law.

[0175] It is understandable that in the above formula, Soil adsorption partition coefficient A quadratic polynomial fit to the nonlinear relationship between pH and the distribution coefficient. The physical meaning of this distribution coefficient is the ratio of the amount of pollutants adsorbed per unit mass of soil to the concentration of pollutants in pore water under equilibrium conditions. This quadratic polynomial fit can express... The nonlinear characteristics of pH variation—for example, for most heavy metal pollutants, The value is higher under neutral to slightly alkaline conditions (strong adsorption) and lower under strongly acidic conditions (weak adsorption). Restriction factor The higher the value, the stronger the soil's adsorption of pollutants, and the slower the actual transport rate of pollutants compared to the groundwater flow rate. In the reverse calculation, the Darcy velocity is used... Divide by The product of these factors effectively accelerates the reverse tracing process, thereby compensating for the velocity loss caused by adsorption during forward migration. This approach avoids the tracing distance shortage error caused by using fixed flow rate parameters in traditional tracing methods.

[0176] In this embodiment, the empirical fitting coefficient , , It can be calibrated separately for different types of specific pollutants. When there are multiple target pollutants in the monitoring area, the system can calculate the blocking factor of each pollutant separately and perform reverse tracing independently.

[0177] Sub-step 4: Solve the reverse convection-diffusion equation within a three-dimensional heterogeneous spatiotemporal grid using the finite volume method.

[0178] After constructing the reverse convection-diffusion equation and all its parameters, this invention employs the finite volume method for numerical solution on a three-dimensional heterogeneous spatiotemporal grid. The choice of the finite volume method is based on its inherent mass conservation property—this method establishes a discrete set of equations by integrating the governing equations over each control volume (i.e., voxel unit), ensuring that the mass flux budget of each voxel unit is strictly conserved during the numerical solution process, which is crucial for the physical rationality of the retrospective calculation.

[0179] In the solution process, the inverse time variable With a preset time step The computation proceeds incrementally (corresponding to a gradual rollback in forward physical time), with each step advancing synchronously across the entire 3D voxel mesh. The output of each step is an updated spatial probability density field. As time progresses backward, the probability density field gradually expands and migrates upstream from the initial deep anomaly detection point, gradually forming a probability distribution peak along a physically feasible migration path.

[0180] Sub-step 5: Execute the boundary space pruning mechanism during the reverse solution process.

[0181] In this embodiment, in order to improve the computational efficiency and physical rationality of the reverse solution, the present invention performs spatial pruning operations in each iteration cycle of the reverse time-stepping calculation.

[0182] The motivation behind this spatial pruning mechanism is that in three-dimensional underground spaces, there are numerous layers of rock or dense clay with extremely low permeability, making it physically impossible for pollutants to penetrate these areas. If these areas are not excluded during the reverse solution process, the probability density field will generate spurious diffusion along the impermeable spatial branches. This wastes computational resources and reduces the accuracy of locating the true source due to the dispersion of probability density caused by spurious diffusion.

[0183] Specifically, the system reads non-contact geophysical exploration boundary data from a three-dimensional heterogeneous spatiotemporal grid and identifies absolute aquitard voxels with permeability (permeability) below a preset lower limit. In each iteration of the reverse time-step calculation, the mass exchange flux of the absolute aquitard voxels is forcibly set to zero.

[0184] For example, in this embodiment, the spatial pruning operator can be represented by the following formula:

[0185]

[0186] in, For spatial location The permeability at that location was obtained by inversion from geophysical exploration data; The minimum permeability threshold is preset; voxels below this value are considered to be absolutely waterproof. The penalty value is a sufficiently large positive number.

[0187] It is understandable that when the permeability of a certain voxel is below a threshold... At that time, by applying extremely large negative values ​​to the source and sink terms of the inverse equation. This causes the probability density of the voxel to decay rapidly to zero during the solution process (because a maximal negative source sink is equivalent to a very strong probability sink, quickly absorbing the probability density). This is equivalent to completely removing physically impenetrable regions from the solution domain, ensuring that the probability density field of the reverse tracing propagates only along physically feasible migration paths. This hard spatial pruning not only shortens the convergence time of the inversion matrix solution (because the size of the effective solution domain is reduced) but also improves the positioning accuracy (because the probability density is not dispersed onto impossible paths). Figure 11 This example shows a comparison of the convergence speed of the inverse inversion residuals before and after spatial pruning in this embodiment. Figure 11 The study compared the slow decrease in residuals without impermeable layer pruning with the rapid convergence trend after the impermeable layer flux was forced to zero in this scheme. This demonstrates that the scheme can shorten the convergence time of the inversion matrix solution by pruning path branches that are physically impossible for pollution transport through the boundary space pruning machine.

[0188] Sub-step 6: Extract the peak cluster center coordinates of the spatial probability density field to confirm the three-dimensional coordinates of the leakage source.

[0189] When the reverse solution is advanced to the depth corresponding to the Earth's surface... Reverse time endpoint At time, spatial probability density field One or more probability density peak regions are formed in three-dimensional space, and these peak regions are the most likely sources of pollution.

[0190] For example, in this embodiment, the source location coordinates can be calculated using the following formula:

[0191]

[0192] in, The estimated spatial coordinates of the pollution source, i.e., the three-dimensional coordinates of the leak source; To find the solution domain; The end time of reverse solution The spatial probability density field; It is a spatial position vector; It is a volumetric infinitesimal element.

[0193] Understandably, the above formula is essentially a weighted spatial integral of the spatial probability density field to calculate its centroid location. When the probability density field exhibits a unimodal distribution (i.e., there is only one possible source), the centroid is the most likely source location. When the probability density field exhibits a multimodal distribution (i.e., there may be multiple leakage sources), further cluster analysis (such as density-based clustering algorithms) can be performed on the probability density field to extract the cluster center coordinates of each peak as multiple candidate source coordinates.

[0194] In this embodiment, the candidate source coordinates can also be combined with surface topography and elevation data for spatial projection mapping, projecting the reverse tracing results in three-dimensional space onto the surface digital elevation model, thereby associating and matching the three-dimensional spatial coordinates of the leak source with identifiable geographical locations on the surface (such as pipelines, storage tanks, landfill areas, etc.), assisting on-site investigators in quickly locating the leaking facility.

[0195] The system ultimately outputs the three-dimensional coordinates of the leak source and a second dataset (containing a three-dimensional heterogeneous spatiotemporal grid and dynamic sampling control commands) for regulatory authorities and on-site engineers to use for subsequent precise investigation and emergency response. Figure 12 The spatial probability distribution and two-dimensional cross-sectional view of the leakage source location in this embodiment are shown.

[0196] Example 2

[0197] This embodiment discloses a soil and groundwater collaborative monitoring and hidden danger investigation system.

[0198] Specifically, the soil and groundwater collaborative monitoring and hazard investigation system can be integrated into an electronic device, which can be a terminal, server, or embedded edge computing gateway. The terminal can be an industrial tablet, laptop, or personal computer; the server can be a single server or a server cluster composed of multiple servers; the edge computing gateway can be an industrial-grade embedded device deployed on-site for performing data preprocessing and command issuance locally. When the electronic device is running, it can implement the soil and groundwater collaborative monitoring and hazard investigation method described in Embodiment 1 of this application.

[0199] In this embodiment, the soil and groundwater collaborative monitoring and hidden danger investigation system can also be integrated into multiple electronic devices. For example, the soil and groundwater collaborative monitoring and hidden danger investigation system can be integrated into multiple servers and multiple edge computing gateways, and the soil and groundwater collaborative monitoring and hidden danger investigation method in Embodiment 1 of this application can be implemented by multiple servers and edge computing gateways working together.

[0200] In this embodiment, the server can also be implemented in the form of a terminal.

[0201] Example 3

[0202] This embodiment discloses a device for coordinated monitoring of soil and groundwater and for identifying potential hazards.

[0203] This device corresponds to the method for coordinated monitoring of soil and groundwater and hazard investigation in Example 1, including:

[0204] The first acquisition device is used to acquire a multimodal environmental dataset as the first dataset based on a pre-constructed three-dimensional spatial monitoring network covering the surface soil layer, vadose zone and groundwater aquifer. The first dataset includes in-situ sensing sequence data, non-contact geophysical exploration boundary data and borehole sampling final value data.

[0205] The first acquisition device may include, but is not limited to, various in-situ sensors (such as temperature sensors, humidity sensors, conductivity sensors, pH sensors, etc.), geophysical exploration equipment (such as high-density resistivity meters, ground penetrating radar, etc.), borehole sampling and testing equipment, and a data acquisition and transmission gateway for aggregating the above-mentioned data.

[0206] The second processing unit is used to process the first dataset by executing deterministic multiphase flow coupling and deep gradient evaluation logic, and generate a second dataset containing a three-dimensional heterogeneous spatiotemporal grid and dynamic sampling control instructions.

[0207] Specifically, the second processing device can be implemented as a software module deployed on a cloud server or a local high-performance computing node. This module integrates the aforementioned spatial interpolation engine based on hydraulic distance, laboratory data calibration engine, Richards equation and Darcy's law coupled solver, cascade early warning judgment logic, adaptive sampling frequency calculation engine, and differential incremental hash flow control module.

[0208] The third source tracing and output device is used to extract hydrogeological parameters from the first dataset and execute reverse hydrogeochemical convection diffusion tracing logic to generate the three-dimensional coordinates of the leakage source when the bottom layer data of the three-dimensional heterogeneous spatiotemporal grid in the second dataset meets the preset deep anomaly threshold, and output the three-dimensional coordinates of the leakage source and the second dataset.

[0209] Specifically, the third source tracing and output device can be implemented as a reverse tracing calculation engine deployed on a cloud server. This engine integrates the aforementioned reverse convection-diffusion equation construction module, finite volume method solver, pH-dependent retardation factor calculation module, spatial pruning engine, and probability density field peak clustering and coordinate output module.

[0210] In this embodiment, the first acquisition device, the second processing device, and the third tracing and output device can be deployed in the same physical device, or they can be distributed and deployed in multiple devices to work together through a communication network.

[0211] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for synergistic monitoring and hazard investigation of soil and groundwater, characterized in that, The methods for identifying potential hazards include: Based on a pre-constructed three-dimensional spatial monitoring network covering the surface soil layer, vadose zone, and groundwater aquifer, the first dataset was obtained. The first dataset includes in-situ sensing sequence data, non-contact geophysical exploration boundary data, and borehole sampling final value data; The first dataset is processed by executing deterministic multiphase flow coupling and deep gradient evaluation logic to generate a second dataset, which contains a three-dimensional heterogeneous spatiotemporal grid and dynamic sampling control instructions. When the bottom data corresponding to the groundwater aquifer in the three-dimensional heterogeneous spatiotemporal grid meets the deep anomaly threshold, hydrogeological parameters are extracted from the first dataset, and reverse hydrogeochemical convection-diffusion tracing logic is executed based on the hydrogeological parameters to generate the three-dimensional coordinates of the leakage source. Output the three-dimensional coordinates of the leakage source and the second dataset; The three-dimensional heterogeneous spatiotemporal mesh is generated through the following steps: Using the non-contact geophysical exploration boundary data as three-dimensional spatial geometric constraints, the three-dimensional spatial monitoring network is divided into multiple three-dimensional voxel units; A hydraulic distance interpolation algorithm based on geological permeation topology maps the in-situ sensing sequence data to the corresponding three-dimensional voxel unit to generate an initial voxel state vector. The final value data of the borehole sampling is extracted as a hard calibration factor, and a weighted residual correction calculation is performed on the initial voxel state vector to generate the three-dimensional heterogeneous spatiotemporal mesh.

2. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 1, characterized in that, The hydraulic distance interpolation algorithm based on geological permeability topology includes: For each of the three-dimensional voxel units, the hydraulic distance from each in-situ sensor to that three-dimensional voxel unit is calculated, wherein, The hydraulic distance is defined as a scalar value obtained by taking the reciprocal of the permeability at each point along the spatial path from the in-situ sensor position to the target three-dimensional voxel unit and then performing a line integral. The permeability is obtained by inversion processing of the non-contact geophysical exploration boundary data. Using the negative power of the hydraulic distance as the contribution weight of each in-situ sensor to the three-dimensional voxel unit, the weighted average of the actual measurement values ​​of each in-situ sensor at the current moment is obtained to obtain the initial voxel state vector of the three-dimensional voxel unit.

3. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 2, characterized in that, The weighted residual correction calculation includes: The deviation between the final value of the borehole sampling at the time of the test sampling and the interpolated state value of the corresponding three-dimensional voxel unit at the same time is calculated and used as the calibration residual. The calibration residuals are weighted by spatial dimension using a spatial confidence tensor matrix, wherein the spatial confidence tensor matrix independently adjusts the calibration intensity in the horizontal and vertical directions; The spatially weighted calibration residual is decayed in the time dimension using an exponential time decay factor, wherein the exponential time decay factor approaches zero as the time interval between the current moment and the test sampling moment increases; The calibration residuals, after spatial weighting and time decay, are superimposed onto the interpolated state value at the current moment to obtain the calibrated environmental state value, which is then written into the three-dimensional heterogeneous spatiotemporal grid.

4. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 1, characterized in that, The deterministic multiphase flow coupling includes: A multiphase flow equation set is established on the three-dimensional heterogeneous spatiotemporal grid, and spatiotemporal evolution dynamics calculations are performed, wherein: For the aforementioned vadose zone region, the Richards equation is used to describe the unsaturated water movement flux; For the aforementioned groundwater aquifer region, Darcy's law is used to calculate the saturated water flow rate; At the physical interface between the vadose zone and the groundwater aquifer, a mass conservation boundary condition is set, which forces the normal mass flux on both sides of the interface to be strictly equal. The specific mass conservation boundary conditions are as follows: On the side above the physical interface, the unsaturated vertical flux calculated by the Richards equation includes the combined contribution of capillary driving force and gravitational components. On the lower side of the physical interface, the saturated vertical flux calculated by Darcy's law is driven by the head gradient. The mass conservation boundary condition forces the unsaturated vertical flux on the upper side to be equal to the saturated vertical flux on the lower side. The unsaturated mass flux output by the Richards equation is used as the inflow flux input to the Darcy law computation grid, thereby achieving seamless coupling of pollutants across phase interfaces.

5. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 1, characterized in that, The depth gradient evaluation logic includes: Extract a first depth data vector corresponding to the depth range of the surface soil layer along the vertical direction from the three-dimensional heterogeneous spatiotemporal grid, and a second depth data vector corresponding to the depth range from the lower part of the vadose zone to the groundwater aquifer. Calculate the rate of change of the first time derivative of the concentration of a specific pollutant in the first depth data vector; The rate of change of the first-order time derivative is compared with the dynamic correction threshold; Based on the comparison results and the current status of the deep monitoring nodes in the second depth data vector, a tiered early warning judgment is executed and the dynamic sampling control command is generated.

6. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 5, characterized in that, The dynamic correction threshold is determined in the following way: Simultaneously acquire meteorological precipitation sequence data, and based on the baseline warning threshold, superimpose the equivalent concentration dilution bias caused by rainfall infiltration to obtain the dynamic correction threshold; wherein... The equivalent concentration dilution bias is calculated as follows: The actual infiltrated water volume is obtained by multiplying the rainfall with time delay by the soil infiltration coefficient, then dividing by the shallow soil porosity to obtain the concentration dilution caused by rainwater injection, and finally scaling by multiplying by the rainfall dilution response coefficient. The time delay reflects the time required for rainfall to infiltrate from the surface to the shallow sensor installation depth, and its value is determined by the infiltration path length and surface soil permeability. The tiered early warning determination adopts a combined shallow-deep criterion: When the absolute value of the rate of change of the first time derivative exceeds the dynamic correction threshold, and the current state value of the deep monitoring node in the second depth data vector is lower than the preset deep safety threshold, the first-level tiered warning is triggered. When anomalies occur simultaneously in both the shallow and deep layers, the first-level tiered early warning will not be triggered.

7. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 1, characterized in that, The logic for executing reverse hydrogeochemical convection-diffusion tracing includes: Using hydrogeological parameters as the tracking boundary parameters, and using the three-dimensional voxel units in the three-dimensional heterogeneous spatiotemporal grid that satisfy the deep anomaly threshold as the initial source and sink terms, a reverse convection-diffusion equation is constructed. The reverse convection-diffusion equation is solved on a three-dimensional heterogeneous spatiotemporal grid to obtain the spatial probability density field; The peak cluster center coordinates of the spatial probability density field are extracted to determine the three-dimensional coordinates of the leakage source.

8. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 7, characterized in that, The reverse convection-diffusion equation is constructed by reversing the time variable into a reverse time variable, which increases in the past direction from the current moment when the anomaly is triggered; the reverse convection-diffusion equation includes a diffusion term, a convection term, and a source-sink term, wherein: The dispersion term describes the diffusion and distribution effect of pollutants caused by medium inhomogeneity and molecular diffusion using the hydrodynamic dispersion coefficient tensor. The convection term describes the propagation of the spatial probability density field in the opposite direction of the groundwater flow field using an effective reverse velocity vector; The source and sink terms include spatial pruning operators used to exclude physically impenetrable regions from the solution domain; The initial condition for the spatial probability density field is to take a high probability density value at the three-dimensional voxel unit corresponding to the initial source and sink terms, and a zero or near-zero value at other positions.

9. The method for coordinated monitoring and hazard investigation of soil and groundwater according to claim 7, characterized in that, The reverse convection diffusion equation is solved using the finite volume method. By integrating the reverse convection diffusion equation on each three-dimensional voxel unit, a discrete set of equations is established. The solution is then gradually advanced with a preset reverse time step to ensure that the mass flux balance of each three-dimensional voxel unit is conserved. The peak cluster center coordinates of the extracted spatial probability density field include: When the spatial probability density field has a single-peak distribution, the spatial probability density field is weighted and spatially integrated in the solution domain to calculate its centroid position, which is used as the three-dimensional coordinates of the leakage source. When the spatial probability density field has a multi-peak distribution, density-based clustering analysis is performed on the spatial probability density field to extract the cluster center coordinates of each peak as the three-dimensional coordinates of multiple candidate leakage sources.

10. A soil and groundwater collaborative monitoring and hidden danger investigation system, characterized in that, The hazard detection system includes: processor; The memory stores a computer program, which, when executed by a processor, implements the method for coordinated monitoring and hazard investigation of soil and groundwater as described in any one of claims 1 to 9.

Images