Method for comprehensive evaluation of karst fissure zone leakage recharge groundwater suitability

Abstract

This invention relates to the field of groundwater recharge technology, specifically a comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones. The method includes the following steps: acquiring fissure node information to construct paths; calculating density and return rate to identify abrupt changes; extracting flow and head sequences to determine unstable zones; combining lithology and integrity to locate abrupt change nodes; analyzing fissures and impermeability to screen interfering paths; integrating overlapping segments to generate weakened patches and mapping them to an evaluation layer to form an influence map. In this invention, by constructing fissure node paths and calculating return rate and density changes, structural abrupt change segments are identified; by superimposing hydraulic response sequences to calculate fluctuation intensity, unstable recharge zones are determined; lithological labels and integrity indices are extracted; microscale structural variations are identified; abrupt change nodes are located; and by combining fissure density and impermeability level, interfering paths are screened. Overlapping paths are integrated into weakened patches and mapped to an evaluation layer, thus constructing a linkage mechanism for identification, response, interference judgment, and evaluation constraints.

Description

Technical Field

[0001] This invention relates to the field of groundwater recharge technology, and in particular to a comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones. Background Technology

[0002] Groundwater recharge technology is an important branch of hydrogeology, primarily studying the processes, mechanisms, and patterns of water recharge from natural or artificial sources to underground aquifers. Core aspects of this field include recharge type identification, recharge path tracing, recharge volume calculation, and recharge suitability analysis. It typically employs a comprehensive approach combining field surveys, hydrogeological mapping, geophysical exploration, hydrochemical analysis, and numerical simulation to achieve effective understanding and scientific assessment of groundwater recharge processes. Traditional karst fissure zone seepage recharge suitability evaluation refers to the feasibility analysis and location selection of surface water recharge processes from fissure systems in karst-developed areas. The key technical challenge is how to scientifically determine the suitability of seepage recharge areas in environments with complex karst fissure structures and variable hydrological conditions. Traditional methods typically employ karst fissure development observation, hydrogeological parameter surveys, groundwater level dynamic monitoring, and seepage path analysis to complete the evaluation.

[0003] Traditional technologies rely heavily on observations of fracture development, hydrogeological parameter surveys, and groundwater level monitoring. While possessing some macroscopic identification capabilities, they struggle to achieve high-precision local identification in complex karst regions, particularly in areas with densely interwoven fracture systems or diverse lithologies. The lack of methods to detect abrupt structural changes at the node level can lead to insufficient identification of important recharge pathways. Existing technologies typically lack dynamic combination analysis of fracture paths, making it difficult to extract anomalies from changes in path network characteristics. Significant gaps exist in spatial matching and hydraulic fluctuation coupling analysis, resulting in delayed and incomplete identification of unstable areas. For example, in areas with frequent lithological changes or dramatic hydraulic head fluctuations, relying solely on static hydrogeological parameters often fails to reveal the dynamic response patterns of the groundwater system, easily leading to misjudgments in recharge suitability assessments. The absence of a closed-loop mechanism linking structural identification, dynamic response, and limiting factors makes it difficult to comprehensively reflect the actual interference of geological conditions on recharge pathways. These problems substantially constrain the scientific selection of recharge engineering sites, increasing the risk of resource waste during implementation. Summary of the Invention

[0004] To address the technical problems existing in the prior art, embodiments of the present invention provide a comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones, comprising the following steps:

[0005] To achieve the above objectives, the present invention employs the following technical solution: a comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones, comprising the following steps:

[0006] S1: Obtain spatial information of the cross nodes of the crack, construct path combinations, calculate the turnback rate and node density changes, identify the location of structural abrupt changes, and generate a set of intersecting abnormal path segments.

[0007] S2: Based on the path coordinates in the set of intersecting abnormal path segments, extract the flow change sequence and head fluctuation sequence of the corresponding monitoring point, calculate the fluctuation intensity, and combine it with the classification standard to generate the range of unstable supply section.

[0008] S3: Within the unstable supply zone, extract the lithological tag sequence and integrity index along the path, calculate the difference between the lithological transformation frequency and the integrity change rate, locate nodes that exceed the discrimination threshold, and generate a set of lithological structure mutation nodes.

[0009] S4: Based on the spatial location of the set of abrupt change nodes in the lithological structure, extract the fracture line density and impermeability grade on both sides of the node, calculate the combination variation characteristics, screen out path segments that meet the interference identification criteria, and generate a set of structural interference paths.

[0010] S5: Based on the set of structural interference paths, match overlapping areas in the unstable supply segment range, extract and integrate overlapping path segments, generate weakened area patches, map them onto the suitability evaluation layer as a constraint factor, and form suitability weakening influence area patches.

[0011] As a further aspect of the present invention, the specific steps of S1 are as follows:

[0012] S101: Based on the three-dimensional coordinate data of the cross nodes of the crack, extract the spatial connection relationship between the nodes, filter the node pairs that meet the continuous connection condition, construct the path set, record the connection number and index information of the nodes in the path, and obtain the node path combination information.

[0013] S102: Call the path number and node sequence in the node path combination information, extract the direction change sequence formed by the nodes in the path, calculate the ratio of the direction change amount to the number of path segments in each path, and obtain the path turnaround rate sequence value;

[0014] S103: Based on the path turnaround rate sequence value and the node coordinates in the node path combination information, extract the density change amplitude of adjacent node segments, filter the location segments that simultaneously meet the conditions of turnaround rate change and density change, and obtain a set of intersecting abnormal path segments.

[0015] As a further aspect of the present invention, the specific steps of S2 are as follows:

[0016] S201: Based on the path coordinates in the set of intersecting abnormal path segments, extract the monitoring point numbers within the coverage area of ​​the path segments, obtain the flow change sequence and head fluctuation sequence for the corresponding time period, and classify them according to the path segment numbers to obtain the path segment flow and head sequence set.

[0017] S202: Based on the flow change value and head fluctuation value of the path segment in the sequence of flow and head, extract the change amplitude and fluctuation trend of the sequence, calculate the fluctuation statistics value corresponding to the path segment, and obtain the fluctuation intensity value of the path segment.

[0018] S203: Call the path segment fluctuation intensity value and path segment number, determine the fluctuation level of the path segment according to the intensity threshold of the fluctuation classification standard, filter the path segment numbers whose fluctuation level is in the unstable range, and obtain the range of the unstable supply segment.

[0019] As a further aspect of the present invention, the specific steps of S3 are as follows:

[0020] S301: Based on the path coordinates within the supply unstable section, extract the lithological label sequence and corresponding integrity index along the path, establish the correspondence between the path location and lithological and integrity data, and obtain the lithological and integrity sequence set along the path;

[0021] S302: Based on the lithology labels and integrity indices in the lithology and integrity sequence set along the path, the lithology transformation frequency and integrity change rate are statistically analyzed, and the difference between the two is obtained to obtain the lithology transformation frequency and integrity change rate difference sequence.

[0022] S303: Based on the difference sequence between the lithological transformation frequency and the integrity change rate, call the discontinuity discrimination threshold to determine whether the difference exceeds the limit, mark the path node positions that meet the conditions, and obtain the set of lithological structure mutation nodes.

[0023] As a further aspect of the present invention, the specific steps of S4 are as follows:

[0024] S401: Based on the spatial location provided by the set of abrupt change nodes in the lithological structure, extract the fracture line density and impermeability grade of the path segments on both sides of the node, establish the pairing relationship between node location and associated attributes, and obtain the sequence of fracture and impermeability attributes of the path node.

[0025] S402: Based on the crack linear density and impermeability grade in the crack and impermeability attribute sequence of the path node, calculate the combination change characteristics, determine the attribute change of the path segments on both sides of the node, and obtain the crack impermeability combination change characteristic sequence.

[0026] S403: Based on the change feature values ​​of the path segments in the fracture anti-seepage combination change feature sequence, call the interference identification standard threshold, filter the path number segments that meet the conditions, and obtain the set of structural interference paths.

[0027] As a further aspect of the present invention, the interference identification standard threshold is determined by preset criteria for the degree of change in crack development and the level of change in impermeability. When the degree of crack development on both sides of the path segment of the path node is significantly enhanced and the level of impermeability changes significantly, it is identified as a structural interference path segment.

[0028] As a further aspect of the present invention, the specific steps of S5 are as follows:

[0029] S501: Based on the spatial boundary of the path in the set of structural interference paths, match the overlapping positions in the boundary of the unstable supply section, extract the overlapping path number and coordinate section, and obtain the overlapping path segment number set;

[0030] S502: Based on the coordinate information corresponding to the overlapping path segment number set, integrate the path boundary and the unstable segment boundary, unify the spatial contour of the overlapping segment, and generate the spatial weakening patch contour value.

[0031] S503: Based on the position and shape parameters of the spatial weakening patch outline values, match the grid distribution in the suitability evaluation layer, map the constraint parameters, and obtain the suitability weakening influence area patch.

[0032] As a further aspect of the present invention, the path combination is a set of paths constructed based on the connectivity between the cross nodes of the crack;

[0033] The turnaround rate is the ratio of the total propagation distance in the path to the straight-line distance between the starting and ending points;

[0034] The node density change is the gradient of the number of fracture nodes per unit distance.

[0035] The intersecting abnormal path segments are crack path regions where the change in turnback rate and node density is abrupt, and these path segments are identified as structural disturbance sensitive areas.

[0036] The flow rate change sequence is a set of flow rate data collected by monitoring points deployed along the fracture path over a continuous time period.

[0037] The hydraulic head fluctuation sequence is a record sequence of elevation values ​​formed by the change of groundwater level over time. The hydraulic head fluctuation reflects the disturbance response capability of the groundwater system.

[0038] The fluctuation intensity is a fluctuation index obtained by combining the amplitude of flow rate change and the amplitude of head fluctuation.

[0039] The grading standard is a classification rule for the wave intensity value based on statistical parameters such as the coefficient of variation, and is determined based on regional hydrological statistical analysis.

[0040] The range of the unstable supply section is the set of spatial boundaries of the path segments corresponding to the fluctuation intensity values ​​exceeding the stability limit.

[0041] As a further aspect of the present invention, the lithological tag sequence is an ordered set composed of lithological codes of rock strata extracted sequentially along the path direction;

[0042] The integrity index is a parameter used to measure the continuity and stability of rock strata structure, derived from the calculation results of rock strata continuity and tectonic disturbance level data;

[0043] The lithology transformation frequency is the number of times the lithology label changes per unit path length;

[0044] The discrimination threshold is a frequency difference limit used to define whether lithological changes have reached the degree of structural abruptness, and is set based on statistical experience values ​​or simulated judgment values.

[0045] The crack line density is the total length of crack line segments within a unit area;

[0046] The impermeability grade is a classification grade based on the formation permeability coefficient and fracture distribution, which determines the ability to resist water permeability.

[0047] The combined variation characteristic is the joint variation trend of fracture linear density and impermeability grade at the location of structural abrupt change;

[0048] The interference identification standard is a multi-parameter discrimination system used to determine whether a path segment is severely disturbed by structural factors, and is set by fracture characteristics, hydraulic properties, and lithological grade.

[0049] The set of structural interference paths is a spatial set of path segments that satisfy the interference identification conditions;

[0050] The overlapping region is the set of locations where two path sets intersect and overlap in spatial coordinates.

[0051] The weakened region patch is a graphic unit formed by the intersection of the interference path and the unstable supply path. The patch is used to limit the suitability evaluation results and is the representation of the unfavorable factor layer.

[0052] The suitability evaluation layer is a spatial evaluation unit layer constructed based on the regional groundwater recharge conditions;

[0053] The suitability reduction effect area patch is an effect area unit formed by overlaying restrictive patches onto the suitability layer.

[0054] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0055] In this invention, by constructing spatial connections and combining paths at fracture intersection nodes, and by integrating the return rate and the magnitude of node density changes, the invention achieves accurate identification of structurally abrupt change zones. It overlays hydraulic response sequences and calculates wave intensity values ​​to effectively identify unstable recharge zones. It extracts lithological labels and integrity indices, and uses difference indicators to identify microscale structural variations, further locating lithological structural abrupt change nodes. It calculates the combined characteristics of fracture line density and impermeability grade, accurately screening out interference paths. Overlapping path segments are integrated into weakened patches and mapped onto the evaluation layer. This constructs a linkage mechanism from identification, response analysis, interference judgment to evaluation constraints, enhancing the scientific rigor and adaptability of recharge feasibility assessment in complex geological environments. Attached Figure Description

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

[0057] Figure 1 This is a schematic diagram of the steps of the present invention;

[0058] Figure 2 This is a detailed schematic diagram of S1 of the present invention;

[0059] Figure 3 This is a detailed schematic diagram of S2 of the present invention;

[0060] Figure 4 This is a detailed schematic diagram of S3 of the present invention;

[0061] Figure 5 This is a detailed schematic diagram of S4 of the present invention;

[0062] Figure 6 This is a detailed schematic diagram of S5 of the present invention. Detailed Implementation

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

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

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

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

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

[0068] Please see Figure 1 This invention provides a comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones, comprising the following steps:

[0069] S1: Obtain spatial information of the cross nodes of the crack, construct path combinations, calculate the path turnaround rate and the change range of node density, identify the location segments where structural abrupt changes occur, and generate a set of intersecting abnormal path segments;

[0070] Path combinations are sets of paths constructed based on the connectivity between fracture intersection nodes, used to simulate the spatial connectivity paths through which water flows.

[0071] The path turnaround rate is the ratio of the total propagation distance in a path to the straight-line distance between the start and end points, and is used to measure the tortuosity of the path;

[0072] The magnitude of the node density variation is the gradient of the number of fracture nodes within a unit distance, used to characterize the uniformity and abrupt changes in the density of fracture distribution.

[0073] Intersecting abnormal path segments are crack path regions where the change in turnback rate and node density is abrupt, and these path segments are identified as structural disturbance sensitive areas.

[0074] S2: Based on the path coordinates in the set of intersecting abnormal path segments, extract the flow change sequence and head fluctuation sequence of the corresponding monitoring points, calculate the fluctuation intensity value of the path segment, and make a judgment in combination with the fluctuation classification standard to generate the range of the unstable supply section.

[0075] The flow change sequence is a set of flow data collected from monitoring points deployed along the fracture path over a continuous time period, used to assess the dynamic response characteristics during the groundwater recharge process.

[0076] Hydraulic head fluctuation sequence is a record sequence of elevation values ​​formed by the change of groundwater level over time. Hydraulic head fluctuation reflects the disturbance response capability of the groundwater system.

[0077] The fluctuation intensity value is a fluctuation index obtained by combining the amplitude of flow change and the amplitude of head fluctuation, and is used to express the degree of dynamic changes in hydrology of the route segment.

[0078] The fluctuation classification standard is a classification rule for the fluctuation intensity value based on statistical parameters such as the coefficient of variation. It is used to classify the replenishment path according to the level of hydrological fluctuation and is determined based on regional hydrological statistical analysis.

[0079] The range of unstable supply sections is the set of spatial boundaries of the path segments where the fluctuation intensity value exceeds the stability limit. It is used to identify path areas where the supply capacity decreases or changes drastically.

[0080] S3: Based on the path range within the unstable supply section, extract the lithological tag sequence and integrity index along the path, calculate the difference between the lithological transformation frequency and the integrity change rate, locate nodes that exceed the discontinuity discrimination threshold, and generate a set of lithological structure abrupt change nodes.

[0081] The lithological tag sequence is an ordered set of lithological codes extracted sequentially along the path direction, used to reflect the changes in lithological types traversed by the water flow path;

[0082] The integrity index is a parameter used to measure the continuity and stability of rock strata structure. It is derived from the calculation results of rock strata continuity and tectonic disturbance level data and is used to identify potential structural abrupt change locations.

[0083] Lithology transformation frequency is the number of times the lithology label changes within a unit path length, used to characterize the degree of change in lithological structure;

[0084] The discontinuity discrimination threshold is a frequency difference limit used to define whether lithological changes have reached the degree of structural abruptness. It is set based on statistical experience values ​​or simulated judgment values.

[0085] S4: Based on the spatial location provided by the set of abrupt change nodes in lithological structure, extract the fracture line density and impermeability grade on both sides of the node, calculate the combination variation characteristics, screen the path segments that meet the interference identification criteria, and generate a set of structural interference paths.

[0086] Crack linear density is the total length of crack segments per unit area;

[0087] The impermeability grade is a classification grade based on the formation's permeability coefficient and fracture distribution, which determines its ability to resist water seepage.

[0088] The combined variation characteristic is the joint variation trend of fracture linear density and impermeability grade at the location of structural abrupt change. This characteristic is used to determine whether there is a complex interference.

[0089] The interference identification standard is a multi-parameter discrimination system used to determine whether a path segment is severely disturbed by structural factors. It is set by fracture characteristics, hydraulic properties, and lithological grade.

[0090] The structural interference path set is a spatial set of path segments that meet the interference identification conditions, and is used as key input for restricted area identification and suitability reduction determination;

[0091] S5: Based on the spatial boundary of the path in the set of structural interference paths, match the overlapping areas in the range of the unstable supply section, extract the overlapping path segments and integrate them into the weakened area patches, map them to the suitability evaluation layer as a constraint factor, and generate the suitability weakening influence area patches.

[0092] Spatial boundaries are the outer boundary contours of path segments constructed in geographic coordinate systems or model coordinate systems. These boundaries are used for spatial intersection analysis.

[0093] Overlapping regions are the sets of locations where two sets of paths intersect in spatial coordinates, used to identify interactive regions between paths affected by interference from the same source;

[0094] The weakened region patch is a graphic unit formed by the intersection of interference path and unstable supply path. The patch is used to limit the suitability evaluation results and is the representation of the unfavorable factor layer.

[0095] The suitability assessment layer is a spatial assessment unit layer constructed based on regional groundwater recharge conditions, used to represent the distribution of suitability levels;

[0096] The suitability reduction influence area patch is an influence area unit formed by overlaying restrictive patches onto the suitability layer, used to apply spatial corrections in the evaluation model.

[0097] The set of intersecting abnormal path segments includes path turnback rate, node density variation range, and structural abrupt change location segments. The range of unstable recharge sections includes path coordinates, flow change sequence, head fluctuation sequence, and fluctuation intensity value. The set of nodes with lithological structural abrupt changes includes lithological label sequence, integrity index, and the difference between lithological transformation frequency and integrity change rate. The set of structural interference paths includes fracture line density, permeability grade, and combination variation characteristics. The suitability weakening influence area patches include the spatial boundary of structural interference path segments, overlapping areas of unstable recharge sections, weakened area patches, and constraint factors in the suitability evaluation layer.

[0098] Please see Figure 2 The specific steps of S1 are as follows:

[0099] S101: Based on the three-dimensional coordinate data of the cross nodes of the crack, extract the spatial connection relationship between the nodes, filter the node pairs that meet the continuous connection condition, construct the path set, record the connection number and index information of the nodes in the path, and obtain the node path combination information.

[0100] The 3D coordinate data of fracture intersection nodes are typically obtained through geological modeling or laser scanning. After acquiring these spatial points, the spatial distance between node pairs needs to be constructed one by one. The calculation method is to obtain the straight-line distance between any two nodes based on their coordinate differences. If the distance is less than the set connection limit, the two nodes are considered to be related. The connection limit can be determined based on the spatial scale of the fractures. For example, if the average fracture spacing is 1.5m, this value can be used as a reference standard for determining the connection. After establishing all node connections, a depth-first search or breadth-first search is used to determine whether the nodes form a continuous path. The condition for a continuous connection is that each intermediate node in the path has a bidirectional connection, and the angle between adjacent connecting lines cannot exceed a certain angle limit, such as 30 degrees, to ensure the stability of the path direction. The angle is calculated based on the vector direction difference formed by three points. The determination method is that if the angle is less than the set angle limit, it is considered a continuous connection. Through the above screening, all node pairs that meet the continuity condition are obtained, and the sequence number and path number of each pair of nodes in the path are recorded, forming a complete path combination record. For example, if a path consists of 4 nodes in sequence, the corresponding path number is set as T001, and the node sequence numbers are 1 to 4. Finally, a path set containing path numbers and corresponding node sequence information is constructed as the basis for subsequent path analysis.

[0101] S102: Call the path number and node sequence in the node path combination information, extract the direction change sequence formed by the nodes in the path, calculate the ratio of the direction change amount to the number of path segments in each path, and obtain the path turnaround rate sequence value;

[0102] The system reads path numbers and node sequences one by one from the existing path set. For each path, a direction vector sequence is constructed based on the node sequence, with the direction vector formed by the coordinate difference between two consecutive nodes. The angle between any two consecutive vector segments reflects the change in path direction; a smaller angle indicates the path tends towards a straight line, while a larger angle indicates a significant bend. All angle values ​​for each path are calculated and summed to obtain the total change in path direction. Then, the average is taken based on the number of path segments to obtain the path's direction change ratio. This ratio is the path turnaround rate. If a path contains 3 connecting segments and the total change in direction is 60 degrees, the turnaround rate is 20 degrees. The turnaround rate of each path is recorded as a sequence to identify cases of large direction jumps within the path. The critical value for the turnaround rate is set according to application requirements; for example, setting it to 15 degrees marks all paths exceeding this value as having a high turnaround rate. If there are error offsets in the nodes within the path, the direction vector will also be affected. Therefore, the accuracy of the path direction can be improved before processing through averaging filtering or error correction methods. Once completed, output the turnaround rate information for all paths and associate it with the path number for subsequent judgment.

[0103] S103: Based on the path turnaround rate sequence value and the node coordinates in the node path combination information, extract the density change amplitude of adjacent node segments, filter the location segments that simultaneously meet the conditions of turnaround rate change and density change, and obtain the set of intersecting abnormal path segments.

[0104] Paths exceeding a critical turnaround value are selected from the existing turnaround rate sequence. If the critical value is 15 degrees, all paths with turnaround rates higher than this value will be marked for subsequent density analysis. Next, all adjacent node segments within the path are extracted, and the ratio of the number of nodes in each segment to the segment length is calculated to obtain the node density per unit length. For example, if a segment is 5 meters long and contains 8 nodes, its density is 1.6 nodes per meter. The density values ​​of all segments are calculated sequentially according to the path nodes to obtain the density variation between consecutive segments. If a segment has a density of 1.6 and the preceding segment has a density of 1.0, the variation is 0.6. A criterion for density variation is set, such as a value of 0.3. Any segment with a density variation greater than this value and a turnaround rate greater than 15 degrees is defined as an abnormal path segment. Such path segments typically exhibit frequent changes in direction and uneven node distribution, which may represent areas of structural disturbance or overlapping fractures. Record all path segments that meet the above conditions as an anomaly set. The record includes information such as the path number, the start and end node numbers and coordinates of the segment, the turnaround rate value and the density change range, etc., to form a list of intersecting anomaly segments for subsequent analysis or annotation.

[0105] Please see Figure 3 The specific steps of S2 are as follows:

[0106] S201: Based on the path coordinates in the set of intersecting abnormal path segments, extract the monitoring point numbers within the coverage area of ​​the path segments, obtain the flow change sequence and head fluctuation sequence for the corresponding time period, and classify them according to the path segment number to obtain the path segment flow and head sequence set;

[0107] Based on the path coordinates in the set of intersecting abnormal path segments, the starting and ending coordinates of each path segment are first used as the basis for spatial query to construct a spatial envelope boundary encompassing the entire length of the path segment. The boundary range is set as the minimum envelope cube of the segment in the three-dimensional coordinate system. Then, the monitoring point coordinate database of the groundwater monitoring system is called to compare whether the position of each monitoring point is within the spatial envelope of the corresponding path segment. If the condition is met, the monitoring point number is associated with the path segment. Subsequently, the flow change data sequence and head fluctuation data sequence within a specified time period are retrieved based on the monitoring point number. This time period is generally set to... For 24 to 48 hours after path identification, continuously recorded data is obtained at 1-hour intervals. For example, for monitoring point M001, the flow data recorded in the 0-24 hour interval may be 23 values. The head data is also recorded synchronously as a sequence of the same time length. All data are classified and collected according to the path segment number. Each path segment may be associated with multiple monitoring points. Each monitoring point corresponds to an independent flow and head data pair. This process is executed line by line until all path segments have completed the pairing of the corresponding monitoring points with their flow and head data, forming a structured data set for subsequent fluctuation intensity analysis.

[0108] S202: Based on the flow change value and head fluctuation value of the path segment in the sequence of flow and head, extract the change amplitude and fluctuation trend of the sequence, calculate the fluctuation statistics value corresponding to the path segment, and obtain the fluctuation intensity value of the path segment.

[0109] Based on the path segment and its associated flow and head data set, the monitoring point data associated with each path segment are processed one by one. First, the maximum and minimum values ​​of the flow sequence for each monitoring point are identified, and the difference between them is calculated and compared with the mean to obtain the percentage change in flow. For example, if the maximum flow of a monitoring point is 3.5, the minimum is 2.7, and the average is 3.1 within a day, then its change is approximately 26%. The head data is processed in a similar manner to obtain its maximum fluctuation value, minimum fluctuation value, and average fluctuation level within that time range. Simultaneously, the fluctuation trend of the data sequence is analyzed, i.e., whether the head value is continuously rising, falling, or oscillating at high frequency, by comparing adjacent... The sign of the data difference can determine the continuity of the direction of change. Then, the standard deviation, range, and maximum instantaneous change of the change sequence are calculated to obtain multiple statistical characteristic values ​​of flow and head. These multiple indicators are weighted according to a set ratio, such as setting the standard deviation weight to 40%, the range weight to 30%, and the maximum change weight to 30%. The multiple statistical values ​​are summed proportionally to obtain the comprehensive fluctuation intensity value. For example, if the fluctuation intensity calculated for a certain monitoring point is 0.52, and there are multiple monitoring points under the path segment, the average value is taken as the overall fluctuation intensity value of the path segment. All path segments are executed in sequence to form a data sequence containing the correspondence between the path segment number and its fluctuation intensity, which is used for the next step of stability classification judgment.

[0110] S203: Call the path segment fluctuation intensity value and path segment number, determine the fluctuation level of the path segment according to the intensity threshold of the fluctuation classification standard, filter the path segment numbers whose fluctuation level is in the unstable range, and obtain the range of the supply unstable segment.

[0111] Based on the generated set of path segment fluctuation intensity values ​​and path segment numbers, the fluctuation intensity is classified into four levels according to the established grading standard: Level 1 is the stable zone, with fluctuation intensity less than or equal to 0.2; Level 2 is the mild wave zone, with fluctuation intensity between 0.2 and 0.4; Level 3 is the medium wave zone, with fluctuation intensity between 0.4 and 0.6; and Level 4 is the strong wave zone, with fluctuation intensity higher than 0.6. For example, if the fluctuation intensity value of path segment T001 is 0.54, it is determined to belong to the medium wave zone. Similarly, if the fluctuation intensity of path segment T002 is 0.67, it is classified into the strong wave zone. All path segments must be classified according to this standard. The grading process is completed by judging the numerical range of each fluctuation intensity value. Then, all path segment numbers with fluctuation levels in the medium wave zone and strong wave zone are selected. These path segment numbers, such as T001, T002, and T005, will form a set of unstable segment ranges, which will be output as the final identification result, forming a data list of the range for further processing.

[0112] Please see Figure 4 The specific steps of S3 are as follows:

[0113] S301: Based on the path coordinates within the unstable supply zone, extract the lithological label sequence and corresponding integrity index along the path, establish the correspondence between the path location and lithological and integrity data, and obtain the lithological and integrity sequence set along the path;

[0114] Based on the path coordinates within the unstable supply zone, all path segments are first numbered and managed. A continuous sequence of path points is constructed according to the start and end coordinates of each path segment. Spatial sampling is performed at equal intervals, such as collecting path node coordinates every 0.5m, forming a complete set of path spatial points. Then, lithological distribution information and integrity index layers from the regional geological database are used to spatially match the coordinates of each path point with the 3D geological model, extracting the corresponding lithological labels, such as sandstone, siltstone, mudstone, etc., and obtaining the integrity index for that location. The integrity index value ranges from 0 to... 1. Generally, a value closer to 1 indicates a more complete structure. For example, if the lithology corresponding to a path point is sandstone and the integrity index is 0.83, it is recorded as the attribute data of that node. Then, the same operation is performed on all nodes of the path to form the lithology sequence and integrity sequence of that path segment. For example, if the lithology label is [sandstone, sandstone, mudstone, mudstone, limestone], the corresponding integrity is [0.83, 0.80, 0.58, 0.55, 0.61]. The lithology and integrity data of all nodes formed by multiple path segments will be collected to construct a data set corresponding to the path segment location, lithology, and integrity, which prepares the foundation for the next step of data analysis.

[0115] S302: Based on the lithology labels and integrity indices in the lithology and integrity sequence set along the route, the lithology transformation frequency and integrity change rate are statistically analyzed, and the difference between the two is obtained to obtain the lithology transformation frequency and integrity change rate difference sequence.

[0116] Based on the lithology and integrity data set formed along the route, node-level analysis is performed on each route segment. First, the frequency of lithology label changes is statistically analyzed. Starting from the first node of the route, the lithology category of each adjacent node is compared to see if it changes. For example, if the lithology sequence of a route is [sandstone, sandstone, mudstone, mudstone, limestone], and sandstone changes to mudstone and mudstone changes to limestone (two lithology transformations), and there are four adjacent node pairs in the route, then the lithology change frequency is 2 to 4, which translates to 50%. Next, the integrity index difference between adjacent nodes in the integrity sequence is extracted, and its absolute value is taken. All differences are then merged. Divide by the number of nodes minus 1. For example, if the integrity sequence is [0.83, 0.80, 0.58, 0.55, 0.61], its variation range is 0.03, 0.22, 0.03, 0.06, totaling 0.34, with an average variation rate of 8.5%. Finally, the difference between the lithology variation frequency and the integrity variation rate is calculated, that is, the lithology frequency minus the integrity variation rate. In the above example, 50% minus 8.5% equals 41.5%. This difference reflects the difference between lithology category mutation and structural response. Each path segment is processed in the above manner to obtain the difference sequence for subsequent mutation point judgment.

[0117] S303: Based on the difference sequence between lithological transformation frequency and integrity change rate, call the discontinuity discrimination threshold to determine whether the difference exceeds the limit, mark the path node positions that meet the conditions, and obtain the set of lithological structure abrupt change nodes;

[0118] Based on the difference sequence calculation results, a discontinuity discrimination threshold is set to identify lithological abrupt change nodes. It is recommended that this threshold be set based on historical geological structure data and monitoring analysis results, generally chosen as 20%. This means that when the difference corresponding to a node is higher than 20%, it is considered that the location exhibits characteristics of drastic lithological transformation but insufficient structural integrity response, belonging to a potential lithological structural abrupt change point. The judgment process involves scanning the difference sequence node by node along each path segment, identifying all node numbers whose differences exceed the set threshold, recording the path segment and node number, and extracting their three-dimensional coordinates to form a node set. For example, if the path segment number is P07 and the difference of the 4th node is 37.2%, exceeding the threshold, it is marked as P07-4, with corresponding coordinates (156.4, 32.7, 9.8). All nodes meeting the criteria are then assigned to the abrupt change node set. The output structure includes the path segment number, node number, difference percentage, and spatial coordinates, ultimately constructing a lithological structural abrupt change node dataset for subsequent geological section assessment and regional response control parameter input.

[0119] Please see Figure 5 The specific steps of S4 are as follows:

[0120] S401: Based on the spatial location provided by the set of abrupt change nodes in lithological structure, extract the fracture line density and impermeability grade of the path segments on both sides of the node, establish the pairing relationship between node location and associated attributes, and obtain the sequence of fracture and impermeability attributes of path nodes.

[0121] Based on the spatial location information provided by the lithological structure abrupt change node set, it is necessary to extract the adjacent path segments before and after each abrupt change node. Using the node as the center point, extend approximately 2 meters forward and backward, dividing the path into a segment before and after the node. Then, using the node number as the retrieval condition, extract the spatial coordinate range of the preceding and following segments from the path database. Next, extract the geological attributes of this range, specifically including extracting the fracture line distribution map and impermeability parameter map of the corresponding path segment from the geotechnical engineering information system. The fracture line density is calculated based on the total fracture length per unit path length. For example, if a total fracture length of 1.6 meters is identified within a 2-meter path segment, the density is 0.8 m / m. The impermeability grade is determined according to national or international standards. Industry-related standards classify the nodes into levels I to V, with levels 1 to 5 corresponding to the strongest to the weakest impermeability. For example, if the impermeability level of the preceding section is level II, it is recorded as 2, and if the subsequent section is level IV, it is recorded as 4. In this way, values ​​are assigned to the preceding and following path segments of each node. Each node is combined with its corresponding two sets of attributes to form a bidirectional pairing relationship. For example, the attribute pairing for node number N15 is fracture density 0.8 and impermeability level 2, and the other side is fracture density 1.3 and impermeability level 4. After extracting the preceding and following information of all abrupt change nodes in sequence, an attribute sequence containing node number, preceding fracture density, preceding impermeability level, subsequent fracture density, and subsequent impermeability level is finally constructed for subsequent change analysis and processing.

[0122] S402: Based on the crack linear density and impermeability grade in the crack and impermeability attribute sequence of the path node, calculate the combination change characteristics, determine the attribute change of the path segments on both sides of the node, and obtain the crack impermeability combination change characteristic sequence.

[0123] Based on the established sequence of fracture and impermeability attributes at each path node, a change characteristic analysis is performed on each node. First, the difference in fracture linear density between the preceding and following segments is calculated, using absolute values ​​to measure the degree of change. For example, if the linear density is 0.8 m / m in the preceding segment and 1.3 m / m in the following segment, the difference is 0.5 m / m. Changes in impermeability grade are measured by the absolute value of the grade difference. For example, if the grade is 2 in the preceding segment and 4 in the following segment, the difference is grade 2. Based on these two indicators, a fracture-impermeability change pair is formed for each node. Then, all nodes are processed uniformly to form a sequence set of fracture change values ​​and impermeability grade change values. For ease of comparison, data with different dimensions can be processed... Normalization is performed. For example, if the maximum difference in fracture density in the sample is 2 m / m, the normalized value is the actual difference divided by 2. If the maximum difference in impermeability grade is 4, the normalized result is the difference in impermeability grade divided by 4. Taking the above example, a change in fracture density of 0.5 corresponds to a normalized value of 0.25, and a change in impermeability grade of 2 corresponds to a normalized value of 0.5. These two values ​​are used to form a feature pair to represent the change characteristics of the node. This method is then extended to all nodes to form a complete sequence of changes in fracture and impermeability combination features. Through this sequence, the change trend of the geological structure on both sides of each node can be visualized and numerically expressed.

[0124] S403: Based on the change feature value of the path segment in the fracture anti-seepage combination change feature sequence, call the interference identification standard threshold, filter the path number segments that meet the conditions, and obtain the set of structural interference paths;

[0125] Based on the obtained characteristic sequence of fracture impermeability combination, a set of standard thresholds for identifying structural interference paths needs to be set. The threshold setting should be based on a comprehensive judgment of engineering experience, historical project data, and actual geological conditions. For example, the identification threshold for the comprehensive change characteristic value can be set to 0.6. The combined characteristic value can be obtained by weighted averaging, assigning a weight of 0.4 to the normalized value of fracture density change and a weight of 0.6 to the normalized value of impermeability grade change. The combined value is the sum of the two values ​​multiplied by their respective weights. For example, if the normalized value of fracture density change at a certain node is 0.25 and the normalized value of impermeability grade is 0.5, the calculated combined change characteristic value is 0.4 multiplied by 0.25 plus 0. Multiplying 6 by 0.5 yields 0.4. If the result is below the threshold, the corresponding path segment is not included in the interference path. If the combined feature value of another node is 0.73, it is greater than 0.6, and its path segment should be marked as a structural interference path. For each node, its combined change feature value is compared with the standard threshold. If any exceeds the threshold, its corresponding path segment number is extracted and recorded. For example, if the path segment numbers corresponding to node N21 are L45 and L46, and the interference condition is met, L45 and L46 should be included in the structural interference path set. Finally, a list of all path numbers that meet the interference identification criteria is formed, which can serve as an important spatial positioning basis for subsequent structural intervention and engineering treatment.

[0126] Please see Figure 6The specific steps of S5 are as follows:

[0127] S501: Based on the spatial boundary of the path in the set of structural interference paths, match the overlapping positions in the boundary of the unstable supply section, extract the overlapping path number and coordinate section, and obtain the overlapping path segment number set;

[0128] Based on the spatial boundaries of the paths in the structural interference path set, the start and end coordinates of each path need to be extracted to form a boundary segment dataset. The boundary information of the unstable section needs to be read and processed synchronously. The two types of boundary data are transformed with the same spatial precision, and coordinate interpolation is performed at 1-meter intervals to ensure spatial comparability between the path boundaries and the unstable section boundaries. Then, each path segment coordinate sequence is matched point-by-point with the unstable section coordinate sequence using spatial comparison. A Euclidean spatial error threshold of 0.5 meters is set for each path segment. The number of coordinate points that coincide with the boundary of the unstable region is counted, and the proportion of this number to the total number of points in the path segment is defined as the coincidence rate. The coincidence rate identification threshold is set to 30%. If a path segment has a total of 40 coordinate points sampled, and 14 of these points are less than 0.5m away from the unstable region points, the calculated coincidence rate is 35%, which is higher than the set threshold. Therefore, the path segment is identified as a coincident segment, and the path segment number and start and end coordinates are recorded. All structural interference paths are processed in sequence using this method, and all path segment numbers and corresponding coordinate segments that meet the conditions are extracted to form a set of coincident path segment numbers.

[0129] S502: Based on the coordinate information corresponding to the overlapping path segment number set, integrate the path boundary and unstable segment boundary, unify the spatial contour of overlapping segments, and generate spatial weakening patch contour values.

[0130] Based on the extracted set of overlapping path segment numbers, the boundaries of each path segment and the boundaries of unstable sections are spatially integrated. First, the boundary point coordinate data is called to construct a continuous line segment from the start to the end of each path segment. At the same time, the boundary coordinates of the unstable sections that overlap with it are read to form polygonal patches. To unify the spatial representation of the boundaries, the path segment boundaries need to be buffered. The width of the buffer zone is set to 2m. By constructing the buffer zone, a strip-shaped spatial region is formed. Then, spatial overlay analysis is performed with the polygonal boundaries of the unstable sections to identify overlapping areas and extract the boundary line segments of the overlapping parts. These overlapping line segments are spliced, topologically cleaned, and redundant points are removed to finally generate a unified and continuous polygonal contour. This contour expresses its boundary shape with a sequence of coordinate points, including the coordinates of the start point, inflection point, and end point. Each contour record is mapped to its corresponding path segment number and overlapping unstable section number. The same processing is performed on all overlapping path segments in turn. Finally, a set of spatial weakening patch contour values ​​is uniformly output. This set records the boundary shape and position distribution of the weakening region formed by the overlap of each path segment and the unstable section in two-dimensional space.

[0131] S503: Based on the location and shape parameters of the spatial weakening patch outline values, match the grid distribution in the suitability evaluation layer, map the constraint parameters, and obtain the suitability weakening influence area patch;

[0132] Based on the set of spatially weakened patch contour values, a suitability impact matching analysis is performed on their location and shape information. The grid distribution information in the suitability evaluation layer is read, with each grid having a side length of 10m. A grid spatial framework is established by reading the grid center point coordinates and numbers. A spatial matching process is performed on each set of patch contour points to determine whether the contour covers a certain grid area. According to the judgment rules, if any point of the contour falls within the grid boundary area, the grid is considered an affected grid. For each covered grid, its number is read and its affected status is recorded. Simultaneously, the corresponding constraint parameter values ​​are extracted from the layer, and the range of constraint parameters is determined. The value is set to a decimal between 0 and 1, representing the suitability level of the area. For example, if the outline of a map patch falls into three grids numbered G17, G18, and G25, and their corresponding constraint parameter values ​​are 0.7, 0.5, and 0.9, then these three grids are defined as map patches in the suitability weakening influence area. After the correspondence is established, the map patch number, grid number, and constraint parameter value are combined into a structured record. All map patch outline values ​​are processed in sequence, and finally, a set of map patches in the suitability weakening influence area is generated. This set expresses the spatial influence information of the map patch outline on the grid suitability distribution and can be used in subsequent engineering spatial planning and parameter adjustment stages.

[0133] 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 variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones, characterized in that, Includes the following steps: S1: Obtain spatial information of the cross nodes of the crack, construct path combinations, calculate the return rate and node density changes, identify the location of structural abrupt changes, and generate a set of intersecting abnormal path segments; The specific steps of S1 are as follows: S101: Based on the three-dimensional coordinate data of the cross nodes of the crack, extract the spatial connection relationship between the nodes, filter the node pairs that meet the continuous connection condition, construct the path set, record the connection number and index information of the nodes in the path, and obtain the node path combination information. S102: Call the path number and node sequence in the node path combination information, extract the direction change sequence formed by the nodes in the path, calculate the ratio of the direction change amount to the number of path segments in each path, and obtain the path turnaround rate sequence value; S103: Based on the path turnaround rate sequence value and the node coordinates in the node path combination information, extract the density change amplitude of adjacent node segments, filter the position segments that simultaneously meet the conditions of turnaround rate change and density change, and obtain the set of intersecting abnormal path segments; S2: Based on the path coordinates in the set of intersecting abnormal path segments, extract the flow change sequence and head fluctuation sequence of the corresponding monitoring point, calculate the fluctuation intensity, and combine it with the classification standard to generate the range of unstable supply section. The specific steps of S2 are as follows: S201: Based on the path coordinates in the set of intersecting abnormal path segments, extract the monitoring point numbers within the coverage area of ​​the path segments, obtain the flow change sequence and head fluctuation sequence for the corresponding time period, and classify them according to the path segment numbers to obtain the path segment flow and head sequence set. S202: Based on the flow change value and head fluctuation value of the path segment in the sequence of flow and head, extract the change amplitude and fluctuation trend of the sequence, calculate the fluctuation statistics value corresponding to the path segment, and obtain the fluctuation intensity value of the path segment. S203: Call the path segment fluctuation intensity value and path segment number, determine the fluctuation level of the path segment according to the intensity threshold of the fluctuation classification standard, filter the path segment numbers whose fluctuation level is in the unstable range, and obtain the range of the supply unstable segment. S3: Within the unstable supply zone, extract the lithological tag sequence and integrity index along the path, calculate the difference between the lithological transformation frequency and the integrity change rate, locate nodes that exceed the discrimination threshold, and generate a set of lithological structure mutation nodes; The specific steps for S3 are as follows: S301: Based on the path coordinates within the unstable supply zone, extract the lithological label sequence and corresponding integrity index along the path, establish the correspondence between the path location and lithological and integrity data, and obtain the lithological and integrity sequence set along the path; S302: Based on the lithology labels and integrity indices in the lithology and integrity sequence set along the path, the lithology transformation frequency and integrity change rate are statistically analyzed, and the difference between the two is obtained to obtain the lithology transformation frequency and integrity change rate difference sequence. S303: Based on the difference sequence between the lithological transformation frequency and the integrity change rate, call the discontinuity discrimination threshold to determine whether the difference exceeds the limit, mark the path node positions that meet the conditions, and obtain the set of lithological structure mutation nodes; S4: Based on the spatial location of the set of abrupt change nodes in the lithological structure, extract the fracture line density and impermeability grade on both sides of the node, calculate the combination variation characteristics, screen out path segments that meet the interference identification criteria, and generate a set of structural interference paths. S5: Based on the set of structural interference paths, match overlapping areas in the unstable supply segment range, extract and integrate overlapping path segments, generate weakened area patches, map them onto the suitability evaluation layer as a constraint factor, and form suitability weakening influence area patches.

2. The comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones according to claim 1, characterized in that, The set of intersecting abnormal path segments includes path turnback rate, node density variation range, and structural abrupt change location segments. The range of the recharge unstable section includes path coordinates, flow change sequence, head fluctuation sequence, and fluctuation intensity value. The set of lithological structural abrupt change nodes includes lithological label sequence, integrity index, and the difference between lithological transformation frequency and integrity change rate. The set of structural interference paths includes fracture line density, permeability grade, and combination variation characteristics. The suitability weakening influence area patches include the spatial boundary of structural interference path segments, overlapping areas of recharge unstable sections, weakening area patches, and constraint factors in the suitability evaluation layer.

3. The comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones according to claim 1, characterized in that, The specific steps of S4 are as follows: S401: Based on the spatial location provided by the set of abrupt change nodes in the lithological structure, extract the fracture line density and impermeability grade of the path segments on both sides of the node, establish the pairing relationship between node location and associated attributes, and obtain the sequence of fracture and impermeability attributes of the path node. S402: Based on the crack linear density and impermeability grade in the crack and impermeability attribute sequence of the path node, calculate the combination change characteristics, determine the attribute change of the path segments on both sides of the node, and obtain the crack impermeability combination change characteristic sequence. S403: Based on the change feature values ​​of the path segments in the fracture anti-seepage combination change feature sequence, call the interference identification standard threshold, filter the path number segments that meet the conditions, and obtain the set of structural interference paths.

4. The comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones according to claim 1, characterized in that, The interference identification standard threshold is determined by preset criteria for the degree of crack development and the level of change in impermeability. When the degree of crack development on both sides of the path node is significantly enhanced and the level of impermeability changes significantly, it is identified as a structural interference path segment.

5. The comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones according to claim 1, characterized in that, The specific steps of S5 are as follows: S501: Based on the spatial boundary of the path in the set of structural interference paths, match the overlapping positions in the boundary of the unstable supply section, extract the overlapping path number and coordinate section, and obtain the overlapping path segment number set; S502: Based on the coordinate information corresponding to the overlapping path segment number set, integrate the path boundary and the unstable segment boundary, unify the spatial contour of the overlapping segment, and generate the spatial weakening patch contour value. S503: Based on the position and shape parameters of the spatial weakening patch outline values, match the grid distribution in the suitability evaluation layer, map the constraint parameters, and obtain the suitability weakening influence area patch.

6. The comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones according to claim 1, wherein the path combination is a set of paths constructed based on the connectivity between fissure intersection nodes; The turnaround rate is the ratio of the total propagation distance in the path to the straight-line distance between the starting and ending points; The node density change is the gradient of the number of fracture nodes per unit distance. The intersecting abnormal path segments are crack path regions where the change in turnback rate and node density is abrupt, and these path segments are identified as structural disturbance sensitive areas. The flow rate change sequence is a set of flow rate data collected by monitoring points deployed along the fracture path over a continuous time period. The hydraulic head fluctuation sequence is a record sequence of elevation values ​​formed by the change of groundwater level over time. The hydraulic head fluctuation reflects the disturbance response capability of the groundwater system. The fluctuation intensity is a fluctuation index obtained by combining the amplitude of flow rate change and the amplitude of head fluctuation. The grading standard is a classification rule for the wave intensity value based on statistical parameters such as the coefficient of variation, and is determined based on regional hydrological statistical analysis. The range of the unstable supply section is the set of spatial boundaries of the path segments corresponding to the fluctuation intensity values ​​exceeding the stability limit.

7. The comprehensive evaluation method for the suitability of groundwater recharge through seepage in karst fissure zones according to claim 1, characterized in that, The lithological tag sequence is an ordered set composed of lithological codes of rock strata extracted sequentially along the path direction; The integrity index is a parameter used to measure the continuity and stability of rock strata structure, derived from the calculation results of rock strata continuity and tectonic disturbance level data; The lithology transformation frequency is the number of times the lithology label changes per unit path length; The discrimination threshold is a frequency difference limit used to define whether lithological changes have reached the degree of structural abruptness, and is set based on statistical experience values ​​or simulated judgment values. The crack line density is the total length of crack line segments within a unit area; The impermeability grade is a classification grade based on the formation permeability coefficient and fracture distribution, which determines the ability to resist water permeability. The combined variation characteristic is the joint variation trend of fracture linear density and impermeability grade at the location of structural abrupt change; The interference identification standard is a multi-parameter discrimination system used to determine whether a path segment is severely disturbed by structural factors, and is set by fracture characteristics, hydraulic properties, and lithological grade. The set of structural interference paths is a spatial set of path segments that satisfy the interference identification conditions; The overlapping region is the set of locations where two path sets intersect and overlap in spatial coordinates. The weakened region patch is a graphic unit formed by the intersection of the interference path and the unstable supply path. The patch is used to limit the suitability evaluation results and is the representation of the unfavorable factor layer. The suitability evaluation layer is a spatial evaluation unit layer constructed based on the regional groundwater recharge conditions; The suitability reduction effect area patch is an effect area unit formed by overlaying restrictive patches onto the suitability layer.

Images