A Dynamic Groundwater Level Assessment Method Based on Multi-Source Data Mapping and Fluid Model Correction

Abstract

This invention relates to the field of groundwater status assessment technology, and discloses a method for dynamic assessment of groundwater levels based on multi-source data mapping and fluid model correction. This invention aims to address the technical problems of low spatial resolution in existing satellite gravity sounding data, and the difficulty in accurately separating deep groundwater change signals from total land water storage signals, resulting in insufficient accuracy and a lack of spatial detail in groundwater level monitoring. This invention unifies the spatiotemporal reference of multi-source heterogeneous data, utilizes high-resolution surface texture features to guide the spatial downscaling reconstruction of the low-resolution gravity field; combines spectral features and the principle of thermal inertia to progressively eliminate surface water and soil water storage signals; introduces measured data to construct a spatial error correction surface for water yield parameters; and uses image connectivity analysis technology to identify groundwater funnels.

Description

Technical Field

[0001] This invention relates to the field of groundwater condition assessment technology, and more specifically, to a method for dynamic assessment of groundwater levels based on multi-source data mapping and fluid model correction. Background Technology

[0002] Groundwater, as a critical freshwater resource, directly impacts regional ecological security and the sustainability of agricultural irrigation due to its dynamic changes in water levels. Currently, monitoring of regional groundwater levels relies primarily on two methods: traditional ground-based observation well networks and remote sensing inversion based on satellite gravity sounding technology (such as the GRACE satellite and its subsequent missions). Ground-based observation wells can provide high-precision data for single points, but due to construction and maintenance costs, the distribution of these sites is often sparse and uneven, making it difficult to cover large areas. In contrast, satellite gravity sounding technology can provide time-varying gravity field data covering the globe. The overall changes in terrestrial water storage derived from this inversion have become an important tool for macro-level water resource assessment.

[0003] However, in practical engineering applications, existing monitoring methods face a substantial technical problem: the spatial resolution of satellite gravity data is generally low (typically on the order of 300 kilometers), and the signals detected are a comprehensive superposition of changes in surface water, soil water, and groundwater quality, making direct differentiation difficult. Existing processing methods often struggle to effectively identify the spatial distribution details of groundwater at small scales (such as county or irrigation district scales) while maintaining accuracy in macroscopic total data. They also struggle to accurately extract deep groundwater change signals from complex surface hydrological noise, resulting in groundwater level monitoring results that are often overly smooth or contain systematic biases, failing to meet the practical needs for refined identification and risk assessment of groundwater over-extraction funnel areas. Summary of the Invention

[0004] In view of the above-mentioned existing problems, the present invention aims to solve the technical problems of low resolution of existing satellite gravity data and difficulty in separating deep groundwater signals from the total land water storage, resulting in insufficient monitoring accuracy. The present invention provides the following technical solution:

[0005] This invention provides a method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction, which includes the following steps:

[0006] S1: Acquire high-resolution optical images, low-resolution gravity field images, radar images, and terrain elevation images of the target monitoring area, perform spatiotemporal benchmark unified processing on each image data, and generate a multi-channel standard image set.

[0007] S2: Using the high-resolution optical image and the terrain elevation image, perform a spatial thinning operation based on texture features on the low-resolution gravity field image to output a high-resolution land water storage distribution image that incorporates surface details.

[0008] S3: Identify the surface water coverage area from the multi-channel standard image set, and remove the corresponding surface water storage signal from the high-resolution land water storage distribution image;

[0009] S4: Soil water storage is inverted based on surface features, and the corresponding soil water storage signals are further removed from the image to generate a groundwater storage change image;

[0010] S5: Based on the measured data, the initial water yield parameters of the region are spatially corrected, and the groundwater storage change image is converted into a dynamic change distribution image of groundwater level using the corrected water yield parameters.

[0011] S6: Perform spatial connectivity analysis on the dynamic distribution image of the groundwater level to identify abnormal areas and generate situation assessment results.

[0012] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: the spatiotemporal reference unification processing of each image data in S1 specifically includes:

[0013] The high-resolution optical image is selected as the reference base map, and the same geographic feature points in each image are identified to construct a coordinate transformation matrix. The pixel grid of the low-resolution gravity field image and the radar image is then resampled using the coordinate transformation matrix to match the reference base map. Figure 1 For invalid pixel areas in the high-resolution optical image that are obscured by clouds, valid pixel values ​​at the same coordinate position in adjacent shooting times are retrieved and filled in to replace them.

[0014] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: the spatial thinning operation based on texture features performed on the low-resolution gravity field image in step S2 specifically includes:

[0015] A data processing network structure including a feature extraction unit and an image reconstruction unit is constructed. The vegetation index features, surface temperature features, and elevation features in the multi-channel standard image set are used as high-frequency guiding information, and the low-resolution gravity field image after spatial resolution alignment is used as low-frequency basic information. These are input into the data processing network structure. The data processing network structure is used to perform nonlinear feature fusion operations to map the surface texture features in the high-frequency guiding information to the low-frequency basic information, thereby generating the high-resolution terrestrial water storage distribution image that incorporates surface texture details.

[0016] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: step S3, removing the corresponding surface water storage signal from the high-resolution terrestrial water storage distribution image, specifically includes:

[0017] The multi-band spectral index features of the optical image and the microwave scattering intensity features of the radar image are extracted; the gray-level distribution histogram of the image pixel values ​​is statistically analyzed to determine the segmentation threshold for distinguishing water and non-water pixels, and a binarized water mask image is generated using the segmentation threshold; the number of pixels representing water in the water mask image is counted, and the volume change corresponding to the pixels representing water in the water mask image is obtained by combining the elevation difference in the topographic elevation image, thereby generating a surface water storage distribution layer; the surface water storage distribution layer is then stripped from the high-resolution land water storage distribution image using pixel-level differential extraction.

[0018] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: the soil water storage inversion based on surface features in step S4 specifically includes:

[0019] A two-dimensional feature space is constructed with surface temperature as the vertical axis and vegetation index as the horizontal axis. The upper and lower boundary lines of the scatter distribution in the two-dimensional feature space are fitted. The relative position coefficient of each pixel between the upper and lower boundary lines is determined to obtain a surface soil moisture image. Numerical recursion processing on the time series is performed to simulate the delayed process of water infiltration from the surface to the deep layers. The root zone soil water storage distribution layer is calculated from the surface soil moisture image.

[0020] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: step S4, which involves further removing the corresponding soil water storage signal from the image, specifically includes:

[0021] The root zone soil water storage distribution layer is stripped pixel-level from the high-resolution terrestrial water storage distribution image after removing surface water signals to generate the groundwater storage change image reflecting the gravity changes of deep underground aquifers.

[0022] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: the spatial correction of the initial water yield parameter of the region based on measured data in step S5 specifically includes:

[0023] A digital lithology classification map and soil texture classification map of the target monitoring area are acquired. A corresponding conversion coefficient is assigned to each pixel based on its soil and rock type to generate an initial water yield parameter distribution map. Water level observation data from discretely distributed measured points within the target monitoring area are acquired, and the true water yield conversion coefficient at each measured point is determined. The deviation of the true water yield conversion coefficient at each measured point relative to the corresponding value in the initial water yield parameter distribution map is determined. The deviation value of the discretely distributed measured points is smoothly diffused across geographic space to construct a continuous spatial error correction surface covering the entire area. The initial water yield parameter distribution map is superimposed on the spatial error correction surface to obtain the corrected water yield parameter distribution map.

[0024] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: step S5 converts the groundwater storage change image into a groundwater level dynamic change distribution image using the corrected water yield parameter, specifically including:

[0025] Using the corrected water yield parameter distribution map, a pixel-by-pixel parametric inversion calculation is performed on the groundwater storage change image to obtain a dynamic distribution image of groundwater level changes.

[0026] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: step S6 involves performing spatial connectivity analysis on the groundwater level dynamic change distribution image, specifically including:

[0027] Obtain the groundwater level dynamic change distribution image relative to the historical reference image of the same period; filter out the set of pixels in the groundwater level change difference image whose values ​​are lower than a preset outlier; search for pixel clusters that are spatially adjacent to each other in the pixel set, mark each independent pixel cluster as a connected region, and calculate the coverage area and average pixel value of each connected region.

[0028] As a preferred embodiment of the groundwater level dynamic assessment method based on multi-source data mapping and fluid model correction described in this invention, wherein: the identification of abnormal areas and generation of situation assessment results in step S6 specifically includes:

[0029] The coverage area and average pixel value of the connected region are compared with a preset risk threshold. Connected regions whose index values ​​exceed the preset risk threshold are identified as groundwater funnel risk zones. The boundary contour of the groundwater funnel risk zone is highlighted in the output map results to generate situation assessment data containing the geographical range information of the groundwater funnel risk zone.

[0030] The beneficial effects of this invention are as follows: it effectively overcomes the limitation of low resolution in satellite gravity data, which makes it impossible to identify local details, and refines the macroscopic water volume change signal to the specific surface texture scale; by removing interference signals from surface water bodies and soil moisture layer by layer, it eliminates the noise influence of non-groundwater factors on the monitoring results, ensuring that the inversion results can truly reflect the changes in deep aquifers; at the same time, by fusing measured data to dynamically correct geological parameters, it significantly improves the accuracy of groundwater level numerical calculation, and realizes precise positioning and automatic early warning of groundwater over-extraction funnel areas. Attached Figure Description

[0031] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the 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.

[0032] Figure 1 The flowchart shows a method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction.

[0033] Figure 2 This is a logic diagram for image fusion based on deep learning.

[0034] Figure 3 Flowchart for layered stripping of hydrological signals;

[0035] Figure 4 Flowchart for spatial error correction of water yield parameters;

[0036] Figure 5 This is a flowchart for anomaly identification and situation assessment. Detailed Implementation

[0037] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0038] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0039] Secondly, the term "one embodiment" or "example" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the invention. The appearance of an embodiment in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that mutually excludes other embodiments.

[0040] Example 1

[0041] Reference Figures 1-5 This is one embodiment of the present invention, which provides a method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction, including the following steps:

[0042] S1. Acquire high-resolution optical images, low-resolution gravity field images, radar images, and terrain elevation images of the target monitoring area. Perform spatiotemporal benchmark unified processing on each image data to generate a multi-channel standard image set.

[0043] High-resolution optical images are selected as the baseline reference map. Common geographic feature points are identified in each image, and a coordinate transformation matrix is ​​constructed. This matrix is ​​then used to resample the pixel grid of the low-resolution gravity field image and radar image to match the baseline reference map. Figure 1 To address invalid pixel regions obscured by clouds in high-resolution optical images, valid pixel values ​​at the same coordinate position in adjacent shooting times are retrieved and used to fill and replace them.

[0044] Addressing the inconsistencies in spatial resolution, imaging time, and coordinate systems among different satellite sensors includes:

[0045] 1. Specific selection of data source:

[0046] High-resolution optical images: Wide-field images (16-meter resolution) from my country's independently developed Gaofen-1 or Gaofen-6 satellites, or Sentinel-2 images, are preferred for extracting surface textures. In this embodiment, the time window for each image data is preferably the same ten-day period to ensure the consistency of hydrological conditions.

[0047] Low-resolution gravity field image: Monthly time-varying gravity field data (raw resolution approximately 300 km) released by the Gravity Inversion and Climate Experiment satellite were used to obtain the total water storage signal.

[0048] Radar imagery: Imagery from my country's Gaofen-3 or Sentinel-1 synthetic aperture radar is preferred, utilizing their C-band sensitivity to water bodies for auxiliary observation.

[0049] Topographic elevation image: Using a digital elevation model, this image records the elevation value of every pixel on the ground. In this embodiment, 90-meter resolution data released by the Space Shuttle Radar Topographic Mapping mission can be selected.

[0050] 2. Geometric fine correction and mesh reconstruction:

[0051] The system uses the WGS84 coordinate system of high-resolution optical images as a reference. For radar images, the system extracts feature points such as road intersections and river bifurcation points, and uses the scale-invariant feature transformation algorithm to calculate the affine transformation matrix to geometrically correct them to the reference coordinate system.

[0052] For gravity field images, the system employs bicubic interpolation. Specifically, the system first constructs a target grid (e.g., a 1000×1000 matrix) with the same resolution as the optical image. For each interpolation point (x, y) in the target grid, the system finds its corresponding floating-point coordinates in the original low-resolution gravity image (e.g., a 10×10 matrix) and selects 16 original pixels within a 4×4 neighborhood of those coordinates. The system then calculates the weighting coefficients for these 16 points using bicubic basis functions, performs a weighted summation operation, and thus calculates the value of the interpolation point, ultimately smoothing the coarse gravity signal onto the high-resolution grid.

[0053] 3. Cloud removal and data completion:

[0054] To address cloud cover obstruction in optical images, the system employs a spatiotemporal filling strategy. First, it retrieves cloud-free observations from the preceding and following 15 days for the pixel and replaces them.

[0055] If clouds are present throughout the time window, then multi-source correlation-based grayscale mapping completion is performed: the system selects nearby cloudless areas in the image as training samples and extracts the grayscale values ​​of the optical image within those areas. and the backscattering coefficient of radar images The least squares method is used to fit a linear regression equation between the two:

[0056]

[0057] in is the regression coefficient.

[0058] For cloud-covered areas, the system reads the corresponding radar backscattering coefficient, substitutes it into the equation, and calculates the predicted optical grayscale value, thereby restoring the texture information of the obscured ground objects.

[0059] S2. Using high-resolution optical images and topographic elevation images, perform spatial thinning operations based on texture features on low-resolution gravity field images to output high-resolution terrestrial water storage distribution images that incorporate surface details.

[0060] High-resolution optical images are selected as the baseline reference map. Common geographic feature points are identified in each image, and a coordinate transformation matrix is ​​constructed. This matrix is ​​then used to resample the pixel grid of the low-resolution gravity field image and radar image to match the baseline reference map. Figure 1 For invalid pixel regions obscured by clouds in high-resolution optical images, valid pixel values ​​at the same coordinate position in adjacent time phases are retrieved and filled / replaced.

[0061] A data processing network structure containing feature extraction and image reconstruction units is constructed. Vegetation index features, surface temperature features, and elevation features from a multi-channel standard image set are used as high-frequency guiding information, and low-resolution gravity field images aligned with spatial resolution are used as low-frequency basic information. These are then input into the data processing network structure. The data processing network structure is used to perform nonlinear feature fusion operations to map the surface texture features in the high-frequency guiding information to the low-frequency basic information, generating a high-resolution terrestrial water storage distribution image that incorporates surface texture details.

[0062] This involves the nonlinear mapping and fusion of high-frequency spatial structure information into low-resolution gravity field images. The specific data processing steps are as follows:

[0063] 1. Dual-stream feature extraction:

[0064] This embodiment constructs a convolutional neural network with dual-stream input.

[0065] Low-frequency feature flow: Input the interpolated gravity field image, and extract large-scale background water storage features through convolutional layers to maintain the total water volume constraint.

[0066] High-frequency feature stream: Input normalized optical, temperature and elevation images, and extract fine surface spatial structure features (such as vegetation cover gradient and topographic water lines) through deep convolutional networks.

[0067] 2. Channel splicing and nonlinear fusion:

[0068] The system performs feature channel concatenation, stacking high-frequency and low-frequency feature maps. The data then enters the residual-dense module. During this process, the network performs pixel-level non-linear weighted operations. Specifically, the network establishes mappings using learned physical correlations: for pixel regions in the input features that exhibit high vegetation index, low surface temperature, and low-lying terrain, the network assigns them higher water storage weights in the gravity field feature map; conversely, it assigns lower weights to regions with low vegetation, high temperature, and high elevation.

[0069] Through this spatial weight redistribution mechanism based on surface environmental factors, high-frequency texture details of the surface are successfully injected and the originally blurred gravity field signal is refined.

[0070] 3. Residual learning and image reconstruction:

[0071] The network employs a global residual learning strategy, meaning it does not directly predict the final water volume but instead predicts a high-frequency detail residual map. The system performs matrix addition between the input original interpolated gravity image (base map) and this high-frequency detail residual map, thereby outputting a high-resolution terrestrial water storage distribution image that retains both the original gravity accuracy and clear texture details.

[0072] S3. Identify the surface water coverage area from the multi-channel standard image set, and remove the corresponding surface water storage signal from the high-resolution land water storage distribution image.

[0073] Multi-band spectral index features of optical images and microwave scattering intensity features of radar images are extracted. The gray-level distribution histogram of image pixel values ​​is statistically analyzed to determine the segmentation threshold for distinguishing water and non-water pixels. A binarized water mask image is generated using the segmentation threshold. The number of pixels representing water in the water mask image is counted, and the volume change corresponding to the pixels representing water in the water mask image is obtained by combining the elevation difference in the topographic elevation image, thus generating a surface water storage distribution layer. The surface water storage distribution layer is then stripped from the high-resolution terrestrial water storage distribution image using pixel-level differential extraction.

[0074] Among them, based on the principle of water balance, the first layer of signal stripping is performed.

[0075] 1. Adaptive water extraction:

[0076] The system first calculates the improved normalized difference water index of the optical image. The formula is:

[0077]

[0078] in, It is in the green light band, and water has a high reflectivity in this band; This is the shortwave infrared band, where water has extremely strong absorption characteristics (i.e., extremely low reflectivity), while land and buildings have high reflectivity. Simultaneously, the backscattering coefficient of the radar image is read. The feature map is fused with the radar feature map, and the gray-level histogram of all pixels in the image is calculated. Since water and non-water bodies exhibit a clear bimodal distribution in the histogram, the system automatically calculates the optimal segmentation threshold using the Otsu's algorithm (maximum inter-class variance method). Pixels larger than this threshold are marked as 1 (water body), and the rest are marked as 0, generating a binarized water body mask.

[0079] 2. Quantification of water volume:

[0080] In order to convert the water surface area into water volume, the system performs volume integration on each connected region identified as a water body.

[0081] The system identifies the edge pixels of the connected water area, extracts the average elevation of these pixels in the topographic elevation image, and uses this as the water level elevation for the current time phase. ).

[0082] Calculate surface water storage using the following integral formula. :

[0083]

[0084] in, This represents the summation of all pixels within a connected region; For the first The bottom elevation of each pixel; (area of ​​a single pixel).

[0085] 3. Signal stripping:

[0086] The generated surface water storage distribution layer is then stripped from the high-resolution terrestrial water storage distribution image using pixel-level differential extraction:

[0087]

[0088] in, It is the total water storage, which includes all forms of water (river and lake water + soil water + groundwater). This is what's left after the subtraction. At this point, the image no longer shows the weight of rivers and lakes; what remains is mainly soil water and groundwater.

[0089] In the processed images, the signals from river and lake areas were removed, and the remaining signals mainly reflect soil water and groundwater.

[0090] S4. Based on surface features, soil water storage is inverted, and the corresponding soil water storage signals are further removed from the image to generate a groundwater storage change image.

[0091] A two-dimensional feature space is constructed with surface temperature as the vertical axis and vegetation index as the horizontal axis. The upper and lower boundary lines of the scatter distribution in the two-dimensional feature space are fitted. The relative position coefficient of each pixel between the upper and lower boundary lines is determined to obtain the surface soil moisture image. Numerical recursion processing on the time series is performed to simulate the delayed process of water infiltration from the surface to the deep layers. The root zone soil water storage distribution layer is inferred from the surface soil moisture image.

[0092] The root zone soil water storage distribution layer is stripped pixel-level from the high-resolution terrestrial water storage distribution image after removing surface water signals to generate a groundwater storage change image reflecting the gravity changes of deep underground aquifers.

[0093] Since satellite sensors cannot directly detect deep soil moisture, the system uses a physical inference model for inversion.

[0094] 1. Surface humidity inversion:

[0095] The system is constructed based on surface temperature ( The two-dimensional feature space is defined by the vertical axis (I) and the horizontal axis (VI).

[0096] Iterate through all pixels in the image and fit the top edge (dry edge) of the scatter plot using linear regression. ) and bottom edge (wet edge) ).

[0097] For any pixel, calculate its Temperature Vegetation Drought Index (TVDI) as a proxy for surface (0-5 cm) soil moisture:

[0098]

[0099] in, This represents the measured surface temperature at that point.

[0100] 2. Deep infiltration extrapolation (exponential filtering):

[0101] To simulate the time lag effect of water infiltration from the surface to deeper layers, the system employs a recursive exponential filtering algorithm. The soil moisture index in the root zone (0-100 cm) is calculated from the surface temperature and vegetation drought index. ):

[0102] in, For time step, This is the gain factor, and the time constant is... Depending on the soil texture, in this embodiment, a period of 20 days is preferably set for loam areas. for The root zone soil moisture index at a given time. for The surface humidity at any given time. This calculation generates a layer showing the distribution of soil water storage in the root zone. ).

[0103] 3. Final stripping:

[0104] The root zone soil water storage distribution layer is then peeled off from the image processed in the previous step using pixel-level differential extraction:

[0105]

[0106] in, Indicates the amount of deep groundwater stored. This indicates the remaining total water storage.

[0107] At this point, the remaining gravity anomaly signals in the image are caused only by changes in the reserves of deep underground aquifers, generating a groundwater reserve change image.

[0108] S5. Based on the measured data, the initial water yield parameters of the region are spatially corrected, and the groundwater storage change image is converted into a dynamic change distribution image of groundwater level using the corrected water yield parameters.

[0109] Digital lithology and soil texture classification maps of the target monitoring area are acquired. A conversion coefficient is assigned to each pixel based on its soil and rock type to generate an initial water yield parameter distribution map. Water level observation data from discretely distributed measured points within the target monitoring area are acquired, and the true water yield conversion coefficient at each measured point is determined. The deviation of the true water yield conversion coefficient at each measured point from the corresponding value in the initial water yield parameter distribution map is determined. The deviation values ​​of the discretely distributed measured points are smoothly diffused across geographic space to construct a continuous spatial error correction surface covering the entire area. The initial water yield parameter distribution map and the spatial error correction surface are superimposed to obtain the corrected water yield parameter distribution map.

[0110] Using the corrected water yield parameter distribution map, a pixel-by-pixel parametric inversion solution is performed on the groundwater storage change image to obtain the dynamic change distribution image of groundwater level.

[0111] Among them, a data assimilation strategy is used to solve the nonlinear problem of mass (weight) to water level (height) conversion.

[0112] 1. Parameter initialization and deviation determination:

[0113] The system first reads a 1:200,000 digital hydrogeological map, and assigns empirical water yield values ​​based on the lithology of each pixel (such as coarse sand, silt, bedrock), generating an initial parameter map. .

[0114] At discretely distributed measured well locations, the actual water yield is calculated backward using the water balance principle. :

[0115] in For changes in reserves retrieved by satellite, This is a measured change in water level.

[0116] Calculate the local deviation value: .

[0117] 2. Spatial error correction:

[0118] The system first bases its data on the deviation values ​​of discrete measured points. The experimental variogram is calculated to quantify the spatial correlation of deviation values ​​(i.e., the closer two points are, the more similar their deviation values). This embodiment preferably uses a spherical model as the theoretical variogram for fitting, determining key structural parameters such as nugget value (representing micro-random error), sill value (representing the maximum degree of variation), and range (representing the maximum distance of spatial correlation). This captures the spatial distribution pattern of groundwater conversion errors in the region.

[0119] Kriging interpolation: Using the ordinary Kriging interpolation algorithm, based on the fitted variogram model, the deviation of discrete points is interpolated. The error is diffused across the entire region. The system calculates weights based on the spatial distance and structural relationship between the pixel to be interpolated and known points, estimates the deviation value of each pixel in the entire region, and generates a continuous spatial error correction surface.

[0120] Parametric overlay: The initial parametric map is matrix-overlaid with the error correction surface to obtain the corrected parametric map. .

[0121] It is worth noting that the corrected water yield here is essentially a generalized regional hydrogeological conversion coefficient. It is not limited to characterizing lithological physical properties, but rather, by introducing measured bias data, it implicitly absorbs and compensates for systematic errors caused by complex physical factors such as differences in aquifer structure (e.g., the elastic effect of confined water), redistribution of crustal mass caused by land subsidence, and crustal isostatic adjustment, thereby significantly improving the adaptability and accuracy of the inversion model.

[0122] 3. Water level mapping:

[0123] Using the corrected water yield parameter distribution map, a pixel-by-pixel parametric inversion solution is performed on the groundwater storage change image:

[0124] in: This represents the final change in groundwater level; Changes in groundwater reserves; This is the corrected water yield parameter. The result is a high-precision groundwater level distribution image output in meters (m).

[0125] S6. Perform spatial connectivity analysis on the dynamic distribution image of groundwater level changes, identify abnormal areas and generate situation assessment results.

[0126] The system acquires images showing the dynamic distribution of groundwater levels relative to historical baseline images of the same period; filters out sets of pixels in the images showing the dynamic distribution of groundwater levels with values ​​lower than preset outliers; searches for spatially adjacent pixel clusters within these pixel clusters, marks each independent pixel cluster as a connected region, and calculates the coverage area and average pixel value of each connected region.

[0127] The coverage area and average pixel value of the connected regions are compared with the preset risk threshold. Connected regions whose index values ​​exceed the preset risk threshold are identified as groundwater funnel risk areas. The boundary contour of the groundwater funnel risk areas is highlighted in the output map results, generating situation assessment data containing the geographical range information of the groundwater funnel risk areas.

[0128] Among these are the intelligent applications of monitoring results and risk warnings.

[0129] 1. Anomaly Detection:

[0130] The system calculates the difference between the current groundwater level distribution map and the average groundwater level map of the same period over the past 5 years, and obtains a map of differences in groundwater level changes.

[0131] Filter out the set of pixels in the difference image whose values ​​are lower than a preset outlier value (e.g., -2 meters) and mark them as candidate outliers.

[0132] 2. Connectivity analysis:

[0133] The system performs 8-neighborhood connectivity analysis on the binarized abnormal image (i.e., determines whether there are similar pixels in the 8 directions around a pixel).

[0134] Spatially connected outlier pixels are merged into an independent connected region, and two key metrics for each connected region are calculated: coverage area (number of pixels × pixel area) and average descent depth (arithmetic mean of pixel values).

[0135] 3. Situation assessment:

[0136] The indicators of the connected regions are compared with the preset risk thresholds.

[0137] Judgment rule: If the coverage area of ​​a connected area is greater than 10 square kilometers and the average drop depth is more than 3 meters, it is judged as a groundwater funnel risk area.

[0138] Output results: The system automatically extracts the vector boundary of the risk area, highlights it in red on the map, and generates a situation assessment data report containing the center coordinates of the area, the scope of influence, and the risk level.

[0139] Example 2

[0140] This is the second embodiment of the present invention. To further illustrate the present invention, the technical solution of the present invention will be described in more detail below using an agricultural irrigation area in the North China Plain (hereinafter referred to as the experimental area) as an example.

[0141] In this embodiment, the local water resources department needs to conduct detailed monitoring of the dynamic changes in groundwater level in the experimental area (approximately 2,000 square kilometers) during the peak irrigation season in the summer of 2023, and identify potential groundwater funnel risk areas.

[0142] S1: Multi-source data acquisition and spatiotemporal normalization preprocessing

[0143] First, we acquired multi-source remote sensing data for the experimental area in June 2023.

[0144] Optical image: A wide-field image from the Gaofen-1 satellite with a resolution of 16 meters clearly records the wheat-corn rotation planting structure in this area.

[0145] Gravity field image: Acquired monthly gravity anomaly data released by the Gravity Inversion and Climate Experiment (GRACE) satellite during the same period. Although the raw resolution is only 300 km, it provides a macroscopic signal of the reduction in total water storage in the region.

[0146] Radar and Topographic Data: Acquire Gaofen-3 synthetic aperture radar imagery and 90-meter digital elevation model data released by the Space Shuttle radar topographic mapping mission.

[0147] The system uses Gaofen-1 imagery as a baseline and employs a scale-invariant feature transform algorithm to register road intersection feature points in radar imagery. For the approximately 15% cloud-occupied area in the optical imagery, the system retrieves cloud-free imagery from the preceding and following 10 days and replaces it with such imagery. It also uses radar texture features from the same period to fill in the remaining gaps, generating a seamless multi-channel standard image set.

[0148] S2: Deep Learning-Based Super-Resolution Reconstruction of Gravity Field Images

[0149] The system constructs a multi-scale residual dense network. Normalized vegetation index (reflecting crop water requirements), surface temperature (reflecting soil drought level), and topographic elevation from the standard image set are used as high-frequency guiding information and input into the network.

[0150] Based on pre-trained mapping rules, the network determined that the farmland area in the central part of the experimental zone had an extremely high vegetation index (vigorous crop growth) and a low surface temperature (strong evaporation due to irrigation). Therefore, the network automatically assigned a higher weight to gravity water storage in this area. Ultimately, the fuzzy gravity signal was intelligently focused onto specific farmland plots, outputting a high-precision terrestrial water storage distribution map with a resolution of 1 kilometer.

[0151] S3: Precise identification and stripping of surface water signals

[0152] The experimental area includes a main river and two small reservoirs. The system automatically determines the segmentation threshold using the Otsu's method and accurately extracts the water mask from the improved normalized difference water index feature map. Combined with digital elevation model data, the system calculates the water level elevations of the two reservoirs in June to be 52 meters and 48 meters, respectively, and then calculates the surface water storage to be approximately 0.5 billion cubic meters through integration. The system performs pixel-level subtraction to remove this surface water signal from the total water storage, avoiding misclassification as groundwater.

[0153] S4: Inversion and Secondary Stripping of Soil Water Signals

[0154] The system constructs a temperature-vegetation feature space, fits dry and wet edge equations, and calculates the average temperature-vegetation drought index of the experimental area as 0.65, indicating that the surface soil is relatively dry.

[0155] Using exponential filtering, the system calculated the changes in soil water storage in the root zone (0-100 cm). Since it was the irrigation season, the soil water storage in the root zone showed an increasing trend. The system then performed subtraction to remove the soil water signal from the image. At this point, the remaining gravity anomaly signal purely reflected changes in deep underground aquifers.

[0156] S5: Spatial Correction of Water Supply Parameters and Water Level Mapping

[0157] The system reads the geological map, which shows that the northern part of the experimental area is coarse sand (theoretical water yield 0.25) and the southern part is silty clay (theoretical water yield 0.05).

[0158] The system incorporated data from five measured wells within the region for correction. Calculations revealed a significant discrepancy between satellite-retrieved and measured values ​​in the southern clay region. The system used ordinary kriging interpolation to generate an error surface, correcting the specific yield parameter in the south to 0.07. Using the corrected parameter map, the system performed a division operation to convert water weight into water level, generating a dynamic distribution map of groundwater level changes in June 2023. The results showed that the overall water level in the experimental area decreased by an average of 1.2 meters.

[0159] S6: Connectivity Component Intelligent Analysis and Risk Assessment

[0160] The system calculates the difference between the water level map and the historical average for the same period (the average of the past 5 years).

[0161] Using a connected component analysis algorithm, the system automatically identified an anomalous connected region in the center of the experimental area on the image. Statistics show that this region covers an area of ​​35 square kilometers, with an average water level drop of 2.8 meters.

[0162] By comparing this indicator with preset thresholds (area > 10 square kilometers, depth > 2 meters), the system determines that the area is a high-risk groundwater funnel zone. The system automatically extracts the vector boundary of the funnel zone and highlights it in red on the map, while generating a report suggesting that agricultural irrigation water use in the area be restricted immediately.

[0163] In summary, this invention effectively overcomes the limitation of low resolution in satellite gravity data, which makes it impossible to identify local details, and refines macroscopic water volume change signals to specific surface texture scales. By eliminating interference signals from surface water bodies and soil moisture layer by layer, it eliminates the noise impact of non-groundwater factors on monitoring results, ensuring that the inversion results can truly reflect changes in deep aquifers. At the same time, by fusing measured data to dynamically correct geological parameters, it significantly improves the accuracy of groundwater level calculations, and achieves precise location and automatic early warning of groundwater over-extraction funnel areas.

[0164] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for groundwater level dynamic assessment based on multi-source data mapping and fluid model correction, characterized in that, Includes the following steps: S1: Acquire high-resolution optical images, low-resolution gravity field images, radar images, and terrain elevation images of the target monitoring area, perform spatiotemporal benchmark unified processing on each image data, and generate a multi-channel standard image set. S2: Construct a data processing network structure that includes a feature extraction unit and an image reconstruction unit; use the vegetation index features, surface temperature features, and elevation features from the multi-channel standard image set as high-frequency guiding information, and use the low-resolution gravity field image after spatial resolution alignment as low-frequency basic information, and input them into the data processing network structure; use the data processing network structure to perform nonlinear feature fusion operations, map the surface texture features in the high-frequency guiding information to the low-frequency basic information, and generate a high-resolution land water storage distribution image that incorporates surface texture details; S3: Identify the surface water coverage area from the multi-channel standard image set, and remove the corresponding surface water storage signal from the high-resolution land water storage distribution image; S4: Soil water storage is inverted based on surface features, and the corresponding soil water storage signals are further removed from the image to generate a groundwater storage change image; S5: Obtain digital lithology and soil texture classification maps of the target monitoring area; assign corresponding conversion coefficients to each pixel based on its soil and rock type to generate an initial water yield parameter distribution map; acquire water level observation data from discretely distributed measured points within the target monitoring area; determine the true water yield conversion coefficient at each measured point; determine the deviation of the true water yield conversion coefficient at each measured point relative to the corresponding value in the initial water yield parameter distribution map; smoothly diffuse the deviation values ​​of the discretely distributed measured points in geographic space to construct a continuous spatial error correction surface covering the entire area; superimpose the initial water yield parameter distribution map with the spatial error correction surface to obtain a corrected water yield parameter distribution map; use the corrected water yield parameter distribution map to perform pixel-by-pixel parametric inversion calculations on the groundwater storage change image to obtain a dynamic distribution image of groundwater level changes; S6: Perform spatial connectivity analysis on the dynamic distribution image of the groundwater level to identify abnormal areas and generate situation assessment results.

2. The method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction according to claim 1, characterized in that, The process of performing spatiotemporal benchmark unification processing on each image data in S1 specifically includes: The high-resolution optical image is selected as the reference base map, and the same geographical feature points in each image are identified to construct a coordinate transformation matrix. The pixel grid of the low-resolution gravity field image and the radar image is resampled to be consistent with the reference base map using the coordinate transformation matrix. For invalid pixel areas in the high-resolution optical image that are obscured by clouds, valid pixel values ​​at the same coordinate position in adjacent shooting times in the time series are retrieved and filled and replaced.

3. The method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction according to claim 1, characterized in that, The step S3, which involves removing the corresponding surface water storage signal from the high-resolution land water storage distribution image, specifically includes: The multi-band spectral index features of the optical image and the microwave scattering intensity features of the radar image are extracted; the gray-level distribution histogram of the image pixel values ​​is statistically analyzed to determine the segmentation threshold for distinguishing water and non-water pixels, and a binarized water mask image is generated using the segmentation threshold; the number of pixels representing water in the water mask image is counted, and the volume change corresponding to the pixels representing water in the water mask image is obtained by combining the elevation difference in the topographic elevation image, thereby generating a surface water storage distribution layer; the surface water storage distribution layer is then stripped from the high-resolution land water storage distribution image using pixel-level differential extraction.

4. The method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction according to claim 1, characterized in that, The soil water storage inversion based on surface features in S4 specifically includes: A two-dimensional feature space is constructed with surface temperature as the vertical axis and vegetation index as the horizontal axis. The upper and lower boundary lines of the scatter distribution in the two-dimensional feature space are fitted. The relative position coefficient of each pixel between the upper and lower boundary lines is determined to obtain a surface soil moisture image. Numerical recursion processing on the time series is performed to simulate the delayed process of water infiltration from the surface to the deep layers. The root zone soil water storage distribution layer is calculated from the surface soil moisture image.

5. The method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction according to claim 4, characterized in that, The step S4, which involves further removing the corresponding soil water storage signal from the image, specifically includes: The root zone soil water storage distribution layer is stripped pixel-level from the high-resolution terrestrial water storage distribution image after removing surface water signals to generate the groundwater storage change image reflecting the gravity changes of deep underground aquifers.

6. The method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction according to claim 1, characterized in that, S6 involves performing spatial connectivity analysis on the dynamic distribution image of the groundwater level, specifically including: Obtain the groundwater level dynamic change distribution image relative to the historical reference image of the same period; filter out the set of pixels in the groundwater level change difference image whose values ​​are lower than a preset outlier; search for pixel clusters that are spatially adjacent to each other in the pixel set, mark each independent pixel cluster as a connected region, and calculate the coverage area and average pixel value of each connected region.

7. The method for dynamic assessment of groundwater level based on multi-source data mapping and fluid model correction according to claim 6, characterized in that, The process of identifying abnormal regions and generating situation assessment results in S6 specifically includes: The coverage area and average pixel value of the connected region are compared with a preset risk threshold. Connected regions whose index values ​​exceed the preset risk threshold are identified as groundwater funnel risk zones. The boundary contour of the groundwater funnel risk zone is highlighted in the output map results to generate situation assessment data containing the geographical range information of the groundwater funnel risk zone.

Images