A method and system for assessing the safety of groundwater source water supply

Abstract

The present invention provides a method and system for assessing the safety of water supply in a groundwater source, and relates to the technical field of water supply safety assessment. The method comprises: obtaining a permeability distribution map, dividing a grid and analyzing directional changes, constructing an edge set to extract a closed form to generate a partition layer, setting measuring points to collect water level data, evaluating the consistency of changes to screen paths, matching boundaries to determine crossing conditions, marking mutation positions to record interference paths, and establishing an identification list. In the present invention, by comparing the consistency of the water level direction per unit time with the response fluctuation, path segments with hydraulic coherence characteristics are identified, the ability to depict actual groundwater flow paths in a non-uniform recharge environment is improved, the position where the path direction mutation and the boundary coincidence are analyzed, the response identification of the structural interference area is strengthened, the key sections of water supply interruption or recharge anomaly are located, and the response capability of actual hydrological and geological changes and the spatial precision of the assessment and judgment are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of water supply safety assessment, and in particular to a method and system for assessing the safety of groundwater source water supply. Background Art

[0002] The technical field of water supply safety assessment involves quantitative or qualitative analysis and evaluation of factors such as water quality and quantity, water recharge, groundwater flow paths, pollutant migration, and the impact of human activities at riverside groundwater sources. The core of this technical field lies in the construction of a systematic assessment framework to scientifically analyze the hydrogeological structure of the water source, its hydrological recharge relationship, its hydraulic connection with surface water bodies, and its susceptibility to pollution risks, in order to determine its water supply stability and safety under different hydrological conditions. This overall technical field utilizes field surveys, hydrogeological monitoring, experimental simulations, and numerical simulations to comprehensively analyze the recharge sources, runoff processes, and response characteristics of groundwater resources to external disturbances, providing a scientific basis for the safe management of riverside water sources in urban water supply systems.

[0003] Among them, the groundwater source water supply safety assessment method refers to the safe water supply capacity of riverside water sources under different conditions. By establishing an analysis method based on the interaction between groundwater and surface water, combined with hydrogeological profiles, hydrological monitoring data and multi-time series groundwater level changes, the reliability and risk of its water supply are systematically evaluated. The subject of this patent covers the calculation of groundwater runoff direction and intensity, the deduction of pollutant migration path and speed, the boundary determination of water source recharge areas, and the dynamic response analysis of rainfall and surface water changes to groundwater systems. The assessment method is specifically completed by means of hydrological monitoring data processing, determination of borehole profile hydrological parameters, hydrodynamic gradient calculation, and groundwater-surface water exchange zone parameter measurement.

[0004] Existing techniques for identifying underground structures primarily rely on interpreting hydrogeological profiles and determining linear trends in geological maps. These techniques ignore localized continuous disturbances under atypical trends or small-scale fault conditions, resulting in reduced accuracy in identifying fault boundaries. Hydraulic path identification often relies on deducing regional water level contours, making it difficult to depict subtle changes in hydraulic path behavior under non-uniform recharge conditions. Recharge capacity assessment often uses a surface mean to represent recharge status, failing to distinguish between variations in recharge capacity under micro-topographic conditions and resulting in a poor match with actual recharge distribution. Existing methods often assess groundwater fluctuations separately from recharge capacity, failing to comprehensively identify sensitive areas with high fluctuation responses and insufficient recharge. In areas where high evaporation zones overlap with shallow water levels, water levels fluctuate frequently but lack sustained recharge capacity, making existing systems difficult to identify such high-risk areas. These shortcomings lead to the omission of sensitive locations in water supply security assessments and misjudgment of the spatial distribution of recharge capacity, compromising the scientific nature of regional water source allocation and scheduling strategies. Summary of the Invention

[0005] In order to solve the problems existing in the prior art, such as neglecting local continuous disturbances under atypical trends or small-scale fracture conditions, resulting in decreased accuracy in fracture boundary identification; in terms of hydraulic path identification, it mostly relies on regional water level contour deduction, which makes it difficult to depict the subtle changes in the hydraulic path under non-uniform recharge conditions; in terms of recharge capacity determination, the surface mean is often used to represent the recharge state, and there is a lack of distinction between the differences in recharge capacity performance under micro-topography undulations, resulting in a poor match with the actual recharge distribution; the existing methods mostly separate the groundwater fluctuation value and the recharge capacity for evaluation, and fail to comprehensively identify sensitive areas with high fluctuation response and insufficient recharge; under the condition of overlapping of strong evaporation zone and shallow buried water level area, although the water level fluctuates frequently, it does not have the ability to continuously recharge. The existing system has difficulty in identifying such high-risk areas, resulting in the omission of sensitive points in water supply safety assessment, misjudgment of the spatial distribution of recharge capacity, and the scientific nature of regional water source configuration and scheduling strategies. The present invention provides a method and system for assessing the water supply safety of groundwater sources. The technical solution is as follows:

[0006] In one aspect, a method for assessing the safety of groundwater source water supply is provided, the method comprising:

[0007] S1: Obtain the permeability distribution map of the area covered by the riverside water source, divide the image into equal-area grids, analyze the permeability changes in the vertical and horizontal directions of the units, construct an edge point set based on the positions of continuous jump points, extract the closed boundary shape, and generate a structural disturbance impact partition layer;

[0008] S2: calling the grid area in the structural disturbance impact partition layer, evaluating the consistency of the water level change direction and fluctuation time per unit time between the measuring points, screening the response path segments according to the degree of coherence, and obtaining the groundwater response path distribution set;

[0009] S3: Based on the groundwater response path distribution set, spatial matching is performed on the overlapping areas between the path segments and the boundary segments to determine whether the path line crosses the boundary or turns along the boundary, and the position of the path direction mutation position and the degree of overlap with the fracture boundary are marked to establish a path structure interference interactive identification list;

[0010] S4: According to the path structure interference interactive identification list, the precipitation data, evaporation data and elevation data corresponding to the unit are extracted, the precipitation and evaporation ratio is calculated, and then the elevation difference between adjacent cells is combined to form an undulation index value, which is divided into levels according to the recharge performance to generate a spatial recharge capacity level layer.

[0011] Optionally, the structural disturbance impact zoning layer includes a closed boundary line network, a disturbance boundary surface unit, and a continuous jump direction aggregation area; the groundwater response path distribution set includes a cross-section path line, a response direction marking point, and a water level continuity record index; the path structure interference interactive identification list includes a crossing path code, a junction change position identifier, and a recharge interruption area index; the spatial recharge capacity level layer includes a undulation correction index grid, a level zoning code map, and a recharge performance distribution layer.

[0012] Optionally, the steps for obtaining the structural disturbance impact partition layer are specifically as follows:

[0013] S101: Obtain a permeability distribution map of the area covered by the riverside water source, divide it into equal-area grids, treat the grid cells as independent spatial units, extract the initial soil permeability values ​​within the cell range, match the soil attribute point coordinates with the corresponding grid numbers, and generate a grid permeability distribution;

[0014] S102: Calling the permeability values ​​of adjacent grid cells in the grid permeability distribution, calculating the permeability change difference between the real-time cell and the adjacent cells in the horizontal and vertical directions, and performing a determination operation on each pair of direction groups to determine whether there is a direction jump, thereby generating a direction jump distribution point set;

[0015] S103: Based on the continuous spatial distribution of the jump points in the directional jump distribution point set, identify the jump point groups that meet the continuous quantity requirement, use the edge coordinates of the qualified point groups as edge candidate points, sort the edge candidate point coordinates by spatial proximity, extract the set of circumscribed boundary lines that constitute the closed area, and generate a structural disturbance impact partition layer.

[0016] Optionally, the step of obtaining the groundwater response path distribution set is specifically:

[0017] S201: calling the grid area in the structural disturbance impact partition layer, laying out a cross-section measuring point array in each group of grids according to a preset cross-section spacing, collecting continuous time water level change values ​​of the measuring points, organizing the water level time series data in the cross section according to the measuring point numbers, statistically analyzing the trend of the water level change direction and fluctuation amplitude between the measuring points in unit time, performing synchronization judgment on whether the direction and trend changes between the measuring point pairs in the cross section are consistent, and generating a consistency assessment result of the measuring point water level change;

[0018] S202: Based on the consistency evaluation results of the water level changes at the measuring points, the measuring point pairs that meet the consistency requirements are screened in the cross section, connected paths are constructed using the connecting lines of the measuring point pairs as line segments, connected distance correction values ​​are calculated, and the spatial direction distributions of the multiple connected paths are integrated to generate a groundwater response path distribution set.

[0019] Optionally, the formula for calculating the connectivity distance correction value is:

[0020] ;

[0021] Among them, L ij Represents the connected distance correction value between measuring point i and measuring point j, x i 、x j Represents the horizontal coordinates of measuring point i and measuring point j, y i 、y j Represent the vertical coordinates of measuring point i and measuring point j respectively, n represents the total number of observation times, It represents the water level difference between measuring point i and measuring point j at the kth moment, It represents the average value of the water level average value of measuring point i and measuring point j during the entire period. They represent the average water level of measuring points i and j during the entire period, and hc represents the centralized reference value of the water level of the cross-section measuring point during the entire period.

[0022] Optionally, the step of obtaining the path structure interference interaction identification list is specifically as follows:

[0023] S301: Calling the spatial coordinates of the path segments in the groundwater response path distribution set, and performing spatial matching with the coordinate set of the boundary segments in the structural disturbance impact partition layer, determining whether the path segments have intersections or overlaps with the boundary segments in the coordinate range, if there is an intersection, marking it as crossing, if there is an overlap, marking it as turning along the boundary, and generating a path line and boundary matching type dataset;

[0024] S302: Extracting the path direction change angle based on the path segments marked as crossing or turning in the path line and boundary matching type dataset, calculating the direction mutation value based on the angle difference between the previous and next segments of the path, and using the boundary segment coordinates to perform spatial overlap judgment to generate a path boundary mutation overlap annotation result;

[0025] S303: Calling the mutation path segment number in the path boundary mutation coincidence degree annotation result, marking the path segment with disturbance signs, summarizing the abnormal marked areas at the partition boundary according to the path number, extracting the corresponding path segment and its location as the identification object, calculating the path segment disturbance identification value, and generating a path structure interference interactive identification list;

[0026] The formula for calculating the path segment disturbance identification value is:

[0027] ;

[0028] Among them, D a is the path segment disturbance identification value, N is the total number of path segments, S ao is the matching score between the oth path segment and the ath path segment, Mo is the average value of path segment o, σ o is the standard deviation of path segment o, α is the adjustment coefficient, W o is the path segment weight coefficient.

[0029] Optionally, the steps for obtaining the space replenishment capability level layer are specifically as follows:

[0030] S401: Based on the unit numbers listed in the path structure interference interactive identification list, precipitation raster data, evaporation raster data, and elevation raster data of the corresponding unit area are extracted, the spatial resolution and projection format of the three types of data are unified, and the raster data are cropped according to the unit boundary range to obtain precipitation-evaporation ratio distribution data;

[0031] S402: Based on the drop-evaporation ratio distribution data, call the elevation value of the corresponding pixel position in the clipped elevation raster data, compare the elevation difference between the real-time pixel and the neighboring pixels, calculate the average neighborhood elevation difference, and multiply the average elevation difference by the drop-evaporation ratio corresponding to the pixel to obtain the relief index distribution value;

[0032] S403: Call the pixel data of the relief index distribution value, divide it into levels according to the numerical range, use the set relief index level threshold range as the classification basis, assign differentiated level values ​​respectively, and generate a level grid layer according to the spatial position of the division result to obtain the spatial replenishment capacity level layer.

[0033] Optionally, the method further comprises step S5:

[0034] S5: calling the spatial recharge capacity level layer and the marked areas in the path structure interference interactive identification list, superimposing the annual groundwater level fluctuation data in the area, extracting regional units with structural interference and recharge instability, analyzing the response amplitude and frequency characteristics, and obtaining a water supply security sensitive distribution area layer;

[0035] The water supply safety sensitive distribution area layer includes response abnormality unit mark, water level fluctuation frequency index, and structural recharge overlapping area code.

[0036] Optionally, the steps for obtaining the water supply safety sensitive distribution area layer are specifically as follows:

[0037] S501: The spatial recharge capacity level layer is called to identify the marked area in the path structure interference interactive identification list. For the recharge units in the marked area, the annual groundwater level monitoring data of the corresponding location is collected. The groundwater level values ​​at the same location are arranged in a monthly sequence. The water level fluctuation amplitude of the recharge unit in the year is calculated, and the recharge units that do not fluctuate are screened out to obtain a groundwater level fluctuation amplitude dataset.

[0038] S502: Based on the groundwater level fluctuation amplitude dataset and in combination with the regional interference level zoning results in the path structure interference interactive identification list, the fluctuation amplitude value change frequency of the recharge units within the interference level area is counted, and the recharge units whose change frequency exceeds the regional annual change benchmark are screened to obtain the interference area unit sequence;

[0039] S503: According to the interference area unit sequence, the spatial position index value of the corresponding unit is extracted, and combined with the grade division of the corresponding position in the spatial replenishment capacity grade layer, the area whose distribution ratio exceeds the internal benchmark level of the grade area to which it belongs is screened through index matching operation to generate a water supply safety sensitive distribution area layer.

[0040] On the other hand, a groundwater source water supply safety assessment system is provided for executing the above-mentioned groundwater source water supply safety assessment method, the system comprising:

[0041] The structural disturbance identification module obtains the permeability distribution map of the area covered by the riverside water source. After dividing the image into equal-area grids, the permeability values ​​within the grids are called in the vertical and horizontal directions. After forming unit pairs, the permeability differences between adjacent grids in the same direction are compared. The closed boundaries are connected in sequence and the spatial contours are extracted to generate the structural disturbance impact partition layer.

[0042] A response path extraction module sets a cross-sectional measurement point array based on the grid area within the boundary covered by the structural disturbance impact partition layer, collects groundwater level time series data at the corresponding locations, and selects line segments with coherent change trends to obtain a groundwater response path distribution set;

[0043] The interference mechanism identification module performs spatial overlap judgment with the boundary line segment of the disturbance limit distribution morphological value based on the coordinates of the corresponding path segment in the groundwater response path distribution set, marks whether the path turns or changes its extension direction at the intersection, and obtains a path structure interference interactive identification list;

[0044] The recharge capacity classification module calls the path structure interference interactive identification list, extracts the precipitation data, evaporation data and elevation data of the unit, divides the recharge level according to the combination level of the calculation results, and generates a spatial recharge capacity level layer;

[0045] The water supply sensitive area identification module, based on the spatial recharge capacity level layer, superimposes the regional annual groundwater level change data, calculates the joint value of the water level response amplitude and fluctuation frequency corresponding to the interference unit in the recharge level area, and screens out the units whose joint value exceeds the distribution mean of the level area to obtain the water supply safety sensitive distribution area layer.

[0046] The beneficial effects brought about by the technical solution provided by the embodiment of the present invention include at least:

[0047] In an embodiment of the present invention, by determining the jump characteristics of permeability change direction and extracting closed boundary morphology based on the spatial clustering of edge points, the structural identification stage no longer relies on geological boundary indices or empirical judgment, thereby enhancing the resolution capability of spatially continuous disturbance identification. In path identification, by comparing the consistency of water level direction per unit time with response fluctuations, path segments with hydraulic coherence are identified, improving the ability to characterize actual groundwater flow paths in non-uniform recharge environments. During the path and boundary interaction stage, spatial superposition is used to analyze the locations where sudden changes in path direction coincide with boundary overlaps, enhancing the response identification of areas affected by structural disturbances and effectively locating key sections of water supply interruptions or recharge anomalies. In recharge capacity identification, fluctuation indicators are used to modify the precipitation-evaporation relationship, linking local terrain characteristics with recharge performance and distinguishing the differences in the ability of different geomorphic units to compensate for water supply. By superimposing annual groundwater level fluctuations, a composite judgment criterion between recharge performance and fluctuation amplitude is established, achieving high recognition accuracy for water supply-sensitive areas. The processing logic emphasizes the four major elements of structure, hydraulics, topography and dynamic response, and sequentially constructs the basis for assessment, integrating the identification of underground structures, identification of hydraulic paths, quantification of recharge stability and spatial expression of risks, thus overall improving the assessment system's ability to respond to actual hydrological and geological changes and the spatial precision of assessment judgments. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 It is a schematic diagram of the workflow of the present invention;

[0049] Figure 2 It is a system flow chart of the present invention. DETAILED DESCRIPTION

[0050] The technical solution of the present invention is described below in conjunction with the accompanying drawings.

[0051] In the embodiments of the present invention, words such as "exemplarily" and "for example" are used to indicate examples, illustrations, or explanations. Any embodiment or design described as an "exemplary" in the present invention should not be interpreted as being preferred or advantageous over other embodiments or designs. Rather, the use of the word "exemplary" is intended to present concepts in a concrete manner. Furthermore, in the embodiments of the present invention, "and / or" can mean both or either of the two.

[0052] In order to make the technical problems, technical solutions and advantages to be solved by the present invention clearer, a detailed description will be given below with reference to the accompanying drawings and specific embodiments.

[0053] See also Figure 1 The embodiment of the present invention provides a method for assessing the safety of groundwater source water supply. The processing flow of the method may include the following steps:

[0054] S1: Obtain the permeability distribution map of the area covered by the riverside water source, divide the image into equal-area grids, analyze the vertical and horizontal permeability changes of each unit, construct direction pairs with adjacent units, and then determine the direction jump characteristics. According to the positions of continuous jump points, construct edge point sets and extract closed boundary shapes to generate a structural disturbance impact partition layer;

[0055] S2: Call the grid area in the structural disturbance impact partition layer, set the cross-section measurement point array and collect water level time series data, evaluate the consistency of the water level change direction and fluctuation time per unit time between the measurement points, filter the response path segments according to the degree of coherence, and obtain the groundwater response path distribution set;

[0056] S3: Based on the distribution set of groundwater response paths, spatially match the overlapping areas between path segments and boundary segments to determine whether the path line crosses the boundary or bends along the boundary. Annotate the location of the path's sudden change in direction and the degree of overlap with the fault boundary. Record the path segments with direction changes as potential interference paths. Combined with the continuity change trend of the path in the intersection area, mark the recharge anomaly area and establish a path structure interference interactive identification list.

[0057] S4: Based on the interactive identification list of path structure interference, the corresponding precipitation data, evaporation data, and elevation data of the unit are extracted. The precipitation-evaporation ratio is calculated and then combined with the height difference of adjacent cells to form an ups and downs index value. The recharge performance is divided into levels to generate a spatial recharge capacity level layer.

[0058] S5: Call the spatial recharge capacity level layer and the marked areas in the path structure interference interactive identification list, superimpose the annual fluctuation data of the groundwater level in the area, extract the regional units with structural interference and unstable recharge, analyze the response amplitude and frequency characteristics, and obtain the water supply security sensitive distribution area layer.

[0059] Among them, the structural disturbance impact zoning layer includes closed boundary line networks, disturbance boundary surface units, and continuous jump direction aggregation areas; the groundwater response path distribution set includes cross-section path lines, response direction marking points, and water level continuity record indexes; the path structure interference interactive identification list includes crossing path codes, intersection change location markers, and recharge interruption area indexes; the spatial recharge capacity level layer includes undulation correction index grids, level zoning coding maps, and recharge performance distribution layers; the water supply safety sensitive distribution area layer includes response anomaly unit markers, water level fluctuation frequency indexes, and structural recharge overlap area codes.

[0060] The specific steps for obtaining the structural disturbance impact partition layer are as follows:

[0061] S101: Obtain a permeability distribution map of the area covered by the riverside water source, divide it into equal-area grids, treat the grid cells as independent spatial units, extract the initial soil permeability values ​​within the cell range, match the soil attribute point coordinates with the corresponding grid numbers, and generate a grid permeability distribution;

[0062] It is necessary to use geographic information software to load the remote sensing image or topographic map of the riverside water source area, determine the boundary range of the area by manual outlining or automatic identification, and use a fixed scale to divide the equal area grid within the boundary range. Set each grid as a 100m x 100m square unit, and the entire area will be automatically divided into several numbered units. Each unit has a unique number and spatial coordinate range. Import the soil permeability survey data points. Each point contains the geographic coordinates and the original permeability value. Use spatial matching to match the points with the grid. If a soil point is located in the unit numbered A25 If a cell contains multiple points, its value is assigned to cell A25. If a cell contains multiple points, it is necessary to extract the original permeability value of the point and perform statistical integration on it. The average method or median method can be used to obtain the representative value of the cell. Suppose a cell contains 3 points, located in the upper left corner, middle and lower right corner, with permeabilities of 2.3, 2.5 and 2.7 respectively. After integration, the representative value of the cell can be set as 2.5. At the same time, excluding cells with no point data, its value can be estimated by spatial interpolation, such as referring to the data of the surrounding adjacent cells. After completion, the cell number and its permeability are integrated into a distribution map to generate the grid permeability distribution.

[0063] S102: Calling the permeability values ​​of adjacent grid cells in the grid permeability distribution, calculating the permeability change difference between the real-time cell and the adjacent cells in the horizontal and vertical directions, and performing a determination operation on each pair of direction groups to determine whether there is a direction jump, thereby generating a direction jump distribution point set;

[0064] The data of adjacent cells of each grid cell are extracted and compared in the horizontal and vertical directions respectively. The permeability values ​​of the cells on the left and right, and above and below are set to be compared. The differences in the two directions are calculated and recorded separately. If a cell shows a large difference in permeability between its right neighbor and it in the horizontal comparison, and a small change in the vertical comparison between its lower neighbor, the horizontal difference value is marked as a key jump data. The jump judgment criteria are further set to determine whether the difference meets the conditions for determining a jump. If it does, the current cell is recorded as a jump point, its spatial location is recorded, and added to the jump point dataset. This process needs to be executed for all grid cells in the entire area. If an area contains 400 grid cells, all grids will be traversed, and the change relationship in the four directions will be analyzed. The permeability difference of each direction pair is identified and classified to obtain the jump point layer data. Each point in this layer represents a jump location, including its grid number and directional difference characteristics. It is used for subsequent jump point group identification and structural perturbation analysis to generate a directional jump distribution point set.

[0065] S103: Based on the continuous spatial distribution of the jump points in the directional jump distribution point set, identify jump point groups that meet the continuous quantity requirement, use the edge coordinates of the qualified point groups as edge candidate points, sort the edge candidate point coordinates by spatial proximity, extract the set of circumscribed boundary lines that constitute the closed area, and generate a structural disturbance impact partition layer;

[0066] The jump point layer is spatially analyzed. Based on the distribution density and continuity of the jump points in space, jump point groups that meet the set conditions are extracted. It is set that a continuous jump point group must contain at least 5 adjacent points. The spatial distance between the jump point and the surrounding points will be checked in turn to see if it meets the proximity condition. If so, they are classified into the same jump point group. After the point group identification is completed, the edge points of each group of point groups will be extracted. The boundary sorting algorithm is used to connect the edge points in the spatial direction to form a closed boundary path. This path represents the maximum circumscribed range of the jump point distribution. The boundary line set of the area enclosed by the closed path is further drawn. If a jump point group is distributed in the southwest area of ​​the river bank, its outer edge points are arranged in the order of southeast-northeast-northwest-southwest, and their connecting lines will form a convex quadrilateral or polygon. The boundary line set is exported as a spatial layer for structural disturbance impact zoning and superimposed on the existing water source map layer to form a composite layer, providing input data for subsequent hydrological simulation and risk analysis, and generating a structural disturbance impact zoning layer.

[0067] The specific steps for obtaining the groundwater response path distribution set are:

[0068] S201: Calling the grid area in the structural disturbance impact partition layer, laying out a cross-section measuring point array in each group of grids according to the preset cross-section spacing, collecting the continuous time water level change values ​​of the measuring points, organizing the water level time series data in the cross section according to the measuring point number, and calculating the trend of the water level change direction and fluctuation amplitude between the measuring points in unit time. Whether the direction and trend changes between the measuring points in the cross section are consistent is synchronously judged, and a consistency assessment result of the measuring point water level change is generated;

[0069] The cross-section spacing is usually set according to the on-site terrain characteristics, hydrogeological conditions and survey accuracy requirements. If a cross-section is laid out every 200 meters, the grid line crossing the disturbance partition will be automatically selected as the cross-section within the spacing interval, and several measuring points will be evenly laid out on each cross-section. The measuring point spacing can be set to 20 meters to ensure that the full width of the cross-section is covered. The name of each measuring point needs to be coded to include the cross-section number and point sequence. The second measuring point on the No. 3 section is named D3P2. Each measuring point is equipped with a groundwater level observation device or a remote sensing node to collect water level change data regularly or continuously. The data is recorded in the form of a time series list. The data of different measuring points are organized in a unified format according to the time dimension. For points such as 3P3, a continuous list of time-water level pair data is generated. The water level change trend between each measuring point is counted according to the unit time step. The difference change between adjacent moments can be used to represent the rising or falling direction. If both measuring points show an upward trend within a 5-minute time period, the recorded direction is consistent. At the same time, the change amplitude between each pair of measuring points is counted to see if it is within the same level range. The fluctuation amplitude level boundary value (such as 2cm, 5cm, 5~10cm, etc.) is set to determine whether it belongs to the same level, and the water level fluctuation synchronization between the measuring point pairs is judged based on this. After the measuring point pair is judged, each pair of measuring points is marked to see if it is consistent, and an attribute table is formed, listing the measuring point pair number, direction consistency mark, and fluctuation amplitude level consistency mark. This is used as the basic data for subsequent path construction to generate the consistency assessment result of the measuring point water level change.

[0070] S202: Based on the consistency evaluation results of the water level changes at the measuring points, select the measuring point pairs that meet the consistency requirements in the cross section, construct connected paths using the connecting lines of the measuring point pairs as line segments, calculate the connected distance correction value, integrate the spatial trend distribution of multiple connected paths, and generate a groundwater response path distribution set;

[0071] The formula for calculating the connectivity distance correction value is:

[0072] ;

[0073] Among them, L ij Represents the connected distance correction value between measuring point i and measuring point j, x i 、x jRepresents the horizontal coordinates of measuring point i and measuring point j, y i 、y j Represent the vertical coordinates of measuring point i and measuring point j respectively, n represents the total number of observation times, It represents the water level difference between measuring point i and measuring point j at the kth moment, It represents the average value of the water level average value of measuring point i and measuring point j during the entire period. They represent the average water level of measuring points i and j during the entire period, and hc represents the centralized reference value of the water level of the cross-section measuring point during the entire period.

[0074] Parameter meaning and formula calculation derivation process:

[0075] Calculation of measuring point coordinate difference:

[0076] ;

[0077] This item calculates the straight-line distance between measuring points i and j in two-dimensional space;

[0078] x i 、y i is the coordinate value of measuring point i, x j 、y j is the coordinate value of measuring point j;

[0079] The coordinates of measuring point i are (10, 5), and the coordinates of measuring point j are (15, 8), then:

[0080] ;

[0081] The straight-line distance between measuring points i and j on the plane is approximately 5.83 units;

[0082] Mean absolute difference of water level differences:

[0083] ;

[0084] This item calculates the average absolute value of the water level difference between measuring point i and measuring point j at each observation time k, reflecting the consistency of water level changes;

[0085] n represents the number of observation points;

[0086] is the water level difference between measuring point i and measuring point j at the kth moment,

[0087] Set n = 5, and the water level difference at each moment is [0.5, 0.7, 0.3, 0.4, 0.6], then:

[0088] Mean absolute difference:

[0089] ;

[0090] This means that the average water level difference between measuring points i and j during the observation period is 0.5 units;

[0091] The difference between the water level difference and the reference value:

[0092] ;

[0093] This item calculates the absolute difference between the mean of the water levels at measuring points i and j and the global reference water level hc;

[0094] is the average water level of measuring point i and measuring point j at all times;

[0095] hc is the centralized reference value of the water level at all measuring points in the cross section during the entire period, such as the median or overall mean;

[0096] The average water levels of measuring points i and j are set to and ,but:

[0097] ;

[0098] Set the reference water level hc = 3.0, then:

[0099] ;

[0100] This shows that the difference between the average of the water levels at measuring points i and j and the global reference water level is 0.35 units;

[0101] Substitute into the formula for calculation:

[0102] ;

[0103] The results show that the correction value of the connectivity path between measuring points i and j is 6.68 units. This value not only considers the linear distance in space, but also combines the consistency of water level changes and the difference from the global water level to correct the connectivity path.

[0104] The specific steps for obtaining the path structure interference interaction identification list are as follows:

[0105] S301: Call the spatial coordinates of the path segments in the groundwater response path distribution set and perform spatial matching with the coordinate set of the boundary segments in the structural disturbance impact partition layer. Determine whether the path segments have intersections or overlaps with the boundary segments in the coordinate range. If there is an intersection, mark it as crossing; if there is an overlap, mark it as turning along the boundary. Generate a path and boundary matching type dataset.

[0106] The path segments are spatially overlapped with the boundary segments in the structural disturbance impact partition layer in turn. The specific operation process is as follows: read the starting and ending coordinates of the path segments, build a continuous segment chain, and number and index the paths one by one, extract the boundary segment set from the disturbance partition boundary layer, establish a spatial index for each boundary line, and use the fast bounding box method of the coordinate range to determine whether there is a potential contact probability between the path and the boundary. After screening out the contact candidate pairs, perform a more accurate segment intersection test to determine whether there is an intersection point. If the test result is that the path segment intersects the boundary segment and the intersection point falls within the segment range, the path is marked as a crossing match. If the test result is that the path segment and the boundary segment partially or completely overlap, especially at multiple points If continuous matching in precise coordinates exceeds the set length threshold, it is marked as a boundary-bending match. A path is numbered P21, with endpoint coordinates from (135.221, 36.447) to (135.245, 36.472). If it intersects with the boundary segment numbered B12, the intersection information is recorded in the P21 attribute and assigned a crossing label. At the same time, the numbers of all paths with intersections or overlapping segments are recorded in a path boundary matching type dataset. This dataset includes key fields such as the path number, matching type (crossing or turning), matching boundary number, and intersection coordinates or the start and end positions of the overlapping segment. This serves as the basis for subsequent trend analysis and interference identification, generating a path line and boundary matching type dataset.

[0107] S302: Extract the path direction change angle based on the path segments marked as crossing or turning in the path line and boundary matching type dataset, calculate the direction mutation value based on the angle difference between the previous and next path segments, and use the boundary segment coordinates to determine the spatial overlap, thereby generating the path boundary mutation overlap annotation result;

[0108] Filter the path segments marked as crossing or turning, extract the direction angles of the path segments before and after the crossing or turning points, and calculate the change in the direction angle to determine the degree of sudden change in the direction. The operation steps include: reading the coordinates of the connecting line segments before and after the boundary intersection point, obtaining their direction vectors respectively, and calculating the sudden change angle of the path direction based on the angle formed by the vectors. If the angle change exceeds the preset threshold (such as 30 degrees), it is regarded as a sudden change point, and its sudden change value is recorded. At the same time, the boundary segment coordinates are read again at this position to determine whether there is spatial overlap between the path and the boundary, and the overlapping length accounts for the boundary length or the path length. The ratio of path lengths can be used as an overlap index. If the ratio exceeds 50%, it is marked as significant overlap. The mutation angle and overlap index are integrated into the annotation results. Each record includes the path number, mutation angle, whether it is a mutation flag, overlap percentage, mutation point coordinates, and the corresponding boundary segment number. If the direction angle difference between the front and rear segments of path P31 is 48 degrees, and the mutation segment has a 70% overlap ratio with boundary B22, then the record will be assigned a dual label of mutation and significant overlap, providing a direct identification basis for the next step of identifying the disturbed path and generating the path boundary mutation overlap annotation result.

[0109] S303: Call the mutation path segment number in the path boundary mutation coincidence annotation result, mark the path segment with disturbance signs, summarize the abnormal annotated areas at the partition boundary according to the path number, extract the corresponding path segment and its location as the identification object, calculate the path segment disturbance identification value, and generate the path structure interference interactive identification list;

[0110] The formula for calculating the path segment disturbance identification value is:

[0111] ;

[0112] Among them, D a is the path segment disturbance identification value, N is the total number of path segments, S ao is the matching score between the oth path segment and the ath path segment, M o is the average value of path segment o, σ o is the standard deviation of path segment o, α is the adjustment coefficient, W o is the path segment weight coefficient.

[0113] Parameter meaning and formula calculation derivation process:

[0114] Path segment matching score (S ao ):

[0115] Calculate the matching score between path segment a and path segment o using a path matching algorithm (such as the dynamic time warping algorithm), and set the matching score between path segment a and path segment o to 0.85;

[0116] The average value of the path segment (Mo ):

[0117] Calculate the average value of path segment o, which is the arithmetic mean of all eigenvalues ​​of the path segment. Set the eigenvalue of path segment o to [0.8, 0.9, 0.85], then its average value M o for:

[0118] ;

[0119] Path segment standard deviation (σ o ):

[0120] Calculate the standard deviation of the path segment o, which represents the discrete degree of the path segment eigenvalue. Set the eigenvalue of the path segment o to [0.8, 0.9, 0.85], then its standard deviation σ o for:

[0121] ;

[0122] Adjustment coefficient (α):

[0123] The adjustment coefficient is used to control the influence of the normalized deviation. It is set according to actual needs and α=2.

[0124] Path segment weight coefficient (W o ):

[0125] Set the weight coefficient according to the importance or confidence of the path segment, and set the weight coefficient W of the path segment o o = 1.5;

[0126] Substituting the above parameters into the formula:

[0127] ;

[0128] For path segments a and o, substitute the known values:

[0129] ;

[0130] Due to S ao = M o , so the disturbance identification value D a = 0, indicating that there is no significant disturbance between path segments a and o;

[0131] This result shows that the matching degree between path segment a and path segment o is high, and the discreteness of their eigenvalues ​​is low, so there is no need to specially mark path segments with signs of disturbance in the path structure.

[0132] The steps to obtain the space supply capacity level layer are as follows:

[0133] S401: Based on the cell numbers listed in the path structure interference interactive identification list, precipitation raster data, evaporation raster data, and elevation raster data of the corresponding cell area are extracted. After unifying the spatial resolution and projection format of the three types of data, the raster data are cropped according to the cell boundary range to obtain the precipitation-evaporation ratio distribution data;

[0134] The spatial unit range corresponding to each number is retrieved in turn, and its boundary outline is located through the spatial index. The precipitation raster data, evaporation raster data and digital elevation model data within the coverage area of ​​the area are extracted from the remote sensing data platform or the preset raster database. The formats are converted with a unified spatial reference (such as WGS84) and the same resolution (such as 30 meters) to ensure that the projection format and pixel size of the three data sources are consistent. After completion, the three types of data are subjected to raster clipping operations according to the unit boundaries, that is, only the pixel area within the identified unit range is retained. A unit number U107 is set to be located in the southwest of the layer, and the corresponding spatial range is 13 5.215 longitude, 36.455 latitude. The clipping operation extracts the precipitation, evaporation, and elevation data pixels within this range to form three sets of unified raster data. The clipped precipitation data and evaporation data are compared pixel by pixel and a ratio operation is performed to form a precipitation-evaporation ratio distribution layer. The ratio represents the water surplus or deficit at each pixel. If the precipitation of a pixel is 85 mm and the evaporation is 65 mm, then its precipitation-evaporation ratio is 1.31. If the evaporation of another pixel is higher than the precipitation, the ratio is less than 1. The basic data for the dynamic evaluation of water in each unit is completed through the ratio calculation, which provides input for the subsequent recharge capacity assessment and obtains the precipitation-evaporation ratio distribution data.

[0135] S402: Based on the drop-evaporation ratio distribution data, the elevation value of the corresponding pixel position in the clipped elevation raster data is called, the elevation difference between the real-time pixel and the neighboring pixels is compared, the average elevation difference of the neighborhood is calculated, and the average elevation difference is multiplied by the drop-evaporation ratio corresponding to the pixel to obtain the relief index distribution value;

[0136] Read the ratio result of each pixel, and perform neighborhood analysis on the position of the pixel in the elevation raster data, extract the elevation values ​​of the adjacent pixels within its eight neighborhoods, calculate the height difference between the current pixel and its neighbors, average the height differences, and obtain the average height difference of the neighborhood of the pixel. The height difference can reflect the degree of terrain undulation. Set the elevation of a pixel to 310 meters, and the elevations of its surrounding pixels are 308, 312, 315, 309, 305, 307, 311, and 313 respectively. The average height difference is the average value of the difference set. After obtaining the result , the pixel's drop-evaporation ratio is multiplied by its average height difference, which is called the relief index distribution value. The comprehensive factors of terrain relief and hydrological surplus and deficit capacity are integrated, and this calculation process is performed on the pixels in sequence to form a complete relief index layer. Set in a certain unit, if the drop-evaporation ratio of a certain pixel is 1.25 and the average height difference of the neighborhood is 6.2 meters, then the relief index value of the pixel is 7.75. After the pixel is executed, the index layer can be used to reveal the spatial recharge capacity differentiation characteristics of the groundwater response path of the interaction between terrain and hydrology, and obtain the relief index distribution value.

[0137] S403: Retrieving pixel data of relief index distribution values, classifying them by numerical range, using the set relief index level threshold range as the classification basis, assigning differentiated level values, and generating a level grid layer based on the division results according to spatial position to obtain a spatial recharge capacity level layer;

[0138] The pixel relief index value is used as the classification object and is graded according to the set grade threshold interval. The grade classification rules can be preset into several grade segments according to the landform type and regional characteristics. The relief index value is set to be divided into: low grade (05), medium-low grade (10), medium grade (15), medium-high grade (20), and high grade (greater than 20). All pixels are assigned corresponding grade identification values ​​according to the interval in which their values ​​fall, such as grade 1 to grade 5. After completing the grade assignment of the pixels, a new grade raster layer is generated according to their spatial position, that is, each pixel corresponds to a grade identification in the layer, forming a spatial recharge capacity grade layer. The layer can be directly loaded and browsed in the GIS platform, and can be superimposed on the path disturbance layer or the structural partition layer to form a linkage analysis result. The layer is set to show that most path interaction sections in a certain disturbance partition are located in the high-grade recharge interval, which means that the area has good water recharge capacity. The grade layer is output for further analysis or visualization expression to obtain the spatial recharge capacity grade layer.

[0139] The specific steps for obtaining the water supply security sensitive distribution area layer are as follows:

[0140] S501: The spatial recharge capacity level layer is used to identify the marked areas in the interactive identification list of path structure interference. For the recharge units in the marked areas, the annual groundwater level monitoring data of the corresponding locations are collected. After the groundwater level values ​​at the same location are arranged in a monthly sequence, the annual water level fluctuation amplitude of the recharge units is calculated. The recharge units that do not fluctuate are screened out to obtain the groundwater level fluctuation amplitude dataset.

[0141] After calling the spatial recharge capacity level layer and the annotated area included in the path structure interference interactive identification list, the two layers are superimposed, and the recharge level unit in the annotated area is extracted through spatial intersection. The unique code and position coordinates of the unit in the spatial layer are identified, and the annual groundwater level monitoring data of the unit are read one by one. The data source can be an automatic water level recorder or a monitoring record table to ensure that the monthly data is complete and time-series continuous. The monthly water level values ​​of each unit are sorted to construct an annual time series table. The monthly water level values ​​of the unit numbered U308 are 12.3, 12.1, 11.9, 11.7, 11.8, 12.0, 12.2, 12.5, 12.6, 12.4, 12.2, and 12.1 meters. The difference between the maximum and minimum values ​​is automatically identified as the annual amplitude. In this case, the amplitude is 0.9 meters. The same operation is performed on all recharge units. The water level amplitude is extracted and recorded in the table. The results are screened and the units with no change in water level values ​​throughout the year or all monthly values ​​are exactly the same are eliminated. Because they do not have fluctuation characteristics, they are not included in the subsequent analysis. Each record contains the unit number, the maximum water level value throughout the year, the minimum water level value, the amplitude, and the spatial location, which lays the foundation for further identification of the intensity of interference impact and obtains the groundwater level fluctuation amplitude dataset.

[0142] S502: Based on the groundwater level fluctuation amplitude dataset and the regional interference level zoning results in the path structure interference interactive identification list, the fluctuation amplitude value change frequency of the recharge units within the interference level area is counted, and the recharge units whose change frequency exceeds the regional annual change benchmark are selected to obtain the interference area unit sequence;

[0143] The data were matched and analyzed with the regional interference levels in the path structure interference interactive identification list to identify the interference level area of ​​each recharge unit. The interference level was divided into low, medium and high level areas. The water level fluctuation amplitude of the recharge units in each level area was counted, and the frequency distribution of each fluctuation amplitude value in different intervals was recorded. For example, the fluctuation amplitude occurred 8 times within the 0.2 meter range, 12 times within the 0.4 meter range, and 5 times above the 0.6 meter range. The annual fluctuation frequency distribution of the recharge units in the area was used as a benchmark, and the annual change benchmark frequency of the area was set. For example, the frequency of occurrence was set to be greater than a certain set threshold (such as 1.5 times the average number of occurrences in the area). The fluctuation frequency of all recharge units was compared with the benchmark, and the recharge units with a change frequency exceeding the benchmark were screened out. The units were considered to have active fluctuations and have the characteristics of being enhanced by interference, including unit number, annual fluctuation amplitude, interference level area, frequency excess multiple, spatial location, etc., which provided a fluctuation basis for identifying water supply sensitive areas and formed a sequence of interference area units.

[0144] S503: Extract the spatial position index value of the corresponding unit based on the interference area unit sequence. Combined with the grade division of the corresponding position in the spatial recharge capacity grade layer, perform an index matching operation to screen areas whose distribution ratio exceeds the internal benchmark level of the grade area to which they belong, and generate a water supply security sensitive distribution area layer.

[0145] The corresponding level categories of the units in the spatial recharge capacity level layer are extracted, and their spatial position index values ​​are called to form a "unit number-recharge level" matching table. The spatial matching algorithm is used to count the number of interference units contained in each recharge level category and their proportion to the total number of units in the corresponding level. It is assumed that there are 50 units in the recharge level 4 area, of which 18 units appear in the interference unit sequence, and the interference distribution ratio is 36%. A benchmark level is set within the level area. For example, the benchmark is 30%. When the actual ratio exceeds this benchmark level, it is considered that the recharge level area has water supply security risk sensitivity. The area is extracted and marked as a sensitive area. The sensitive area layer is reconstructed using the spatial index value. Each identified block in this layer contains information such as its level category, interference unit ratio, excess multiple, and spatial boundary. The area numbered R42 is set as a level 4 recharge area with an interference ratio of 38%, which exceeds the benchmark by 30%, and is therefore marked as a sensitive area. This layer will serve as an important spatial reference for groundwater resource safety management and generate a water supply security sensitive distribution area layer.

[0146] See also Figure 2 , a groundwater source water supply safety assessment system, the system includes:

[0147] The structural disturbance identification module obtains the permeability distribution map of the area covered by the riverside water source. After dividing the image into equal-area grids, the permeability values ​​within the grids are called in the vertical and horizontal directions. After forming unit pairs, the permeability differences between adjacent grids in the same direction are compared. The closed boundaries are connected in sequence and the spatial contours are extracted to generate the structural disturbance impact partition layer.

[0148] The response path extraction module sets up a cross-sectional measurement point array based on the grid area within the boundary covered by the structural disturbance impact partition layer and collects groundwater level time series data at the corresponding locations. It then selects line segments with consistent change trends to obtain a groundwater response path distribution set.

[0149] The interference mechanism identification module determines the spatial overlap between the coordinates of the corresponding path segments in the groundwater response path distribution set and the boundary segments of the disturbance limit distribution morphological values, marks whether the path turns or changes its extension direction at the intersection, and obtains a path structure interference interactive identification list;

[0150] The recharge capacity classification module calls the path structure interference interactive identification list to extract the precipitation data, evaporation data and elevation data of the unit, divides the recharge level according to the combination level of the calculated results, and generates a spatial recharge capacity level layer;

[0151] The water supply sensitive area identification module superimposes the regional annual groundwater level change data based on the spatial recharge capacity level layer, calculates the joint value of the water level response amplitude and fluctuation frequency corresponding to the interference unit in the recharge level area, and screens out units whose joint value exceeds the distribution mean of the level area to obtain the water supply safety sensitive distribution area layer.

[0152] The system of this embodiment can be used to perform Figure 1 The technical solution of the method embodiment shown has similar implementation principles and technical effects, which will not be repeated here.

[0153] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. A method for assessing the safety of groundwater source water supply, characterized in that: The following steps are involved: S1: Obtain the permeability distribution map of the area covered by the riverside water source, divide the image into equal-area grids, analyze the permeability changes in the vertical and horizontal directions of the units, construct an edge point set based on the positions of continuous jump points, extract the closed boundary shape, and generate a structural disturbance impact partition layer; S2: calling the grid area in the structural disturbance impact partition layer, evaluating the consistency of the water level change direction and fluctuation time per unit time between the measuring points, screening the response path segments according to the degree of coherence, and obtaining the groundwater response path distribution set; S3: Based on the groundwater response path distribution set, spatial matching is performed on the overlapping areas between the path segments and the boundary segments to determine whether the path line crosses the boundary or turns along the boundary, and the position of the path direction mutation position and the degree of overlap with the fracture boundary are marked to establish a path structure interference interactive identification list; S4: Extracting precipitation data, evaporation data, and elevation data corresponding to the unit based on the path structure interference interactive identification list, calculating the precipitation-evaporation ratio and combining it with the neighborhood height difference to form an ups and downs index value, classifying the values ​​according to the recharge performance, and generating a spatial recharge capacity level layer; wherein the neighborhood height difference is the neighborhood average height difference; the steps for calculating the neighborhood average height difference include: According to the distribution data of the drop-evaporation ratio, the elevation value of the corresponding pixel position in the clipped elevation raster data is called, the elevation difference between the real-time pixel and the neighboring pixel is compared, and the average elevation difference of the neighborhood is calculated; S5: calling the spatial recharge capacity level layer and the marked areas in the path structure interference interactive identification list, superimposing the annual groundwater level fluctuation data in the area, extracting regional units with structural interference and recharge instability, analyzing the response amplitude and frequency characteristics, and obtaining a water supply security sensitive distribution area layer; The water supply safety sensitive distribution area layer includes response abnormality unit mark, water level fluctuation frequency index, and structural recharge overlapping area code.

2. The groundwater source water supply safety assessment method according to claim 1, characterized in that: The structural disturbance impact zoning layer includes a closed boundary line network, a disturbance boundary surface unit, and a continuous jump direction aggregation area; the groundwater response path distribution set includes a cross-section path line, a response direction marking point, and a water level continuity record index; the path structure interference interactive identification list includes a crossing path code, a junction change position identifier, and a recharge interruption area index; the spatial recharge capacity level layer includes a undulation correction index grid, a level zoning code map, and a recharge performance distribution layer.

3. The groundwater source water supply safety assessment method according to claim 1, characterized in that: The steps for obtaining the structural disturbance impact partition layer are specifically as follows: S101: Obtain a permeability distribution map of the area covered by the riverside water source, divide it into equal-area grids, treat the grid cells as independent spatial units, extract the initial soil permeability values ​​within the cell range, match the soil attribute point coordinates with the corresponding grid numbers, and generate a grid permeability distribution; S102: Calling the permeability values ​​of adjacent grid cells in the grid permeability distribution, calculating the permeability change difference between the real-time cell and the adjacent cells in the horizontal and vertical directions, and performing a determination operation on each pair of direction groups to determine whether there is a direction jump, thereby generating a direction jump distribution point set; S103: Based on the continuous spatial distribution of the jump points in the directional jump distribution point set, identify the jump point groups that meet the continuous quantity requirement, use the edge coordinates of the qualified point groups as edge candidate points, sort the edge candidate point coordinates by spatial proximity, extract the set of circumscribed boundary lines that constitute the closed area, and generate a structural disturbance impact partition layer.

4. The groundwater source water supply safety assessment method according to claim 3, characterized in that: The steps for obtaining the groundwater response path distribution set are specifically as follows: S201: calling the grid area in the structural disturbance impact partition layer, laying out a cross-section measuring point array in each group of grids according to a preset cross-section spacing, collecting continuous time water level change values ​​of the measuring points, organizing the water level time series data in the cross section according to the measuring point numbers, statistically analyzing the trend of the water level change direction and fluctuation amplitude between the measuring points in unit time, performing synchronization judgment on whether the direction and trend changes between the measuring point pairs in the cross section are consistent, and generating a consistency assessment result of the measuring point water level change; S202: Based on the consistency evaluation results of the water level changes at the measuring points, the measuring point pairs that meet the consistency requirements are screened in the cross section, connected paths are constructed using the connecting lines of the measuring point pairs as line segments, connected distance correction values ​​are calculated, and the spatial direction distributions of the multiple connected paths are integrated to generate a groundwater response path distribution set.

5. The method for assessing the safety of groundwater source water supply according to claim 4, characterized in that: The formula for calculating the connected distance correction value is: ； Among them, L ij Represents the connected distance correction value between measuring point i and measuring point j, x i 、x j Represents the horizontal coordinates of measuring point i and measuring point j, y i 、y j They represent the vertical coordinates of measuring points i and j respectively, n represents the total number of observation times, Δh ijk It represents the water level difference between measuring point i and measuring point j at the kth moment, It represents the average value of the water level at measuring point i and measuring point j during the entire period, and hc represents the centralized reference value of the water level at the cross-section measuring point during the entire period.

6. The groundwater source water supply safety assessment method according to claim 4, characterized in that: The steps for obtaining the path structure interference interaction identification list are specifically as follows: S301: Calling the spatial coordinates of the path segments in the groundwater response path distribution set, and performing spatial matching with the coordinate set of the boundary segments in the structural disturbance impact partition layer, determining whether the path segments have intersections or overlaps with the boundary segments in the coordinate range, if there is an intersection, marking it as crossing, if there is an overlap, marking it as turning along the boundary, and generating a path line and boundary matching type dataset; S302: Extracting the path direction change angle based on the path segments marked as crossing or turning in the path line and boundary matching type dataset, calculating the direction mutation value based on the angle difference between the previous and next segments of the path, and using the boundary segment coordinates to perform spatial overlap judgment to generate a path boundary mutation overlap annotation result; S303: Calling the mutation path segment number in the path boundary mutation coincidence degree annotation result, marking the path segment with disturbance signs, summarizing the abnormal marked areas at the partition boundary according to the path number, extracting the corresponding path segment and its location as the identification object, calculating the path segment disturbance identification value, and generating a path structure interference interactive identification list; The formula for calculating the path segment disturbance identification value is: ； Among them, D a is the path segment disturbance identification value, N is the total number of path segments, S ao is the matching score between the oth path segment and the ath path segment, M o is the average value of path segment o, σ o is the standard deviation of path segment o, α is the adjustment coefficient, W o is the path segment weight coefficient.

7. The groundwater source water supply safety assessment method according to claim 6, characterized in that: The steps for obtaining the space replenishment capability level layer are as follows: S401: Based on the unit numbers listed in the path structure interference interactive identification list, precipitation raster data, evaporation raster data, and elevation raster data of the corresponding unit area are extracted, the spatial resolution and projection format of the three types of data are unified, and the raster data are cropped according to the unit boundary range to obtain precipitation-evaporation ratio distribution data; S402: Based on the drop-evaporation ratio distribution data, call the elevation value of the corresponding pixel position in the clipped elevation raster data, compare the elevation difference between the real-time pixel and the neighboring pixels, calculate the average neighborhood elevation difference, and multiply the average elevation difference by the drop-evaporation ratio corresponding to the pixel to obtain the relief index distribution value; S403: Call the pixel data of the relief index distribution value, divide it into levels according to the numerical range, use the set relief index level threshold range as the classification basis, assign differentiated level values ​​respectively, and generate a level grid layer according to the spatial position of the division result to obtain the spatial replenishment capacity level layer.

8. The groundwater source water supply safety assessment method according to claim 1, characterized in that: The steps for obtaining the water supply safety sensitive distribution area layer are as follows: S501: The spatial recharge capacity level layer is called to identify the marked area in the path structure interference interactive identification list. For the recharge units in the marked area, the annual groundwater level monitoring data of the corresponding location is collected. The groundwater level values ​​at the same location are arranged in a monthly sequence. The water level fluctuation amplitude of the recharge unit in the year is calculated, and the recharge units that do not fluctuate are screened out to obtain a groundwater level fluctuation amplitude dataset. S502: Based on the groundwater level fluctuation amplitude dataset and in combination with the regional interference level zoning results in the path structure interference interactive identification list, the fluctuation amplitude value change frequency of the recharge units within the interference level area is counted, and the recharge units whose change frequency exceeds the regional annual change benchmark are screened to obtain the interference area unit sequence; S503: According to the interference area unit sequence, the spatial position index value of the corresponding unit is extracted, and combined with the grade division of the corresponding position in the spatial replenishment capacity grade layer, the area whose distribution ratio exceeds the internal benchmark level of the grade area to which it belongs is screened through index matching operation to generate a water supply safety sensitive distribution area layer.

9. A groundwater source water supply safety assessment system, characterized in that: The system is used to implement the groundwater source water supply safety assessment method according to any one of claims 1 to 8, and the system includes: The structural disturbance identification module obtains the permeability distribution map of the area covered by the riverside water source. After dividing the image into equal-area grids, the permeability values ​​within the grids are called in the vertical and horizontal directions. After forming unit pairs, the permeability differences between adjacent grids in the same direction are compared. The closed boundaries are connected in sequence and the spatial contours are extracted to generate the structural disturbance impact partition layer. A response path extraction module sets a cross-sectional measurement point array based on the grid area within the boundary covered by the structural disturbance impact partition layer, collects groundwater level time series data at the corresponding locations, and selects line segments with coherent change trends to obtain a groundwater response path distribution set; The interference mechanism identification module performs spatial overlap judgment with the boundary line segment of the disturbance limit distribution morphological value based on the coordinates of the corresponding path segment in the groundwater response path distribution set, marks whether the path turns or changes its extension direction at the intersection, and obtains a path structure interference interactive identification list; The recharge capacity classification module calls the path structure interference interactive identification list, extracts the precipitation data, evaporation data and elevation data of the unit, divides the recharge level according to the combination level of the calculation results, and generates a spatial recharge capacity level layer; The water supply sensitive area identification module, based on the spatial recharge capacity level layer, superimposes the regional annual groundwater level change data, calculates the joint value of the water level response amplitude and fluctuation frequency corresponding to the interference unit in the recharge level area, and screens out the units whose joint value exceeds the distribution mean of the level area to obtain the water supply safety sensitive distribution area layer.

Images