A method and system for analyzing lead pollution in paddy field soil

Abstract

The application discloses a kind of rice field soil lead pollution analysis method and system, method includes: obtaining the current soil attribute data set of target rice field and soil moisture content distribution data, combine future expected water management parameters, simulate moisture content space-time variation sequence by soil water migration model;Dissolution-precipitation equilibrium tendency state index of each time step Pb is calculated based on moisture content change and soil attribute, the advancing track of Pb active front is identified;The front track is overlaid with the current rice Pb overproof risk spatial distribution, and the potential Pb activation amplification area is determined;Along the front advancing direction, strip-shaped analysis unit is divided, and Pb pollution risk evolution index is calculated by combining front arrival time and current risk value;Finally, the analysis report containing amplification area position, risk evolution order and leading factor diagnosis result is output.The technical leap of rice field Pb pollution risk from "static status evaluation" to "dynamic evolution prediction" is realized.

Description

Technical Field

[0001] This invention belongs to the field of soil environmental monitoring technology, and in particular relates to a method and system for analyzing lead pollution in paddy field soil. Background Technology

[0002] Rice, as a major food crop in my country, has a strong ability to accumulate phosphorus (Pb) in the soil. Excessive Pb content in rice directly threatens food security and public health. After being absorbed by the plant's roots, Pb in rice paddies can be transported and accumulated in the rice grains. Long-term consumption of rice with excessive Pb levels can cause chronic damage to multiple systems, including the nervous, hematopoietic, and digestive systems.

[0003] A prominent problem with current technologies for the safe utilization and risk assessment of phosphorus-contaminated paddy fields is that they focus only on the total phosphorus content in the soil and a few static physicochemical indicators, lacking a systematic understanding of the intrinsic relationship between the dynamic water management unique to paddy fields and changes in phosphorus activity. Paddy fields experience alternating periods of flooding and drying during the planting cycle, and the soil redox state fluctuates with water changes. Pb exhibits significant dissolution-precipitation equilibrium migration under different redox conditions. However, existing technologies cannot predict how the phosphorus activity front will advance under a specific water management plan, or which areas will experience a significant increase in phosphorus pollution risk driven by water. This limits risk analysis to a descriptive level of current static information, making it difficult to proactively predict future pollution risks.

[0004] Therefore, there is an urgent need to develop an analytical method that can combine future water management plans for paddy fields, predict the dynamic advancement of Pb active fronts, identify potential Pb activation and amplification areas, and assess the evolution trend of pollution risk in each area, so as to provide a scientific basis for the accurate monitoring and early warning of Pb-contaminated paddy fields. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for analyzing Pb pollution in paddy field soil, which solves the technical problems of existing technologies being unable to predict the advancement process of Pb active fronts in conjunction with future water management plans and being unable to assess the evolution trend of Pb pollution risk in various regions.

[0006] In a first aspect, the present invention provides a method for analyzing lead pollution in paddy field soil, comprising: Obtain the current soil property dataset and current soil moisture content distribution data of the target paddy field. The current soil property dataset includes at least the total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity of the soil. Obtain the expected water management parameters for the target paddy field within a future preset time period; The current soil moisture content distribution data and the expected water management parameters are input into a preset soil moisture transport model to simulate the spatiotemporal change sequence of soil moisture content within a preset future time period. Based on the spatiotemporal variation sequence of soil moisture content and the current soil property dataset, calculate the spatial distribution of state indicators characterizing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step; Based on the spatial distribution of the state index, the advancing trajectory of the Pb active front is identified. The Pb active front is the projection of the critical isosurface of the state index from precipitation tendency to dissolution tendency in horizontal space. The advancing trajectory of the Pb active front is overlaid with the spatial distribution of the current Pb exceedance risk in the target paddy field to determine the area that the Pb active front is expected to sweep through and where the current predicted Pb value in the paddy exceeds the safety threshold, as the potential Pb activation and amplification area. Within the potential Pb activation and amplification region, the region is divided into multiple strip-shaped analysis units perpendicular to the advancing direction of the Pb active front. Based on the expected arrival time sequence of the Pb active fronts in each analysis unit and the average Pb exceedance risk value of each unit, the Pb pollution risk evolution index of each analysis unit is calculated. The output includes an analysis report containing the spatial location of the potential Pb activation and amplification region, the ranking of Pb pollution risk evolution indices for each analysis unit, and the diagnostic results of key soil factors that dominate Pb activity changes within each analysis unit.

[0007] Secondly, the present invention provides a paddy field soil lead pollution analysis system, comprising: The first acquisition module is configured to acquire the current soil attribute dataset and the current soil moisture content distribution data of the target paddy field. The current soil attribute dataset includes at least the total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity of the soil. The second acquisition module is configured to acquire the expected water management parameters of the target paddy field in a future preset period. The simulation module is configured to input the current soil moisture content distribution data and the expected water management parameters into a preset soil moisture transport model to simulate the spatiotemporal change sequence of soil moisture content within a preset future time period. The calculation module is configured to calculate the spatial distribution of state indicators characterizing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step, based on the spatiotemporal variation sequence of the soil moisture content and the current soil attribute dataset. The identification module is configured to identify the advancing trajectory of the Pb active front based on the spatial distribution of the state index, wherein the Pb active front is the projection of the critical isosurface of the state index from precipitation tendency to dissolution tendency in horizontal space. The determination module is configured to overlay the advancement trajectory of the Pb active front with the spatial distribution of the current Pb exceedance risk of the target paddy field to determine the area that the Pb active front is expected to sweep through and the current predicted Pb value of the paddy exceeds the safety threshold, as the potential Pb activation and amplification area. The partitioning module is configured to divide the potential Pb activation and amplification region into multiple strip-shaped analysis units perpendicular to the advancing direction of the Pb active front. The calculation module is configured to calculate the Pb pollution risk evolution index of each analysis unit based on the expected arrival time sequence of the Pb active fronts in each analysis unit and the average Pb exceedance risk value of each unit. The output module is configured to output an analysis report containing the spatial location of the potential Pb activation and amplification region, the ranking of Pb pollution risk evolution indices for each analysis unit, and the diagnostic results of key soil factors that dominate Pb activity changes within each analysis unit.

[0008] Thirdly, an electronic device is provided, comprising: at least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the steps of the paddy field soil lead pollution analysis method according to any embodiment of the present invention.

[0009] Fourthly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the program instructions are executed by a processor, the processor performs the steps of the method for analyzing lead pollution in paddy field soil according to any embodiment of the present invention.

[0010] This application presents a method and system for analyzing lead pollution in paddy field soil. By combining future water management plans for paddy fields with a dynamic model of soil Pb dissolution-precipitation balance, it achieves advanced prediction of the trajectory of Pb active fronts. By overlaying the trajectory of Pb active fronts with the current spatial distribution of Pb exceedance risk in rice, it accurately identifies potential Pb activation and amplification areas most likely to experience future Pb pollution escalation. By dividing these potential Pb activation and amplification areas into analytical units along the frontal advancement direction and calculating the Pb pollution risk evolution index, it achieves quantitative ranking and graded early warning of pollution risks in each area. The final analysis report clearly provides diagnostic results of the dominant driving factors for each area, providing a scientific basis for subsequent precise monitoring and targeted regulation. This method represents a technological leap from "static status evaluation" to "dynamic evolution prediction," significantly improving the predictability and scientific rigor of Pb pollution risk analysis in paddy fields. Attached Figure Description

[0011] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0012] Figure 1 A flowchart of a method for analyzing lead pollution in paddy field soil provided in an embodiment of the present invention; Figure 2 A structural block diagram of a paddy field soil lead pollution analysis system provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0013] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, 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.

[0014] Please see Figure 1 The diagram shows a flowchart of a method for analyzing lead pollution in paddy field soil according to this application.

[0015] like Figure 1 As shown, the specific steps in the analysis method for lead pollution in paddy field soil include: Step S101: Obtain the current soil attribute dataset and current soil moisture content distribution data of the target paddy field. The current soil attribute dataset includes at least the total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity.

[0016] In this step, multiple monitoring points are set up in the target paddy field according to a preset grid, surface soil samples are collected at each monitoring point, and various indicators in the current soil attribute dataset are measured. At each monitoring point, in-situ moisture content sensors were used to acquire the current soil volumetric moisture content data at each point. Spatial interpolation is performed on the current soil volumetric moisture content data at each monitoring point to obtain the current soil moisture content distribution data.

[0017] In one specific embodiment, a systematic grid method is used to deploy monitoring points within the target paddy field. The grid spacing is determined based on the field area and terrain complexity, generally set at 50 meters × 50 meters; for fields with significant terrain variations or uneven expected pollution distribution, the spacing can be increased to 25 meters × 25 meters. The starting point of the grid is selected at a clearly marked location on the field boundary, and the grid direction is consistent with the main cultivation direction of the field. Using high-precision differential GPS equipment or real-time dynamic differential positioning equipment, the actual geographic coordinates of each grid node are located on-site, and the latitude, longitude, and elevation information are recorded, with the error controlled within sub-meter level. For cases where the actual location of a grid node is inaccessible due to field obstacles, an alternative point representing the soil conditions of that grid point is selected within a 2-meter radius, and the offset coordinates are recorded. Finally, a spatial distribution layer file containing several monitoring points is generated, with each monitoring point having a unique number and spatial coordinates.

[0018] At each established monitoring point, soil samples from the top 0-20cm depth were collected using a stainless steel soil auger. Five samples were collected from each point, one from the center and one from each of the surrounding areas. The soil samples from these five augers were thoroughly mixed on a clean plastic sheet and then divided into quarters. Approximately 1kg of the final mixed sample was placed in a clean polyethylene self-sealing bag and labeled with the site number, sampling date, and field name. Samples were protected from compression and high temperatures during transport and delivered to the laboratory promptly. In the laboratory, the soil samples were spread evenly on a clean enamel plate or plastic film and allowed to air dry naturally in a well-ventilated environment away from direct sunlight. During air drying, the samples were turned occasionally, and large clumps were gently broken up along natural cracks. The dried samples were then crushed with a wooden stick and passed through a 2mm nylon sieve. The sieved material was thoroughly mixed and then sealed in a wide-mouthed glass bottle for later testing. Some of the samples were further ground using an agate mortar and pestle and passed through a 100-mesh nylon sieve for testing projects requiring finer-grained soil.

[0019] For the treated soil samples, the following five key attribute indicators were measured according to currently effective agricultural industry standard methods: Total Pb content in soil: Approximately 0.2 g of soil sample passed through a 100-mesh sieve was weighed and digested using a nitric acid-hydrofluoric acid-perchloric acid triacid system. After acid removal and volume adjustment, the digestion solution was determined using a graphite furnace atomic absorption spectrophotometer or inductively coupled plasma mass spectrometer. The method detection limit should be below 0.1 mg / kg. Results are expressed in mg / kg.

[0020] Soil pH value: Weigh 10.0g of air-dried soil sample that has passed through a 2mm sieve, place it in a 50mL beaker, add 25mL of distilled water to remove carbon dioxide, with a water-to-soil ratio of 2.5:1, stir vigorously with a glass rod for 1-2 minutes, and let stand for 30 minutes. Use a pH meter glass electrode calibrated with a standard buffer solution to measure the pH value of the supernatant, and record the reading after it stabilizes.

[0021] Soil organic matter content: Weigh approximately 0.2g of soil sample that has passed through a 100-mesh sieve, place it in a hard glass test tube, add potassium dichromate-sulfuric acid solution, heat to boiling in an oil bath at 170-180℃ for 5 minutes, cool, and titrate with ferrous sulfate standard solution. Using o-phenanthroline as an indicator, calculate the organic carbon content based on the amount of potassium dichromate consumed, multiply by the conversion factor of 1.724 to obtain the organic matter content, and express the result in g / kg.

[0022] Available Si content in soil: Weigh 5.0 g of air-dried soil sample that has passed through a 2 mm sieve, add acetic acid buffer for extraction, shake for 1 hour and then filter. The silicon in the filtrate is determined by the silicomolybdenum blue colorimetric method, with the absorbance measured at a wavelength of 810 nm. The available silicon content is calculated according to the standard curve and the result is expressed in mg / kg.

[0023] Soil cation exchange capacity: Weigh 5.0 g of air-dried soil sample that has passed through a 2 mm sieve, rinse it multiple times with ammonium acetate solution to remove excess ammonium ions, and then wash it with ethanol. The amount of ammonium ions adsorbed by the soil is determined by the Kjeldahl method, and the cation exchange capacity is calculated. The result is expressed as cmol / kg.

[0024] For each batch of samples, blank samples, standard substance samples, and no less than 10% parallel samples are simultaneously set up for quality control. The measured values ​​of the standard substances should be within the allowable error range of the standard values, and the relative deviation of the parallel samples should be controlled within 15%.

[0025] At each monitoring point, portable time-domain reflectometers or frequency-domain reflectometers were used to acquire the current soil volumetric moisture content data in situ. During measurement, the probe was vertically inserted into the top 0-20cm of soil, ensuring close contact between the probe and the surrounding soil. The reading was recorded after it stabilized. Data was taken three times at each point, including the center and three times at different locations around the perimeter. The arithmetic mean was taken as the representative value of the current soil volumetric moisture content at that point. The measurement time should be as consistent as possible with the soil sample collection time to ensure temporal matching between soil moisture content and soil property data. For fields where high salinity may affect the accuracy of dielectric measurements, a small amount of undisturbed soil sample could be collected simultaneously, and the sensor readings could be calibrated using a drying and weighing method.

[0026] The current soil volumetric moisture content data of each monitoring point was correlated with its spatial coordinates and imported into the spatial analysis module of the geographic information system software. Ordinary Kriging interpolation was used for spatial interpolation. First, exploratory statistical analysis was performed on the moisture content data to check if it conformed to a normal distribution. If the deviation was large, necessary logarithmic transformation or normal score transformation was performed. Then, semivariance function modeling was performed, and a spherical model, exponential model, or Gaussian model was selected for fitting. The optimal model and its parameters were determined based on the goodness-of-fit index. During the interpolation process, the cell size of the output raster was set, generally to one-fifth to one-tenth of the grid spacing. For example, when the grid spacing was 50 meters, the raster cell size was set to 5-10 meters. After interpolation, the raster area exceeding the boundary of the study area was masked to obtain raster data of the current soil moisture content distribution covering the entire target paddy field. Each raster cell in this data has a current soil volumetric moisture content value, forming complete current soil moisture content distribution data, which serves as the initial input for subsequent soil moisture transport simulations.

[0027] Step S102: Obtain the expected water management parameters for the target paddy field in the future within a preset time period.

[0028] In one specific embodiment, the determination of the future preset time period must be consistent with the planting system of the target rice field and the time span of subsequent soil moisture transport simulations. Typically, the current sampling date is used as the starting point, covering the period from the next complete rice growing season to harvest. For single-season rice planting areas, the preset time period is generally from the current date to 100 to 140 days in the future; for double-season rice planting areas, the preset time period is from the current date to 60 to 90 days before the harvest of the current season's rice. The start and end dates of the time period are clearly recorded as the time domain boundaries for subsequent simulations. When determining the time period, it is necessary to consult farmers or check local agricultural management records to confirm key growth period time nodes such as rice variety, transplanting date, expected start and end times of field drying, and expected harvest date, ensuring that the time period has sufficient overlap with the rice's most sensitive growth period to soil moisture.

[0029] The expected methods for obtaining water management parameters include: Historical data from IoT-based water monitoring devices installed in the fields can be used to infer water management practices. For fields without written irrigation plans but with deployed field water monitoring devices, continuous monitoring data on surface water depth or soil moisture content from the same period of the previous year can be retrieved. Analysis of historical data identifies characteristic parameters such as the alternation cycle of flooding and drying, the approximate time intervals between irrigations, and the typical range of flood depths. Combined with an investigation of the field's irrigation and drainage facilities, the feasibility of the water management model derived from historical data in the current year can be verified. Once confirmed, this model can be used as the expected water management parameters.

[0030] The acquired expected water management parameters are converted into a structured digital format recognizable by subsequent soil moisture transport models. The data structure of the expected water management parameters is a two-dimensional table or an equivalent structured file, containing the following three fields: The first field is the time series, recording the date or relative time of each water management event. It increments daily, starting from the current sampling day as time zero, until the end of a preset future time period. The resolution of the time series is generally set to 1 day, but can be increased to 12 hours or 6 hours for periods of frequent water management changes.

[0031] The second field is the flood depth value, recording the depth of the flooded layer on the field surface at the corresponding time point, in millimeters. A flood depth of zero indicates that there is no water layer on the field surface, and the soil is in an unsaturated state; a positive flood depth indicates that there is water of a corresponding depth on the field surface, and its value is equal to the height of the water level above the soil surface. Between two irrigation events, if there is no rainfall replenishment and drainage operation, the flood depth decreases linearly according to the preset daily average evaporation and seepage rate. The evaporation and seepage rate is determined based on local meteorological data or existing research literature, and is generally taken in the range of 5 to 15 mm / day.

[0032] The third field is the irrigation identifier, recording whether an irrigation event occurred at that time. 0 indicates no irrigation that day, and 1 indicates irrigation that day. For irrigation times, the irrigation water volume, in millimeters, should also be noted. The irrigation water volume is the quota for this irrigation, and its value is equal to the difference between the flood depth after irrigation and the flood depth before irrigation, plus the evaporation and seepage losses that occurred during the irrigation process.

[0033] The three fields above are arranged in a time series to form a complete curve showing the expected change in flood depth from the current moment to the end of a preset future period. This curve serves as the upper boundary condition for soil moisture transport simulation in subsequent steps. The mathematical expression of this curve is h(t), a function of time t, where h is the flood depth at the corresponding time t. The function shows a step increase in flood depth at the instant of irrigation, and an approximately linear decrease between two irrigations. This pattern accurately reflects the water management characteristics of alternating flooding and drying in paddy fields.

[0034] Taking a specific scenario as an example: In a single-season rice field, the fixed irrigation interval is 7 days, and the target flooding depth after each irrigation is 40 mm. The average daily evaporation and seepage loss rate between irrigations is approximately 6 mm / day. The expected water management parameters for this field are as follows: On each irrigation day, the flooding depth jumps from the lower pre-irrigation level to 40 mm; thereafter, the flooding depth decreases by approximately 6 mm each day, reaching approximately 4 mm on the 6th day; on the 7th day, irrigation resumes and the water depth returns to 40 mm. During the planned field drying period, irrigation is stopped on the first day of drying, and the flooding depth continues to decrease to zero and remains there for a period until irrigation resumes after drying. Recording this information completely as a flooding depth change curve completes the digital representation of the expected water management parameters.

[0035] Step S103: Input the current soil moisture content distribution data and the expected water management parameters into a preset soil moisture transport model to simulate the spatiotemporal change sequence of soil moisture content within a preset future time period.

[0036] In this step, the preset soil moisture transport model is a numerical model of unsaturated soil moisture movement based on the Richards equation. The simulated spatiotemporal variation sequence of soil moisture content within a preset future time period includes: Using the current soil moisture content distribution data as the initial condition, the expected flooding depth change curve as the upper boundary condition, and free drainage or a set groundwater level as the lower boundary condition, the Richards equation is iteratively solved at a preset time step to obtain a raster map of the spatial distribution of soil moisture content at each time step.

[0037] In one specific embodiment, the preset soil moisture transport model is a numerical model of unsaturated soil moisture movement based on the Richards equation. The Richards equation is the fundamental governing equation describing soil moisture movement in the unsaturated zone, and its one-dimensional vertical flow form is: , Where θ is the soil volumetric water content; t is time, in days; z is the vertical coordinate, in meters, with the soil surface as the origin and downward as positive; h is the soil matrix potential, in meters; K(h) is the unsaturated hydraulic conductivity, in meters per day, a function of matrix potential; and S(h) is the root water uptake term, in days. -1 For fallow periods without crop growth or simulation scenarios where root water absorption is not precisely considered, the value can be zero.

[0038] This equation describes the change in soil volumetric water content over time as equal to the difference between the vertical gradient of soil water flux and root water uptake.

[0039] The governing equations were discretized spatially and temporally using the finite difference method. Spatially, the soil profile was divided into several node layers, with closer node spacing at the surface (e.g., 1-2 cm) and gradually increasing spacing downwards (e.g., 5-10 cm). Each horizontal spatial location corresponded to an independent soil column, with no lateral water exchange between the columns. Temporally, an implicit difference scheme was used to ensure numerical stability. The time step was adaptively adjusted based on the computational convergence, with an initial time step of 0.01 days and a maximum allowable step size of no more than 1 day. The nonlinear algebraic equations within each time step were solved using the Picard iteration method, with the iteration convergence criterion set at a maximum absolute value of the matric potential change between two iterations being less than 0.01 m. Besides the finite difference method, the finite element method could also be used for discretization, employing mature open-source code such as HYDRUS-1D or SWAP, or a self-developed equivalent numerical program.

[0040] The numerical solution process is carried out independently on soil columns at each horizontal spatial location, and the differences in water content between soil columns are only reflected in the spatial variability given the initial conditions. The solution result for each soil column is a time series volumetric water content profile. By recombining the solution results of all soil columns according to their horizontal spatial locations, a three-dimensional spatiotemporal variation sequence of soil water content for the entire target paddy field is obtained.

[0041] The model requires parameters including soil moisture characteristic curve parameters and saturated hydraulic conductivity parameters. The soil moisture characteristic curve describes the nonlinear relationship between soil matric potential and volumetric water content, expressed using the van Genuchten equation: θ(h) = θr + (θs - θr) / [1 + |αh| n ] m Where θr is the residual volumetric water content, θs is the saturated volumetric water content, α is the reciprocal parameter of the inlet value, n is the pore size distribution index parameter, and m is related to n through the relationship m = 1 - 1 / n. The unsaturated hydraulic conductivity function K(h) is described by the Mualem-van Genuchten model.

[0042] The above parameters are obtained as follows: In the target paddy field, undisturbed soil samples were collected. Soil water holding capacity data under multiple pressure levels were measured in the laboratory using a pressure membrane apparatus. The van Genuchten equation was fitted to obtain the values ​​of each parameter. The saturated hydraulic conductivity was determined through a constant head permeability experiment. This method has the highest parameter accuracy.

[0043] If the target paddy field contains multiple different soil texture types, values ​​should be assigned according to the spatial distribution of each texture type, and different sets of hydraulic parameters should be used for soil columns in different texture types. The spatial distribution information of soil texture comes from field soil survey maps or remote sensing interpretation data.

[0044] The initial conditions are the current soil moisture content distribution data obtained in step S101. For each soil column at a horizontal spatial location, the initial moisture content within the 0-20cm depth range of the surface layer is set according to the current soil volumetric moisture content value at that location. For moisture content profiles below 20cm, if there is no measured data, it is set to 80% of the field capacity based on experience or estimated according to the capillary rise equation based on the groundwater level depth. A model "warm-up" period is set in the early stage of the simulation to allow the deep moisture content profile to reach natural equilibrium.

[0045] The upper boundary condition is the expected flood depth variation curve generated in step S102. For periods where the flood depth is greater than zero, the upper boundary is set as a constant head boundary, where the pressure head value equals the flood depth; for periods where the flood depth is zero and the surface soil moisture content is below saturation, the upper boundary is set as an atmospheric boundary, allowing for evaporative water loss. The evaporation rate is set as the daily average evaporation rate based on local meteorological data and crop water requirements.

[0046] The lower boundary conditions are set based on the hydrogeological conditions of the target paddy field. When the groundwater depth is large and free drainage at the bottom of the soil profile is unobstructed, the lower boundary is set as a free drainage boundary, i.e., a unit gradient boundary, with gravity drainage as the main factor. When the groundwater level is shallow and its seasonal fluctuation pattern is known, the lower boundary is set as a variable head boundary, with the time function of the groundwater level depth as the lower boundary head value. When there is an impermeable layer underground, the lower boundary is set as a zero flux boundary.

[0047] After assigning model parameters and setting initial and boundary conditions, the numerical simulation program is started. The simulation begins at the current time and progresses along the time axis until the preset time period ends. At each time step, the program performs the following operations: The first step is to load the current upper boundary conditions (flood depth) into the surface nodes of the model; The second step is to calculate the unsaturated hydraulic conductivity and specific water capacity of the current iteration step at each spatial node of each soil column based on the matrix potential value of the previous time step. The third step is to assemble the coefficient matrix, solve the discretized Richards equations, and obtain the matrix potential and volumetric water content values ​​at each node at the current time step. The fourth step is to check the mass conservation error. If the cumulative mass error exceeds the preset allowable value, the time step is reduced and the calculation is repeated. Fifth, after convergence is achieved, the volumetric water content values ​​of each node at that time step are written to the output file. The time is advanced by one step, and the process returns to the first step to continue the calculation.

[0048] During the simulation, the volumetric water content values ​​for each time step, each grid position, and each vertical node are continuously recorded. After the simulation, the volumetric water content within a depth range of 0-20cm in the surface layer is extracted. The average volumetric water content of the surface soil at each grid point at each output time step is obtained using a depth-weighted average method or by taking representative node values. All grid points and the average volumetric water content of the surface layer at all output time steps are organized into a time series to form a spatiotemporal variation sequence of soil water content within the predetermined future time period. The physical format of this sequence is as follows: for each time step, a soil water content spatial distribution raster map with the same spatial range and grid resolution as the current soil water content distribution data is output; the raster maps for all time steps are arranged in chronological order to form a complete time series dataset.

[0049] The output time step is set to 1 day, meaning one raster map of the spatial distribution of moisture content is output daily. For periods of frequent changes in water management, the output interval can be increased to 12 hours. All output raster maps are stored in GeoTIFF format with embedded geospatial reference information, facilitating overlay analysis with other spatial data in subsequent steps.

[0050] Step S104: Based on the spatiotemporal variation sequence of soil moisture content and the current soil attribute dataset, calculate the spatial distribution of state indicators characterizing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step.

[0051] In this step, for each grid point at each time step, the simulated values ​​of soil moisture content, soil pH, and organic matter content of the grid point are obtained; Calculate the soil water-filling porosity based on the simulated soil moisture content; The soil water-filled porosity, soil pH, and organic matter content are input into a preset Pb dissolution-precipitation equilibrium tendency function to calculate the state index values ​​of the grid points at the time step. The preset Pb dissolution-precipitation equilibrium tendency function is a multivariate nonlinear function used to characterize the Pb dissolution-precipitation equilibrium tendency under given water-filled porosity, pH, and organic matter content. 2+ Thermodynamic stability tendency of PbCO3 and Pb(OH)2 precipitates in soil solution relative to solid PbCO3 and Pb(OH)2 precipitates.

[0052] In one specific embodiment, three types of input data are required: the first type is the spatiotemporal variation sequence of soil moisture content output in step S103, i.e., the spatial distribution raster map of soil moisture content at each time step; the second type is the soil pH value and organic matter content in the current soil attribute dataset obtained in step S101; and the third type is soil bulk density data, which is used to calculate total porosity.

[0053] First, spatial matching of the three types of data is completed. Soil pH and organic matter content are point measurements at each monitoring point, which need to be converted into continuous raster surface data with the same spatial range and raster resolution as the soil moisture content raster map. Using the same ordinary kriging interpolation method as in step S101, spatial interpolation is performed on the pH and organic matter content at each monitoring point to generate spatial distribution raster maps for pH and organic matter content. The semivariogram function model parameters in the interpolation process can be kept consistent with those in the soil moisture content interpolation to ensure that the three raster layers are completely aligned in space and that each raster cell corresponds one-to-one. The soil bulk density data is processed in the same way as pH and organic matter content. Spatial interpolation is performed on the soil bulk density values ​​measured at each monitoring point to generate a bulk density spatial distribution raster map. If the bulk density difference between monitoring points is small, the arithmetic mean of the bulk density of all points can be taken as the representative value of the entire area, and a bulk density layer can be established in the form of a constant value raster map.

[0054] After spatial matching is completed, for each time step, the data record for each raster point includes the following five attributes: simulated soil moisture content, soil pH, organic matter content, soil bulk density, and soil texture classification code. These attribute data are stored in a unified spatial database, indexed by time step and raster row and column number, for quick retrieval in subsequent calculations.

[0055] Soil water-filled porosity is defined as the ratio of current soil volumetric water content to total soil porosity. It is a dimensionless index, ranging from 0 to 1, and characterizes the degree to which soil pores are filled with water. This index directly controls the oxygen diffusion rate and redox potential in the soil and is one of the key driving variables for the Pb dissolution-precipitation equilibrium.

[0056] The calculation process consists of three steps. The first step is to calculate the total soil porosity at each grid point based on its soil bulk density and soil particle density. The calculation formula is: φ=1-(ρb / ρs), Where φ is the total porosity of the soil, which is dimensionless; ρb is the soil bulk density, in g / cm³. 3 ρs represents the soil particle density, which is typically taken as 2.65 g / cm³ for paddy soil. 3When the soil organic matter content is high, the soil particle density can be appropriately adjusted based on the organic matter content.

[0057] The second step is to read the simulated soil volumetric water content θ of the grid point from the spatial distribution raster map of soil moisture content at that time step. This value is the average volumetric water content of the surface soil at that grid point at that time step, calculated in step S103, and its value is equal to the percentage of the current water content of the grid point to the total soil volume.

[0058] The third step is to divide the soil volumetric water content by the total porosity to obtain the soil water-filled porosity: WFPS=θ / φ, Where WFPS is the soil water-filled porosity, dimensionless; θ is the current soil volumetric water content; and φ is the total soil porosity, dimensionless. When WFPS is close to 0, it indicates that there is very little water in the soil pores, oxygen supply is sufficient, and oxidation reactions dominate; when WFPS is close to 1, it indicates that the soil pores are almost completely filled with water, oxygen diffusion is hindered, the soil is in a hypoxic or anaerobic state, and reduction reactions dominate.

[0059] The preset Pb dissolution-precipitation equilibrium tendency function is a multivariate nonlinear function used to quantify the Pb dissolution-precipitation tendency under given soil water-filled porosity, pH, and organic matter content. 2+ Thermodynamic stability tendency of PbCO3 and Pb(OH)2 precipitates in soil solution relative to solid PbCO3 and Pb(OH)2 precipitates.

[0060] The theoretical basis for this function is that the dissolution-precipitation equilibrium of Pb in paddy soil is controlled by the coupling effect of the solubility product equilibrium of Pb-containing minerals and the adsorption-desorption equilibrium of soil colloid surfaces. In a typical flooded paddy soil environment, PbCO3 and Pb(OH)2 are the controlling factors for Pb concentration. 2+ The main solid-phase precipitation forms of activity. When the Pb in the soil solution... 2+ When the product of the activity of Pb and the activities of other anions exceeds the solubility product of the corresponding mineral, Pb 2+ Pb tends to precipitate and be removed from solution, thus reducing its bioavailability; conversely, precipitates tend to dissolve, thus reducing Pb bioavailability. 2+ When released into the solution, the activity of Pb increases.

[0061] The function has the following form: SI=f(WFPS,pH,SOM)=k1×ln(WFPS+δ)+k2×(pH-pH ref ) 2 -k3×ln(SOM+ε)+k0; Where SI is the state index value, which is dimensionless; a higher value indicates higher Pb. 2+The stronger the solubility, the higher the activity; the lower the value, the stronger the precipitation tendency and the higher the degree of Pb fixation; WFPS is soil water-filled porosity; pH ​​is soil pH value; SOM is soil organic matter content, unit is g / kg; pH ref The reference pH value is set to the pH value corresponding to the lowest solubility of Pb(OH)2 precipitate, approximately between 9 and 10; δ is the water filling porosity adjustment parameter, set to 0.01 to prevent the logarithm from being meaningless when WFPS is zero; ε is the organic matter adjustment parameter, set to 1.0 to prevent the logarithm from being meaningless when SOM is zero; k1, k2, and k3 are the influence weight coefficients of each factor, and k0 is the intercept constant. All of the above coefficients were pre-calibrated using experimental data.

[0062] The four parameters k1, k2, k3, and k0 in the function were determined through the following experimental calibration procedure: Multiple sets of undisturbed soil samples of different types and different levels of Pb contamination were collected from the target paddy field area. Under controlled laboratory conditions, multiple soil moisture content gradients were set (e.g., WFPS of 0.2, 0.4, 0.6, 0.8, and 1.0). After equilibration for 48 hours under each moisture content gradient, the Pb content in the soil solution was measured. 2+ Activity; pH ​​and organic matter content of each soil sample were measured simultaneously; all experimental data were compiled into a training dataset, and nonlinear least squares method was used to measure Pb. 2+ Using soil activity as a reference target, the four parameters in the function are optimally fitted. The calibrated function is added to the model parameter library for subsequent calculations. For different soil type regions, appropriate parameter sets can be calibrated separately to improve the regional applicability of the model.

[0063] After completing function calibration and data preparation, the batch calculation phase begins. The calculation process iterates through each time step and each grid point.

[0064] For the current time step and the current grid point, perform the following operations: read the current simulated soil moisture content of the grid point from the spatial database, calculate the soil water-filled porosity (WFPS) for the current time step according to the aforementioned steps; read the soil pH and organic matter content of the grid point; substitute the three values ​​of WFPS, pH and SOM into the calibrated Pb dissolution-precipitation equilibrium tendency function to calculate the state index value SI.

[0065] After the calculation is completed, the SI value is assigned to the current grid point and stored in association with the grid point's spatial row and column number and time step index.

[0066] After performing the above calculations on all grid points within the current time step, a complete spatial distribution raster map of SI values ​​covering the entire target paddy field is obtained. The above process is repeated for all time steps to obtain a sequence of spatial distribution raster maps of SI values ​​for all time steps. These raster maps are arranged in chronological order, and the output is a spatial distribution dataset of state indicators representing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step. The format of this dataset is consistent with the input spatiotemporal variation sequence of soil moisture content. Each time step corresponds to one GeoTIFF format raster map, and the raster resolution, spatial extent, and projected coordinate system are completely consistent with the input data, ensuring that pixel-level overlay operations can be directly performed in subsequent steps.

[0067] Step S105: Based on the spatial distribution of the state index, identify the advancing trajectory of the Pb active front, wherein the Pb active front is the projection of the critical isosurface of the state index from precipitation tendency to dissolution tendency onto the horizontal space.

[0068] In this step, the state index value of each grid point at a certain time step is compared with a preset state critical threshold, which is the boundary value between precipitation tendency and dissolution tendency. Connect the grid points whose state index values ​​are exactly equal to the state critical threshold to form a dissolution-precipitation boundary line at a certain time step; The dissolution-precipitation boundary lines of all time steps are superimposed in chronological order to form a set of trajectory lines that continuously advance from the initial position into or out of the target paddy field. The set of trajectory lines constitutes the advancing trajectory of the Pb active front.

[0069] In one specific embodiment, the preset critical threshold is the boundary between precipitation tendency and dissolution tendency. Its physical meaning is: when the state index value equals this threshold, the Pb concentration in the soil solution... 2+ The activity of Pb reaches thermodynamic equilibrium with the solubility product of the PbCO3 or Pb(OH)2 solid phase; when the state index value is higher than this threshold, Pb 2+ When the activity exceeds the equilibrium level, the system tends to shift towards dissolution; when the state index value is below this threshold, Pb... 2+ When the activity is below the equilibrium level, the system tends to move towards precipitation.

[0070] The method for determining the critical threshold of a state is as follows: Determine the solubility product constant (Ksp) of the main Pb-containing solid minerals in the target paddy field soil. The solubility product constant of PbCO3 is Ksp(PbCO3) = 7.4 × 10⁻⁶. -14The corresponding negative logarithm of the solubility product constant is pKsp = 13.13; the solubility product constant of Pb(OH)₂ has two crystal forms, among which the solubility product constant of Pb(OH)₂ formed by the hydration of orthorhombic PbO is Ksp(Pb(OH)₂) = 1.43 × 10⁻⁶. -20 The corresponding negative logarithm of the solubility product constant is pKsp = 19.85. The literature standard value at 25℃ (normal temperature) is chosen as the basis for calculation. If actual soil temperature variations need to be considered, Ksp can be corrected for temperature using the van der Hoff equation.

[0071] Determine the equilibrium anion activity in the soil solution. The activity of Pb in the soil solution... 2+ The paired anions are mainly CO3. 2- and OH - CO3 2- Soil CO2 activity is controlled by soil CO2 partial pressure and pH. In flooded paddy fields, soil CO2 partial pressure is typically 1 × 10⁻⁶. -2 Up to 5×10 -2 Based on this, combined with Henry's law and the second-order dissociation equilibrium of carbonic acid, CO3 can be calculated. 2- Activity; OH - The activity is directly calculated from pH, aOH - =Kw / aH + Where Kw is the ion product constant of water, and aOH - aH represents the activity of hydroxide ions in the soil solution. + This represents the activity of hydrogen ions in the soil solution. The calculation is based on the median value within the typical pH range of paddy fields.

[0072] Calculate equilibrium Pb using the solubility product formula. 2+ Activity. For PbCO3, the activity of lead ions aPb in the soil solution at equilibrium. 2+ = Ksp(PbCO3) / aCO3 2- aCO3 2- This represents the activity of carbonate ions in the soil solution; for Pb(OH)₂, at equilibrium, aPb 2+ =Ksp(Pb(OH)2) / (aOH - ) 2 Under given pH and CO2 partial pressure conditions, the lower aPb calculated from two minerals was selected. 2+ The equilibrium activity threshold for Pb precipitation under these conditions, i.e., the mineral phase that first reaches saturation and precipitates under these conditions, controls the Pb precipitation. 2+ The upper limit of activity.

[0073] Establishing a balance of Pb 2+The correspondence between activity and state index values. Using the measurement data of a series of soil samples from the function calibration experiment in step S104 under known WFPS and pH conditions as training data, the equilibrium Pb under each condition was simultaneously measured during the calibration experiment. 2+ Activity, through inverse regression to balance Pb 2+ By mapping the activity value to the dimensions of the state index value, the state index value SI corresponding to the Pb dissolution-precipitation equilibrium critical point is determined. critical The SI critical This is the preset critical threshold value, whose value is specifically calibrated through the above steps and stored in the model parameter library for use in each subsequent frontal identification operation.

[0074] For each time step of the state index spatial distribution raster map output in step S104, the dissolution-precipitation boundary line is extracted. The extraction process is as follows: The first step is to load the spatial distribution raster map of the state indicators for the current time step. Each cell in this raster map stores a state indicator value SI, and the cell row and column numbers correspond to the horizontal spatial coordinates.

[0075] The second step is to compare the SI value of each pixel in the raster image with the preset state threshold SI. critical Compare the values. For each cell, perform the following judgment: if SI is greater than SI... critical The pixel is marked as a region prone to dissolution; if SI is less than SI critical This cell is marked as a precipitation-prone region; if SI is exactly equal to SI critical This pixel is the candidate pixel for the boundary.

[0076] The third step is to use the contour tracing algorithm to extract the SI value, which equals SI. critical The precise spatial location. Because raster data is discretely sampled, the SI value is exactly equal to SI. critical The number of pixels may be small, requiring linear interpolation to determine the precise boundary location. Specifically, this involves iterating through all adjacent pixel pairs in the raster image; if the SI value of one pixel in an adjacent pixel pair is greater than SI... critical The SI value of another pixel is less than SI. critical If a segment with a dissolution-precipitation boundary exists between adjacent pixel pairs, the SI value is determined by linear interpolation on the connection line between the two pixels. critical The precise coordinates of these points are determined. Connecting all these precise coordinates according to their spatial adjacency forms one or more continuous contour lines. These contour lines are defined as the dissolution-precipitation boundary at the current time step.

[0077] The fourth step involves vectorizing and post-processing the extracted dissolution-precipitation boundary lines. The contour lines are converted into GIS vector line feature format, assigned a time attribute, and labeled with their respective time step number and corresponding date. Extremely short, fragmented lines caused by raster noise are filtered out by setting a minimum segment length threshold. Secondary closed loops formed by local spatial anomalies are determined based on their area and topological relationship with the major boundary line, deciding whether to retain them as independent secondary boundaries or remove them as noise.

[0078] Fifth, repeat steps one through four for all time steps to form a vector dataset of dissolution-precipitation boundaries for each time step. Each boundary line in this dataset corresponds to a time step and records the spatial boundary position between the precipitation-prone area and the dissolution-prone area at that moment.

[0079] The dissolution-precipitation boundary lines of all time steps generated in the previous step are spatially superimposed in chronological order to synthesize the propagation trajectory of the Pb active front. The specific implementation is as follows: Dissolution-precipitation boundaries for all time steps are loaded into a unified geographic information system environment, using the same spatial reference coordinate system. Each boundary has a unique attribute identifier, including the time step number and the corresponding simulation date.

[0080] Spatial tracking is performed on the boundary lines between adjacent time steps according to chronological order. A shortest distance matching algorithm is used to establish the correspondence between nodes on the boundary lines between two adjacent time steps. For each sampling point on the boundary line of an earlier time step, the nearest spatially distant corresponding point is searched on the boundary line of the next time step; the vector connecting the two corresponding points represents the local advancing direction and distance of that point within a time step.

[0081] All sampling points and their local advance vectors along the time step boundaries are concatenated in chronological order to form a set of trajectory lines describing the continuous advance of the Pb active front from the initial moment to the final moment. Specifically, starting from the sampling point on the first time step boundary, corresponding points are sequentially connected to the second, third, and final time step boundaries, forming several spatial trajectory lines extending continuously from the starting position into or out of the target paddy field. These trajectory lines reflect the movement path of the Pb active front during the simulation period, pointing in the direction of expansion or contraction of the Pb dissolution tendency region.

[0082] Statistical analysis was performed on the trajectory set to extract the macroscopic characteristics of frontal advance. The total advance length, average advance rate, and advance direction angle of each trajectory were calculated. Spatial clustering was performed on the average advance direction and advance rate of each trajectory to identify the rapid advance zone and stagnation zone of the Pb active front.

[0083] The trajectory set is stored and output as a vector line feature format. Each trajectory line is accompanied by attribute fields, including the start time step, end time step, total advance distance, average advance rate, and average advance direction angle. This trajectory set constitutes the advance trajectory of the Pb active front, providing direct input for the overlay analysis in step S106.

[0084] Step S106: Overlay analysis is performed on the advancing trajectory of the Pb active front and the spatial distribution of the current Pb exceedance risk in the target paddy field to determine the area that the Pb active front is expected to sweep through and where the current predicted Pb value in the paddy exceeds the safety threshold, as the potential Pb activation and amplification area.

[0085] In this step, the spatial distribution of the current risk of Pb exceeding the standard in the target paddy field is obtained through the following steps: Input the current soil attribute dataset of each monitoring point into the preset rice Pb prediction model to calculate the predicted rice Pb value of each monitoring point. The rice Pb prediction model is a function model that is pre-established based on the multiple linear regression relationship between total soil Pb, pH, organic matter, available Si and cation exchange capacity and rice Pb content. Spatial interpolation was performed on the predicted Pb values ​​of rice at each monitoring point to generate a continuous distribution map of the predicted Pb values ​​of rice. In the continuous distribution map of predicted Pb values ​​in rice, the regions where the predicted Pb values ​​of rice are greater than a preset safety threshold are extracted as the current spatial distribution of Pb exceeding the standard risk in rice.

[0086] In one specific embodiment, the spatial distribution of the risk of excessive phosphorus (Pb) in rice is obtained by calculating the predicted Pb values ​​for each monitoring point using a preset rice Pb prediction model, then generating a continuous distribution map through spatial interpolation, and finally extracting the areas exceeding the standard based on a safety threshold to obtain raster data. The generation process is as follows: Sub-step S1061: Construction of the rice Pb prediction model. The rice Pb prediction model is a pre-established function model based on the multiple linear regression relationship between five key factors—total soil Pb content, pH value, organic matter content, available Si content, and cation exchange capacity—and rice Pb content. The function form of this model is: Pb grain =a1×Pb soil +a2×pH+a3×SOM+a4×Si avail +a5×CEC+b, Among them, Pb grain The predicted Pb values ​​for rice are in mg / kg. soilTotal phosphorus (Pb) content in soil, in mg / kg; pH is the soil pH value, dimensionless; SOM is the soil organic matter content, in g / kg; Si avail 1 represents the available Si content in the soil, in mg / kg; CEC represents the soil cation exchange capacity, in cmol / kg; a1 to a5 are the regression coefficients of each factor; b is the intercept constant.

[0087] The regression coefficients a1 to a5 and the intercept constant b of the model were determined by fitting multiple linear regression analysis to soil-rice co-sampling data from the target paddy field area. The co-sampling data needed to cover paddy fields with different soil types, different Pb pollution levels, and different planting systems in the area, and the sample size should meet the statistical requirements of multiple regression analysis, typically no less than 100 groups. In the regression analysis, five key soil factors at each monitoring point were used as independent variables, and the measured Pb content of harvested rice at the corresponding point was used as the dependent variable. The least squares method was used to solve for the regression coefficients, and significance tests were performed to ensure that each factor coefficient was statistically significant. The overall model passed the F-test.

[0088] Sub-step S1062: Calculation of predicted Pb values ​​for rice at each monitoring point. The current soil attribute dataset (total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity) obtained in step S101 for each monitoring point in the target paddy field is substituted into the aforementioned rice Pb prediction model to calculate the predicted Pb value for each monitoring point. This predicted value is a numerical value in mg / kg, representing the theoretical value of the expected Pb content of rice in the corresponding area under the current soil conditions.

[0089] Sub-step S1063: Generation of a continuous distribution map of predicted Pb values ​​for rice. The predicted Pb values ​​and their spatial coordinates for each monitoring point are correlated and imported into the spatial analysis module of the geographic information system software. Spatial interpolation is performed using the same ordinary kriging interpolation method as in step S101. Before interpolation, exploratory statistical analysis is conducted on the predicted data to verify whether it conforms to a normal distribution; data transformation is performed if necessary. The semivariogram function model type and parameters are consistent with those of the previous interpolation, and the cell size of the output raster is completely consistent with the soil moisture content distribution raster map generated in step S101, ensuring precise spatial alignment of the two raster layers during subsequent overlay analysis. After interpolation, a raster map of predicted Pb values ​​for rice covering the entire continuous surface of the target paddy field is obtained, with each raster cell storing one predicted Pb value.

[0090] Sub-step S1064: Extraction of the spatial distribution of current rice Pb exceeding the standard risk. The preset safety threshold is the Pb content limit for rice specified in the National Food Safety Standard for Maximum Levels of Contaminants in Food (GB 2762-2022), which is 0.2 mg / kg. In the continuous distribution map of predicted Pb values ​​in rice, a binarization classification operation is performed: raster pixels with predicted Pb values ​​greater than 0.2 mg / kg are marked as 1, representing areas with exceeding the standard risk; raster pixels with predicted Pb values ​​less than or equal to 0.2 mg / kg are marked as 0, representing areas that meet the standard. The generated binarized raster map is the spatial distribution of the current rice Pb exceeding the standard risk in the target paddy field. The spatial range covered by the pixels marked as 1 in this raster map represents the area where the Pb content in rice may exceed the national food safety standard under the current soil conditions.

[0091] The set of Pb active front propulsion trajectory lines output in step S105 is converted into spatial coverage area data that can be used for overlay analysis.

[0092] Load all trajectory vector data output in step S105. Each trajectory line records the complete motion path of the Pb active front from the initial time to the final time.

[0093] Spatial aggregation is performed on all trajectory lines. Using the smallest convex polygon bounded by the start and end points of all trajectory lines as the base contour, a buffer extension is made outwards by a certain distance. This buffer distance can be set to twice the size of the raster cell in step S103 to ensure that the frontal influence range is not missed. The buffered polygonal region serves as the total influence range of the Pb active front.

[0094] The total impact area is converted into raster surface data with the same raster resolution, spatial range, and projected coordinate system as the current spatial distribution of Pb exceedance risk in rice. Raster cells falling within the total impact area are marked as 1, indicating that the cell is located within the region swept by the Pb active front; cells outside the range are marked as 0. This binarized raster map is the raster representation of the region swept by the Pb active front.

[0095] The raster layers generated in the above two steps are overlaid pixel by pixel to define the potential Pb activation amplification region.

[0096] Load the "Pb active front swept area" raster layer and the "Spatial distribution of current rice Pb exceedance risk" raster layer into the same geographic information system environment. Ensure that the cell size, spatial range, and number of rows and columns of the two layers are completely consistent, and that all corresponding cells overlap spatially.

[0097] Perform a cell-by-cell logical AND operation. For each raster cell, read its attribute values ​​from the two raster maps mentioned above. A cell is marked as a potential Pb-activated amplification region cell only if it meets both of the following conditions: Condition 1: The value of this pixel in the "Pb active front swept area" grid is 1, that is, this pixel is located within the expected advance path of the Pb active front and will be swept by the Pb active front in the future. Condition 2: The value of this pixel in the "Spatial Distribution of Current Rice Pb Exceedance Risk" raster map is 1, that is, the current predicted value of rice Pb in this pixel is greater than the safety threshold of 0.2 mg / kg, and there is already a risk of exceeding the standard.

[0098] A new binary raster image is generated from the result of the logical AND operation. Pixels that simultaneously satisfy both of the above conditions are assigned a value of 1, indicating that the pixel belongs to the potential Pb activation and amplification region; the remaining pixels are assigned a value of 0. The continuous spatial range formed by the pixels marked with 1 is the desired potential Pb activation and amplification region.

[0099] Spatial filtering is performed on the raster map of potential Pb activation and amplification areas. For scattered and isolated pixels caused by raster noise or edge effects, small patches below a minimum patch area threshold are merged into adjacent main patches or removed. The minimum patch area threshold can be set to a value equivalent to the area of ​​the grid cell at the monitoring point, such as 2500 square meters, to ensure that the identified potential Pb activation and amplification areas have practical management significance.

[0100] The above raster-format overlay analysis results are converted into a vector format that facilitates subsequent calculations, and relevant attribute information is added.

[0101] The raster image of the potential Pb-activated amplification region obtained in step three is vectorized. A boundary tracing algorithm is used to extract the boundaries of the continuous set of cells labeled 1, generating one or more closed polygonal features. Each polygonal feature represents a spatially continuous potential Pb-activated amplification region.

[0102] Calculate and assign attribute fields to each polygon. Attribute fields include: polygon area, in square meters or hectares; polygon perimeter, in meters; average predicted Pb value of rice in raster cells inside the polygon; estimated time range of the polygon's interior swept by the Pb active front; and distribution of the main exceedance risk levels inside the polygon.

[0103] Each polygon is output and stored as a vector polygon file, serving as the final spatial data result for the potential Pb activation and amplification regions. This result, along with the description field, constitutes the direct input for the region division in step S107, allowing subsequent operations to focus solely on these clearly defined regions. Simultaneously, in the subsequent output analysis report, this vector polygon file will serve as the underlying data for the first part, "Spatial Location Map of Potential Pb Activation and Amplification Regions," clearly displaying the specific location and extent of these regions within the target paddy field in the form of a spatial distribution map.

[0104] Step S107: Within the potential Pb activation and amplification region, the region is divided into multiple strip-shaped analysis units perpendicular to the advancing direction of the Pb active front.

[0105] In this step, the trajectory line of the advancing direction of the Pb active front in the potential Pb activation amplification region is extracted; Multiple dividing points are set along the trajectory line of the propulsion direction with a preset spatial step size; At each segmentation point, a segmentation line perpendicular to the propulsion direction trajectory line is generated, and each segmentation line extends to both sides until it intersects the boundary of the potential Pb activation amplification region; The region enclosed by two adjacent dividing lines and the two side boundaries of the potential Pb activation amplification region is defined as a strip-shaped analysis unit.

[0106] In one specific embodiment, from the set of Pb active front propulsion trajectory lines output in step S105, trajectory segments that fall into the potential Pb activation and amplification region determined in step S106 are extracted, and their dominant propulsion direction is determined.

[0107] Load the complete set of Pb-active front advancement trajectories output in step S105. This dataset is in vector line feature format, with each trajectory line accompanied by attribute fields for advancement direction and advancement rate. Simultaneously, load the vector surface file of the potential Pb-activated amplification region output in step S106.

[0108] Using the boundaries of the potential Pb activation and amplification region as the clipping range, a spatial clipping operation is performed on the set of Pb active front advancement trajectories. The spatial relationship between each trajectory and the polygon of the potential Pb activation and amplification region is determined: if all or part of the trajectory segment falls inside the region polygon, the part falling inside is retained; if the trajectory line is completely outside the region polygon, it is discarded. After clipping, a dataset containing only trajectory segments within the target region is obtained.

[0109] Principal direction analysis was performed on all trajectory segments falling within the same potential Pb activation and amplification region. The azimuth angle of each trajectory segment was calculated (with true north as 0° and expressed as clockwise angles). The azimuth angles of all trajectory segments were weighted and averaged according to their lengths to obtain the frontal weighted average advancing azimuth angle for that region. The weighting method was as follows: the length of each trajectory segment was used as the weight, multiplied by the azimuth angle of that segment, summed, and then divided by the total length of all segments.

[0110] Starting from the geometric center of the potential Pb activation and amplification region, a main propagation direction trajectory line is generated, extending outwards along the weighted average propagation direction and traversing the entire polygonal region. This trajectory line runs through the entire potential Pb activation and amplification region, and its two ends intersect with the boundaries of the polygonal region, forming the starting and ending points of the main propagation direction trajectory line.

[0111] Dividing points are set up along the main propulsion trajectory line according to the preset spatial step length.

[0112] Determine the preset spatial step size. The value of the spatial step size is set according to the size of the target paddy field and the required management precision, generally ranging from 10 meters to 30 meters. For smaller fields or scenarios requiring higher analysis precision, a smaller step size should be used; for larger fields, a larger step size can be used to reduce subsequent computational load. The preset spatial step size should be uniformly determined before processing and kept constant within the same potential Pb activation and amplification region.

[0113] Starting from the beginning of the trajectory line in the main propulsion direction, segmentation points are sequentially set at equal intervals along the trajectory line towards the end point, according to a preset spatial step size. Let the total length of the trajectory line be L, and the preset spatial step size be Δs. Then the number of segmentation points is n = floor(L / Δs) + 1, where floor represents rounding down. The i-th segmentation point is located on the trajectory line at a distance of (i-1) × Δs from the starting point. All segmentation points are arranged along the trajectory line, covering the complete segment from the starting point to the end point.

[0114] Record the spatial coordinates and mileage position of each segment point on the trajectory line. Each segment point is assigned a unique number, numbered sequentially from the starting point to the ending point of the trajectory line as P1, P2, ..., Pn. At the same time, record the proportion of the segment point in the total length.

[0115] At each segmentation point, a segmentation line perpendicular to the main propulsion direction trajectory line is generated and extends to both sides until it intersects the boundary of the potential Pb activation amplification region.

[0116] Calculate the local tangent direction of the trajectory line at each dividing point. For a dividing point located in the middle of the trajectory line, take the direction of the line connecting the two adjacent dividing points before and after the dividing point as the local tangent direction; for dividing points located at the starting point and the ending point of the trajectory line, take the direction of the line connecting the starting point and the second dividing point, and the direction of the line connecting the ending point and the penultimate dividing point as the local tangent directions of the starting point and the ending point, respectively.

[0117] Rotate the local tangent direction by 90° (clockwise or counterclockwise, but it must be consistent throughout the calculations) to obtain the partitioning direction at that partitioning point. The partitioning direction is perpendicular to the local tangent direction of the trajectory line at that partitioning point, representing the analytical direction that laterally crosses the potential Pb activation and amplification region.

[0118] Starting from each segmentation point, rays are generated in two opposite directions along the segmentation direction, extending to both sides. The extension distance is set to the diagonal length of the maximum bounding rectangle of the potential Pb activation amplification region, ensuring that the rays completely penetrate the region and terminate outside the region.

[0119] Perform an intersection operation between each ray and the boundary polygon of the potential Pb-activated amplification region. Each ray intersects the region boundary polygon at two points, located on either side of the dividing point. Record the spatial coordinates of these two intersection points. The line segment from intersection point A to intersection point B is the complete dividing line corresponding to that dividing point. This dividing line runs through the potential Pb-activated amplification region, and its two endpoints fall precisely on the region boundary.

[0120] Repeat the above operation for all split points to generate a series of parallel split lines. The distance between two adjacent split lines, measured along the main propulsion direction trajectory, is approximately equal to the preset spatial step size Δs.

[0121] By utilizing adjacent partition lines and the boundaries of potential Pb-activated amplification regions, closed strip-shaped analysis unit polygons are constructed and attribute information is assigned to them.

[0122] The first step is to construct the polygon geometry of the strip-shaped analysis unit. For the i-th subdivision line L... i and the (i+1)th subdivision line L i+1 (i from 1 to n-1), with the dividing line L i The line segment extending to the left from the dividing point to the boundary of the region is the upper segment of the left boundary, L. i The line segment extending to the right from the dividing point to the boundary of the region is the lower segment of the left boundary; the dividing line L i+1 The two corresponding line segments are the upper right boundary segment and the lower right boundary segment; the boundary of the potential Pb activation amplification region is at L i With L i+1 The left arc segment is the upper boundary, and the right arc segment is the lower boundary. Connecting these four boundary segments end to end forms a closed polygon, which is the i-th strip-shaped analysis unit.

[0123] Handle the special cases at the beginning and end. For the first analysis unit (i=1), one side of its outer boundary is the region boundary at the starting point of the trajectory line in the main propulsion direction; for the last analysis unit (i=n-1), one side of its outer boundary is the region boundary at the ending point of the trajectory line in the main propulsion direction. The polygons of these two analysis units are closed by the natural region boundary on one side along the trajectory line direction.

[0124] Perform a topology consistency check on all constructed analysis cells. Check for gaps or overlaps between adjacent analysis cells. If any exist, eliminate them by fine-tuning the endpoints of the partition lines to ensure that all analysis cells seamlessly cover the entire potential Pb activation amplification region and do not overlap with each other.

[0125] Assign attribute information to each analysis unit. The attribute information includes: analysis unit number, in order from the starting point to the ending point of the trajectory line in the advancement direction; analysis unit area, in square meters, calculated from polygon geometry; the range of subdivision line numbers contained within the analysis unit; and the sorting position of the analysis unit in the advancement direction.

[0126] All analysis cells are output in vector polygon file format, with each cell being an independent polygon feature and including the aforementioned attribute fields. These analysis cells are arranged sequentially along the advance direction of the Pb active front, forming the basic calculation unit for calculating the Pb pollution risk evolution index in step S108. The spatial width of each analysis cell reflects the regional characteristics within a range of Δs along the advance direction, and its spatial arrangement is perpendicular to the front advance direction, enabling it to accurately capture the risk change gradient during the front advance process.

[0127] Step S108: Based on the expected arrival time sequence of the Pb active fronts in each analysis unit and the average Pb exceedance risk value of each unit, calculate the Pb pollution risk evolution index of each analysis unit.

[0128] In this step, the arithmetic mean of the predicted Pb values ​​of rice for all grid points within a certain analysis unit is calculated as the current average risk value for that analysis unit. Based on the advancing trajectory of the Pb active front, the estimated time when a certain analytical unit is swept by the Pb active front is determined; The current average risk value is weighted and combined with the expected time to calculate the Pb pollution risk evolution index.

[0129] In one specific embodiment, the time when the Pb active front first enters each analytical unit is determined using the dissolution-precipitation boundary vector data for each time step generated in step S105. Specifically: Load the dissolution-precipitation boundary vector dataset for each time step output in step S105. Each boundary line includes a time step number and the corresponding simulation date attribute. Simultaneously load the analysis cell vector surface file.

[0130] For each analytical unit, determine the estimated time when it will be swept by the Pb active front. Specifically, in chronological order of time steps, spatially intersect the dissolution-precipitation boundary line of each time step with the polygon of the analytical unit. The time step in which the boundary line first intersects the polygon or falls within it is the time when the Pb active front first enters that analytical unit.

[0131] Mathematically, consider a set of dissolution-precipitation boundaries {L1, L2, ..., Lm} arranged in chronological order, where L1 is the initial time boundary and Lm is the final time boundary. For the analysis unit Ui, find the minimum value of k such that the spatial intersection relationship between Lk and Ui is true. Then, the date Tk corresponding to the k-th time step is the estimated time T when the analysis unit is swept by the Pb active front. arrival(i) .

[0132] If an analytical cell does not intersect any dissolution-precipitation boundary at any time step, it indicates that the Pb active front failed to advance to that cell during the entire simulation period. Such analytical cells are marked as front-unreached areas, and their estimated time is set to the simulation termination time plus a maximum value, or they are treated separately and assigned a specific low-urgency assignment in the subsequent risk evolution index calculation.

[0133] The determined estimated time is recorded in days, counting the number of days from the current moment, and assigned to each analysis unit as the estimated time. Simultaneously, the duration of the frontal sweep across the analysis unit is recorded, i.e., the time span from the first entry to the last exit.

[0134] Step S109: Output an analysis report containing the spatial location of the potential Pb activation and amplification region, the ranking of Pb pollution risk evolution indices for each analysis unit, and the diagnostic results of key soil factors that dominate Pb activity changes within each analysis unit.

[0135] In this step, for each analysis unit, the correlation coefficients between soil pH, organic matter content, available Si content, cation exchange capacity, and total Pb content and the state indicators are calculated. The soil factor with the largest absolute value of the correlation coefficient was identified as the dominant factor of the analysis unit; The analysis report presents the ranking of Pb pollution risk evolution indices for each analysis unit in the form of a spatial distribution map, and lists the dominant factors and the average deviation of the dominant factors for each analysis unit in the form of tables or annotations.

[0136] In summary, the method of this application acquires the current soil attribute dataset and soil moisture content distribution data of the target paddy field, combines it with expected future water management parameters, and simulates the spatiotemporal change sequence of moisture content through a soil moisture transport model; based on the moisture content change and soil attributes, it calculates the Pb dissolution-precipitation equilibrium tendency index at each time step to identify the advancement trajectory of the Pb active front; it overlays the front trajectory with the current spatial distribution of Pb exceedance risk in rice to determine potential Pb activation and amplification areas; it divides the analysis into strip-shaped units along the front advancement direction, and calculates the Pb pollution risk evolution index by combining the front arrival time and the current risk value; finally, it outputs an analysis report including the location of the amplification area, the risk evolution ranking, and the diagnostic results of the dominant factors; thus achieving a technological leap from "static status evaluation" to "dynamic evolution prediction" of Pb pollution risk in paddy fields.

[0137] Please see Figure 2 The diagram shows a structural block diagram of a paddy field soil lead pollution analysis system according to this application.

[0138] like Figure 2 As shown, the paddy field soil lead pollution analysis system 200 includes a first acquisition module 210, a second acquisition module 220, a simulation module 230, a calculation module 240, an identification module 250, a determination module 260, a division module 270, a calculation module 280, and an output module 290.

[0139] The system comprises the following modules: a first acquisition module 210, configured to acquire current soil property datasets and current soil moisture content distribution data for the target paddy field, wherein the current soil property dataset includes at least total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity; a second acquisition module 220, configured to acquire expected water management parameters for the target paddy field within a preset future time period; a simulation module 230, configured to input the current soil moisture content distribution data and the expected water management parameters into a preset soil moisture transport model to simulate the spatiotemporal variation sequence of soil moisture content within a preset future time period; a calculation module 240, configured to calculate the spatial distribution of state indicators characterizing the Pb dissolution-precipitation equilibrium tendency in the soil at each time step based on the spatiotemporal variation sequence of soil moisture content and the current soil property dataset; and an identification module 250, configured to identify the advancing trajectory of the Pb active front based on the spatial distribution of the state indicators, wherein the Pb active front is the state indicator transitioning from precipitation tendency to dissolution tendency. The transformation is defined by the projection of the critical isosurface in horizontal space; the determination module 260 is configured to overlay the advancement trajectory of the Pb active front with the spatial distribution of the current Pb exceedance risk of the target paddy field to determine the area that the Pb active front is expected to sweep through and where the current predicted Pb value of the paddy exceeds the safety threshold, as the potential Pb activation and amplification area; the division module 270 is configured to divide the potential Pb activation and amplification area into multiple strip-shaped analysis units perpendicular to the advancement direction of the Pb active front; the calculation module 280 is configured to calculate the Pb pollution risk evolution index of each analysis unit based on the expected arrival time sequence of the Pb active front in each analysis unit and the average Pb exceedance risk value of each unit; and the output module 290 is configured to output an analysis report containing the spatial location of the potential Pb activation and amplification area, the ranking of the Pb pollution risk evolution index of each analysis unit, and the diagnostic results of the key soil factors that dominate the Pb activity changes in each analysis unit.

[0140] It should be understood that Figure 2 The modules and references described in the document Figure 1 The steps described in the text correspond to those in the method described above. Therefore, the operations, features, and corresponding technical effects described above also apply to the method described in the text. Figure 2 The various modules in the document will not be described in detail here.

[0141] In other embodiments, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the program instructions are executed by a processor, the processor performs the method for analyzing lead pollution in paddy field soil in any of the above method embodiments. In one embodiment, the computer-readable storage medium of the present invention stores computer-executable instructions, which are configured as follows: Obtain the current soil property dataset and current soil moisture content distribution data of the target paddy field. The current soil property dataset includes at least the total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity of the soil. Obtain the expected water management parameters for the target paddy field within a future preset time period; The current soil moisture content distribution data and the expected water management parameters are input into a preset soil moisture transport model to simulate the spatiotemporal change sequence of soil moisture content within a preset future time period. Based on the spatiotemporal variation sequence of soil moisture content and the current soil property dataset, calculate the spatial distribution of state indicators characterizing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step; Based on the spatial distribution of the state index, the advancing trajectory of the Pb active front is identified. The Pb active front is the projection of the critical isosurface of the state index from precipitation tendency to dissolution tendency in horizontal space. The advancing trajectory of the Pb active front is overlaid with the spatial distribution of the current Pb exceedance risk in the target paddy field to determine the area that the Pb active front is expected to sweep through and where the current predicted Pb value in the paddy exceeds the safety threshold, as the potential Pb activation and amplification area. Within the potential Pb activation and amplification region, the region is divided into multiple strip-shaped analysis units perpendicular to the advancing direction of the Pb active front. Based on the expected arrival time sequence of the Pb active fronts in each analysis unit and the average Pb exceedance risk value of each unit, the Pb pollution risk evolution index of each analysis unit is calculated. The output includes an analysis report containing the spatial location of the potential Pb activation and amplification region, the ranking of Pb pollution risk evolution indices for each analysis unit, and the diagnostic results of key soil factors that dominate Pb activity changes within each analysis unit.

[0142] Computer-readable storage media may include a stored program area and a stored data area, wherein the stored program area may store an operating system and an application program required for at least one function; the stored data area may store data created based on the use of the paddy field soil lead pollution analysis system, etc. Furthermore, the computer-readable storage medium may include high-speed random access memory, and may also include memory, such as at least one disk storage device, flash memory device, or other non-volatile solid-state storage device. In some embodiments, the computer-readable storage medium may optionally include memory remotely configured relative to a processor, which can be connected to the paddy field soil lead pollution analysis system via a network. Examples of such networks include, but are not limited to, the Internet, corporate intranets, local area networks, mobile communication networks, and combinations thereof.

[0143] Figure 3 This is a schematic diagram of the structure of the electronic device provided in the embodiment of the present invention, such as... Figure 3 As shown, the device includes a processor 310 and a memory 320. The electronic device may also include an input device 330 and an output device 340. The processor 310, memory 320, input device 330, and output device 340 can be connected via a bus or other means. Figure 3 Taking a bus connection as an example, the memory 320 is the computer-readable storage medium described above. The processor 310 executes various server functions and data processing by running non-volatile software programs, instructions, and modules stored in the memory 320, thereby implementing the paddy field soil lead pollution analysis method described in the above embodiment. The input device 330 can receive input digital or character information and generate key signal inputs related to user settings and function control of the paddy field soil lead pollution analysis system. The output device 340 may include a display screen or other display device.

[0144] The aforementioned electronic device can execute the method provided in the embodiments of the present invention, and has the corresponding functional modules and beneficial effects for executing the method. Technical details not described in detail in this embodiment can be found in the method provided in the embodiments of the present invention.

[0145] In one implementation, the above-described electronic device is used in a paddy field soil lead pollution analysis system as a client, comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to: Obtain the current soil property dataset and current soil moisture content distribution data of the target paddy field. The current soil property dataset includes at least the total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity of the soil. Obtain the expected water management parameters for the target paddy field within a future preset time period; The current soil moisture content distribution data and the expected water management parameters are input into a preset soil moisture transport model to simulate the spatiotemporal change sequence of soil moisture content within a preset future time period. Based on the spatiotemporal variation sequence of soil moisture content and the current soil property dataset, calculate the spatial distribution of state indicators characterizing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step; Based on the spatial distribution of the state index, the advancing trajectory of the Pb active front is identified. The Pb active front is the projection of the critical isosurface of the state index from precipitation tendency to dissolution tendency in horizontal space. The advancing trajectory of the Pb active front is overlaid with the spatial distribution of the current Pb exceedance risk in the target paddy field to determine the area that the Pb active front is expected to sweep through and where the current predicted Pb value in the paddy exceeds the safety threshold, as the potential Pb activation and amplification area. Within the potential Pb activation and amplification region, the region is divided into multiple strip-shaped analysis units perpendicular to the advancing direction of the Pb active front. Based on the expected arrival time sequence of the Pb active fronts in each analysis unit and the average Pb exceedance risk value of each unit, the Pb pollution risk evolution index of each analysis unit is calculated. The output includes an analysis report containing the spatial location of the potential Pb activation and amplification region, the ranking of Pb pollution risk evolution indices for each analysis unit, and the diagnostic results of key soil factors that dominate Pb activity changes within each analysis unit.

[0146] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods of various embodiments or some parts of embodiments.

[0147] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for analyzing lead pollution in paddy field soil, characterized in that, include: Obtain the current soil property dataset and current soil moisture content distribution data of the target paddy field. The current soil property dataset includes at least the total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity of the soil. Obtain the expected water management parameters for the target paddy field within a future preset time period; The current soil moisture content distribution data and the expected water management parameters are input into a preset soil moisture transport model to simulate the spatiotemporal change sequence of soil moisture content within a preset future time period. Based on the spatiotemporal variation sequence of soil moisture content and the current soil property dataset, the spatial distribution of state indicators characterizing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step is calculated, including: For each grid point at each time step, obtain the simulated soil moisture content, soil pH, and organic matter content of the grid point; Calculate the soil water-filling porosity based on the simulated soil moisture content; The soil water-filled porosity, soil pH, and organic matter content are input into a preset Pb dissolution-precipitation equilibrium tendency function to calculate the state index values ​​of the grid points at the time step. The preset Pb dissolution-precipitation equilibrium tendency function is a multivariate nonlinear function used to characterize the Pb dissolution-precipitation equilibrium tendency under given water-filled porosity, pH, and organic matter content. 2+ Thermodynamic stability tendency of PbCO3 and Pb(OH)2 precipitates in soil solution relative to solid PbCO3 and Pb(OH)2 precipitates; Based on the spatial distribution of the state index, the advancing trajectory of the Pb active front is identified. The Pb active front is the projection of the critical isosurface of the state index transitioning from precipitation tendency to dissolution tendency onto the horizontal space, specifically including: The state index values ​​of each grid point at a certain time step are compared with a preset state critical threshold, which is the boundary value between precipitation tendency and dissolution tendency. Connect the grid points whose state index values ​​are exactly equal to the state critical threshold to form a dissolution-precipitation boundary line at a certain time step; The dissolution-precipitation boundary lines of all time steps are superimposed in chronological order to form a set of trajectory lines that continuously advance from the initial position into or out of the target paddy field. The set of trajectory lines constitutes the advancing trajectory of the Pb active front. The advancing trajectory of the Pb active front is overlaid with the spatial distribution of the current Pb exceedance risk in the target paddy field to determine the area that the Pb active front is expected to sweep through and where the current predicted Pb value in the paddy exceeds the safety threshold, as the potential Pb activation and amplification area. Within the potential Pb activation and amplification region, the region is divided into multiple strip-shaped analysis units perpendicular to the advancing direction of the Pb active front. Based on the expected arrival time sequence of the Pb active fronts in each analysis unit and the average Pb exceedance risk value of each unit, the Pb pollution risk evolution index of each analysis unit is calculated. The output includes an analysis report containing the spatial location of the potential Pb activation and amplification region, the ranking of Pb pollution risk evolution indices for each analysis unit, and the diagnostic results of key soil factors that dominate Pb activity changes within each analysis unit.

2. The method for analyzing lead pollution in paddy field soil according to claim 1, characterized in that, The acquisition of the current soil attribute dataset and current soil moisture content distribution data of the target paddy field includes: In the target paddy field, multiple monitoring points are set up according to a preset grid, surface soil samples are collected at each monitoring point, and various indicators in the current soil attribute dataset are measured. At each monitoring point, in-situ moisture content sensors were used to obtain the current soil volumetric moisture content data at each point. Spatial interpolation is performed on the current soil volumetric moisture content data at each monitoring point to obtain the current soil moisture content distribution data.

3. The method for analyzing lead pollution in paddy field soil according to claim 1, characterized in that, The expected water management parameters include the timing of irrigation events, the amount of water used in a single irrigation, and the expected flood depth variation curve. The preset soil moisture transport model is a numerical model of unsaturated soil moisture movement based on the Richards equation; The simulated spatiotemporal variation sequence of soil moisture content within a preset future time period includes: Using the current soil moisture content distribution data as the initial condition, the expected flooding depth change curve as the upper boundary condition, and free drainage or a set groundwater level as the lower boundary condition, the Richards equation is iteratively solved at a preset time step to obtain a raster map of the spatial distribution of soil moisture content at each time step.

4. The method for analyzing lead pollution in paddy field soil according to claim 1, characterized in that, The spatial distribution of the current risk of Pb exceeding the standard in the target paddy field was obtained through the following steps: Input the current soil attribute dataset of each monitoring point into the preset rice Pb prediction model to calculate the predicted rice Pb value of each monitoring point. The rice Pb prediction model is a function model that is pre-established based on the multiple linear regression relationship between total soil Pb, pH, organic matter, available Si and cation exchange capacity and rice Pb content. Spatial interpolation was performed on the predicted Pb values ​​of rice at each monitoring point to generate a continuous distribution map of the predicted Pb values ​​of rice. In the continuous distribution map of predicted Pb values ​​in rice, the regions where the predicted Pb values ​​of rice are greater than a preset safety threshold are extracted as the current spatial distribution of Pb exceeding the standard risk in rice.

5. The method for analyzing lead pollution in paddy field soil according to claim 1, characterized in that, The process involves dividing the potential Pb activation and amplification region into multiple strip-shaped analytical units perpendicular to the advancing direction of the Pb active front, along the advancing direction of the Pb active front. Extract the trajectory line of the advancing direction of the Pb active front within the potential Pb activation and amplification region; Multiple dividing points are set along the trajectory line of the propulsion direction with a preset spatial step size; At each segmentation point, a segmentation line perpendicular to the propulsion direction trajectory line is generated, and each segmentation line extends to both sides until it intersects the boundary of the potential Pb activation amplification region; The region enclosed by two adjacent dividing lines and the two side boundaries of the potential Pb activation amplification region is defined as a strip-shaped analysis unit.

6. The method for analyzing lead pollution in paddy field soil according to claim 1, characterized in that, Based on the expected arrival time sequence of the Pb active fronts within each analysis unit and the average Pb exceedance risk value of each unit, the Pb pollution risk evolution index of each analysis unit is calculated, including: Calculate the arithmetic mean of the predicted Pb values ​​of rice for all grid points within a certain analysis unit, and use it as the current average risk value for that analysis unit; Based on the advancing trajectory of the Pb active front, the estimated time when a certain analytical unit is swept by the Pb active front is determined; The current average risk value is weighted and combined with the expected time to calculate the Pb pollution risk evolution index.

7. The method for analyzing lead pollution in paddy field soil according to claim 1, characterized in that, Output an analytical report containing diagnostic results of key soil factors dominating Pb activity changes within each analysis unit, including: For each analysis unit, the correlation coefficients between soil pH, organic matter content, available Si content, cation exchange capacity, and total Pb content and the state indicators within the analysis unit were calculated. The soil factor with the largest absolute value of the correlation coefficient was identified as the dominant factor of the analysis unit; The analysis report presents the ranking of Pb pollution risk evolution indices for each analysis unit in the form of a spatial distribution map, and lists the dominant factors and the average deviation of the dominant factors for each analysis unit in the form of tables or annotations.

8. A system for analyzing lead pollution in paddy field soil, characterized in that, include: The first acquisition module is configured to acquire the current soil attribute dataset and the current soil moisture content distribution data of the target paddy field. The current soil attribute dataset includes at least the total Pb content, pH value, organic matter content, available Si content, and cation exchange capacity of the soil. The second acquisition module is configured to acquire the expected water management parameters of the target paddy field in a future preset period. The simulation module is configured to input the current soil moisture content distribution data and the expected water management parameters into a preset soil moisture transport model to simulate the spatiotemporal change sequence of soil moisture content within a preset future time period. The calculation module is configured to calculate the spatial distribution of state indicators characterizing the tendency of Pb dissolution-precipitation equilibrium in the soil at each time step, based on the spatiotemporal variation sequence of soil moisture content and the current soil attribute dataset, including: For each grid point at each time step, obtain the simulated soil moisture content, soil pH, and organic matter content of the grid point; Calculate the soil water-filling porosity based on the simulated soil moisture content; The soil water-filled porosity, soil pH, and organic matter content are input into a preset Pb dissolution-precipitation equilibrium tendency function to calculate the state index values ​​of the grid points at the time step. The preset Pb dissolution-precipitation equilibrium tendency function is a multivariate nonlinear function used to characterize the Pb dissolution-precipitation equilibrium tendency under given water-filled porosity, pH, and organic matter content. 2+ Thermodynamic stability tendency of PbCO3 and Pb(OH)2 precipitates in soil solution relative to solid PbCO3 and Pb(OH)2 precipitates; The identification module is configured to identify the propagation trajectory of the Pb active front based on the spatial distribution of the state index. The Pb active front is the projection of the critical isosurface of the state index transitioning from precipitation tendency to dissolution tendency onto a horizontal surface in space, specifically including: The state index values ​​of each grid point at a certain time step are compared with a preset state critical threshold, which is the boundary value between precipitation tendency and dissolution tendency. Connect the grid points whose state index values ​​are exactly equal to the state critical threshold to form a dissolution-precipitation boundary line at a certain time step; The dissolution-precipitation boundary lines of all time steps are superimposed in chronological order to form a set of trajectory lines that continuously advance from the initial position into or out of the target paddy field. The set of trajectory lines constitutes the advancing trajectory of the Pb active front. The determination module is configured to overlay the advancement trajectory of the Pb active front with the spatial distribution of the current Pb exceedance risk of the target paddy field to determine the area that the Pb active front is expected to sweep through and the current predicted Pb value of the paddy exceeds the safety threshold, as the potential Pb activation and amplification area. The partitioning module is configured to divide the potential Pb activation and amplification region into multiple strip-shaped analysis units perpendicular to the advancing direction of the Pb active front. The calculation module is configured to calculate the Pb pollution risk evolution index of each analysis unit based on the expected arrival time sequence of the Pb active fronts in each analysis unit and the average Pb exceedance risk value of each unit. The output module is configured to output an analysis report containing the spatial location of the potential Pb activation and amplification region, the ranking of Pb pollution risk evolution indices for each analysis unit, and the diagnostic results of key soil factors that dominate Pb activity changes within each analysis unit.

Images