Underground water vulnerability spatial interpolation method based on deep learning

Abstract

The invention discloses an underground water vulnerability spatial interpolation method based on deep learning, and relates to the technical field of underground water environmental protection and evaluation, and the method comprises the steps: collecting and preprocessing multi-source data; constructing a composite spatial weight matrix based on a spatial topological relation and a distance attenuation principle; a deep learning model is constructed, the weight matrix and the hydrogeological parameters are fused into input features, and a composite loss function is adopted for training; by means of prediction error distribution, weight parameters are dynamically adjusted through Bayesian optimization, the model is retrained, and a vulnerability partition boundary is optimized in combination with morphological operation and filtering processing; according to the method, dynamic calibration and boundary optimization of the spatial weight are realized, spatial geographic constraints and attribute features can be effectively fused, and the spatial consistency, boundary rationality and prediction precision of underground water vulnerability evaluation are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of groundwater environmental protection and assessment technology, specifically a deep learning-based spatial interpolation method for groundwater vulnerability. Background Technology

[0002] Groundwater vulnerability assessment is a key technology for identifying the potential risk of pollution to groundwater systems and serves as an important scientific basis for groundwater resource protection and management. Traditional assessment methods, such as the DRASTIC model, generate vulnerability zoning maps by weighting and overlaying hydrogeological parameters, and have been widely used. However, with increasing demands for assessment accuracy and spatial detail, the limitations of traditional methods, such as fixed weight allocation and coarse boundary characterization, are becoming increasingly apparent.

[0003] In existing technologies, researchers have explored various improvement approaches to enhance the objectivity and accuracy of assessments. For example, some existing techniques utilize the Analytic Hierarchy Process (AHP) to determine the weights of evaluation indicators and combine this with GIS for spatial overlay analysis. While this method introduces a more scientific approach to weight determination, the weights remain fixed across the entire region once determined, failing to reflect spatial heterogeneity. Other existing techniques propose using machine learning algorithms such as random forests to make predictions by learning the relationships between monitoring data and various environmental factors. This method improves the model's nonlinear fitting ability, but it remains essentially an attribute-driven prediction, failing to explicitly incorporate and optimize the geographical dependencies and interactions between spatial units, and lacking effective modeling and post-processing mechanisms for the spatial continuity and smoothness of vulnerable zone boundaries.

[0004] Therefore, existing technologies, whether based on traditional exponential superposition models with fixed weights or currently introduced machine learning prediction methods, all share a common drawback: the failure to achieve adaptive dynamic calibration of spatial weight parameters. This prevents the model from optimizing and adjusting the weight matrix representing spatial relationships during training based on the spatiotemporal distribution characteristics of prediction errors. This deficiency leaves room for improvement in the spatial consistency and boundary rationality of the final vulnerability zoning map, making it difficult to meet the needs of high-precision, refined management of groundwater resources.

[0005] In summary, existing groundwater vulnerability assessment methods still lack effective solutions for how to deeply integrate spatial geographic constraints with machine learning models and construct a closed-loop assessment system that can dynamically calibrate spatial weights and iteratively optimize zoning boundaries. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a deep learning-based spatial interpolation method for groundwater vulnerability. By constructing a dynamically calibrated spatial weight matrix and deeply integrating it into the model training and optimization process, it can achieve adaptive learning of spatial dependencies and refined characterization of vulnerability partition boundaries, thereby overcoming the shortcomings of existing methods in spatial heterogeneity modeling and boundary optimization.

[0007] To solve the above-mentioned technical problems, this invention provides the following technical solution: a deep learning-based spatial interpolation method for groundwater vulnerability, comprising the following steps:

[0008] Step 1: Data Acquisition and Preprocessing: Collect multi-source data of the study area, including water quality monitoring data, geological parameters, land use information and meteorological data. Perform spatial coordinate unification, resolution adjustment, missing value filling and standardization on the collected multi-source data.

[0009] Step 2: Construction of GIS spatial weight matrix: Based on the spatial topological relationship of the study area, adjacency relationships are defined, and distance decay weights are calculated in combination with the distance decay principle. The adjacency relationships and distance decay weights are then integrated to construct a composite spatial weight matrix.

[0010] Step 3, Deep Learning Model Training: Construct a deep learning model based on the U-Net architecture. The composite spatial weight matrix obtained in Step 2 is fused with the multi-source data obtained in Step 1 after preprocessing to form input features. The deep learning model is trained using the input features. The loss function used in the training is a composite loss function, which includes a cross-entropy loss term, a spatial structure similarity loss term, and a spatial autocorrelation loss term.

[0011] Step 4, Dynamic Calibration and Boundary Optimization: Vulnerability prediction is performed using the deep learning model trained in Step 3. The error between the prediction result and the actual monitoring value is calculated. Based on the error distribution, the parameters of the composite spatial weight matrix are dynamically adjusted using the Bayesian optimization method. The input features are re-fused using the adjusted spatial weight matrix, and the deep learning model is retrained. Morphological operations and smoothing filtering are performed on the prediction results of the retrained deep learning model to optimize the vulnerability partition boundary.

[0012] Step 5: Vulnerability partition map generation: Calculate the vulnerability index of each region based on the prediction results of the optimized deep learning model, classify the vulnerability index according to the preset classification threshold, generate a raster format vulnerability partition map, and post-process the raster partition map to convert it into a vector boundary map.

[0013] Furthermore, in step one, the data acquisition and preprocessing further includes:

[0014] The multi-source data is spatially aligned and resampled to a uniform resolution;

[0015] Use interpolation algorithms to fill in missing values ​​in the data;

[0016] Hydrogeological parameters are scored, graded, and standardized.

[0017] The training dataset is augmented by applying random perturbations.

[0018] Furthermore, in step two, constructing the composite spatial weight matrix further includes:

[0019] The adjacency relationship of spatial units is defined based on the Queen adjacency or Rook adjacency rules;

[0020] The distance attenuation weights between spatial units are calculated using the inverse power function of distance or the Gaussian kernel function.

[0021] The adjacency matrix and the distance decay weight matrix are linearly weighted and fused to generate the composite spatial weight matrix;

[0022] Spatial autocorrelation of the generated composite spatial weight matrix is ​​verified by calculating the Moran index.

[0023] Furthermore, in step three, the deep learning model is a U-Net model with an encoder-decoder structure;

[0024] The input features are constructed by concatenating or multiplying each element of the composite spatial weight matrix with the hydrogeological parameter vector of the corresponding spatial unit.

[0025] The hydrogeological parameters include at least the scores for deep groundwater level, net recharge, aquifer medium, soil medium, topography, vadose zone medium, hydraulic conductivity, and land use type, based on the DRASTIC-LU framework.

[0026] Furthermore, the composite loss function The expression is:

[0027]

[0028] in, This represents the cross-entropy loss term, used to measure the difference between the model's predicted classification result and the true label. The spatial structure similarity loss term is used to constrain the spatial structure fidelity of the predicted vulnerability map. This represents the spatial autocorrelation loss term, used to ensure that the prediction results conform to the spatial dependence described by the first law of geography. , , These are the weighting coefficients for each loss, and all are adjustable hyperparameters with a positive value.

[0029] Furthermore, in step four, the error between the calculated prediction result and the actual monitored value further includes:

[0030] Calculate the root mean square error (RMSE) for each spatial cell;

[0031] Generate a heat map of the spatial distribution of prediction errors on a GIS platform;

[0032] A clustering algorithm based on maximum-minimum distance is used to perform spatial clustering analysis on error data to identify the spatial distribution pattern of errors;

[0033] The distance threshold parameter in the clustering algorithm based on maximum-minimum distance is denoted as: .

[0034] Furthermore, in step four, the step of dynamically adjusting the parameters of the composite spatial weight matrix based on the error distribution using a Bayesian optimization method specifically involves:

[0035] The decay coefficient of the distance decay function involved in the construction of the composite spatial weight matrix, and the fusion ratio coefficient of the adjacency weight and the distance weight are used as parameters to be optimized.

[0036] With the objective function of minimizing the overall prediction error of the model on the validation set, the optimal combination of the parameters to be optimized is iteratively searched using a Bayesian optimization framework.

[0037] Furthermore, in step four, the morphological operations and smoothing filtering of the prediction results further include:

[0038] A Gaussian filtering algorithm is used to smooth the vulnerability index raster image to suppress local noise;

[0039] Using opening and closing operations from mathematical morphology to handle vulnerable partition boundaries;

[0040] Set an area threshold to remove vulnerable patches from the partition map whose connected region area is smaller than the threshold.

[0041] A spatial consistency constraint term is introduced into the loss function during the model training phase to penalize unreasonable abrupt changes in vulnerability levels between adjacent spatial units.

[0042] Furthermore, in step four, the dynamic calibration and boundary optimization process is an iterative loop process:

[0043] In each iteration, error analysis, spatial weight parameter adjustment, model retraining, and boundary optimization steps are performed sequentially.

[0044] The iteration termination condition is that the model's performance evaluation metric on the independent test set reaches the preset standard or the number of iterations reaches the upper limit.

[0045] After each iteration, save the updated model and the corresponding spatial weight parameters.

[0046] Furthermore, in step five, generating the raster-formatted vulnerability partition map further includes:

[0047] Based on the output logic value of the optimized model, calculate the vulnerability index DI for each cell;

[0048] The vulnerability index of continuous distribution is divided into high, medium and low vulnerability levels by using the natural breakpoint method, the equal interval method or the quantile method.

[0049] Different levels of vulnerability are assigned colors and rendered to generate vulnerability-level thematic maps;

[0050] The vector boundaries of vulnerable regions at each level are extracted from the raster partition map using a raster-to-vector conversion algorithm.

[0051] Compared with existing technologies, this deep learning-based spatial interpolation method for groundwater vulnerability has the following advantages:

[0052] I. This invention constructs a composite weight matrix that integrates spatial adjacency and distance decay, and uses Bayesian optimization to dynamically calibrate its parameters based on prediction error. This enables the spatial weights to be adaptively adjusted during training, allowing the model to continuously optimize the spatial dependency representation based on error distribution characteristics, thus significantly improving the spatial consistency and boundary rationality of the vulnerability partitioning map.

[0053] Second, this invention integrates a composite spatial weight matrix with multi-source hydrogeological parameters as input features for a deep learning model, and uses a composite loss function including cross-entropy loss, spatial structure similarity loss, and spatial autocorrelation loss for training. This enables the model to learn attribute features and geospatial constraints simultaneously, thereby enhancing the model's ability to model spatial structure fidelity and spatial dependency patterns. While ensuring classification accuracy, it effectively improves the spatial continuity and smoothness of vulnerability prediction results, and enhances the practicality and interpretability of the assessment results.

[0054] Other advantages, objectives and features of the invention will be set forth in part in the description which follows, and in part will be apparent to those skilled in the art from the following examination or study, or may be learned from the practice of the invention. Attached Figure Description

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

[0056] Figure 1 This is a diagram illustrating the method steps of the present invention;

[0057] Figure 2 This is a flowchart of the spatial weight dynamic calibration and model iterative optimization process of the present invention;

[0058] Figure 3 This is a structural diagram of the deep learning model with fused spatial weights according to the present invention. Detailed Implementation

[0059] To further illustrate the technical means and effects of the present invention in achieving its intended purpose, the following detailed description of the specific implementation methods, structures, features, and effects of the present invention, in conjunction with the accompanying drawings and preferred embodiments, is provided below.

[0060] Example

[0061] like Figures 1 to 3 As shown, to clearly present the technical solution, implementation process, and application effects of this invention, the following detailed explanation of the deep learning-based spatial interpolation method for groundwater vulnerability, combined with the actual research area and standardized operating procedures, is provided. This embodiment selects a plain agricultural area as the research area. This area has a gentle terrain, uniform aquifer distribution, and land use types mainly consisting of cultivated land, forest land, and residential areas. Groundwater extraction activities are frequent, making it a typical area requiring accurate groundwater vulnerability assessment.

[0062] This embodiment sequentially executes five core steps: data acquisition and preprocessing, GIS spatial weight matrix construction, deep learning model training, dynamic calibration and boundary optimization, and vulnerability partition map generation, fully disclosing the implementation details of the technical solution.

[0063] The specific implementation steps in this embodiment are as follows:

[0064] Step 1: Data Acquisition and Preprocessing

[0065] The core objective of data acquisition and preprocessing is to obtain high-quality, standardized multi-source data to provide reliable input for subsequent model training and interpolation analysis. The specific implementation is as follows:

[0066] Multi-source data acquisition: Collect multi-source data of the study area. The specific content and acquisition methods of each type of data are as follows:

[0067] Water quality monitoring data: Annual water quality monitoring data from 30 groundwater monitoring wells in the study area were collected, including the concentrations of pollutants such as nitrate, ammonia nitrogen, and total phosphorus. The data came from long-term monitoring records of the local hydrological and environmental monitoring department, with a monitoring frequency of once a month. The annual average value was taken as the valid data.

[0068] Geological parameters: obtained through regional geological survey reports and borehole data, including aquifer media types such as sandstone and gravel, soil media types such as loam, sandy soil and clay, vadose zone thickness and lithology, and hydraulic conductivity coefficient, of which the hydraulic conductivity coefficient was calculated from pumping test data.

[0069] Land use information: The latest Landsat 8 satellite remote sensing imagery of the study area was used to obtain land use type data through supervised classification, which was divided into five categories: cultivated land, forest land, grassland, settlements and water areas, with a spatial resolution of 30 meters.

[0070] Meteorological data: Annual precipitation and evaporation data were collected from three meteorological stations in the study area and surrounding areas. The net replenishment data for the entire region was obtained by spatial interpolation. The net replenishment is equal to precipitation minus evaporation minus surface runoff.

[0071] Spatial Alignment and Resampling: To ensure spatial consistency across various data types, all data were first unified in spatial coordinates using the National Geodetic Coordinate System 2000 and the Gauss-Kruger projection with 3-degree zoning. All data were then resampled to a uniform resolution; in this embodiment, a 30m × 30m raster resolution was selected. The resampling method employed was bilinear interpolation, which effectively preserves the spatial continuity of the data and is suitable for adjusting the resolution of continuous data.

[0072] Missing value imputation: Due to occasional missing data at some monitoring points, ordinary kriging interpolation is used to impute the missing values. Ordinary kriging interpolation is based on the principle of spatial autocorrelation and assumes that the spatial distribution of the data conforms to second-order stationarity. Its calculation formula is as follows:

[0073]

[0074] in These are estimates of the missing points. These are the measured values ​​from known monitoring points. For the weighting coefficients, satisfying The weighting coefficients are obtained by solving the variogram equation. In this embodiment, a spherical variogram model is selected, with a nugget value of 0.02, a sill value of 0.15, and a range of 5000 meters. Cross-validation is used to determine the rationality of the model parameters, ensuring that the data error after padding is within an acceptable range. The nugget value reflects measurement error and micro-variation, the sill value reflects the total degree of data variation, and the range represents the effective range of spatial autocorrelation.

[0075] Feature engineering processing: Scoring, grading, and standardization of hydrogeological parameters:

[0076] Scoring and Grading: Based on the DRASTIC-LU framework, eight parameters were scored: depth, groundwater level, net recharge, aquifer medium, soil medium, topography, vadose zone medium, hydraulic conductivity coefficient, and land use type. The scoring criteria referenced relevant industry standards: deeper groundwater levels indicated lower vulnerability (scoring range 1-10); larger net recharge indicated higher vulnerability (scoring range 1-10); gravel in the aquifer medium received higher scores than sandstone; sandy soil received higher scores than loam in the soil medium; gentler topographic slopes resulted in higher scores; poorer permeability in the vadose zone medium resulted in lower scores; higher hydraulic conductivity coefficients resulted in higher scores; and cultivated land received higher scores than forest land and residential areas in land use type. Specific scoring results were fine-tuned based on the actual hydrogeological conditions of the study area to ensure that the scores reflected the degree of influence of the parameters on groundwater vulnerability.

[0077] Standardization: To eliminate differences in the units and ranges of different parameters, the min-max standardization method is used to process all scored parameters. The formula is as follows:

[0078]

[0079] in The values ​​are standardized. The original score. The minimum score for this parameter. This represents the maximum score value for this parameter. The standardized data ranges from 0 to 1, facilitating feature fusion and training of the model.

[0080] Data augmentation: Data augmentation techniques expand the training dataset by applying random perturbations to improve the model's generalization ability. Specifically, this involves adding a normally distributed random perturbation term to each parameter of the standardized raster data. The perturbation magnitude is controlled within 5% of the original data value, meaning the perturbed data becomes... ,in The dataset follows a normal distribution with a mean of 0 and a variance of 0.05². Simultaneously, the data is horizontally and vertically flipped to generate new sample data. Through data augmentation, the number of samples in the original dataset is increased to three times, ensuring sufficient data support for model training.

[0081] Step 2: Constructing the GIS Spatial Weight Matrix

[0082] The GIS spatial weight matrix is ​​used to quantify the interdependencies between spatial units within the research area, providing a foundation for incorporating spatial geographic constraints into the model. The specific implementation of this embodiment is as follows:

[0083] Spatial cell division: First, the study area was divided into several spatial cells with a grid resolution of 30 meters × 30 meters. Each spatial cell corresponds to one grid pixel, and the entire study area was divided into approximately 556,000 spatial cells. Each spatial cell serves as an independent evaluation unit, storing the multi-source data characteristics and vulnerability assessment results of that unit.

[0084] Adjacency Definition: Adjacency relationships between spatial units are defined based on the Queen adjacency rule. The Queen adjacency rule states that two spatial units are considered adjacent if they share an edge or a vertex. This embodiment chooses the Queen adjacency rule because it can more comprehensively reflect the geographical relationships between spatial units, and is particularly suitable for study areas with flat terrain and strong spatial continuity. Adjacency relationships are implemented using the spatial analysis tools in ArcGIS software, generating an adjacency matrix. ,in 1 represents a unit and unit Adjacency, 0 indicates non-adjacency.

[0085] Distance attenuation weight calculation: The distance attenuation weight between spatial units is calculated using the inverse power function of distance. The inverse power function of distance can better reflect the basic geographical law that the closer the distance, the greater the mutual influence. The calculation formula is as follows:

[0086]

[0087] in spatial unit and unit Distance decay weight between them Let be the Euclidean distance between the center points of the two spatial units. This is the attenuation coefficient, used to control the degree to which distance affects the weighting. In this embodiment... The value of is 2, which was determined through multiple experiments. This value can highlight the influence of nearby units while avoiding excessive concentration of weights on adjacent units.

[0088] In some alternative implementations, a Gaussian kernel function can also be used to calculate the distance attenuation weight, and the calculation formula is as follows:

[0089]

[0090] in This is a bandwidth parameter used to adjust the rate of distance attenuation. In this embodiment, if this function is used... The value was set at 100 meters, and its rationality was determined through spatial autocorrelation testing.

[0091] Construction of composite weight matrix: The adjacency relation matrix and the distance decay weight matrix are linearly weighted and fused to generate a composite spatial weight matrix. The fusion formula is:

[0092]

[0093] in This is the fusion ratio coefficient of adjacency weight and distance weight, with a value ranging from 0 to 1. In this embodiment... The initial value is 0.4, and it will be dynamically adjusted through Bayesian optimization. After fusion, the weights of each spatial unit are normalized to ensure that the sum of the weights corresponding to each unit is 1. Normalization can prevent excessively large or small weight values ​​from affecting model training.

[0094] Weight Matrix Validation: The spatial autocorrelation of the generated composite spatial weight matrix is ​​validated by calculating the Moran index. The Moran index is a commonly used indicator for measuring the correlation of spatial data, and its calculation formula is as follows:

[0095]

[0096] in This represents the total number of spatial units. and They are spatial units and The standardized value of a certain hydrogeological parameter, This is the average value of this parameter across all spatial units. The Moran index ranges from -1 to 1. When the data shows a positive spatial correlation, it means that units with similar attributes are clustered in space, indicating that the weight matrix can effectively characterize spatial dependence; when The data is randomly distributed; when The data exhibits a spatial negative correlation. In this embodiment, the Moran index is calculated to be 0.68, indicating that the composite spatial weight matrix has a good spatial positive correlation and can meet the model's spatial constraint requirements.

[0097] Step 3: Deep learning model training:

[0098] Deep learning model training is the core of achieving accurate interpolation of groundwater vulnerability. By fusing spatial weights and features from multiple data sources, a predictive model capable of capturing spatial dependence and nonlinear relationships is constructed. The specific implementation is as follows:

[0099] Model Architecture Design: A deep learning model based on the U-Net architecture is constructed. The U-Net model adopts an encoder-decoder structure, which has strong feature extraction and spatial reconstruction capabilities, and is suitable for spatial interpolation tasks on raster data. The specific structure of the model is as follows:

[0100] The encoder section contains four downsampling blocks, each consisting of two convolutional layers, one batch normalization layer, and one activation function. The convolutional layers have a 3×3 kernel size, a stride of 1, and SAME padding to ensure the feature map size remains unchanged after convolution. The batch normalization layer accelerates model convergence and prevents gradient vanishing. The activation function is ReLU, expressed as follows: This allows for the introduction of non-linear features. A max-pooling layer is set after each downsampling block, with a pooling kernel size of 2×2 and a stride of 2, which halves the size of the feature map and increases the receptive field of the model.

[0101] The decoder section consists of four upsampling blocks. Each upsampling block uses a transposed convolutional layer to double the size of the feature map. The transposed convolutional kernel size is 2×2, and the stride is 2. After upsampling, the corresponding encoder feature map is concatenated with the current decoder feature map through skip connections to supplement low-level spatial details. The concatenated feature map is then processed by two convolutional batch normalization layers and a ReLU activation function to further extract fused features.

[0102] Output layer: A 1×1 convolutional layer is used to map the feature map output by the decoder into a feature vector corresponding to the number of input spatial units. Then, the Softmax activation function is used to convert the feature vector into a probability distribution of each vulnerability level. The number of channels in the output layer is consistent with the number of vulnerability levels. In this embodiment, there are 3 categories, namely high, low and medium vulnerability.

[0103] Input feature construction: The composite spatial weight matrix is ​​fused with the DRASTIC-LU parameters to form the model's input features. Specifically, element-wise multiplication is used; that is, for each spatial unit i, its input feature vector... Normalize the DRASTIC-LU parameter vector of this unit. The eight parameters correspond to eight dimensions, which are the weight vectors corresponding to the units in the composite space weight matrix. Perform element-wise multiplication to obtain , This represents element-wise multiplication. In this way, the interdependencies between spatial units are incorporated into the feature representation of each unit, enabling the model to consider both its own attribute features and spatial correlation features during training.

[0104] In some alternative implementations, input features can also be constructed using a concatenation method, i.e. By concatenating the parameter vector and the weight vector into a 16-dimensional feature vector, this method can also achieve feature fusion and is suitable for scenarios that place equal emphasis on spatial correlation and attribute features.

[0105] Dataset partitioning: The preprocessed dataset is divided into training, validation, and test sets in a 7:2:1 ratio. The training set is used for learning model parameters, the validation set is used to monitor overfitting during model training and adjust hyperparameters, and the test set is used to finally evaluate the model's generalization performance. Stratified sampling is used during partitioning to ensure that the proportion of samples at each vulnerability level in the training, validation, and test sets is consistent with the original dataset, avoiding model training bias due to uneven sample distribution.

[0106] Model training process: The model is trained using supervised learning, and the specific process is as follows:

[0107] Initialize model parameters: The weights of all convolutional and transposed convolutional layers are initialized using a He normal distribution, and the bias parameters are initialized to 0.

[0108] Forward propagation: The input features of the training set are fed into the model, and after feature extraction and fusion by the encoder and decoder, the vulnerability level probability distribution of each sample is output.

[0109] Loss Calculation: The error between the model's prediction and the true label is calculated using a composite loss function. The expression for the composite loss function is:

[0110]

[0111] in The cross-entropy loss term measures the difference between the model's predicted classification result and the true label. The calculation formula is as follows: , The one-hot encoded vector of the true label, if the sample belongs to the first... Class Otherwise, it is 0. The model predicts the first Class probability, In this embodiment, the number of vulnerability levels is [number]. .

[0112] The spatial structure similarity loss term is used to constrain the spatial structure fidelity of the predicted vulnerability map, and its calculation formula is as follows: ,in This is a structural similarity index, with values ​​ranging from 0 to 1. For the true vulnerability map, To predict vulnerability maps. The calculation considers three dimensions: brightness, contrast, and structure. The formula is:

[0113]

[0114] and They are respectively and The mean, and They are variance, For covariance, and The constant is used to avoid a denominator of 0. In this embodiment... , .

[0115] The spatial autocorrelation loss term is used to ensure that the prediction results conform to the spatial dependence described by the first law of geography. The calculation formula is as follows: , To predict the Moran index of the vulnerability map, the calculation method is the same as the Moran index formula in step 2, point 5.

[0116] , , The weighting coefficients for each loss are all adjustable hyperparameters with a positive value, used to balance the importance of each loss. In this embodiment, they are determined through grid search. , , This combination can ensure classification accuracy while taking into account spatial structure and spatial autocorrelation.

[0117] Backpropagation: An adaptive moment estimation optimizer is used to update the model parameters based on the gradient of the loss function. The learning rate of the Adam optimizer is initially set to 0.001. As the training epochs increase, a cosine annealing strategy is used to adjust the learning rate to avoid the model getting trapped in local optima. During training, the batch size is set to 32 and the number of training epochs is set to 100. Training stops when the loss function value on the validation set no longer decreases for 10 consecutive epochs to avoid overfitting.

[0118] Model Evaluation: Multiple evaluation metrics were used to assess model performance, including overall classification accuracy, kappa coefficient, mean intersection-over-union ratio (Intersection over Union), and RMSE. Overall classification accuracy is the proportion of correctly classified samples out of the total number of samples. The kappa coefficient measures the deviation of the classification result from the random classification, ranging from -1 to 1; a value closer to 1 indicates better classification performance. The mean intersection-over-union ratio is the average of the intersection ratios of each category, and it is the ratio of the area of ​​intersection between the predicted and true categories to the area of ​​union. RMSE measures the root mean square error between the predicted vulnerability index and the true index; a smaller value indicates higher interpolation accuracy. In this embodiment, the trained model achieved an overall classification accuracy of 89.2%, a kappa coefficient of 0.83, a mean intersection-over-union ratio of 0.81, and an RMSE of 0.06 on the test set, indicating good predictive performance.

[0119] Step 4: Dynamic calibration and boundary optimization

[0120] Dynamic calibration and boundary optimization improve the model's prediction accuracy and spatial consistency by iteratively adjusting spatial weight parameters and optimizing partition boundaries. The specific implementation is as follows:

[0121] Error Analysis: The trained deep learning model is used to predict vulnerabilities on the test set, and the error between the prediction results and the actual monitored values ​​is calculated. The specific analysis process is as follows:

[0122] Calculate RMSE: For each spatial unit, calculate the RMSE between its predicted vulnerability index and the actual monitored value. The calculation formula is as follows:

[0123]

[0124] in The number of times this unit is monitored, in this embodiment That is, monthly monitoring data, For the first The actual value of the monitoring. For the first The predicted value for the next time.

[0125] Generating an error heatmap: In the ArcGIS platform, the RMSE value of each spatial unit is used as an attribute value to generate a spatial distribution heatmap of the prediction error. The heatmap uses a hierarchical color scheme, dividing the RMSE values ​​into three levels: low error (0 to 0.05), medium error (0.05 to 0.10), and high error (greater than 0.10), corresponding to green, yellow, and red respectively. The heatmap allows for intuitive location of high-error areas, providing targeted targets for parameter adjustments.

[0126] Error Spatial Clustering Analysis: A maximum-minimum distance-based clustering algorithm is used to perform spatial clustering analysis on error data to identify spatial distribution patterns of errors. The specific steps of the algorithm are as follows: randomly select a high-error spatial unit as the initial cluster center; calculate the minimum distance from each of the remaining spatial units to all selected cluster centers; if this minimum distance is greater than a distance decision threshold parameter... If a new cluster center is selected, that unit is chosen as the new cluster center; repeat the above two steps until no new cluster centers are generated; then assign each spatial unit to the cluster containing the nearest cluster center. In this embodiment... The value was set to 0.07, which is determined based on twice the standard deviation of the error data (0.035), and can effectively distinguish different error clustering areas. Through cluster analysis, three main high error clustering areas were identified, located in the northeastern concentrated farmland area, the southwestern densely populated area, and the central area with high groundwater extraction intensity, respectively.

[0127] Spatial error statistics: For each cluster, i.e., the error region, the mean, variance, and standard deviation of the error are calculated. For example, the mean error of the high-error cluster in the northeast is 0.12, the variance is 0.003, and the standard deviation is 0.055; the mean error of the cluster in the southwest is 0.11, the variance is 0.0028, and the standard deviation is 0.053; and the mean error of the cluster in the central region is 0.13, the variance is 0.0032, and the standard deviation is 0.057. These statistical results provide data support for the subsequent targeted adjustment of spatial weight parameters.

[0128] Spatial weight parameter adjustment: A Bayesian optimization method is used to dynamically adjust the parameters of the composite spatial weight matrix based on the error distribution. The specific implementation is as follows:

[0129] Determine the parameters to be optimized: the attenuation coefficients of the distance attenuation function involved in the construction of the composite spatial weight matrix. That is, the fusion ratio coefficient of the inverse power function of distance, adjacency weight and distance weight. As parameters to be optimized. The search space for the parameters to be optimized is set to... Between 1 and 3, Between 0.2 and 0.8.

[0130] Define the objective function: Minimize the overall prediction error of the model on the validation set, i.e., the objective function is... ,in This represents the average RMSE of the model on the validation set.

[0131] Bayesian Optimization Search: Bayesian optimization constructs a probabilistic model of the objective function based on Gaussian process regression. In this embodiment, the expected improvement (EI) function is used to select the next parameter combination to be evaluated, iteratively searching for the optimal parameters. The specific steps are as follows: Initialization: Randomly select 5 parameter combinations within the search space and calculate the objective function value corresponding to each parameter combination; Construct a Gaussian process model with the parameter combination as input and the objective function value as output, and train the Gaussian process model; Select the next set of parameters and calculate the expected improvement value of each candidate parameter combination using the EI function, selecting the parameter combination with the largest expected improvement value; Evaluate the objective function: Substitute this parameter combination into the model and calculate the objective function value; Update the model: Add the new parameter combination and objective function value to the training set and update the Gaussian process model; Repeat the above three steps until the number of iterations reaches a preset upper limit (20 in this embodiment) or the objective function value converges, and the difference in the objective function value between two adjacent iterations is less than 0.001. Through Bayesian optimization, the final optimal parameter combination is determined as follows: , .

[0132] Model retraining: Reconstruct the composite spatial weight matrix using the adjusted spatial weight parameters, then re-integrate the input features, and retrain the deep learning model according to the model architecture and training parameters from step three. The specific implementation is as follows:

[0133] Reconstructing input features: based on optimized and The distance decay weight and composite spatial weight matrix are recalculated and then multiplied element-wise with the DRASTIC-LU parameter normalization vector to construct new input features.

[0134] Model retraining: The model architecture, batch size, training epochs, and optimizer are the same as in the initial training, with only the input features changed. During training, a custom callback function updates the spatial weight parameters in real time, achieving co-optimization of the parameters and the model. After each training epoch, the custom callback function calculates the error on the validation set. If the error decreases, the current spatial weight parameters and model parameters are saved; if the error increases, the parameters from the previous epoch are restored. The model is saved after each iteration for subsequent analysis and validation.

[0135] Boundary optimization: Morphological operations and smoothing filters are performed on the prediction results of the retrained model to optimize the boundaries of vulnerable partitions. The specific implementation is as follows:

[0136] Gaussian filtering smoothing: A Gaussian filtering algorithm is used to smooth the vulnerability index raster image to suppress local noise. The convolution kernel size of the Gaussian filter is 3×3, and the standard deviation of the Gaussian function is [not specified]. The value is set to 1.0, and each element of the convolution kernel is processed by a Gaussian function:

[0137]

[0138] The calculated values ​​for the 3×3 Gaussian convolution kernel are as follows:

[0139]

[0140] By using convolution operations, the vulnerability index value of each spatial unit is influenced by its neighboring units, reducing local fluctuations and making the index distribution smoother.

[0141] Morphological operations: Opening and closing operations from mathematical morphology are used to process the boundaries of vulnerable partitions. The opening operation first performs erosion followed by dilation to eliminate small white noise, i.e., small, highly vulnerable patches. The calculation formula is as follows: ,in This is a vulnerability partitioning map. As a structural element, this embodiment uses a 3×3 square structural element. This indicates a corrosion operation. This indicates an expansion operation. The closing operation performs expansion before erosion, used to fill small black voids, i.e., small, low-vulnerability patches. The calculation formula is as follows: By using opening and closing operations, the boundaries of vulnerable partitions are made more regular.

[0142] Small Area Removal: An area threshold is set to remove vulnerable patches with connected regions smaller than the threshold from the zoning map. In this embodiment, the area threshold is set to 900 square meters, which is the area of ​​10 30m × 30m grid cells. This threshold is determined based on the actual situation of the study area and the readability of the zoning map, and can effectively remove small patches with no practical significance while retaining important small-scale highly vulnerable areas.

[0143] Spatial consistency constraint: A spatial consistency constraint term is introduced into the loss function during model training to penalize unreasonable abrupt changes in vulnerability levels between adjacent spatial units. The formula for calculating the spatial consistency constraint term is:

[0144]

[0145] in spatial unit The set of adjacent units, The average number of adjacent cells per cell. and Units and The vulnerability level is represented by values ​​1, 2, and 3, indicating low, medium, and high vulnerability. This constraint is added to the composite loss function, resulting in the new loss function:

[0146]

[0147] in The constraint weight is set to 0.1. This constraint aims to minimize the difference in vulnerability levels between adjacent cells, thereby improving the spatial consistency of the partition map.

[0148] Iterative optimization: The dynamic calibration and boundary optimization process is an iterative loop process, and the specific iterative process is as follows:

[0149] First iteration: using initial spatial weight parameters , Train the model, perform error analysis, and obtain the optimized parameters through Bayesian optimization. , The model was retrained and the boundaries were optimized. The model performance was evaluated. The RMSE on the test set was 0.048, and the overall classification accuracy was 91.5%.

[0150] Second iteration: Based on the model prediction results of the first iteration, error analysis is performed, and Bayesian optimization is repeated to obtain new parameters. The model was retrained and the boundaries were optimized. The performance was evaluated, and the RMSE on the test set was 0.042, with an overall classification accuracy of 93.2%.

[0151] Third iteration: Repeat the above process to obtain the parameters. , The model's RMSE on the test set was 0.041, and the overall classification accuracy was 93.5%.

[0152] Iteration Termination: In the fourth iteration, the RMSE on the test set is 0.040, the difference from the third iteration is less than 0.001, and the overall classification accuracy reaches 93.6%, meeting the preset performance criteria: RMSE ≤ 0.05 and overall classification accuracy ≥ 90%. Therefore, the iteration terminates. After each iteration, the updated model and corresponding spatial weight parameters are saved, and the model from the fourth iteration is finally selected as the optimal model.

[0153] Step 5: Vulnerability Partition Map Generation

[0154] Vulnerability partitioning map generation transforms the optimized model predictions into visual and applicable evaluation results. The specific implementation is as follows:

[0155] Vulnerability Index Calculation: Based on the prediction results of the optimized deep learning model, the vulnerability index DI for each spatial unit is calculated. The model output is the probability distribution of each vulnerability level. These correspond to the probabilities of low, medium, and high vulnerability, respectively. The vulnerability index DI is calculated using a weighted summation method, and the formula is as follows: 1, 2, and 3 represent the weights for low, medium, and high vulnerability levels, respectively. The value of DI ranges from 1 to 3. The closer the value is to 1, the lower the vulnerability, and the closer it is to 3, the higher the vulnerability.

[0156] Vulnerability Classification: The natural break method is used to classify the continuously distributed vulnerability index into three vulnerability levels: high, low, and medium. Based on the natural distribution characteristics of the data, the natural break method classifies levels by finding natural breaks in the data values, minimizing the data differences within each level and maximizing the differences between levels. Specifically, this is implemented using the classification tool in ArcGIS software. Inputting DI raster data, selecting the natural break method, and setting the number of classifications to 3, the software automatically calculates two breakpoint values: 1.8 and 2.5. Therefore, the vulnerability classification standards are: low vulnerability DI between 1 and 1.8, medium vulnerability DI between 1.8 and 2.5, and high vulnerability DI between 2.5 and 3.

[0157] In some alternative implementations, the equidistant method or the quantile method can also be used for classification. The equidistant method divides the DI value range of 1 to 3 into three equal intervals, each interval having a length of 0.667. The classification criteria are low vulnerability (1 to 1.667), medium vulnerability (1.667 to 2.333), and high vulnerability (2.333 to 3). The quantile method divides the data according to the quantiles, so that the number of samples in each level is approximately equal. In this embodiment, if the quantile method is used, the classification criteria are low vulnerability (1 to 1.7), medium vulnerability (1.7 to 2.4), and high vulnerability (2.4 to 3).

[0158] Raster zoning map generation: Based on the vulnerability classification results, a raster-format vulnerability zoning map is generated. In ArcGIS software, the DI raster data is reclassified according to the classification criteria, assigning different numerical codes to different DI values: low vulnerability is 1, medium vulnerability is 2, and high vulnerability is 3. Then, different colors are assigned to different levels: low vulnerability is green, medium vulnerability is yellow, and high vulnerability is red. Color rendering is then performed to generate a vulnerability classification thematic map. The resolution of the raster zoning map is consistent with the spatial cell resolution, at 30 meters × 30 meters, which can clearly present the spatial distribution characteristics of groundwater vulnerability in the study area.

[0159] Post-processing and Vector Boundary Extraction: Post-processing of the raster partition map mainly includes removing isolated small patches to ensure consistency with the small area removal in step 4, point 4, and smoothing the boundaries to ensure the standardization and readability of the partition map. Then, the vector boundaries of vulnerable areas at each level are extracted from the raster partition map using a raster-to-vector algorithm. In this embodiment, the Douglas-Pock algorithm is used for vector boundary extraction. The core idea of ​​this algorithm is to simplify the boundary by reducing the number of points on the vector boundary while maintaining the basic shape of the boundary. The specific steps are as follows: connect the start and end points of the vector boundary to form a straight line; calculate the distance from each intermediate point on the boundary to the straight line; if the maximum distance is greater than a preset threshold (5 meters in this embodiment), the point is retained, and the boundary is divided into two segments using this point as the boundary, and the above steps are repeated; if the maximum distance is less than or equal to the threshold, all intermediate points are deleted, and the original boundary segment is replaced with a straight line. The vector boundaries extracted by the Douglas-Pock algorithm are concise and smooth, which is convenient for subsequent spatial analysis and groundwater resource protection planning applications.

[0160] In some optional embodiments, the deep learning model architecture and spatial weight parameter optimization method are replaced to verify the flexibility and feasibility of the technical solution of the present invention.

[0161] Model architecture replacement: Replace the U-Net model with the ResU-Net model. The ResU-Net model introduces residual connections on top of U-Net, effectively solving the gradient vanishing problem in deep networks and improving the model's feature extraction capabilities. The formula for calculating residual connections is as follows: ,in For the residual function, Using the input features as inputs, the model can directly learn the residuals between the input and output features through residual connections, thus accelerating the training convergence speed.

[0162] Optimization method replacement: The Bayesian optimization method is replaced with a grid search method. The grid search method selects several parameter combinations evenly within the search space of the parameters to be optimized, evaluates the model performance of each combination, and selects the optimal parameter combination. In this replacement example, the parameters... The search space is 1 to 3, and five candidate values ​​are selected: 1, 1.5, 2, 2.5, and 3; parameters The search space is from 0.2 to 0.8. Seven candidate values ​​of 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 and 0.8 are selected, forming a total of 35 parameter combinations. The model is trained and its performance is evaluated one by one.

[0163] After the replaced model was trained, the overall classification accuracy of the test set was 92.1% and the RMSE was 0.045, which is close to the overall classification accuracy of 93.6% and the RMSE of 0.040 of the original embodiment. This shows that even if the model architecture and optimization method are replaced, the technical solution of the present invention can still achieve high prediction accuracy, verifying the flexibility and reliability of the solution.

[0164] Through the detailed implementation of this embodiment, a deep learning-based spatial interpolation method for groundwater vulnerability was successfully constructed, generating a high-precision groundwater vulnerability zoning map of the study area. Compared with the traditional DRASTIC model, this method improves both prediction accuracy and spatial consistency, produces smoother and more regular vulnerability zoning boundaries, and accurately identifies the spatial distribution of highly vulnerable areas. It provides reliable technical support for groundwater resource protection and management and has already been applied in the "Yunnan Province Groundwater Pollution Prevention and Control Key Zone Delineation Scheme (Trial)" with good results.

[0165] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A deep learning-based spatial interpolation method for groundwater vulnerability, characterized in that, Includes the following steps: Step 1: Data Acquisition and Preprocessing: Collect multi-source data of the study area, including water quality monitoring data, geological parameters, land use information and meteorological data. Perform spatial coordinate unification, resolution adjustment, missing value filling and standardization on the collected multi-source data. Step 2: Construction of GIS spatial weight matrix: Based on the spatial topological relationship of the study area, adjacency relationships are defined, and distance decay weights are calculated in combination with the distance decay principle. The adjacency relationships and distance decay weights are then integrated to construct a composite spatial weight matrix. Step 3, Deep Learning Model Training: Construct a deep learning model based on the U-Net architecture. The composite spatial weight matrix obtained in Step 2 is fused with the multi-source data obtained in Step 1 after preprocessing to form input features. The deep learning model is trained using the input features. The loss function used in the training is a composite loss function, which includes a cross-entropy loss term, a spatial structure similarity loss term, and a spatial autocorrelation loss term. Step 4, Dynamic Calibration and Boundary Optimization: Vulnerability prediction is performed using the deep learning model trained in Step 3. The error between the prediction result and the actual monitoring value is calculated. Based on the error distribution, the parameters of the composite spatial weight matrix are dynamically adjusted using the Bayesian optimization method. The input features are re-fused using the adjusted spatial weight matrix, and the deep learning model is retrained. Morphological operations and smoothing filtering are performed on the prediction results of the retrained deep learning model to optimize the vulnerability partition boundary. Step 5: Vulnerability partition map generation: Calculate the vulnerability index of each region based on the prediction results of the optimized deep learning model, classify the vulnerability index according to the preset classification threshold, generate a raster format vulnerability partition map, and post-process the raster partition map to convert it into a vector boundary map.

2. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, In step one, the data acquisition and preprocessing further includes: The multi-source data is spatially aligned and resampled to a uniform resolution; Use interpolation algorithms to fill in missing values ​​in the data; Hydrogeological parameters are scored, graded, and standardized. The training dataset is augmented by applying random perturbations.

3. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, In step two, constructing the composite spatial weight matrix further includes: The adjacency relationship of spatial units is defined based on the Queen adjacency or Rook adjacency rules; The distance attenuation weights between spatial units are calculated using the inverse power function of distance or the Gaussian kernel function. The adjacency matrix and the distance decay weight matrix are linearly weighted and fused to generate the composite spatial weight matrix; Spatial autocorrelation of the generated composite spatial weight matrix is ​​verified by calculating the Moran index.

4. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, In step three, the deep learning model is a U-Net model with an encoder-decoder structure; The input features are constructed by concatenating or multiplying each element of the composite spatial weight matrix with the hydrogeological parameter vector of the corresponding spatial unit. The hydrogeological parameters include at least the scores for deep groundwater level, net recharge, aquifer medium, soil medium, topography, vadose zone medium, hydraulic conductivity, and land use type, based on the DRASTIC-LU framework.

5. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, The composite loss function The expression is: in, This represents the cross-entropy loss term, used to measure the difference between the model's predicted classification result and the true label. The spatial structure similarity loss term is used to constrain the spatial structure fidelity of the predicted vulnerability map. This represents the spatial autocorrelation loss term, used to ensure that the prediction results conform to the spatial dependence described by the first law of geography. , , These are the weighting coefficients for each loss, and all are adjustable hyperparameters with a positive value.

6. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, In step four, the error between the calculated prediction result and the actual monitored value further includes: Calculate the root mean square error (RMSE) for each spatial cell; Generate a heat map of the spatial distribution of prediction errors on a GIS platform; A clustering algorithm based on maximum-minimum distance is used to perform spatial clustering analysis on error data to identify the spatial distribution pattern of errors; The distance threshold parameter in the clustering algorithm based on maximum-minimum distance is denoted as: .

7. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, In step four, the step of dynamically adjusting the parameters of the composite space weight matrix based on the error distribution using a Bayesian optimization method specifically involves: The decay coefficient of the distance decay function involved in the construction of the composite spatial weight matrix, and the fusion ratio coefficient of the adjacency weight and the distance weight are used as parameters to be optimized. With the objective function of minimizing the overall prediction error of the model on the validation set, the optimal combination of the parameters to be optimized is iteratively searched using a Bayesian optimization framework.

8. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, Step four involves performing morphological operations and smoothing filtering on the prediction results, further including: A Gaussian filtering algorithm is used to smooth the vulnerability index raster image to suppress local noise; Using opening and closing operations from mathematical morphology to handle vulnerable partition boundaries; Set an area threshold to remove vulnerable patches from the partition map whose connected region area is smaller than the threshold. A spatial consistency constraint term is introduced into the loss function during the model training phase to penalize unreasonable abrupt changes in vulnerability levels between adjacent spatial units.

9. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, In step four, the dynamic calibration and boundary optimization process is an iterative loop process: In each iteration, error analysis, spatial weight parameter adjustment, model retraining, and boundary optimization steps are performed sequentially. The iteration termination condition is that the model's performance evaluation metric on the independent test set reaches the preset standard or the number of iterations reaches the upper limit. After each iteration, save the updated model and the corresponding spatial weight parameters.

10. The deep learning-based spatial interpolation method for groundwater vulnerability according to claim 1, characterized in that, Step five, generating a raster-formatted vulnerability partition map, further includes: Based on the output logic value of the optimized model, calculate the vulnerability index DI for each cell; The vulnerability index of continuous distribution is divided into high, medium and low vulnerability levels by using the natural breakpoint method, the equal interval method or the quantile method. Different levels of vulnerability are assigned colors and rendered to generate vulnerability-level thematic maps; The vector boundaries of vulnerable regions at each level are extracted from the raster partition map using a raster-to-vector conversion algorithm.

Images