Underground water source water supply safety assessment method and system

Abstract

The invention provides an underground water source water supply safety assessment method and system, and relates to the technical field of water supply safety assessment. The method comprises the steps of obtaining a permeability distribution diagram, dividing grids and analyzing direction changes, constructing an edge set to extract a closed form to generate a partition layer, setting measuring points to collect water level data, evaluating a change consistency screening path, matching a boundary to judge a crossing condition, marking a sudden change position, recording an interference path and establishing a recognition list. According to the method, by comparing the consistency of the water level direction in unit time and the response fluctuation, the path line segment with the hydraulic coherence characteristic is identified, the depicting ability of the actual underground water flowing path in the non-uniform supply environment is improved, the path trend sudden change and boundary overlapping position is analyzed, the response identification of the structure interference action area is enhanced, and the accuracy of the structure interference action area is improved. A key section with water supply interruption or supply abnormity is positioned, and the response capability of actual hydrological and geological changes and the spatial fineness of evaluation and judgment are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of water supply safety assessment, and particularly 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 includes quantitative or qualitative analysis and assessment of elements such as water quality and quantity, water source recharge, groundwater flow path, pollutant migration, and human activity impact of riverside groundwater source areas. The core content of this technical field lies in constructing a systematic assessment framework to scientifically analyze the hydrogeological structure, hydrogeological recharge relationship, hydraulic connection with surface water bodies, and susceptibility to pollution risk of the water source area, so as to determine its water supply stability and safety under different hydrological conditions. The overall technical field uses means such as on-site investigation, hydrogeological monitoring, experimental simulation, and numerical simulation to comprehensively analyze the recharge source, runoff process of groundwater resources, and their response characteristics to external disturbances, providing a scientific basis for the safety management of riverside water sources in urban water supply systems.

[0003] Among them, the method for assessing the safety of groundwater source water supply refers to the ability to supply water safely under different conditions for riverside water source areas. By establishing an analysis method based on the interaction relationship between groundwater and surface water, combined with hydrogeological profiles, hydrological monitoring data, and multi-temporal groundwater level changes, the reliability and risk of its water supply are systematically evaluated. This patent theme covers the calculation of groundwater runoff direction and intensity, the deduction of pollutant migration paths and speeds, the determination of the boundaries of water source recharge areas, and the dynamic response analysis of rainfall and surface water changes to the groundwater system. The assessment method is specifically completed by means such as hydrological monitoring data processing, determination of hydrogeological parameters in borehole profiles, calculation of hydrodynamic gradients, and determination of groundwater-surface water exchange zone parameters.

[0004] In the prior art, the identification of underground structures mainly relies on the interpretation of hydrogeological profiles and the judgment of the linear trend of geological maps, ignoring local continuity disturbances under atypical strike or small-scale fracture conditions, resulting in a decrease in the accuracy of fracture boundary identification. In terms of hydraulic path identification, it mostly relies on the deduction of regional water level contour lines, making it difficult to depict the micro-variation trend of hydraulic paths under non-uniform recharge conditions. In the determination of recharge capacity, the recharge state is often represented by planar means, lacking the distinction of the differences in recharge capacity performance under micro-topographic undulations, resulting in a low match with the actual recharge distribution. Existing methods mostly separately evaluate the groundwater fluctuation value and recharge capacity, failing to comprehensively identify sensitive areas with high fluctuation response and insufficient recharge at the same time. Under the condition where the strong evaporation zone overlaps with the shallow groundwater level area, although the water level fluctuates frequently, it does not have continuous recharge capacity, and existing systems are difficult to identify such high-risk areas. These deficiencies lead to the omission of sensitive points in water supply safety assessment and misjudgment of the spatial distribution of recharge capacity, affecting the scientific nature of regional water source allocation and scheduling strategies. Summary of the Invention

[0005] In order to solve the technical problems existing in the prior art, such as the neglect of local continuity disturbances under non - typical trends or small - scale fracture conditions, resulting in a decline in the accuracy of fracture boundary identification; in the aspect of hydraulic path identification, relying mostly on the deduction of regional water level contour lines, it is difficult to depict the micro - variable trend of hydraulic paths under non - uniform recharge conditions; in the determination of recharge capacity, often using the planar average value to represent the recharge state, lacking the distinction of the differences in recharge capacity performance under micro - topographic undulations, resulting in a low matching of the actual recharge distribution; most existing methods separately evaluate the groundwater fluctuation value and recharge capacity, failing to comprehensively identify sensitive areas with high - fluctuation response and insufficient recharge at the same time; under the condition of the overlap of the strong evaporation zone and the shallow groundwater level zone, although the water level fluctuates frequently but does not have continuous recharge capacity, the existing system is difficult to identify such high - risk areas, leading to the omission of sensitive points in the water supply safety assessment and misjudgment of the spatial distribution of recharge capacity, affecting the scientific nature of regional water source allocation and scheduling strategies, the present invention provides a method and system for evaluating the water supply safety of groundwater source areas. The technical solutions are as follows:

[0006] On the one hand, a method for evaluating the water supply safety of a groundwater source area is provided, and the method includes:

[0007] S1: Obtain the permeability distribution map of the area covered by the river - side water source area, divide the image into equal - area grids, analyze the permeability change information in the longitudinal and transverse directions of the unit, construct an edge point set according to the position of continuous jump points and extract the closed - boundary shape, and generate a structural disturbance influence zoning layer;

[0008] S2: Call the grid area in the structural disturbance influence zoning layer, evaluate the consistency of the water level change direction and fluctuation time per unit time between measurement points, screen the response path line segments according to the coherence degree, and obtain the groundwater response path distribution set;

[0009] S3: According to the groundwater response path distribution set, perform spatial matching on the overlapping area between the path line segment and the boundary line segment, judge whether the path line crosses the boundary or turns along the boundary, mark the overlapping degree between the path - trend mutation position and the fracture boundary, and establish a path - structure interference interaction identification list;

[0010] S4: According to the path - structure interference interaction identification list, extract the corresponding precipitation data, evaporation data and elevation data of the unit, calculate the precipitation - evaporation ratio and combine it with the elevation difference of adjacent grids to form a relief index value, divide the grades according to the recharge performance, and generate a spatial recharge capacity grade layer.

[0011] Optionally, the structure disturbance impact zoning layer includes a closed boundary line network, disturbance boundary surface units, and continuous jump direction aggregation areas. The groundwater response path distribution set includes cross-sectional path lines, response direction marking points, and water level coherence record indexes. The path structure interference interaction identification list includes crossing path codes, junction direction change position identifiers, 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.

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

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

[0014] S102: Call the permeability values of adjacent grid cells in the grid permeability distribution, calculate the difference in permeability change between the real-time cell and the adjacent cell in the horizontal and vertical directions respectively, and perform a judgment operation on whether there is a direction jump for each pair of direction pairs to generate a set of direction jump distribution points;

[0015] S103: According to the continuous spatial distribution of the jump points in the direction jump distribution point set, identify the jump point groups that meet the continuous quantity requirements, use the edge coordinates of the qualified point groups as edge candidate points, perform a closure sorting on the edge candidate point coordinates according to spatial adjacency, and extract the set of circumscribed boundary lines that form a closed area to generate the structure disturbance impact zoning layer.

[0016] Optionally, the steps for obtaining the groundwater response path distribution set are specifically as follows:

[0017] S201: Call the grid areas in the structure disturbance impact zoning layer, arrange cross-sectional measurement point arrays in each group of grids at a preset cross-sectional spacing, collect the continuous time water level change values of the measurement points, organize the water level time series data within the cross-section according to the measurement point numbers, count the change trends of the water level change direction and fluctuation amplitude between the measurement points within the unit time, and perform a synchronization judgment on whether the direction and trend changes between the measurement point pairs within the cross-section are consistent to generate an evaluation result of the consistency of the water level change of the measurement points;

[0018] S202: According to the evaluation result of the consistency of the water level change of the measurement points, screen the measurement point pairs that meet the consistency requirements within the cross-section, construct a connected path with the connection line of the measurement point pairs as a line segment, calculate the corrected value of the connected distance, and integrate the spatial trend distributions of multiple connected paths to generate the groundwater response path distribution set.

[0019] Optionally, the formula for calculating the corrected value of the connected distance is as follows:

[0020] ;

[0021] where L ij represents the corrected value of the connected distance between measuring point i and measuring point j, x i , x j respectively represent the abscissas of measuring point i and measuring point j, y i , y j respectively represent the ordinates of measuring point i and measuring point j, n represents the total number of observation times, represents the water level difference between measuring point i and measuring point j at the kth moment, represents the average value of the average water levels of measuring point i and measuring point j over the entire period, respectively represent the average water levels of measuring point i and j over the entire period, and hc represents the centralized reference value of the water level of the cross-section measuring points over the entire period.

[0022] Optionally, the steps for obtaining the list of path structure interference interaction identifications are specifically as follows:

[0023] S301: Call the spatial coordinates of the path segments in the groundwater response path distribution set, perform spatial matching with the coordinate set of the boundary segments in the structure disturbance influence zoning layer, determine whether there are intersections or overlapping segments between the path segments and the boundary segments in terms of coordinate range. If there are intersections, mark them as crossing; if there are overlapping segments, mark them as turning along the boundary, and generate a path line and boundary matching type data set;

[0024] S302: According to the path segments marked as crossing or turning in the path line and boundary matching type data set, extract the path direction change angle, calculate the degree of sudden change in the direction according to the angle difference between the front and rear segments of the path, and call the boundary segment coordinates for spatial coincidence degree judgment to generate a path boundary sudden change coincidence degree annotation result;

[0025] S303: Call the numbers of the sudden change path segments in the path boundary sudden change coincidence degree annotation result, mark the path segments with signs of disturbance, summarize the abnormal annotation areas of the partition boundaries where they are located according to the path numbers, extract the corresponding path segments and their locations as the identification objects, calculate the disturbance identification value of the path segments, and generate a list of path structure interference interaction identifications;

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

[0027] ;

[0028] where D a is the disturbance identification value of the path segment, N is the total number of path segments, S ao is the matching score between the oth path segment and the a-th 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 spatial recharge capacity level layer are specifically as follows:

[0030] S401: Based on the unit numbers listed in the path structure interference interaction recognition list, extract the precipitation grid data, evaporation grid data, and elevation grid data of the corresponding unit areas. After unifying the spatial resolutions and projection formats of the three types of data, perform grid clipping according to the unit boundary range to obtain the precipitation-evaporation ratio distribution data;

[0031] S402: According to the precipitation-evaporation ratio distribution data, call the elevation values at the corresponding pixel positions in the clipped elevation grid data, compare the elevation differences between the real-time pixel and the neighborhood pixels, calculate the average elevation difference of the neighborhood, and perform a product calculation on the average elevation difference value and the precipitation-evaporation ratio corresponding to the pixel to obtain the undulation index distribution value;

[0032] S403: Call the pixel data of the undulation index distribution value, perform grade division according to the numerical interval, use the set undulation index grade threshold interval as the classification benchmark, assign different grade values respectively, and generate a grade grid layer according to the spatial position of the division result to obtain the spatial recharge capacity level layer.

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

[0034] S5: Call the spatial recharge capacity level layer and the marked area in the path structure interference interaction recognition list, overlay the annual fluctuation data of the groundwater level within the area, extract the area units with structural interference and unstable recharge, analyze the response amplitude and frequency characteristics, and obtain the water supply security sensitive distribution area layer;

[0035] The water supply security sensitive distribution area layer includes response anomaly unit marks, water level fluctuation frequency indexes, and structural recharge overlay area codes.

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

[0037] S501: Call the spatial recharge capacity level layer and the marked area in the path structure interference interaction recognition list. For the recharge units within the marked area, collect the annual groundwater level monitoring data at the corresponding positions, arrange the groundwater level values at the same position in monthly sequence, calculate the water level change amplitude of the recharge units within the year, and screen out the recharge units without fluctuations to obtain the groundwater level fluctuation amplitude data set;

[0038] S502: Based on the dataset of the fluctuation amplitude of the groundwater level, combined with the regional interference level zoning results in the list of path structure interference interaction identification, count the change frequency of the fluctuation amplitude values of recharge units within the interference level area, and screen out the recharge units whose change frequency exceeds the annual change benchmark of the area to obtain a sequence of interference area units;

[0039] S503: According to the sequence of interference area units, extract the spatial position index values of the corresponding units, combine with the grade zoning at the corresponding positions in the spatial recharge capacity level layer, and through index matching operations, screen out the areas where the distribution ratio exceeds the internal benchmark level of the corresponding grade area 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 performing the above-mentioned groundwater source water supply safety assessment method, and the system includes:

[0041] A structural disturbance identification module, which obtains the permeability distribution map of the area covered by the river-side water source, divides the image into equal-area grids, then calls the permeability values within the grids in the vertical and horizontal directions respectively, constructs unit pairs, compares the permeability differences between adjacent grids in the same direction, connects them in sequence to form a closed boundary and extracts the spatial contour graph, and generates a structural disturbance impact zoning layer;

[0042] A response path extraction module, based on the grid area within the boundary covered by the structural disturbance impact zoning layer, sets up a cross-section measuring point array and collects the groundwater level time series data at the corresponding positions, and screens out the line segments with consistent change trends to obtain a set of groundwater response path distributions;

[0043] An interference mechanism identification module, according to the coordinates of the corresponding path line segments in the set of groundwater response path distributions, makes a spatial overlap judgment with the boundary line segments of the disturbance limit distribution form value, marks whether the path turns or changes the extension direction at the junction, and obtains a list of path structure interference interaction identifications;

[0044] A recharge capacity division module, calls the list of path structure interference interaction identifications, extracts the precipitation data, evaporation data and elevation data of the corresponding units, divides the recharge grade according to the combined level of the calculation results, and generates a spatial recharge capacity level layer;

[0045] A water supply sensitive area identification module, according to the spatial recharge capacity level layer, superimposes the annual groundwater level change data of the area, calculates the combined value of the water level response amplitude and fluctuation frequency corresponding to the interference units in the recharge grade area, and screens out the units whose combined value exceeds the distribution mean of the grade area to obtain a water supply safety sensitive distribution area layer.

[0046] The beneficial effects brought by the technical solutions provided in the embodiments of the present invention at least include:

[0047] In the embodiment of the present invention, by judging the jump characteristics of the permeability change direction, the closed boundary morphology is extracted with the help of the spatial aggregation of edge points, and the geological boundary indication or empirical judgment is no longer relied on in the structural identification stage, so as to enhance the resolution ability of spatial continuous disturbance identification. In path identification, the path line segment with hydraulic coherence characteristics is identified by comparing the consistency of the water level direction per unit time with the response fluctuation, so as to improve the ability to describe the actual groundwater flow path under the non-uniform recharge environment. In the path and boundary interaction stage, the sudden change of the path direction and the coincidence position of the boundary are analyzed by spatial superposition, and the response identification of the structural interference area is strengthened, which can effectively locate the key section of water supply interruption or recharge anomaly. In the identification of recharge capacity, the fluctuation index participates in the correction of the precipitation-evaporation relationship, and the local terrain characteristics are linked with the recharge performance to distinguish the differences in the ability of different geomorphic units to compensate for water supply. By superimposing the annual fluctuation of groundwater level, a composite judgment standard between recharge performance and fluctuation amplitude is established to achieve high recognition accuracy of 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, connecting each link of underground structure identification, hydraulic path identification, recharge stability quantification and risk spatial expression, thus improving the overall response capability of the assessment system to actual hydrological and geological changes and the spatial precision of assessment judgment. 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 "example" in the present invention should not be interpreted as being more preferred or more advantageous than other embodiments or designs. Specifically, the use of the word "example" is intended to present the concept in a specific way. In addition, in the embodiments of the present invention, the meaning expressed by "and / or" can be both, or it can be either of the two.

[0052] In order to make the technical problems, technical solutions and advantages to be solved by the present invention more clear, 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 river-side water source area, divide the image into equal-area grids, analyze the permeability change information in the longitudinal and transverse directions of the analysis unit, construct direction pairs with adjacent units and then judge the direction jump characteristics, construct an edge point set according to the positions of continuous jump points and extract the closed boundary form, and generate a structural disturbance impact zoning layer;

[0055] S2: Call the grid areas in the structural disturbance impact zoning layer, set up cross-section measuring point arrays and collect water level time series data, evaluate the consistency of the water level change direction and the fluctuation time per unit time between measuring points, screen the response path line segments according to the coherence degree, and obtain the groundwater response path distribution set;

[0056] S3: According to the groundwater response path distribution set, perform spatial matching on the overlapping areas between the path line segments and the boundary line segments, judge whether the path line crosses the boundary or turns along the boundary, mark the overlapping degree between the path direction mutation position and the fracture boundary, record the path line segments with direction changes as potential interference paths, and mark the recharge anomaly areas in combination with the continuous change trend of the path in the junction area, and establish a path structure interference interaction recognition list;

[0057] S4: According to the path structure interference interaction recognition list, extract the corresponding precipitation data, evaporation data and elevation data of the unit, calculate the precipitation-evaporation ratio and then form a relief index value in combination with the elevation difference between adjacent grids, divide into grades according to the recharge performance, and generate a spatial recharge capacity grade layer;

[0058] S5: Call the areas marked in the spatial recharge capacity grade layer and the path structure interference interaction recognition list, overlay the annual fluctuations of the groundwater level in the area, extract the area units with structural interference and unstable recharge, analyze the response amplitude and frequency characteristics, and obtain the water supply safety sensitive distribution area layer.

[0059] Among them, the structural disturbance impact zoning layer includes a closed boundary line network, a disturbed boundary surface unit, and a continuous jump direction aggregation area; the groundwater response path distribution set includes cross-section path lines, response direction marking points, and water level coherence record indexes; the path structure interference interaction recognition list includes crossing path codes, junction direction change position identifiers, and recharge interruption area indexes; the spatial recharge capacity grade layer includes relief correction index grids, grade zoning coding maps, and recharge performance distribution layers; the water supply safety sensitive distribution area layer includes response anomaly unit marks, water level fluctuation frequency indexes, and structural recharge overlay area coding.

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

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

[0062] It is necessary to use geographic information software to load the remote sensing image or topographic map of the river-side water source area, determine the boundary range of the area through manual delineation or automatic recognition, divide it into equal-area grids at a fixed scale within this boundary range, set each grid as a square unit of 100 meters × 100 meters, 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 of which contains geographical coordinates and the original permeability values. Match the points with the grids through spatial matching. If a certain soil point is within the unit numbered A25, then assign its value to the A25 unit. If a unit contains multiple points, it is necessary to extract the original permeability values of the points, perform statistical integration on them, and the representative value of the unit can be obtained by using the average method or the median method. Suppose there are 3 points in a unit, located in the upper left corner, middle, and lower right corner respectively, with permeabilities of 2.3, 2.5, and 2.7. After integration, the representative value of this unit can be set as 2.5. At the same time, exclude the units without point data, and the value can be estimated by using spatial interpolation methods, such as referring to the data of adjacent units around. After completion, integrate the unit numbers and their permeabilities into a distribution map to generate the grid permeability distribution.

[0063] S102: Call the permeability values of adjacent grid cells in the grid permeability distribution, calculate the difference in permeability change between the real-time unit and the adjacent unit in the horizontal and vertical directions respectively, and perform a judgment operation on whether there is a direction jump for each pair of direction pairs to generate a set of direction jump distribution points;

[0064] Extract the adjacent cell data of each grid cell, compare them separately in the horizontal and vertical directions, set the permeability values for comparing the left and right, upper and lower cells, calculate and record the differences in the two directions separately. If a certain cell shows a large difference in permeability between its right adjacent cell in the horizontal comparison and a small change in its lower adjacent cell in the vertical comparison, the horizontal difference value will be marked as key jump data. Further, by setting a jump judgment criterion, determine whether this difference meets the condition for judging a jump. If it does, record the current cell as a jump point, record its spatial position, and add it to the jump point dataset. This process needs to be executed for all grid cells in the entire region. Suppose a region contains 400 grid cells, then all the grids will be traversed, the change relationships in its four directions will be analyzed, and the permeability differences for each pair of directions will be identified and classified to obtain the jump point layer data. Each point in this layer represents a jump position, including its grid number and direction difference characteristics, which are used for subsequent jump point group identification and structural disturbance analysis to generate a direction jump distribution point set.

[0065] S103: According to the continuous spatial distribution of the jump points in the direction jump distribution point set, identify the jump point groups that meet the continuous quantity requirements, use the edge coordinates of the qualified point groups as the edge candidate points, perform a closed sorting on the edge candidate point coordinates according to spatial adjacency, extract the set of circumscribed boundary lines forming the closed region, and generate a structural disturbance impact zoning layer;

[0066] Perform a spatial analysis on the jump point layer. According to the distribution density and continuity of the jump points in space, extract the jump point groups that meet the set conditions. Set that the set continuous jump point group should contain at least 5 adjacent points. Successively check whether the spatial distances between the jump points and their surrounding points meet the proximity conditions. If they do, classify them into the same jump point group. After the point group identification is completed, extract the edge points of each group of point groups, and use the boundary sorting algorithm 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. Further, draw the set of boundary lines of the region enclosed by this closed path. Suppose a certain jump point group is distributed in the southwestern region of the riverbank, and its outer edge points are arranged in the order of southeast - northeast - northwest - southwest, then their connection will form a convex quadrilateral or polygon. Export this set of boundary lines as the spatial layer of the structural disturbance impact zoning and overlay it 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 as follows:

[0068] S201: Call the grid area in the structure perturbation impact partition layer, arrange a cross-section measurement point array in each group of grids at a preset cross-section spacing, collect the continuous time water level change values of the measurement points, organize the water level time series data in the cross-section according to the measurement point numbers, statistically analyze the change trends of the water level change directions and fluctuation amplitudes between the measurement points per unit time, perform a synchronization judgment on whether the directions and trend changes between the measurement points in the cross-section are consistent, and generate an evaluation result of the consistency of the water level changes of the measurement points;

[0069] Perform the arrangement operation according to the preset cross-section spacing. The cross-section spacing is usually set according to the on-site terrain features, hydrogeological conditions, and survey accuracy requirements. If it is set to arrange a cross-section every 200 meters, then the grid lines crossing the perturbation partition will be automatically selected as cross-sections within this spacing interval, and several measurement points will be evenly arranged on each cross-section. The measurement point spacing can be set to 20 meters to ensure coverage of the entire width of the cross-section. The naming of each measurement point needs to be encoded to include the cross-section number and the point position sequence. Set the second measurement point on the 3rd cross-section to be named D3P2. An underground water level observation device or a remote sensing induction node is arranged at each measurement point, and the water level change data is collected regularly or continuously. The data recording form is a time series list. According to the time dimension, the data of different measurement points is sorted in a unified format. For points such as D3P2 and D3P3, a continuous time-water level pair data list is generated. The water level change trend between each measurement point is statistically analyzed per unit time step. The difference change between adjacent moments can be used to represent the rising or falling direction. If both measurement points show an upward trend within a 5-minute time period, it is recorded that the directions are the same. At the same time, it is statistically analyzed whether the change amplitudes between each pair of measurement points are within the same level interval. Set the boundary values of the fluctuation amplitude levels (such as 2 cm, 5 cm, 5 - 10 cm, etc.), and judge whether they belong to the same level. Based on this, the water level fluctuation synchronization between the measurement point pairs is judged. After the judgment of each measurement point pair, each pair of measurement points is marked whether they are consistent, and an attribute table is formed, listing the measurement point pair number, the direction consistency flag, and the fluctuation amplitude level consistency flag, as the basic data for subsequent path construction, and generating an evaluation result of the consistency of the water level changes of the measurement points.

[0070] S202: According to the evaluation result of the consistency of the water level changes of the measurement points, screen the measurement point pairs in the cross-section that meet the consistency requirements, construct a connected path with the connection line of the measurement point pairs as a line segment, calculate the correction value of the connected distance, and integrate the spatial orientation distributions of multiple connected paths to generate a set of groundwater response path distributions;

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

[0072] ;

[0073] Among them, L ij represents the correction value of the connected distance between measurement point i and measurement point j, x i , x jRepresent the abscissas of measurement points i and j respectively, y i and y j Represent the ordinates of measurement points i and j respectively, and n represents the total number of observation times. Indicates the water level difference between measurement points i and j at the k-th moment. Represents the average of the average water levels of measurement points i and j over the entire period. Represent the average water levels of measurement points i and j over the entire period respectively, and hc represents the centralized reference value of the water level of the cross-section measurement points over the entire period.

[0074] Meaning of parameters and derivation process of formula calculation:

[0075] Calculation of the coordinate difference of measurement points:

[0076] ;

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

[0078] x i and y i Are the coordinate values of measurement point i, and x j and y j Are the coordinate values of measurement point j;

[0079] The coordinate of measurement point i is (10, 5), and the coordinate of measurement point j is (15, 8), then:

[0080] ;

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

[0082] Average absolute difference of water level difference:

[0083] ;

[0084] This item calculates the average absolute value of the water level difference between measurement points i and 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 measurement points i and j at the k-th moment,

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

[0088] Average absolute difference:

[0089] ;

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

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

[0092] ;

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

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

[0095] hc is the centralized reference value of the water levels of all measuring points in this cross-section over the entire period, such as the median or the overall mean;

[0096] Set the average water levels of measuring point i and measuring point j to be and , respectively, then:

[0097] ;

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

[0099] ;

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

[0101] Substitute into the formula for calculation:

[0102] ;

[0103] This result indicates that the corrected value of the connection path between measuring point i and measuring point j is 6.68 units. This value not only considers the straight-line distance in space but also combines the consistency of water level changes and the difference from the global water level to correct the connection path.

[0104] The specific steps for obtaining the list of path structure interference interactions 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 structure disturbance impact zoning layer. Determine whether there are intersections or overlapping segments between the path segments and the boundary segments in terms of coordinate range. If there are intersections, mark them as crossing; if there are overlapping segments, mark them as turning along the boundary to generate a data set of the matching type between the path line and the boundary;

[0106] The path segments are successively subjected to spatial overlapping matching with the boundary segments in the structural disturbance influence zoning layer. The specific operation process is as follows: Read the starting and ending coordinates of the path segments, construct a continuous line segment chain, number and index each path one by one, extract the boundary segment set from the disturbance zoning boundary layer, establish a spatial index for each boundary line, judge whether there is a potential contact probability between the path and the boundary through the fast bounding box method of the coordinate range. After screening out the contact candidate pairs, perform a more accurate line segment intersection detection to judge whether there is an intersection point. If the detection result is that the path segment intersects with the boundary segment and the intersection point falls within the line segment range, mark the path as a crossing type match. If the detection result is that the path segment partially or completely coincides with the boundary segment, especially when the continuous match exceeds the set length threshold under multi-point precision coordinates, mark it as a boundary-folding type match. Set a path number P21, whose endpoint coordinates range from (135.221, 36.447) to (135.245, 36.472). If it generates an intersection point when passing through the boundary segment numbered B12, record the intersection point information in the P21 attribute and assign a crossing label. At the same time, record the path numbers of all paths with intersection points or overlapping segments into the path boundary matching type data set. This data set includes key fields such as path number, matching type (crossing or folding), matching boundary number, intersection point coordinates or start and end positions of the overlapping segment, etc., as the basis for subsequent trend analysis and interference recognition, and generate the path line and boundary matching type data set.

[0107] S302: According to the path segments marked as crossing or folding in the path line and boundary matching type data set, extract the path trend change angle, calculate the trend mutation degree value according to the angle difference between the front and rear segments of the path, and call the boundary segment coordinates to judge the spatial coincidence degree, and generate the path boundary mutation coincidence degree 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 respectively, calculate the change amplitude of the direction angles to judge the degree of sudden change in the trend. The operation steps include: read the coordinates of the connecting line segments before and after the boundary intersection points of the path segments, obtain their direction vectors respectively, calculate the sudden change angle of the path trend according to the included angle formed by the vectors. If the angle change exceeds the preset threshold (such as 30 degrees), it is regarded as a mutation point, record its mutation degree value, and at the same time read the boundary line segment coordinates again at this position to judge whether there is a spatial coincidence between the path and the boundary. The proportion of the coincidence length to the boundary length or the path length can be used as the coincidence degree index. If the proportion exceeds 50%, it is marked as significant coincidence. Integrate the mutation angle and the coincidence degree index into the annotation result. Each record includes the path number, mutation angle, mutation flag, coincidence percentage value, mutation point coordinates and the corresponding boundary line segment number, etc. Set the path P31. If the direction angle difference between its front and rear segments is 48 degrees, and there is a 70% coincidence ratio between this mutation segment and the boundary B22, then this record will be given double labels of mutation and significant coincidence, providing a direct identification basis for the next step of identifying the disturbed path, and generating the annotation result of the path boundary mutation coincidence degree.

[0109] S303: Call the mutation path segment numbers in the annotation result of the path boundary mutation coincidence degree, mark the path segments with signs of disturbance, summarize the abnormal annotation areas of the boundaries of the corresponding partitions according to the path numbers, extract the corresponding path segments and their locations as the recognition objects, calculate the disturbance recognition value of the path segments, and generate the interactive recognition list of path structure interference;

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

[0111] ;

[0112] where, D a is the disturbance recognition value of the path segment, N is the total number of path segments, S ao is the matching score between the o-th path segment and the a-th path segment, M o is the average value of the path segment o, σ o is the standard deviation of the path segment o, α is the adjustment coefficient, W o is the weight coefficient of the path segment.

[0113] Meaning of parameters and derivation process of formula calculation:

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

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

[0116] Path segment average value (Mo ):

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

[0118] ;

[0119] Standard deviation of path segment (σ o ):

[0120] Calculate the standard deviation of path segment o, which represents the degree of dispersion of the eigenvalues of the path segment. Set the eigenvalues of path segment o as [0.8, 0.9, 0.85], then its standard deviation σ o is:

[0121] ;

[0122] Adjustment coefficient (α):

[0123] The adjustment coefficient is used to control the influence degree of the deviation after standardization, and is set according to actual requirements. Set α = 2;

[0124] Weight coefficient of path segment (W o ):

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

[0126] Substitute the above parameters into the formula:

[0127] ;

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

[0129] ;

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

[0131] This result shows that the matching degree between path segment a and path segment o is high, and the degree of dispersion of its eigenvalues is low. Therefore, there is no need to specifically mark the path segment with disturbance signs in the path structure.

[0132] The specific steps for obtaining the spatial supply capacity level layer are as follows:

[0133] S401: Based on the unit numbers listed in the path structure interference interaction recognition list, extract the precipitation grid data, evaporation grid data, and elevation grid data of the corresponding unit area. After unifying the spatial resolutions and projection formats of the three types of data, perform raster clipping according to the unit boundary range to obtain the distribution data of the precipitation-evaporation ratio;

[0134] Retrieve the spatial unit ranges corresponding to each number in sequence, locate its boundary contour through spatial indexing, and extract the precipitation grid data, evaporation grid data, and digital elevation model data within the coverage of this area from the remote sensing data platform or the pre-set raster database. Perform format conversion with a unified spatial reference (such as WGS84) and the same resolution (such as 30 meters) respectively to ensure that the projection formats and pixel sizes of the three data sources are consistent. After completion, perform raster clipping operations on the three types of data according to the unit boundary, that is, only retain the pixel area within the identified unit range. Set a unit number U107 to be located in the southwestern part of the layer, with the corresponding spatial range of 135.215 longitude and 36.455 latitude. The clipping operation extracts the precipitation, evaporation, and elevation data pixels within this range to form three groups of unified raster data. Perform ratio operations on the clipped precipitation data and evaporation data pixel by pixel to form a precipitation-evaporation ratio distribution layer. This ratio represents the water surplus and deficit situation at each pixel. Set the precipitation of a certain pixel to be 85 mm and the evaporation to be 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. Through ratio calculation, the basic data for the water dynamics evaluation of each unit is completed, providing input for the subsequent recharge capacity assessment, and obtaining the distribution data of the precipitation-evaporation ratio.

[0135] S402: According to the distribution data of the precipitation-evaporation ratio, call the elevation value at 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 height difference of the neighborhood, and perform a product calculation on the average height difference value and the precipitation-evaporation ratio corresponding to the pixel to obtain the distribution value of the undulation index;

[0136] Read the ratio results of each pixel, perform neighborhood analysis on the position of the pixel in the elevation raster data, extract the elevation values of adjacent pixels within its eight-neighborhood range, calculate the elevation difference between the current pixel and its neighborhood, average the elevation differences, and obtain the average neighborhood elevation difference of the pixel. This elevation difference can reflect the degree of terrain undulation. Set the elevation of a certain pixel to 310 meters, and the elevations of its surrounding pixels are 308, 312, 315, 309, 305, 307, 311, and 313 respectively. Then the average elevation difference is the average value of this set of differences. After obtaining the result, multiply the evapotranspiration reduction ratio of the pixel by its average elevation difference, which is called the undulation index distribution value, a comprehensive factor that combines terrain undulation and hydrological profit and loss ability. Perform this calculation process on each pixel in turn to form a complete undulation index layer. Set within a certain unit, the evapotranspiration reduction ratio of a certain pixel is 1.25, and the average neighborhood elevation difference is 6.2 meters. Then the undulation index value of this pixel is 7.75. After the pixel is processed, this index layer can be used to reveal the spatial recharge ability differentiation characteristics of the response path of terrain and hydrology interaction to groundwater, and obtain the undulation index distribution value.

[0137] S403: Call the pixel data of the undulation index distribution value, perform grade division according to the numerical interval, use the set undulation index grade threshold interval as the classification benchmark, assign different grade values respectively, and generate a grade raster layer according to the spatial position of the division result to obtain the spatial recharge ability grade layer;

[0138] Take the undulation index value of the pixel as the classification object, perform grade division on it according to the set grade threshold interval. The grade division rules can be preset as several grade sections according to the landform type and regional characteristics. Set the undulation index value to be divided into: low grade (05), medium-low grade (10), medium grade (15), medium-high grade (20), high grade (greater than 20). Assign corresponding grade identification values to all pixels according to the interval into which their values fall, such as grade 1 to grade 5. After completing the grade assignment of the pixels, generate a new grade raster layer according to their spatial positions, that is, each pixel corresponds to a grade identification in the layer, forming the spatial recharge ability grade layer. This layer can be directly loaded and browsed on the GIS platform, and can be superimposed on the path disturbance layer or the structure partition layer to form a linkage analysis result. Set that most path interaction paragraphs in a certain disturbance partition in the layer are located in the high-grade recharge interval, which indicates that this area has good water recharge ability. Output this grade layer for further analysis or visualization expression to obtain the spatial recharge ability grade layer.

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

[0140] S501: Call the marked areas in the list of interactive identification of the spatial supply capacity level layer and the path structure interference, collect the annual groundwater level monitoring data at the corresponding positions for the supply units within the marked areas, arrange the groundwater level values at the same position in monthly sequence, calculate the water level variation range of the supply units within the year, and screen out the supply units without fluctuations to obtain the groundwater level fluctuation range data set;

[0141] After calling the marked areas included in the list of interactive identification of the spatial supply capacity level layer and the path structure interference, overlay the two layers, extract the supply level units within the marked areas through spatial intersection, identify the unique code and position coordinates of the units in the spatial layer, read the annual groundwater level monitoring data one by one for the positions where the units are located. The data source can be an automatic water level recorder or a monitoring record form to ensure that the monthly data is complete and the time series is continuous. Sort the monthly water level values of each unit to construct an annual time series table. Set the monthly water level values of unit U308 to be 12.3, 12.1, 11.9, 11.7, 11.8, 12.0, 12.2, 12.5, 12.6, 12.4, 12.2, 12.1 meters respectively. Automatically identify the difference between the maximum and minimum values as the annual variation range. In this example, the variation range is 0.9 meters. Perform the same operation for all supply units, extract the water level variation values and record them in the table. Screen the results, and eliminate the units with no change in the annual water level value or all monthly values being exactly the same, because they do not have the fluctuation characteristics and are not included in the subsequent analysis. Each record includes the unit number, the maximum annual water level value, the minimum water level value, the variation value, and the spatial position, laying a foundation for further identifying the interference influence intensity, and obtaining the groundwater level fluctuation range data set.

[0142] S502: Based on the groundwater level fluctuation range data set, combined with the regional interference level zoning results in the list of interactive identification of the path structure interference, count the change frequency of the fluctuation range values of the supply units within the interference level area, and screen out the supply units whose change frequency exceeds the annual change benchmark within the area to obtain the interference area unit sequence;

[0143] Match and analyze the data with the regional interference levels in the list of path structure interference interactions to identify the interference level areas where each recharge unit is located. Set the interference levels to be divided into low, medium, and high-level areas. In each level area, count the amplitude change values of the water levels of the recharge units in that area, and record the frequency distribution of the amplitude change values in different intervals. For example, the amplitude change occurs 8 times within a range of 0.2 meters, 12 times at 0.4 meters, and 5 times above 0.6 meters. Based on the annual amplitude change frequency distribution of the recharge units in this area, set the annual change reference frequency for this area. For example, set it to occur at a frequency greater than a certain set threshold (such as 1.5 times the average number of occurrences in the area). Compare the amplitude change frequencies of all recharge units with this reference, and screen out the recharge units whose change frequencies exceed this reference. The units are regarded as active fluctuation units and have the characteristics of being enhanced by interference, including information such as unit number, annual amplitude change value, the interference level area to which they belong, the multiple by which the frequency exceeds, and spatial location, etc., to provide a basis for fluctuations in the identification of water supply sensitive areas and form a sequence of interference area units.

[0144] S503: According to the sequence of interference area units, extract the spatial position index values of the corresponding units, combine with the grade zoning at the corresponding positions in the spatial recharge capacity level layer, and through index matching operations, screen out the areas where the distribution ratio exceeds the internal reference level of the corresponding grade area to generate a water supply safety sensitive distribution area layer;

[0145] Extract the corresponding grade categories of the units in the spatial recharge capacity level layer, and call their spatial position index values to form a "unit number - recharge grade" matching table. Through spatial matching algorithms, count the number of interference units included in each recharge grade category and their proportion in the total number of units in the corresponding grade. Set that there are 50 units in the area with a recharge grade of 4, and 18 of them appear in the interference unit sequence, then the interference distribution ratio is 36%. Set the reference level within the grade area, such as the reference is 30%. When the actual ratio exceeds this reference level, it is considered that there is a water supply safety risk sensitive performance in this recharge grade area. Extract the area and mark it as a sensitive area. Reconstruct the sensitive area layer through the spatial index value. Each identification block in this layer contains information such as its grade category, interference unit ratio, exceeding multiple, and spatial boundary, etc. Set the area numbered R42 as the recharge area of grade 4, with an interference ratio of 38%, exceeding the reference of 30%, so it is marked as a sensitive area. This layer will be used as an important spatial reference basis for the safety management of groundwater resources to generate a water supply safety sensitive distribution area layer.

[0146] Please refer to 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 area. After dividing the image into equal-area grids, it calls the permeability values within the grids in the vertical and horizontal directions respectively, constructs unit pairs, compares the permeability differences between adjacent grids in the same direction, connects them in sequence to form a closed boundary, extracts the spatial contour graph, and generates a structural disturbance impact zoning layer;

[0148] The response path extraction module, based on the grid area within the boundary covered by the structural disturbance impact zoning layer, sets up a cross-section measuring point array and collects the time series data of the groundwater level at the corresponding positions, and screens the line segments with consistent change trends to obtain the groundwater response path distribution set;

[0149] The interference mechanism identification module judges the spatial overlap between the coordinates of the corresponding path segments in the groundwater response path distribution set and the boundary line segments of the disturbance limit distribution form value, marks whether the path turns or changes the extension direction at the junction, and obtains the path structure interference interaction identification list;

[0150] The recharge capacity division module calls the path structure interference interaction identification list, extracts the precipitation data, evaporation data and elevation data of the affiliated unit, divides the recharge level according to the combined level of the calculation results, and generates a spatial recharge capacity level layer;

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

[0152] The system of this embodiment can be used to execute Figure 1 the technical solutions of the method embodiment shown, and its implementation principle and technical effects are similar, so they will not be elaborated here.

[0153] As mentioned above, the above are only the specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art can easily think of changes or substitutions within the technical scope disclosed by the present invention, and all should be covered by the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims

1. A method for evaluating the water supply security of a groundwater source area, characterized in that, It includes the following steps: S1: Obtain the permeability distribution map of the area covered by the river-side water source area, divide the image into equal-area grids, analyze the permeability change information in the longitudinal and transverse directions of the unit, construct an edge point set according to the position of continuous jump points, extract the closed boundary shape, and generate a structural disturbance impact zoning layer; S2: Call the grid area in the structural disturbance impact zoning layer, evaluate the consistency between the water level change direction and the fluctuation time per unit time between measurement points, screen the response path line segments according to the coherence degree, and obtain the groundwater response path distribution set; S3: According to the groundwater response path distribution set, perform spatial matching on the overlapping area between the path line segment and the boundary line segment, judge whether the path line crosses the boundary or turns along the boundary, mark the overlapping degree between the path direction mutation position and the fracture boundary, and establish a path structure interference interaction recognition list; S4: According to the path structure interference interaction recognition list, extract the corresponding precipitation data, evaporation data and elevation data of the unit, calculate the precipitation-evaporation ratio, combine it with the elevation difference of adjacent grids to form a relief index value, divide the grades according to the recharge performance, and generate a spatial recharge capacity grade layer.

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 disturbed boundary surface unit, and a continuous jump direction aggregation area. The groundwater response path distribution set includes a cross-sectional path line, a response direction marking point, and a water level coherence record index. The path structure interference interaction recognition list includes a crossing path code, a boundary change direction position identifier, and a recharge interruption area index. The spatial recharge capacity grade layer includes a relief correction index grid, a grade zoning coding map, and a recharge performance distribution layer.

3. The groundwater source water supply safety assessment method according to claim 1, characterized in that The specific steps for obtaining the structural disturbance impact zoning layer are as follows: S101: Obtain the permeability distribution map of the area covered by the river-side water source area, divide it into equal-area grids, regard the grid unit as an independent spatial unit, extract the initial value of the soil permeability within the unit range, match the soil property point coordinates with the corresponding grid numbers, and generate the grid permeability distribution; S102: Call the permeability values of adjacent grid units in the grid permeability distribution, calculate the permeability change difference between the real-time unit and the adjacent unit in the horizontal and vertical directions respectively, and perform a judgment operation on whether there is a direction jump for each pair of direction pairs to generate a direction jump distribution point set; S103: According to the continuous spatial distribution of the jump points in the direction jump distribution point set, identify the jump point groups that meet the continuous quantity requirements, use the edge coordinates of the qualified point groups as edge candidate points, sort the edge candidate point coordinates according to spatial adjacency, extract the set of circumscribed boundary lines forming the closed area, and generate the structural disturbance impact zoning layer.

4. The groundwater source water supply safety assessment method according to claim 3, characterized in that The specific steps for obtaining the groundwater response path distribution set are as follows: S201: Invoke the grid areas in the structural disturbance impact zoning layer, arrange cross-section measurement point arrays in each group of grids at a preset cross-section spacing, collect the continuous-time water level change values of the measurement points, organize the water level time series data within the cross-section according to the measurement point numbers, statistically analyze the change trends of the water level change directions and fluctuation amplitudes between the measurement points per unit time, judge the synchronization of whether the directions and trend changes between the measurement point pairs within the cross-section are consistent, and generate an evaluation result of the consistency of the water level changes of the measurement points; S202: According to the evaluation result of the consistency of the water level changes of the measurement points, screen the measurement point pairs within the cross-section that meet the consistency requirements, construct a connected path with the connection line of the measurement point pairs as the line segment, calculate the corrected value of the connected distance, and integrate the spatial orientation distributions of multiple connected paths to generate a set of groundwater response path distributions.

5. The groundwater source water supply safety assessment method according to claim 4, wherein, The formula for calculating the corrected value of the connected distance is: ； Among them, L ij represents the correction value of the connected distance between measurement points i and j, x i , x j respectively represent the abscissas of measurement points i and j, y i , y j respectively represent the ordinates of measurement points i and j, n represents the total number of observation times, represents the water level difference between measurement points i and j at the k-th moment, represents the average value of the average water levels of measurement points i and j over the entire period, respectively represent the average water levels of measurement points i and j over the entire period, hc represents the centralized reference value of the water level of the cross-section measurement points over the entire period.

6. The groundwater source water supply safety assessment method according to claim 4, characterized in that The specific steps for obtaining the list of path structure interference interaction identifications are as follows: S301: Invoke the spatial coordinates of the path segments in the set of groundwater response path distributions, perform spatial matching with the coordinate set of the boundary segments in the structural disturbance impact zoning layer, judge whether there are intersection points or overlapping segments between the path segments and the boundary segments in terms of coordinate range. If there are intersection points, mark them as crossing; if there are overlapping segments, mark them as turning along the boundary, and generate a dataset of path line and boundary matching types; S302: According to the path segments marked as crossing or turning in the dataset of path line and boundary matching types, extract the path orientation change angles, calculate the orientation mutation degree value according to the angle difference between the front and rear segments of the path, and call the boundary segment coordinates to perform a spatial coincidence degree judgment to generate a path boundary mutation coincidence degree annotation result; S303: Invoke the numbers of the mutation path segments in the path boundary mutation coincidence degree annotation result, mark the path segments with signs of disturbance, summarize the abnormal annotation areas of the partition boundaries where they are located according to the path numbers, extract the corresponding path segments and their locations as the identification objects, calculate the path segment disturbance identification value, and generate a list of path structure interference interaction identifications; The formula for calculating the path segment disturbance identification value is: ； Among them, D a is the path segment perturbation recognition value, N is the total number of path segments, S ao is the matching score between the o-th path segment and the a-th 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 specific steps for obtaining the spatial recharge capacity level layer are as follows: S401: Based on the unit numbers listed in the list of path structure interference interaction identifications, extract the precipitation grid data, evaporation grid data, and elevation grid data of the corresponding unit areas. After unifying the spatial resolutions and projection formats of the three types of data, perform grid clipping according to the unit boundary ranges to obtain the distribution data of the precipitation-evaporation ratio; S402: According to the distribution data of the precipitation-evaporation ratio, call the elevation values of the corresponding pixel positions in the clipped elevation grid data, compare the elevation differences between the real-time pixels and the neighboring pixels, calculate the average height difference of the neighborhood, and perform a product calculation of the average height difference value and the precipitation-evaporation ratio corresponding to the pixel to obtain the distribution value of the undulation index; S403: Invoke the pixel data of the distribution value of the undulation index, perform grade division according to the numerical interval, use the set undulation index grade threshold interval as the classification benchmark, assign different grade values respectively, and generate a grade grid layer according to the spatial positions of the division results to obtain the spatial recharge capacity level layer.

8. The groundwater source area water supply safety assessment method according to claim 1, characterized in that, The method further includes step S5: S5: Call the marked areas in the list of interactive identification of the spatial recharge capacity level layer and path structure interference, overlay the annual fluctuation data of the groundwater level within the overlay area, extract the area units with structural interference and unstable recharge, analyze the response amplitude and frequency characteristics, and obtain the water supply security sensitive distribution area layer; The water supply security sensitive distribution area layer includes response anomaly unit marks, water level fluctuation frequency indexes, and structural recharge overlay area codes.

9. The groundwater source water supply safety assessment method according to claim 8, wherein, The specific steps for obtaining the water supply security sensitive distribution area layer are as follows: S501: Call the marked areas in the list of interactive identification of the spatial recharge capacity level layer and path structure interference. For the recharge units within the marked areas, collect the annual groundwater level monitoring data at the corresponding positions. After arranging the groundwater level values at the same position in monthly sequences, calculate the annual water level change amplitude of the recharge units, and screen out the recharge units without fluctuations to obtain the groundwater level fluctuation amplitude data set; S502: Based on the groundwater level fluctuation amplitude data set, combined with the regional interference level zoning results in the list of interactive identification of path structure interference, count the change frequency of the fluctuation amplitude values of the recharge units within the interference level areas, and screen out the recharge units whose change frequency exceeds the annual change benchmark of the area to obtain the interference area unit sequence; S503: According to the interference area unit sequence, extract the spatial position index values of the corresponding units, combined with the grade zoning at the corresponding positions in the spatial recharge capacity level layer, and through index matching operations, screen out the areas where the distribution ratio exceeds the internal benchmark level of the corresponding grade area to generate the water supply security sensitive distribution area layer.

10. A groundwater source area water supply safety assessment system, characterized in that, The system is used to implement the groundwater source water supply security assessment method described in any one of claims 1-9. The system includes: A structural disturbance identification module, which obtains the permeability distribution map of the area covered by the river-side water source, divides the image into equal-area grids, then calls the permeability values within the grids in the vertical and horizontal directions respectively, constructs unit pairs to compare the permeability differences between adjacent grids in the same direction, connects them in sequence to form a closed boundary and extracts the spatial contour graph, and generates the structural disturbance impact zoning layer; A response path extraction module, based on the grid area within the boundary covered by the structural disturbance impact zoning layer, sets up a cross-section measurement point array and collects the groundwater level time series data at the corresponding positions, and screens out the line segments with consistent change trends to obtain the groundwater response path distribution set; An interference mechanism identification module, according to the coordinates of the corresponding path line segments in the groundwater response path distribution set, makes a spatial overlap judgment with the boundary line segments of the disturbance limit distribution form value, marks whether the path turns or changes the extension direction at the junction, and obtains the list of interactive identification of path structure interference; A recharge capacity division module, which calls the list of interactive identification of path structure interference, extracts the precipitation data, evaporation data and elevation data of the corresponding units, divides the recharge grade according to the combined level of the calculation results, and generates the spatial recharge capacity level layer; Water supply sensitive area identification module, based on the spatial recharge capacity level layer, superimpose the annual groundwater level change data in the area, calculate the combined value of the water level response amplitude and fluctuation frequency corresponding to the interference unit in the recharge grade area, screen the units with the combined value exceeding the distribution mean of the grade area, and obtain the water supply safety sensitive distribution area layer.

Images