Mine heavy metal pollution simulation and prediction method and system based on digital twinning

Abstract

The application provides a mine heavy metal pollution simulation prediction method and system based on digital twinning, relates to the technical field of digital twinning, and first acquires a multi-source environmental monitoring data set containing heavy metal concentration sampling values and spatial coordinate information of a target mine, injects the multi-source environmental monitoring data set into a pre-constructed mine digital twinning body, generates a discrete distribution of heavy metal concentration spatial scatter point sets in the digital twinning body according to spatial coordinate information mapping, then performs three-dimensional space continuous field reconstruction processing, generates a heavy metal three-dimensional distribution field, further calls a heavy metal migration and diffusion simulation model based on a physical process drive, performs multi-medium coupling migration and diffusion simulation on the three-dimensional distribution field in the digital twinning body, and generates a space-time evolution prediction result containing heavy metal space distribution prediction fields and pollution diffusion front advancing trajectories of future time periods. The application realizes deep fusion of monitoring data and physical simulation in the digital twinning body, and improves the accuracy of mine heavy metal pollution prediction.

Description

Technical Field

[0001] This invention relates to the field of digital twin technology, and more specifically, to a method and system for simulating and predicting heavy metal pollution in mining areas based on digital twins. Background Technology

[0002] Heavy metal pollution in mining areas is a major environmental problem that restricts the ecological environment protection and sustainable development of mining areas. Heavy metal pollutants migrate and spread in the soil, water and sediment of mining areas through various pathways such as surface runoff, underground infiltration and atmospheric deposition, forming a complex spatiotemporal evolution pattern, which poses a serious threat to the surrounding ecosystem and human health.

[0003] In existing technologies, methods for simulating and predicting heavy metal pollution in mining areas mainly fall into the following categories: The first category is extrapolation prediction methods based on statistical regression models. These methods utilize historical monitoring data to establish statistical regression relationships between heavy metal concentrations and time or spatial variables, and then extrapolate and predict future pollution conditions. However, these methods treat the heavy metal migration process as a purely statistical data pattern, completely ignoring the physicochemical processes that pollutants undergo during migration and transformation between different media. This leads to a sharp decline in the reliability of the prediction results when pollution source conditions change or the prediction time span is large.

[0004] The second category is the physical model prediction method based on numerical simulation. This type of method uses partial differential equations such as the convection-diffusion equation to describe the migration and diffusion process of heavy metals in a single medium, and performs numerical solutions using the finite element method or finite difference method. Although this type of method has a certain physical mechanism basis, it usually simplifies the mining area as a homogeneous single-medium environment, and does not fully consider the coupled migration and transformation process of heavy metals in multiple media such as soil, surface water, groundwater and sediment. Moreover, the construction of the initial field of the model often relies on manual interpolation or empirical assignment, and cannot effectively integrate the spatial heterogeneity information contained in the actual discrete monitoring data, resulting in a large deviation between the simulation results and the actual pollution distribution.

[0005] The third category is spatial interpolation prediction methods based on geographic information systems. These methods utilize spatial interpolation algorithms such as Kriging interpolation and inverse distance weighting to estimate heavy metal concentrations in unmonitored areas based on existing monitoring data. While these methods can generate continuous spatial distribution prediction results, they are essentially purely data-driven spatial extrapolation methods. They lack any characterization of the physical processes of heavy metal migration and diffusion, cannot predict the dynamic advancement of pollution diffusion fronts, and cannot generate pollution evolution sequences with temporal continuity. Summary of the Invention

[0006] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a method for simulating and predicting heavy metal pollution in mining areas based on digital twins, the method comprising: Obtain a multi-source environmental monitoring data set of the target mining area, wherein the multi-source environmental monitoring data set includes heavy metal concentration sampling values ​​of multiple monitoring points and spatial coordinate information of each monitoring point; The multi-source environmental monitoring data set is injected into a pre-constructed digital twin of the mining area. Based on the spatial coordinate information of each monitoring point, the heavy metal concentration sampling value of each monitoring point is mapped to the corresponding spatial position of the digital twin of the mining area, generating a spatially scattered set of heavy metal concentration points that are discretely distributed within the digital twin of the mining area. The three-dimensional spatial continuous field reconstruction process is performed on the spatial scatter set of heavy metal concentration to generate the three-dimensional distribution field of heavy metals inside the digital twin of the mining area. A physical process-driven heavy metal migration and diffusion simulation model is invoked to perform multi-media coupled migration and diffusion simulation of the three-dimensional distribution field of the heavy metal in the digital twin of the mining area, generating a three-dimensional evolution sequence of the heavy metal in the digital twin of the mining area in the future time period. Based on the three-dimensional evolution sequence, a spatiotemporal evolution prediction result of heavy metal pollution in the target mining area is generated. The spatiotemporal evolution prediction result of heavy metal pollution includes the predicted field of heavy metal spatial distribution and the advancing trajectory of the pollution diffusion front in each future time period.

[0007] Furthermore, embodiments of the present invention also provide a digital twin-based simulation and prediction system for heavy metal pollution in mining areas, comprising: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to execute the above-described digital twin-based method for simulating and predicting heavy metal pollution in mining areas by executing the machine-executable instructions.

[0008] Based on the above, firstly, a multi-source environmental monitoring data set containing spatial coordinate information is injected into a pre-constructed digital twin of the mining area. According to the spatial coordinate information, the heavy metal concentration sampling values ​​of each monitoring point are accurately mapped to the corresponding spatial location of the digital twin, generating a discretely distributed spatial scatter set of heavy metal concentrations. On this basis, the spatial scatter set of heavy metal concentrations is reconstructed into a three-dimensional continuous field to generate a three-dimensional distribution field of heavy metals inside the digital twin. This step transforms the discrete monitoring sampling points into a continuous three-dimensional concentration field. Then, a heavy metal migration and diffusion simulation model based on physical processes is invoked to simulate the multi-media coupling migration and diffusion of this three-dimensional distribution field within the digital twin, generating a three-dimensional evolution sequence of heavy metals in the future time period. The above simulation process is not carried out in an abstract mathematical space, but is executed in a digital twin that is completely consistent with the actual geometric structure and media distribution of the mining area. Moreover, the initial field comes from real monitoring data reconstructed through continuous field rather than being manually set, so that the physical process of multi-media coupling migration and diffusion can be unfolded on a real spatial skeleton, thereby ensuring the physical authenticity and spatial accuracy of the simulation results. Finally, based on the three-dimensional evolution sequence, a spatiotemporal evolution prediction result is generated, which includes the predicted field of heavy metal spatial distribution and the trajectory of pollution diffusion front in each future time period. This spatiotemporal evolution prediction result has both a continuous evolution sequence in the time dimension and the ability to dynamically track the diffusion front in the spatial dimension. This allows decision-makers not only to know the concentration distribution of heavy metals in various spatial locations in the mining area at any future time, but also to clearly grasp the direction and speed of the pollution diffusion front. This achieves a panoramic and accurate description of the spatiotemporal evolution process of heavy metal pollution in the mining area, significantly improving the effectiveness of the prediction results for decision-making on environmental governance in the mining area. Attached Figure Description

[0009] Figure 1 This is a schematic diagram of the execution flow of the digital twin-based simulation and prediction method for heavy metal pollution in mining areas provided in this embodiment of the invention.

[0010] Figure 2 This is a logical schematic diagram of the method for simulating and predicting heavy metal pollution in mining areas based on digital twins provided in this embodiment of the invention.

[0011] Figure 3 This is a schematic diagram of exemplary hardware and software components of the digital twin-based heavy metal pollution simulation and prediction system for mining areas provided in an embodiment of the present invention. Detailed Implementation

[0012] Figure 1 This is a flowchart illustrating a digital twin-based method for simulating and predicting heavy metal pollution in mining areas, provided in one embodiment of the present invention. The following is a detailed description in conjunction with the attached diagram. Figure 2This method is described in detail. It can be applied to the field of mining area environmental monitoring and pollution control, and is particularly suitable for using digital twin technology to perform three-dimensional spatiotemporal dynamic simulation and prediction of the migration and diffusion process of heavy metal pollution in multi-media environments within mining areas. An example of a non-ferrous metal mining area is used for illustration, but this application scenario is merely illustrative and does not constitute a limitation on the scope of protection of this application.

[0013] Step S110: Obtain a multi-source environmental monitoring data set of the target mining area. The multi-source environmental monitoring data set includes heavy metal concentration sampling values ​​of multiple monitoring points and spatial coordinate information of each monitoring point.

[0014] Specifically, the data interface of the mining area environmental monitoring data management system reads various environmental monitoring data collected from the target mining area within a specified monitoring period. This multi-source environmental monitoring data set is stored in a standardized data table structure, with each record representing a single sampling of a monitoring point. The fields in the data table include: a monitoring point identifier, which uniquely identifies a sampling point; a sampling timestamp, recording the specific time of sampling; spatial coordinate information, including the X, Y, and Z coordinates of the monitoring point in the unified coordinate system of the mining area; a media category identifier, used to distinguish different environmental media such as soil, groundwater, surface water, or atmospheric deposition; and a heavy metal concentration sampling value, recording the concentration of a specific heavy metal element (such as lead or cadmium) in the sample. In this embodiment, the spatial coordinate information of the monitoring point is represented in the form of three-dimensional coordinates (x, y, z), where the z-coordinate reflects the depth or height of the sampling point. For soil and groundwater monitoring, the z-coordinate is the depth below the surface; for surface water and atmospheric deposition monitoring, the z-coordinate is the height above or below the surface. The digital twin of the mining area is a three-dimensional spatial gridded model whose spatial coordinate system is consistent with the coordinate system of the monitoring points, so as to facilitate subsequent spatial mapping.

[0015] Step S120: Inject the multi-source environmental monitoring data set into the pre-constructed digital twin of the mining area, and map the heavy metal concentration sampling value of each monitoring point to the corresponding spatial location of the digital twin of the mining area based on the spatial coordinate information of each monitoring point, thereby generating a spatially distributed set of heavy metal concentration points within the digital twin of the mining area.

[0016] Through the data injection interface, the multi-source environmental monitoring data set obtained in step S110 is input into the pre-constructed digital twin model of the mining area. This digital twin is a three-dimensional spatial gridded model, where each grid node corresponds to a spatial coordinate position. Internally, the model is divided into soil, groundwater, surface water, and atmospheric media layers according to the physical stratification characteristics of the environmental media, each layer having its own grid and attributes. Based on the spatial coordinate information of each monitoring point, the specific location of that point in the three-dimensional spatial coordinate system of the mining area's digital twin is located, and the media layer type corresponding to that point is further identified. The heavy metal concentration sampling value (incremental value after background correction) and its media category identifier are assigned to a scatter object at that spatial location. By traversing all monitoring points, all injected scatter objects are aggregated into a discretely distributed set, namely the spatial scatter set of heavy metal concentrations, where each element contains spatial coordinates, a media category identifier, and a concentration increment value.

[0017] Step S121: Extract the medium category identifier corresponding to the heavy metal concentration sampling value of each monitoring point in the multi-source environmental monitoring data set. The medium category identifier is used to distinguish between soil monitoring media, groundwater monitoring media, surface water monitoring media, and atmospheric deposition monitoring media.

[0018] The media category identifier field in each record of the multi-source environmental monitoring dataset is parsed to extract its corresponding media category. This media category identifier can be an integer code, such as code 1 for soil, code 2 for groundwater, code 3 for surface water, and code 4 for atmospheric deposition. This identifier allows for the differentiation of heavy metal data in different environmental media, facilitating subsequent stratification and media parameter configuration.

[0019] Step S122: Perform media stratification processing on the multi-source environmental monitoring data set according to the media category identifier to generate a soil media sampling record subset, a groundwater media sampling record subset, a surface water media sampling record subset, and an atmospheric deposition media sampling record subset.

[0020] Based on media category identification, the records in the multi-source environmental monitoring dataset were grouped. Records categorized as media category 1 were extracted to form a soil media sampling record subset; records categorized as media category 2 were extracted to form a groundwater media sampling record subset; records categorized as media category 3 were extracted to form a surface water media sampling record subset; and records categorized as media category 4 were extracted to form an atmospheric deposition media sampling record subset. Data records within each subset were sorted according to sampling timestamps to maintain chronological order.

[0021] Step S123: Perform soil background correction processing on the heavy metal concentration sampling values ​​of each monitoring point in the soil medium sampling record subset, and perform a difference operation between the heavy metal concentration sampling value of each monitoring point and the preset soil background concentration benchmark value to generate a set of heavy metal concentration increments in the soil medium. Each element in the set of heavy metal concentration increments in the soil medium corresponds to the heavy metal concentration deviation of a monitoring point relative to the soil background.

[0022] For each record in the subset of soil media sampling records, its heavy metal concentration sampling value is read. Simultaneously, the baseline soil concentration value corresponding to that monitoring point is retrieved from the mining area environmental baseline database. This baseline value represents the average background level of heavy metals in the soil of that area under conditions unaffected by mining pollution. The sampling value is subtracted from this baseline value to obtain the heavy metal concentration increment at that point relative to the baseline. If the sampling value is lower than the baseline value, the increment is negative; however, considering that pollution is usually characterized by an increase in concentration, negative increments are set to 0. The calculated increment values ​​are associated with the spatial coordinates, media category identifier, and monitoring point identifier of that point to construct a set of heavy metal concentration increments in the soil media.

[0023] Step S124: Perform groundwater background correction processing on the heavy metal concentration sampling values ​​of each monitoring point in the groundwater medium sampling record subset, and perform differential calculation between the heavy metal concentration sampling values ​​of each monitoring point and the preset groundwater background concentration benchmark value to generate a set of groundwater medium heavy metal concentration increments.

[0024] The groundwater medium data is processed using the same logic as in step S123. Sampled values ​​are read, and the corresponding groundwater background concentration benchmark value is subtracted to generate the groundwater medium heavy metal concentration increment. If the difference is less than 0, it is set to 0. The increment is associated with the spatial coordinate information and medium category identifier of that location to construct a set of groundwater medium heavy metal concentration increments.

[0025] Step S125: Perform surface water background correction processing on the heavy metal concentration sampling values ​​of each monitoring point in the surface water medium sampling record subset, and perform differential calculation between the heavy metal concentration sampling values ​​of each monitoring point and the preset surface water background concentration benchmark value to generate a set of surface water medium heavy metal concentration increments.

[0026] The same logic is used to process surface water media data. The difference between the sampled values ​​and the baseline concentration of surface water is calculated to generate a set of heavy metal concentration increments in surface water media.

[0027] Step S126: Perform atmospheric background correction processing on the heavy metal concentration sampling values ​​of each monitoring point in the atmospheric deposition medium sampling record subset, and perform differential calculation between the heavy metal concentration sampling values ​​of each monitoring point and the preset atmospheric background concentration benchmark value to generate a set of heavy metal concentration increments in the atmospheric deposition medium.

[0028] The same logic is used to process atmospheric deposition media data. The difference between the sampled values ​​and the baseline atmospheric concentration is calculated to generate a set of heavy metal concentration increments in atmospheric deposition media.

[0029] Step S127: Map the sets of heavy metal concentration increments in the soil medium, the groundwater medium, the surface water medium, and the atmospheric deposition medium to the corresponding spatial locations of the soil medium layer, groundwater medium layer, surface water medium layer, and atmospheric medium layer in the digital twin of the mining area, according to their respective medium category identifiers and the spatial coordinate information of each monitoring point. This generates a spatial scatter set of heavy metal concentration increments for each of the four medium layers in the digital twin of the mining area. Each scatter point contains spatial coordinates, a medium category identifier, and a concentration increment value.

[0030] In the digital twin of the mining area, soil, groundwater, surface water, and atmosphere are modeled as independent media layers, each with its own three-dimensional spatial grid structure. Each scatter point in the four incremental sets generated in steps S123 to S126 is mapped to the corresponding grid position in the corresponding media layer based on its spatial coordinates and media category identifier. After mapping, a series of scatter points marked with concentration increment values ​​are distributed in the soil media layer, and corresponding scatter points are also distributed in the groundwater, surface water, and atmospheric media layers, forming four independent subsets, namely, the spatial scatter point sets of heavy metal concentration increments.

[0031] Step S128: Align the heavy metal concentration increments corresponding to the four media layers in the spatial scatter set of heavy metal concentration increments with the media layer positions in the unified three-dimensional spatial coordinate system of the mining area digital twin. Record the heavy metal concentration increments of different media layers at the same spatial coordinate projection position together with their corresponding media category identifiers to generate a spatial scatter set of heavy metal concentrations that is discretely distributed in the mining area digital twin and marked with the media category.

[0032] The four subsets generated in step S127 are spatially aligned in a unified three-dimensional coordinate system. For a specific spatial location, if there are scatter points in both the soil and groundwater media layers, the two scatter points at that location are merged into a single entry, and the concentration increment values ​​of the soil layer, the groundwater layer, and their respective media category identifiers are recorded simultaneously. Through the above alignment operation, a single, multi-dimensional spatial scatter point set of heavy metal concentrations containing information from multiple media layers is finally generated. Each scatter point contains the heavy metal concentration increment information corresponding to all media layers at that specific spatial location.

[0033] Step S130: Perform three-dimensional spatial continuous field reconstruction processing on the spatial scatter set of heavy metal concentration to generate a three-dimensional distribution field of heavy metals inside the digital twin of the mining area.

[0034] The spatial scatter set of heavy metal concentrations generated in step S128 is converted into a three-dimensional continuous concentration distribution field covering the entire digital twin space. Since the monitoring points (scatter points) are sparse and discretely distributed in space, it is necessary to estimate the concentration variation patterns between scatter points and in unsampled areas through spatial interpolation and reconstruction methods, thereby generating a smooth, complete, and physically consistent three-dimensional continuous distribution field.

[0035] Step S131: Perform spatial clustering segmentation on the spatial scatter point set of heavy metal concentrations. Divide the spatial scatter point set of heavy metal concentrations into multiple spatial clusters based on the spatial distribution density. Each spatial cluster contains a group of scatter points with spatially adjacent locations and similar heavy metal concentration values.

[0036] A density-based spatial clustering algorithm, such as the DBSCAN algorithm, is used to process the scatter set generated in step S128. A spatial neighborhood radius and a threshold for the minimum number of scatter points within the neighborhood are set. Scatter points whose spatial distance is less than this radius and whose concentration values ​​are within a certain range of difference are clustered into the same cluster. Through clustering operations, discrete scatter points are divided into multiple clusters with similar internal data point characteristics. These clusters correspond to areas in the mining area with similar pollution characteristics and spatial continuity, such as strong pollution diffusion centers or various branches of pollution plumes.

[0037] Step S132: Extract the local concentration gradient direction for the concentration scatter points within each spatial cluster, calculate the rate of change of heavy metal concentration values ​​within each spatial cluster along each coordinate axis in the spatial coordinate system, generate the local concentration gradient vector corresponding to each spatial cluster, construct the concentration gradient field for each spatial cluster based on the local concentration gradient vector corresponding to each spatial cluster, and estimate the concentration value at any spatial location within the spatial cluster using the spatial coordinates, concentration values, and local concentration gradient vector of the concentration scatter points within each spatial cluster through spatial interpolation methods, thereby generating the local concentration gradient field corresponding to each spatial cluster.

[0038] For each spatial cluster generated in step S131, local concentration gradient analysis is performed. This step specifically includes: first, extracting the spatial coordinates and concentration values ​​of all scattered points within the cluster; then, fitting a local linear plane using the least squares regression method to obtain the spatial trend of concentration variation, i.e., calculating the partial differential approximations of the concentration values ​​relative to the X, Y, and Z coordinates, forming a local concentration gradient vector composed of three components; finally, using the discrete scattered points within the cluster as control points and the calculated local concentration gradient vector as a directional constraint, using a radial basis function interpolation algorithm, estimating the concentration value at any grid position within the cluster, thereby generating the local concentration gradient field corresponding to the cluster.

[0039] Step S133: Perform concentration field smoothing processing on adjacent spatial clusters in the cluster boundary region, extract the concentration value difference of the local concentration gradient field corresponding to each adjacent spatial cluster in the boundary region, and use a weighted fusion method based on spatial distance to gradually eliminate the concentration difference in the boundary region, generating a smooth concentration transition zone between adjacent spatial clusters.

[0040] When there is a significant concentration jump between adjacent clusters generated in step S132 at the boundary, smoothing processing is required to ensure the physical continuity of the distribution field. This step specifically includes: extracting a set of scattered points near the boundary of adjacent clusters and calculating their respective average concentrations; constructing a transition zone in the boundary region, the width of which is determined by preset parameters; within the transition zone, employing an inverse distance weighted average algorithm, assigning higher weights to clusters closer to the point and lower weights to clusters farther away; through this weighted fusion method, the concentration distribution on both sides of the boundary achieves a smooth transition, eliminating unnatural abrupt changes.

[0041] Step S134: Spatial splicing and integration of the local concentration gradient fields corresponding to all spatial clusters and all smooth concentration transition zones are performed to generate an initial three-dimensional distribution field of heavy metals covering the entire internal space of the digital twin of the mining area. Non-negative constraint correction processing is performed on the initial three-dimensional distribution field of heavy metals, and the concentration values ​​of spatial locations with concentration values ​​less than 0 in the initial three-dimensional distribution field of heavy metals are set to 0 to generate a non-negative three-dimensional distribution field of heavy metals.

[0042] The local gradient fields of all clusters and all transition zones are stitched together to form an initial three-dimensional distribution field of heavy metals covering the entire digital twin space. Since the interpolation process may produce negative values, which is physically unrealistic, a non-negativity constraint correction is needed for the distribution field. Specifically, the concentration values ​​of all grid nodes in the distribution field are traversed. If a node's concentration value is found to be less than 0, it is forcibly modified to 0 to ensure that the concentration field is physically realizable and meaningful, thus generating a non-negative three-dimensional distribution field of heavy metals.

[0043] Step S135: Perform physical consistency correction processing on the non-heavy metal three-dimensional distribution field, adjust the concentration value of spatial location points in the non-heavy metal three-dimensional distribution field that exceed the preset theoretical upper limit of heavy metal concentration in the mining area to the theoretical upper limit of heavy metal concentration in the mining area, and generate the heavy metal three-dimensional distribution field inside the digital twin of the mining area.

[0044] In nature, pollutant concentrations do not increase indefinitely, thus requiring constraints on the extreme values ​​of the distribution field. This step involves: obtaining the theoretical upper limit of the target heavy metal element concentration from a mining area environmental database. This theoretical upper limit is typically the highest concentration value obtained through actual measurements of pollution sources in the mining area; traversing all grid nodes in the non-negative three-dimensional heavy metal distribution field, if the concentration value of a node exceeds the upper limit, forcibly modifying the node's concentration value to the upper limit; after this physical consistency correction, a physically reasonable three-dimensional heavy metal distribution field that satisfies both the non-negativity constraint and the physical upper limit constraint is finally generated.

[0045] Step S140: Call the physical process-driven heavy metal migration and diffusion simulation model to perform multi-media coupled migration and diffusion simulation of the three-dimensional distribution field of the heavy metal in the digital twin of the mining area, and generate a three-dimensional evolution sequence of the heavy metal in the digital twin of the mining area in the future time period.

[0046] Using the three-dimensional heavy metal distribution field generated in step S135 as the initial concentration condition, the heavy metal migration and diffusion simulation model is invoked. This heavy metal migration and diffusion simulation model is based on a physics-driven numerical simulation engine. Based on the convection-diffusion equation, it comprehensively considers the physicochemical processes of heavy metals in four media layers—soil, groundwater, surface water, and atmosphere—including convection, diffusion, adsorption, desorption, sedimentation, and resuspension, as well as the interfacial material exchange processes between these media layers. Through iterative solutions along the time axis, the dynamic changes of heavy metal concentration in three-dimensional space are simulated, outputting a three-dimensional concentration distribution field for multiple consecutive time steps, forming a three-dimensional evolution sequence. This step includes several sub-steps, detailing the preparation of model parameters and the specific implementation methods for coupled solutions of each media layer.

[0047] Step S141: Extract spatial distribution data of soil moisture content and spatial distribution data of soil porosity from the soil medium layer of the digital twin of the mining area, and input the spatial distribution data of soil moisture content and spatial distribution data of soil porosity into the soil medium parameter module of the heavy metal migration and diffusion simulation model to generate a spatial distribution field of heavy metal convection diffusion coefficient in the soil medium.

[0048] Two key physical parameters are read from each grid node of the soil medium layer grid in the digital twin: soil moisture content, denoted as Wc, representing the volume percentage of water in the soil; and soil porosity, denoted as Ps, representing the volume percentage of pores in the soil. After receiving these two parameters, the soil medium parameter module first calculates the effective diffusion coefficient De of soil moisture using a soil moisture characteristic curve model (such as the van Genuchten model). Then, combining the soil bulk density and the distribution coefficient of heavy metals in the soil Kd, the soil medium retardation factor R = 1 + (Dd*Kd) / Wc is calculated, where Dd is the soil dry density. Finally, the convective diffusion coefficient Dsoil is calculated as Dsoil = De / R. This calculation is performed on all grid nodes to generate a spatial distribution field Dsoil(x, y, z) of the heavy metal convective diffusion coefficient covering the entire soil medium layer.

[0049] Step S142: Extract groundwater velocity vector field and aquifer permeability coefficient spatial distribution data from the groundwater medium layer of the digital twin of the mining area, and input the groundwater velocity vector field and aquifer permeability coefficient spatial distribution data into the groundwater medium parameter module of the heavy metal migration and diffusion simulation model to generate the spatial distribution field of heavy metal convection diffusion coefficient in the groundwater medium.

[0050] The groundwater velocity vector field Vw(x, y, z) (containing Ux, Uy, Uz components) and the spatial distribution of aquifer permeability Kw(x, y, z) are read from the groundwater medium layer grid of the digital twin. The groundwater medium parameter module first calculates the dispersion based on the permeability and aquifer thickness, then calculates the mechanical dispersion coefficient of heavy metals in the groundwater: Dm = A*Vw + B*Vw, where A and B are the longitudinal and transverse dispersions, respectively. Molecular diffusion Dmol_water is also considered. The effective diffusion coefficient in the groundwater is Dw = Dm + Dmol_water. For the adsorption process, the groundwater retardation factor Rw = 1 + (Dw*Kdw) / nw is calculated, where nw is the effective porosity and Kdw is the distribution coefficient in the groundwater. Finally, the spatial distribution field of the heavy metal convective diffusion coefficient in the groundwater medium is Deff_w = Dw / Rw. This output covers the entire groundwater medium layer.

[0051] Step S143: Extract the surface runoff velocity vector field and the spatial distribution data of heavy metal adsorption coefficient in riverbed sediment from the surface water medium layer of the digital twin of the mining area. Input the surface runoff velocity vector field and the spatial distribution data of heavy metal adsorption coefficient in riverbed sediment into the surface water medium parameter module of the heavy metal migration and diffusion simulation model to generate the spatial distribution field of heavy metal convection diffusion coefficient in the surface water medium.

[0052] The surface runoff velocity vector field Vs(x, y, z) and the heavy metal adsorption coefficient distribution Ksed(x, y, z) in the riverbed sediment heavy metal medium layer grid are read from the digital twin. The surface water medium parameter module calculates the turbulent diffusion coefficient Dturb=At*Vs based on the flow velocity. The longitudinal dispersion coefficient Dlong=Al*Vs is calculated. The total diffusion coefficient in surface water Dsurf=Dturb+Dlong+Dmol_water. The adsorption of riverbed sediment is corrected by an interfacial exchange coefficient Kexch=Ksed / (1+Ksed*Ased), where Ased is the sediment surface area. Finally, the spatial distribution field Dsurf_eff of the heavy metal convective diffusion coefficient in the surface water medium is generated.

[0053] Step S144: Extract near-surface wind speed vector field and spatial distribution data of atmospheric boundary layer turbulent diffusion coefficient from the atmospheric medium layer of the digital twin of the mining area, and input the near-surface wind speed vector field and spatial distribution data of atmospheric boundary layer turbulent diffusion coefficient into the atmospheric medium parameter module of the heavy metal migration and diffusion simulation model to generate the spatial distribution field of heavy metal convection diffusion coefficient in the atmospheric medium.

[0054] The near-surface wind speed vector field Va(x, y, z) and the atmospheric boundary layer turbulent diffusion coefficient distribution Kturb(x, y, z) are read from the atmospheric medium layer mesh of the digital twin. The atmospheric medium parameter module calculates wind shear based on wind speed to determine the intensity of atmospheric mechanical turbulence. Assuming that the diffusion of heavy metals in the atmosphere is mainly controlled by turbulence, the effective atmospheric diffusion coefficient Dair = Kturb + Dmol_air. Simultaneously, for settleable particulate matter, the settling velocity Vdep = (Pp*g*dp^2) / (18*Mu), where Pp is the particulate density, dp is the particle diameter, and Mu is the aerodynamic viscosity. Finally, the atmospheric convective diffusion coefficient field Dair_eff = Dair / (1 + Vdep / Va). This output covers the entire atmospheric medium layer.

[0055] Step S145: Call the media coupling module of the heavy metal migration and diffusion simulation model to perform heavy metal flux exchange processing between soil and groundwater at the interface between soil and groundwater, between soil and surface water at the interface between soil and surface water, between soil and atmosphere at the interface between soil and atmosphere, and between surface water and groundwater at the interface between surface water and groundwater, generating a multi-media interface heavy metal coupling flux set.

[0056] The media coupling module is invoked. This module handles the mass exchange between media layers, specifically as follows: First, at the soil-groundwater interface (usually the water table), the flux Jsg = Ksg * (Cs - Cg) is calculated based on the concentration difference and exchange coefficient between the two media layers at the interface, where Ksg is the exchange coefficient. Similarly, for the soil-surface water interface, the flux Jss = Kss * (Cs - Csw); for the surface water-groundwater interface, the flux Jsw = Ksw * (Csw - Cg); and for the atmosphere-soil interface, the flux Jas = Kas * (Ca - Cs). The case of soil dust entering the atmosphere and returning through sedimentation will also be considered. Finally, a set containing the fluxes of all interfaces is generated.

[0057] Step S146: Using the three-dimensional distribution field of heavy metals as the initial concentration field, the spatial distribution fields of the heavy metal convection diffusion coefficients in the soil medium, the groundwater medium, the surface water medium, the atmospheric medium, and the multi-media interface heavy metal coupling flux set are input into the convection diffusion solution module of the heavy metal migration and diffusion simulation model. The multi-media convection diffusion equations are coupled and iteratively solved within the digital twin of the mining area according to a preset time step. The spatial distribution of heavy metal concentrations in the soil medium layer, groundwater medium layer, surface water medium layer, and atmospheric medium layer is updated within each time step to generate a three-dimensional evolution sequence of heavy metals in the digital twin of the mining area within the future time period.

[0058] The convection-diffusion solution module discretizes the equations on a three-dimensional regular grid using the finite volume method. For the soil medium layer, its concentration update follows the convection-diffusion equation. In the iteration, the soil layer concentration is updated as C_new = C_old + delta_t * (convection term + diffusion term + interface flux term - attenuation term). The groundwater layer, surface water layer, and atmosphere layer use the same iterative logic. At the end of each time step, the updated concentration data of the four layers are written back to the digital twin, forming a new distribution field. This iterative calculation continues until the preset total simulation time is reached, and the distribution field of each time step is stored in chronological order to form the final three-dimensional evolution sequence. This three-dimensional evolution sequence covers all medium layers, all grid nodes, and all time steps.

[0059] Step S150: Generate the spatiotemporal evolution prediction results of heavy metal pollution in the target mining area based on the three-dimensional evolution sequence. The spatiotemporal evolution prediction results of heavy metal pollution include the predicted field of heavy metal spatial distribution and the advancing trajectory of the pollution diffusion front in each future time period.

[0060] The three-dimensional evolution sequence generated in step S146 is used as input data, and a series of post-processing operations are performed to generate the final prediction result. This post-processing includes: first, superimposing the spatial distribution of each medium layer within each time step to form a comprehensive pollution distribution; second, dividing the superimposed field into regions of different pollution levels according to a preset threshold to form a labeled field; then, extracting and stitching the pollution diffusion front by comparing the differences between the labeled fields of adjacent time steps; further, calculating the front's advancing velocity and direction, and performing spatial interpolation to obtain a continuous velocity field; next, tracing streamlines from the pollution source to obtain the pollution diffusion path; and finally, generating the front's advancing trajectory through time integration and envelope surface construction.

[0061] Step S151: Extract the spatial distribution of heavy metal concentration in soil, groundwater, surface water, and atmosphere corresponding to each time step from the three-dimensional evolution sequence. Spatially superimpose the spatial distribution of heavy metal concentration in the four media at the same time step in the unified three-dimensional spatial coordinate system of the mining area digital twin to generate a heavy metal spatial distribution superposition field corresponding to each time step.

[0062] The data records for each time step in the 3D evolution sequence are traversed. For each time step, four independent 3D arrays are extracted from the record, corresponding to the spatial distribution of heavy metal concentrations in the soil, groundwater, surface water, and atmospheric media layers, respectively. These four 3D arrays have identical dimensions, and each array element corresponds to a grid node in the digital twin. The concentration values ​​of the soil, groundwater, surface water, and atmospheric layers within the time step are summed at the same grid node to obtain a scalar value representing the total heavy metal concentration at that node. The calculation result is stored as a new 3D array, i.e., the superimposed field of heavy metal spatial distribution for that time step. All grid nodes are traversed to complete the construction of the superimposed field for that step. The above process is repeated until all time steps in the 3D evolution sequence have been processed.

[0063] Step S152: Perform pollution level classification on the superimposed field of heavy metal spatial distribution corresponding to each time step, mark the spatial area where the heavy metal concentration value exceeds the preset heavy metal pollution risk threshold as the pollution risk area, mark the spatial area where the heavy metal concentration value exceeds the preset heavy metal pollution warning threshold as the pollution warning area, generate the heavy metal pollution zoning label field corresponding to each time step, arrange the heavy metal pollution zoning label fields corresponding to all time steps in chronological order, and generate the heavy metal spatial distribution prediction field for each future time period.

[0064] For each overlay field generated in step S151, a pollution level classification operation is performed. Two threshold parameters are set: the first is the heavy metal pollution risk threshold Tr, and the second is the heavy metal pollution warning threshold Ta, where Ta is greater than Tr. Each grid node in the overlay field is traversed, and the total concentration value C of that node is obtained (this value is the summation result, with units consistent with the input). If C is greater than or equal to Ta, the region of that node is marked as a pollution warning region and assigned a corresponding code (e.g., number 3); if C is greater than or equal to Tr but less than Ta, the region of that node is marked as a pollution risk region and assigned a corresponding code (e.g., number 2); if C is less than Tr, the region of that node is marked as a pollution-free region and assigned a corresponding code (e.g., number 1). All node marking results are stored in a new three-dimensional array as the heavy metal pollution zoning labeling field for that time step. All zoning labeling fields corresponding to all time steps are arranged in chronological order to form a time series, which is the predicted field for the spatial distribution of heavy metals in future time periods.

[0065] Step S153: Extract the heavy metal pollution zone labeling fields of two adjacent time steps from the predicted heavy metal spatial distribution field of each future time period. Take the area boundary where the spatial location of the same pollution zone labeling category changes in two adjacent time steps as pollution diffusion front segments. Perform spatial continuity splicing on all pollution diffusion front segments between two adjacent time steps to generate pollution diffusion fronts between adjacent time steps.

[0066] From the predicted field sequence generated in step S152, two consecutive time steps T and T+1 are extracted and denoted as labeled fields PT and PT+1, respectively. Using an edge detection algorithm (such as the Canny algorithm) in image processing, the boundary polygon sets of the pollution risk area and pollution warning area in PT and PT+1 are extracted, respectively. For each boundary polygon in PT, a similarly shaped boundary polygon with a displaced position is found in PT+1, and the displacement vector between the boundary polygons is calculated. The displaced boundary polygons in PT are then concatenated with their corresponding boundary polygons in PT+1 to form a continuous spatial surface patch, which is the pollution diffusion front segment. All front segments between adjacent time steps are summarized to form a continuous pollution diffusion front.

[0067] Step S154: Arrange the pollution diffusion fronts between all adjacent time steps in time series, extract the spatial displacement vector of each spatial location point on each pollution diffusion front between adjacent time steps, and normalize the spatial displacement vector of each spatial location point on each pollution diffusion front according to the time interval between the two time steps corresponding to the pollution diffusion front to generate the diffusion propulsion velocity vector of each spatial location point on each pollution diffusion front.

[0068] Iterate through all pollution diffusion fronts generated in step S153. For each front, uniformly sample a series of spatial locations. For each sampling point, calculate its spatial coordinates X(T) at time step T and X(T+1) at time step T+1. Calculate the displacement vector V = X(T+1) - X(T) from X(T) to X(T+1). Divide the displacement vector V by the time step value (e.g., T_dt) to obtain the diffusion propulsion velocity vector V_speed = V / T_dt for that sampling point. Store all these velocity vectors for subsequent analysis.

[0069] Step S155: The diffusion propulsion velocity vectors of all spatial locations on the entire pollution diffusion front are spatially interpolated in the three-dimensional spatial coordinate system of the digital twin of the mining area to generate a continuous diffusion velocity vector field covering the entire spatial region where the pollution diffusion front is located.

[0070] Using the discrete velocity vector points obtained in step S154 as control points, a continuous vector field is generated on the three-dimensional grid of the mining area digital twin using a three-dimensional kriging interpolation algorithm or an inverse distance weighted interpolation algorithm. Each grid node in this vector field has a velocity vector value, representing the trend direction of pollution diffusion at that location, and forms the basis of the entire frontal advance velocity field.

[0071] Step S156: In the continuous vector field of diffusion velocity, starting from the preset initial spatial location point of the pollution source, perform three-dimensional streamline tracking processing along the vector direction of each spatial location point in the continuous vector field of diffusion velocity to generate multiple pollution diffusion streamlines. Each pollution diffusion streamline represents the movement path of heavy metal pollution in space from the initial spatial location point of the pollution source along the diffusion direction.

[0072] One or more preset initial spatial locations of pollution sources are selected as the starting points for streamline tracing. In the continuous velocity vector field generated in step S155, a fourth-order Runge-Kutta numerical integration method is used to advance step-by-step along the velocity vector direction from the starting point, with the step size controlled within a certain range. This advancement process generates a curved trajectory, which is a pollution diffusion streamline. Applying different random perturbations from different starting points or the same starting point can generate multiple streamlines, covering possible pollution diffusion paths.

[0073] Step S157: Perform time dimension integration processing on all pollution diffusion streamlines in the three-dimensional space of the digital twin of the mining area. Along the path direction of each pollution diffusion streamline, accumulate the product of the diffusion propagation velocity vector of each spatial location point and the time integration step with a preset time integration step to obtain the spatial coordinates corresponding to each time integration node on each pollution diffusion streamline. Perform spatial envelope surface construction processing on the spatial coordinates corresponding to the same time integration node on all pollution diffusion streamlines to generate the propagation trajectory of the pollution diffusion front.

[0074] For each pollution diffusion streamline generated in step S156, discretization sampling is performed along the streamline path to obtain a series of discrete nodes. For each node, the total path length from the starting point to the current node is calculated, and the time required to reach the node is roughly estimated by dividing the total length by the average speed. Nodes with the same arrival time in all streamlines are collected to form a point set. A closed polygonal surface is constructed by applying a 3D convex hull algorithm or an alpha-shape algorithm to the above point set. This polygonal surface represents the spatial location of the pollution front at the corresponding time point. Arranging the front positions at all time points in chronological order constitutes the advancement trajectory of the pollution diffusion front.

[0075] Step S210: Obtain the time-series forecast data of rainfall in the target mining area and the surface runoff generation and confluence model of the mining area. Input the time-series forecast data of rainfall in the mining area into the surface runoff generation and confluence model of the mining area to generate the time-series forecast curve of surface runoff flow and the time-series evolution data of surface runoff inundation range of the target mining area in the future time period.

[0076] Rainfall time-series forecast data for the target mining area is obtained from meteorological departments or hydrological monitoring stations. This data is presented in a list format, recording the predicted rainfall intensity values ​​for multiple future time points (e.g., recorded at hourly or 6-hour intervals). This data is then input into a pre-established surface runoff model for the mining area. Based on parameters such as topographic slope, soil type, and land use type, this model analyzes the relationship between rainfall and surface runoff, calculates the curve of surface runoff flow over time, and uses data from a digital elevation model to calculate the spatial inundation extent vector map or raster data for runoff at various time points.

[0077] Step S220: Input the time-series prediction curve of surface runoff flow and the time-series evolution data of surface runoff inundation range into the surface water medium layer of the mining area digital twin. Combined with the digital elevation model, calculate the spatial distribution data of surface runoff water depth and the surface runoff velocity vector field for each time step, and update the corresponding data in the surface water medium layer. Based on the updated surface runoff velocity vector field and surface runoff water depth spatial distribution data, recalculate the spatial distribution field of heavy metal convection diffusion coefficient in the surface water medium in the heavy metal migration and diffusion simulation model.

[0078] The surface runoff flow curve and inundation range data generated in step S210 are imported into the surface water medium layer of the mining area's digital twin. Using GIS spatial analysis functions, the inundation range data is overlaid with the digital elevation model to determine the specific water depth changes and temporal distribution in the inundated area. Using Manning's formula or other hydraulic models, combined with surface topography and flow data, the surface runoff velocity vector field is calculated. These updated water depth and velocity data replace the original corresponding parameters in the surface water layer. Then, based on the new velocity vector and water depth data, the calculation process in step S143, i.e., the surface water medium parameter module calculation, is re-executed to obtain a new spatial distribution field of heavy metal convection diffusion coefficients in the surface water medium.

[0079] Step S230: Using the three-dimensional distribution field of heavy metals at the current time step as the initial concentration field, input the recalculated spatial distribution field of the heavy metal convection diffusion coefficient in the surface water medium and the updated surface runoff velocity vector field into the convection diffusion solution module of the heavy metal migration and diffusion simulation model, and re-simulate the migration and diffusion process of heavy metals in the digital twin of the mining area under rainfall conditions to generate a corrected three-dimensional evolution sequence of heavy metals under rainfall conditions.

[0080] The three-dimensional heavy metal distribution field generated in step S135 is used as the initial concentration condition. The updated spatial distribution field of the surface water medium convection diffusion coefficient and the updated surface runoff velocity vector field (especially the velocity vector) obtained in step S220 are used as input parameters and then re-aggregated into the multi-media coupled convection diffusion solution module. The iterative calculation process described in step S146 is then executed. In this simulation, since the parameters of the surface water layer have been updated (considering the influence of rainfall), the final generated three-dimensional evolution sequence will more accurately reflect the heavy metal dynamics during the rainfall and runoff periods, i.e., it is the heavy metal three-dimensional evolution sequence corrected under rainfall conditions.

[0081] Step S240: Based on the three-dimensional evolution sequence of heavy metals corrected under the rainfall conditions, update the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period, and generate the prediction results of the spatiotemporal evolution of heavy metal pollution under the rainfall scenario.

[0082] Using the corrected three-dimensional evolution sequence generated in step S230 as input data, the entire calculation process from steps S151 to S157 is repeated. Since the input data already reflects the impact of rainfall, the newly generated pollution distribution field and frontal advancement trajectory naturally take into account runoff and scouring effects under rainfall. Through this series of operations, the final result is a prediction of the spatiotemporal evolution of heavy metal pollution under a rainfall scenario. Users can use this result to assess the extent and risk of pollution diffusion during the rainy season.

[0083] Step S310: Obtain the heavy metal adsorption characteristics data of the target mining area soil, which includes the spatial distribution data of soil cation exchange capacity, the spatial distribution data of soil organic matter content, and the spatial distribution data of soil pH.

[0084] By reviewing the soil survey report of the mining area or through laboratory sampling and analysis, data on the cation exchange capacity (CEC), soil organic matter content, and soil pH of the target mining area are obtained. This data is typically provided in tabular form and linked to the specific spatial location or sampling unit of the mining area (sometimes provided in raster or vector polygon format). This data is then coordinate-registered and formatted for later use in the digital twin.

[0085] Step S320: Input the heavy metal adsorption characteristic data of the mining area soil into the soil medium parameter module of the heavy metal migration and diffusion simulation model, perform adsorption correction processing on the heavy metal migration process in the soil medium, generate a spatial distribution field of adsorption inhibition factors for heavy metals in the soil medium based on the spatial distribution data of soil cation exchange capacity and soil organic matter content; generate a spatial distribution field of reaction rate parameters for heavy metals in the soil medium based on the spatial distribution data of soil pH; use the spatial distribution field of adsorption inhibition factors to correct the convective transport velocity of heavy metals in the soil medium, and use the spatial distribution field of reaction rate parameters to characterize the adsorption reaction terms in the migration and diffusion simulation, thereby realizing the adsorption correction of the heavy metal migration and diffusion process in the soil medium.

[0086] The soil adsorption characteristic data obtained in step S310 is input into the soil medium parameter module to calculate two key parameters: the adsorption retardation factor and the reaction rate constant. The adsorption retardation factor is usually calculated using empirical formulas or linear isothermal adsorption models. It determines the ratio of the migration rate of heavy metals in the soil medium to the water flow rate, i.e., R = 1 + (rho_d * Kd) / theta, where rho_d is the soil dry density and Kd is the partition coefficient. This partition coefficient Kd depends on the soil's cation exchange capacity and organic matter content. Similarly, the hydrolysis reaction rate or adsorption rate constant of heavy metals can be determined based on the pH value, which can be used to construct the reaction rate parameter field. These parameters are then applied to the convection-diffusion equation in the soil medium to perform adsorption corrections to the original physical processes.

[0087] Step S330: Using the three-dimensional distribution field of heavy metals at the current time step as the initial concentration field, replace the original spatial distribution field of heavy metal convection diffusion coefficient in the soil medium with the spatial distribution field of heavy metal convection diffusion coefficient in the soil medium after adsorption correction, and input it into the convection diffusion solution module of the heavy metal migration and diffusion simulation model to re-simulate the migration and diffusion process of heavy metals in the soil medium within the digital twin of the mining area, and generate the three-dimensional evolution sequence of heavy metals after soil adsorption correction.

[0088] The soil medium convective diffusion coefficient field (including adsorption correction information) calculated in step S320 is used. The heavy metal distribution field generated in step S135 is used as the initial condition, and the adsorption-corrected soil medium convective diffusion coefficient replaces the original soil medium parameters. This is input into the convective diffusion solution module to re-simulate the migration and diffusion process of heavy metals in the soil medium. Because adsorption is considered, the model can more accurately reflect the occurrence state and mobility of heavy metals in the soil. Finally, a three-dimensional evolution sequence of heavy metals after soil adsorption correction is generated, which can be used to evaluate the natural decay law.

[0089] Step S340: Based on the soil adsorption-corrected three-dimensional evolution sequence of heavy metals, update the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period, and generate a prediction result of the spatiotemporal evolution of heavy metal pollution combined with soil adsorption.

[0090] Using the new evolution sequence generated in step S330, the post-processing procedures described in steps S150 to S157 are executed again. The resulting pollution prediction results reflect the influence of heavy metal adsorption / desorption processes in the soil; for example, the migration speed of the pollution plume may slow down or exhibit lag, and the advancement speed of the pollution front may also decrease. The output is a prediction of the spatiotemporal evolution of heavy metal pollution that incorporates the effects of soil adsorption.

[0091] Step S410: Obtain the time-series monitoring data of the groundwater level in the target mining area and the groundwater recharge and discharge parameters of the mining area. Based on the time-series monitoring data of the groundwater level in the mining area and the groundwater recharge and discharge parameters of the mining area, determine the groundwater level fluctuation curve and groundwater flow field change trend of the target mining area in the future time period.

[0092] Continuous historical groundwater level data was collected from various groundwater monitoring wells within the mining area, along with recharge and discharge parameters of the groundwater system (such as precipitation infiltration coefficient, river seepage, and artificial drainage volume). Using a groundwater flow dynamics model (such as MODFLOW or FEFLOW) and combining this data, groundwater level changes over future periods can be predicted, including the trend of rise or fall and the magnitude of fluctuations. Simultaneously, this groundwater flow dynamics model can also calculate the vector field changes in the direction of groundwater flow (i.e., streamlines) during the predicted period. This trend prediction forms the basis for subsequent calculations of heavy metal migration in groundwater.

[0093] Step S420: Input the groundwater level fluctuation curve and groundwater flow field change trend into the groundwater medium layer of the mining area digital twin. Combined with the spatial distribution data of aquifer permeability coefficient, calculate the groundwater velocity vector field for each time step through the groundwater flow model, and update the aquifer head spatial distribution data and groundwater velocity vector field. Based on the updated groundwater velocity vector field and aquifer head spatial distribution data, recalculate the spatial distribution field of heavy metal convection diffusion coefficient in the groundwater medium in the heavy metal migration and diffusion simulation model.

[0094] The groundwater level fluctuation curve and flow field trend obtained in step S410 are input into the groundwater flow model. Combined with parameters such as aquifer permeability coefficient and storage coefficient in the digital twin, the groundwater head distribution and velocity vector at each time step are solved. This data is then updated in the groundwater medium layer of the digital twin. Finally, using the updated flow field and aquifer head data, the spatial distribution field of the heavy metal convective diffusion coefficient in the groundwater medium is recalculated.

[0095] Step S430: Using the three-dimensional distribution field of heavy metals at the current time step as the initial concentration field, input the recalculated spatial distribution field of the heavy metal convection diffusion coefficient in the groundwater medium and the updated groundwater velocity vector field into the convection diffusion solution module of the heavy metal migration and diffusion simulation model, and re-simulate the migration and diffusion process of heavy metals in the groundwater medium within the digital twin of the mining area to generate a three-dimensional evolution sequence of heavy metals after groundwater level change correction.

[0096] The updated parameters of the groundwater medium calculated in step S420 are applied to the heavy metal migration model. The current heavy metal distribution is used as the initial condition, and the updated groundwater flow velocity and convection-diffusion coefficient field are used for simulation calculation. Finally, a three-dimensional evolution sequence of heavy metals after groundwater level change correction is generated (such as reflecting the impact of groundwater flow changes during dry or wet seasons on the morphology of the pollution plume).

[0097] Step S440: Based on the three-dimensional evolution sequence of heavy metals corrected by the groundwater level change, update the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period, and generate a prediction result of the spatiotemporal evolution of heavy metal pollution combined with the groundwater level change.

[0098] Similar to step S340, the new evolution sequence is post-processed to generate the final prediction result, taking into account the impact of groundwater level changes on groundwater pollution migration.

[0099] Step S510: Obtain the time-series prediction data of the atmospheric wind field and the topographic roughness distribution data of the target mining area, input the time-series prediction data of the atmospheric wind field and the topographic roughness distribution data of the mining area into the atmospheric boundary layer flow field calculation model, and generate the corrected wind speed vector field and spatial distribution data of the turbulence diffusion coefficient near the surface of the target mining area in the future time period.

[0100] Hourly atmospheric wind field time-series prediction data for the target mining area are extracted from the output of the meteorological forecasting service system or meteorological model (such as the Weather Research and Forecasting Model WRF) for the target mining area. This dataset is stored in a structured three-dimensional array, containing three-dimensional wind speed vectors (including components in the U, V, and W directions) and wind direction data covering the spatial extent of the mining area at multiple future time points. Simultaneously, topographic roughness distribution data for the mining area is obtained from a surface feature database. This topographic roughness distribution data is stored in raster form, where the value of each grid cell represents the surface roughness length at that location, which directly affects the intensity of near-surface wind field and turbulence.

[0101] The two sets of data are simultaneously input into a pre-established atmospheric boundary layer flow field calculation model (such as the CALMET diagnostic wind field model). This model uses terrain roughness, initial wind field, and thermodynamic parameters as boundary conditions. By solving the mass continuity equation and momentum equation, it diagnoses and adjusts the input wind field data, generating more refined and physically consistent near-surface wind speed vector field and turbulent diffusion coefficient spatial distribution data. The output of the atmospheric boundary layer flow field calculation model is a three-dimensional gridded dataset covering the atmospheric medium layer of the mining area's digital twin, containing the corrected wind speed vector (U_c, V_c, W_c) and turbulent diffusion coefficient at each grid node.

[0102] Step S520: Load the corrected wind speed vector field and the spatial distribution data of the turbulent diffusion coefficient into the atmospheric medium layer of the mining area digital twin, and update the near-surface wind speed vector field and the spatial distribution data of the atmospheric boundary layer turbulent diffusion coefficient for each time step in the atmospheric medium layer.

[0103] Through the data integration interface, the output results of the atmospheric boundary layer flow field calculation model generated in step S510 (corrected wind speed vector field and spatial distribution data of turbulent diffusion coefficient) are loaded into the data storage area corresponding to the atmospheric medium layer of the mining area digital twin. The loading process includes matching and replacing the calculation results with the atmospheric medium layer grid in the digital twin node by node and time step according to the time step (e.g., per hour) and the position coordinates of each spatial node. After the operation is completed, the original old wind field and turbulent diffusion coefficient data in the atmospheric medium layer are overwritten by the updated data obtained from this prediction.

[0104] Step S530: Based on the updated near-surface wind speed vector field and the spatial distribution data of atmospheric boundary layer turbulent diffusion coefficient, recalculate the spatial distribution field of heavy metal convective diffusion coefficient in the atmospheric medium in the heavy metal migration and diffusion simulation model.

[0105] The heavy metal migration and diffusion simulation model reads updated near-surface wind speed vector fields (U_c, V_c, W_c) and spatial distribution data of turbulent diffusion coefficients from the atmospheric medium layer of the digital twin. The atmospheric medium parameter module in the model converts the above data into convective diffusion coefficient parameters of heavy metals in the atmospheric medium. Specifically, the convective diffusion coefficient vector field is defined as D_v = D_0 + (U_c, V_c, W_c), where D_0 is the basic molecular diffusion coefficient. The turbulent diffusion coefficient is directly used to correct the dispersion term in the diffusion equation. Through the above calculations, a spatial distribution field of heavy metal convective diffusion coefficients covering all grid nodes of the atmospheric medium layer is generated.

[0106] Step S540: Using the three-dimensional distribution field of heavy metals at the current time step as the initial concentration field, input the recalculated spatial distribution field of heavy metal convection diffusion coefficient in the atmospheric medium and the updated near-surface wind speed vector field into the convection diffusion solution module of the heavy metal migration and diffusion simulation model, and re-simulate the migration and diffusion process of heavy metals in the atmospheric medium within the digital twin of the mining area to generate a three-dimensional evolution sequence of heavy metals after atmospheric wind field correction.

[0107] The three-dimensional distribution field of heavy metals generated in step S135 (for the initial time step) or the distribution field of the previous time step (for cyclic simulation) is used as the initial concentration condition input. The spatial distribution field of the heavy metal convective diffusion coefficient in the atmospheric medium calculated in step S530 and the near-surface wind speed vector field updated in step S520 are input into the convective diffusion solution module of the heavy metal migration and diffusion simulation model. The solution module, based on the discretized form of the convection-diffusion equation (such as the finite volume method or the finite difference method), iteratively performs time-step integral solutions on the three-dimensional grid of the digital twin atmospheric medium layer. In each iteration, the model uses the new wind speed vector field to update the convection term and the new turbulent diffusion coefficient field to update the diffusion term, thereby updating the atmospheric medium heavy metal concentration value of each grid node. After multiple time-step iterations, a new set of three-dimensional heavy metal concentration distribution sequences reflecting the influence of atmospheric wind field changes is generated, namely, the atmospheric wind field-corrected three-dimensional heavy metal evolution sequence.

[0108] Step S550: Based on the heavy metal three-dimensional evolution sequence corrected by the atmospheric wind field, update the predicted field of heavy metal spatial distribution and the advancing trajectory of the pollution diffusion front for each future time period, and generate a prediction result of the spatiotemporal evolution of heavy metal pollution combined with changes in the atmospheric wind field.

[0109] Using the corrected 3D evolution sequence of the atmospheric wind field generated in step S540 as input data, the pollution level classification, front tracking, and streamline trajectory generation operations described in steps S151 to S157 are executed again. When calculating the frontal advancement trajectory, because the new evolution sequence has different concentration distribution patterns (e.g., pollution plumes deflecting towards a specific wind direction or intensity changes), the starting point and direction of the streamline tracking in step S156 will change accordingly. The generated pollution source diffusion path and frontal advancement trajectory will also reflect the influence of the wind field. Finally, the updated prediction field and trajectory information are visualized and output to form a prediction result of the spatiotemporal evolution of heavy metal pollution combined with changes in the atmospheric wind field. The pollution distribution topographic map, contour map, and frontal advancement animation included will accurately reflect the driving effect of wind and aerodynamics on the diffusion of pollution in the mining area.

[0110] For example, the method may further include: step S610: obtaining spatial distribution data of surface cover type and spatial distribution data of root depth of vegetation in the target mining area, and generating spatial distribution field of surface cover permeability coefficient and spatial distribution field of root absorption depth of vegetation in the mining area based on the spatial distribution data of surface cover type and spatial distribution data of root depth of vegetation in the mining area.

[0111] Spatial distribution data of land cover types in the target mining area are obtained through high-resolution remote sensing image interpretation or field surveys. This data is stored in raster form, with the value of each grid cell representing the land cover type at that point (e.g., exposed soil and rock, grassland, shrubland, built-up area, water body, etc.). Simultaneously, spatial distribution data of root depth of vegetation in the mining area is obtained through vegetation ecological survey data. This is also a raster layer, with the value of each grid cell representing the maximum root depth or effective absorption depth of the main vegetation type at that point.

[0112] The spatial distribution data of land cover types were digitized, mapping each type to its corresponding permeability coefficient value. For example, the permeability coefficient of exposed soil and rock was set as K1, grassland as K2, shrubland as K3, and built-up areas as K4. Through this mapping, the land cover type raster was converted into a spatial distribution field of permeability coefficients for mining area land cover, denoted as Perm_surf(x, y, z=land surface). For the spatial distribution data of vegetation root absorption depth, the root depth was estimated as a depth range based on vegetation type and soil type, constructing a vertical absorption depth field TopRoot(x, y, z). These two distribution fields describe the spatial influence of surface conditions and root systems on the migration and absorption processes of heavy metals in the soil. The legality, accuracy, and relevant usage permits of the above calculation data have been compliant with relevant national and industry laws and regulations, and all data collection involved have obtained explicit authorization.

[0113] Step S620: Load the spatial distribution field of the surface cover permeability coefficient of the mining area onto the interface between the soil medium layer and the surface water medium layer of the digital twin of the mining area, and update the spatial distribution of the interface permeability coefficient of the heavy metal flux exchange treatment between the soil medium and the surface water medium.

[0114] Through the digital twin interface, the spatial distribution field of surface cover permeability coefficients, Perm_surf(x, y), generated in step S610, is assigned to the interface between the soil layer and the surface water layer in the digital twin. This interface is a two-dimensional manifold surface defined in three-dimensional space (approximately the surface elevation surface). For each grid node on the interface, its original interface permeability coefficient parameter is replaced with the value of Perm_surf(x, y). After this update, the flux exchange calculation between the soil medium and the surface water medium in the simulation of heavy metal migration and diffusion will use the new permeability coefficient, which conforms to the surface cover type, thereby more accurately simulating the surface runoff and soil water exchange process.

[0115] Step S630: Load the spatial distribution field of the vegetation root absorption depth into the soil medium layer of the digital twin of the mining area, add the vegetation root absorption sink term to the heavy metal migration and diffusion simulation model, and determine the spatial range and vertical depth of the vegetation root absorption sink term in the soil medium layer based on the spatial distribution field of the vegetation root absorption depth.

[0116] The spatial distribution field of vegetation root absorption depth, TopRoot(x, y, z), generated in step S610, is loaded into the soil medium layer. When simulating heavy metal migration and diffusion, a sink term S_veg(x, y, z, t) is added to the convection-diffusion equation. The value of this sink term is determined by the root absorption depth field TopRoot and the current heavy metal concentration in the soil. Specifically, if a node (x, y, z) in the soil is within the effective absorption depth range indicated by TopRoot, the model sets a negative source term (i.e., sink term) at that node, the strength of which is related to the concentration and vegetation absorption rate parameters at that point; if the node is outside the absorption depth range, the sink term is 0. In this way, the process of vegetation roots reducing heavy metals in the soil through absorption is characterized in the mathematical model.

[0117] Step S640: Using the current time step's three-dimensional distribution field of heavy metals as the initial concentration field, input the updated spatial distribution of interface permeability coefficients and the sink term of vegetation root absorption into the convection-diffusion solution module of the heavy metal migration and diffusion simulation model, and re-simulate the interface exchange process of heavy metals in the digital twin of the mining area between soil medium and surface water medium, as well as the absorption process of vegetation roots in the soil medium, to generate a three-dimensional evolution sequence of heavy metals after surface cover and vegetation correction.

[0118] The three-dimensional heavy metal distribution field (initial time) generated in step S135 is used as the initial concentration condition input. The interface permeability coefficient distribution field updated in step S620 and the vegetation root absorption sink term set in step S630 are substituted into the convection-diffusion solution module of the heavy metal migration and diffusion simulation model. During the iteration process, the solution module not only considers convection and diffusion, but also calculates the heavy metal flux exchange at the soil-surface water interface based on the new interface permeability coefficient, and at each node of the soil layer, the mass balance equation of that node is modified according to the presence or absence of the sink term S_veg. Through iterative advancement step by step, a new set of three-dimensional heavy metal concentration distribution sequences is obtained, namely, the three-dimensional evolution sequence of heavy metals after land cover and vegetation correction. This three-dimensional evolution sequence of heavy metals reflects the enhancement or weakening of soil-surface water exchange caused by land cover, and the fixation effect of vegetation absorption on soil heavy metals.

[0119] Step S650: Based on the three-dimensional evolution sequence of heavy metals corrected by land cover and vegetation, update the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period, and generate a prediction result of the spatiotemporal evolution of heavy metal pollution combining the effects of land cover and vegetation.

[0120] Using the three-dimensional evolution sequence of land cover and vegetation corrected in step S640 as input data, the post-processing procedures described in steps S151 to S157 are executed again. When updating the pollution spatial distribution prediction field, the concentration distribution field corrected by vegetation absorption and interface exchange will exhibit local concentration reduction characteristics, and the pollution level classification results will also be adjusted accordingly. When tracking pollution diffusion fronts, the advance speed of the front may be slowed down due to the influence of vegetation on the soil surface, and the surface diffusion path will also be guided and blocked by the distribution of land cover types. Finally, a spatiotemporal evolution prediction result of heavy metal pollution combining the effects of land cover and vegetation is generated, which can be used to evaluate the mitigation effect of control measures such as ecological restoration on pollution diffusion.

[0121] Step S710: Obtain the historical sequence of the three-dimensional distribution field of heavy metals and the corresponding historical sequence of environmental driving factors of the mining area in the historical period of the digital twin of the mining area. Use the historical sequence of the three-dimensional distribution field of heavy metals as training labels and the historical sequence of environmental driving factors of the mining area as training inputs to form a set of historical sample pairs.

[0122] Historical data from the mining area's digital twin stored in the data warehouse system over several time periods (e.g., the past year) are retrieved. The heavy metal concentration distribution sequence covering the entire 3D grid of the digital twin is extracted and denoted as C_hist(t, x, y, z), which serves as the training label. Simultaneously, historical sequences of environmental driving factors corresponding to the mining area for that time period are extracted, including rainfall time-series data, groundwater level data, wind speed and direction data, and temperature data, denoted as D_hist(t, x, y, z), which serve as the training input. The driving factors at each time step are paired with the corresponding heavy metal concentration distribution, forming a historical sample pair set {(D_i, C_i)} containing N samples, where i = 1, 2, ..., N. This historical sample pair set is used to train the machine learning model.

[0123] Step S720: Discretize the three-dimensional space within the digital twin of the mining area into multiple spatial cube units, using each spatial cube unit as a node in the graph structure and the spatial adjacency relationship between adjacent spatial cube units as an edge in the graph structure to construct a spatial topology graph structure for the mining area.

[0124] The continuous three-dimensional space of the mining area's digital twin is discretized into a set of regular spatial cube units, each with the same size. Each unit occupies a fixed coordinate range in space. Each cube unit is treated as a node in a graph, with node numbers corresponding one-to-one with the unit's spatial index. For each node, its adjacent nodes (upper, lower, front, back, left, and right) are connected by edges, forming an undirected graph structure. The number of nodes in this graph equals the total number of units in the three-dimensional space, and the number of edges equals the number of adjacent unit pairs. This spatial topology graph structure of the mining area describes the local spatial dependence of heavy metal concentration.

[0125] Step S730: Perform node feature assignment processing on the historical sequence of mining area environmental driving factors for each time step in the historical sample pair set, assign the environmental driving factor value of each spatial cube unit to the corresponding node in the mining area spatial topology graph structure, and generate a node feature map corresponding to each time step.

[0126] For each time step of the historical sample pair set, the environmental driving factor data D_i is assigned to the corresponding node in the spatial topology graph structure based on its spatial coordinates. For example, if the driving factor is rainfall, the rainfall intensity or cumulative value of a cell is assigned to that node; if the driving factor is groundwater level, the water level elevation of the node is assigned to that cell. After the assignment, each node in the graph structure has one or more feature values ​​(determined by the number of driving factors). In this way, the environmental state at each time step is represented as a node feature graph, which contains the spatial information of all nodes and their corresponding environmental states.

[0127] Step S740: Perform spatial neighborhood feature aggregation processing on the node feature map corresponding to each time step. For each node in the spatial topology map structure of the mining area, extract the node features of all neighboring nodes of the node, and perform a weighted summation of the node features of all neighboring nodes and the node features of the node itself to generate the spatial aggregation feature vector of each node at that time step.

[0128] Construct a graph neural network layer (e.g., a graph convolutional network layer) that receives node feature maps as input. For each node v, collect the node feature vectors F_u of all its neighboring nodes u, and sum them with the node v's own feature vector F_v, i.e., F_v'=sum_{uinN(v)}w_{vu}*F_u+b. Here, the weighted edges w_{vu} are learnable parameter matrices, and b is a bias term. Through this operation, each node aggregates the environmental driving information of its surrounding neighborhood, and the output node feature vector F' contains the environmental features of the node and its local region, serving as the basis for further analysis.

[0129] Step S750: Perform temporal evolution feature extraction processing on the spatial aggregation feature vector of each node in all time steps. Input the spatial aggregation feature vector of each node in different time steps into the recursive calculation module along the time axis, and calculate the recursive change of the spatial aggregation feature vector between adjacent time steps in turn to generate the spatiotemporal recursive feature sequence of each node.

[0130] The spatial aggregated feature vector sequence generated in step S740 is arranged in chronological order. For each node, its spatial aggregated feature vectors at different time steps are extracted to form a time series. This time series is input into a recurrent neural network (e.g., a Long Short-Term Memory network LSTM) or a gated recurrent unit. The recurrent neural network processes the input vectors sequentially in chronological order and maintains an internal hidden state. At each step, the recurrent neural network calculates the recursive change of the current step's features with respect to the previous step's hidden state, delta_h=f(h_{t-1}, X_t), and outputs a new hidden state. Finally, each node will obtain a set of spatiotemporal recursive feature sequences characterizing the evolution of its environmental features over time, denoted as H_seq(v)=[h_1, h_2, ..., h_T].

[0131] Step S760: Based on the feature state of the last time step in the spatiotemporal recursive feature sequence, perform recursive prediction processing on the heavy metal concentration of each node in the next time step, fuse the spatial aggregation feature vector of the last time step with the recursive change amount, and generate the predicted value of heavy metal concentration of each node in the next time step.

[0132] For each node, take the state h_T of the last time step of the spatiotemporal recursive feature sequence generated in step S750 and the spatial aggregated feature vector F'_T of that time step. Fuse these two values ​​using a fully connected prediction head to generate the predicted value formula: y_{pred}=W*concat(h_T, F'_T)+b, where W is the weight matrix and b is the bias term. This output value represents the predicted heavy metal concentration for that node in the next time step. Repeat this operation for all nodes to obtain a predicted three-dimensional distribution map containing the concentration distribution of the entire region in the next time step.

[0133] Step S770: Compare the predicted heavy metal concentration of each node at the next time step with the corresponding heavy metal concentration label value in the historical sequence of the three-dimensional distribution field of heavy metals to construct a prediction deviation feedback quantity. Iteratively adjust the weighting weight in the spatial neighborhood feature aggregation processing and the recursive parameters in the recursive calculation module through the prediction deviation feedback quantity to obtain the trained spatial neighborhood aggregation parameters and temporal recursive parameters.

[0134] The predicted concentration value C_pred generated in step S760 for the next time step is compared with the label concentration value C_true for the corresponding time step, and the loss function L = |C_pred - C_true|^2 is calculated. The gradient of the loss function with respect to the graph convolutional layer weight matrix W and the recursive parameters of the recurrent neural network is calculated using the backpropagation algorithm. An optimizer (such as Adam) is used to update the parameters iteratively to reduce the gap between the prediction and the label. Training is repeated multiple times until the loss decreases to a preset value or convergence occurs. At this point, the model training is complete, and a set of optimal values ​​is obtained as the spatial neighborhood aggregation parameters and temporal recursive parameters for the completed training.

[0135] Step S780: Obtain the real-time mining area environment driving factor sequence of the target mining area, perform the node feature assignment processing, the spatial neighborhood feature aggregation processing, and the temporal evolution feature extraction processing on the real-time mining area environment driving factor sequence, generate the real-time spatiotemporal recursive feature sequence of each node using the trained spatial neighborhood aggregation parameters and temporal recursive parameters, and generate the surrogate prediction three-dimensional evolution sequence of heavy metals based on the real-time spatiotemporal recursive feature sequence.

[0136] The system reads the sequence of environmental driving factors in the mining area (such as rainfall and wind fields from numerical weather prediction) in real time after the current moment. Following the logic of steps S730 to S750, the real-time driving factors are assigned to nodes in the topology graph, and spatial neighborhood feature aggregation is performed to obtain a spatial aggregated feature vector. This feature vector is then input into the recursive calculation module trained in step S770 to generate a real-time spatiotemporal recursive feature sequence. After processing by the prediction head, this real-time spatiotemporal recursive feature sequence outputs predicted heavy metal concentrations for several future time steps starting from the current moment. These predicted sequences constitute the proxy-predicted three-dimensional evolution sequence of heavy metals. Because the driving factor sequence data is legally sourced, compliant, and authorized, the prediction results can be continuously updated based on real-time data.

[0137] Step S790: Perform weighted fusion processing on the heavy metal concentration values ​​of the same spatial cubic unit at the same time step in the heavy metal three-dimensional evolution sequence predicted by the proxy and the heavy metal migration and diffusion simulation model generated by the heavy metal three-dimensional evolution sequence to generate a fused and corrected heavy metal three-dimensional evolution sequence.

[0138] For each time step and each spatial cube cell, the surrogate model prediction value C_agent(x, y, z, t) generated in step S780 is weighted and fused with the heavy metal migration and diffusion simulation model prediction value C_physics(x, y, z, t) generated in step S146. The specific fusion formula is: C_fused = alpha * C_physics + (1-alpha) * C_agent, where alpha is a weighting coefficient (e.g., 0.7). This fusion combines the accurate results driven by physical processes with the advantages of rapid prediction from graph neural networks, generating more efficient and reliable prediction results. Performing the above operation on all cells and time steps generates a fused and corrected three-dimensional evolution sequence of heavy metals.

[0139] Step S7100: Update the predicted field of heavy metal spatial distribution and the trajectory of pollution diffusion front for each future time period according to the fused and corrected three-dimensional evolution sequence of heavy metals, and generate the spatiotemporal evolution prediction result of heavy metal pollution assisted by proxy prediction.

[0140] Using the fused and corrected three-dimensional evolution sequence of heavy metals generated in step S790 as input, the pollution level classification and front trajectory tracking operations in steps S151 to S157 are performed again. Since this is a set of comprehensive prediction data that integrates physical models and surrogate models, the generated prediction field and front trajectory will more comprehensively reflect the actual pollution diffusion trend of the mining area, and the computational efficiency and data real-time performance are improved. The final generated surrogate prediction-assisted spatiotemporal evolution prediction results of heavy metal pollution are obtained.

[0141] Step S810: Obtain the real-time monitoring values ​​of heavy metal concentrations and the corresponding real-time monitoring values ​​of environmental driving factors in the target mining area from multiple monitoring points. Combine the real-time monitoring values ​​of heavy metal concentrations and environmental driving factors in the same time stamp into a real-time monitoring sample. Arrange the real-time monitoring samples from multiple consecutive time stamps in chronological order to form a real-time monitoring sequence.

[0142] Data is continuously collected through online monitoring equipment installed at various key locations in the mining area. This equipment includes soil heavy metal sensors, groundwater quality monitoring wells, automatic surface water monitoring stations, and atmospheric particulate matter samplers. Simultaneously, real-time data on environmental driving factors, including rainfall and wind speed, are obtained from meteorological and hydrological stations. Heavy metal concentration data collected at the same time stamp are paired with driving factor data to form a sample, denoted as sample S_t={C_t, D_t}, where C_t contains the heavy metal concentration at that time, and D_t contains the driving factor at that time. Samples from multiple consecutive time stamps are arranged chronologically to form a real-time monitoring sequence.

[0143] Step S820: Discretize the three-dimensional space within the digital twin of the mining area into multiple spatial cube units, and determine the target spatial cube unit to which each monitoring point belongs among the multiple spatial cube units.

[0144] The three-dimensional space of the digital twin of the mining area is divided into a regular three-dimensional grid, and the size of each cubic unit is set as needed. A coordinate system is established, and by calculating the spatial coordinates (x, y, z) of each monitoring point and the spatial range of the grid unit, the unique target cubic unit to which that point belongs is determined, denoted as the grid index (gx, gy, gz). Subsequent data assimilation and correction are all performed within this index range.

[0145] Step S830: Perform time window sliding slicing processing on the real-time monitoring sequence, and slide and cut along the time axis of the real-time monitoring sequence with a preset time window length to generate multiple time window real-time monitoring segments.

[0146] Set a fixed time window length (e.g., 6 hours of continuous monitoring data). Starting from the beginning of the real-time monitoring sequence, extract continuous data segments of length L to form a fragment; then slide the window backward by a preset step size to extract the next fragment. Repeat this process until the entire real-time monitoring sequence has been processed.

[0147] Step S840: Perform multi-dimensional feature separation processing on each time window real-time monitoring segment, and split the real-time monitoring value of heavy metal concentration and the real-time monitoring value of environmental driving factors in the time window real-time monitoring segment along the feature dimension to generate a time series segment of heavy metal concentration and a time series segment of environmental driving factors.

[0148] Take a time window segment generated in step S830, separate the heavy metal concentration data sequence and the driving factor data sequence in it, and generate two independent sub-segments: one consisting of the time series of concentration values, and the other consisting of the time series of driving factors (such as rainfall, wind speed, temperature and other parameters).

[0149] Step S850: Perform multi-scale temporal feature extraction processing on the time series segments of the environmental driving factors, and smooth the time series segments of the environmental driving factors by using long-term span moving average operation and short-term span moving average operation respectively to generate long-term environmental driving features and short-term environmental driving features.

[0150] For the time series segments of environmental driving factors obtained in step S840, a long-span moving average filter (e.g., 10 time steps) is applied to each segment to obtain long-term environmental driving characteristics. Then, a short-span moving average filter (e.g., 2 time steps) is applied to obtain short-term environmental driving characteristics. These two characteristics reflect the long-term trend and short-term fluctuations of the driving factors, respectively.

[0151] Step S860: Input the long-term environmental driving features and short-term environmental driving features into the parameter offset generation module, perform time-step splicing of the long-term environmental driving features and short-term environmental driving features, and generate a time-step offset correction sequence of the multi-medium convection diffusion coefficient in the heavy metal migration and diffusion simulation model based on the spliced ​​features.

[0152] The long-term and short-term environment-driven feature sequences from step S850 are concatenated time-step by time to form a new multidimensional feature vector sequence. This sequence is then input into a parameter offset generation module (e.g., a neural network based on a multilayer perceptron). The network outputs correction coefficients for the multi-medium convection-diffusion coefficients at the corresponding time steps, generating a complete sequence representing the magnitude of the adjustment to the current model parameters at each time step.

[0153] Step S870: Using the time-step offset correction sequence, perform time-step offset correction on the convective diffusion coefficient values ​​of the heavy metal convective diffusion coefficients in the soil medium, groundwater medium, surface water medium, and atmospheric medium in the heavy metal migration and diffusion simulation model at the target spatial cube unit, and generate the online corrected spatial distribution field of the four-layer medium convective diffusion coefficients.

[0154] Based on the correction amount of each medium convection diffusion coefficient generated in step S860, calculate the new coefficient at the target grid cell position: K_input_new=K_input_old+K_mod, and write the corrected convection diffusion coefficient back into the grid cell parameters of the corresponding medium layer for online updating.

[0155] Step S880: Using the real-time three-dimensional distribution field of heavy metals at the current time step as the initial concentration field, the spatial distribution field of the online corrected four-layer medium convection diffusion coefficient is input into the convection diffusion solution module of the heavy metal migration and diffusion simulation model to perform rolling re-simulation processing on the subsequent migration and diffusion process of heavy metals in the digital twin of the mining area, and generate an online calibrated three-dimensional evolution sequence of heavy metals.

[0156] A simulation calculation is performed using the updated convection-diffusion coefficient from step S870. Using the current concentration distribution as the initial condition and the corrected coefficient as the model parameters, the solution process in step S146 is executed to obtain a new prediction result. This prediction result reflects the impact of real-time environmental changes on diffusion.

[0157] Step S890: Extract the predicted heavy metal concentration at the target spatial cube cell in the next time step from the online calibrated three-dimensional evolution sequence of heavy metals, calculate the deviation between the predicted heavy metal concentration and the actual real-time monitoring value of heavy metal concentration obtained in the next time step, generate an online calibration deviation feedback quantity, and use the online calibration deviation feedback quantity to perform online iterative adjustment of the calculation parameters in the parameter offset generation module.

[0158] The predicted concentration of the target grid for the next time step is extracted from the heavy metal 3D evolution sequence generated in step S880 and compared with the actual monitored concentration at the same location of the real-time monitoring node to calculate the deviation value. The deviation value is fed back to the parameter offset generation module to adjust its internal parameters (such as the weights and biases of neurons), so that the output correction is gradually optimized and the deviation between prediction and measurement is reduced.

[0159] Step S8100: Update the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period based on the online calibrated three-dimensional evolution sequence of heavy metals, and generate the spatiotemporal evolution prediction results of heavy metal pollution with online data assimilation calibration.

[0160] Using the evolution sequence updated in step S880, steps S151 to S157 are re-executed to generate a predicted field and frontal trajectory that reflect the actual monitoring data, thus obtaining online assimilation results.

[0161] Step S910: Obtain the geological structure data of the target mining area, which includes the spatial location data of rock strata interfaces, the spatial distribution data of fault strikes, and the spatial distribution data of fracture development zones.

[0162] Surface coordinates of rock strata interfaces, spatial distribution data of fault lines, and distribution area data of fracture zones were obtained using geological exploration data, borehole data, and geological models. This data was then imported into the mining area's environmental geological database to ensure spatial coordinate consistency.

[0163] Step S911: Extract soil media units and groundwater media units from the soil media layer and groundwater media layer of the digital twin of the mining area, respectively. Compare the spatial location data of the rock strata interface with the spatial locations of the soil media units and groundwater media units. Find the spatial surface morphology of the rock strata interface at the soil media unit and groundwater media unit. On the spatial surface morphology of the rock strata interface, for each spatial surface point, determine the anisotropic conduction coefficient tensor used to describe the solute exchange near the interface based on the interface normal vector and the anisotropic parameters of the media on both sides of the interface.

[0164] The spatial location data of the rock strata interface are projected onto a digital twin mesh. For each mesh node intersecting the interface, a local surface is fitted using the interface shape function, and its normal vector is calculated. Based on the anisotropic parameters of the media on both sides of the interface, a 3x3 anisotropic conduction coefficient tensor is constructed, which replaces the original isotropic coefficients.

[0165] Step S912: Overlay the spatial distribution data of the fault strike with the spatial location of the groundwater medium unit, extract the groundwater medium unit that has spatial overlap with the spatial distribution data of the fault strike, mark it as the fault channel medium unit, and store the three-dimensional spatial location and interconnection relationship of all fault channel medium units as the fault channel priority flow network.

[0166] The fault line vector is projected onto the groundwater grid of the digital twin, and an influence radius Rf is set according to the width of the fault. Grid cells within this influence radius are marked as fault channel medium cells, and a connectivity graph of these cells is constructed.

[0167] Step S913: Overlay the spatial distribution data of the fracture development zone with the spatial location of the soil medium unit, extract the soil medium unit that has spatial overlap with the spatial distribution data of the fracture development zone, mark it as fracture channel medium unit, and store the spatial location and interconnection relationship of all fracture channel medium units as fracture channel priority flow network.

[0168] The polygonal region of the fracture zone is projected onto the soil layer, and the grid cells within the polygon are marked as fracture channel medium cells to establish a fracture channel preferential flow network.

[0169] Step S914: In the heavy metal migration and diffusion simulation model, the heavy metal flux exchange treatment between the soil medium and the groundwater medium is modified. At the location of the spatial curved surface morphology of the rock layer interface, the original isotropic convection diffusion coefficient in the heavy metal flux exchange treatment is replaced with the anisotropic convection diffusion tensor at each spatial curved surface point to generate a multi-media interface heavy metal coupling flux combined with the geological interface.

[0170] During the model solution process, when encountering rock strata interfaces, the isotropic exchange coefficients are no longer used. Instead, the flux is calculated using the anisotropic tensor calculated in step S911 to correct the interface flux.

[0171] Step S915: In the convection diffusion solution module of the heavy metal migration and diffusion simulation model, the relevant parameters of the groundwater medium at the fault channel medium unit are locally enhanced. The hydraulic conductivity coefficient of the fault channel medium unit is multiplied by the preset fault preferential flow enhancement factor to obtain the enhanced hydraulic conductivity coefficient. Then, the spatial distribution field of the heavy metal convection diffusion coefficient of the groundwater medium enhanced by the fault is calculated.

[0172] For the unit marked as a fault channel, its permeability coefficient is multiplied by a preset enhancement factor to increase the convection diffusion coefficient of the unit, simulating the effect of fault preferential flow rapidly transporting heavy metals.

[0173] Step S916: In the convection diffusion solution module of the heavy metal migration and diffusion simulation model, the soil medium-related parameters at the fracture channel medium unit are locally enhanced. The seepage velocity of the fracture channel medium unit is multiplied by the preset fracture preferential flow enhancement factor to obtain the enhanced seepage velocity. Then, the spatial distribution field of the heavy metal convection diffusion coefficient of the fracture-enhanced soil medium is calculated.

[0174] For fractured channel media units, the soil water infiltration velocity is increased to simulate the fracture preferential flow effect.

[0175] Step S917: Using the current time step's three-dimensional distribution field of heavy metals as the initial concentration field, input the heavy metal coupling flux of the multi-media interface combined with the geological interface, the spatial distribution field of the heavy metal convection diffusion coefficient of the fault-enhanced groundwater medium, and the spatial distribution field of the heavy metal convection diffusion coefficient of the fracture-enhanced soil medium into the convection diffusion solution module of the heavy metal migration and diffusion simulation model, simulate the rapid migration and diffusion process of heavy metals in the preferential flow channel within the digital twin of the mining area, and generate a three-dimensional evolution sequence of heavy metals combined with geological structures.

[0176] Simulations were performed using the updated coupling parameters and priority flow parameters from steps S914, S915, and S916 to generate a new evolution sequence.

[0177] Step S918: Update the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period based on the three-dimensional evolution sequence of heavy metals combined with geological structure, and generate the prediction result of the spatiotemporal evolution of heavy metal pollution combined with the influence of geological structure.

[0178] The process of steps S151 to S157 is performed using the new evolution sequence described above, and the prediction results that take into account the influence of geological structures are output.

[0179] In one exemplary embodiment, a digital twin-based simulation and prediction system for heavy metal pollution in mining areas is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 3 As shown, it includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements a digital twin-based method for simulating and predicting heavy metal pollution in mining areas. The display unit is used to form a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the casing of the digital twin-based mining heavy metal pollution simulation and prediction system, or an external keyboard, touchpad, or mouse, etc.

[0180] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A method for simulating and predicting heavy metal pollution in mining areas based on digital twins, characterized in that, The method includes: Obtain a multi-source environmental monitoring data set of the target mining area, wherein the multi-source environmental monitoring data set includes heavy metal concentration sampling values ​​of multiple monitoring points and spatial coordinate information of each monitoring point; The multi-source environmental monitoring data set is injected into a pre-constructed digital twin of the mining area. Based on the spatial coordinate information of each monitoring point, the heavy metal concentration sampling value of each monitoring point is mapped to the corresponding spatial position of the digital twin of the mining area, generating a spatially scattered set of heavy metal concentration points that are discretely distributed within the digital twin of the mining area. The three-dimensional spatial continuous field reconstruction process is performed on the spatial scatter set of heavy metal concentration to generate the three-dimensional distribution field of heavy metals inside the digital twin of the mining area. A physical process-driven heavy metal migration and diffusion simulation model is invoked to perform multi-media coupled migration and diffusion simulation of the three-dimensional distribution field of the heavy metal in the digital twin of the mining area, generating a three-dimensional evolution sequence of the heavy metal in the digital twin of the mining area in the future time period. Based on the three-dimensional evolution sequence, a spatiotemporal evolution prediction result of heavy metal pollution in the target mining area is generated. The spatiotemporal evolution prediction result of heavy metal pollution includes the predicted field of heavy metal spatial distribution and the advancing trajectory of the pollution diffusion front in each future time period.

2. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The process involves injecting the multi-source environmental monitoring data set into a pre-constructed digital twin of the mining area, mapping the heavy metal concentration sampling value of each monitoring point to the corresponding spatial location within the digital twin of the mining area based on the spatial coordinate information of each monitoring point, and generating a discretely distributed spatial scatter set of heavy metal concentration points within the digital twin of the mining area, including: Extract the medium category identifier corresponding to the heavy metal concentration sampling value of each monitoring point in the multi-source environmental monitoring data set. The medium category identifier is used to distinguish between soil monitoring media, groundwater monitoring media, surface water monitoring media and atmospheric deposition monitoring media. The multi-source environmental monitoring data set is divided into media layers according to the media category identifier, generating a soil media sampling record subset, a groundwater media sampling record subset, a surface water media sampling record subset, and an atmospheric deposition media sampling record subset. The heavy metal concentration sampling values ​​of each monitoring point in the soil medium sampling record subset are subjected to soil background correction processing. The heavy metal concentration sampling values ​​of each monitoring point are compared with the preset soil background concentration benchmark value to generate a set of heavy metal concentration increments in the soil medium. Each element in the set of heavy metal concentration increments in the soil medium corresponds to the heavy metal concentration deviation of a monitoring point relative to the soil background. The heavy metal concentration sampling values ​​of each monitoring point in the groundwater medium sampling record subset are subjected to groundwater background correction processing. The heavy metal concentration sampling values ​​of each monitoring point are compared with the preset groundwater background concentration benchmark value to generate a set of groundwater medium heavy metal concentration increments. The heavy metal concentration sampling values ​​of each monitoring point in the surface water medium sampling record subset are subjected to surface water background correction processing. The heavy metal concentration sampling values ​​of each monitoring point are compared with the preset surface water background concentration benchmark value to generate a set of surface water medium heavy metal concentration increments. The atmospheric background correction process is performed on the heavy metal concentration sampling values ​​of each monitoring point in the atmospheric deposition medium sampling record subset. The heavy metal concentration sampling values ​​of each monitoring point are then compared with the preset atmospheric background concentration benchmark value to generate a set of heavy metal concentration increments in the atmospheric deposition medium. The sets of heavy metal concentration increments in soil, groundwater, surface water, and atmospheric deposition media are mapped to the corresponding spatial locations of the soil, groundwater, surface water, and atmospheric media layers in the digital twin of the mining area, according to their respective media category identifiers and the spatial coordinate information of each monitoring point. This generates a set of spatial scatter points for the heavy metal concentration increments of the four media layers in the digital twin of the mining area, with each scatter point containing spatial coordinates, media category identifiers, and concentration increment values. The heavy metal concentration increments corresponding to the four media layers in the spatial scatter set of heavy metal concentration increments are aligned between media layers in the unified three-dimensional spatial coordinate system of the mining area digital twin. The heavy metal concentration increments of different media layers at the same spatial coordinate projection position are recorded together with their corresponding media category identifiers to generate a spatial scatter set of heavy metal concentrations that is discretely distributed in the mining area digital twin and marked with the media category.

3. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The step of performing three-dimensional spatial continuous field reconstruction processing on the spatial scatter set of heavy metal concentrations to generate a three-dimensional distribution field of heavy metals inside the digital twin of the mining area includes: The spatial scatter point set of heavy metal concentration is subjected to spatial clustering segmentation. Based on the spatial distribution density, the spatial scatter point set of heavy metal concentration is divided into multiple spatial clusters. Each spatial cluster contains a group of scatter points with spatially adjacent heavy metal concentration values. For each spatial cluster, the local concentration gradient direction is extracted from the scatter points of concentration. The rate of change of heavy metal concentration values ​​in each spatial cluster along each coordinate axis in the spatial coordinate system is calculated to generate a local concentration gradient vector corresponding to each spatial cluster. Based on the local concentration gradient vector corresponding to each spatial cluster, a concentration gradient field is constructed for each spatial cluster. Using the spatial coordinates, concentration values, and local concentration gradient vectors of the scatter points of concentration in each spatial cluster, the concentration values ​​at any spatial location within the spatial cluster are estimated by spatial interpolation methods to generate a local concentration gradient field corresponding to each spatial cluster. The concentration field of adjacent spatial clusters is smoothed in the boundary region. The concentration difference of the local concentration gradient field of each adjacent spatial cluster in the boundary region is extracted. The concentration difference in the boundary region is gradually eliminated by a weighted fusion method based on spatial distance, and a smooth concentration transition band between adjacent spatial clusters is generated. The local concentration gradient fields corresponding to all spatial clusters and all smooth concentration transition zones are spatially spliced ​​and integrated to generate an initial three-dimensional distribution field of heavy metals covering the internal space of the digital twin of the entire mining area. The initial three-dimensional distribution field of heavy metals is then subjected to non-negative constraint correction processing, and the concentration values ​​of spatial locations with concentration values ​​less than 0 in the initial three-dimensional distribution field of heavy metals are set to 0 to generate a non-negative three-dimensional distribution field of heavy metals. The physical consistency correction process is performed on the non-heavy metal three-dimensional distribution field. The concentration values ​​of spatial locations in the non-heavy metal three-dimensional distribution field that exceed the preset theoretical upper limit of heavy metal concentration in the mining area are adjusted to the theoretical upper limit of heavy metal concentration in the mining area, thereby generating a heavy metal three-dimensional distribution field inside the digital twin of the mining area.

4. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The process involves invoking a physics-driven heavy metal migration and diffusion simulation model to perform multi-media coupled migration and diffusion simulations of the three-dimensional distribution field of heavy metals within the digital twin of the mining area, generating a three-dimensional evolution sequence of heavy metals in the digital twin of the mining area over a future time period, including: Spatial distribution data of soil moisture content and spatial distribution data of soil porosity are extracted from the soil medium layer of the digital twin of the mining area. The spatial distribution data of soil moisture content and spatial distribution data of soil porosity are input into the soil medium parameter module of the heavy metal migration and diffusion simulation model to generate a spatial distribution field of heavy metal convection and diffusion coefficient in the soil medium. The groundwater velocity vector field and aquifer permeability coefficient spatial distribution data are extracted from the groundwater medium layer of the digital twin of the mining area. The groundwater velocity vector field and aquifer permeability coefficient spatial distribution data are input into the groundwater medium parameter module of the heavy metal migration and diffusion simulation model to generate the spatial distribution field of heavy metal convection diffusion coefficient in the groundwater medium. The surface runoff velocity vector field and the spatial distribution data of heavy metal adsorption coefficient in riverbed sediment are extracted from the surface water medium layer of the digital twin of the mining area. The surface runoff velocity vector field and the spatial distribution data of heavy metal adsorption coefficient in riverbed sediment are input into the surface water medium parameter module of the heavy metal migration and diffusion simulation model to generate the spatial distribution field of heavy metal convection diffusion coefficient in the surface water medium. The near-surface wind speed vector field and the spatial distribution data of the atmospheric boundary layer turbulent diffusion coefficient are extracted from the atmospheric medium layer of the digital twin of the mining area. The near-surface wind speed vector field and the spatial distribution data of the atmospheric boundary layer turbulent diffusion coefficient are input into the atmospheric medium parameter module of the heavy metal migration and diffusion simulation model to generate the spatial distribution field of the heavy metal convection diffusion coefficient in the atmospheric medium. The medium coupling module of the heavy metal migration and diffusion simulation model is invoked to perform heavy metal flux exchange processing between soil and groundwater at the interface between soil and groundwater, between soil and surface water at the interface between soil and atmospheric, between soil and groundwater at the interface between surface water and groundwater, generating a multi-media interface heavy metal coupling flux set. Using the three-dimensional distribution field of heavy metals as the initial concentration field, the spatial distribution fields of the heavy metal convection diffusion coefficients in the soil medium, the groundwater medium, the surface water medium, the atmospheric medium, and the multi-media interface heavy metal coupled flux set are input into the convection diffusion solution module of the heavy metal migration and diffusion simulation model. The multi-media convection diffusion equations are coupled and iteratively solved within the digital twin of the mining area according to a preset time step. The spatial distribution of heavy metal concentrations in the soil medium layer, groundwater medium layer, surface water medium layer, and atmospheric medium layer is updated within each time step, generating a three-dimensional evolution sequence of heavy metals in the digital twin of the mining area within the future time period.

5. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The method of generating spatiotemporal evolution prediction results of heavy metal pollution in the target mining area based on the three-dimensional evolution sequence, wherein the spatiotemporal evolution prediction results of heavy metal pollution include the predicted spatial distribution field of heavy metals in future time periods and the advancing trajectory of the pollution diffusion front, including: The spatial distribution of heavy metal concentration in soil, groundwater, surface water, and atmosphere at each time step is extracted from the three-dimensional evolution sequence. The spatial distribution of heavy metal concentration in the four media at the same time step is then spatially superimposed in the unified three-dimensional spatial coordinate system of the digital twin of the mining area to generate a spatial distribution superposition field of heavy metal at each time step. The pollution level is classified for the superimposed field of heavy metal spatial distribution corresponding to each time step. Spatial areas with heavy metal concentration values ​​exceeding the preset heavy metal pollution risk threshold are marked as pollution risk areas, and spatial areas with heavy metal concentration values ​​exceeding the preset heavy metal pollution warning threshold are marked as pollution warning areas. A heavy metal pollution zoning label field corresponding to each time step is generated. The heavy metal pollution zoning label fields corresponding to all time steps are arranged in chronological order to generate a prediction field of heavy metal spatial distribution for future time periods. Extract the heavy metal pollution zone labeling field of two adjacent time steps from the predicted heavy metal spatial distribution field of each future time period. Take the area boundary where the spatial location of the same pollution zone labeling category changes in two adjacent time steps as pollution diffusion front segment. Perform spatial continuity splicing on all pollution diffusion front segments between two adjacent time steps to generate pollution diffusion front between adjacent time steps. The pollution diffusion fronts between all adjacent time steps are arranged in time series. The spatial displacement vectors of each spatial location point on each pollution diffusion front between adjacent time steps are extracted. The spatial displacement vectors of each spatial location point on each pollution diffusion front are time normalized according to the time interval between the two time steps corresponding to the pollution diffusion front, and the diffusion propulsion velocity vectors of each spatial location point on each pollution diffusion front are generated. The diffusion velocity vectors of all spatial locations on the entire pollution diffusion front are spatially interpolated in the three-dimensional spatial coordinate system of the digital twin of the mining area to generate a continuous diffusion velocity vector field covering the entire spatial region where the pollution diffusion front is located. In the continuous vector field of diffusion velocity, starting from the preset initial spatial location of the pollution source, three-dimensional streamline tracking is performed along the vector direction of each spatial location in the continuous vector field of diffusion velocity to generate multiple pollution diffusion streamlines. Each pollution diffusion streamline represents the movement path of heavy metal pollution in space from the initial spatial location of the pollution source along the diffusion direction. All pollution diffusion streamlines are integrated in the three-dimensional space of the digital twin of the mining area in terms of time dimension. Along the path direction of each pollution diffusion streamline, the product of the diffusion propagation velocity vector of each spatial location point and the time integration step is accumulated segment by segment with a preset time integration step to obtain the spatial coordinates corresponding to each time integration node on each pollution diffusion streamline. The spatial envelope surface is constructed by the spatial coordinates corresponding to the same time integration node on all pollution diffusion streamlines to generate the propagation trajectory of the pollution diffusion front.

6. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The method further includes: Obtain the time-series forecast data of rainfall in the target mining area and the surface runoff generation and runoff model of the mining area. Input the time-series forecast data of rainfall in the target mining area into the surface runoff generation and runoff model of the mining area to generate the time-series forecast curve of surface runoff flow and the time-series evolution data of surface runoff inundation range of the target mining area in the future time period. The time-series prediction curve of surface runoff and the time-series evolution data of surface runoff inundation range are input into the surface water medium layer of the digital twin of the mining area. Combined with the digital elevation model, the spatial distribution data of surface runoff water depth and the surface runoff velocity vector field at each time step are calculated, and the corresponding data in the surface water medium layer are updated. Based on the updated surface runoff velocity vector field and surface runoff water depth spatial distribution data, the spatial distribution field of heavy metal convection diffusion coefficient in the surface water medium is recalculated in the heavy metal migration and diffusion simulation model. Using the current time step of the three-dimensional distribution field of heavy metals as the initial concentration field, the recalculated spatial distribution field of the heavy metal convection diffusion coefficient in the surface water medium and the updated surface runoff velocity vector field are input into the convection diffusion solution module of the heavy metal migration and diffusion simulation model to re-simulate the migration and diffusion process of heavy metals in the digital twin of the mining area under rainfall conditions, and generate a corrected three-dimensional evolution sequence of heavy metals under rainfall conditions. Based on the modified three-dimensional evolution sequence of heavy metals under the rainfall conditions, the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period are updated to generate the prediction results of the spatiotemporal evolution of heavy metal pollution under the rainfall scenario.

7. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The method further includes: Obtain heavy metal adsorption characteristic data of the target mining area soil, which includes spatial distribution data of soil cation exchange capacity, spatial distribution data of soil organic matter content, and spatial distribution data of soil pH. The heavy metal adsorption characteristics data of the mining area soil are input into the soil medium parameter module of the heavy metal migration and diffusion simulation model to perform adsorption correction processing on the heavy metal migration process in the soil medium. Based on the spatial distribution data of soil cation exchange capacity and soil organic matter content, a spatial distribution field of adsorption inhibition factors for heavy metals in the soil medium is generated; based on the spatial distribution data of soil pH, a spatial distribution field of reaction rate parameters for heavy metals in the soil medium is generated; the spatial distribution field of adsorption inhibition factors is used to correct the convective transport velocity of heavy metals in the soil medium, and the spatial distribution field of reaction rate parameters is used to characterize the adsorption reaction terms in the migration and diffusion simulation, thereby realizing the adsorption correction of the heavy metal migration and diffusion process in the soil medium. Using the current time step of the three-dimensional distribution field of heavy metals as the initial concentration field, the spatial distribution field of the convective diffusion coefficient of heavy metals in the soil medium after adsorption correction is replaced with the original spatial distribution field of the convective diffusion coefficient of heavy metals in the soil medium. The convective diffusion solution module of the heavy metal migration and diffusion simulation model is input to re-simulate the migration and diffusion process of heavy metals in the soil medium in the digital twin of the mining area, and generate the three-dimensional evolution sequence of heavy metals after soil adsorption correction. Based on the three-dimensional evolution sequence of heavy metals corrected by soil adsorption, the predicted spatial distribution field of heavy metals and the trajectory of the pollution diffusion front for each future time period are updated, generating a prediction result of the spatiotemporal evolution of heavy metal pollution combined with soil adsorption.

8. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The method further includes: Obtain the time-series monitoring data of the groundwater level in the target mining area and the groundwater recharge and discharge parameters of the mining area. Based on the time-series monitoring data of the groundwater level in the mining area and the groundwater recharge and discharge parameters of the mining area, determine the groundwater level fluctuation curve and groundwater flow field change trend of the target mining area in the future time period. The groundwater level fluctuation curve and groundwater flow field change trend are input into the groundwater medium layer of the digital twin of the mining area. Combined with the spatial distribution data of aquifer permeability coefficient, the groundwater velocity vector field at each time step is calculated through the groundwater flow model. The spatial distribution data of aquifer head and the groundwater velocity vector field are updated. Based on the updated groundwater velocity vector field and aquifer head spatial distribution data, the spatial distribution field of heavy metal convection diffusion coefficient in the groundwater medium is recalculated in the heavy metal migration and diffusion simulation model. Using the current time step of the three-dimensional distribution field of heavy metals as the initial concentration field, the recalculated spatial distribution field of the heavy metal convection diffusion coefficient in the groundwater medium and the updated groundwater velocity vector field are input into the convection diffusion solution module of the heavy metal migration and diffusion simulation model to re-simulate the migration and diffusion process of heavy metals in the groundwater medium within the digital twin of the mining area, and generate a three-dimensional evolution sequence of heavy metals after groundwater level change correction. Based on the three-dimensional evolution sequence of heavy metals corrected by the groundwater level change, the predicted spatial distribution field of heavy metals and the trajectory of the pollution diffusion front for each future time period are updated, generating a prediction result of the spatiotemporal evolution of heavy metal pollution combined with groundwater level changes.

9. The method for simulating and predicting heavy metal pollution in mining areas based on digital twins according to claim 1, characterized in that, The method further includes: Acquire the time-series prediction data of atmospheric wind field and the topographic roughness distribution data of the target mining area, input the time-series prediction data of atmospheric wind field and the topographic roughness distribution data of the mining area into the atmospheric boundary layer flow field calculation model, and generate the corrected wind speed vector field and spatial distribution data of turbulence diffusion coefficient near the surface of the target mining area in the future time period; The modified wind speed vector field and the spatial distribution data of the turbulent diffusion coefficient are loaded into the atmospheric medium layer of the digital twin of the mining area, and the near-surface wind speed vector field and the spatial distribution data of the atmospheric boundary layer turbulent diffusion coefficient at each time step in the atmospheric medium layer are updated. Based on the updated near-surface wind speed vector field and the spatial distribution data of atmospheric boundary layer turbulent diffusion coefficient, the spatial distribution field of heavy metal convective diffusion coefficient in the atmospheric medium is recalculated in the heavy metal migration and diffusion simulation model. Using the current time step of the three-dimensional distribution field of heavy metals as the initial concentration field, the recalculated spatial distribution field of the heavy metal convection diffusion coefficient in the atmospheric medium and the updated near-surface wind speed vector field are input into the convection diffusion solution module of the heavy metal migration and diffusion simulation model to re-simulate the migration and diffusion process of heavy metals in the atmospheric medium within the digital twin of the mining area, and generate a three-dimensional evolution sequence of heavy metals after atmospheric wind field correction. Based on the three-dimensional evolution sequence of heavy metals corrected by the atmospheric wind field, the predicted spatial distribution field of heavy metals and the advancing trajectory of the pollution diffusion front for each future time period are updated, generating a prediction result of the spatiotemporal evolution of heavy metal pollution that combines changes in the atmospheric wind field.

10. A simulation and prediction system for heavy metal pollution in mining areas based on digital twins, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the digital twin-based method for simulating and predicting heavy metal pollution in mining areas as described in any one of claims 1 to 9 by executing the machine-executable instructions.

Images