Regional groundwater reserve change dynamic monitoring method based on gravity satellite time series data analysis

Abstract

The invention relates to the technical field of hydrological monitoring, and discloses a regional groundwater reserve change dynamic monitoring method based on gravity satellite time series data analysis. According to the method, the multi-source data scale effect and observation noise are effectively eliminated, the physical reliability and resolution of groundwater signal extraction are improved, and the recognition precision and early warning perspectiveness of the abnormal loss trend are enhanced. A full-chain monitoring closed loop of multi-source satellite observation, high-precision physical inversion and dynamic statistical benchmark is constructed, multi-dimensional environment data is aligned, a component separation strategy based on multi-element physical mechanism constraint is established based on a water balance principle, and anomaly diagnosis is performed in combination with a self-adaptively updated long-time-sequence historical database.

Description

Technical Field

[0001] This invention relates to the field of hydrological monitoring technology, and more specifically, to a method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis. Background Technology

[0002] Groundwater, a crucial component of the terrestrial water cycle, is a vital strategic resource for maintaining global ecosystem balance, ensuring agricultural irrigation, and supporting urban water supply security. Long-term, comprehensive dynamic monitoring of regional groundwater storage changes can reveal the evolutionary patterns of water resources driven by both natural climate fluctuations and human activities, providing indispensable basic data support for the optimal allocation of regional water resources, the management of groundwater over-extraction, and the prevention of geological disasters. Especially against the backdrop of increasingly severe climate change, accurately understanding the state of groundwater abundance and deficiency has profound practical significance for maintaining wetland functions, preventing land subsidence, and ensuring sustainable socio-economic development.

[0003] However, existing groundwater monitoring technologies still face numerous technical bottlenecks in practical applications. While traditional surface well observations offer high accuracy, they are limited by sparse station distribution, high construction and maintenance costs, and insufficient spatial representativeness, making it difficult to achieve continuous area coverage over large regions. In contrast, while satellite remote sensing technology achieves wide-area coverage and can macroscopically detect changes in surface mass, its observation principle dictates that it acquires total mass changes including surface water, soil water, snowmelt, and groundwater. This mixed signal characteristic makes it difficult to directly distinguish between different water components. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis, aiming to solve the problem of groundwater independent components that are difficult to calculate from multi-dimensional remote sensing data.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0007] This invention provides a method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis, which includes the following steps:

[0008] S1. Obtain time-varying gravity field spherical harmonic data of the target area through gravity satellite, and invert the time-varying gravity field spherical harmonic data to obtain the land water storage data of the target area;

[0009] S2. Acquire multi-source environmental data of the target area, and resample the multi-source environmental data to the same spatial resolution as the gravity satellite data;

[0010] S3. Obtain the non-groundwater reserves of the target area by parsing the multi-source environmental data, and then calculate the groundwater reserve distribution data of the target area;

[0011] S4. Continuously observe the target area, construct regional groundwater time series distribution data, compare and analyze the regional groundwater time series distribution data with the groundwater dynamic evolution database of the target area, and generate a dynamic monitoring report after identifying abnormal groundwater depletion.

[0012] As a preferred embodiment of the regional groundwater storage change dynamic monitoring method based on gravity satellite time-series data analysis according to the present invention, wherein: before inverting the time-varying gravity field spherical harmonic data in step S1, a preprocessing step of the time-varying gravity field spherical harmonic data is further included, the preprocessing step specifically including:

[0013] Obtain the C20 coefficients from the satellite laser ranging data, replace the corresponding terms in the time-varying gravity field spherical harmonic data with the C20 coefficients, and correct the reference frame origin drift by introducing a geocentric correction term.

[0014] Then, polynomial fitting is performed to remove the correlation between odd and even order coefficients to eliminate strip error, and Gaussian smoothing operator is applied to suppress noise in higher order terms.

[0015] The signal leakage scale factor is determined by the signal strength attenuation ratio before and after combined filtering, and the amplitude of the filtered data is corrected using the scale factor.

[0016] As a preferred embodiment of the regional groundwater storage change dynamic monitoring method based on gravity satellite time-series data analysis described in this invention, the process of inverting the terrestrial water storage data of the target area in step S1 specifically includes: introducing the Earth load Love number to weight and correct the preprocessed spherical harmonic coefficients to eliminate the interference of Earth's elastic deformation on the gravity signal; substituting the corrected spherical harmonic coefficients into a physical conversion formula containing the Earth's average radius, Earth's average density, and water density to convert the dimensionless spherical harmonic coefficients into mass coefficients with equivalent water height physical meaning; using spherical harmonic synthesis processing, mapping the mass coefficients to the spatial domain, calculating the values ​​of each latitude and longitude grid point in the target area point by point, and then generating terrestrial water storage data reflecting the target area.

[0017] As a preferred embodiment of the regional groundwater storage dynamic monitoring method based on gravity satellite time-series data analysis described in this invention, the multi-source environmental data in step S2 includes: normalized vegetation index data reflecting the growth status of surface vegetation, surface temperature data reflecting the thermal status of the surface, and rainfall data reflecting the regional water cycle flux; the multi-source environmental data also includes digital elevation model data, used to provide a fine three-dimensional topographic spatial benchmark for the target area.

[0018] As a preferred embodiment of the regional groundwater storage change dynamic monitoring method based on gravity satellite time-series data analysis according to the present invention, the step S2 of resampling the multi-source environmental data to the same spatial resolution as the gravity satellite data specifically includes: extracting the spatial grid parameters of the land water storage data as the target projection reference; and uniformly projecting the normalized vegetation index data, surface temperature data, rainfall data, and digital elevation model data in the multi-source environmental data to the coordinate system where the target projection reference is located.

[0019] For normalized vegetation index, surface temperature and rainfall data with continuous values, the pixel averaging method is used to calculate the arithmetic mean falling within each target grid; for the digital elevation model data, the bilinear interpolation method is used to extract the grid center elevation and simultaneously calculate the standard deviation of topographic relief within the grid.

[0020] The processed data is uniformly mapped to the grid cells of the target projection reference to generate an environmental dataset that is spatially aligned with the land water storage data.

[0021] As a preferred embodiment of the regional groundwater storage change dynamic monitoring method based on gravity satellite time-series data analysis described in this invention, the process of obtaining the non-groundwater storage in the target area in step S3 specifically includes: calculating the regional cumulative precipitation recharge using the rainfall data, and estimating the surface evapotranspiration loss by combining the surface temperature data and the normalized vegetation index data; performing preliminary correction on the basic soil moisture data in the area based on the difference between the regional cumulative precipitation recharge and the surface evapotranspiration loss; using the normalized vegetation index data as a biomass proxy index, calculating the vegetation biological water content per unit area through nonlinear mapping; and obtaining the vegetation biological water content record of the previous monitoring period, calculating the difference in biological water content between the current period and the previous period, and defining it as the vegetation canopy water storage change value.

[0022] Then, lateral runoff correction based on terrain constraints is performed. The digital elevation model is used to set a lateral flow trigger slope threshold, and lateral outflow grids and their downstream receiving grids with slopes exceeding the slope threshold are identified. The soil moisture content of the lateral outflow grids exceeds the field capacity. The lateral water migration is calculated based on the terrain gradient and soil hydraulic conductivity, and the basic soil moisture data and snow water equivalent data of the corresponding grids and their downstream receiving grids are redistributed and updated using the lateral migration. Excess water that exceeds the soil water-holding capacity and does not form lateral outflow is marked as a seepage component and removed from the non-groundwater storage.

[0023] The corrected soil moisture data, snow water equivalent data, and vegetation canopy water storage change values ​​are spatially superimposed to generate the non-groundwater storage that reflects the distribution of the comprehensive non-groundwater body in the target area.

[0024] As a preferred embodiment of the regional groundwater storage change dynamic monitoring method based on gravity satellite time series data analysis according to the present invention, the process of calculating the groundwater storage distribution data of the target area in step S3 specifically includes: performing spatial grid operation based on the water balance principle in the spatial coordinate system defined by the target projection reference; deducting the non-groundwater storage from the terrestrial water storage data in the spatial grid, obtaining the water surplus in each grid cell as the state parameter of the underground aquifer in the spatial grid, and forming the groundwater storage distribution data of the target area.

[0025] As a preferred embodiment of the regional groundwater storage change dynamic monitoring method based on gravity satellite time-series data analysis described in this invention, the construction and updating process of the groundwater dynamic evolution database specifically includes: collecting groundwater storage change records of the target area based on historical gravity satellite observation data inversion and archiving, and establishing an initial groundwater dynamic evolution database after time-dimensional standardization processing; forming time-series increments by standardizing the regional groundwater time-series distribution data according to the time dimension, and adding them to the groundwater dynamic evolution database; based on the updated groundwater dynamic evolution database, traversing the full-time series records of each grid unit in the groundwater dynamic evolution database, recalculating the full-time water level mean and standard deviation covering historical and current data, and constructing a dynamic comparison benchmark for identifying abnormal losses.

[0026] As a preferred embodiment of the regional groundwater storage dynamic monitoring method based on gravity satellite time series data analysis of the present invention, the comparative analysis process in step S4 specifically includes: extracting the long-term trend component from the regional groundwater time series distribution data, calculating the first derivative of the long-term trend component with time, and obtaining the groundwater storage change rate.

[0027] The rate of change of groundwater storage is compared with the historical average rate of change in the groundwater dynamic evolution database. If the current rate of change is negative and the absolute value exceeds a preset multiple of the historical average rate of change, it is determined that there is an accelerated depletion trend of groundwater in the area.

[0028] The statistical benchmark is updated using all data in the groundwater dynamic evolution database. The deviation of the current groundwater storage data from the updated statistical benchmark is calculated. If the deviation exceeds a preset multiple of the historical standard deviation, the area is determined to be in an abnormal water level state.

[0029] As a preferred embodiment of the regional groundwater storage change dynamic monitoring method based on gravity satellite time-series data analysis described in this invention, the dynamic monitoring report specifically includes: a spatial distribution map of abnormal areas composed of geographical grid areas identified as having an accelerated groundwater depletion trend and abnormal water level status, in the form of a visual layer; a table listing the groundwater storage change rate values ​​of grid points within the abnormal areas and their deviation from historical benchmarks as a depletion rate quantification table; and assigning graded warning color labels to different abnormal areas based on the degree to which the change rate exceeds a preset multiple and the standard deviation level of the deviation, along with targeted water resource management recommendations.

[0030] The beneficial effects of this invention are as follows: By constructing a full-chain monitoring system from satellite observation to physical inversion and dynamic early warning, this invention effectively solves the contradiction between spatial coverage and the separation of multi-source water components in traditional groundwater monitoring methods. By introducing satellite laser ranging coefficients and combining them with combined filtering and signal recovery techniques, the invention eliminates striping errors and high-frequency noise while restoring the true signal to the greatest extent, ensuring high signal-to-noise ratio and reliability of large-scale inversion benchmark data. It accurately reconstructs the distribution of non-groundwater using multi-source environmental data and deducts it from the spatially aligned total storage benchmark based on the water balance principle, thereby avoiding scale mismatch problems while achieving the extraction and physical separation of groundwater signals. By constructing a dynamic database containing long-term historical records as a statistical benchmark, and utilizing dual-dimensional diagnosis of rate and state, the invention keenly identifies abnormal loss trends masked by natural fluctuations and generates dynamic reports integrating spatial perception, quantitative statistics, and decision-making suggestions, providing strong decision support for the precise management of regional water resources. 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 A flowchart of a method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis; Figure 2 Here is a flowchart of the preprocessing and inversion process; Figure 3 Here is a flowchart of the component separation and calculation process; Figure 4 This is a flowchart for monitoring and early warning. Detailed Implementation

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

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

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

[0036] Example 1

[0037] Reference Figures 1-4 This is one embodiment of the present invention, which provides a method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis, including the following steps:

[0038] S1. Obtain time-varying gravity field spherical harmonic data of the target area through gravity satellite, and invert the land water storage data of the target area from the time-varying gravity field spherical harmonic data;

[0039] Before inverting the time-varying gravity field spherical harmonic data, a preprocessing step is also included to obtain the C20 coefficients in the satellite laser ranging data, replace the corresponding terms in the time-varying gravity field spherical harmonic data with the C20 coefficients, and correct the reference frame origin drift by introducing a geocentric correction term.

[0040] Then, polynomial fitting is performed to remove the correlation between odd and even order coefficients to eliminate strip error, and Gaussian smoothing operator is applied to suppress noise in higher order terms.

[0041] The signal leakage scale factor is determined by the signal strength attenuation ratio before and after combined filtering, and the amplitude of the filtered data is corrected using the scale factor.

[0042] The process of inverting to obtain terrestrial water storage data for the target area involves introducing the Earth Load Love number to weight and correct the preprocessed spherical harmonic coefficients, eliminating the interference of Earth's elastic deformation on the gravity signal. The corrected spherical harmonic coefficients are then substituted into a physical conversion formula that includes the Earth's average radius, average density, and water density, transforming the dimensionless spherical harmonic coefficients into mass coefficients with equivalent physical meaning of water height. Using spherical harmonic synthesis, the mass coefficients are mapped to the spatial domain, and the values ​​at each latitude and longitude grid point within the target area are calculated point by point, thus generating terrestrial water storage data reflecting the target area.

[0043] The process begins by accessing Level-2 monthly spherical harmonic coefficient data of the Earth's gravity field, released by the Gravity Recovery and Climate Experiment or its subsequent missions, via a dedicated data interface. This data represents the long-wavelength components of the Earth's gravity field in the form of spherical harmonic coefficients truncated to a specific order. Before data processing, a crucial coefficient replacement operation is performed. Higher-precision C20 coefficients from the International Laser Ranging Service are obtained and directly replace the corresponding terms in the original data. Geocentric correction data for the first-order coefficients are also introduced to align the reference frame origin from the Earth's center of mass to the Earth's geometric center, thus eliminating seasonal drift errors caused by geocentric motion.

[0044] To address the north-south high-frequency stripe error caused by satellite orbital resonance in the original data, a polynomial fitting decorrelation algorithm was executed. Within the coefficient sequences of odd and even categories, the background trend of coefficient changes with order was fitted and subtracted, thus removing the generation mechanism of the stripe error at the spectral level. Subsequently, to suppress residual isotropic high-frequency noise in higher-order terms, a Gaussian smoothing operator was applied. A preset smoothing radius was used to assign exponentially decaying weights to the higher-order coefficients as the order increases, thereby smoothing out shortwave noise. Since the above filtering operation attenuates the true signal amplitude while removing noise, an amplitude correction step is then required. By calculating the ratio of signal intensity before and after the same filtering process for historical average water storage simulation data, a spatially distributed signal leakage scale factor is quantified. This factor is then used to perform point-by-point multiplicative correction on the filtered actual observation data to restore the attenuated true signal amplitude.

[0045] After data preprocessing, the Love number of Earth's load is introduced to weight and correct each spherical harmonic coefficient, eliminating the influence of elastic deformation of the solid Earth under surface load on the gravitational potential coefficient, and purifying the surface mass change signal. Then, the corrected dimensionless spherical harmonic coefficients are substituted into a physical conversion equation that includes the Earth's average radius, Earth's average density, and liquid water density, transforming them into mass coefficients with physical meaning of equivalent water content values.

[0046]

[0047] in, Represents the calculated grid points The equivalent water content value at the location; The average radius of the Earth; This represents the Earth's average density. The average density of liquid water; These represent the order and degree of the spherical harmonic coefficients, respectively. The maximum truncation order; It is a fully normalized Legendre polynomial; Indicates the order is The Earth's load Love number at that time; These represent the spherical harmonic coefficients after noise reduction.

[0048] Finally, by using the spherical harmonic synthesis algorithm and weighted summation of Legendre polynomials, the quality coefficients in the frequency domain are mapped to the spatial domain, and the values ​​at each latitude and longitude grid point in the target area are calculated point by point, generating terrestrial water storage data that reflects the water storage surplus and deficit status of the region at a macro scale.

[0049] S2. Acquire multi-source environmental data of the target area and resample the multi-source environmental data to the same spatial resolution as the gravity satellite data;

[0050] Multi-source environmental data includes: normalized vegetation index data reflecting the growth status of surface vegetation, surface temperature data reflecting the thermal state of the surface, and rainfall data reflecting the regional water cycle flux; multi-source environmental data also includes digital elevation model data, which is used to provide a fine three-dimensional topographic spatial benchmark for the target area.

[0051] Multi-source environmental data is resampled to the same spatial resolution as gravity satellite data, and spatial grid parameters of terrestrial water storage data are extracted as the target projection reference. Normalized vegetation index data, surface temperature data, rainfall data and digital elevation model data from multi-source environmental data are uniformly projected to the coordinate system where the target projection reference is located.

[0052] For normalized vegetation index, surface temperature and rainfall data with continuous values, the pixel averaging method is used to calculate the arithmetic mean falling within each target grid; for digital elevation model data, the bilinear interpolation method is used to extract the grid center elevation and simultaneously calculate the standard deviation of topographic relief within the grid.

[0053] The processed data is uniformly mapped to the grid cells of the target projection reference to generate an environmental dataset that is spatially aligned with the terrestrial water storage data.

[0054] This involves parallel acquisition of various environmental element data for the region from distributed environmental monitoring networks and remote sensing databases or GLDAS. Specifically, the collected data includes: Normalized Difference Vegetation Index (NDI) data obtained through inversion from optical remote sensing satellites, which quantifies surface vegetation cover and growth vitality through the spectral response characteristics of vegetation chlorophyll, indirectly indicating the intensity of root water absorption and transpiration; surface temperature data reflecting the thermodynamic state and potential evaporation capacity of the surface, acquired using thermal infrared sensors; and rainfall data directly characterizing the water replenishment flux in the regional water cycle system, extracted from meteorological monitoring stations or precipitation radar grids. Furthermore, high-precision digital elevation model (DEM) data is required to provide a detailed three-dimensional topographic spatial benchmark for subsequent hydrological analysis, determining the convergence direction of surface runoff and the potential flow direction of groundwater.

[0055] After aggregating the heterogeneous data, a rigorous spatial alignment process is required to eliminate differences in spatial resolution and coordinate system definition between different data sources. First, the terrestrial water storage data is analyzed to extract parameters such as spatial resolution (e.g., 1 degree or 0.5 degrees), latitude and longitude range, and projection method, constructing a standardized target projection reference grid. Then, a geospatial processing engine is invoked to uniformly reproject and transform the normalized vegetation index data, surface temperature data, rainfall data, and digital elevation model data to the coordinate system of this target projection reference, eliminating geometric distortion.

[0056] To address the scale mismatch issue where environmental data generally has high spatial resolution (e.g., 30 meters to 1 kilometer) while terrestrial water storage data has lower resolution (e.g., 300 kilometers), spatial aggregation is required. For continuously varying physical quantities (e.g., temperature, rainfall, vegetation index), a regional averaging algorithm is used to calculate the arithmetic mean of all high-resolution pixel values ​​falling within each target grid cell. For topographic data, the average elevation and elevation standard deviation (characterizing topographic relief) within the grid are calculated simultaneously. Through spatial resampling and numerical aggregation, surface environmental data are uniformly mapped spatially to the grid system of gravity inversion data, generating an environmental dataset that spatially corresponds one-to-one with terrestrial water storage data and retains the environmental statistical characteristics within the grid. This lays a unified data foundation for establishing a quantitative relationship between environmental factors and water storage.

[0057] S3. Obtain the non-groundwater reserves of the target area by analyzing multi-source environmental data, and then calculate the groundwater reserve distribution data of the target area;

[0058] The process of obtaining non-groundwater reserves in the target area involves calculating the cumulative precipitation recharge in the area using rainfall data, and estimating the surface evapotranspiration loss by combining surface temperature data and normalized vegetation index data. Based on the difference between the cumulative precipitation recharge and the surface evapotranspiration loss in the area, the basic soil moisture data in the area is preliminarily corrected.

[0059] Using the normalized vegetation index data as a biomass proxy, the vegetation biological water content per unit area is calculated through nonlinear mapping; and the vegetation biological water content record of the previous monitoring period is obtained, the difference in biological water content between the current period and the previous period is calculated, and it is defined as the change value of vegetation canopy water storage.

[0060] Then, terrain-constrained lateral runoff correction is performed. A slope threshold for lateral flow triggering is set using a digital elevation model. Lateral outflow grids and their downstream receiving grids with slopes exceeding the slope threshold are identified. The soil moisture content of the lateral outflow grids exceeds the field capacity. The lateral migration of water is calculated based on the terrain gradient and soil hydraulic conductivity. The lateral migration is then used to redistribute and update the basic soil moisture data and snow water equivalent data of the corresponding grids and their downstream receiving grids. Excess water that exceeds the soil water-holding capacity but does not form lateral outflow is marked as a seepage component and removed from the non-groundwater storage.

[0061] The corrected soil moisture data, snow water equivalent data and vegetation canopy water storage change values ​​are spatially superimposed to generate non-groundwater storage that reflects the distribution of non-groundwater comprehensive water bodies in the target area.

[0062] The process of calculating the groundwater storage distribution data of the target area involves performing a spatial grid operation based on the water balance principle in the spatial coordinate system defined by the target projection reference; subtracting non-groundwater storage from the terrestrial water storage data within the spatial grid, and obtaining the water surplus in each grid cell as the state parameter of the groundwater aquifer in the spatial grid to form the groundwater storage distribution data of the target area.

[0063] The process of constructing and updating the groundwater dynamic evolution database involves collecting groundwater storage change records of the target area based on historical gravity satellite observation data. After time-dimensional standardization, an initial groundwater dynamic evolution database is established. The time-series distribution data of the regional groundwater is then standardized according to the time dimension to form time-series increments, which are added to the groundwater dynamic evolution database. Based on the updated groundwater dynamic evolution database, the full time-series records of each grid unit in the database are traversed, and the mean and standard deviation of the groundwater level covering both historical and current data are recalculated to construct a dynamic comparison benchmark for identifying abnormal deficits.

[0064] The process begins with in-depth analysis of non-surface features. For the input rainfall data, it's not simply a matter of accumulation; instead, an effective rainfall calculation based on a time-sliding window is performed. By analyzing rainfall intensity and duration, the proportion of precipitation that can infiltrate the soil and provide effective replenishment is estimated, thus obtaining a more accurate cumulative rainfall replenishment figure that aligns with hydrophysical processes. In estimating evapotranspiration loss, surface temperature data is used as a thermal driver, and the normalized difference in vegetation index (NDVI) is used as a biophysical constraint. Physical evaporation from the soil surface and physiological transpiration from plant leaves are calculated separately, and these two calculations are combined to obtain the total surface evapotranspiration loss.

[0065] The Normalized Difference Vegetation Index (NDVI) is used as a dynamic proxy indicator for vegetation biomass. By establishing a nonlinear mapping relationship between the NDVI and the water content per unit area of ​​vegetation, the amount of biological water absorbed and stored within the tissues of plants (including leaves and stems) due to growth activities is quantified. The system calculates the difference in vegetation water content between the current period and the previous monitoring period, defining it as the change in vegetation canopy water storage, thereby accurately capturing water storage fluctuation signals caused by changes in vegetation phenology.

[0066] To address the accuracy issues of traditional assimilation data in complex terrain, the system introduces digital elevation model (DEM) data as a spatial regulator. The system first calculates the slope, aspect, and flow direction matrices of the DEM data to identify catchment areas and watersheds. Based on the physical law of water flow migrating from high to low elevations, and combined with the previously calculated difference between precipitation and evapotranspiration, a lateral runoff redistribution mechanism is established. This mechanism simulates the lateral migration of excess water from high-slope grids to low-lying grids during periods of abundant rainfall and low evaporation.

[0067] The lateral runoff redistribution mechanism utilizes a digital elevation model and employs the D8 unidirectional flow algorithm to traverse the entire region's grid. A lateral flow trigger slope threshold is set, and grids with slopes exceeding this threshold are marked as potential lateral outflow sources. The adjacent grid with the largest elevation drop within its 8-neighborhood is then identified as the downstream receiving grid. For grids marked as outflow sources, the field capacity of the soil type is set as the water content trigger threshold. The current soil volumetric water content of the grid is monitored in real time. If it is lower than or equal to the field capacity, the water is considered to be bound by capillary forces and not undergoing lateral migration; if it is higher than the field capacity, movable gravity water is considered to be present.

[0068]

[0069] in, This represents the lateral water flux flowing out of the target grid; It is the anisotropic attenuation coefficient, used to correct for grid scale effects; The saturated hydraulic conductivity of the soil in this grid is given. It is the terrain slope angle of the target grid; This represents the soil volumetric water content of the grid in the current time period; This refers to the field water holding capacity of the soil type in this grid. It is the effective soil layer thickness that participates in lateral flow; Width of the grid cell

[0070] Lateral outflow flux was calculated based on a modified Darcy's law. This calculation process comprehensively considered soil saturated hydraulic conductivity, the tangent of topographic gradient, effective gravitational water content exceeding field capacity, effective soil layer thickness, and grid width. An anisotropic attenuation coefficient negatively correlated with vegetation root density was introduced to correct for velocity overestimation. Each grid was synchronously updated based on the calculated lateral outflow flux for the entire region. The corrected final soil moisture content was calculated by subtracting the outflow to adjacent lower areas from the original moisture content and adding the sum of inflow from all upstream adjacent higher areas. This physically simulated the process of water converging from ridges to valleys, achieving a refined reconstruction of non-groundwater storage.

[0071]

[0072] in, This is the final soil volumetric water content after the grid update; This represents the soil volumetric water content of the grid in the current time period; This represents the lateral water flux flowing out of the target grid; This represents the effective soil volume within the grid cell; For all upstream neighbor grids pointing to the current grid The sum of the lateral outflow fluxes.

[0073] This mechanism is used to perform spatial fine-tuning of the basic soil moisture data (typically including four soil layers: 0-10cm, 10-40cm, 40-100cm, and 100-200cm) and snow water equivalent data provided by the Global Land Surface Assimilation System in this region. Specifically, for low-lying grids identified as catchment areas, their soil moisture content readings are appropriately increased based on the upstream runoff; for sunny grids with steep slopes, their snow water equivalent readings are appropriately decreased based on the accelerated snowmelt effect caused by high radiation.

[0074] After completing the lateral runoff redistribution calculation, for grid cells where the soil moisture content still exceeds the field capacity after lateral recharge, this excess gravity water is determined to no longer remain in the shallow soil but has undergone downward vertical infiltration. This portion of water is marked as the infiltration component and removed from the current non-groundwater storage statistics. This ensures that the shallow water data only includes water bodies truly present in the surface and rhizosphere soils, preventing the erroneous deduction of infiltration water already recharged to the groundwater aquifer in subsequent subtraction calculations.

[0075] After completing these physical mechanism-based corrections, the corrected soil moisture data at all depths, the corrected snow water equivalent data, and the canopy interception water volume inferred from the vegetation leaf area index are numerically superimposed grid by grid to generate a non-groundwater storage map that can truly reflect the distribution of non-groundwater comprehensive water bodies in the region.

[0076] After acquiring non-groundwater reserves, groundwater signal separation is performed in a unified spatial coordinate system. Using terrestrial water reserve data as the benchmark for total water reserves, a grid-by-grid interpolation calculation is performed based on the water balance principle. Specifically, the system traverses each grid cell and directly subtracts the corresponding non-groundwater reserves from the macroscopic total water reserves value of that cell.

[0077]

[0078] in, This represents the groundwater storage value within grid cell (i,j); The total water storage value of the grid obtained by inversion; This is the corrected shallow soil moisture value; This is the corrected snow water equivalent value; This represents the water storage capacity intercepted by the vegetation canopy.

[0079] Non-groundwater signals such as surface snow cover, canopy intercepted water, and shallow soil water are separated, and the remaining water volume residuals are locked as changes in groundwater aquifers. To ensure the spatial continuity of the data, the calculated residual values ​​of each grid are combined and necessary smoothing is performed to finally generate groundwater storage distribution data for the target area.

[0080] The construction and updating process of the groundwater dynamic evolution database is executed synchronously. During the initialization phase, the system automatically connects to the historical data archive center to trace back and aggregate all historical groundwater storage records for the target area over the past few decades (e.g., 20 years of data since the launch of the GRACE satellite). Since historical data may come from different processing versions or time sampling intervals, strict time dimension standardization processing is also required. All non-standard timestamp data is uniformly resampled or interpolated to standard monthly time nodes, and obvious outliers are removed to establish an initial groundwater dynamic evolution database with controllable quality and continuous time series.

[0081] After calculating the latest regional groundwater time-series distribution data, it is used as a new time-series increment. First, the completeness and physical validity of this increment data are checked. Once confirmed to be correct, it is appended to the end of the database according to a standard time format, automatically extending the database's time span to the current time. Whenever new data is added, the background analysis engine immediately triggers a full database scan. The engine extracts long-term time-series sample sets for each spatial grid by month. Based on this sample set, the historical mean water level (representing the normal water level) and standard deviation (representing the natural fluctuation range) for that grid in the current season are recalculated. These two sets of statistical parameters constitute the latest dynamic comparison benchmark, characterizing the seasonal fluctuations of groundwater in the region with seasonal and climatic changes, providing a solid statistical basis for accurately identifying abnormal deficits exceeding the natural fluctuation range.

[0082] S4. Continuously monitor the target area, construct regional groundwater time series distribution data, compare and analyze the regional groundwater time series distribution data with the groundwater dynamic evolution database of the target area, and generate a dynamic monitoring report after identifying abnormal groundwater depletion.

[0083] The comparative analysis process involves extracting the long-term trend component from the time-series distribution data of regional groundwater, calculating the first derivative of the long-term trend component over time, and obtaining the rate of change of groundwater storage.

[0084] The rate of change of groundwater storage is compared with the historical average rate of change in the groundwater dynamic evolution database. If the current rate of change is negative and the absolute value exceeds the preset multiple of the historical average rate of change, it is determined that there is an accelerated depletion trend of groundwater in the area.

[0085] The statistical benchmark is updated using the full data in the groundwater dynamic evolution database. The deviation of the current groundwater storage data from the updated statistical benchmark is calculated. If the deviation exceeds a preset multiple of the historical standard deviation, the area is determined to be in an abnormal water level state.

[0086] The dynamic monitoring report uses a visual layer to identify the spatial distribution of anomalous areas, which are composed of geographic grid areas that are determined to have an accelerated trend of groundwater depletion and abnormal water level status. It lists the groundwater storage change rate values ​​of grid points in the anomalous area and the degree of deviation from the historical benchmark as a quantitative table of loss rate. Based on the degree to which the change rate exceeds a preset multiple and the standard deviation level of the deviation, different anomalous areas are assigned graded warning color labels, and targeted water resource management recommendations are attached.

[0087] Once in continuous monitoring mode, groundwater storage data for discrete time periods are serialized and stored according to timestamps to construct a spatiotemporal cube covering the target area, i.e., regional groundwater time-series distribution data. To understand the hydrological driving mechanisms behind the data, the system introduces a time-series decomposition algorithm to decouple the original observation signal into three orthogonal components: a seasonal term representing the natural hydrological cycle, a residual term representing random climate disturbances, and a trend term representing the long-term surplus and deficit state of the groundwater aquifer. Based on this, focusing on the trend term, a time-dimensional differential operation is performed on it to calculate the first derivative of this component with time, thereby obtaining the groundwater storage change rate, which has a clear physical orientation and quantifies the rate of consumption or recovery of groundwater resources.

[0088]

[0089]

[0090] in, express Raw time-series data of groundwater storage observed during the period; The separated long-term trend component; The separated seasonal periodic components; These are the separated random residual components;

[0091] express The rate of change of groundwater storage over a period of time; The time sampling interval is typically 1 month.

[0092] After obtaining the current state parameters, the groundwater dynamic evolution database is used as a reference to perform anomaly diagnosis in two dimensions: rate and state. The historical average rate of change stored in the database is extracted as a baseline, and the currently calculated rate of change in groundwater reserves is compared with this baseline. If the current rate is negative and its absolute value significantly exceeds a preset multiple of the historical baseline, it indicates that groundwater consumption in the area is deviating from historical norms, thus indicating an accelerating trend of groundwater depletion. The deviation of the current groundwater reserve value from the historical average is calculated, and the historical standard deviation is introduced as a unit of measurement for volatility. If the deviation exceeds a preset standard deviation threshold in the negative direction, it means that the current water level has fallen below the safe lower limit allowed by natural climate fluctuations in the area, thus indicating that the area is in an abnormal water level state.

[0093]

[0094] in, express Abnormal water level index for a given period; This represents the measured groundwater storage value during that period; This represents the average groundwater storage for the same historical period (same month) stored in the database. It is the standard deviation of historical groundwater reserves for the same period stored in the database.

[0095] Based on the above judgment results, a dynamic monitoring report integrating spatial perception, quantitative statistics, and decision-making recommendations is synthesized. The report first generates a spatial distribution map of the abnormal areas, mapping the identified abnormal grids onto a geographic information base map. Using a spatial clustering algorithm, scattered abnormal grids are aggregated into continuous loss areas, visually presenting the spatial spread of the groundwater crisis. The report also outputs a loss rate quantification table, listing the measured change rate of each key node within the abnormal area, its historical average, and the percentage deviation between the two, providing water resource management departments with grid-level quantitative evidence. The report also integrates trend warning level indicators, establishing a color-coding system based on severity. According to the rate exceeding the standard multiple and the standard deviation of water level deviation, abnormal areas are assigned graded warning colors, and targeted water resource management recommendations are automatically matched and attached based on the warning level, thus achieving closed-loop support from data monitoring to governance decision-making.

[0096] In this invention, gravity satellite data is accessed through publicly available platforms, obtained via data service ports of the National Earth System Science Data Center, the Space Research Center, the Natural Resources Airborne Geophysical and Remote Sensing Center of the Geological Survey, or the Geological Science Research Center.

[0097] Preferably, the gravity field model used in the conception and verification of this invention is the EGM2008 or EIGEN-6C4 Earth gravity field model.

[0098] In summary, this invention effectively resolves the contradiction between spatial coverage and the separation of multiple water source components in traditional groundwater monitoring methods by constructing a full-chain monitoring system from satellite observation to physical inversion and dynamic early warning. By introducing satellite laser ranging coefficients and combining them with combined filtering and signal recovery techniques, the true signal is restored to the greatest extent while eliminating striping errors and high-frequency noise, ensuring high signal-to-noise ratio and reliability of large-scale inversion benchmark data. Non-groundwater distribution is accurately reconstructed using multi-source environmental data and subtracted from the spatially aligned total storage benchmark based on the water balance principle, thus avoiding scale mismatch problems while achieving the extraction and physical separation of groundwater signals. By constructing a dynamic database containing long-term historical records as a statistical benchmark, and utilizing dual-dimensional diagnostics of rate and state, abnormal deficit trends masked by natural fluctuations are keenly identified, and a dynamic report integrating spatial perception, quantitative statistics, and decision-making recommendations is generated, providing strong decision support for the precise management of regional water resources.

[0099] 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 dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis, characterized in that, Performed by a computer device, including the following steps: S1. Obtain time-varying gravity field spherical harmonic data of the target area through gravity satellite, and invert the time-varying gravity field spherical harmonic data to obtain the land water storage data of the target area; S2. Acquire multi-source environmental data of the target area, and resample the multi-source environmental data to the same spatial resolution as the gravity satellite data; S3. Obtain the non-groundwater reserves of the target area by parsing the multi-source environmental data, and then calculate the groundwater reserve distribution data of the target area; S4. Continuously observe the target area, construct regional groundwater time series distribution data, compare and analyze the regional groundwater time series distribution data with the groundwater dynamic evolution database of the target area, and generate a dynamic monitoring report after identifying abnormal groundwater depletion.

2. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 1, characterized in that, Before inverting the time-varying gravity field spherical harmonic data in step S1, a preprocessing step for the time-varying gravity field spherical harmonic data is also included. The preprocessing step specifically includes: Obtain the C20 coefficients from the satellite laser ranging data, replace the corresponding terms in the time-varying gravity field spherical harmonic data with the C20 coefficients, and correct the reference frame origin drift by introducing a geocentric correction term. Then, polynomial fitting is performed to remove the correlation between odd and even order coefficients to eliminate strip error, and Gaussian smoothing operator is applied to suppress noise in higher order terms. The signal leakage scale factor is determined by the signal strength attenuation ratio before and after combined filtering, and the amplitude of the filtered data is corrected using the scale factor.

3. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 2, characterized in that, The process of inverting the terrestrial water storage data of the target area in step S1 specifically includes: introducing the Earth Load Love number to weight and correct the preprocessed spherical harmonic coefficients to eliminate the interference of Earth's elastic deformation on the gravity signal; substituting the corrected spherical harmonic coefficients into a physical conversion formula that includes the Earth's average radius, Earth's average density, and water density to convert the dimensionless spherical harmonic coefficients into mass coefficients with equivalent physical meaning of water height; using spherical harmonic synthesis processing to map the mass coefficients to the spatial domain, calculating the values ​​of each latitude and longitude grid point in the target area point by point, and then generating terrestrial water storage data reflecting the target area.

4. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 3, characterized in that, The multi-source environmental data in step S2 includes: normalized vegetation index data reflecting the growth status of surface vegetation, surface temperature data reflecting the thermal status of the surface, and rainfall data reflecting the regional water cycle flux; the multi-source environmental data also includes digital elevation model data, which is used to provide a fine three-dimensional topographic spatial benchmark for the target area.

5. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 4, characterized in that, The step S2, which involves resampling the multi-source environmental data to the same spatial resolution as the gravity satellite data, specifically includes: extracting the spatial grid parameters of the land water storage data as the target projection reference; and uniformly projecting the normalized vegetation index data, surface temperature data, rainfall data, and digital elevation model data from the multi-source environmental data onto the coordinate system where the target projection reference is located. For normalized vegetation index, surface temperature and rainfall data with continuous values, the pixel averaging method is used to calculate the arithmetic mean falling within each target grid; for the digital elevation model data, the bilinear interpolation method is used to extract the grid center elevation and simultaneously calculate the standard deviation of topographic relief within the grid. The processed data is uniformly mapped to the grid cells of the target projection reference to generate an environmental dataset that is spatially aligned with the land water storage data.

6. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 5, characterized in that, The process of obtaining the non-groundwater reserves in the target area in step S3 specifically includes: calculating the cumulative precipitation replenishment of the area using the rainfall data, and estimating the surface evapotranspiration loss by combining the surface temperature data and the normalized vegetation index data, and performing preliminary correction on the basic soil moisture data in the area based on the difference between the cumulative precipitation replenishment of the area and the surface evapotranspiration loss. Using the normalized vegetation index data as a biomass proxy, the vegetation biological water content per unit area is calculated through nonlinear mapping; and the vegetation biological water content record of the previous monitoring period is obtained, the difference in biological water content between the current period and the previous period is calculated, and it is defined as the change value of vegetation canopy water storage. Then, lateral runoff correction based on terrain constraints is performed. The digital elevation model is used to set a lateral flow trigger slope threshold, and lateral outflow grids and their downstream receiving grids with slopes exceeding the slope threshold are identified. The soil moisture content of the lateral outflow grids exceeds the field capacity. The lateral water migration is calculated based on the terrain gradient and soil hydraulic conductivity, and the basic soil moisture data and snow water equivalent data of the corresponding grids and their downstream receiving grids are redistributed and updated using the lateral migration. Excess water that exceeds the soil water-holding capacity and does not form lateral outflow is marked as a seepage component and removed from the non-groundwater storage. The corrected soil moisture data, snow water equivalent data, and vegetation canopy water storage change values ​​are spatially superimposed to generate the non-groundwater storage that reflects the distribution of the comprehensive non-groundwater body in the target area.

7. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 6, characterized in that, The process of calculating the groundwater storage distribution data of the target area in step S3 specifically includes: performing a spatial grid operation based on the water balance principle in the spatial coordinate system defined by the target projection reference; subtracting the non-groundwater storage from the terrestrial water storage data in the spatial grid, obtaining the water surplus in each grid cell as the state parameter of the underground aquifer in the spatial grid, and forming the groundwater storage distribution data of the target area.

8. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 1, characterized in that, The construction and updating process of the groundwater dynamic evolution database specifically includes: collecting groundwater storage change records of the target area based on historical gravity satellite observation data, and establishing an initial groundwater dynamic evolution database after time-dimensional standardization; adding time-series increments to the groundwater dynamic evolution database after standardizing the time-series distribution data of the area according to the time dimension; and based on the updated groundwater dynamic evolution database, traversing the full-time series records of each grid unit in the groundwater dynamic evolution database, recalculating the full-time water level mean and standard deviation covering historical and current data, and constructing a dynamic comparison benchmark for identifying abnormal deficits.

9. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 8, characterized in that, The comparative analysis process in step S4 specifically includes: extracting the long-term trend component from the time-series distribution data of groundwater in the region, calculating the first derivative of the long-term trend component with time, and obtaining the rate of change of groundwater storage. The rate of change of groundwater storage is compared with the historical average rate of change in the groundwater dynamic evolution database. If the current rate of change is negative and the absolute value exceeds a preset multiple of the historical average rate of change, it is determined that there is an accelerated depletion trend of groundwater in the area. The statistical benchmark is updated using all data in the groundwater dynamic evolution database. The deviation of the current groundwater storage data from the updated statistical benchmark is calculated. If the deviation exceeds a preset multiple of the historical standard deviation, the area is determined to be in an abnormal water level state.

10. The method for dynamic monitoring of regional groundwater storage changes based on gravity satellite time-series data analysis according to claim 9, characterized in that, The dynamic monitoring report specifically includes: a spatial distribution map of abnormal areas, which is identified in the form of a visual layer, consisting of geographic grid areas that are determined to have an accelerated trend of groundwater depletion and abnormal water level status; a table listing the groundwater storage change rate values ​​of the grid points in the abnormal areas and the degree of deviation from historical benchmarks as a depletion rate quantification table; and assigning graded warning color labels to different abnormal areas according to the degree to which the change rate exceeds a preset multiple and the standard deviation level of the deviation, along with targeted water resource management recommendations.

Images