Soil remediation target dynamic formulation method based on protection of underground water conservation function

Abstract

The invention relates to the technical field of soil remediation, in particular to a soil remediation target dynamic formulation method based on an underground water conservation protection function, which comprises the following steps: acquiring a remediation starting point and water level data, identifying a lagging time period, adjusting a boundary to generate a dynamic mapsheet, and extracting a water level gradient and seepage direction sorting intervention sequence. Fluctuation and stagnant water capacity are judged to mark a regulation and control section, and concentration and water level trends are superposed to set a restoration time history table. According to the method, an intervention lag relationship is determined by extracting a water level inflection point and an infiltration difference value, a boundary response direction is adjusted according to a lag time period, the adaptability of a control range is enhanced, an intervention priority sequence is constructed by combining a hydraulic gradient and a seepage trend, the repair efficiency of a continuous segment is highlighted, and a restricted area is judged by utilizing water layer fluctuation and water retention capacity. The concentration regulation and control position is accurately delimited, the water level trend and pollution distribution are superimposed, the repair frequency and time history are set, and dynamic cooperation of intervention recognition, boundary adjustment and the repair rhythm is achieved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of soil remediation, in particular to a soil remediation target dynamic formulation method based on protection of groundwater source conservation function. BACKGROUND

[0002] The technical field of soil remediation involves the governance and restoration of contaminated soil to reduce or eliminate the risk of harmful substances to the ecological environment and human health, which is a key component of environmental engineering. The core matters of this technical field include identification and assessment of contaminated soil, setting of remediation targets, selection and implementation of remediation technology, and post-monitoring and management, etc. It usually covers physical remediation, chemical remediation, biological remediation and other means, and a comprehensive remediation plan is developed according to the type of pollutants, site characteristics and environmental function requirements. In particular, in the protection of regional water and soil ecosystems, the coordination of soil remediation and groundwater conservation function is particularly critical, and has become an important direction of soil environmental governance in recent years. Among them, the traditional soil remediation target dynamic formulation method refers to formulating remediation targets according to the type and degree of pollution of contaminated soil and dynamically adjusting the targets during the implementation process to adapt to changing environmental conditions and remediation progress. The technical matter targeted is how to scientifically set and dynamically adjust the soil remediation target to effectively protect the groundwater conservation function. Traditionally, the method of revising the remediation target regularly by combining pollutant concentration detection data, hydrogeological investigation results and national or local soil environmental quality standards, pollutant migration trend simulation results and risk assessment parameters is commonly used. Common means include regional groundwater recharge intensity-based plot zoning strategy, pollution source response analysis under hydrological boundary conditions, and reference to experience model parameters of similar geological conditions to develop remediation path and threshold standard.

[0003] The existing technology relies on the static data structure of the stage detection results and the geological investigation parameters in the adjustment process of the remediation target, and fails to build a continuous response discrimination system. In the case of frequent dynamic changes of groundwater or unclear human intervention points, it is difficult to determine whether the intervention behavior produces effective response in time, and there are problems of target adjustment lag and control interface solidification. For example, when the hydrological boundary changes, the original remediation area range lacks dynamic correction ability, causing the remediation intervention time sequence and the pollution diffusion path to be dislocated, affecting the accurate allocation of remediation resources and the coordinated maintenance of regional hydrological function. SUMMARY

[0004] In order to solve the technical problems existing in the prior art, the present application provides a soil remediation target dynamic formulation method based on protection of groundwater source conservation function, comprising the following steps: In order to achieve the above purpose, the present application adopts the following technical scheme: the soil remediation target dynamic formulation method based on protection of groundwater source conservation function comprises the following steps: S1: Obtain groundwater source conservation area repair starting point and water level data, extract three-day water level rising trend and locate change inflection point, record interval between repair time point, detect evaporation and infiltration rate difference, obtain repair intervention effective lag period; S2: Call the repair intervention effective lag period, select boundary expansion line, judge whether the period is shorter than the response time limit, if short, move the boundary out, if long, freeze the boundary and mark the core area, obtain the dynamic boundary map of the pollution control area; S3: Call the dynamic boundary map of the pollution control area, extract water level gradient change and seepage direction consistency, screen water power continuous area and sort intervention order, obtain repair intervention priority sequence list; S4: Call the repair intervention priority sequence list, extract the difference between the upper and lower boundaries of the capillary water layer and calculate the fluctuation, detect the particle size and organic matter content, compare the fluctuation and the limit value of the lag water capacity, obtain the intervention concentration control section marker map; S5: Call the intervention concentration control section marker map, extract the pollution distribution map and mark the high value section, superimpose the underground water level rising trend map, judge the consistency of the water content rate trend, obtain the dynamic soil repair target development scheme.

[0005] As a further scheme of the application, the repair intervention effective lag period includes water level change inflection point time sequence, lag interval length, rate difference amplitude, lag type classification, the dynamic boundary map of the pollution control area includes boundary adjustment state, core reserved area position, map revision range, response time limit comparison result, the repair intervention priority sequence list includes gradient change frequency, change rate sorting, seepage consistency ratio, water power continuous segment, intervention priority level, the intervention concentration control section marker map includes fluctuation amplitude index, particle size proportion estimate, organic matter content level, lag water capacity judgment result, limit section identification, and the dynamic soil repair target development scheme includes concentration high value section distribution, water content rate trend direction, repair frequency setting, time cycle arrangement.

[0006] As a further scheme of the application, the response time limit refers to the shortest time threshold required for hydrological changes to respond to repair intervention measures after boundary expansion.

[0007] As a further scheme of the application, the screening of water power continuous area and the sorting of intervention order refers to identifying water power coherent segments according to water level gradient change frequency and seepage direction consistency, and determining intervention priority according to water level mutation level.

[0008] As a further scheme of the application, the specific steps of S1 are: S101: Obtain the groundwater level monitoring data frame, extract the water level sequence of three consecutive days, calculate the daily increment, screen the time period of continuous water level rise, pair the end date with the corresponding water level, and obtain the continuous rising water level marker set; S102: Call the continuous rising water level marker set, calculate the increment change rate, locate the time point corresponding to the maximum increment, judge the change of the direction of the front and rear increments, screen the time of trend turning point, and obtain the water level change inflection point time node; S103: Call the water level change inflection point time node and the repair starting time stamp, calculate the time difference between the two, eliminate the inflection points earlier than the repair starting point, retain the time difference value and aggregate, obtain the continuous rising water level marker set, and detect the difference between the evaporation rate and the infiltration rate, and generate the repair intervention lag period.

[0009] As a further scheme of the present application, the specific steps of S2 are: S201: Call the repair intervention lag period, select the original extension line position of the pollution control area boundary, perform difference judgment on the lag period and the boundary response time limit, screen the boundary section whose lag period is shorter than the response time limit, and obtain the boundary extension adjustment section set; S202: According to the boundary extension adjustment section set, perform outward moving operation on the corresponding boundary, retain the boundary whose lag period is longer than the response time limit and demarcate the core reserved area, integrate the boundary change, and generate the repair boundary adjustment layer; S203: Call the repair boundary adjustment layer, superimpose the original boundary graph, identify the boundary change area, redraw the overall boundary line and label the area attribute, extract the layer shape data, and generate the pollution control area dynamic boundary map.

[0010] As a further scheme of the present application, the specific steps of S3 are: S301: Call the pollution control area dynamic boundary map, extract the water level data sequence between monitoring points, calculate the water level difference between adjacent monitoring points and count the positive and negative direction change frequency, and simultaneously calculate the change rate of the absolute value of the continuous difference value, aggregate the frequency and rate data, and obtain the water level gradient change index set; S302: According to the water level gradient change index set, extract the seepage direction scalar value of adjacent monitoring segments, judge the direction consistency between adjacent segments, aggregate the consistent segments and arrange them in order, extract the continuous segment set, and generate the water power continuity area list; S303: Call the water power continuity area list, extract the water level mutation difference value in the area, mark the difference value domain by level, sort the area number from high to low according to the level, establish the area and number mapping table and record the sorting position, and establish the repair intervention priority sequence list.

[0011] As a further scheme of the present application, the specific steps of S4 are: S401: Call the calibrated segment in the repair intervention priority sequence list, extract the upper and lower boundary elevation values of the corresponding capillary water layer, calculate the upper and lower boundary difference value, and calculate the difference between the maximum and minimum amplitude of the continuous difference value sequence to obtain the capillary water layer fluctuation amplitude value; S402: According to the capillary water layer fluctuation amplitude value, extract the soil particle size distribution and organic matter content data frame of the same segment, calculate the coarse grain component proportion, and estimate the water storage capacity index combined with the organic matter content. Compare the index with the water storage retention capacity threshold to generate a set of restricted area identification results; S403: Call the restricted area identification result set, locate the corresponding segment position, mark the area range in the pollution control zone layer, and establish a mapping combined with the segment number and control zone segment number. Output the vector boundary and label attribute field to generate the intervention concentration control segment label map.

[0012] As a further scheme of the present application, the specific steps of S5 are: S501: Call the area identified by the intervention concentration control segment label map, extract the pollutant concentration sequence in the soil layer depth direction of the corresponding area, identify the concentration peak position and mark the space, construct the atlas structure combined with the profile depth data, and generate the vertical pollution concentration distribution atlas; S502: Call the vertical pollution concentration distribution atlas, superimpose the groundwater level recovery trend change graph, locate the concentration change segment in the atlas and extract the corresponding spatial position, and judge whether the trend direction of the position is consistent with the pollution trend direction to obtain the trend matching relationship determination result; S503: According to the trend matching relationship determination result, set the corresponding repair period for the trend consistent area and the inconsistent area respectively, establish the repair frequency and time corresponding table according to the spatial segment number, arrange the intervention segment operation in time sequence, and establish the dynamic development scheme of soil repair target.

[0013] As a further scheme of the present application, the water storage retention capacity threshold is used to determine the minimum water storage performance standard value of whether the soil layer has water storage capacity.

[0014] Compared with the prior art, the present application has the following advantages and positive effects: In the present application, the intervention lag relationship is determined by extracting the water level inflection point and the infiltration difference, the boundary response direction is adjusted according to the lag period, the adaptability of the control range is enhanced, the intervention priority order is constructed combined with the hydraulic gradient and seepage trend, the repair efficiency of continuous segments is highlighted, the limited area is determined by water layer fluctuation and water retention capacity, the concentration control position is accurately determined, the repair frequency and time are set by superimposing the water level trend and pollution distribution, and the dynamic cooperation of intervention identification, boundary adjustment and repair rhythm is realized. Attached Figure Description

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

[0016] Figure 1 This is a schematic diagram of the steps of the present invention; Figure 2 This is a detailed schematic diagram of S1 of the present invention; Figure 3 This is a detailed schematic diagram of S2 of the present invention; Figure 4 This is a detailed schematic diagram of S3 of the present invention; Figure 5 This is a detailed schematic diagram of S4 of the present invention; Figure 6 This is a detailed schematic diagram of S5 of the present invention. Detailed Implementation

[0017] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0018] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0019] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning.

[0020] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0021] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0022] Please refer to Figure 1 The embodiment of the present application provides a soil remediation target dynamic formulation method based on groundwater source conservation function, comprising the following steps: S1: Obtain the groundwater source conservation area repair starting point and groundwater level monitoring data, extract the three-day water level rising trend and locate the maximum change inflection point, record the time interval between the inflection point and the repair time point, detect the difference between the evaporation amount and the vertical infiltration rate at the same period, compare the change amplitude between the time interval and the rate difference, divide the lag relationship between the intervention behavior and the hydrological response, and obtain the repair intervention effective lag period; S2: Call the repair intervention effective lag period, select the original set of pollution control area boundary extension line position, judge whether the lag period is shorter than the boundary response time limit, if it is short, move the repair boundary outward, if it is long, freeze the boundary unchanged and mark the core reserved area, correct the repair processing area corresponding to the regional boundary, and obtain the pollution control area dynamic boundary map; S3: Call the pollution control area dynamic boundary map, extract the water level gradient change frequency and change rate between monitoring points, calculate the consistency ratio of adjacent segments in the seepage direction and determine the trend direction, screen the hydrodynamic continuity area and sort the intervention order according to the water level mutation grade, establish the intervention operation sequence, and obtain the repair intervention priority sequence list; S4: Call the marked segment in the repair intervention priority sequence list, extract the difference between the upper and lower boundaries of the capillary water layer and calculate the fluctuation amplitude, detect the lag water retention capacity of the particle size proportion and organic matter content in the same segment, compare the fluctuation amplitude with the retention capacity threshold and determine the restricted area, and obtain the intervention concentration control section marker map; S5: Call the intervention concentration control section marker map identified area, extract the vertical distribution map of pollutants in the soil layer and mark the high concentration section, superimpose the groundwater level rising trend change map, determine whether the trend direction of the water content in the overlapping position is consistent, set the repair frequency and time according to the trend matching situation, and obtain the soil remediation target dynamic formulation scheme.

[0023] The repair intervention effective lag period includes the water level change inflection point time sequence, the lag interval length, the rate difference amplitude, and the lag type classification. The pollution control area dynamic boundary map includes the boundary adjustment state, the core reserved area position, the map correction range, and the response time limit comparison result. The repair intervention priority sequence list includes the gradient change frequency, the change rate sorting, the seepage consistency ratio, the hydrodynamic continuity segment, and the intervention priority grade. The intervention concentration control section marker map includes the fluctuation amplitude index, the particle size proportion estimate, the organic matter content level, the lag water capacity judgment result, and the restricted section identification. The soil remediation target dynamic formulation scheme includes the high concentration section distribution, the water content trend direction, the repair frequency setting, and the time cycle arrangement.

[0024] Please refer to Figure 2 The specific steps of S1 are: S101: Obtain the groundwater level monitoring data frame, extract the continuous three-day water level sequence, calculate the daily increment, screen the time period of continuous water level rise, pair the endpoint date with the corresponding water level, and obtain the continuous rising water level marker set; The groundwater level monitoring data frame is usually obtained by deploying automatic water level monitoring equipment in underground water observation wells for regular collection. Each record includes a timestamp, a water level value, and a site number. The data frame structure is arranged in ascending order of timestamp to ensure sequence continuity. Then, a time window of 3 days is divided, and a sliding window method with a step size of 1 day is used to extract the continuous three-day water level sequence. The water level values of the three days are recorded for each three-day sequence. The water level changes of the first day and the second day, and the second day and the third day are calculated. The water level difference between the first day and the second day, and the water level difference between the second day and the third day are extracted. The increment calculation is completed by direct subtraction. This calculation operation is performed in groups. For example, at a certain observation point, the first day water level is 3.12 meters, the second day is 3.24 meters, and the third day is 3.41 meters. The daily increment is 0.12 meters and 0.17 meters, respectively. When both increments are positive, it indicates that the water level is rising for two consecutive days. Further, the timestamp and water level value of the third day are combined into a group as the continuous rising endpoint record in the current time window. After continuously traversing the entire monitoring data, multiple continuous rising water level time nodes and water level value pairing results can be formed. These paired values constitute the continuous rising water level marker set, which is used for subsequent judgment of the turning point time distribution of water level change trend.

[0025] S102: Call the continuous rising water level marker set, calculate the increment change rate, locate the time point corresponding to the maximum increment, judge the change of the direction of the previous and subsequent increments, screen the time of trend turning point, and obtain the water level change turning point time node; After arranging the continuous rising water level mark set in chronological order, the water level values and time stamps of the adjacent two records are extracted, the water level increment between the two records is calculated, the current increment value is obtained, and then the adjacent two increments are compared to determine whether the increment has an upward or downward trend. By continuously judging the difference trend of adjacent increments, it can be determined whether the increment change shows a decreasing or increasing trend. For example, in the records of a certain observation station, the first continuous rising water level is 3.50 meters, the second is 3.68 meters, and the interval is two days, with an increment of 0.18 meters. The next one is 3.79 meters, with an interval of one day, and the increment is 0.11 meters. The increment decreases from 0.18 to 0.11, so it can be determined that the increment has a downward trend. By detecting this trend change process, the time node where the maximum increment value is located is locked, and the increment records before and after it are compared. If it changes from continuous growth to decrease, it means that the turning point of water level rising trend occurs here. The necessary condition for judging the turning point is that the direction of the increment change changes, i.e. from increasing to decreasing or reverse change. A fixed interval increment comparison method is used to identify the trend inflection point. In order to avoid misjudgment, a minimum change threshold is set. It is recommended to set the change amplitude not less than 0.05 meters. Changes less than this amplitude will be considered as invalid fluctuations. In actual operation, the value can be set according to historical data statistics. For example, if the standard deviation of annual water level change is 0.12 meters, 0.05 meters can be set as the lower limit of effective trend change. Finally, the time nodes that meet the above change characteristics are extracted to form the water level change inflection point time nodes.

[0026] S103: Call the water level change inflection point time node and the repair starting time stamp, calculate the time difference between the two, eliminate the inflection points earlier than the repair starting point, retain the time difference value and aggregate, obtain the continuous rising water level mark set, and detect the difference between the evaporation rate and the infiltration rate to generate the repair intervention lag period of effectiveness; With the water level change inflection point time node set as the basis, compare each inflection point time with the repair intervention start time stamp, calculate the time interval between the two, if the inflection point time is earlier than the repair start time, it is judged as a non-intervention related inflection point and is eliminated, the remaining inflection point time is the response time node potentially affected by the repair behavior, for each inflection point time that meets the condition, calculate the specific interval days between it and the repair start point, and record it as the basic data of the lag response, for example, the repair start time is May 1st, and the inflection point is May 3rd, then the corresponding lag days are 2 days, if the inflection point is April 29th, it is discarded because it is earlier than the intervention start point, in order to avoid calculation errors, the date needs to be standardized to ensure that the interval calculation is in a unified time format, further, all lag time points that meet the requirements are collected to form a lag period set, in order to ensure that the identified lag response is significant, a lag threshold range can be set, for example, only identify the lag inflection points within 1 to 10 days, and those more than 10 days are considered as non-direct response relationship and are eliminated, the setting of this threshold can refer to the usual effective time of the intervention measure, combined with the setting of the historical response record, if the average response time is 7 days and the maximum response time is 10 days, the lag analysis interval can be set between 1 to 10 days, finally output the repair intervention effective lag period as the basic input data for subsequent intervention lag effect analysis.

[0027] Please refer to Figure 3 , the specific steps of S2 are: S201: Call the repair intervention effective lag period, select the original set of pollution control area boundary expansion line positions, difference judge the lag period and the boundary response time limit, select the boundary section whose lag period is shorter than the response time limit, and obtain the boundary extension adjustment section set; When calling the lag period during which the repair intervention takes effect, it is necessary to traverse the data of the lag days, extract the spatial coordinate information corresponding to each lag period one by one, ensure that this period can be accurately mapped to the boundary area of the pollution control zone geographically, select the position of the originally set extended line in the boundary layer. The originally set extended line is a buffer line segment reserved for the outer edge of the initial pollution control range. Each extended line segment has a unique number and its corresponding spatial coordinate information, and the corresponding boundary response time limit is marked. The boundary response time limit refers to the theoretical minimum time for the pollutant impact to spread to this boundary segment set during the original assessment process. This response time limit can be determined comprehensively according to parameters such as the source release intensity, regional groundwater flow velocity, and formation permeability. It is recommended to set the interval to 6 to 10 days. Among them, when more than 80% of the data in the frequency distribution of the historical response time of a certain region is concentrated between 6 and 8 days, the upper bound can be taken as the response threshold. For example, when setting the boundary response time limit to 8 days, judge the difference between each lag period T_delay and the response time limit T_response of the boundary segment where it is located. When T_delay is less than T_response, that is, T_delay < T_response, it means that the pollution water level reaction speed at this point is relatively fast and has not reached the response time that the boundary segment should have, which is regarded as the boundary segment not fully responding to the pollution diffusion and meets the conditions for extension adjustment. Otherwise, it is retained unchanged. For example, if the lag time of a certain point is 5 days and the response time limit of the boundary segment where it is located is 8 days, the difference is days, meeting the adjustment conditions; if the lag time is 9 days, it does not meet the screening conditions and is discarded. After repeating this process to traverse all lag points and boundary segments, finally, all boundary segments that meet T_delay < T_response are extracted to generate a boundary extension adjustment section set.

[0028] S202: According to the boundary extension adjustment section set, perform an outward movement operation on the corresponding boundary, retain the boundary with a lag period longer than the response time limit and mark the core retention area, integrate the boundary changes, and generate a repair boundary adjustment layer; According to the number and spatial position provided by the boundary extension adjustment section set, each boundary section is subjected to an outward operation, and the outward direction should be along the pollution source diffusion trend direction. Generally, the boundary normal direction is taken for extension. The outward distance needs to comprehensively consider the actual migration ability of the pollutant, the groundwater flow rate, the terrain slope and other factors. It is recommended to set the offset interval to 50 meters to 300 meters. If in a high permeability area, it is recommended to take the upper limit of the offset value. If in a low permeability clay area, the lower limit can be taken to ensure that the boundary extrapolation meets the geological conditions. For example, if a high permeability section is set to an outward distance of 250 meters, the original boundary coordinate point is translated along the normal direction by 250 meters to generate a new boundary vertex. After constructing a new boundary section, the next one is processed. During the processing, it is judged whether the lag time of the boundary section exceeds the response time limit. If the lag time is greater than or equal to the response time limit, the boundary section is not adjusted, and is directly marked as a core reserved area. The identification field "reserved area = yes" is set in the attribute table. In addition, in order to facilitate subsequent layer integration, the "adjustment area = yes" field is added to all boundary sections subjected to the outward operation. After completing all boundary section operations, the boundary sections subjected to the outward operation and the boundary sections not changed are combined. According to the connection relationship, a new continuous boundary line is reconstructed to realize seamless splicing of the layer boundary line sections. The overall result is output as a repaired boundary adjustment layer.

[0029] S203: Call the repaired boundary adjustment layer, superimpose the original boundary map, identify the boundary change area, redraw the overall boundary line and mark the area attribute, extract the layer shape data, and generate a pollution control area dynamic boundary map; After calling the repaired boundary adjustment layer, it needs to be spatially superimposed with the original boundary map. Through spatial joint analysis of the two layers on the GIS platform, the boundary line sections with non-coinciding coordinates in the two layers are identified as changed areas. These changed sections and adjacent unchanged sections are further spliced to construct a complete new boundary line in order. The redrawn boundary line needs to be displayed in different colors or line types on the graphical interface to facilitate manual identification. At the same time, the attributes of all boundary sections are marked, including "whether to adjust", "whether to reserve", "offset distance", "adjustment direction" and other fields. These attribute data are bound with the layer spatial data and stored in the boundary layer data table. When extracting the layer shape data, the start and end point coordinate information of each boundary section needs to be extracted to ensure that the original boundary section geometry can be restored. Finally, all the marked boundary lines are combined and output as a pollution control area dynamic boundary map. The map should be in a standard GIS vector data format, such as shp or geojson format, and contain metadata information such as layer timestamp and version number. Each boundary section in the map contains complete attribute fields and geometric information, which can be used for subsequent supervision, analysis and display.

[0030] Please refer to Figure 4 , and the specific steps of S3 are as follows: S301: Call the pollution control area dynamic boundary map, extract the water level data sequence between monitoring points, calculate the water level difference between adjacent monitoring points and count the positive and negative direction change frequency, and calculate the change rate of the absolute value of the continuous difference value, aggregate the frequency and rate data, and obtain the water level gradient change index set; First, the geographic coordinates and corresponding water level time series data of all monitoring points within the coverage of the map need to be extracted. Each monitoring point should have a uniform time step, such as daily or hourly, to ensure comparability in the time dimension. The extracted data is constructed into a two-dimensional sequence matrix, with rows representing time and columns representing different monitoring points. On this basis, pairs of adjacent monitoring points are selected according to the principle that the spatial distance is not more than 200 meters or there is an underground water channel along the actual connection line. For example, if monitoring point P1 is 120 meters away from P2, they form a pair of effective comparison units. For each pair of monitoring points, the water level difference is calculated at each time step, and the water level of the downstream point is subtracted from the water level of the upstream point to form a difference sequence. The positive and negative directions of the difference value are recorded. If the difference value is positive, it indicates that the water flow direction is P1 to P2, and if it is negative, the direction is opposite. By traversing all time steps, the frequency of positive and negative directions for each pair of monitoring points is recorded. If the positive direction appears 8 times and the negative direction appears 2 times, the positive frequency is recorded as 80% and the negative frequency is recorded as 20%. The absolute value of the difference value at each time step is subtracted from the absolute value of the difference value at the adjacent time step to obtain the water level change rate. The water level change rate is the change amplitude of the absolute value of the water level difference between adjacent time steps. All rate values are averaged, maximized, and standardized to aggregate the characteristic values of the change rate of the monitoring point pair. Finally, the positive and negative direction frequencies and change rate characteristics of all monitoring point pairs are integrated and summarized, and the output is the water level gradient change index set.

[0031] S302: According to the water level gradient change index set, extract the seepage direction scalar value of adjacent monitoring segments, judge the direction consistency between adjacent segments, aggregate and arrange the consistent segments in order, extract the continuous segment set, and generate the water power continuity area list; According to the water level difference and direction change information of each pair of monitoring points extracted from the water level gradient change index set, the seepage direction scalar value of the adjacent point pairs constituting the monitoring segment is induced. The seepage direction scalar value is set according to the water level sequence of the monitoring points. If the water level of the upstream point is greater than that of the downstream point, the direction scalar is positive, otherwise it is negative. The direction scalars of all monitoring segments are compared. If the direction scalars of adjacent segments are consistent in sign, it is judged that the directions are consistent. If the scalars of three adjacent segments are all positive, it is considered that the direction of the segment is continuously consistent, and the segment is aggregated into a continuous segment. If the direction changes suddenly at a certain place, such as the scalar changes from positive to negative, it is considered as a discontinuity and is not included in the current aggregated segment. A new continuous segment is constructed from the beginning. For example, if the direction of segment A-B is positive and the direction of B-C is also positive, A-B-C constitutes a continuous segment. If the direction of C-D changes to negative, a new segment starts from D. The aggregation of the direction continuous segment is completed by traversing the entire index set. The start point, end point and intermediate monitoring points of each aggregated segment are recorded. The segments are arranged in spatial connection order to form a traceable path segment. The segment set that meets the continuous direction condition is integrated to generate a water dynamic continuity area list. The list should include the spatial position, length, start and end monitoring point number and direction scalar of each continuous segment.

[0032] S303: Call the water dynamic continuity area list, extract the water level mutation difference in the area, mark the difference value range by level, sort the area number from high to low, establish the area and number mapping table and record the sorting position, and establish the repair intervention priority sequence list; After calling the water dynamic continuity area list, the water level time series data in each continuity area needs to be extracted again, especially whether there is a mutation difference in the continuous segment, that is, the water level difference of two adjacent monitoring points at a certain time step is much higher than the average value of the area. For this, the average value and standard deviation of the water level difference in the area need to be calculated first, and then the water level mutation point is screened. The threshold value of mutation determination is recommended to be set as the average difference value of the area plus twice the standard deviation. For example, if the average value of the water level difference in a continuous segment is 0.2 meters and the standard deviation is 0.05 meters, the threshold value is set to 0.3 meters. The point pair exceeding the value is considered as a mutation point. After extracting all the mutation points, the difference value range is processed by level. The division level is recommended to be grouped according to 0.3-0.5 meters, 0.5-0.8 meters and 0.8 meters and above. The corresponding labels are first, second and third mutation levels. Each level corresponds to different intervention intensity. According to the mutation level from high to low, the area is numbered and sorted. The number priority can be set by multiplying the mutation level by the length weighting method, that is, the area with long length and severe mutation is sorted earlier. The sorting result is bound to the corresponding area to establish a one-to-one mapping table of the area and the number, and the position of each area in the whole sorting is recorded. Finally, it is output as a repair intervention priority sequence list, which can be used for subsequent repair task scheduling.

[0033] Referring to Figure 5 The specific steps of S4 are as follows: S401: Call the calibrated segment in the repair intervention priority sequence list, extract the upper and lower boundary elevation values of the capillary water layer in the corresponding area, calculate the difference between the upper and lower boundary values, and calculate the difference between the maximum amplitude and the minimum amplitude of the continuous difference value sequence to obtain the fluctuation amplitude value of the capillary water layer; When calling the calibrated segment in the repair intervention priority sequence list, the spatial number of the corresponding segment should be retrieved one by one in the area map or database, and the capillary water layer upper and lower boundary elevation value data within the coverage of the segment should be extracted as the retrieval condition. This data comes from borehole stratification records or geological radar profile data, and is generally recorded in segments with a spacing of one meter to ensure that the upper and lower boundary elevations are consistent in units and reference surfaces. Then, the upper and lower boundary elevations at all sampling points are subjected to one-to-one difference operation to calculate the capillary water layer thickness at each point. Then, according to time or space sorting, a continuous thickness difference sequence is constructed. For this sequence, the maximum and minimum values are extracted, which represent the maximum and minimum thicknesses of the capillary water layer in the time or space scale of the segment, respectively. Then, the difference between the two values is calculated to obtain the fluctuation amplitude value of the capillary water layer, which reflects the degree of vertical variation of the capillary water layer within a certain period or within a certain section. For example, if the upper boundary elevation range of a segment is 63.0 to 64.2 meters and the lower boundary is 60.5 to 62.8 meters, the minimum thickness is 1.2 meters, the maximum thickness is 3.0 meters, and the fluctuation amplitude value is 1.8 meters. This provides key basic parameters for subsequent limit area identification and control strategy formulation.

[0034] S402: According to the fluctuation amplitude value of the capillary water layer, extract the soil particle size distribution and organic matter content data frame of the same segment, calculate the proportion of coarse-grained components, and estimate the water storage capacity index combined with the organic matter content. Compare the index with the water storage retention capacity threshold to generate a limit area identification result set; Extract the soil particle size distribution data and organic matter content data frames in the spatial range of the same section. The data comes from the in-situ sampling particle size analysis report or the soil test result file. The particle size distribution data is usually expressed in percentage of grading mass, covering coarse sand, medium sand, fine sand, silt, clay and other particle grades. The coarse particle component ratio is obtained by adding the proportions of coarse sand and medium sand. Combined with the organic matter content data of the region, the organic matter content is generally expressed in percentage of total organic carbon. The two are combined to estimate the water storage capacity index. The estimation process is to multiply the coarse particle ratio by the organic matter content to obtain a dimensionless value. For example, the coarse particle ratio of a certain section is 55%, and the organic matter content is 2.5%. The water storage capacity index is 1.375. To make a judgment, the water storage capacity index needs to be compared with the water storage retention capacity threshold. The threshold value should be set by statistical analysis of the historical soil water storage performance of the pollution control area. It is recommended to select typical representative values from historical data and set them by percentile. For example, if the 25th percentile water storage capacity index is 1.5, then 1.5 can be set as the threshold, indicating that sections below 1.5 generally have weak water storage capacity and need to be included in the restricted range. The index of all sections is iterated to determine whether it is below the threshold. The section number and attribute information of the records below the threshold are extracted to form the restricted area identification result set, which includes section number, index value, and whether it is below the threshold.

[0035] S403: Call the restricted area identification result set, locate the corresponding section position, mark the area range in the pollution control area layer, and establish a mapping between the section number and the control area section number. Output the vector boundary and label attribute field to generate the intervention concentration control section label map. After calling the restricted area identification result set, the spatial range of the section should be located in the pollution control area layer, and its spatial geometry boundary in the layer should be determined. The extracted spatial range exists in the form of polygon or polyline, and the coordinate system of the layer should be consistent to ensure accurate projection. At the same time, write the restriction label field value into the layer attribute table, and set the label value to "1" to indicate that the area is identified as a restricted section. Sort all marked areas by section number and map them with the original boundary section number of the control area. The mapping method is to correspond each restricted section to its corresponding pollution control section, for example, a restricted section number L-P009 corresponds to a control section number Z3-7, which is recorded in the mapping table as After the mapping table is established, it is attached to the layer attribute table as an external association field, and finally a layer file containing vector boundary graphics and attribute data is formed. The layer attribute field should at least contain the fields of segment number, restriction mark, belonging control section number, storage capacity index value, etc. The output format is a standard GIS supported vector data format such as SHP or GeoJSON, and can be used for spatial calling and area identification in the downstream intervention area regulation strategy module, and finally an intervention concentration regulation section mark map is generated.

[0036] Please refer to Figure 6 The specific steps of S5 are: S501: Call the identified area of the intervention concentration regulation section mark map, extract the pollutant concentration sequence in the corresponding area in the soil depth direction, identify the concentration peak position and mark it in space, construct the atlas structure combined with the profile depth data, and generate the vertical pollution concentration distribution atlas; When calling the identified area of the intervention concentration regulation section mark map, the vector boundary of the marked intervention area in geographic space should be located, the pollution monitoring drill hole or sample point number within the boundary should be extracted, and the soil depth direction pollutant concentration sequence of the corresponding drill hole should be retrieved accordingly. The sequence records the pollutant content value at different horizons from the ground surface to the survey final hole depth, usually in mg / kg or mg / L. The concentration value is recorded at intervals of every 0.5 meters or 1 meter. A depth-concentration one-dimensional data curve is constructed. The maximum value of each point on the curve is extracted as the concentration peak position of the drill hole. The depth value and corresponding concentration of the position are recorded, and the position is marked in space. The horizontal coordinates of the drill hole and the peak depth are combined to form a three-dimensional spatial coordinate, and a visual point set is generated. After extracting the peak positions of all drill holes, further superimpose the geological profile data of each drill hole to obtain the burial depth range, lithology layering structure and other information of each layer in the profile. The relative position of the concentration peak in the profile is located, and the three-dimensional distribution information of the drill holes in the whole area is combined to construct a spatial interpolation structure. The point concentration is interpolated to generate a surface atlas, which is layered according to the depth dimension, forming a vertical pollution concentration distribution atlas. The atlas is output in 2D profile and 3D stereogram formats, which is convenient for subsequent superposition analysis with other environmental factors.

[0037] S502: Call the vertical pollution concentration distribution atlas, superimpose the groundwater level rising trend change map, locate the concentration change section in the atlas and extract the corresponding spatial position, judge whether the trend direction of the position is consistent with the pollution trend direction, and obtain the trend matching relationship determination result; After calling the vertical pollution concentration distribution map, it needs to be spatially superimposed with the groundwater level recovery trend change map. The groundwater level change map records the groundwater level elevation of each point at different monitoring periods, and calculates the trend change direction according to the time sequence. If the water level data of a certain area in four consecutive periods shows an upward trend, the trend direction is upward, marked as positive, otherwise, if it is downward, it is marked as negative. Match the depth and horizontal position of the high concentration section in the map with the same position in the groundwater level change map, extract the groundwater level trend identifier at that position, and judge the consistency of the direction with the pollution trend direction. The pollution trend direction is determined according to the change of concentration value in the depth direction. For example, if the concentration value of a certain concentration section increases layer by layer from top to bottom, it is determined to be a downward migration trend. If the groundwater level is rising and the pollution is migrating downward, the trend direction is consistent, otherwise it is inconsistent. Repeat the above steps to traverse all high concentration sections in the map to form a position-trend matching result set. In the result set, mark the spatial coordinates, pollution migration trend, groundwater trend, and trend consistency judgment result of each position. If the trend is consistent, mark it as "1", and if it is inconsistent, mark it as "0". Finally, generate the trend matching relationship judgment result.

[0038] S503: According to the trend matching relationship judgment result, set the corresponding repair period for the trend consistent area and the inconsistent area respectively, establish the repair frequency and time corresponding table according to the spatial segment number, arrange the intervention section operation according to the time sequence, and establish the dynamic development scheme of soil repair target; According to the trend matching relationship judgment result, set the intervention frequency for each identified point area. If the trend judgment result is consistent, that is, the trend mark is "1", then it is divided into a response synchronous area, and a higher frequency repair period is allocated to this area, for example, set to repair once every 15 days. If the trend is inconsistent, that is, the trend mark is "0", then it is divided into a response lag area, and a delayed intervention period is set, for example, extended to repair once every 30 days. The setting of frequency value refers to the pollution response period of similar sections in historical repair projects. If the average repair period of synchronous area in historical data is 14 to 17 days, the middle value 15 days can be taken as the setting value, and the lag area is doubled to 30 days. Then generate the repair frequency and repair period record corresponding to each segment number. Take the segment number as the primary key and the repair frequency (such as "15d") as the value item to establish the repair frequency and time corresponding table. The table is in structured data format, and the fields include segment number, trend judgment, repair period days, repair frequency level, etc. Finally, sort the repair arrangement of different segments according to the time axis, group the segments with the same frequency into the same execution batch, and mark the execution order in the table. Output the complete intervention section operation arrangement table to establish the dynamic development scheme of soil repair target. The scheme should cover all identified segments and associate their intervention time plan.

[0039] The above merely illustrates the specific embodiments of the present application, but the protection scope of the present application is not limited thereto, any person skilled in the art can easily think of the changes or replacements within the technical range disclosed by the present application, which should be covered in the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A method for dynamically establishing a soil restoration target based on groundwater source conservation, characterized by, The method comprises the following steps: S1: obtaining groundwater source conservation area repair starting point and water level data, extracting three-day water level rising trend and locating change inflection point, recording interval between repair time point, detecting evaporation and infiltration rate difference, obtaining repair intervention effective lag period; S2: calling the repair intervention effective lag period, selecting boundary expansion line, judging whether the period is shorter than the response time limit, if so, moving the boundary out, if not, freezing the boundary and marking the core area, obtaining the dynamic boundary map of the pollution control area; S3: calling the dynamic boundary map of the pollution control area, extracting water level gradient change and seepage direction consistency, screening water power continuous area and sorting intervention order, obtaining repair intervention priority sequence list; S4: calling the repair intervention priority sequence list, extracting capillary water layer upper and lower boundary difference and calculating fluctuation, detecting particle size and organic matter content, comparing fluctuation and lag water capacity limit value, obtaining intervention concentration control section marker map; S5: calling the intervention concentration control section marker map, extracting pollutant distribution map and marking high value section, superimposing groundwater level rising trend map, judging water content rate trend consistency, obtaining soil repair target dynamic development scheme.

2. The method according to claim 1, wherein, The repair intervention effective lag period includes water level change inflection point time sequence, lag interval length, rate difference amplitude, lag type classification, the dynamic boundary map of the pollution control area includes boundary adjustment state, core reserved area position, map revision range, response time limit comparison result, the repair intervention priority sequence list includes gradient change frequency, change rate sorting, seepage consistency ratio, water power continuous segment, intervention priority level, the intervention concentration control section marker map includes fluctuation amplitude index, particle size proportion estimate, organic matter content level, lag water capacity judgment result, limited section identification, the soil repair target dynamic development scheme includes concentration high value section distribution, water content rate trend direction, repair frequency setting, time cycle arrangement.

3. The method according to claim 1, wherein the method is characterized by, The response time limit refers to the shortest time threshold required for hydrological changes to respond to repair intervention measures after boundary expansion.

4. The method according to claim 1, wherein, The screening water power continuous area and sorting intervention order refers to identifying water power coherent segments and determining intervention priority according to water level gradient change frequency and seepage direction consistency.

5. The method according to claim 1, wherein, The specific steps of S1 are: S101: obtaining groundwater level monitoring data frame, extracting continuous three-day water level sequence, calculating daily increment, screening water level continuous rising time period, pairing end date with corresponding water level, obtaining continuous rising water level marker set; S102: calling the continuous rising water level marker set, calculating increment change rate, locating the time point corresponding to the maximum increment, judging the change of increment direction before and after, screening the time when the trend turns, obtaining water level change inflection point time node; S103: calling the water level change inflection point time node and repair starting point timestamp, calculating the time difference between them, eliminating inflection points earlier than the repair starting point, retaining time difference and aggregating, obtaining continuous rising water level marker set, and detecting evaporation rate and infiltration rate difference, generating repair intervention effective lag period.

6. The method according to claim 1, wherein, The specific steps of S2 are: S201: Call the repair intervention lag period, select the original extension line position of the pollution control area boundary, judge the difference between the lag period and the boundary response time limit, filter the boundary section with a shorter lag period than the response time limit, and obtain the boundary expansion adjustment section set; S202: According to the boundary expansion adjustment section set, perform the outward moving operation on the corresponding boundary, retain the boundary with a longer lag period than the response time limit and mark the core retention area, integrate the boundary changes, and generate the repair boundary adjustment layer; S203: Call the repair boundary adjustment layer, superimpose the original boundary map, identify the boundary change area, redraw the overall boundary line and mark the area attribute, extract the layer shape data, and generate the pollution control area dynamic boundary map.

7. The method according to claim 1, wherein the method is characterized by, The specific steps of S3 are: S301: Call the pollution control area dynamic boundary map, extract the water level data sequence between monitoring points, calculate the water level difference between adjacent monitoring points and count the positive and negative direction change frequency, and calculate the change rate of the absolute value of the continuous difference value, aggregate the frequency and rate data, and obtain the water level gradient change index set; S302: According to the water level gradient change index set, extract the seepage direction scalar value of adjacent monitoring segments, judge the direction consistency between adjacent segments, aggregate consistent segments and arrange them in order, extract the continuous segment set, and generate the water power continuity area list; S303: Call the water power continuity area list, extract the water level mutation difference in the area, mark the difference value range by level, sort the area number from high to low, establish the area and number mapping table and record the sorting position, and establish the repair intervention priority sequence list.

8. The method according to claim 1, wherein the method is characterized by, The specific steps of S4 are: S401: Call the marked segment in the repair intervention priority sequence list, extract the upper and lower boundary elevation values of the capillary water layer in the corresponding area, calculate the difference between the upper and lower boundaries, and calculate the difference between the maximum and minimum amplitude of the continuous difference sequence, and obtain the capillary water layer fluctuation amplitude value; S402: According to the capillary water layer fluctuation amplitude value, extract the soil particle size distribution and organic matter content data frame of the same segment, calculate the coarse grain component ratio and estimate the water storage capacity index combined with the organic matter content, compare the index with the water storage retention capacity threshold, and generate the restriction area identification result set; S403: Call the restriction area identification result set, locate the corresponding segment position, mark the area range in the pollution control area layer, and establish the mapping combined with the segment number and control section number, output the vector boundary and mark attribute field, and generate the intervention concentration control section mark graph.

9. The method according to claim 1, wherein the method is characterized by, The specific steps of S5 are: S501: Call the area identified by the intervention concentration control section mark graph, extract the pollution concentration sequence in the soil layer depth direction of the corresponding area, identify the concentration peak position and mark the space, construct the atlas structure combined with the profile depth data, and generate the vertical pollution concentration distribution atlas; S502: Call the vertical pollution concentration distribution atlas, superimpose the groundwater level recovery trend change graph, locate the concentration change section in the atlas and extract the corresponding spatial position, judge whether the trend direction of the water content rate at the position is consistent with the pollution trend direction, and obtain the trend matching relationship judgment result; S503: According to the trend matching relationship judgment result, the corresponding repair period is set for the trend consistent region and the inconsistent region respectively, the repair frequency and time corresponding table is established according to the space segment number, the intervention section operation arrangement is arranged in time sequence, and the soil repair target dynamic development scheme is established.

10. The method according to claim 8, wherein the method is characterized by, The water storage retention capacity threshold is used to determine the minimum water storage performance standard value of whether the soil layer has water retention capacity.

Images