A method and system for digital simulation of soil evolution processes

Abstract

The application relates to the technical field of soil water-heat-salt evolution simulation, and discloses a soil evolution process digital simulation method and system.The method comprises the following steps: acquiring soil monitoring data of a monitored soil area, wherein the soil monitoring data comprises salt observation data and water-heat-pore observation data; performing space-time gridding on the soil monitoring data to obtain three-dimensional soil grid data comprising surface layer grid values and deep layer grid values; generating surface layer constraints according to the surface layer grid values, generating deep layer constraints according to the deep layer grid values, and performing continuous field estimation to obtain an initial water-heat-salt field; identifying a salt abnormal area; performing water-salt coupling state updating according to a last period salt state, a current salt observation set and a water-heat-pore state to obtain a local salt state; writing the local salt state into the initial water-heat-salt field to obtain a fused water-heat-salt field, and performing continuous correction to obtain a target soil evolution simulation field.The method can recursively update and continuously correct the salt abnormal area.

Description

Technical Field

[0001] This invention relates to the field of soil hydrothermal-salt evolution simulation technology, and in particular to a digital simulation method and system for soil evolution processes. Background Technology

[0002] Currently, soil internal state parameters such as temperature, moisture content, salinity, and pore structure continuously change over time and space, and these parameters are interconnected. Continuous simulation of soil evolution processes helps analyze the changing trends of water migration, heat transfer, and salt accumulation in different soil layers. With the development of remote sensing, soil sensing, and spatial interpolation technologies, soil condition monitoring has gradually shifted from single-point sampling to multi-source data collaborative analysis. Related methods typically require mapping surface bedding data and deep point data into a unified space for model fusion.

[0003] In a current technology, surface remote sensing data and deep sensor data are typically acquired first, and then a soil state distribution map is generated through spatial interpolation or statistical regression. This method can supplement data gaps in areas without sensor deployment to some extent, but it is mostly static reconstruction, focusing more on the soil state distribution at a single moment. For areas with abnormally high salinity concentrations, existing methods often rely directly on the current interpolation results, lacking a joint update of the relationship between the previous period's salinity state, current salinity observation data, and hydrothermal porosity. This makes it difficult to continuously correct for local salinity anomalies as soil evolves. Furthermore, surface remote sensing data has a large spatial coverage but mainly reflects surface conditions; deep sensor data can reflect local conditions at different depths, but the observation points are discretely distributed. Simply stitching together or interpolating the two types of data can easily weaken the coupling relationship between deep salinity changes and moisture, heat, and porosity. When the salinity anomaly area is located in a deep or localized area, it is difficult to correct the local salinity state in a timely manner by simply relying on the global continuous field reconstruction results, and it is also difficult to ensure that the updated local area maintains spatial continuity with the surrounding soil state field.

[0004] Existing technologies have the problem that it is difficult to couple and continuously correct local salinity anomaly areas in soil hydrothermal-salt evolution simulation fields. Summary of the Invention

[0005] This invention provides a digital simulation method and system for soil evolution processes to solve the problem in existing soil hydrothermal-salt evolution simulation fields that are difficult to couple, update, and continuously correct local salinity anomaly areas.

[0006] Firstly, in order to solve the above-mentioned technical problems, the present invention provides a digital simulation method for soil evolution processes, comprising: Soil monitoring data, including salinity observation data and hydrothermal porosity observation data, is acquired from the monitored soil area. The soil monitoring data is then spatiotemporally gridded to obtain three-dimensional soil grid data, including surface grid values ​​and deep grid values. Surface constraints are generated based on the surface grid values, deep constraints are generated based on the deep grid values, and continuous field estimation is performed on the three-dimensional soil grid nodes based on the surface constraints and the deep constraints to obtain the initial hydrothermal salinity field. From the initial hydrothermal salt field, identify spatially continuous grid nodes where the salt concentration exceeds the preset soil salinity level limit to obtain salinity anomaly zones. Obtain the previous cycle salinity state corresponding to the salinity anomaly area, extract the current salinity observation set from the soil monitoring data, extract the hydrothermal porosity state from the initial hydrothermal salinity field, and update the water-salt coupling state based on the previous cycle salinity state, the current salinity observation set, and the hydrothermal porosity state to obtain the local salinity state. The local salinity state is written into the salinity field of the grid node corresponding to the salinity anomaly zone in the initial hydrothermal salinity field to obtain the fused hydrothermal salinity field. The fused hydrothermal salinity field is then continuously corrected according to the boundary of the salinity anomaly zone to obtain the target soil evolution simulation field.

[0007] Secondly, the present invention provides a digital simulation system for soil evolution processes, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method described above.

[0008] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention acquires soil monitoring data from a monitored soil area and performs spatiotemporal gridding on the soil monitoring data to obtain three-dimensional soil grid data containing surface grid values ​​and deep grid values; further, it generates surface constraints based on the surface grid values ​​and deep constraints based on the deep grid values, and performs continuous field estimation on the three-dimensional soil grid nodes based on the surface constraints and the deep constraints to obtain an initial hydrothermal salinity field. This process converts surface bedding observation information and deep point observation information into a unified three-dimensional grid representation, enabling soil monitoring data from different spatial scales and depths to participate in the same continuous field reconstruction process, thereby reducing the problem of incomplete expression of the soil's internal state caused by relying on only a single observation source.

[0009] (2) This invention identifies spatially continuous grid nodes in the initial hydrothermal salt field where the salt concentration exceeds a preset soil salinity level limit, thus obtaining a salinity anomaly zone. It then acquires the previous cycle salinity state corresponding to the salinity anomaly zone, extracts the current salinity observation set from the soil monitoring data, extracts the hydrothermal porosity state from the initial hydrothermal salt field, and updates the water-salt coupling state based on the previous cycle salinity state, the current salinity observation set, and the hydrothermal porosity state to obtain the local salinity state. This process does not rely solely on the current salinity interpolation result for static judgment; instead, it incorporates historical salinity states, current salinity observations, and hydrothermal porosity conditions into the local state update, enabling the state correction of the salinity anomaly zone to reflect the temporal changes and water-salt coupling relationship during soil evolution.

[0010] (3) This invention obtains a fused hydrothermal salt field by writing the local salinity state into the salinity field of the grid node corresponding to the salinity anomaly area in the initial hydrothermal salt field, and then continuously corrects the fused hydrothermal salt field according to the boundary of the salinity anomaly area to obtain the target soil evolution simulation field. This process writes the local update results back to the global three-dimensional hydrothermal salt field and corrects the boundary state between the local replacement area and the surrounding area, so that the updated salinity anomaly area maintains a spatial transition relationship with the surrounding soil state field, thereby reducing the risk of boundary abrupt changes after local state update. Attached Figure Description

[0011] Figure 1 This is a schematic diagram of the digital simulation method for soil evolution provided in the first embodiment of the present invention. Detailed Implementation

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

[0013] Reference Figure 1 The first embodiment of the present invention provides a digital simulation method for soil evolution processes, including the following steps: S1. Acquire soil monitoring data in the monitored soil area, including salinity observation data and hydrothermal porosity observation data. Perform spatiotemporal gridding on the soil monitoring data to obtain three-dimensional soil grid data containing surface grid values ​​and deep grid values. S2, generate surface constraints based on the surface grid values, generate deep constraints based on the deep grid values, and perform continuous field estimation on the three-dimensional soil grid nodes based on the surface constraints and the deep constraints to obtain the initial hydrothermal saline field; S3, Identify spatially continuous grid nodes in the initial hydrothermal salt field where the salt concentration exceeds the preset soil salinity level limit to obtain salinity anomaly zones; S4, obtain the previous cycle salinity state corresponding to the salinity anomaly area, extract the current salinity observation set from the soil monitoring data, extract the hydrothermal pore state from the initial hydrothermal salt field, update the water-salt coupling state according to the previous cycle salinity state, the current salinity observation set and the hydrothermal pore state, and obtain the local salinity state; S5, the local salinity state is written into the salinity field of the grid node corresponding to the salinity anomaly area in the initial hydrothermal salinity field to obtain the fused hydrothermal salinity field. The fused hydrothermal salinity field is continuously corrected according to the boundary of the salinity anomaly area to obtain the target soil evolution simulation field.

[0014] In step S1, soil monitoring data including salinity observation data and hydrothermal porosity observation data of the monitored soil area are acquired, and the soil monitoring data is spatiotemporally gridded to obtain three-dimensional soil grid data including surface grid values ​​and deep grid values.

[0015] The soil monitoring data includes surface remote sensing temperature data, surface remote sensing moisture content data, deep sensing temperature data, deep sensing moisture content data, deep sensing salinity data, and deep sensing porosity data.

[0016] Specifically, the soil monitoring data is spatiotemporally gridded to obtain three-dimensional soil grid data containing surface grid values ​​and deep grid values, including: Extract the sampling time, spatial coordinates, and observation values ​​from the soil monitoring data; The soil monitoring data is resampled according to the sampling time to obtain synchronous monitoring data; The synchronous monitoring data is mapped to a preset three-dimensional soil grid based on the spatial coordinates to obtain three-dimensional soil grid data containing surface grid values ​​and deep grid values.

[0017] In one implementation, surface remote sensing temperature data is acquired using a thermal infrared imaging sensor, surface remote sensing moisture content data is acquired using a near-infrared imaging sensor, deep sensing temperature and moisture content data are acquired using soil temperature and humidity sensors buried at different depths, deep sensing salinity data is acquired using a soil conductivity sensor, and deep sensing porosity data is acquired using a soil porosity sensor. The sampling period for surface remote sensing data is determined by the time difference between adjacent sampling in the output records of the thermal infrared and near-infrared imaging sensors, and the sampling period for deep sensing data is determined by the time difference between adjacent sampling in the output records of the soil temperature and humidity, soil conductivity, and soil porosity sensors. For example, the sampling period for surface remote sensing data is 30 minutes, and the sampling period for deep sensing data is 10 minutes.

[0018] It should be noted that each soil monitoring data entry includes a data source identifier, sampling time, spatial coordinates, depth identifier, observation field name, and observation value. Observation field names include temperature, moisture content, salinity, and porosity. The salinity observation data is generated from deep-sensor salinity data, and the hydrothermal porosity observation data is generated from surface remote-sensor temperature data, surface remote-sensor moisture content data, deep-sensor temperature data, deep-sensor moisture content data, and deep-sensor porosity data. If multiple observation values ​​exist under the same data source identifier, sampling time, spatial coordinates, and observation field name, the arithmetic mean of these multiple observation values ​​is calculated, and the result is used as the observation value for the corresponding data record.

[0019] It should be noted that time resampling adopts a reference time axis alignment method. The reference time axis consists of sampling times where valid observation values ​​exist simultaneously in surface remote sensing temperature data and surface remote sensing water content data. For deep sensing temperature data, deep sensing water content data, deep sensing salinity data, and deep sensing porosity data, if the target sampling time on the reference time axis is the same as a certain sampling time in the deep sensing data, the observation value corresponding to that sampling time is read; if the target sampling time is between two adjacent sampling times of the same sensing point and the same observation field name, linear interpolation is performed according to the time interval ratio between the two adjacent sampling times and the target sampling time to obtain the observation value corresponding to the target sampling time; if the target sampling time is not between two adjacent sampling times, the corresponding observation field is marked as null. The above processing results constitute synchronous monitoring data.

[0020] It is worth noting that the preset 3D soil mesh is generated during the system initialization phase. The coverage of the horizontal grid nodes is determined by the minimum circumscribed boundary of the spatial coordinates in the soil monitoring data. The horizontal grid spacing is taken as the median of the pixel spatial spacing in the surface remote sensing temperature data and the surface remote sensing moisture content data. The depth grid nodes include surface depth nodes and deep depth nodes. The interlayer spacing of the deep depth nodes is taken as the median of the differences between adjacent depth markers in the deep sensing temperature data, deep sensing moisture content data, deep sensing salinity data, and deep sensing porosity data.

[0021] In one implementation, surface remote sensing temperature and moisture content data are mapped to surface depth nodes in a 3D soil mesh. If the spatial coordinates of the surface cell coincide with the coordinates of the surface depth node, the corresponding observed value is directly written into the temperature or moisture content field of the surface depth node. If the spatial coordinates of the surface cell do not coincide with the coordinates of the surface depth node, four adjacent cells surrounding the surface depth node are selected, and the temperature or moisture content field of the surface depth node is calculated using bilinear interpolation.

[0022] It should be noted that when mapping deep-sensed temperature, moisture, salinity, and porosity data to a 3D soil mesh, the nearest 3D soil mesh node is located based on the spatial coordinates and depth identifier of the sensing point. The corresponding observation value is then written to the corresponding field of that 3D soil mesh node. If a 3D soil mesh node receives multiple observation values ​​for the same observation field name at the same sampling time, the multiple observation values ​​are arithmetically averaged. If a 3D soil mesh node does not receive any observation values ​​for a certain observation field name at the same sampling time, that field remains empty.

[0023] For example, if the median pixel spacing of surface remote sensing temperature data and surface remote sensing moisture content data is 0.5 meters, and the median difference between adjacent depth markers of deep sensing data is 0.2 meters, then the horizontal grid spacing of the three-dimensional soil layer grid is set to 0.5 meters, and the depth grid spacing is set to 0.2 meters. If a certain deep sensing salinity data is recorded as 3.0 g / kg at the 20th minute and 4.0 g / kg at the 40th minute, then the corresponding salinity observation value at the 30th minute is 3.5 g / kg.

[0024] It is worth noting that the 3D soil grid data is organized according to sampling time and 3D soil layer grid nodes. Each 3D soil layer grid node includes a grid node identifier, sampling time, spatial coordinates, depth identifier, temperature field, moisture content field, salinity field, porosity field, and data source identifier. Grid node data located at surface depth nodes and containing surface remote sensing temperature data or surface remote sensing moisture content data constitutes the surface grid value; grid node data located at deep depth nodes and containing deep sensing temperature data, deep sensing moisture content data, deep sensing salinity data, or deep sensing porosity data constitutes the deep grid value. Fields with empty values ​​in the 3D soil grid data are not filled in step S1.

[0025] In step S2, surface constraints are generated based on the surface grid values, deep constraints are generated based on the deep grid values, and continuous field estimation is performed on the three-dimensional soil grid nodes based on the surface constraints and the deep constraints to obtain the initial hydrothermal saline field.

[0026] The process of generating surface constraints based on the surface mesh values ​​and generating deep constraints based on the deep mesh values ​​includes: Extract the surface mesh values ​​and corresponding mesh node coordinates, and associate the surface mesh values ​​and corresponding mesh node coordinates with surface boundary data to obtain surface constraints; The deep mesh values ​​and corresponding mesh node coordinates are extracted, and the deep mesh values ​​and corresponding mesh node coordinates are associated with deep correction data to obtain deep constraints.

[0027] In one implementation, the surface grid values ​​are the temperature and moisture content fields corresponding to the surface depth nodes in the 3D soil grid data. The surface grid values ​​and corresponding grid node coordinates are extracted, and the observed values, grid node coordinates, sampling time, and observation field names are combined to form surface boundary data, thus obtaining surface constraints. The deep grid values ​​are the temperature, moisture content, salinity, and porosity fields corresponding to the deep depth nodes in the 3D soil grid data. The deep grid values ​​and corresponding grid node coordinates are extracted, and the observed values, grid node coordinates, sampling time, and observation field names are combined to form deep correction data, thus obtaining deep constraints.

[0028] It should be noted that surface and deep constraints are used as observational data with defined three-dimensional coordinates in continuous field estimation. Temperature, water content, salinity, and porosity fields under the same grid node are stored independently, without numerical mixing between different physical quantities.

[0029] The initial hydrothermal-salinity field is obtained by performing continuous field estimation on the three-dimensional soil grid nodes based on the surface constraints and the deep constraints, including: Based on the surface constraints and the deep constraints, determine the known grid nodes with observed values; The grid nodes in the three-dimensional soil grid that lack observation values ​​are identified as grid nodes to be estimated. Based on the known grid nodes, perform three-dimensional spatial interpolation on the grid nodes to be estimated to obtain continuous grid values; An initial hydrothermal salt field is generated based on the continuous grid values.

[0030] The process of interpolating the known grid nodes to be estimated in three-dimensional space to obtain continuous grid values ​​includes: The observation values ​​and three-dimensional coordinates corresponding to the target fields are extracted from the known grid nodes. The target fields include temperature, water content, salinity, and porosity. Construct a three-dimensional semivariance function based on the observed values ​​and the three-dimensional coordinates; The Kriging weights corresponding to the grid nodes to be estimated are solved based on the three-dimensional semivariance function. The observations are weighted and summed according to the Kriging weights to obtain the continuous grid values ​​of the grid nodes to be estimated.

[0031] In one implementation, known grid nodes with observed values ​​are determined based on the surface constraints and the deep constraints. Specifically, each grid node in the 3D soil grid data is traversed. If a grid node has an observed value in any target field at the current sampling time, then the grid node is determined as a known grid node for that target field. If a grid node lacks an observed value for a target field, then the grid node is determined as a grid node to be estimated for that target field.

[0032] It should be noted that three-dimensional spatial interpolation is performed separately for the temperature, moisture content, salinity, and porosity fields. For any target field, known grid nodes with observed values ​​for that target field at the current sampling time are extracted, and grid nodes lacking observed values ​​for that target field at the current sampling time are extracted. Based on the three-dimensional coordinates and observed values ​​of the known grid nodes, the continuous grid values ​​of the grid nodes to be estimated under that target field are estimated.

[0033] It is worth noting that the 3D spatial interpolation uses ordinary kriging interpolation. The inputs to ordinary kriging interpolation are the 3D coordinates of known grid nodes, the observed values ​​of known grid nodes, and the 3D coordinates of the grid node to be estimated. The 3D coordinates include horizontal, vertical, and depth coordinates. Ordinary kriging interpolation does not project surface grid values ​​to deep grid values; instead, it estimates the field values ​​of the grid node to be estimated based on spatial correlation within the same 3D coordinate system.

[0034] In one implementation, when constructing the three-dimensional semivariogram function, the three-dimensional spatial distance between any two known grid nodes is first calculated, and then the distances are grouped in ascending order of distance. The group width is taken as the larger value between the horizontal grid spacing and the depth grid spacing in the three-dimensional soil layer grid. For each distance group, half the square of the difference between the observed values ​​of the known grid nodes within that group is calculated, and this result is used as the empirical semivariogram value for that group. The three-dimensional semivariogram function is then fitted based on the median distance and the empirical semivariogram value of each distance group.

[0035] It should be noted that the three-dimensional semivariance function adopts a spherical semivariance function. The spherical semivariance function includes three parameters: nugget value, sill value, and range. The nugget value is the empirical semivariance value corresponding to the minimum distance group; the sill value is the median of all empirical semivariance values; and the range is the median distance of the distance group from which the empirical semivariance value first reaches the sill value. If the empirical semivariance value does not reach the sill value, the range is half of the maximum spatial distance between known grid nodes in the three-dimensional soil mesh. All the above parameters are calculated from the known grid nodes and observed values ​​at the current sampling time.

[0036] It should be noted that the Kriging weights corresponding to the grid nodes to be estimated are obtained by solving the three-dimensional semi-variance function. Specifically, a system of Kriging linear equations is constructed based on the three-dimensional semi-variance function values ​​between known grid nodes, and the right-hand side of the equations is constructed based on the three-dimensional semi-variance function values ​​between the grid node to be estimated and each known grid node. Solving the system of Kriging linear equations yields the Kriging weights of each known grid node relative to the grid node to be estimated. Subsequently, the observed values ​​of the corresponding known grid nodes are weighted and summed according to the Kriging weights to obtain the continuous grid values ​​of the grid node to be estimated in the target field.

[0037] It is worth noting that if the number of known grid nodes for a certain target field at the current sampling time is less than 4, ordinary Kriging interpolation is not performed; instead, inverse distance weighted interpolation is used. The distance exponent for inverse distance weighted interpolation is set to 2. Known grid nodes that are closest to the grid node to be estimated and whose number does not exceed 6 are selected. The distances are normalized based on the inverse square of the 3D spatial distance from each known grid node to the grid node to be estimated, yielding corresponding interpolation weights. The observed values ​​are then weighted and summed according to these interpolation weights to obtain continuous grid values. The distance exponent of 2 ensures that the weight of nearby known grid nodes is greater than that of distant known grid nodes. The use of inverse distance weighted interpolation when the number of known grid nodes is less than 4 is because 3D spatial interpolation requires at least several known grid nodes in different spatial locations to participate in spatial relationship estimation. The selection of no more than 6 known grid nodes corresponds to the local neighborhood range of the six adjacent directions in the 3D soil layer grid.

[0038] In one implementation, an initial hydrothermal saline field is generated based on the continuous grid values. Specifically, the temperature, moisture content, salinity, and porosity fields of each three-dimensional soil grid node at each sampling time are aggregated to form a three-dimensional state field for the corresponding sampling time. The initial hydrothermal saline field includes three-dimensional state fields for multiple sampling times, and each three-dimensional state field includes a grid node identifier, three-dimensional coordinates, a temperature field, a moisture content field, a salinity field, and a porosity field.

[0039] For example, at a certain sampling time, a certain 3D soil grid node, after continuous field estimation, yields a temperature field of 22.8 degrees Celsius, a moisture content field of 19.2%, a salinity field of 3.6 g / kg, and a porosity field of 40%. The grid node identifier, 3D coordinates, and four fields of the 3D soil grid node are written into the same 3D state field record. The same writing process is performed on all 3D soil grid nodes at that sampling time to obtain the 3D state field corresponding to that sampling time.

[0040] In step S3, spatially continuous grid nodes with salt concentrations exceeding preset soil salinity level limits are identified from the initial hydrothermal salt field to obtain salinity anomaly zones, including: The salt concentration of each grid node is extracted from the initial hydrothermal salt field; Grid nodes whose salt concentration exceeds the preset soil salinity level limit are identified as abnormal nodes; The abnormal nodes are connected based on the grid adjacency relationship to obtain the salt anomaly region.

[0041] In one implementation, the initial hydrothermal salinity field is a set of three-dimensional state fields arranged according to sampling time. Each three-dimensional state field includes multiple three-dimensional soil grid nodes, and each three-dimensional soil grid node includes a grid node identifier, three-dimensional coordinates, a temperature field, a moisture content field, a salinity field, and a porosity field. Step S3 is executed within the three-dimensional state field corresponding to a single sampling time, and the same processing is performed on the three-dimensional state fields corresponding to different sampling times according to the order of sampling time.

[0042] It should be noted that the salinity concentration refers to the numerical value corresponding to the salinity field of each grid node in the initial hydrothermal salt field. Specifically, the process involves traversing all three-dimensional soil grid nodes at the current sampling time, reading the salinity field of each three-dimensional soil grid node, and determining the value in the salinity field as the salinity concentration of the corresponding grid node. If the salinity field of a certain grid node is empty, that grid node is skipped, and it is marked as a node not involved in the judgment in the abnormal node identification record.

[0043] It should be noted that the preset soil salinity level limit is determined by the historical salinity field during the system initialization phase. Specifically, the salinity field is extracted from the historical three-dimensional state field of the same monitored soil area, and null values ​​and salinity values ​​less than 0 g / kg are deleted to obtain historical salinity samples. These historical salinity samples are then arranged in ascending order of value, and the 25th and 75th percentiles are extracted. The difference between the 75th and 25th percentiles is calculated to obtain the interquartile range. The 75th percentile is then added to the interquartile range by 1.5 times to obtain the preset soil salinity level limit. 1.5 times is a known multiplier for the high value limit of the box plot, and this multiplier varies with the dispersion of the historical salinity samples.

[0044] It should be noted that if no historical three-dimensional state field for the same monitored soil area exists during the initial system run, a non-empty salt field with a value not less than 0 grams per kilogram is extracted from the deep sensing salinity data of the current cycle and the salinity field of the initial hydrothermal salinity field to obtain an initial salinity sample. The preset soil salinity level limit is then calculated based on this initial salinity sample. From the second cycle onwards, the salinity field from the target soil evolution simulation field generated in the previous cycle is incorporated into the historical salinity sample.

[0045] For example, the 25th percentile of historical salinity samples is 2.2 g / kg, the 75th percentile is 3.4 g / kg, and the interquartile range is 1.2 g / kg. Multiplying 1.2 g / kg by 1.5 gives 1.8 g / kg, and adding 1.8 g / kg to 3.4 g / kg gives 5.2 g / kg, which is then determined as the preset soil salinity level limit.

[0046] In one implementation, the salt concentration of each three-dimensional soil grid node at the current sampling time is numerically compared with a preset soil salinity level limit. If the salt concentration is greater than the preset soil salinity level limit, the corresponding three-dimensional soil grid node is marked as an abnormal node; if the salt concentration is less than or equal to the preset soil salinity level limit, the corresponding three-dimensional soil grid node is marked as a normal node.

[0047] It is worth noting that the mesh adjacency relationship is determined by the node coordinates of the 3D soil layer mesh. If two mesh nodes have a difference equal to the mesh spacing in only one of the horizontal, vertical, and depth coordinates, and the difference in the other two coordinate directions is 0, then the two mesh nodes are determined to be adjacent mesh nodes. This adjacency relationship corresponds to the six-sided adjacency relationship in the 3D soil layer mesh; two mesh nodes that only share edges or corners are not determined to be adjacent mesh nodes.

[0048] In one implementation, the abnormal nodes are marked as connected based on grid adjacency relationships. Specifically, the following steps are taken: First, all abnormal nodes are marked as unvisited. Then, one abnormal node is selected from the unvisited abnormal nodes as a starting node, an abnormal region identifier is assigned to the starting node, and the starting node is added to a node queue. Nodes are sequentially removed from the node queue, and abnormal nodes that satisfy the grid adjacency relationship with the removed nodes and have an unvisited access mark are found. These found abnormal nodes are assigned the same abnormal region identifier and added to the node queue. When the node queue is empty, all abnormal nodes with the same abnormal region identifier form a connected node set. This process is repeated until all abnormal nodes are assigned abnormal region identifiers, resulting in a salt anomaly region.

[0049] It should be noted that the salinity anomaly region includes one or more sets of connected nodes. Each set of connected nodes includes an anomaly region identifier, sampling time, a set of anomaly nodes, a set of three-dimensional coordinates corresponding to the anomaly nodes, and a set of salinity concentrations corresponding to the anomaly nodes. Subsequent steps read the salinity state of the previous period based on the set of three-dimensional coordinates corresponding to the anomaly nodes, and write the local salinity state based on the set of three-dimensional coordinates corresponding to the anomaly nodes.

[0050] In step S4, the previous cycle salinity state corresponding to the salinity anomaly area is obtained, the current salinity observation set is extracted from the soil monitoring data, the hydrothermal porosity state is extracted from the initial hydrothermal salt field, and the water-salt coupling state is updated according to the previous cycle salinity state, the current salinity observation set and the hydrothermal porosity state to obtain the local salinity state.

[0051] This includes obtaining the salinity status of the previous period corresponding to the salinity anomaly area and extracting the current salinity observation set from the soil monitoring data, including: Based on the grid node coordinates of the salinity anomaly zone, the salinity field and timestamp field of the corresponding grid node are read from the target soil evolution simulation field generated in the previous cycle to obtain the salinity status and timestamp of the previous cycle. Extract the coordinates of the current period's salinity observation point, the corresponding salinity observation value, and the current period's sampling time from the soil monitoring data; Calculate the spatial distance from the coordinates of each current period's salinity observation point to the coordinates of each grid node in the salinity anomaly zone, and determine the minimum spatial distance as the anomaly zone distance for the corresponding current period's salinity observation point; The current salinity observation set is obtained by selecting a preset number of current period salinity observation point coordinates, corresponding salinity observation values, and current period sampling time according to the distance of the anomaly zone from smallest to largest.

[0052] Specifically, the local salinity state is obtained by updating the water-salt coupling state based on the previous salinity state, the current salinity observation set, and the hydrothermal pore state, including: Extract the adjacent salinity state of the salinity anomaly area and surrounding grid nodes from the target soil evolution simulation field generated in the previous cycle, and extract the hydrothermal porosity state of the salinity anomaly area and surrounding grid nodes from the initial hydrothermal salinity field. The monitoring time step is calculated based on the timestamp of the previous period and the sampling time of the current period in the current salinity observation set; Based on the salinity state of the previous cycle, the adjacent salinity state, the hydrothermal porosity state, and the monitoring time step, a water-salt coupling migration recursion is performed to obtain the prior salinity state; The current grid observation value is obtained by performing a distance-weighted mapping based on the current salinity observation set and the grid node coordinates of the salinity anomaly area; Based on the observation residual between the prior salinity state and the current grid observation, the prior salinity state is corrected by Kalman filtering to obtain the local salinity state.

[0053] In one implementation, the previous cycle salinity state is the set of salinity fields of grid nodes in the target soil evolution simulation field generated in the previous cycle that have the same grid node coordinates as the salinity anomaly zone. When reading the previous cycle salinity state, firstly, all grid node coordinates in the salinity anomaly zone are extracted, and then grid nodes with the same coordinates are searched in the target soil evolution simulation field generated in the previous cycle. If a grid node with the same coordinates is found, then the salinity field and timestamp field of that grid node are read. If a grid node in a salinity anomaly zone does not have a grid node with the same coordinates in the target soil evolution simulation field generated in the previous cycle, then the salinity field of the grid node closest to that salinity anomaly zone grid node is read, and the timestamp field of that closest grid node is used as the corresponding timestamp field.

[0054] In one implementation, the current salinity observation set is a set of salinity observation records obtained by filtering soil monitoring data. Specifically, the process involves first filtering data records from the soil monitoring data whose observation field name is "saltiness," and deleting data records where the salinity field is null or less than 0 grams per kilogram. Then, the spatial coordinates, salinity field, and sampling time from the remaining data records are read and used as the coordinates of the salinity observation point for the current period, the corresponding salinity observation value, and the sampling time for the current period, respectively.

[0055] It should be noted that the anomaly zone distance represents the minimum three-dimensional spatial distance from the current period's salinity observation point to the salinity anomaly zone. For any current period's salinity observation point coordinates, the three-dimensional spatial distance between the current period's salinity observation point coordinates and the coordinates of each grid node in the salinity anomaly zone is calculated. Then, the minimum value is extracted from all three-dimensional spatial distances, and this minimum value is determined as the anomaly zone distance for the current period's salinity observation point.

[0056] It should be noted that the preset number is determined based on the spatial scale of the salinity anomaly region during the system initialization phase. Specifically, the number of grid nodes in the salinity anomaly region is counted, and the square root of this number is calculated and rounded up to obtain the candidate number. If the candidate number is less than 3, the preset number is set to 3; if the candidate number is greater than 8, the preset number is set to 8; if the candidate number is between 3 and 8, the candidate number is determined as the preset number. A lower limit of 3 ensures that the current salinity observation set includes observation points in multiple spatial directions, while an upper limit of 8 restricts the participation range of distant salinity observation points in local state updates.

[0057] In one implementation, the hydrothermal porosity state is a set of temperature, water content, and porosity fields corresponding to the salinity anomaly zone and surrounding grid nodes in the initial hydrothermal salt field. The surrounding grid nodes are obtained by extending the grid nodes of the salinity anomaly zone by one layer through grid adjacency relationships. Specifically, grid nodes that satisfy a six-sided adjacency relationship with any grid node in the salinity anomaly zone but do not belong to the salinity anomaly zone are extracted, and the extracted grid nodes are combined with the grid nodes in the salinity anomaly zone to form the salinity anomaly zone and surrounding grid nodes.

[0058] In one implementation, when calculating the monitoring time step, if there are multiple current period sampling times within the current salinity observation set, these multiple current period sampling times are arranged in chronological order, and the median sampling time is taken as the current period sampling time. Then, the time difference between the current period sampling time and the timestamp of the previous period is calculated, and this time difference is determined as the monitoring time step. If the time difference is less than 0, the current salinity observation set is discarded; if the time difference is equal to 0, the monitoring time step is set to one sampling period. The sampling period is determined by the difference between two adjacent sampling times on the reference time axis in step S1.

[0059] In one implementation, during the water-salt coupled migration recursion, for each grid node in the salinity anomaly zone, the salinity state of the previous cycle for that grid node is read, the salinity states of adjacent grid nodes that satisfy a six-sided adjacency relationship with that grid node are read, and the temperature, moisture content, and porosity fields of that grid node are also read. Subsequently, a moisture driving coefficient is calculated based on the moisture content field, a pore passage coefficient is calculated based on the pore field, an evaporation correction coefficient is calculated based on the temperature field, and the salinity migration increment is calculated based on the moisture driving coefficient, pore passage coefficient, evaporation correction coefficient, adjacent salinity states, and monitoring time step. This salinity migration increment is then superimposed on the salinity state of the previous cycle to obtain the prior salinity state of that grid node.

[0060] It should be noted that the water-salt coupled migration recursion is performed one by one according to the grid nodes in the salinity anomaly zone. For any grid node, the arithmetic mean of the salinity states of the six adjacent grid nodes of that grid node is subtracted from the salinity state of the previous cycle of that grid node to obtain the salinity state difference; the salinity state difference is then multiplied by the water driving coefficient, the pore passage coefficient, the evaporation correction coefficient, and the time step multiplier in sequence to obtain the salinity migration increment; the salinity migration increment is then added to the salinity state of the previous cycle to obtain the prior salinity state.

[0061] It should be noted that the moisture driving coefficient is determined based on the difference between the arithmetic mean of the moisture content fields of adjacent grid nodes and the moisture content field of the current grid node. The moisture content field is used in the calculation as a percentage value; that is, the value used in the calculation is 18 when the moisture content field is 18%, and the value used in the calculation is 22 when the moisture content field is 22%. If the arithmetic mean of the percentage values ​​of the moisture content fields of adjacent grid nodes is greater than the percentage value of the current grid node's moisture content field, the difference between the two values ​​is divided by 100 to obtain the moisture driving coefficient; if the arithmetic mean of the percentage values ​​of the moisture content fields of adjacent grid nodes is less than or equal to the percentage value of the current grid node's moisture content field, the moisture driving coefficient is 0. The pore passage coefficient is obtained by normalizing the pore field, specifically by dividing the percentage value of the pore field by 100. The evaporation correction coefficient is determined based on the local relative deviation of the temperature field. Specifically, the maximum, minimum, and median values ​​of the temperature field of the salinity anomaly area and surrounding grid nodes are extracted at the same sampling time. If the maximum and minimum values ​​are the same, the evaporation correction coefficient is set to 1. If the current grid node temperature field is greater than the median, the difference between the current grid node temperature field and the median is divided by the difference between the maximum and the minimum values, and the result is incremented by 1 to obtain the evaporation correction coefficient. If the current grid node temperature field is less than or equal to the median, the evaporation correction coefficient is set to 1. The time step multiplier is the monitoring time step divided by the sampling period in step S1, where the sampling period is determined by the difference between two adjacent sampling times on the reference time axis.

[0062] For example, the salinity of a grid node in a salinity anomaly zone in the previous cycle is 5.0 g / kg, and the arithmetic mean of the salinity of adjacent grid nodes is 4.4 g / kg; the moisture content field of this grid node is 18%, the arithmetic mean of the moisture content fields of adjacent grid nodes is 22%, and the porosity field is 40%; at the same sampling time, the maximum value of the temperature field of the salinity anomaly zone and surrounding grid nodes is 30 degrees Celsius, the minimum value is 20 degrees Celsius, the median is 25 degrees Celsius, and the temperature field of this grid node is 26 degrees Celsius. The monitoring time step is equal to the sampling period in step S1. Dividing the difference between 22 and 18 by 100 yields a moisture driving coefficient of 0.04; dividing 40 by 100 yields a porosity passage coefficient of 0.4; dividing the difference between 26 degrees Celsius and 25 degrees Celsius by the difference between 30 degrees Celsius and 20 degrees Celsius, and adding 1 to the result, yields an evaporation correction coefficient of 1.1; the time step multiplier is 1. Multiplying the difference between 4.4 g / kg and 5.0 g / kg by 0.04, 0.4, 1.1 and 1, we get a salt migration increment of -0.01056 g / kg; adding -0.01056 g / kg to 5.0 g / kg, we get a priori salt state of 4.98944 g / kg.

[0063] In one implementation, the current grid observation value is obtained by performing a distance-weighted mapping based on the current salinity observation set and the grid node coordinates of the salinity anomaly region. The specific operation is as follows: For any grid node in the salinity anomaly zone, read the horizontal, vertical, and depth coordinates of the grid node's coordinates, and read the horizontal, vertical, and depth coordinates of each salinity observation point in the current period of the current salinity observation set; calculate the differences between the grid node coordinates and each current period salinity observation point coordinates in the horizontal, vertical, and depth coordinates respectively, square the differences in the three directions respectively, add them together, and take the square root of the sum to obtain the three-dimensional spatial distance between the grid node coordinates and the corresponding current period salinity observation point coordinates; if the three-dimensional spatial distance is 0, the corresponding salinity observation value is directly determined as the current grid observation value of the grid node; if the three-dimensional spatial distances are not all 0, take the inverse square of each three-dimensional spatial distance, and normalize all the inverse squares to obtain the observation mapping weights corresponding to each current period salinity observation point; perform a weighted summation of the corresponding salinity observation values ​​according to the observation mapping weights to obtain the current grid observation value of the grid node.

[0064] It should be noted that the number of grid nodes in the current grid observation is the same as the number of grid nodes in the prior salinity state. The observation residual is calculated under the same grid node coordinates, specifically the difference between the current grid observation and the prior salinity state.

[0065] In one implementation, when performing Kalman filtering correction on the prior salinity state, the observation residuals are first calculated, and then the Kalman gain is calculated based on the process variance and the observation variance. The process variance is taken as the sample variance of the salinity field of the grid node and its six adjacent grid nodes in the target soil evolution simulation field generated in the previous cycle. The observation variance is taken as the sample variance of the corresponding salinity observation value in the current salinity observation set. The process variance uses local spatial salinity differences to characterize the state uncertainty in the recursive process, and the observation variance uses the dispersion of the current salinity observation set to characterize the uncertainty of the observation data. If the current salinity observation set contains only one corresponding salinity observation value, then the observation variance is taken as the sample variance of all salinity observation values ​​in the current salinity observation set at the same sampling time. The Kalman gain is the process variance divided by the sum of the process variance and the observation variance. If the sum of the process variance and the observation variance is 0, then the current grid observation value is determined as the local salinity state. If the sum of the process variance and the observation variance is not zero, the observation residual is multiplied by the Kalman gain, and the product is superimposed on the prior salinity state to obtain the local salinity state.

[0066] For example, the prior salinity of a grid node in a salinity anomaly area is 5.4 g / kg, the current grid observation value obtained by distance-weighted mapping is 5.8 g / kg, and the observation residual is 0.4 g / kg. The sample variance of the salinity field of this grid node and its six adjacent grid nodes in the previous period's target soil evolution simulation field is 0.30, and the sample variance of the salinity observation value in the current salinity observation set is 0.20. The Kalman gain is 0.30 divided by 0.50, resulting in 0.6. Multiplying 0.4 g / kg by 0.6 yields 0.24 g / kg, and adding 0.24 g / kg to 5.4 g / kg gives 5.64 g / kg, which is then determined as the local salinity state of this grid node.

[0067] It is worth noting that the local salinity status includes the grid node identifier, grid node coordinates, corrected salinity field, and current period sampling time for each grid node in the salinity anomaly zone. The local salinity status only updates the salinity field; the temperature, water content, and porosity fields are still retained from the initial hydrothermal salt field.

[0068] In step S5, the local salinity state is written into the salinity field of the grid node corresponding to the salinity anomaly zone in the initial hydrothermal salinity field to obtain the fused hydrothermal salinity field. The fused hydrothermal salinity field is then continuously corrected according to the boundary of the salinity anomaly zone to obtain the target soil evolution simulation field.

[0069] In one implementation, the local salinity state is the set of corrected salinity fields corresponding to each grid node within the salinity anomaly zone generated in step S4. The local salinity state includes grid node identifiers, grid node coordinates, corrected salinity fields, and the current period sampling time. The initial hydrothermal salt field is the three-dimensional state field generated in step S2, and each grid node in the initial hydrothermal salt field includes a grid node identifier, three-dimensional coordinates, a temperature field, a water content field, a salinity field, and a porosity field.

[0070] It should be noted that writing the local salinity state into the salinity field of the grid node corresponding to the salinity anomaly zone in the initial hydrothermal salt field specifically involves searching for grid nodes with the same coordinates in the initial hydrothermal salt field according to the grid node coordinates in the local salinity state. If a grid node with the same coordinates is found, the corrected salinity field from the local salinity state is written into the salinity field of that grid node; if no grid node with the same coordinates is found, the grid node closest to the grid node coordinates in the local salinity state is selected, and the corrected salinity field is written into the salinity field of that nearest grid node. During the writing process, the temperature, water content, and porosity fields retain their original values ​​from the initial hydrothermal salt field.

[0071] It is worth noting that the fused hydrothermal salt field is a three-dimensional state field after the local salinity state has been written. The fused hydrothermal salt field includes grid nodes within the salinity anomaly zone whose salinity fields have been updated, and grid nodes outside the salinity anomaly zone whose salinity fields have not been updated. Each grid node in the fused hydrothermal salt field is still stored according to its grid node identifier, three-dimensional coordinates, temperature field, water content field, salinity field, and porosity field.

[0072] In one implementation, the boundary of the salinity anomaly region is determined by the coordinates of the grid nodes and their adjacency relationships. Specifically, all grid nodes in the salinity anomaly region are traversed. If a grid node in the salinity anomaly region has at least one six-sided adjacent grid node that does not belong to the salinity anomaly region, then that grid node is designated as an inner boundary grid node. Grid nodes that satisfy the six-sided adjacency relationship with the inner boundary grid node but do not belong to the salinity anomaly region are designated as outer boundary grid nodes. The inner and outer boundary grid nodes together constitute the boundary of the salinity anomaly region.

[0073] It should be noted that the continuity correction of the fused hydrothermal salt field based on the boundary of the salinity anomaly zone specifically involves boundary smoothing of the salinity fields of the inner and outer boundary grid nodes. During boundary smoothing, the salinity field of the current boundary grid node is read, along with the salinity fields of adjacent grid nodes that satisfy a six-sided adjacency relationship with the current boundary grid node; adjacent grid nodes with null salinity fields are deleted; the arithmetic mean of the salinity field of the current boundary grid node and the remaining adjacent grid nodes is calculated, and the arithmetic mean result is written into the salinity field of the current boundary grid node.

[0074] It should be noted that if a boundary grid node is located on the outer edge of the three-dimensional soil layer grid, and the number of actual six adjacent grid nodes is less than six, only the actual adjacent grid nodes are read. If the salinity fields of all adjacent grid nodes of a boundary grid node are null, the salinity field of the boundary grid node in the fused hydrothermal salt field is retained.

[0075] For example, the salt content field of an inner boundary grid node is 5.64 g / kg, and the salt content fields of two valid adjacent grid nodes that satisfy a six-sided adjacency relationship with the inner boundary grid node are 5.10 g / kg and 4.80 g / kg, respectively. During boundary smoothing, the values ​​of 5.64 g / kg, 5.10 g / kg, and 4.80 g / kg are added together and divided by 3 to obtain 5.18 g / kg, which is then written into the salt content field of the inner boundary grid node.

[0076] It is worth noting that boundary smoothing only applies to the salinity field corresponding to the boundary of the salinity anomaly zone. The salinity field of non-boundary grid nodes inside the salinity anomaly zone retains the value after the local salinity state is written; the salinity field of non-boundary grid nodes outside the salinity anomaly zone retains the value in the initial hydrothermal salt field; the temperature field, water content field, and porosity field are not changed by boundary smoothing.

[0077] In one implementation, after boundary smoothing, the three-dimensional state field is re-encapsulated according to the current periodic sampling time, grid node identifier, three-dimensional coordinates, temperature field, moisture content field, salinity field, and porosity field, and this three-dimensional state field is determined as the target soil evolution simulation field. The target soil evolution simulation field also records the salinity anomaly zone identifier, the number of local salinity state records involved in the writing, and the number of boundary grid nodes that have completed continuous correction.

[0078] In summary, this invention performs spatiotemporal gridding processing on surface and deep sensing observation data to form a three-dimensional soil grid data containing both surface and deep grid values. Based on surface and deep constraints, continuous field estimation is performed on the three-dimensional soil grid nodes to obtain an initial hydrothermal-salinity field. Further, salinity anomaly areas are identified, and water-salinity coupling state updates are performed by combining the previous period's salinity state, the current salinity observation set, and the hydrothermal porosity state to obtain the local salinity state. Finally, the local salinity state is written into the initial hydrothermal-salinity field and boundary continuity correction is applied to obtain the target soil evolution simulation field. Thus, a three-dimensional continuous representation of soil hydrothermal-salinity state is achieved, and local salinity anomaly areas can be recursively updated based on current observation data and the previous period's state, reducing abrupt boundary changes between local update results and the surrounding soil state field.

[0079] A second embodiment of the present invention provides a digital simulation system for soil evolution processes, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method described above.

[0080] It should be noted that the soil evolution process digital simulation system provided in this embodiment of the invention is used to execute all the process steps of the soil evolution process digital simulation method in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0081] It should be noted that the device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0082] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A digital simulation method for soil evolution processes, characterized in that, include: Soil monitoring data, including salinity observation data and hydrothermal porosity observation data, is acquired from the monitored soil area. The soil monitoring data is then spatiotemporally gridded to obtain three-dimensional soil grid data, including surface grid values ​​and deep grid values. Surface constraints are generated based on the surface grid values, deep constraints are generated based on the deep grid values, and continuous field estimation is performed on the three-dimensional soil grid nodes based on the surface constraints and the deep constraints to obtain the initial hydrothermal salinity field. From the initial hydrothermal salt field, identify spatially continuous grid nodes where the salt concentration exceeds the preset soil salinity level limit to obtain salinity anomaly zones. Obtain the previous cycle salinity state corresponding to the salinity anomaly area, extract the current salinity observation set from the soil monitoring data, extract the hydrothermal porosity state from the initial hydrothermal salinity field, and update the water-salt coupling state based on the previous cycle salinity state, the current salinity observation set, and the hydrothermal porosity state to obtain the local salinity state. The local salinity state is written into the salinity field of the grid node corresponding to the salinity anomaly zone in the initial hydrothermal salinity field to obtain the fused hydrothermal salinity field. The fused hydrothermal salinity field is then continuously corrected according to the boundary of the salinity anomaly zone to obtain the target soil evolution simulation field.

2. The digital simulation method for soil evolution processes according to claim 1, characterized in that, The soil monitoring data includes surface remote sensing temperature data, surface remote sensing moisture content data, deep sensing temperature data, deep sensing moisture content data, deep sensing salinity data, and deep sensing porosity data.

3. The digital simulation method for soil evolution processes according to claim 1, characterized in that, The process of spatiotemporally gridding the soil monitoring data to obtain three-dimensional soil grid data containing surface grid values ​​and deep grid values ​​includes: Extract the sampling time, spatial coordinates, and observation values ​​from the soil monitoring data; The soil monitoring data is resampled according to the sampling time to obtain synchronous monitoring data; The synchronous monitoring data is mapped to a preset three-dimensional soil grid based on the spatial coordinates to obtain three-dimensional soil grid data containing surface grid values ​​and deep grid values.

4. The digital simulation method for soil evolution processes according to claim 1, characterized in that, The step of generating surface constraints based on the surface mesh values ​​and generating deep constraints based on the deep mesh values ​​includes: Extract the surface mesh values ​​and corresponding mesh node coordinates, and associate the surface mesh values ​​and corresponding mesh node coordinates with surface boundary data to obtain surface constraints; The deep mesh values ​​and corresponding mesh node coordinates are extracted, and the deep mesh values ​​and corresponding mesh node coordinates are associated with deep correction data to obtain deep constraints.

5. The digital simulation method for soil evolution processes according to claim 1, characterized in that, The step of performing continuous field estimation on the three-dimensional soil grid nodes based on the surface constraints and the deep constraints to obtain the initial hydrothermal saline field includes: Based on the surface constraints and the deep constraints, determine the known grid nodes with observed values; The grid nodes in the three-dimensional soil grid that lack observation values ​​are identified as grid nodes to be estimated. Based on the known grid nodes, perform three-dimensional spatial interpolation on the grid nodes to be estimated to obtain continuous grid values; An initial hydrothermal salt field is generated based on the continuous grid values.

6. The digital simulation method for soil evolution processes according to claim 5, characterized in that, The step of performing three-dimensional spatial interpolation on the grid nodes to be estimated based on the known grid nodes to obtain continuous grid values ​​includes: The observation values ​​and three-dimensional coordinates corresponding to the target fields are extracted from the known grid nodes. The target fields include temperature, water content, salinity, and porosity. Construct a three-dimensional semivariance function based on the observed values ​​and the three-dimensional coordinates; The Kriging weights corresponding to the grid nodes to be estimated are solved based on the three-dimensional semivariance function. The observations are weighted and summed according to the Kriging weights to obtain the continuous grid values ​​of the grid nodes to be estimated.

7. The digital simulation method for soil evolution processes according to claim 1, characterized in that, The process of identifying spatially continuous grid nodes in the initial hydrothermal salt field where the salt concentration exceeds a preset soil salinity level limit to obtain salinity anomaly zones includes: The salt concentration of each grid node is extracted from the initial hydrothermal salt field; Grid nodes whose salt concentration exceeds the preset soil salinity level limit are identified as abnormal nodes; The abnormal nodes are connected based on the grid adjacency relationship to obtain the salt anomaly region.

8. The digital simulation method for soil evolution processes according to claim 1, characterized in that, The step of obtaining the previous period's salinity status corresponding to the salinity anomaly area and extracting the current salinity observation set from the soil monitoring data includes: Based on the grid node coordinates of the salinity anomaly zone, the salinity field and timestamp field of the corresponding grid node are read from the target soil evolution simulation field generated in the previous cycle to obtain the salinity status and timestamp of the previous cycle. Extract the coordinates of the current period's salinity observation point, the corresponding salinity observation value, and the current period's sampling time from the soil monitoring data; Calculate the spatial distance from the coordinates of each current period's salinity observation point to the coordinates of each grid node in the salinity anomaly zone, and determine the minimum spatial distance as the anomaly zone distance for the corresponding current period's salinity observation point; The current salinity observation set is obtained by selecting a preset number of current period salinity observation point coordinates, corresponding salinity observation values, and current period sampling time according to the distance of the anomaly zone from smallest to largest.

9. The digital simulation method for soil evolution processes according to claim 8, characterized in that, The step of updating the water-salt coupling state based on the previous salinity state, the current salinity observation set, and the hydrothermal pore state to obtain the local salinity state includes: Extract the adjacent salinity state of the salinity anomaly area and surrounding grid nodes from the target soil evolution simulation field generated in the previous cycle, and extract the hydrothermal porosity state of the salinity anomaly area and surrounding grid nodes from the initial hydrothermal salinity field. The monitoring time step is calculated based on the timestamp of the previous period and the sampling time of the current period in the current salinity observation set; Based on the salinity state of the previous cycle, the adjacent salinity state, the hydrothermal porosity state, and the monitoring time step, a water-salt coupling migration recursion is performed to obtain the prior salinity state; The current grid observation value is obtained by performing a distance-weighted mapping based on the current salinity observation set and the grid node coordinates of the salinity anomaly area; Based on the observation residual between the prior salinity state and the current grid observation, the prior salinity state is corrected by Kalman filtering to obtain the local salinity state.

10. A digital simulation system for soil evolution processes, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the method according to any one of claims 1 to 9.

Images