Coastal wetland groundwater development suitability grading evaluation method and system

Abstract

The invention relates to the technical field of hydrological exploration, in particular to a coastal wetland groundwater development suitability grading evaluation method and system.The method comprises the following steps that groundwater and seawater pressure monitoring devices are arranged to obtain data, gradients are calculated to generate a pressure field distribution diagram, intrusion paths are recognized, peak time difference is extracted, conduction delay is calculated, and sensitivity indexes are obtained; threshold values are compared to delimit risk areas and generate a risk partition map, mining intensity is matched to calculate weights, the weights are input into a fuzzy comprehensive evaluation model to generate evaluation indexes, and mining quantity and permeability coefficients are combined to be input into an analytic hierarchy process for operation to generate a suitability grading result. Underground water and seawater pressure field distribution is revealed through multi-depth and multi-azimuth pressure monitoring, quantitative conduction delay is extracted in combination with time difference, sensitivity indexes reflect intrusion path activity and regional differences, dynamic signals are introduced into risk partitions to match evolution trends, and weight balance mining and permeation are fused in a grading manner. The result layering highlights high-risk warning and low-risk potential.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of hydrological exploration, in particular to a coastal wetland groundwater development suitability grading evaluation method and system. BACKGROUND

[0002] The technical field of hydrological exploration relates to the research and detection of the distribution, recharge, runoff and discharge rules of surface water bodies and groundwater bodies. The core matters of this technical field include groundwater resource exploration, determination of groundwater aquifer characteristics, groundwater dynamic monitoring and analysis of the mutual transformation relationship between surface water and groundwater. The overall technical system of this technical field covers methods such as geophysical exploration, hydrogeological surveying and mapping, water chemical sampling and analysis and numerical simulation, and is used to reveal the occurrence conditions and distribution characteristics of groundwater.

[0003] Among them, the traditional coastal wetland groundwater development suitability grading evaluation method refers to the use of hydrogeological exploration data as the basis, taking aquifer depth, aquifer thickness, water quality type and salinity, recharge and discharge conditions and seawater intrusion sensitivity as the main basis, and combining with water chemical testing, drilling sampling and surface water dynamic observation to conduct grading evaluation when the coastal wetland area is divided into zones and grades for groundwater resource development and utilization.

[0004] The existing technology relies on static geological and water chemical parameters for grading, lacks real-time tracking of the pressure conduction relationship between groundwater and seawater, and in the scenario where seawater intrusion is progressive expansion and latent advancement, parameter changes often lag behind the actual intrusion process, making it difficult to reveal hidden risks in time, and some areas are in a potential invasion state before water quality deterioration occurs, so that the suitability evaluation is disconnected with the actual evolution, making it difficult to accurately grasp the safety limit of the development plan, increasing the risk of overexploitation and salinization. SUMMARY

[0005] In order to solve the technical problem that the existing technology relies on static geological and water chemical parameters for grading, lacks real-time tracking of the pressure conduction relationship between groundwater and seawater, and in the scenario where seawater intrusion is progressive expansion and latent advancement, parameter changes often lag behind the actual intrusion process, making it difficult to reveal hidden risks in time, and some areas are in a potential invasion state before water quality deterioration occurs, so that the suitability evaluation is disconnected with the actual evolution, making it difficult to accurately grasp the safety limit of the development plan, increasing the risk of overexploitation and salinization, the present application provides a coastal wetland groundwater development suitability grading evaluation method, comprising the following steps: In order to achieve the above purpose, the technical scheme adopted by the present application is as follows: a coastal wetland groundwater development suitability grading evaluation method, comprising the following steps: S1: Obtain groundwater pressure by laying multiple-depth groundwater pressure monitoring devices, obtain seawater surface pressure by laying seawater surface pressure monitoring devices, calculate the gradient of adjacent depth monitoring points, calculate the second-order gradient of adjacent horizontal monitoring points, and generate a pressure field distribution map; S2: Identify the seawater intrusion advancing path based on the pressure field distribution map, extract the pressure change peak time difference of the monitoring points, calculate the pressure conduction delay characteristics, and obtain the sensitivity index; S3: Compare the sensitivity index with a preset critical threshold value, and if the sensitivity index exceeds the threshold value, it is determined as a high-risk area, and if the sensitivity index is below the threshold value, it is determined as a safe area; spatially classify according to the numerical interval, and generate a risk zoning map; S4: Match the risk zoning with the allowable exploitation intensity standard based on the risk zoning map, calculate the development suitability evaluation weight coefficient by calculating the difference of seawater intrusion risk in each zone, input the weight coefficient into the fuzzy comprehensive evaluation model for comprehensive evaluation operation, and generate a suitability evaluation index; S5: Obtain the exploitation amount and the permeability coefficient based on the suitability evaluation index, input the suitability evaluation index, the exploitation amount, and the permeability coefficient into the analytic hierarchy process for weight distribution operation, and generate a suitability classification result.

[0006] As a further scheme of the present application, the pressure field distribution map includes spatial distribution pattern, abnormal change area and distribution boundary characteristics, the sensitivity index includes pressure conduction coefficient, invasion path difference and pressure response time difference, the risk zoning map includes grade division unit, regional boundary range and zoning attribute characteristics, the suitability evaluation index includes resource carrying capacity, environmental stability and utilization feasibility, and the suitability classification result includes priority exploitation area, restricted exploitation area and prohibited exploitation area.

[0007] As a further scheme of the present application, the specific steps of S1 are as follows: S101: Perform index sorting and difference operation on the pressure values of adjacent depth monitoring points based on the collected data of the groundwater pressure monitoring devices and the seawater surface pressure monitoring devices, aggregate the difference results in the time sequence dimension, and generate a depth pressure gradient sequence; S102: Call the depth pressure gradient sequence, perform first-order difference processing on the same layer data values of adjacent horizontal monitoring points, then perform amplitude operation on the difference results to obtain change trend values, and combine them into continuous distribution in a two-dimensional coordinate system to generate a horizontal pressure gradient sequence; S103: According to the corresponding relationship between the horizontal pressure gradient sequence and the depth pressure gradient sequence, perform coordinate mapping and numerical superposition, and discretize the superposition results in the spatial grid to obtain the pressure field distribution map.

[0008] As a further scheme of the present application, the specific steps of S2 are as follows: S201: Based on the pressure field distribution map, the pressure data sequence of the monitoring point is retrieved point by point, the local peak amplitude and time index are calibrated, the peak value difference of adjacent points is logically judged, and the advancing path sequence is generated; S202: The pressure curve of the monitoring point in the advancing path sequence is called, the peak value occurrence time is extracted, and the peak value time of adjacent points is compared. According to the comparison result, the time difference value is calculated, and the difference value is combined into a data sequence to generate a time difference sequence; S203: According to the numerical distribution of the time difference sequence, the sequence elements are weighted and aggregated, and the result is subtracted from the delay reference value. According to the difference value, the normalized coefficient is calculated and assigned to the monitoring point level to obtain the sensitivity index.

[0009] As a further scheme of the present application, the specific steps of S3 are: S301: Based on the sensitivity index and the preset critical threshold parameter, the sensitivity index value is retrieved, the difference operation is performed on the value and the critical threshold value, and if the value is greater than the threshold value, the corresponding space unit is marked as a risk unit, and if the value is less than the threshold value, the space safety unit is marked as a space safety unit, and the space unit classification identification is obtained; S302: The space unit classification identification is called, the risk unit and the safety unit index are aggregated, the sensitivity index distribution interval is mapped to the level range, and the numerical interval is sorted in ascending order to generate a space classification index sequence; S303: According to the space classification index sequence, the corresponding space unit is located on the two-dimensional space coordinates, the boundary of the unit range of different numerical intervals is drawn, and the interval level index is filled into the corresponding position to generate a risk partition map.

[0010] As a further scheme of the present application, the preset critical threshold is a numerical limit preset according to the sensitivity index value range; The space classification index sequence is a level index sequence formed by dividing the space unit according to the sensitivity index interval and arranging it in ascending order according to the value.

[0011] As a further scheme of the present application, the specific steps of S4 are: S401: Based on the distribution data in the risk partition map, the risk level value is compared with the allowed mining intensity threshold value level by level, the difference between the level and the threshold value is retrieved in the comparison process, and the difference is converted into a regional mining intensity corresponding relationship set to generate a risk intensity matching result; S402: The matching parameters in the risk intensity matching result are called, the partition level value and the seawater intrusion risk difference coefficient are weighted, the weighted coefficient is calculated according to the numerical deviation of the difference coefficient and the partition risk reference value, and the development suitability evaluation weight coefficient is obtained. S403: calling the development suitability evaluation weight coefficient, inputting the partition weight into the fuzzy comprehensive evaluation model, performing numerical aggregation operation according to the distribution proportion of the weight on the multi-factor index, and obtaining the suitability evaluation index.

[0012] As a further scheme of the present application, the allowable exploitation intensity threshold refers to the allowable exploitation limit of groundwater in a specific area according to the hydrogeological conditions and resource carrying characteristics. The development suitability evaluation weight coefficient refers to a partition evaluation weight numerical coefficient calculated by the partition grade value, the risk difference coefficient and the reference value deviation.

[0013] As a further scheme of the present application, the specific steps of S5 are: S501: based on the suitability evaluation index, comparing the index parameter value with the corresponding threshold parameter one by one, recording the index value exceeding the threshold parameter, and sequentially integrating the recorded index value according to the belonging category to generate index sequence data; S502: calling the index sequence data and aggregating with the exploitation amount data and the permeability coefficient data, scaling and converting the three types of data in the same numerical interval into a vector structure according to the relative difference between the data, and analyzing the allocation proportion factor to obtain a weight coefficient vector; S503: according to the weight coefficient vector, weighting and summing the corresponding elements of the index sequence data, exploitation amount data and permeability coefficient data, and then segmenting and sorting according to the numerical interval of the operation result to generate suitability classification results.

[0014] The coastal wetland groundwater development suitability classification evaluation system comprises: A groundwater pressure monitoring module obtains groundwater pressure by arranging multi-depth groundwater pressure monitoring devices and obtains seawater surface pressure by arranging seawater surface pressure monitoring devices, performs gradient calculation on adjacent depth monitoring points and quadratic gradient calculation on adjacent horizontal monitoring points to generate a pressure field distribution map, which is transmitted to a seawater intrusion identification module; The seawater intrusion identification module identifies the seawater intrusion advancing path based on the pressure field distribution map, extracts the pressure change peak time difference of the monitoring points, calculates the pressure conduction delay characteristics, obtains a sensitivity index, and transmits it to a risk assessment classification module; The risk assessment classification module compares the sensitivity index with a preset critical threshold value, and if the threshold value is exceeded, it is classified as a high-risk area, and if the threshold value is lower, it is classified as a low-risk area, and performs spatial grid classification operation according to the sensitivity coefficient numerical interval to generate a risk distribution map; The suitability comprehensive evaluation module matches the risk zones with the allowable mining intensity standards based on the risk zoning map, calculates the development suitability evaluation weight coefficients through the differences in seawater intrusion risk between zones, inputs the weight coefficients into the fuzzy comprehensive evaluation model for comprehensive evaluation calculation, generates suitability evaluation indicators, and transmits them to the development intensity grading module. The development intensity grading module obtains the extraction volume and permeability coefficient based on the suitability evaluation index, performs weight allocation calculation based on the suitability evaluation index data, recommended extraction volume, and aquifer permeability coefficient, and divides the development intensity level according to the numerical range of the comprehensive evaluation results to generate development suitability grading results.

[0015] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, the pressure field distribution between groundwater and seawater is obtained through multi-depth and multi-directional pressure monitoring, and the pressure transmission delay characteristics are quantified by combining time difference extraction. This allows the sensitivity index to intuitively reflect the activity level and regional differences of the intrusion path. Dynamic response signals are added when partitioning the risk space to ensure that the partitioning results are consistent with the actual evolution trend. Furthermore, risk weights and partition intensity matching are introduced in the grading process, so that the suitability evaluation can be comprehensively balanced by combining the amount of extraction and permeability conditions. This results in a stratified characteristic, which can highlight the warning of high-risk areas and also highlight the utilization potential of low-risk areas. Attached Figure Description

[0016] 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.

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

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

[0019] 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.

[0020] 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.

[0021] 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.

[0022] 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.

[0023] Please see Figure 1 This invention provides a method for grading and evaluating the suitability of groundwater development in coastal wetlands, comprising the following steps: S1: Groundwater pressure is obtained by deploying multi-depth groundwater pressure monitoring devices, sea surface pressure is obtained by deploying sea surface pressure monitoring devices, gradient calculation is performed on adjacent depth monitoring points, secondary gradient calculation is performed on adjacent horizontal monitoring points, and pressure field distribution map is generated. S2: Identify the seawater intrusion propagation path based on the pressure field distribution map, extract the peak time difference of pressure changes at monitoring points, calculate the pressure transmission delay characteristics, and obtain the sensitivity index; S3: The sensitivity index is compared with the preset critical threshold. Areas exceeding the threshold are classified as high-risk areas, and areas below the threshold are classified as low-risk areas. Spatial classification is performed according to the numerical range to generate a risk zoning map. S4: Based on the risk zoning map, the risk zones are matched with the allowable mining intensity standard. The development suitability evaluation weight coefficient is calculated by the difference in seawater intrusion risk between zones. The weight coefficient is input into the fuzzy comprehensive evaluation model for comprehensive evaluation and calculation to generate suitability evaluation index. S5: Based on the suitability evaluation index, obtain the extraction volume and permeability coefficient. Input the suitability evaluation index, extraction volume and permeability coefficient into the analytic hierarchy process for weight allocation calculation to generate the suitability classification result.

[0024] The pressure field distribution map includes spatial distribution patterns, areas of abnormal changes, and distribution boundary characteristics. The sensitivity index includes pressure transmission coefficient, intrusion path differences, and pressure response time difference. The risk zoning map includes grade division units, regional boundary ranges, and zoning attribute characteristics. The suitability evaluation indicators include resource carrying capacity, environmental stability, and utilization feasibility. The suitability classification results include priority mining areas, restricted mining areas, and prohibited mining areas.

[0025] Please see Figure 2 The specific steps of S1 are as follows: S101: Using data collected from groundwater pressure monitoring devices and sea surface pressure monitoring devices, index sorting and difference calculation are performed on the pressure values ​​of adjacent depth monitoring points, and the difference results are aggregated in the time series dimension to generate a depth pressure gradient sequence. Using data collected from groundwater pressure monitoring devices and sea surface pressure monitoring devices, multiple vertically deployed groundwater pressure monitoring points were first selected in the coastal wetland area. Each monitoring point was located at a different depth, such as 5m in the shallow layer, 15m in the middle layer, and 30m in the deep layer. Sea surface pressure monitoring points were also deployed at adjacent locations to simultaneously record sea surface pressure values. The pressure values ​​were continuously recorded at a frequency of 1 hour per data point, and the groundwater pressure data at the same time point were matched using timestamps. Sea surface pressure data Establish a one-to-one correspondence and sort the data by depth from shallow to deep, for example, at the same time... The groundwater pressures are as follows: , , ; Corresponding sea surface pressure When calculating the pressure difference between adjacent depth monitoring points, a subtraction operation is performed: , ; Simultaneously, the difference between groundwater and sea level is directly subtracted, such as the difference in shallow layers. Mid-layer difference Deep interpolation After obtaining the difference at a single moment, the difference results from multiple consecutive moments (e.g., 24 sets within a day) are aggregated according to time series. Specifically, at the same depth, the pressure difference over time is stored as a time series array: ; It also performs statistical operations such as mean, extreme values, and variance on the sequence. For example, the mean of the deep data corresponding to 24 hours is calculated as follows: ; Extreme value The variance is The statistical values ​​can provide direct input when generating depth pressure gradient sequences. Determining the depth pressure gradient requires inputting the mean of adjacent depth difference sequences into the gradient calculation logic. For example, in the 5-15m range, the pressure gradient value... In the 15-30m range, the pressure gradient value At the same time, it is necessary to determine the gradient level by setting a numerical range based on the special characteristics of the interaction between groundwater and sea surface in coastal wetlands. For example, the gradient should be set within a certain range. Defined as low gradient, Defined as the meso gradient, Defined as a high gradient, according to this standard, the gradients of the two intervals mentioned above can be determined to fall within the medium gradient range. The calculated depth pressure gradient sequence is sorted and stored with time as the horizontal axis and depth as the vertical axis, and is ultimately used as the basic data for spatial pressure field calculation, forming a depth pressure gradient sequence that can be used directly.

[0026] Table 1: Initial pressure gauges at monitoring points

[0027] As shown in Table 1, by using monitoring points at different vertical depths and a unified reference sea surface pressure, the pressure difference between adjacent depths and the difference between depth and sea surface can be obtained. The results can be aggregated into a complete depth pressure gradient sequence after being arranged in time series.

[0028] S102: Call the depth pressure gradient sequence, perform first-order difference processing on the data values ​​of adjacent horizontal monitoring points in the same layer, then perform amplitude calculation on the difference results to obtain the trend value, and combine them into a continuous distribution in the two-dimensional coordinate system to generate the horizontal pressure gradient sequence. Using a depth pressure gradient sequence, adjacent monitoring points distributed horizontally are selected within the same coastal wetland area. Each monitoring point is positioned at the same depth layer. For example, three horizontal monitoring points, A, B, and C, are positioned at a depth of 15m. First, the pressure gradient values ​​of multiple points at the same time are obtained, such as point A. Point B Point C The first-order difference process is performed according to the spatial layout order, that is, the difference value is obtained by directly subtracting two adjacent points. , After obtaining the difference value, to reflect the actual range of change in the difference, the absolute value of the difference result is taken to obtain the range of change value. , Then, the change magnitude values ​​at the same time are grouped into a horizontal change trend array. To facilitate time tracking, the trend arrays of changes at different time points (e.g., 24 monitoring times per day) are stacked and stored in chronological order to form a two-dimensional matrix. Rows in the matrix correspond to the time series, and columns correspond to adjacent monitoring point segments. During matrix filling, if the change in a certain segment exceeds a preset threshold (e.g., the threshold is set to...), the matrix will be updated accordingly. If the value is 1, it is marked as 1 in the matrix; otherwise, it is marked as 0. This binarization can quickly determine the spatial location of significant horizontal changes. After obtaining the horizontal change trend, it is associated with the XY coordinates of the monitoring point in the real geographic coordinate system, such as point A (0, 0), point B (1, 0), and point C (2, 0). The change trend value is then filled into the corresponding spatial location in the two-dimensional coordinate system. The trend value of the discrete monitoring points is used to generate a continuous distribution field through interpolation, thus obtaining the spatial distribution of the horizontal pressure gradient. The resulting horizontal pressure gradient sequence not only includes the change amplitude of the monitoring segment and its temporal evolution, but is also arranged according to spatial relationships to facilitate direct superposition calculation with the depth pressure gradient data.

[0029] S103: Based on the correspondence between the horizontal pressure gradient sequence and the depth pressure gradient sequence, coordinate mapping and numerical superposition are performed, and the superposition result is discretized in a spatial grid to obtain a pressure field distribution map. Based on the correspondence between horizontal pressure gradient sequences and depth pressure gradient sequences, the two types of data are first time-stamped together under the same time reference to ensure that spatial locations have paired horizontal gradient values ​​at the same time. With depth gradient value Taking a monitoring point M in a coastal wetland as an example, let's assume that in The horizontal gradient value at time t is The depth gradient value is When performing coordinate mapping, the planar coordinates of monitoring point M are first called. with depth coordinates ,For example , , The horizontal and depth gradients are mapped to the corresponding 3D coordinates and then numerically superimposed to obtain the combined gradient value. Perform the same operation on the monitoring points, for example, point N at the same time. , The combined gradient value is After obtaining the comprehensive gradient values ​​of all monitoring points at this moment, they are discretized, that is, according to the preset spatial grid units (such as...). The continuous area is divided into multiple grids, and the combined gradient values ​​of the monitoring points falling into each grid are spatially interpolated. If a cell has only one monitoring point, the combined gradient value of that point is directly assigned. If a cell has multiple monitoring points, the arithmetic mean is calculated. For example, if the combined gradients of two points in a certain cell are... and Then the gradient value of that grid is The calculations are performed sequentially on the grid to form a comprehensive pressure field data matrix covering the study area. The matrix is ​​then plotted as a two-dimensional or three-dimensional spatial distribution map, so that the pressure differences in different areas are fully mapped in the spatial grid. The resulting pressure field distribution map can show the spatial organization structure of the pressure field at multiple locations at the same time, and can be directly used as the input basic data for groundwater development suitability classification.

[0030] Please see Figure 3 The specific steps of S2 are as follows: S201: Based on the pressure data sequence of monitoring points in the pressure field distribution map, the curve shape of the monitoring points is retrieved point by point, the local peak amplitude and time index are calibrated, the difference of peak values ​​between adjacent points is logically judged, and the propulsion path sequence is generated. Based on the pressure data sequence of monitoring points in the pressure field distribution map, the spatial coordinates of each monitoring point in the map are first established to correspond to its corresponding pressure data curve. For example, monitoring point M is located at coordinates... Its 24-hour pressure data series is as follows: ; The data is arranged in chronological order and includes hourly timestamps. In the process of point-by-point retrieval of curve shapes, the current point of each time series is selected and directly compared with a fixed number of data points before and after it. In this embodiment, the comparison window is one point before and after. If the pressure value of the current point is greater than the value of the points before and after, it is identified as a local peak. For example, in Constant pressure The previous value was The value after that is ,because and Therefore, this point is marked as the peak point. When calibrating the local peak amplitude, the peak value is directly subtracted from the pressure value of the previous monitoring point on the same curve. For example, the peak amplitude... At the same time, a time index number is appended to the peak record. As a calibration record of peak occurrence time, if multiple peaks occur at the same monitoring point, this calibration process is repeated to obtain multiple peak amplitudes and their associated time indices. After obtaining all peak calibration data for a single monitoring point, the logical judgment stage of peak differences between adjacent points is entered. The peak amplitudes of two adjacent geographically close monitoring points (e.g., M and N, 200m apart) under the same time index are directly subtracted. Assuming that in... The peak amplitude at point M is The peak amplitude at point N is The difference is If the difference is greater than the set difference threshold The direction is then marked as a valid propulsion direction on the propulsion path determination matrix; otherwise, it is marked as no propulsion direction. This threshold is set with reference to the range of normal pressure difference fluctuations between different monitoring points in the coastal wetland area, and its reasonable range, derived from raw monitoring data, is 0.3–0.6. The median value was selected as 0.4. As a judgment boundary, the process traverses monitoring point pairs that have spatial adjacency, accumulates effective advancement direction segments one by one in the time dimension, forms an array of advancement paths arranged in time order, and finally outputs a sequence of monitoring point groups sorted by time as the advancement path sequence.

[0031] S202: Call the pressure curves of monitoring points in the advancement path sequence, extract the time of peak occurrence, compare the peak times of adjacent points, calculate the time difference based on the comparison results, and combine the differences in order into a data sequence to generate a time difference sequence; The pressure curves of monitoring points in the propulsion path sequence are retrieved. First, the peak information of the already calibrated points is retrieved. Each propulsion path unit consists of two or more points, such as a path segment. The pressure curve includes points M and N, and the peak value array at point M is... The N-point peak array is When extracting the peak occurrence time, the time index is converted to precise time (hours), for example... , , Similarly, for point N, we obtain... Then, the two points are paired and compared according to the peak order. The matching rule is that the first peak of point M is paired with the peak of point N that is closest in time and does not exceed a certain time limit after it (the time limit is set to 3 hours based on the response characteristics of the monitoring area). The time difference of the pairing is calculated. For example, point M Corresponding to N points When, then the time difference is Point M Corresponding to N points When, the difference is Point M Corresponding to N points When, the difference is If, during the pairing process, a peak occurs at point M and then no corresponding peak occurs at point N within the time limit, the pairing is discarded and will not participate in subsequent calculations. After obtaining all valid time differences for a path segment, they are arranged into an array in chronological order, as in this example. This array represents the time difference sequence of the path segments. The same steps are then performed sequentially on the remaining path segments in the advancement path sequence, for example, path segment... The obtained time difference sequence is Once the path segment completes this processing flow, the path segment time difference sequences are merged and output as a complete time difference sequence for the next step.

[0032] S203: Based on the numerical distribution of the time difference sequence, the sequence elements are weighted and aggregated, and the result is compared with the delay benchmark value. The normalization coefficient is calculated based on the difference result and the monitoring point is assigned a level to obtain the sensitivity index. Based on the time difference series, its numerical distribution is first analyzed. For example, the time difference data of the entire region is merged as follows: In the weighted aggregation process, a weight coefficient is assigned to each time difference. The weighting is set with reference to the influence of time difference on the regional hydraulic conduction characteristics. For example, the weighting is set for time differences greater than or equal to 0 and less than or equal to 1. The weight of the interval is set as The time difference is greater than 1 and less than or equal to 2. The weight of the interval is set as The time difference is greater than 2 and less than or equal to 3. The weight of the interval is set as Time difference greater than The weight is set to Then, the time difference is multiplied by its corresponding weight and summed. For example, the weighted sum of this sequence is... The weighted aggregation result is obtained. Compare the result with the delay baseline value The difference calculation is performed, with the delay reference value set based on the comparison results of the original simulation and actual measurements, and the range is [range missing]. This example is set up The difference is When calculating the normalization coefficient, the normalization denominator is first determined to be the range of the time difference series of the entire sample in the region. Then, take the absolute value of the difference and normalize the coefficient. This coefficient is directly used for grade classification, and the grade classification is based on: Level 1 Level 2 Level 3 It is level four, therefore this example If it belongs to Level 1, assign this level to the corresponding monitoring point and calculate its sensitivity index. The sensitivity index value can be obtained by... Multiply by 100 and divide by percentage to obtain the sensitivity index in this example. .

[0033] Table 2: Calculation Table of Sensitivity Index of Monitoring Points

[0034] As shown in Table 2, the weighted aggregate value can be obtained by calculating the time difference and weight coefficient one by one. After the delay benchmark value difference calculation and normalization processing, the sensitivity index and level of each monitoring path segment and the whole area are obtained, which is the direct basis for the final sensitivity index result.

[0035] Please see Figure 4 The specific steps of S3 are as follows: S301: Based on the sensitivity index and the preset critical threshold parameter, retrieve the sensitivity index value, perform a difference operation between the value and the critical threshold, if the value is greater than the threshold, mark the corresponding spatial unit as a risk unit, if it is less than the threshold, mark it as a spatial safe unit, and obtain the spatial unit classification identifier. Based on the sensitivity index and a preset critical threshold parameter, a set of sensitivity index values ​​corresponding to the spatial unit is retrieved. Each element in this set is associated with its spatial unit ID and planar coordinate position. For example, the sensitivity index of unit A is... The sensitivity index of unit B is The sensitivity index of unit C is And ensure that the units are dimensionless percentage results. During the retrieval process, the sensitivity index value of each unit must be explicitly extracted one by one. and with the set critical threshold When performing direct subtraction, the determination of the critical threshold needs to refer to the risk definition standards for coastal wetland areas. For example, based on the analysis of the original data, the unit exhibits a high-risk response tendency when the sensitivity index is greater than 150. Therefore, [the following is a possible interpretation:] As a judgment criterion, in this embodiment, differential processing is performed for unit A. For unit B, perform differential operation. For unit C, perform differential operation. Then, in the logical judgment stage, each difference result is evaluated for its numerical range. If the difference result is greater than 0, the unit is determined to be a risk unit; if the difference result is less than or equal to 0, the unit is determined to be a safe space unit. For example, if the difference of unit A is negative, it is determined to be a safe space unit; if the difference of units B and C is positive, they are marked as risk units respectively. When constructing the spatial unit identification data structure, the classification label of each unit is bound to its spatial coordinates and number to form a structured array, for example... , , This array can be stored sequentially or sorted by geographical location, making it easy to read and call in subsequent steps, and finally generating spatial unit classification labels with "risk" or "safety" classification tags.

[0036] S302: Call the spatial unit classification identifier, perform aggregation processing on the risk unit and safety unit indexes, map the index to the level range according to the sensitivity index distribution range, sort in ascending order according to the numerical range, and generate a spatial hierarchical index sequence; To retrieve the spatial unit classification identifier, first extract the spatial unit IDs of risk units and safe units according to the classification label, forming two independent index sets, such as the risk unit set. Security Unit Set For each unit within a set, its corresponding sensitivity index value is further retrieved and stored in a corresponding array, such as a risk unit sensitivity index array. Security Unit Sensitivity Index Array Then, an aggregation operation is performed on the indexes in the risk unit set, that is, the IDs of the risk units are grouped according to the sensitivity index distribution range. The sensitivity index distribution range is set according to the risk level mapping rules. For example, the range 0-50 is mapped to level 1, 50-150 is mapped to level 2, 150-250 is mapped to level 3, and greater than 250 is mapped to level 4. In this embodiment, risk unit B (226.7) is mapped to level 3, and risk unit C (316.7) is mapped to level 4. Similarly, the same range mapping is performed on the safety unit set. For example, unit A (16.7) is mapped to level 1, and unit D (85.4) is mapped to level 2. The above mapping results are recombined into a spatial unit index sequence arranged in order of level. The sorting is done in ascending order, that is, first list the level 1 unit IDs, then list the level 2, level 3, and level 4 unit IDs, for example, forming a sequence. The system maintains the level label associated with each ID in its storage, and the resulting spatial hierarchical index sequence can reflect the level distribution after unit classification and can also be directly used for two-dimensional spatial coordinate positioning and drawing.

[0037] S303: Based on the spatial hierarchical index sequence, locate the corresponding spatial unit on the two-dimensional spatial coordinates, perform boundary drawing on the unit range of different numerical intervals, and then fill the interval level index into the corresponding position to generate a risk zoning map. Based on the spatial hierarchical index sequence, the coordinates of the geometric center point of each spatial unit are located one by one in the two-dimensional spatial coordinate system. For example, unit A is located at... Unit D is located Unit B is located in Unit C is located in During the localization process, it is necessary to retrieve the set of polygon vertex coordinates corresponding to the unit ID from the spatial base map data, and use the geometric center point as the filling reference position. Units belonging to the same sensitivity level are aggregated, that is, it is determined whether their geometric boundaries are adjacent or touching. If they are adjacent, they are merged into a continuous polygon block. Then, closed polylines are drawn on the outer boundary of the block according to the actual coordinate order to form the boundary graphic. For discontinuously distributed units, independent boundaries are drawn separately. After the boundary drawing is completed, the level value of each unit in the spatial hierarchical index sequence is filled into the corresponding spatial range according to the mapped color code or mark symbol. For example, level 1 is numbered "1" in light green, level 2 is numbered "2" in light yellow, level 3 is numbered "3" in orange, and level 4 is numbered "4" in red, so as to distinguish the regional distribution of different levels on the two-dimensional graphic. After the drawing is completed, it is stored in a readable vector or raster data format to form a spatial partition map dataset containing hierarchical information.

[0038] Table 3: Two-Dimensional Positioning and Hierarchical Allocation of Spatial Units

[0039] As shown in Table 3, when drawing spatial partitions, the geometric position, boundary information and level markers of each unit are matched and encoded one by one. The generated table can be directly used as important input data in the risk partitioning map generation process.

[0040] Please see Figure 5 The specific steps of S4 are as follows: S401: Based on the distribution data in the risk zoning map, the risk level values ​​are compared with the allowable mining intensity thresholds level by level. During the comparison process, the differences between the level and the threshold are retrieved and the differences are converted into a set of regional mining intensity correspondences to generate risk intensity matching results. Based on the distribution data in the risk zoning map, the risk level values ​​of spatial units are first extracted one by one from the zoning map. This value is a level index; for example, unit A is level 1, unit B is level 3, unit C is level 4, and unit D is level 2. The level value range is... The set of integers, followed by retrieval of the set of allowed mining intensity threshold parameters. Each threshold corresponds one-to-one with a risk level. Parameter settings refer to the regional extraction control standards for groundwater management. For example, Level 1 threshold... Level 2 threshold Level 3 threshold Level 4 threshold During the comparison process, for each unit Obtain its risk level Corresponding allowable mining intensity threshold Then, in conjunction with the proposed mining intensity of the unit. Perform differential calculations; for example, the current proposed mining intensity for unit A (level 1) is... ,difference Unit B (Level 3) Proposed Mining Intensity ,difference Unit C (Level 4) Proposed Mining Intensity ,difference Unit D (Level 2) Proposed Mining Intensity ,difference Then, the above difference results are stored in the difference result list according to the unit number. At the same time, a corresponding regional mining intensity relationship identifier is defined for each difference value. When the difference result is greater than 0, it is marked as "exceeding the threshold"; when it is equal to 0, it is marked as "equal to the threshold"; and when it is less than 0, it is marked as "below the threshold". For example, units B and C are marked as exceeding the threshold, and units A and D are marked as below the threshold. Finally, the unit ID and risk level are entered into the list. Permissible mining intensity threshold Proposed mining intensity Difference The combination of relationship identifiers forms a set of relationships corresponding to regional mining intensity, resulting in risk intensity matching results.

[0041] S402: Call the matching parameters in the risk intensity matching results, calculate the weighted difference coefficient between the zoning level value and the seawater intrusion risk difference coefficient, and calculate the weighted coefficient based on the numerical deviation between the difference coefficient and the zoning risk benchmark value to obtain the development suitability evaluation weight coefficient. To retrieve the risk intensity matching results, first read the partition level value of each spatial unit one by one. and its corresponding seawater intrusion risk difference coefficient in the matching results. ,in The value of risk is defined as the ratio of the absolute value of the difference between the proposed mining intensity and the allowable mining intensity threshold of a unit to the allowable mining intensity threshold. It is expressed in dimensionless form as the relative difference in risk. For example, for unit A: Threshold ,but Unit B: Threshold , Unit C: Threshold , Unit D: Threshold , Then, the partition level values ​​and their corresponding difference coefficients are weighted and calculated. The setting of reference partition levels depends on their importance in resource utilization; for example, level 1 weight. Level 2 weights 3-level weight 4-level weight Calculate the weighted results sequentially by unit. Unit A is obtained: Unit B: Unit C: Unit D: Subsequently, a risk benchmark value for each level was introduced. This value represents the acceptable difference coefficient threshold under conditions where significant seawater intrusion does not occur, for example: Level 1. Level 2 Level 3 Level 4 Then calculate the weighting coefficient for each unit. For example, unit A: Unit B: Unit C: Unit D: , will the unit The dataset is associated with its spatial ID, grade value, difference coefficient, and weight value to obtain the development suitability evaluation weight coefficient dataset.

[0042] S403: Call the development suitability evaluation weight coefficient, input the partition weight into the fuzzy comprehensive evaluation model, and perform numerical aggregation operation according to the weight distribution ratio on the multi-factor indicators to obtain the suitability evaluation indicators. By accessing the suitability assessment weight coefficient dataset, the first step is to retrieve the corresponding set of multi-factor evaluation indicators for each spatial unit, such as the hydrogeological condition score. Groundwater recharge conditions score Ecological sensitivity score The scores are all quantized within the range of 0–1 and are the results after dimensionless processing. In this embodiment, the multi-factor scores of unit A are set as follows: , , Unit B is , , Unit C is , , Unit D is , , Then the obtained evaluation weight coefficients The scores are aggregated proportionally to the indicator scores. The calculation method is to calculate the product of each factor score and its weight coefficient in turn, and then sum all the products. For example, in unit A: Unit B: Unit C: Unit D: The sum obtained from the aggregation operation is used as the suitability evaluation index. Store and retain the ID, level, weight, and factor score of the spatial unit to form a complete suitability evaluation index data table.

[0043] Table 4: Development Suitability Evaluation Indicators

[0044] As shown in Table 4, the suitability evaluation index for each spatial unit is obtained by proportionally calculating the weight coefficients and multi-factor scores. This data table can be directly used as the basis for subsequent comprehensive analysis of regional development suitability.

[0045] Please see Figure 6 The specific steps of S5 are as follows: S501: Based on suitability evaluation indicators, compare the indicator parameter values ​​with the corresponding threshold parameters one by one, record the indicator values ​​that exceed the threshold parameters, and serialize and integrate the recorded indicator values ​​according to their categories to generate indicator sequence data; Based on suitability evaluation indicators, the suitability evaluation indicator values ​​for each spatial unit are first extracted from the dataset. The index values ​​are currently in dimensionless percentage form, with a range of [0, 1]. For example, unit A is 0.17888, unit B is 0.03800, unit C is 0.16500, and unit D is 0.18034. Then, the evaluation threshold parameters set in the configuration file are retrieved. This parameter is determined with reference to the lower limit standard of regional development suitability, for example, in the assessment, when Then we can proceed to the next step of resource allocation assessment, therefore we set... Then for each unit and The comparison is performed by direct subtraction. And record the units with positive differences and their index values. For example, unit A: (Greater than 0, record, index value 0.17888), Unit B: (Not recorded), Unit C: (Record, index value 0.16500), Unit D: (Record, index value 0.18034), then the recorded threshold-exceeding indicators are serialized and integrated according to their categories. The category division is based on the suitability assessment type definition. In this embodiment, there are four categories: hydrogeological conditions, groundwater recharge conditions, ecological sensitivity, and comprehensive. The comprehensive category is used to store the summary of all threshold-exceeding indicators. During the serialization process, for the index values ​​of units A, C, and D, the corresponding classification arrays are inserted according to their respective categories, and the arrays are arranged in ascending order by unit ID. For example, the comprehensive category sequence is... For hydrogeological condition sequences, only the threshold values ​​that meet the conditions of that type are retained. For example, if A and D meet the conditions, then the sequence... Groundwater recharge conditions and ecological sensitivity are processed in the same way, ultimately generating index sequence data containing ordered arrays of multiple categories and their over-threshold indicators.

[0046] S502: Call the indicator sequence data and aggregate it with the mining volume data and permeability coefficient data. Scale the aggregated three types of data within the same numerical range and convert them into a vector structure. Analyze and allocate the proportional factor based on the relative differences between the data to obtain the weight coefficient vector. To retrieve the indicator sequence data, first read the indicator values ​​from the comprehensive sequence, hydrogeological condition sequence, groundwater recharge condition sequence, and ecological sensitivity sequence one by one, and then match them with the extraction volume data of that unit according to the unit ID. and permeability coefficient data When performing paired aggregation, ensure that the numerical order of the three types of data is completely consistent during pairing. For example, the index value of unit A is... Mining volume Permeability coefficient The index value of unit C is Mining volume Permeability coefficient The index value of unit D is Mining volume Permeability coefficient After obtaining the paired dataset, the three data categories are mapped to the same numerical interval [0, 1]. The scaling method involves calculating the minimum and maximum values ​​for each data category and then performing interval normalization. For example, for the index value category: minimum value... Maximum value The scaling result for cell A is: ; Unit C is 0.000, and unit D is 1.000; for the extraction volume category: minimum value 600, maximum value 1500, unit A result is 1.000, unit C is 0.000, and unit D is 0.444; for the permeability coefficient category: minimum value 38, maximum value 50, unit A result is 0.583, unit C is 0.000, and unit D is 1.000. After scaling, the three types of normalized data for each unit are combined in vector order to form a three-dimensional vector structure. For example, the vector for unit A is [0.905, 1.000, 0.583], unit C is [0.000, 0.000, 0.000], and unit D is [1.000, 0.444, 1.000]. Then, based on the relative differences between the components within the vectors of different units, the proportion of each vector component is calculated and assigned a scaling factor. For example, the scaling factor group for unit A is... The same steps are performed on the units. When the vector sum is detected to be close to 0, weights are assigned according to the relative performance of the unit in the original data. For example, the original data of unit C is: index value 0.16500, mining volume 600, permeability coefficient 38. Among the three dimensions, the index value is relatively high (close to the threshold), so the index value is given a higher weight, and the weights are assigned as [0.600, 0.200, 0.200], resulting in a set of weight coefficient vectors arranged by unit ID.

[0047] S503: Based on the weight coefficient vector, the corresponding elements of the indicator sequence data, extraction volume data and permeability coefficient data are weighted and summed, and then the results are sorted in segments according to the numerical range of the calculation results to generate the suitability classification results. Based on the set of weighted coefficient vectors, the indicator sequence data, extraction volume data, and permeability coefficient data are matched sequentially according to unit ID. Then, a weighted multiplication operation is performed on the three elements within the same unit, followed by summation. These are the scaling factors for the three elements. For the normalized values ​​of multiple data types, taking unit A as an example: the normalized index value is 0.905, the normalized extraction volume is 1.000, the normalized permeability coefficient is 0.583, and the scale factors are 0.364, 0.402, and 0.234 respectively. Unit C: If all values ​​are normalized to 0, then... Unit D: Normalized values ​​are 1.000, 0.444, and 1.000; scale factors are assumed to be 0.400, 0.178, and 0.422. After completing the weighted summation of the units, the results are appropriately segmented and sorted according to their numerical ranges. For example, the ranges can be divided into: 0.00–0.30 for level four, 0.30–0.60 for level three, 0.60–0.80 for level two, and 0.80–1.00 for level one. The units are then sorted according to... The corresponding intervals are used to determine the levels, and the results are obtained by sorting them in ascending order from level one to level four to obtain the suitability classification results.

[0048] Table 5: Calculation Results of Suitability Classification

[0049] As shown in Table 5, after weighted summation and interval division, a distribution of suitability levels for multiple units is formed, which can provide a direct ranking basis for the final regional comprehensive development suitability assessment.

[0050] Please see Figure 7 A grading and evaluation system for the suitability of groundwater development in coastal wetlands, including: The groundwater pressure monitoring module obtains groundwater pressure by deploying groundwater pressure monitoring devices at multiple depths and sea surface pressure by deploying sea surface pressure monitoring devices. It performs gradient calculations on adjacent depth monitoring points and secondary gradient calculations on adjacent horizontal monitoring points to generate a pressure field distribution map, which is then transmitted to the seawater intrusion identification module. The seawater intrusion identification module identifies the seawater intrusion propagation path based on the pressure field distribution map, extracts the peak time difference of pressure changes at monitoring points, calculates the pressure transmission delay characteristics, obtains the sensitivity index, and transmits it to the risk assessment and classification module. The risk assessment and grading module compares the sensitivity index with a preset critical threshold. Areas exceeding the threshold are classified as high-risk areas, while areas below the threshold are classified as low-risk areas. Spatial grid grading calculations are performed according to the sensitivity coefficient value range to generate a risk distribution map. The suitability comprehensive evaluation module matches risk zones with permissible mining intensity standards based on the risk zoning map. It calculates development suitability evaluation weight coefficients by the differences in seawater intrusion risk between zones, inputs the weight coefficients into the fuzzy comprehensive evaluation model for comprehensive evaluation calculation, generates suitability evaluation indicators, and transmits them to the development intensity grading module. The development intensity grading module obtains the extraction volume and permeability coefficient based on the suitability evaluation index. It performs weight allocation calculations based on the suitability evaluation index data, recommended extraction volume, and aquifer permeability coefficient, and divides the development intensity level according to the numerical range of the comprehensive evaluation results, generating development suitability grading results.

[0051] 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 method for classifying and evaluating the suitability of groundwater development in coastal wetlands, characterized in that, Includes the following steps: S1: Groundwater pressure is obtained by deploying multi-depth groundwater pressure monitoring devices, sea surface pressure is obtained by deploying sea surface pressure monitoring devices, gradient calculation is performed on adjacent depth monitoring points, secondary gradient calculation is performed on adjacent horizontal monitoring points, and pressure field distribution map is generated. S2: Based on the pressure field distribution map, identify the seawater intrusion propagation path, extract the peak time difference of pressure changes at monitoring points, calculate the pressure transmission delay characteristics, and obtain the sensitivity index; S3: The sensitivity index is compared with a preset critical threshold. If the index exceeds the threshold, it is designated as a high-risk area; if it is below the threshold, it is designated as a safe area. Spatial classification is performed according to the numerical range to generate a risk zoning map. S4: Based on the risk zoning map, match the risk zones with the allowable mining intensity standard, calculate the development suitability evaluation weight coefficient through the difference in seawater intrusion risk between zones, input the weight coefficient into the fuzzy comprehensive evaluation model for comprehensive evaluation calculation, and generate suitability evaluation index; S5: Based on the suitability evaluation index, obtain the extraction volume and permeability coefficient, input the suitability evaluation index, extraction volume and permeability coefficient into the analytic hierarchy process for weight allocation calculation, and generate the suitability classification result.

2. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 1, characterized in that, The pressure field distribution map includes spatial distribution pattern, abnormal change areas and distribution boundary characteristics. The sensitivity index includes pressure transmission coefficient, intrusion path difference and pressure response time difference. The risk zoning map includes level division units, regional boundary range and zoning attribute characteristics. The suitability evaluation index includes resource carrying capacity, environmental stability and utilization feasibility. The suitability classification results include priority mining areas, restricted mining areas and prohibited mining areas.

3. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Using data collected from groundwater pressure monitoring devices and sea surface pressure monitoring devices, index sorting and difference calculation are performed on the pressure values ​​of adjacent depth monitoring points, and the difference results are aggregated in the time series dimension to generate a depth pressure gradient sequence. S102: Call the depth pressure gradient sequence, perform first-order difference processing on the data values ​​of adjacent horizontal monitoring points in the same layer, then perform amplitude calculation on the difference results to obtain the trend value, and combine them into a continuous distribution in the two-dimensional coordinate system to generate the horizontal pressure gradient sequence. S103: Based on the correspondence between the horizontal pressure gradient sequence and the depth pressure gradient sequence, perform coordinate mapping and numerical superposition, and discretize the superposition result in a spatial grid to obtain a pressure field distribution map.

4. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 3, characterized in that, The specific steps of S2 are as follows: S201: Based on the pressure data sequence of the monitoring points in the pressure field distribution map, the curve shape of the monitoring points is retrieved point by point, the local peak amplitude and time index are calibrated, the peak difference between adjacent points is logically judged, and the propulsion path sequence is generated. S202: Call the pressure curves of the monitoring points in the advancement path sequence, extract the time of peak occurrence, compare the peak times of adjacent points, calculate the time difference based on the comparison results, and combine the differences in order into a data sequence to generate a time difference sequence; S203: Based on the numerical distribution of the time difference sequence, the sequence elements are weighted and aggregated, and the result is compared with the delay benchmark value. The normalization coefficient is calculated based on the difference result and the monitoring point is assigned a level to obtain the sensitivity index.

5. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 4, characterized in that, The specific steps for S3 are as follows: S301: Based on the sensitivity index and the preset critical threshold parameter, retrieve the sensitivity index value, perform a difference operation between the value and the critical threshold, if the value is greater than the threshold, mark the corresponding spatial unit as a risk unit, if it is less than the threshold, mark it as a spatial safe unit, and obtain the spatial unit classification identifier. S302: Call the spatial unit classification identifier, perform aggregation processing on the risk unit and safety unit indexes, map the index to the level range according to the sensitivity index distribution interval, sort in ascending order according to the numerical interval, and generate a spatial hierarchical index sequence; S303: Based on the spatial hierarchical index sequence, locate the corresponding spatial unit on the two-dimensional spatial coordinates, perform boundary drawing on the unit range of different numerical intervals, and then fill the interval level index into the corresponding position to generate a risk zoning map.

6. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 5, characterized in that, The preset critical threshold is a numerical limit pre-set according to the range of the sensitivity index. The spatial hierarchical index sequence is a hierarchical index sequence formed by dividing spatial units according to sensitivity index intervals and arranging them in ascending order of values.

7. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 5, characterized in that, The specific steps of S4 are as follows: S401: Based on the distribution data in the risk zoning map, the risk level values ​​are compared with the allowable mining intensity thresholds level by level. During the comparison process, the differences between the level and the threshold are retrieved, and the differences are converted into a set of regional mining intensity correspondences to generate risk intensity matching results. S402: Call the matching parameters in the risk intensity matching result, calculate the weighted difference coefficient between the zoning level value and the seawater intrusion risk difference coefficient, and calculate the weighted coefficient based on the numerical deviation between the difference coefficient and the zoning risk benchmark value to obtain the development suitability evaluation weight coefficient. S403: Call the development suitability evaluation weight coefficient, input the partition weight into the fuzzy comprehensive evaluation model, and perform numerical aggregation operation according to the weight distribution ratio on the multi-factor indicators to obtain the suitability evaluation indicators.

8. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 1, characterized in that, The allowable extraction intensity threshold refers to the limit on the amount of groundwater that can be extracted in a specific area, based on hydrogeological conditions and resource carrying capacity. The development suitability evaluation weight coefficient refers to the zoning evaluation weight coefficient calculated from the zoning level value, risk difference coefficient, and benchmark value deviation.

9. The method for grading and evaluating the suitability of groundwater development in coastal wetlands according to claim 8, characterized in that, The specific steps of S5 are as follows: S501: Based on the suitability evaluation index, compare the index parameter values ​​with the corresponding threshold parameters one by one, record the index values ​​that exceed the threshold parameters, and serialize and integrate the recorded index values ​​according to their categories to generate index sequence data; S502: Call the index sequence data and aggregate it with the mining volume data and permeability coefficient data. Scale the aggregated three types of data within the same numerical range and convert them into a vector structure. Analyze and allocate proportional factors based on the relative differences between the data to obtain a weight coefficient vector. S503: Based on the weight coefficient vector, the corresponding elements of the index sequence data, mining volume data and permeability coefficient data are weighted and summed, and then sorted in segments according to the numerical intervals in which the calculation results are located to generate suitability classification results.

10. A grading and evaluation system for the suitability of groundwater development in coastal wetlands, characterized in that, The system is used to implement the method for grading and evaluating the suitability of groundwater development in coastal wetlands as described in any one of claims 1-9, and the system comprises: The groundwater pressure monitoring module obtains groundwater pressure by deploying groundwater pressure monitoring devices at multiple depths and sea surface pressure by deploying sea surface pressure monitoring devices. It performs gradient calculations on adjacent depth monitoring points and secondary gradient calculations on adjacent horizontal monitoring points to generate a pressure field distribution map, which is then transmitted to the seawater intrusion identification module. The seawater intrusion identification module identifies the seawater intrusion propagation path based on the pressure field distribution map, extracts the peak time difference of pressure changes at monitoring points, calculates the pressure transmission delay characteristics, obtains the sensitivity index, and transmits it to the risk assessment and classification module. The risk assessment and grading module compares the sensitivity index with a preset critical threshold. If the index exceeds the threshold, it is classified as a high-risk area; if it is below the threshold, it is classified as a low-risk area. The module performs spatial grid grading calculations according to the sensitivity coefficient value range to generate a risk distribution map. The suitability comprehensive evaluation module matches the risk zones with the allowable mining intensity standards based on the risk zoning map, calculates the development suitability evaluation weight coefficients through the differences in seawater intrusion risk between zones, inputs the weight coefficients into the fuzzy comprehensive evaluation model for comprehensive evaluation calculation, generates suitability evaluation indicators, and transmits them to the development intensity grading module. The development intensity grading module obtains the extraction volume and permeability coefficient based on the suitability evaluation index, performs weight allocation calculation based on the suitability evaluation index data, recommended extraction volume, and aquifer permeability coefficient, and divides the development intensity level according to the numerical range of the comprehensive evaluation results to generate development suitability grading results.

Images