Groundwater level spatio-temporal dynamic monitoring and visualization analysis method and system thereof

Abstract

The application discloses a groundwater level space-time dynamic monitoring and visual analysis method and system, relates to the technical field of groundwater resource monitoring and visualization, and comprises the following six steps: multi-source data fusion preprocessing, mixed space interpolation of the Kriging and machine learning, isopachous surface generation and optimization, time sequence animation synthesis, abnormality identification and classification, and early warning heat map visualization. The system comprises a data acquisition module, a space interpolation module, an isopachous surface generation module, a time sequence animation module, an abnormality identification module and an early warning visualization module. The application improves the space interpolation precision through multi-source data fusion and mixed interpolation of the Kriging and machine learning, directly presents the space-time evolution of the water level through the isopachous surface and the time sequence animation, identifies abnormal areas such as over-exploitation funnel through space-time gradient analysis, and assists in management decision through the early warning heat map, so as to solve the technical problems of sparse monitoring sites and the difficulty of directly presenting the space-time evolution law in traditional reports.

Description

Technical Field

[0001] This invention relates to the field of groundwater resource monitoring and visualization technology, specifically to a method and system for spatiotemporal dynamic monitoring and visualization analysis of groundwater levels. Background Technology

[0002] Groundwater resources, as an important freshwater resource, play an irreplaceable role in agricultural irrigation, urban water supply, and ecological maintenance. With rapid economic and social development, groundwater extraction has continuously increased, leading to problems such as persistently declining water levels and the expansion of over-extraction funnels in some areas, seriously threatening regional water resource security and ecological environment stability. Therefore, achieving accurate monitoring and intuitive visualization analysis of the spatiotemporal dynamic changes of regional groundwater levels is of great significance for the scientific management and rational allocation of groundwater resources.

[0003] Currently, groundwater level monitoring mainly relies on automatic monitoring stations and manual measurement stations distributed throughout the monitoring area. However, constrained by factors such as economic costs, terrain conditions, and technological levels, the spatial distribution of monitoring stations is often sparse and uneven, making it difficult to comprehensively reflect the spatial distribution characteristics of regional groundwater levels. Furthermore, traditional groundwater level data presentation methods often use tables and graphs, which, while showing the water level change trends of individual or a small number of stations, fail to intuitively present the spatiotemporal evolution patterns of regional water levels, limiting the scientific rigor and timeliness of management decisions. Existing technologies primarily focus on visualizing linear elements of surface water rivers, employing simple color gradient interpolation methods that do not consider spatial autocorrelation and the influence of auxiliary variables, making it difficult to accurately reflect the spatial distribution characteristics of groundwater levels as a continuous surface variable. In addition, the lack of multi-source heterogeneous data fusion mechanisms, time-series animation display functions, and anomaly area identification and early warning capabilities fails to meet the practical needs of spatiotemporal dynamic monitoring and visualization analysis of groundwater levels.

[0004] Therefore, it is necessary to provide a method and system for dynamic monitoring and visualization analysis of groundwater levels in the spatiotemporal region that can integrate multi-source monitoring data, achieve high-precision spatial interpolation, intuitively present the spatiotemporal evolution pattern, and have anomaly identification and early warning functions. Summary of the Invention

[0005] This invention provides a method and system for spatiotemporal dynamic monitoring and visualization analysis of groundwater levels, which solves the technical problems in the prior art, such as incomplete spatial distribution information caused by the sparse groundwater level monitoring stations, the difficulty of intuitively presenting spatiotemporal evolution patterns in traditional report formats, and the lack of anomaly identification and early warning functions. It achieves the technical effects of multi-source data fusion, high-precision spatial interpolation, spatiotemporal dynamic visualization, and anomaly early warning.

[0006] According to a first aspect of the present invention, the present invention provides a method for spatiotemporal dynamic monitoring and visualization analysis of groundwater levels. The method includes: a multi-source data fusion preprocessing step, acquiring multi-source groundwater level monitoring data for a region, performing time-series alignment processing on the multi-source groundwater level monitoring data to obtain a time-series alignment processing result, and performing quality-weighted fusion of the time-series alignment processing result based on the measurement uncertainty of each data source to obtain a fused groundwater level dataset; a cokriging and machine learning hybrid spatial interpolation step, acquiring auxiliary variable data related to the spatial distribution of groundwater levels, constructing a cokriging interpolation model based on the auxiliary variable data and the fused groundwater level dataset and performing spatial interpolation to obtain a preliminary cokriging interpolation result, using a gradient boosting regression model to correct the residuals of the preliminary cokriging interpolation result to obtain a residual correction result, and superimposing the preliminary cokriging interpolation result and the residual correction result to obtain a spatial interpolation result for groundwater levels; and isosurface generation and optimization. The process involves several steps: First, a contour extraction algorithm is used to generate water level contour data based on spatial interpolation results. This data is then simplified and smoothed to create a water level contour map. Second, a time-series animation synthesis step involves serializing the water level contour maps from multiple time points to obtain a multi-temporal water level distribution map sequence. This sequence is then subjected to inter-frame interpolation and synthesized into time-series animation data. Third, anomaly identification and classification steps involve calculating the spatiotemporal gradient features of water level changes based on the multi-temporal water level distribution map sequence. Anomaly indices are calculated based on these features and compared with a preset threshold for anomaly detection. Anomaly regions are then classified based on the anomaly detection results to obtain anomaly region identification results. Finally, an early warning heatmap visualization step involves generating an early warning heatmap based on the anomaly region identification results and the water level decline rate calculated based on the spatiotemporal gradient features. Areas exceeding the warning water level are then marked on the early warning heatmap.

[0007] According to a second aspect of the present invention, the present invention also provides a groundwater level spatiotemporal dynamic monitoring and visualization analysis system, the system comprising: a data acquisition module for acquiring multi-source groundwater level monitoring data in a region and performing time-series alignment and quality-weighted fusion; a spatial interpolation module for constructing a cokriging interpolation model based on auxiliary variable data and performing residual correction by combining gradient boosting regression; an isosurface generation module for generating a water level isosurface map based on the spatial interpolation results; a temporal animation module for performing inter-frame interpolation on multi-temporal water level distribution maps and synthesizing a temporal animation; an anomaly identification module for identifying and classifying abnormal areas based on spatiotemporal gradient features; and an early warning visualization module for generating an early warning heat map and marking areas exceeding the warning water level.

[0008] The beneficial effects of this invention are as follows: First, through multi-source data fusion preprocessing, it effectively integrates heterogeneous data from multiple sources, including automatic monitoring stations, manual measurements, and remote sensing inversion; quality-weighted fusion based on measurement uncertainty improves data reliability. Second, the hybrid interpolation method combining cokriging and machine learning fully utilizes the spatial correlation of auxiliary variables such as topography and land use, and further improves interpolation accuracy through residual correction. Third, the isosurface generation and temporal animation synthesis functions enable intuitive visualization of the spatiotemporal evolution of groundwater levels, facilitating managers to quickly grasp the overall picture of dynamic changes in water levels. Fourth, the anomaly identification algorithm based on spatiotemporal gradient tensor analysis can effectively identify typical abnormal areas such as over-extraction funnel expansion, abnormal water level decline, and enhanced recharge. Fifth, the early warning heat map visualization presents the rate of water level decline and the distribution of areas exceeding the warning level with intuitive color coding, providing decision support for dynamic monitoring and scientific scheduling of groundwater resources. Attached Figure Description

[0009] Figure 1 This is a flowchart illustrating the spatiotemporal dynamic monitoring and visualization analysis method for groundwater levels provided by this invention.

[0010] Figure 2 This is a schematic diagram of the structure of the groundwater level spatiotemporal dynamic monitoring and visualization analysis system provided by the present invention. Detailed Implementation

[0011] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0012] Reference Figure 1 As shown in the figure, this invention provides a method for spatiotemporal dynamic monitoring and visualization analysis of groundwater levels. This method mainly includes six core steps, forming a deeply coupled data processing link between each step. The output of the previous step serves as the key input for the next step, and the analysis results of subsequent steps can inversely influence the parameter configuration of the preceding steps, forming a closed-loop optimization mechanism. The following will describe each step in detail with reference to specific embodiments.

[0013] Step S1: Multi-source data fusion preprocessing step.

[0014] In one embodiment of the present invention, the multi-source data fusion preprocessing step is the data foundation of the entire method. Its core objective is to integrate groundwater level observation data from different monitoring methods and generate a unified and reliable fused water level dataset through time-series alignment and quality-weighted fusion mechanisms, thus laying the data foundation for subsequent spatial interpolation analysis.

[0015] Specifically, this step first acquires multi-source groundwater level monitoring data for the region. Preferably, the multi-source groundwater level monitoring data includes three main data sources: real-time groundwater level data from automatic monitoring stations, manually measured groundwater level data, and groundwater storage change data retrieved via remote sensing. Automatic monitoring stations typically use pressure-type or float-type water level sensors to continuously record groundwater level changes at sampling intervals of 15 minutes to 1 hour. This data has high temporal resolution but limited spatial coverage. Manually measured groundwater level data comes from periodic inspections, with sampling frequencies typically weekly or monthly. Although the temporal resolution is lower, it can cover remote areas where automatic stations are not deployed. Remote sensing data retrieval of groundwater storage changes mainly comes from the inversion results of GRACE satellite gravity measurement data. This data has a temporal resolution of approximately one month and a spatial resolution of approximately 300 km × 300 km, providing macroscopic trend information on groundwater storage changes at a regional scale.

[0016] Because the sampling frequencies of different data sources vary significantly, direct fusion will lead to data mismatch in the time dimension. Therefore, this step performs time-series alignment processing on the multi-source water level monitoring data. Preferably, the time-series alignment processing includes unifying data with different sampling frequencies to the same time reference, where the time interval can be set from 1 hour to 24 hours according to actual application requirements. In one embodiment of the invention, the time reference is set to 6 hours, meaning that water level observation data for four time nodes are generated daily. For automatic monitoring station data with sampling frequencies higher than the time reference, downsampling processing is performed using the time window averaging method, for example, taking the arithmetic mean of multiple observations within a continuous 6 hours as the representative value for that time node. For manually measured data and remote sensing inversion data with sampling frequencies lower than the time reference, linear interpolation or spline interpolation methods are used to densify the time dimension, generating an interpolation sequence that matches the time reference.

[0017] After time-series alignment, water level observations from different data sources need to be fused. Considering the differences in measurement accuracy and reliability among different monitoring methods, a simple arithmetic average fusion method cannot guarantee the accuracy of the fusion result. Therefore, this step employs a mass-weighted fusion method based on measurement uncertainty. Preferably, the weight values ​​of the mass-weighted fusion are proportional to the reciprocal of the measurement uncertainty of each data source. Specifically, let the first... The data source is located in the spatial position and time Water level observation value Its measurement uncertainty is Then merge water level value Calculate using the following formula: ,in: The water level value after merging is in meters (m). For the first Data sources in location and time The observed values ​​are in meters (m). For the first Quality weights for each data source, dimensionless; In this embodiment, the total number of data sources participating in the fusion is [number]. ; Subscript Indicates the data source number. Corresponding automatic monitoring station, Corresponding to manual measurement, Corresponding remote sensing inversion; subscript This indicates the result of fusion.

[0018] The quality weight Calculate using the following formula: ,in: For the first The measurement uncertainty of each data source is expressed in meters (m). This uncertainty is determined based on the sensor's nominal accuracy and historical error statistics. Preferably, the measurement uncertainty of automatic monitoring stations ranges from 0.01m to 0.05m, the uncertainty of manual measurement ranges from 0.05m to 0.15m, and the uncertainty of remote sensing inversion ranges from 0.5m to 2.0m. Through the aforementioned weighting formula, data sources with lower measurement uncertainty receive a larger weight, thus dominating the fusion results and effectively improving the overall accuracy and reliability of the fused water level dataset.

[0019] After time-series alignment and quality-weighted fusion processing, a fused water level dataset is obtained. This dataset has a unified time reference in the time dimension and retains the original observation location information of each data source in the spatial dimension, providing input data for subsequent spatial interpolation steps.

[0020] In one specific embodiment of the present invention, a groundwater monitoring area in the North China Plain is used as an example for illustration. This area covers approximately 2500 km². 2The system comprises 32 automatic monitoring stations and 68 manual measurement stations, and can simultaneously acquire GRACE satellite inversion data for the region. The automatic monitoring stations use pressure-type water level sensors with a sampling interval of 30 minutes and a measurement accuracy of ±0.02m. Manual measurements use electrical water level gauges, measured weekly, with a measurement accuracy of ±0.10m. The GRACE inversion data has a spatial resolution of approximately 300km, a temporal resolution of 1 month, and an uncertainty of approximately 1.5m. Following the method described in this step, the three types of data are first unified to a 6-hour time base. Automatic monitoring station data is downsampled using the arithmetic mean of 12 consecutive observations, while manual measurement and GRACE data are time-encrypted using cubic spline interpolation. During quality-weighted fusion, the weight of automatic monitoring station data is approximately 0.75, the weight of manual measurement data is approximately 0.20, and the weight of GRACE data is approximately 0.05. The fused water level dataset contains water level time series from approximately 2800 spatial sampling points.

[0021] Preferably, the data acquisition module also includes a data quality control function. For data from automatic monitoring stations, a sliding window midpoint filtering method is used to remove obvious outliers, with the window width preferably being 24 sampling points (corresponding to 12 hours). For manually measured data, recording errors are identified through spatial consistency checks with data from adjacent stations. For remote sensing inversion data, data reliability is improved based on signal leakage correction and destriping. After quality control, the data enters a time-series alignment and fusion process to ensure the overall quality of the fused water level dataset.

[0022] Step S2: Hybrid spatial interpolation step combining cokriging and machine learning.

[0023] In one embodiment of the present invention, the hybrid spatial interpolation step combining cokriging and machine learning is a core step in reconstructing the spatial distribution of regional groundwater levels. Because the spatial distribution of monitoring stations is sparse and uneven, spatial interpolation methods are needed to estimate the water level in areas without monitoring stations. Traditional spatial interpolation methods, such as inverse distance weighted interpolation and spline interpolation, fail to fully consider the spatial autocorrelation of groundwater levels and the influence of auxiliary variables, resulting in limited interpolation accuracy. This step innovatively combines cokriging interpolation with gradient boosting regression, fully utilizing auxiliary variable information and correcting interpolation residuals through a machine learning model to achieve high-precision spatial interpolation.

[0024] First, auxiliary variable data related to the spatial distribution of groundwater level are acquired. In one embodiment of the present invention, the auxiliary variable data includes digital elevation model (DEM) data, land use type data, and aquifer hydraulic conductivity data. The DEM data reflects topographic relief characteristics; groundwater level is generally positively correlated with topographic elevation, with areas at higher elevations having relatively shallower groundwater depths. Preferably, the spatial resolution of the DEM data is 30m × 30m or higher, and the data source can be SRTM, ASTER GDEM, or domestic high-resolution satellite data products. Land use type data reflects surface cover and human activity intensity; different land use types correspond to different groundwater extraction intensities and recharge conditions. For example, agricultural irrigation areas typically have higher extraction volumes, while urban construction areas experience reduced recharge due to surface hardening. Aquifer hydraulic conductivity data reflects the hydraulic conductivity of the aquifer; areas with high hydraulic conductivity have strong groundwater flow, and water level changes are more sensitive to extraction activities.

[0025] A cokriging interpolation model was constructed based on auxiliary variable data and a fused water level dataset. Cokriging is a multivariate extension of the ordinary kriging method, capable of simultaneously utilizing the spatial correlation of the main variable (groundwater level) and auxiliary variables for interpolation. Let the point to be estimated be... water level The cokriging estimation formula is: ,in: Points to be estimated The cokriging interpolation results are in m. For the first The measured water level values ​​at each water level monitoring point are in meters (m). ; For the first The observed values ​​of each auxiliary variable observation point are dimensionless or have corresponding physical units; For the first The weighting coefficients of the observation points of each main variable are dimensionless. For the first The weighting coefficients of the auxiliary variable observation points are dimensionless. The number of sample points of the main variable participating in the interpolation; The number of sample points is an auxiliary variable used in interpolation.

[0026] Weighting coefficient and Solving the Cokrigin equations revealed that the equations involve the semivariogram of the main variable. Auxiliary variable semivariogram and the cross-half-variogram of the main and auxiliary variables. Preferably, the semivariogram model is fitted using a spherical model or an exponential model. The model parameters include nugget value, sill value, and range, which are determined based on measured data using the least squares method or weighted least squares method. In one embodiment of the invention, the range of the main variable semivariogram is approximately 15 km to 30 km, reflecting the spatial correlation scale of the groundwater level.

[0027] While Cokriging interpolation can improve interpolation accuracy by utilizing auxiliary variable information, it is essentially still a linear interpolation method, with limited ability to characterize complex nonlinear spatial variations. Therefore, this step further utilizes a gradient boosting regression model to correct the residuals of the initial Cokriging interpolation results. Preferably, the input features of the gradient boosting regression model include spatial coordinates (longitude, latitude), auxiliary variable values ​​(elevation, land use type code, hydraulic conductivity), and Cokriging interpolation variance. Input feature vector It can be represented as: ,in: These are longitude coordinates, in degrees (°). These are latitude coordinates, in degrees (°). These are values ​​from the digital elevation model, in meters (m). The code represents the land use type and is a dimensionless integer. The aquifer's hydraulic conductivity is expressed in meters (m). 2 / d; The variance of the cokriging interpolation is expressed in m. 2 .

[0028] Gradient boosting regression models iteratively fit the residuals of preceding models, ultimately outputting corrected residual values. The model training samples are paired data of cokriging interpolation residuals and corresponding location features, where the residual is defined as the difference between the measured value of the monitoring point and the predicted value of the cokriging cross-validation. Preferably, the hyperparameters of the gradient boosting regression model are set as follows: learning rate of 0.05 to 0.15, maximum tree depth of 4 to 8, subsample ratio of 0.7 to 0.9, and number of iterations of 100 to 500.

[0029] Finally, the preliminary cokriging interpolation results are superimposed with the residual correction results to obtain the water level spatial interpolation results: ,in: The final result of the hybrid interpolation is in meters. The results are preliminary interpolation results for cokriging, in meters (m). The value represents the correction for the gradient boosting regression residual, expressed in meters. Through the aforementioned hybrid interpolation method, cokriging captures the large-scale spatial structure of the water level field, while gradient boosting regression corrects for local nonlinear variations. The two methods work synergistically to achieve high-precision spatial interpolation. In a verification experiment of one embodiment of this invention, the root mean square error (RMSE) of the cross-validation using the hybrid interpolation method was reduced by 18% to 25% compared to using cokriging alone, and by 12% to 18% compared to using gradient boosting regression alone.

[0030] In a specific embodiment of the present invention, the spatial interpolation process is further illustrated using the North China Plain monitoring area as an example. Regarding auxiliary variable data, the DEM data uses the SRTM 30m data product, with regional elevation ranging from 15m to 85m. Land use type data uses the GlobeLand30 product, which mainly includes three types of land within the region: cultivated land, urban land, and rural settlements, with cultivated land accounting for approximately 72% of the area. Aquifer hydraulic conductivity data is derived from regional hydrogeological survey results, with values ​​ranging from 50m. 2 / d to 300m 2 / d, with an average of approximately 150m 2 / d. In the construction of the cokriging interpolation model, the semivariogram of the main variable is fitted using a spherical model, and the nugget value is 0.15m. 2 The sill value is 2.8m. 2 The range of the cross-variogram of the main variable and elevation is 22 km, indicating that the spatial correlation scale of the groundwater level is approximately 22 km. The range of the cross-variogram of the main variable and elevation is 18 km, and the cross-correlation coefficient is 0.65, indicating that topography has a strong control effect on the distribution of groundwater level. The gradient boosting regression model was optimized for hyperparameters using 5-fold cross-validation. The optimal parameter combination was a learning rate of 0.08, a maximum tree depth of 6, a subsample ratio of 0.8, and 280 iterations. The final cross-validation RMSE of the hybrid interpolation result was 0.42 m, which is 20.8% lower than the 0.53 m of cokriging alone, verifying the effectiveness of the hybrid interpolation method.

[0031] Furthermore, the spatial interpolation module also supports interpolation uncertainty assessment. The cokriging method naturally provides interpolation variance as a measure of uncertainty; regions with large interpolation variance indicate lower reliability of the estimation results, typically corresponding to areas with sparse monitoring stations or insufficient coverage of auxiliary variables. This uncertainty information can be passed to subsequent contour surface generation and early warning visualization stages, prompting users to pay attention to areas of high uncertainty in the visualization output through dashed lines or changes in transparency.

[0032] Step S3: Isosurface generation and optimization steps.

[0033] In one embodiment of the present invention, the core objective of the isosurface generation and optimization step is to convert the water level spatial interpolation results obtained in step S2 into intuitive isoline and isosurface graphic representations, so that users can quickly understand the spatial distribution pattern of regional groundwater levels.

[0034] Water level contour data is generated using an improved contour extraction algorithm based on water level spatial interpolation results. Preferably, the improved contour extraction algorithm is a contour extraction algorithm based on Marching Squares. The Marching Squares algorithm divides a regular grid into several cells and determines the direction of the contour line within that cell based on the height relationship between the four vertices of each cell and the contour line threshold. This invention improves the standard Marching Squares algorithm by introducing a saddle point processing mechanism and boundary continuity constraints. For cells with saddle points (i.e., water level values ​​at opposite diagonal vertices are similar but the other opposite diagonal vertex has a large difference), the contour line direction is determined by calculating the water level value at the center of the cell, avoiding breaks or intersections in the contour lines. Boundary continuity constraints ensure a smooth transition of contour lines between adjacent cells, avoiding jagged boundaries.

[0035] In one embodiment of the present invention, the contour interval is adaptively determined based on the magnitude of regional water level changes. Let the regional water level range be... Preferred contour intervals Calculate using the following formula:

[0036] ,in: These are contour intervals, in meters (m). The regional water level range is equal to the difference between the highest and lowest water levels, and is expressed in meters (m). The desired number of contour lines is preferably between 8 and 15. For example, if the water level in the area varies from 10m to 30m, ,Pick ,but That is, draw a contour line every 2m.

[0037] The extracted raw contour data often contains too many redundant nodes and local jitter, requiring curve simplification and smoothing to improve visualization. Preferably, the curve simplification employs the Douglas-Peucker algorithm. This algorithm recursively removes intermediate nodes that contribute little to the curve's shape, reducing the number of nodes while preserving the overall curve shape. Simplification tolerance parameters are also considered. The value is set according to the map display scale and accuracy requirements, preferably ranging from 5m to 50m.

[0038] After curve simplification, the contour lines are smoothed using the Catmull-Rom spline interpolation algorithm. The Catmull-Rom spline is a cubic spline curve passing through control points, capable of generating smooth, continuous curves. Preferably, 4 to 10 smoothing points are inserted between adjacent simplified nodes to give the contour lines a natural and smooth curve shape, avoiding a broken line appearance.

[0039] After contour line extraction, curve simplification, and smoothing, the processed contour line data is rendered to generate a water level isosurface map. The contour lines use gradient color coding to represent different water level intervals. Preferably, warm colors (such as red and orange) are used for low water level areas, and cool colors (such as blue and green) are used for high water level areas, visually reflecting the distribution of water levels. The isosurface map supports overlay display on a base map, which can be satellite imagery, topographic maps, or administrative division maps.

[0040] In a specific embodiment of the present invention, the process of isosurface generation is further illustrated using the North China Plain monitoring area as an example. The grid resolution of the water level spatial interpolation result is 200m × 200m, with approximately 62,500 grid nodes. The regional water level ranges from 12.5m to 38.2m, with a range of 25.7m. Ten isosurfaces are set with an interval of 2.57m, rounded down to 2.5m. An improved Marching Squares algorithm is used for isosurface extraction. For the saddle point area in the southwest of the region, the isosurface direction is determined by calculating the water level value at the center of each cell. The generated original isosurfaces contain approximately 28,000 nodes. Curve simplification uses the Douglas-Peucker algorithm with a simplification tolerance of 20m. After simplification, the number of nodes is reduced to approximately 4,500, with a retention rate of approximately 16%. Smoothing is achieved using Catmull-Rom splines, inserting six smoothing points between adjacent simplified nodes, resulting in a smooth and continuous curve shape for the final isosurfaces. The color coding uses a warm-cool color gradient scheme. The water level range of 12.5m to 20m uses a red-to-orange gradient, the range of 20m to 30m uses a yellow-to-green gradient, and the range of 30m to 38.2m uses a green-to-blue gradient. The generated isosurface map is overlaid on the satellite imagery base map, clearly showing the spatial distribution pattern of the regional groundwater level. A distinct low-water-level area is formed in the central part of the region, corresponding to the location of the historical over-extraction funnel.

[0041] Preferably, the isosurface generation module also supports multiple visualization style configurations. In addition to isoline display, it also supports an isosurface fill display mode, which fills adjacent isolines with gradient colors to form a continuous color distribution map. It supports isoline annotation, allowing water level values ​​to be annotated at appropriate locations on the isolines, with annotation density and font size adaptively adjusted according to the map scale. It supports map symbol overlay, allowing monitoring station locations, administrative boundaries, water system distribution, and other elements to be overlaid on the isosurface map, providing rich spatial reference information.

[0042] Step S4: Timing animation compositing steps.

[0043] In one embodiment of the present invention, the goal of the time-series animation synthesis step is to serialize and synthesize the water level isosurface maps of multiple time points into a dynamic evolution animation, which intuitively presents the seasonal fluctuations and interannual variation trends of groundwater levels, making it easier for managers to grasp the dynamic change patterns of water levels from a time perspective.

[0044] First, the water level isosurface maps at multiple time points are serialized to obtain a multi-temporal water level distribution map sequence. Assume there are a total of [number missing] monitoring points within the monitoring period. Given a time point and a corresponding isosurface map of water level, the multi-temporal water level distribution map sequence can be represented as follows: ,in Indicates the first Water level isosurface maps at specific time points. Preferably, the frequency of time point selection is determined according to the analysis requirements; short-term analysis can use daily or weekly scales, while long-term analysis can use monthly or quarterly scales. In one embodiment of the present invention, monthly scale time points are used to generate 36 isosurface maps from water level data of a certain region for 3 consecutive years (36 months).

[0045] Because there is a time interval between adjacent time points, directly playing the isosurface map sequence sequentially may create a visual jumpiness, affecting the smoothness of the animation. Therefore, inter-frame interpolation processing is performed on the multi-temporal water level distribution map sequence. Preferably, the inter-frame interpolation processing uses linear interpolation or spline interpolation methods to generate transition frames between adjacent time points. Let two adjacent time points... and Insertion is required between them The nth transition frame, then the nth The water level field for each transition frame is calculated using the following formula: ,in: For the first A transition frame at position The water level interpolation results are in meters. Time node The spatial interpolation results of the water level are in meters. Time node The spatial interpolation results of the water level are in meters. These are linear interpolation coefficients. By using inter-frame interpolation, a smooth transition is achieved between adjacent isosurface maps, eliminating visual jarring.

[0046] After inter-frame interpolation is completed, the original isosurface map and transition frames are synthesized into temporal animation data in chronological order. Preferably, the frame rate of the temporal animation data is 10fps to 30fps, that is, 10 to 30 frames of images are played per second. In one embodiment of the present invention, the frame rate is set to 15fps, and 14 transition frames are inserted between adjacent months, so that each month corresponds to an animation duration of 1 second. The animation encoding format can be GIF, MP4, or WebM, supporting smooth playback on web or mobile devices.

[0047] The timing animation synthesis also supports various visualization enhancement features. Preferably, time labels are overlaid on the animation interface to indicate the date and time corresponding to the current frame in real time. Playback control functions are supported, including play / pause, fast forward / rewind, and jumping to a specified time node. A loop playback mode is supported, automatically jumping from the last frame to the first frame to continue playback, making it easy to observe periodic changes.

[0048] In a specific embodiment of the present invention, the time-series animation synthesis process is further illustrated using the North China Plain monitoring area as an example. Water level data from January 2021 to December 2023 (36 months in total) are selected to generate 36 monthly-scale water level isosurface maps. Fourteen transition frames are inserted between every two adjacent months, generated using a linear interpolation method. Adding the original 36 frames, the total number of frames is 36 + 35 × 14 = 526 frames. The animation frame rate is set to 15fps, with a total duration of approximately 35 seconds. The animation is encoded in MP4 format, with a video resolution of 1920 × 1080 pixels, a bitrate of 8Mbps, and a file size of approximately 35MB. A time label is overlaid in the upper left corner of the animation, in the format YYYY year MM month, with a font size of 24pt and a white font color with a black border to ensure readability under different background colors. A color legend is overlaid in the lower right corner of the animation, indicating the correspondence between water level values ​​and colors.

[0049] The time-series animation visually illustrates the seasonal fluctuations in groundwater levels in the region. During the peak irrigation season from March to May, the overall water level declines, and the warm-colored areas (low water levels) on the isosurface map expand significantly. From July to September, the rainy season replenishes the area, causing the water level to gradually rise, and the cool-colored areas (high water levels) increase. Simultaneously, the over-extraction funnel in the central part of the region shows a slow expansion trend over the three years, with the water level at the center of the funnel decreasing from 15.2m at the beginning of 2021 to 13.8m at the end of 2023, a cumulative decrease of 1.4m. The time-series animation presents this gradual process dynamically, making it more intuitive and impactful than static charts.

[0050] Preferably, the time-series animation module also supports comparative playback. Time-series animations from different years within the same region can be displayed side-by-side, facilitating comparison of interannual variations. Time-series animations from different regions can also be displayed side-by-side, facilitating comparison of spatial differences. During comparative playback, multiple animation windows employ a synchronized playback mechanism to ensure timeline consistency, enabling users to perform cross-temporal and spatial comparative analysis.

[0051] Step S5: Anomaly identification and classification steps.

[0052] In one embodiment of the present invention, the anomaly identification and classification step is the core analysis link for realizing groundwater level anomaly early warning. By calculating the spatiotemporal gradient characteristics of water level changes and performing anomaly discrimination, typical anomaly types such as over-extraction funnel expansion area, water level abnormal decline area, and groundwater recharge enhancement area are identified.

[0053] First, the spatiotemporal gradient characteristics of water level changes are calculated based on a multi-temporal water level distribution map sequence. Preferably, the spatiotemporal gradient characteristics include a spatial gradient component and a temporal gradient component. The spatial gradient reflects the spatial rate of change of the water level field and is defined as the directional derivative of the water level in the horizontal direction. Let the location... In time water level value The spatial gradient components are calculated using the following formula: ,in: The magnitude of the spatial gradient component is expressed in m / km. This is the partial derivative of the water level in the east-west direction, with units of m / km; This represents the partial derivative of the water level in the north-south direction, expressed in m / km. The partial derivative is calculated using the central difference method based on the water level values ​​of adjacent grid points.

[0054] The water level time gradient reflects the rate of change of water level over time and is defined as the derivative of water level in the time direction: ,in: This represents the time gradient component, with units of m / month; The time interval between adjacent time nodes is preferably one month.

[0055] In one embodiment of the present invention, a spatiotemporal gradient tensor is further constructed. Comprehensive characterization of the spatiotemporal variation characteristics of the water level field: An anomaly index is calculated based on the spatiotemporal gradient features and compared with a preset threshold for anomaly detection. Anomaly Index The degree of anomalousness of the current water level change relative to historical statistical characteristics is calculated using the following formula: ,in: An anomaly index, dimensionless; This represents the time gradient value at the current moment, in meters per month. This represents the historical average time gradient for the same period, expressed in m / month. The standard deviation of the time gradient for the same period in history is expressed in m / month. The anomaly index is essentially a standardized time gradient deviation, reflecting the degree to which the current rate of water level change deviates from the historical normal fluctuation range.

[0056] Let the preset anomaly detection threshold be... Preferably The value range is from 2.0 to 3.0. When If the location is determined to be in an abnormal state, further classification of the abnormality type is required.

[0057] The abnormal regions are classified according to the anomaly discrimination results to obtain the anomaly region identification results. Preferably, the anomaly region classification includes three main types: over-extraction funnel expansion region, abnormal water level drop region, and enhanced groundwater recharge region. The discrimination conditions for each type are as follows:

[0058] The criteria for identifying the expansion zone of the over-extraction funnel are a continuous decline in water level and a spatial gradient pointing towards the center of the zone. Specifically, when the time gradient of multiple consecutive time points (preferably not less than 3 months) If the spatial gradient direction exhibits a centripetal distribution characteristic (i.e., the gradient vector points from the periphery to the center of the low water level), then it is determined to be an over-extraction funnel expansion area. The centripetal distribution characteristic is quantified by calculating the gradient divergence; when the divergence... hour( The gradient is assumed to be centripetal distribution, with a preferred value of 0.1 / km as the divergence threshold.

[0059] The criterion for identifying areas of abnormal water level decline is that the rate of water level decline exceeds a preset multiple of the historical average for the same period. Specifically, when... hour( The factor (preferably a value of 2 to 3) is used to classify areas as having abnormally low water levels. These areas experience water level drops significantly exceeding the normal seasonal fluctuation range, possibly caused by factors such as localized over-extraction, reduced recharge, or changes in geological conditions.

[0060] The criterion for identifying areas of enhanced groundwater recharge is that the rate of water level rise exceeds a preset multiple of the historical average for the same period. Specifically, when At that time, it was identified as an area of ​​enhanced groundwater recharge. The rise in water level in this type of area significantly exceeds the normal seasonal fluctuation range, which may be caused by factors such as abnormally high precipitation, irrigation recharge, and increased surface water recharge.

[0061] Through the above anomaly identification and classification process, the results of anomaly area identification are obtained, including the spatial range, anomaly type and anomaly intensity of each anomaly area, providing data support for subsequent early warning visualization.

[0062] In a specific embodiment of the present invention, the anomaly identification process is further illustrated using the North China Plain monitoring area as an example. Based on the water level time series for 36 months from January 2021 to December 2023, the temporal gradient sequence of each pixel is calculated. The historical mean and standard deviation of the temporal gradient for the same period are calculated based on historical data from 2015 to 2020 (6 years). The anomaly discrimination threshold is set to 2.5, meaning that an anomaly index exceeding 2.5 is considered an abnormal state. The anomaly identification results for December 2023 show that three anomaly areas were identified within the region: the first is located in the central part of the region, with an area of ​​approximately 85 km². 2 The first location was identified as an over-extraction funnel expansion area. This area had a negative temporal gradient for 12 consecutive months, with an average decrease rate of 0.12 m / month. The spatial gradient showed a clear centripetal distribution, with a gradient divergence of -0.18 / km. The second location was located in the northeastern part of the region, covering an area of ​​approximately 32 km². 2 The first area was identified as having an abnormally low water level. The time gradient for this area from November to December 2023 was -0.35 m / month, significantly exceeding the historical average of -0.08 m / month and three standard deviations for the same period. This may be related to increased water consumption in newly built industrial parks in the area. The third area is located in the southwest of the region, covering approximately 18 km². 2 The area was identified as a region with enhanced groundwater recharge. The temporal gradient of this area from September to October 2023 was +0.28 m / month, which significantly exceeded the historical average of +0.05 m / month and three times the standard deviation, and is related to the abnormally high precipitation in the summer of 2023.

[0063] Preferably, the anomaly identification module also supports anomaly evolution tracking. For persistent anomaly areas, the system automatically records the time of their first appearance, the history of area changes, and the trend of anomaly intensity changes, generating an anomaly evolution report. This report can serve as an important reference for groundwater resource management decisions, helping managers determine whether the anomaly is intensifying or mitigating, thereby enabling them to take appropriate control measures.

[0064] Step S6: Visualization of the early warning heatmap.

[0065] In one embodiment of the present invention, the early warning heat map visualization step is the final output of the entire method. It presents the anomaly identification results and water level decline rate information in the form of a heat map, and marks the distribution of areas exceeding the warning water level, providing visualization support for groundwater resource management decisions.

[0066] A warning heatmap is generated based on the anomaly area identification results and the water level drop rate calculated based on the spatiotemporal gradient features. The heatmap uses a two-dimensional color matrix to represent spatial distribution information, with the color intensity of each pixel reflecting the warning level at that location. Let the heatmap pixels... The corresponding spatial location is Heat map intensity value Calculate using the following formula: ,in: The value represents the pixel intensity of the heatmap, ranging from 0 to 1, and is dimensionless. An anomaly index, dimensionless; The rate of water level decline is expressed in meters per month. This represents the maximum rate of water level decline within the region, expressed in m / month. and For the weighting coefficients, satisfying Preferably , .

[0067] The intensity values ​​of the heatmap are converted into visible colors through color mapping. Preferably, a multi-threshold progressive rendering strategy is used to divide the intensity value range into several warning levels, each corresponding to a different color: Corresponds to green (normal state); Corresponding to yellow (attention status); Corresponding to orange (warning status); Corresponding to red (alert status). Multi-threshold progressive rendering ensures a natural transition between alert levels, avoiding abrupt color changes at the boundaries.

[0068] The distribution of areas exceeding the warning water level is marked on the warning heat map. Let the regional warning water level value be... When the real-time water level at a certain location When the water level at a location is below the warning level, it is determined that the water level is below the warning level. The warning level is usually determined comprehensively based on the sustainable exploitation capacity of the aquifer, ecological water demand requirements, and land subsidence control requirements, and is preferably verified and issued by the water resources authority. Areas exceeding the warning level (i.e., below the warning level value) are marked with special symbols, such as filled shading lines, thickened boundaries, or flashing effects, to remind managers to pay close attention.

[0069] The early warning heatmap also supports various interactive functions. Preferably, it supports layer overlay display, allowing the early warning heatmap to be superimposed on satellite imagery, administrative divisions, or hydrogeological maps. It supports a click-to-query function; clicking anywhere on the heatmap displays detailed information for that location, including the current water level, water level change trend, anomaly type, and warning level. It also supports a time slider function; dragging the slider allows viewing early warning heatmaps at any historical point in time, tracing the evolution of anomalies.

[0070] In a specific embodiment of the present invention, the process of generating the early warning heatmap is further illustrated using the North China Plain monitoring area as an example. Based on the anomaly identification results and water level decline rate data from December 2023, the heatmap intensity value is calculated according to the above formula. The distribution of heatmap intensity values ​​within the region is as follows: approximately 68% of the region has an intensity value less than 0.3, corresponding to a green normal state; approximately 22% of the region has an intensity value between 0.3 and 0.5, corresponding to a yellow concern state; approximately 7% of the region has an intensity value between 0.5 and 0.7, corresponding to an orange warning state; and approximately 3% of the region has an intensity value greater than or equal to 0.7, corresponding to a red alert state. The red alert areas are mainly distributed in the core area of ​​the over-extraction funnel in the central part of the region and the abnormal decline area in the northeast. The warning water level in this area is set at 8.0m by the water resources department (based on the 1985 National Elevation Datum). When the real-time water level is lower than 8.0m, it is determined to be above the warning water level. Monitoring results in December 2023 showed that water levels in about 12% of the area were below the warning level, mainly located in the central area of ​​the over-extraction funnel, which was specially marked with diagonal shading on the early warning heat map.

[0071] The heatmap's color rendering employs bilinear interpolation to achieve smooth spatial transitions, avoiding noticeable color jumps at pixel boundaries. The heatmap resolution adaptively matches the base map's display scale, using a higher resolution at large scales to reveal local details and a lower resolution at small scales to improve rendering performance. The heatmap supports transparency adjustment, with a default transparency of 0.7. Users can adjust the transparency as needed to better observe the base map information.

[0072] Preferably, the early warning visualization module also supports automatic early warning report generation. Based on the early warning heatmap and anomaly identification results, the system automatically generates early warning reports containing text descriptions, data tables, and charts. The report content includes the overall regional water level situation, area statistics for each early warning level, detailed information on abnormal areas, comparative analysis with historical data for the same period, trend forecasts, and management recommendations. The report can be exported to PDF or Word format for easy archiving and reporting by managers.

[0073] Furthermore, this method supports a closed-loop feedback optimization mechanism. The anomaly identification result in step S5 can be fed back to step S2 to adjust the spatial interpolation model parameters. When an anomaly region is identified, the weight of the auxiliary variable in that region can be increased or local anisotropic parameters can be introduced to improve the interpolation accuracy of the anomaly region. The early warning result in step S6 can be fed back to step S1 to dynamically adjust the data acquisition frequency. When the early warning level of a certain region increases, the automatic monitoring station in that region can be automatically triggered to increase the sampling frequency (e.g., from once per hour to once every 15 minutes) to improve the ability to capture abnormal changes. This closed-loop feedback mechanism enables the system to adaptively optimize operating parameters based on monitoring results, continuously improving the quality of monitoring and analysis.

[0074] Through the coordinated processing of the above six steps, this invention realizes the entire process of spatiotemporal dynamic monitoring of groundwater level from multi-source data acquisition, spatial interpolation analysis, visualization presentation to anomaly early warning, providing effective technical support for the scientific management and dynamic scheduling of groundwater resources.

[0075] Reference Figure 2 As shown in the figure, this embodiment of the invention also provides a groundwater level spatiotemporal dynamic monitoring and visualization analysis system. This system corresponds to the aforementioned method embodiment and includes a data acquisition module 1, a spatial interpolation module 2, an isosurface generation module 3, a temporal animation module 4, an anomaly identification module 5, and an early warning visualization module 6. The functions of each module are consistent with the corresponding steps in the method embodiment, and are briefly described below.

[0076] Data acquisition module 1 is used to acquire multi-source groundwater level monitoring data in the region, perform time-series alignment processing on the multi-source groundwater level monitoring data to obtain a time-series alignment result, and perform quality-weighted fusion based on the measurement uncertainty of each data source to obtain a fused groundwater level dataset. In one embodiment, the data acquisition module 1 includes a data interface submodule, a time-series alignment submodule, and a quality-weighted fusion submodule. The data interface submodule is responsible for connecting to the automatic monitoring station data transmission system, the manual measurement data entry system, and the remote sensing data service platform to acquire multi-source raw monitoring data. The time-series alignment submodule performs time-dimensional alignment processing on data with different sampling frequencies. The quality-weighted fusion submodule performs data fusion processing according to the quality-weighted fusion formula described in the method embodiment.

[0077] Spatial interpolation module 2 is used to acquire auxiliary variable data related to the spatial distribution of groundwater level. Based on the auxiliary variable data and the fused groundwater level dataset, a cokriging interpolation model is constructed, and spatial interpolation is performed to obtain preliminary cokriging interpolation results. A gradient boosting regression model is used to correct the residuals of the preliminary cokriging interpolation results to obtain a residual correction result. The preliminary cokriging interpolation results and the residual correction result are superimposed to obtain the spatial interpolation result of the groundwater level. In one embodiment, spatial interpolation module 2 includes an auxiliary variable management submodule, a cokriging interpolation submodule, and a residual correction submodule. The auxiliary variable management submodule is responsible for storing and managing DEM data, land use data, and aquifer parameter data. The cokriging interpolation submodule implements the cokriging interpolation algorithm and supports automatic fitting of the semivariogram and calculation of interpolation weights. The residual correction submodule implements the training and prediction functions of the gradient boosting regression model.

[0078] The isosurface generation module 3 is used to generate water level isosurface data based on the water level spatial interpolation results using an improved isosurface extraction algorithm, and to perform curve simplification and smoothing processing on the water level isosurface data to generate a water level isosurface map. In one embodiment, the isosurface generation module includes an isosurface extraction submodule, a curve optimization submodule, and a rendering submodule. The isosurface extraction submodule implements the improved Marching Squares algorithm. The curve optimization submodule implements Douglas-Peucker simplification and Catmull-Rom spline smoothing. The rendering submodule renders the isosurface data into a color-coded isosurface map.

[0079] The temporal animation module 4 is used to serialize water level isosurface maps at multiple time points to obtain a multi-temporal water level distribution map sequence, perform inter-frame interpolation on the multi-temporal water level distribution map sequence, and synthesize temporal animation data. In one embodiment, the temporal animation module 4 includes a sequence management submodule, an inter-frame interpolation submodule, and an animation encoding submodule. The sequence management submodule is responsible for managing the multi-temporal isosurface map sequence. The inter-frame interpolation submodule generates transition frames between adjacent time points. The animation encoding submodule encodes the image sequence into a video format for output.

[0080] Anomaly identification module 5 is used to calculate the spatiotemporal gradient features of water level changes based on the multi-temporal water level distribution map sequence, calculate anomaly indices based on the spatiotemporal gradient features and compare them with preset thresholds to perform anomaly discrimination, and classify the anomaly regions according to the anomaly discrimination results to obtain anomaly region identification results. In one embodiment, the anomaly identification module 5 includes a gradient calculation submodule, an anomaly discrimination submodule, and a classification submodule. The gradient calculation submodule implements the calculation of the spatiotemporal gradient tensor. The anomaly discrimination submodule implements threshold-based anomaly index discrimination. The classification submodule classifies the anomaly regions into three categories according to discrimination rules: over-extraction funnel, abnormal decline, and enhanced recharge.

[0081] The early warning visualization module 6 is used to generate an early warning heatmap based on the abnormal area identification results and the water level drop rate calculated based on the spatiotemporal gradient features, and to mark the distribution of areas exceeding the warning water level on the early warning heatmap. In one embodiment, the early warning visualization module 6 includes a heatmap generation submodule, a layer overlay submodule, and an interaction submodule. The heatmap generation submodule generates heatmaps with multi-threshold progressive rendering. The layer overlay submodule supports the overlay display of heatmaps with multiple base maps. The interaction submodule provides interactive functions such as click query, time sliding, and layer control.

[0082] The modules of the aforementioned system can be implemented through software, hardware, or a combination of both. In one embodiment, the system is deployed on a cloud server, employing a microservice architecture to enable independent deployment and elastic scaling of each module. In another embodiment, the system is deployed on an edge computing device to achieve on-site processing and rapid response of monitoring data.

[0083] The embodiments of the present invention are not limited to the specific embodiments described above. Those skilled in the art can make various equivalent changes or substitutions based on the technical solutions of the present invention, and all such changes or substitutions should be included within the protection scope of the present invention.

Claims

1. A method for spatiotemporal dynamic monitoring and visualization analysis of groundwater levels, characterized in that, Includes the following steps: The multi-source data fusion preprocessing step involves acquiring multi-source groundwater level monitoring data for the region, performing time-series alignment processing on the multi-source groundwater level monitoring data to obtain a time-series alignment processing result, and performing quality-weighted fusion on the time-series alignment processing result based on the measurement uncertainty of each data source to obtain a fused groundwater level dataset; the weight value of the quality-weighted fusion is proportional to the reciprocal of the measurement uncertainty of each data source, and the measurement uncertainty is determined based on the sensor accuracy and historical error statistics of each data source; The process involves a hybrid spatial interpolation step combining cokriging and machine learning. First, auxiliary variable data related to the spatial distribution of groundwater level are acquired. This auxiliary variable data includes digital elevation model data, land use type data, and aquifer hydraulic conductivity data. Based on this auxiliary variable data and the fused groundwater level dataset, a cokriging interpolation model is constructed, and spatial interpolation is performed to obtain preliminary cokriging interpolation results. The cross-semivariogram of the auxiliary and main variables in the cokriging interpolation model is determined by fitting measured data. Then, a gradient boosting regression model is used to correct the residuals of the preliminary cokriging interpolation results to obtain a corrected residual result. The input features of the gradient boosting regression model include spatial coordinates, auxiliary variable values, and cokriging interpolation variance. The training samples of the gradient boosting regression model are paired data of cokriging interpolation residuals and corresponding location features. Finally, the preliminary cokriging interpolation results and the corrected residual result are superimposed to obtain the spatial interpolation result of the groundwater level. The isosurface generation and optimization steps involve generating water level isosurface data based on the water level spatial interpolation results using a Marching Squares-based isosurface extraction algorithm. This algorithm incorporates a saddle point processing mechanism and boundary continuity constraints. The Douglas-Peucker algorithm is used to simplify the curves in the water level isosurface data, and the Catmull-Rom spline interpolation algorithm is used for smoothing to generate a water level isosurface map. The temporal animation synthesis step involves serializing the water level isosurface maps at multiple time points to obtain a multi-temporal water level distribution map sequence, performing inter-frame interpolation on the multi-temporal water level distribution map sequence, and synthesizing the temporal animation data. The anomaly identification and classification steps involve calculating the spatiotemporal gradient features of water level changes based on the multi-temporal water level distribution map sequence. These features include spatial and temporal gradient components of the water level. An anomaly index is calculated based on these features and compared with a preset threshold for anomaly identification. The anomaly index is calculated based on the deviation of the magnitude of the spatiotemporal gradient features from the statistical mean. Anomaly areas are classified according to the anomaly identification results to obtain anomaly area identification results. The anomaly area classification includes over-extraction funnel expansion areas, abnormal water level decline areas, and enhanced groundwater recharge areas. The identification condition for an over-extraction funnel expansion area is a continuous decline in water level with the spatial gradient pointing towards the center of the area. The identification condition for an abnormal water level decline area is a water level decline rate exceeding a preset multiple of the historical average for the same period. The identification condition for an enhanced groundwater recharge area is a water level rise rate exceeding a preset multiple of the historical average for the same period. The early warning heatmap visualization step involves generating an early warning heatmap based on the abnormal area identification results and the water level drop rate calculated based on the spatiotemporal gradient features. The intensity value of the early warning heatmap is calculated based on a weighted combination of the anomaly index and the normalized water level drop rate. A multi-threshold progressive rendering strategy is used to divide the intensity value range into four early warning levels: normal state, attention state, early warning state, and alarm state. The distribution of areas exceeding the warning water level is marked on the early warning heatmap.

2. The method for spatiotemporal dynamic monitoring and visualization analysis of groundwater level according to claim 1, characterized in that, The multi-source water level monitoring data includes real-time water level data from automatic monitoring stations, manually measured water level data, and groundwater storage change data retrieved by remote sensing; the time-series alignment processing includes unifying data from different sampling frequencies to the same time reference, wherein the time interval of the time reference is from 1 hour to 24 hours.

3. The method for spatiotemporal dynamic monitoring and visualization analysis of groundwater level according to claim 1, characterized in that, The inter-frame interpolation process uses linear interpolation or spline interpolation methods to generate transition frames between adjacent time nodes, and the frame rate of the time-series animation data is 10fps to 30fps.

4. A groundwater level spatiotemporal dynamic monitoring and visualization analysis system, used to implement the groundwater level spatiotemporal dynamic monitoring and visualization analysis method according to any one of claims 1-3, characterized in that, include: The data acquisition module is used to acquire multi-source groundwater level monitoring data in the region, perform time-series alignment processing on the multi-source groundwater level monitoring data to obtain time-series alignment processing results, and perform quality-weighted fusion on the time-series alignment processing results based on the measurement uncertainty of each data source to obtain a fused groundwater level dataset. The weight value of the quality-weighted fusion is proportional to the reciprocal of the measurement uncertainty of each data source. The spatial interpolation module is used to acquire auxiliary variable data related to the spatial distribution of groundwater level. The auxiliary variable data includes digital elevation model data, land use type data, and aquifer hydraulic conductivity data. Based on the auxiliary variable data and the fused groundwater level dataset, a cokriging interpolation model is constructed and spatial interpolation is performed to obtain preliminary cokriging interpolation results. The cross-semivariogram function of the auxiliary variables and the main variables in the cokriging interpolation model is determined by fitting the measured data. The residuals of the preliminary cokriging interpolation results are corrected using a gradient boosting regression model to obtain a residual correction result. The input features of the gradient boosting regression model include spatial coordinates, auxiliary variable values, and cokriging interpolation variance. The preliminary cokriging interpolation results and the residual correction results are superimposed to obtain the spatial interpolation result of groundwater level. The isosurface generation module is used to generate water level isosurface data based on the water level spatial interpolation results using a Marching Squares-based isosurface extraction algorithm. The Marching Squares-based isosurface extraction algorithm introduces a saddle point processing mechanism and boundary continuity constraints. The Douglas-Peucker algorithm is used to simplify the curves of the water level isosurface data, and the Catmull-Rom spline interpolation algorithm is used for smoothing to generate a water level isosurface map. The temporal animation module is used to serialize the water level isosurface maps at multiple time points to obtain a multi-temporal water level distribution map sequence, perform inter-frame interpolation on the multi-temporal water level distribution map sequence, and synthesize temporal animation data. An anomaly identification module is used to calculate the spatiotemporal gradient features of water level changes based on the multi-temporal water level distribution map sequence. The spatiotemporal gradient features include water level spatial gradient components and water level temporal gradient components. An anomaly index is calculated based on the spatiotemporal gradient features and compared with a preset threshold to identify anomalies. Anomaly regions are classified according to the anomaly identification results to obtain anomaly region identification results. The anomaly region classification includes over-extraction funnel expansion regions, water level anomaly decline regions, and groundwater recharge enhancement regions. The early warning visualization module is used to generate an early warning heat map based on the abnormal area identification results and the water level drop rate calculated based on the spatiotemporal gradient features. The intensity value of the early warning heat map is calculated based on the weighted combination of the abnormal index and the normalized water level drop rate. A multi-threshold progressive rendering strategy is used to divide the early warning into four levels, and the distribution of areas exceeding the warning water level is marked on the early warning heat map.

Images