A method and system for synergistic inversion of groundwater storage integrating SWOT and GRACE

Abstract

This invention provides a method and system for collaborative inversion of groundwater storage that integrates SWOT and GRACE methods, belonging to the field of hydrological data processing technology. The method includes: classifying surface water bodies within a target area based on multi-source data; converting water level data into surface water storage observation sequences using water level-storage conversion models corresponding to each water body type; performing temporal reconstruction and spatial aggregation on the surface water storage observation sequences to obtain a gridded surface water storage sequence that matches the total water storage data in terms of spatiotemporal scale; constructing a groundwater storage inversion model, solving the groundwater storage inversion model, and obtaining the groundwater storage sequence for the target area. This invention, through classifying surface water bodies and reconstructing the spatiotemporal data, achieves refined calculation of surface water storage and alignment of observation data in the spatiotemporal dimensions, thereby accurately reflecting the actual changes in surface water under complex hydrological environments and realizing high-precision separation and inversion of groundwater storage.

Description

Technical Field

[0001] This invention relates to the field of hydrological data processing technology, and in particular to a method and system for synergistic inversion of groundwater storage by integrating SWOT and GRACE. Background Technology

[0002] Groundwater storage inversion refers to the process of separating and calculating the groundwater components in terrestrial water storage using observational data such as satellite remote sensing, and it is an important part of water resource management.

[0003] Currently, groundwater storage inversion is usually carried out using a simple separation strategy, which involves subtracting surface water storage from the total observed terrestrial water storage item by item to obtain the groundwater storage.

[0004] However, this method is difficult to accurately reflect the actual changes in surface water under complex hydrological environments and lacks a refined calculation strategy for surface water storage, which makes it difficult to achieve high-precision separation and inversion of groundwater storage. Summary of the Invention

[0005] This invention provides a groundwater storage collaborative inversion method and system that integrates SWOT and GRACE, in order to solve the shortcomings of existing technologies that lack a refined calculation strategy for surface water storage, making it difficult to accurately reflect the actual changes in surface water under complex hydrological environments, thus making it difficult to achieve high-precision separation and inversion of groundwater storage.

[0006] This invention provides a method for synergistic inversion of groundwater storage by integrating SWOT and GRACE, comprising the following steps: Acquire water level data, total water storage data, and multi-source data for the target area; the multi-source data includes optical image data and topographic data of the target area. Based on the multi-source data, the surface water bodies in the target area are classified into multiple water body types; using the water level-storage conversion model corresponding to each water body type, the water level data is converted into a surface water storage observation sequence. Temporal reconstruction and spatial aggregation are performed on the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of temporal and spatial scale; Based on the gridded surface water storage sequence and the total water storage data, a groundwater storage inversion model is constructed. The groundwater storage inversion model is solved to obtain the groundwater storage sequence of the target area.

[0007] According to the present invention, a groundwater storage synergistic inversion method integrating SWOT and GRACE is provided, wherein the groundwater storage inversion model is constructed based on the gridded surface water storage sequence and the total water storage data, comprising: Based on the gridded surface water storage sequence, the hysteresis transfer function of surface water and groundwater is constructed; Based on the hysteresis response transfer function and the total water storage data, water balance constraints are determined. Obtain measured data on groundwater level changes, and determine groundwater level observation constraints based on the measured data on groundwater level changes; Using the groundwater storage sequence as the variable to be solved, the groundwater storage inversion model is constructed based on the water balance constraint, the groundwater level observation constraint, the time series smoothing constraint, and the physical process constraint. The temporal smoothing constraint term is used to constrain the fluctuation range of groundwater storage changes at adjacent time points, and the physical process constraint term is used to constrain the increase and decrease of groundwater storage.

[0008] According to the groundwater storage synergistic inversion method integrating SWOT and GRACE provided by the present invention, the water balance constraint term is calculated based on the following mathematical model: ; in, This is a water balance constraint term. t For time variables, The total water storage data is as follows. This is the groundwater storage sequence. For the gridded surface water storage sequence, For changes in soil water storage, For changes in snow water equivalent, This represents the amount of water exchanged between surface water and groundwater.

[0009] According to the groundwater storage synergistic inversion method integrating SWOT and GRACE provided by the present invention, the step of performing temporal reconstruction and spatial aggregation on the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of temporal and spatial scale includes: A water balance differential equation is established based on the precipitation, inflow, outflow and evaporation of the target area; The surface water storage observation sequence is input into the water balance differential equation to obtain the daily-scale surface water storage sequence of each surface water body. The daily-scale surface water storage sequence is then averaged monthly to obtain the monthly averaged surface water storage of the surface water body. The target area is divided into multiple grid cells based on the spatial resolution of the total water storage data; Within each of the grid cells, structural weights are calculated based on the area and hydraulic connectivity of each of the surface water bodies, and response weights are calculated based on the correlation between the changes in the storage of each of the surface water bodies and the overall storage changes of the grid. Based on the structural weights and the response weights, the aggregate weights of each surface water body are determined. The monthly average surface water storage of each surface water body and the aggregate weights are then weighted and fused to obtain the gridded surface water storage sequence.

[0010] According to the groundwater storage synergistic inversion method integrating SWOT and GRACE provided by the present invention, the determination of the aggregate weight of each surface water body is calculated based on the following mathematical model: ; in, For the first i Aggregate weight of each surface water body For the first i Structural weights of individual surface water bodies For the first i Response weights of individual surface water bodies The fusion coefficient is the coefficient of fusion. , is the coefficient of variation of the water area within the grid cell.

[0011] According to the present invention, a groundwater storage synergistic inversion method integrating SWOT and GRACE is provided, wherein the water balance differential equation is established based on the precipitation, inflow, outflow, and evaporation of the target area, including: Climate zones are defined based on the annual precipitation of the target area, identifying arid and humid zones. The evaporation rate is determined using a first evaporation model for the arid region and a second evaporation model for the humid region. If the target area has a snow cover period, a water balance differential equation is established based on the precipitation, inflow, outflow, snowmelt replenishment, and evaporation of the target area. If the target area belongs to the monsoon climate zone, a water balance differential equation is established based on the precipitation, inflow, outflow, seasonal adjustment term, and evaporation of the target area.

[0012] According to the present invention, a groundwater storage synergistic inversion method integrating SWOT and GRACE is provided. The step of solving the groundwater storage inversion model to obtain the groundwater storage sequence of the target area includes: Obtain the water balance calculation results and use them as the current solution vector; The model optimization and update steps are iteratively executed until the relative change of the solution vectors obtained in two adjacent iterations is less than the preset convergence threshold. The final output solution vector is then determined as the groundwater storage sequence of the target area. The model optimization and update steps include: Calculate the gradient of the groundwater storage inversion model at the current solution vector; The inverse approximation of the Hessian matrix is ​​constructed using the gradient difference and displacement vector stored in the historical iterations to determine the search direction; Perform a line search along the search direction to determine the update step size; Based on the current solution vector, the search direction, and the update step size, an updated solution vector is determined, and the updated solution vector is projected into a preset feasible region as the solution vector for the next iteration.

[0013] According to the present invention, a groundwater storage synergistic inversion method integrating SWOT and GRACE is provided, wherein the classification of surface water bodies within the target area based on the multi-source data includes: Based on the optical image data, the boundaries of each surface water body in the target area are extracted, and the area, aspect ratio, shape index and stability coefficient of each surface water body are determined based on the boundaries of each surface water body. When the aspect ratio of the surface water body is greater than a first threshold and the shape index is greater than a second threshold, the surface water body is identified as a river. When the stability coefficient of the surface water body is greater than the third threshold, the surface water body is identified as a seasonal water body; Based on the geological data of the target area, the topographic features of the dam body are obtained. When the elevation change gradient of the topographic features of the dam body is greater than the fourth threshold and the extension length is greater than the fifth threshold, the water body is identified as a reservoir; otherwise, the water body is classified as a lake.

[0014] According to the present invention, a groundwater storage synergistic inversion method integrating SWOT and GRACE is provided. The water level storage conversion models corresponding to each water body type include: lake storage conversion model, river storage conversion model, reservoir storage conversion model, and seasonal water body storage conversion model.

[0015] This invention also provides a groundwater storage collaborative inversion system integrating SWOT and GRACE, comprising the following modules: The data acquisition module is used to acquire water level data, total water storage data, and multi-source data of the target area; the multi-source data includes optical image data and topographic data of the target area. The classification and conversion module is used to classify the surface water bodies in the target area based on the multi-source data to obtain multiple water body types; and to convert the water level data into a surface water storage observation sequence using a water level storage conversion model corresponding to each water body type. The spatiotemporal reconstruction module is used to perform temporal reconstruction and spatial aggregation on the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of spatiotemporal scale. The collaborative inversion module is used to construct a groundwater storage inversion model based on the gridded surface water storage sequence and the total water storage data, solve the groundwater storage inversion model, and obtain the groundwater storage sequence of the target area.

[0016] The groundwater storage collaborative inversion method and system integrating SWOT and GRACE provided by this invention achieves refined calculation of surface water storage and alignment of observation data in the spatiotemporal dimension by classifying surface water bodies and reconstructing observation data. This enables accurate reflection of the actual changes in surface water under complex hydrological environments and achieves high-precision separation and inversion of groundwater storage. Attached Figure Description

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

[0018] Figure 1 This is a flowchart illustrating the groundwater storage synergistic inversion method integrating SWOT and GRACE provided by the present invention.

[0019] Figure 2 This is a schematic diagram of the process for constructing a groundwater storage inversion model provided by the present invention.

[0020] Figure 3 This is a schematic diagram of the process for temporal reconstruction and spatial aggregation of surface water storage observation sequences provided by the present invention.

[0021] Figure 4 This is an architecture diagram of the groundwater storage synergistic inversion system that integrates SWOT and GRACE provided by the present invention. Detailed Implementation

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

[0023] It should be noted that, in the description of this invention, the terms "comprising," "including," or any other variations thereof are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0024] The terms "first," "second," etc., used in this invention are used to distinguish similar objects and not to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention can be implemented in orders other than those illustrated or described herein, and the objects distinguished by "first," "second," etc., are generally of the same class and the number of objects is not limited; for example, a first object can be one or more.

[0025] To facilitate a full understanding of the technical solution of this application, the following content is hereby introduced: With the escalating global water crisis and the increasingly significant impacts of climate change, monitoring changes in groundwater storage has become a core scientific issue in water resource management and ecological protection. The Gravity Recovery and Climate Experiment (GRACE) satellite, by measuring the time-varying characteristics of the Earth's gravity field to invert changes in terrestrial water storage, has been widely applied in global drought monitoring and groundwater over-extraction identification. Particularly in groundwater over-extraction areas such as the North China Plain and the Ganges River Basin in India, decision-making bodies need to integrate GRACE total water storage observations, high-precision surface water level data from Surface Water and Ocean Topography (SWOT), soil moisture estimation from land surface models, and measured water levels from ground monitoring wells. However, due to significant differences in observation scale, temporal resolution, and physical dimensions, these data result in large errors in the separate calculation of surface water and groundwater storage, severely affecting the accuracy of over-extraction area identification and the scientific validity of water resource regulation decisions.

[0026] Traditional GRACE methods for retrieving groundwater reserves rely primarily on a simple subtraction separation strategy, deducting surface water, soil water, and snowmelt components from the total terrestrial water reserves. However, this approach has significant shortcomings: First, surface water reserves estimation is too coarse. Existing methods use the same area-volume empirical formula for all types of water bodies, ignoring the fundamental differences in morphological characteristics and water level-storage relationships among different water body types, resulting in significant variations in accuracy across different surface water body types. Second, the 21-day revisit period of the SWOT satellite and the monthly scale data of GRACE are mismatched in terms of time scale. Existing interpolation methods are based solely on mathematical fitting and do not incorporate the physical mechanisms of hydrological processes, leading to larger estimation errors during extreme weather events. Third, the point-to-line water bodies observed by SWOT differ in spatial scale from the grid scale of GRACE, and simple area-weighted aggregation ignores the hydraulic connectivity between water bodies. Fourth, traditional subtraction assumes that each water component changes independently, neglecting the mutual replenishment and discharge relationships between surface water and groundwater and the hysteresis response mechanism, which cannot meet the accuracy requirements for precise identification of over-extraction areas.

[0027] Therefore, this invention provides a groundwater storage collaborative inversion method and system that integrates SWOT and GRACE, which can accurately reflect the actual changes of surface water under complex hydrological environments and achieve high-precision separation and inversion of groundwater storage.

[0028] The following is combined with Figures 1-4 This invention describes the groundwater storage synergistic inversion method and system that integrates SWOT and GRACE.

[0029] Figure 1 This is a flowchart illustrating the groundwater storage synergistic inversion method integrating SWOT and GRACE provided by the present invention, as shown below. Figure 1 As shown, the execution subject of the groundwater storage synergistic inversion method integrating SWOT and GRACE provided by the present invention can be a server, a cloud computing platform, or a computer capable of executing the method of the present invention, etc. Unless otherwise specified, the following embodiments will be described using a server as an example.

[0030] As an optional embodiment, this groundwater storage synergistic inversion method integrating SWOT and GRACE mainly includes, but is not limited to, the following steps: Step 110: Obtain water level data, total water storage data, and multi-source data for the target area; the multi-source data includes optical image data and topographic data of the target area.

[0031] The target area refers to the geographical region where groundwater storage inversion and monitoring are required. This target area typically involves hydrological processes involving the transformation of surface water and groundwater. For example, the target area can be a region globally with needs for groundwater resource management, drought monitoring, or identification of groundwater over-extraction, especially typical areas with severe groundwater over-extraction such as the North China Plain, or a complex hydrological basin containing various surface water body types such as lakes, rivers, and reservoirs.

[0032] Water level data refers to observational data that reflects changes in the surface elevation of water bodies within a target area. For example, water level data can be high-precision satellite altimetry data with a certain temporal resolution and spatial coverage.

[0033] Water level data can be obtained through SWOT satellites. For example, it can be obtained from the L2_HR_RiverSP river products and L2_HR_LakeSP lake products published by SWOT satellites. These products contain information such as the water surface elevation of river nodes, river width, and lake water level observations. Their time resolution usually corresponds to the satellite's revisit period, such as a 21-day revisit period, and the water level measurement accuracy is usually better than 10 centimeters.

[0034] Total water storage data refers to data reflecting abnormal changes in terrestrial water storage within a target area, and it includes the sum of various components such as surface water, soil water, groundwater, and snow water equivalent. For example, total water storage data can be obtained by inverting data based on the time-varying characteristics of the Earth's gravity field observed by gravity satellites.

[0035] Total water storage data can be obtained through the GRACE satellite and its subsequent Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) satellites. Specifically, it can be obtained from the GRACE-FO Level-2 RL06 spherical harmonic coefficient product released by data centers such as the Center for Space Research (CSR), Jet Propulsion Laboratory (JPL), and the German Research Centre for Geosciences (GFZ). This product typically includes a Gravity Static Field Model (GSM) and a Global Atmosphere and Ocean De-aliasing Product (GAC), with a temporal resolution typically on a monthly scale and a spatial resolution of approximately 1 degree by 1 degree.

[0036] Optical imagery data of a target area refers to remote sensing images that can reflect land cover, extract water body boundaries, and identify the dynamic range of water bodies. For example, optical imagery data can be multispectral satellite remote sensing imagery.

[0037] Optical imagery data of the target area can be obtained by downloading publicly available satellite imagery products. Specifically, this could involve acquiring multi-spectral instrument (MSI) images from the Sentinel-2 satellite and operational land imager (OLI) images from the Landsat-8 or Landsat-9 satellites. These images will be used for subsequent dynamic identification and morphological feature parameter extraction for different types of water bodies.

[0038] Topographic data refers to Digital Elevation Model (DEM) data that describes the distribution of surface elevation and topographic relief characteristics of a target area. For example, topographic data can be a high-precision digital elevation model, either global or regional.

[0039] Topographic data can be obtained from publicly available mapping data sources, such as DEM data from the Shuttle Radar Topography Mission (SRTM) or data from the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM). This data is mainly used to assist in topographic analysis, reservoir capacity inversion, and determination of the inundation range of water bodies under topographic constraints.

[0040] As an optional implementation, the multi-source data also includes hydro-meteorological reanalysis data and groundwater monitoring well measured data, which will help to construct hydrophysical constraints and verify the inversion results.

[0041] Hydrometeorological reanalysis data refers to datasets of historical meteorological and hydrological variables generated by assimilating observational data through numerical weather prediction models. These datasets provide the boundary conditions required for water balance calculations. For example, hydrometeorological reanalysis data may include precipitation, evaporation, and temperature data acquired from the ECMWF Reanalysis v5 Land (ERA5-Land) product, and soil moisture, snow water equivalent, and surface runoff data acquired from the Noah land surface model product in the GlobalLand Data Assimilation System (GLDAS). This hydrometeorological reanalysis data will be used to quantify regional precipitation recharge, evapotranspiration loss, and changes in soil and snow cover water volume.

[0042] Groundwater monitoring well data refers to measured records directly obtained from a network of surface observation stations, reflecting the true water level state of groundwater aquifers. These records provide ground truth constraints for the inversion model. For example, groundwater monitoring well data can include data collected from the groundwater monitoring network within the study area, specifically including the geographical coordinates of the monitoring wells, the aquifer type at their location, high-frequency water level observation time series, and regional water yield parameters. This groundwater monitoring well data can correct grid-scale inversion results through point-scale observations, improving the accuracy and reliability of the inversion.

[0043] Step 120: Based on multi-source data, the surface water bodies in the target area are classified into multiple water body types; using the water level-storage conversion model corresponding to each water body type, the water level data is converted into a surface water storage observation sequence.

[0044] Surface water bodies within the target area refer to all surface water-covered areas that exist within the geographical scope of the target area and can be observed through remote sensing. For example, surface water bodies within the target area can include various forms of water objects such as naturally formed lakes, perennially flowing river systems, artificially constructed reservoirs, and temporary water bodies that appear with seasonal precipitation or flood events.

[0045] Water body type refers to the specific category classified according to the morphological and geometric characteristics, hydrological connectivity, and temporal dynamic changes of water bodies. For example, water body types can be classified into four main categories: lakes, rivers, reservoirs, and seasonal water bodies. This refined classification aims to address the problem that traditional methods neglect the differences in water bodies, so as to adopt the most physically consistent calculation strategies for water bodies of different forms.

[0046] A water level-storage conversion model refers to a mathematical model or physical formula that describes the quantitative mapping relationship between changes in water level height and changes in water storage for a specific water body type. For example, a water level-storage conversion model can be an adaptive lake model considering the nonlinear characteristics of lake basin morphology for lakes, a river storage model based on cross-sectional morphology index and river segment accumulation method for rivers, a reservoir capacity curve model incorporating a time-varying correction coefficient for reservoirs, and a volume calculation model that distinguishes between the main channel and the floodplain inundation range for seasonal water bodies.

[0047] Surface water storage observation series refers to time series data that reflects the storage status of specific water bodies at the time of satellite observation, calculated using a water level-storage conversion model combined with measured water level data. For example, a surface water storage observation series can be a non-equidistant sequence of surface water volume changes generated based on the observation time points of surface water and ocean topography satellites, with a time resolution corresponding to the satellite revisit period (e.g., 21 days). This surface water storage observation series accurately records the water volume status of each independent water body at the time of observation and will serve as the basic input data for subsequent time series reconstruction and spatiotemporal fusion.

[0048] Step 130: Perform temporal reconstruction and spatial aggregation on the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of temporal and spatial scales.

[0049] Temporal series reconstruction refers to the process of filling in sparse, non-uniformly spaced satellite observation data into a continuous, high-frequency data sequence using mathematical interpolation or physical model assimilation. For example, temporal series reconstruction can combine climate zone adaptive water balance differential equations with ensemble Kalman filtering algorithms to transform sparse observations of surface water and ocean topography satellites with a 21-day revisit cycle into a continuous daily-scale sequence of storage changes, and further calculate the monthly average value, thereby solving the problem of temporal discontinuity in satellite observations.

[0050] Spatial aggregation refers to the process of summarizing the changes in the reserves of geographically dispersed point-like or linear water bodies into a unified, regular grid cell according to certain physical rules. For example, spatial aggregation can be based on a hydrological response correlation mechanism, calculating structural weights including area factors and hydraulic connectivity factors, as well as response weights reflecting the synchronicity of reserve changes. Using these two-factor weights, the reserves of scattered lakes, reservoirs, and linear rivers are weighted and summed into a 1-degree multiplied-1-degree grid defined by the Gravity Recovery and Climate Experiment Satellite.

[0051] Spatiotemporal matching refers to processing remote sensing data from different sources to ensure consistency in sampling frequency and resolution in both time and space, facilitating algebraic operations or collaborative analysis. For example, spatiotemporal matching could ensure that processed surface water storage data is represented as a monthly average over time and as a 1-degree multiplied grid value spatially, thus perfectly aligning with the spatiotemporal specifications of gravity field spherical harmonic coefficient inversion products published by the Gravity Recovery and Climate Experiment Satellite, eliminating calculation errors caused by scale differences.

[0052] A gridded surface water storage sequence refers to a surface water storage dataset with a regular grid structure and continuous time, generated after the aforementioned temporal reconstruction and spatial aggregation processes. For example, a gridded surface water storage sequence can be a time series representing the total volume change of all surface water bodies within each 1-degree multiplied 1-degree grid cell relative to the reference surface in each month. This gridded surface water storage sequence will be input as a known surface water component into the subsequent groundwater storage inversion model to subtract the influence of surface water from the total storage.

[0053] Step 140: Construct a groundwater storage inversion model based on the gridded surface water storage sequence and total water storage data, solve the groundwater storage inversion model, and obtain the groundwater storage sequence of the target area.

[0054] The groundwater storage inversion model refers to a computational framework based on the physical mechanism of the interaction between surface water and groundwater, used to mathematically separate the groundwater storage component from total terrestrial water storage observations.

[0055] The groundwater storage sequence of a target area refers to a dataset obtained by solving a groundwater storage inversion model, reflecting the dynamic changes in aquifer water volume over time within the target area. For example, the groundwater storage sequence of a target area can be a time series of groundwater storage anomalies with a monthly temporal resolution and a grid-scale spatial resolution (e.g., 1 degree by 1 degree). This time series of groundwater storage anomalies can quantitatively characterize the groundwater surplus and deficit status of the target area during the study period, revealing the long-term groundwater deficit trend or seasonal fluctuation characteristics of the target area (e.g., the North China Plain or arid and semi-arid regions). Furthermore, based on the annual change rate threshold of this sequence, groundwater over-extraction areas and their severity classification can be accurately identified.

[0056] The groundwater storage synergistic inversion method integrating SWOT and GRACE provided by this invention achieves refined calculation of surface water storage and alignment of observation data in the spatiotemporal dimension by classifying surface water bodies and reconstructing observation data. This enables it to accurately reflect the actual changes of surface water under complex hydrological environments and realize high-precision separation and inversion of groundwater storage.

[0057] As an optional embodiment, after acquiring water level data, total water storage data, and multi-source data, it is also necessary to preprocess the data to ensure the quality, reliability, and spatiotemporal consistency of data from different sources.

[0058] First, SWOT analysis of water level data quality control and reference surface conversion were performed. Specifically, low-quality observation data affected by cloud contamination or geometric distortion were removed based on the product's built-in quality indicators. To unify the elevation datum, the elevation values ​​in the SWOT data were converted from the reference ellipsoid to the geoid defined by the Earth Gravitational Model 2008 (EGM2008). Simultaneously, an outlier detection method using an adaptive window based on water body type was employed to clean the water level data. A 90-day sliding window was used for lakes, and a 30-day sliding window was used for rivers. Anomalies and jumps in the time series were removed based on a 3-standard-deviation criterion, thus obtaining a continuous and smooth water level observation sequence.

[0059] Next, GRACE spherical harmonic coefficient processing and spatial filtering were performed. To correct low-order term errors, satellite laser ranging (SLR) data was used to replace the original C20 coefficients, and Degree-1 coefficients were added for geocentric motion correction. To suppress high-frequency noise and north-south stripe errors related to satellite orbits, a Gaussian filter with a 300 km radius combined with P4M6 decorrelation filtering was used. Furthermore, a scaling factor method was employed to correct signal leakage and recover signal attenuation caused by filtering. To reduce the systematic bias of a single data center, the arithmetic mean of the product results published by the three data centers (CSR, JPL, and GFZ) was taken as the final gravity field data used.

[0060] Optical remote sensing and hydrometeorological data were processed again. For optical remote sensing imagery, after atmospheric correction, the Modified Normalized Difference Water Index (MNDWI) was calculated, and a minimum cloud cover synthesis strategy was used to generate monthly cloudless imagery for subsequent water body boundary extraction. For ERA5-Land and GLDAS products, based on measured data from ground meteorological stations, the reanalysis data was biased using the quantile mapping method, and soil water storage changes were calculated based on the integrated corrected soil moisture data as a subsequent deduction item.

[0061] Next, the groundwater monitoring well data was standardized. The well coordinates of all monitoring wells were unified to the 1984 World Geodetic System (WGS84), and the elevation datum was unified to the 1985 National Elevation Datum. After removing outlier records and filling in missing values, based on the specific yield parameters of the regional aquifer, the measured water level changes were converted into equivalent water column height changes using a conversion formula. This conversion formula is: ;in, This represents the change in equivalent water column height. The aquifer yield; This is to monitor the actual water level changes in the monitoring well, thereby ensuring the consistency of physical dimensions between the ground-based measured data and the satellite-retrieved data.

[0062] Finally, spatiotemporal registration and standardized dataset construction were performed. A unified spatiotemporal framework was established using the spatial resolution of the GRACE data (i.e., 1 degree multiplied by 1 degree) as the base grid and the month as the basic time unit. Various types of water level, gravity, optical, and hydrometeorological data, processed as described above, were resampled or aggregated to be registered under this unified framework. The data was then stored in the Network Common Data Form - Climate and Forecast (NetCDF-CF) format to facilitate the reading and collaborative computation of multi-source data in subsequent steps.

[0063] Figure 2 This is a schematic diagram of the process for constructing a groundwater storage inversion model provided by the present invention, as shown below. Figure 2 As shown, as another optional embodiment provided by the present invention, a groundwater storage inversion model is constructed based on gridded surface water storage sequence and total water storage data, including but not limited to the following steps: Step 210: Construct the hysteresis transfer function of surface water and groundwater based on the gridded surface water storage sequence.

[0064] Considering the differences in geological structure, the hydrological response rates vary across different regions. Therefore, this invention sets the lag time for different zones based on the hydrogeological conditions of the study area. Specifically, for alluvial plains with good permeability but long infiltration paths, the lag time is set to 15 to 30 days; for hilly areas with significant topographic relief and slow infiltration processes, the lag time is set to 30 to 60 days; and for bedrock mountainous areas with well-developed fissures and rapid infiltration, the lag time is set to 60 to 90 days. Accurate lag time setting is fundamental to constructing a consistent physical process model.

[0065] Secondly, the exchange rate between surface water and groundwater is calculated based on the set lag time. In this process, a lag transfer function model is used to describe the dynamic relationship between the two, as shown below: ; in, express The amount of surface water exchanged with groundwater at any given time, expressed in millimeters of equivalent water column height. A positive value indicates that surface water is replenishing groundwater; it reflects the replenishment efficiency, and its preferred value range is 0.05 to 0.25. Changes in surface water storage at any given time, expressed in millimeters of equivalent water column height; The time lag is the factor. The physical meaning of this model is that the amount of exchange at the current moment depends on... The state of surface water storage before a certain time. When surface water storage is high, it replenishes groundwater; conversely, when surface water storage is low, it causes groundwater to drain into the surface.

[0066] Furthermore, for regions requiring consideration of more complex response processes, the aforementioned linear hysteresis model can be extended to an exponentially decaying transfer function model. In this model, the exchange volume at the current moment is represented as a weighted integral of the changes in surface water storage over a past period. The exponentially decaying transfer function model is shown below: ; in, express The amount of surface water exchanged with groundwater at any given time. Represents the commutation coefficient. The variable is the integral variable, representing the time delay; It is an exponentially decaying weighting function, which reflects the decay characteristics of the replenishment effect with time delay; The normalization coefficient is used to ensure that the integral of the weighting function over the integral domain is 1. The exponential decay transfer function model can more realistically reflect that during the process of rainfall or surface water infiltration, the replenishment effect is not concentrated at a specific moment, but rather exhibits a gradual decay over time.

[0067] Step 220: Determine the water balance constraint terms based on the hysteresis response transfer function and total water storage data.

[0068] Specifically, the water balance constraint is an optimization objective component constructed based on the closure principle of terrestrial water storage, aiming to ensure that the sum of the retrieved components physically approximates the total storage observed by satellite as closely as possible. For example, a least-squares form water balance objective function can be constructed, which calculates the sum of squared residuals between the observed total water storage and the sum of each independent water component at each time step. The specific calculation formula is as follows: ; in, This represents the value of the water balance constraint term. The smaller the value, the more the inversion result conforms to the water balance principle. t For time variables, This represents the change in total water storage observed by gravity satellites (such as GRACE). This represents the change in groundwater storage that needs to be solved, which is the core unknown in the groundwater storage inversion model. This represents the gridded changes in surface water storage. Changes in soil water storage are typically represented by land surface process models. Indicates the change in snow water equivalent. This indicates the amount of water exchanged between surface water and groundwater.

[0069] Step 230: Obtain measured groundwater level change data to determine groundwater level observation constraints based on the measured groundwater level change data.

[0070] Specifically, by fusing point-scale ground-based measured information with area-scale remote sensing inversion information, high-precision measured water level data is used to constrain the model solution process, preventing the inversion results from deviating from actual observation trends. For example, a groundwater level observation constraint term based on monitoring well data can be constructed. This constraint term measures the weighted difference between the change in groundwater storage in the area to be solved and the change in measured storage in a single well after yield conversion. The specific calculation formula is as follows: ; in, This represents the value of the groundwater level monitoring constraint term. t Represents a time variable. k Indicates the monitoring well number, This represents the change in groundwater storage to be solved. Indicates the first k Monitoring wells at t The measured change in water level at any given time. Indicates the first k The specific yield parameters of the aquifer at the location of the monitoring well are used to convert changes in water level height into equivalent changes in reservoir height. Indicates the first k The formula calculates the observation error variance of a monitoring well. By dividing by the observation error variance, a weighting mechanism is implemented: monitoring wells with smaller observation errors and higher data quality are assigned greater constraint weights, while those with smaller errors are assigned less weight.

[0071] Step 240: Using the groundwater storage sequence as the variable to be solved, a groundwater storage inversion model is constructed based on the water balance constraint, groundwater level observation constraint, time series smoothing constraint, and physical process constraint. The time series smoothing constraint is used to constrain the fluctuation range of groundwater storage changes at adjacent times, and the physical process constraint is used to constrain the increase and decrease of groundwater storage.

[0072] Specifically, the optimal groundwater storage change sequence is found by minimizing a joint optimization objective function (i.e., a groundwater storage inversion model) that includes four constraints. For example, the joint optimization objective function is shown below: ; in, To minimize the total value of the objective function, This is a water balance constraint term. This is a constraint term for groundwater level monitoring. For time series smoothing constraints, These are constraints on the physical process.

[0073] It should be noted that in the process of solving the joint optimization objective function, all observation constraints must be reasonably satisfied and no drastic fluctuations should occur. When the ratio of the minimum value of the joint optimization objective function to the initial value is less than 0.1, it can be considered that the convergence requirement is met.

[0074] The purpose of the time-series smoothing constraint is to utilize the inertial characteristics of the groundwater system to limit non-physical anomalous jumps in the storage change sequence. For example, the time-series smoothing constraint can be expressed as the smoothing regularization coefficient multiplied by the sum of squares of the differences in groundwater storage changes at adjacent time points. The calculation formula is as follows: ; in, Represents the time-series smoothing constraint term. To smooth out the regularization coefficients, and These represent the changes in groundwater storage at the current moment and the previous moment, respectively.

[0075] The purpose of the physical process constraint is to ensure that the retrieved recharge and discharge volumes do not exceed the physical limits allowed by the regional hydrogeological conditions. For example, the physical process constraint can be constructed as a function that penalizes the portion exceeding the maximum recharge and discharge limits, and its calculation formula is as follows: ; in, Represents physical process constraints. For physical constraint regularization coefficients, for t The increase in reserves at any given time, i.e., the amount of replenishment. for t The decrease in reserves at any given time, i.e., the amount of waste excreted. This is the maximum monthly supply limit. This represents the upper limit of the maximum monthly excretory capacity. This is the function that maximizes the value. The meaning of this formula is that when the absolute value of the supply or discharge exceeds a set physical threshold, the objective function value will increase significantly, thus forcing the solution to revert to a reasonable physical range.

[0076] The groundwater storage co-inversion method integrating SWOT and GRACE provided by this invention can significantly improve the physical consistency and reliability of groundwater storage separation by constructing a groundwater inversion model that includes multiple constraints such as hysteresis transfer function, water balance, measured water level, time series smoothing and physical processes. This ensures that the inverted recharge and discharge volumes are consistent with the actual carrying capacity of regional hydrogeology, thereby achieving deep coupling of complex hydrophysical mechanisms in the mathematical inversion process.

[0077] Figure 3 This is a schematic diagram of the process for temporal reconstruction and spatial aggregation of surface water storage observation sequences provided by the present invention, as shown below. Figure 3 As shown, as another optional embodiment provided by the present invention, temporal reconstruction and spatial aggregation are performed on the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of spatiotemporal scale, including but not limited to the following steps: Step 310: Establish a water balance differential equation based on the precipitation, inflow, outflow and evaporation of the target area.

[0078] Specifically, the principle of water balance driven by physics can be used to connect discrete satellite observation points in series, establishing a dynamic equation describing the continuous evolution of surface water storage over time. For example, a water balance differential equation of the following form can be established: ; in, The rate of change of surface water storage is usually obtained by differential calculation from the surface water storage observation series, and serves as an observational constraint for the water balance differential equation. For inbound traffic, For the outflow rate, measured flow data is used when there is a control section, and Manning's formula is used when there is no measured section. This refers to daily precipitation. This represents the actual evaporation rate. For water surface area, This represents the amount of water exchanged between surface water and groundwater. This refers to hydrological processes specific to a climate zone.

[0079] Inflow can be obtained through flow calculation, and the formula is as follows: ; in, For inbound traffic, For the first j Surface runoff depth of each sub-basin For the first j The area of ​​the sub-basin, For confluence delay.

[0080] Step 320: Input the surface water storage observation sequence into the water balance differential equation to obtain the daily-scale surface water storage sequence of each surface water body. Then, perform monthly average processing on the daily-scale surface water storage sequence to obtain the monthly averaged surface water storage of the surface water body.

[0081] Specifically, data assimilation techniques can be used to fuse sparse satellite observation data into a continuous physical model, thereby generating a high temporal resolution storage sequence and performing monthly aggregation. For example, at the observation time of the SWOT satellite, an ensemble Kalman filter algorithm can be used to assimilate the surface water storage observation sequence into the state variables of the water balance differential equation.

[0082] Assume the water balance differential equation model predicts the storage capacity at the current time as follows: The actual reserves observed by the satellite are The updated analytical reserves Kalman gain , The variance of the model prediction error. This represents the variance of the observation error. Taking into account both satellite observation errors and model conversion errors, the calculation formula is as follows: ,in This is due to SWOT error in water level observation. m, Sensitivity to water level-storage conversion. This represents the model error.

[0083] Through numerical solution and assimilation processes, the originally discrete-time observation data can be reconstructed into a continuous daily-scale reserve sequence. Furthermore, the daily-scale reserve changes can be calculated. Its calculation formula is ,in For reference time reserves (e.g., average value in the first month of the study period), This represents the area of ​​the corresponding grid. Finally, the monthly average reserve change is calculated by aggregating data over the calendar month. In the formula The number of days in the current month. For the first iThe water body was in the first of the month j The daily-scale changes in reserves were analyzed to obtain monthly averaged results consistent with the time scale of gravity satellite data.

[0084] Step 330: Divide the target area into multiple grid cells according to the spatial resolution of the total water storage data.

[0085] Specifically, to achieve spatial scale uniformity for data from different sources, this invention discretizes the study area into a regular grid system and constructs the water body topology within the grid. For example, the target area can be divided into several regular grid cells according to the spatial resolution of GRACE satellite data (typically approximately 1 degree by 1 degree). For large water bodies that cross grid boundaries, their reserves are allocated to adjacent grid cells according to their area proportion within each grid. Within each divided grid cell, a water body topology map describing the connectivity between water bodies is further constructed. ,in It is a collection of nodes (i.e., the individual water bodies within the grid). It is a set of edges (i.e., the hydraulic connections between water bodies).

[0086] For each water body node in the topology graph, calculate its hydraulic connectivity. This indicator is defined as the number of other water bodies directly connected to the water body, used to quantify the pivotal role of the water body in the water system. Furthermore, based on the monthly average storage change sequence obtained from the preceding steps, the time-series Pearson correlation coefficients between each pair of water body nodes within the grid are calculated. When the correlation coefficient is greater than a set threshold (e.g., 0.7), the water body is considered to be... and water bodies It has synchronous response characteristics.

[0087] Step 340: Within each grid cell, calculate the structural weight based on the area and hydraulic connectivity of each surface water body, and calculate the response weight based on the correlation between the changes in the storage of each surface water body and the overall storage changes of the grid.

[0088] Specifically, this invention quantifies the contribution of each water body to the overall storage change of its grid, considering both static geometric topology and dynamic hydrological response characteristics. For example, the structural weight can be calculated using the following formula: ; in, For the first i The structural weight of each water body For the first i The area of ​​each water body For hydraulic connectivity, The connectivity gain coefficient is preferably set between 0.1 and 0.3. This structural weighting design assigns greater weight to water bodies with high connectivity and those located at key points in the water system.

[0089] For the response weights, the average storage change within the grid is first calculated. ,in, N This represents the total number of water bodies within the grid. For the first i The water body t Changes in reserves over time.

[0090] Then the i Correlation coefficient between changes in water storage in individual bodies and the average value of the grid. Construct response weights , For the first j The correlation coefficient between the changes in water storage in individual bodies and the grid average value. By taking Ensure that water bodies that are negatively correlated with the overall trend (i.e., change in the opposite direction) do not participate in positive weighting, avoid interference from abnormal water bodies on the aggregation results, and ensure that the aggregation results reflect the dominant hydrological trend of the region.

[0091] Step 350: Based on structural weights and response weights, determine the aggregation weights of each surface water body, and perform weighted fusion based on the monthly average surface water storage of each surface water body and the aggregation weights to obtain a gridded surface water storage sequence.

[0092] Specifically, this invention adaptively fuses the two weights mentioned above to generate a final aggregated weight, and calculates the total reserves at the grid scale accordingly. For example, the formula for calculating the aggregated weight is: ,in, For the first i Aggregate weight of each surface water body For the first i Structural weights of individual surface water bodies For the first i Response weights of individual surface water bodies The fusion coefficient is the coefficient of fusion. The calculation formula is based on the adaptive determination of the uniformity of water spatial distribution within the grid. In the formula Let be the coefficient of variation of the water area within the grid, and , The standard deviation of the water area within the grid. This represents the average area of ​​the water body within the grid.

[0093] The physical significance of this mechanism lies in the fact that when the water body is relatively uniformly distributed ( When it is smaller, Approaching 0.7, the model emphasizes structural weights; when the water distribution is uneven ( When it is larger Approaching 0.3, the model emphasizes response weights. After determining the overall weight for each water body, a weighted aggregation calculation is performed: ,in, for t Gridded surface water storage at any given time For the first i Aggregate weight of individual water bodies For the first i The water body t Changes in reserves over time.

[0094] Finally, to ensure that the generated gridded surface water storage sequence is consistent with the inversion results from the Gravity Satellite (GRACE) in terms of spatial spectral characteristics, after completing the above weighted aggregation calculations, it is also necessary to perform spatial dimensional analysis. A Gaussian filter with a radius of 300 km is applied for smoothing to eliminate high-frequency noise and match the spatial smoothing characteristics of GRACE data, ultimately yielding a gridded surface water storage sequence that can be used for collaborative inversion.

[0095] The groundwater storage co-inversion method integrating SWOT and GRACE provided by this invention effectively solves the problem of mismatch between multi-source remote sensing data in terms of spatiotemporal scales by combining the water balance differential equation and the two-factor weighted aggregation model. Specifically, the water balance differential equation achieves the physical reconstruction from sparse observations to a continuous diurnal sequence, ensuring accurate synchronization in the temporal dimension. Meanwhile, the two-factor weighted aggregation strategy, which integrates area, hydraulic connectivity, and hydrological response correlation, overcomes the shortcomings of traditional simple area weighting methods that ignore the differences in spatial structure and dynamic response of water bodies, ensuring the physical rationality of aggregating point / linear water bodies to a grid scale. This provides high-precision spatiotemporally matched input data for groundwater inversion.

[0096] In another embodiment of the present invention, a water balance differential equation is established based on the precipitation, inflow, outflow, and evaporation of the target area, including: Climate zones are defined based on the annual precipitation of the target area, identifying arid and humid regions. The first evaporation model is used to determine evaporation in arid regions, while the second evaporation model is used to determine evaporation in humid regions.

[0097] Specifically, to more accurately simulate the water deficit process under different climatic backgrounds, the target region is first climatically zoned based on annual precipitation. For example, areas with annual precipitation less than 400 mm can be classified as arid zones, while areas with annual precipitation greater than or equal to 400 mm can be classified as humid zones. For arid zones, due to the significant influence of aerodynamic factors, the FAO Penman-Monteith model is used as the primary evaporation model to calculate actual evaporation. For humid zones, which are mainly controlled by energy, the Priestley-Taylor model is used as the secondary evaporation model to calculate actual evaporation. This differentiated calculation strategy can significantly improve the estimation accuracy of the evaporation term in the water balance equation.

[0098] If the target area has a snow cover period, a water balance differential equation is established based on the precipitation, inflow, outflow, snowmelt replenishment, and evaporation of the target area. If the target area belongs to a monsoon climate zone, a water balance differential equation is established based on the precipitation, inflow, outflow, seasonal adjustment, and evaporation of the target area.

[0099] Specifically, this invention introduces additional hydrological process terms for specific climate zones. To improve the water balance equation. For example, for cold regions with an average annual temperature below 0 degrees Celsius or a significant snow cover period, a snowmelt recharge term can be introduced as... Its calculation formula is in, The day factor is typically around 2.56 millimeters per degree Celsius per day. This is the critical temperature for snow melting, usually taken as 0 degrees Celsius. Indicates only when the temperature Higher than The difference is taken as a positive value if it is positive, otherwise it is 0. This represents the area covered by snow.

[0100] For monsoon regions significantly influenced by monsoons, a seasonal adjustment term is introduced as... Its calculation formula is In the formula, This is the monsoon response coefficient, which typically ranges from 0.1 to 0.3. This represents the current precipitation. This represents the average precipitation for the same period over many years. The watershed area is represented by this term. This seasonal moderating term reflects the nonlinear regulatory effect of monsoon precipitation anomalies on water storage.

[0101] As an optional implementation, after establishing the water balance differential equation containing the above terms, a fourth-order Runge-Kutta method is used for numerical integration to obtain a high-precision numerical solution, with a time step set to one day. During the solution process, special attention needs to be paid to the outflow rate. water surface area A It is not a fixed value, but rather a reserve. V The function reflects the dynamic physical relationship between water level, area and flow rate as storage changes, ensuring that the simulation process conforms to the laws of hydraulics.

[0102] The groundwater storage synergistic inversion method integrating SWOT and GRACE provided by this invention significantly improves the physical realism and environmental adaptability of surface water storage simulation by constructing climate-adaptive water balance differential equations. Specifically, differentiated evaporation models are used for arid and humid regions, effectively capturing the main controlling mechanisms of water loss under different climatic backgrounds. Simultaneously, snowmelt recharge terms are introduced for cold regions and seasonal adjustment terms are introduced for monsoon regions, accurately characterizing the specific contributions of snowmelt and monsoon precipitation anomalies to regional water storage. This allows for a detailed reconstruction of special hydrological processes under complex climatic conditions within the model, significantly reducing estimation errors caused by the lack of physical mechanisms.

[0103] In another embodiment of the present invention, solving the groundwater storage inversion model to obtain the groundwater storage sequence of the target area includes: obtaining the water balance calculation result and using the water balance calculation result as the current solution vector; iteratively executing the model optimization and update steps until the relative change of the solution vector obtained in two adjacent iterations is less than a preset convergence threshold, and determining the final output solution vector as the groundwater storage sequence of the target area.

[0104] The water balance calculation result refers to the change in groundwater storage estimated using the traditional subtraction principle, which is the remaining value after directly deducting the equivalent of surface water, soil water and snow water from the total terrestrial water storage.

[0105] The current solution vector refers to the vector composed of the decision variables to be solved in each iteration of the optimization algorithm. For example, the current solution vector could be a column vector composed of the monthly changes in groundwater storage during the study period, with its dimension equal to the total number of months in the study period. The preset convergence threshold refers to the numerical standard used to determine whether the optimization algorithm has reached a stable state. For example, the preset convergence threshold could be set to 10. -6 When the relative change of the solution vectors in two consecutive iterations (i.e., the ratio of norms) is less than this value, the algorithm is considered to have converged and the optimal solution has been obtained.

[0106] Specifically, this invention employs the L-BFGS-B (Limited-memory Broyden–Fletcher–Goldfarb–Shanno with Box constraints) algorithm to solve the aforementioned large-scale nonlinear optimization problem. This algorithm first initializes the solution vector using the results of the traditional water balance subtraction method, and then enters an iterative loop to continuously refine the solution vector to minimize the aforementioned quadruple-constraint objective function.

[0107] The model optimization and update steps include: calculating the gradient of the groundwater storage inversion model at the current solution vector; constructing an inverse approximation of the Hessian matrix using the gradient difference and displacement vector stored in historical iterations to determine the search direction; performing a line search along the search direction to determine the update step size; determining the updated solution vector based on the current solution vector, the search direction, and the update step size, and projecting the updated solution vector into a preset feasible region as the solution vector for the next iteration.

[0108] The gradient difference stored in historical iterations refers to the change in the gradient vector of the objective function over the most recent several iterations. For example, algorithms typically store the gradient difference from the most recent iterations. m The displacement vector and gradient difference vector are obtained from 10 iterations (typically 10) to approximate second-order curvature information.

[0109] The inverse approximation of the Hessian matrix refers to the implicit construction of the inverse matrix of the objective function's Hessian matrix using a double-loop recursive algorithm with limited historical gradient differences and displacement information. This avoids directly calculating and storing the huge Hessian matrix, greatly reducing computational complexity.

[0110] The search direction refers to the direction in which the objective function value decreases the most at the current iteration point. For example, the search direction... p (n) It is calculated by left-multiplying the negative value of the current gradient vector by the inverse approximation of the Hessian matrix.

[0111] Line search along a search direction refers to finding a suitable step size factor along a defined search direction so that the objective function value decreases to a certain extent in that direction (e.g., satisfying the Wolf condition). For example, the line search process is used to determine the scalar step size. α (n) .

[0112] The update step size refers to the distance the solution vector moves along the search direction in one iteration. The updated solution vector can be determined using vector addition; for example, the initial update solution equals the current solution vector plus the step size multiplied by the search direction vector.

[0113] Projecting to a predefined feasible region means forcibly restricting each component of the initially updated solution vector to a physically permissible upper and lower bound. For example, using a projection operator... P Constrain each element of the solution vector to a lower bound. x lb and the Upper Realm x ub Between these two boundaries, the physical range of changes in historical groundwater storage in the region (such as historical extreme values) is usually determined to ensure that the inversion results do not contain non-physical numerical overflows.

[0114] Specifically, the iterative formula can be expressed as: ;in, For the first n The solution vector after +1 iterations For projection operators, x lb The lower realm x ub For the upper realm, For the first n The solution vector of the nth iteration. For the first n The update step size for each iteration. For the first n The search direction for the next iteration. The optimal solution vector obtained after the algorithm converges. x ∗ The individual components represent the changes in groundwater storage for each month. ,in, for t Changes in groundwater storage at any given time. The first in the optimal solution vector t Each component constitutes a complete monthly time series.

[0115] Furthermore, to further improve accuracy, this result can be substituted into the hysteresis transfer function to update the surface water-groundwater exchange rate, and fed back into the model building steps to form an outer iteration until the overall system reaches stable convergence. Finally, the linear error propagation formula can also be used: ,in, for t Variance of groundwater storage estimation at any given time. This represents the change in groundwater storage. For each input component, To account for the standard error, the uncertainty of the groundwater storage estimate is calculated, and a 95% confidence interval is given accordingly to assess the reliability of the inversion results.

[0116] The groundwater storage collaborative inversion method integrating SWOT and GRACE provided by this invention can efficiently and stably handle large-scale nonlinear inversion problems by adopting an iterative optimization solution strategy. It ensures that the solution vector after each iteration is always within the feasible domain that conforms to the regional hydrophysical characteristics, effectively avoiding the risk of divergence and non-physical numerical overflow in the mathematical solution process, thereby ensuring the numerical convergence and physical rationality of the groundwater storage inversion results.

[0117] In another embodiment of the present invention, the surface water bodies in the target area are classified based on multi-source data, including: extracting the boundaries of each surface water body in the target area based on optical image data, and determining the area, aspect ratio, shape index and stability coefficient of each surface water body based on the boundaries of each surface water body.

[0118] The boundaries of surface water bodies refer to the geographical boundaries in remote sensing imagery that distinguish water bodies from non-water-covered types, defining the spatial extent of the water body. The boundaries of surface water bodies can be obtained by thresholding or edge detection extraction from preprocessed optical remote sensing images (such as MNDWI water index images).

[0119] The area of ​​a surface water body refers to the planar area of ​​the enclosed region bounded by the water body's boundary.

[0120] Aspect Ratio It refers to the ratio of the major axis length to the minor axis length of the smallest circumscribed rectangle of a water body, used to reflect the slenderness of the water body.

[0121] The shape index is a dimensionless parameter characterizing the complexity of a water body's boundary. Its calculation formula is typically the water body's circumference divided by the product of the square root of twice pi and the square root of its area. In the formula Shape index (circular water body) ); The circumference of the water body; The area of ​​the water body.

[0122] The stability coefficient refers to the coefficient of variation of a multi-temporal water body area sequence. It reflects the degree of drastic change in the water body over time. The formula is the standard deviation of the multi-temporal water body area divided by the average area. In the formula The stability coefficient; Standard deviation of water body area in multiple time phases; This is the area average. This stability coefficient addresses the problem of seasonal water body misclassification caused by traditional methods relying solely on single-phase imagery.

[0123] When the length-to-width ratio of a surface water body is greater than the first threshold and the shape index is greater than the second threshold, the surface water body is identified as a river; when the stability coefficient of a surface water body is greater than the third threshold, the surface water body is identified as a seasonal water body.

[0124] The first threshold refers to the aspect ratio limit used to distinguish rivers from other bodies of water; for example, it can be set to 10.

[0125] The second threshold refers to the shape index limit used to help confirm the characteristics of the river; for example, it can be set to 2.5.

[0126] The third threshold refers to the coefficient of variation limit used to identify seasonal water bodies; for example, it can be set to 0.4.

[0127] Seasonal water bodies refer to those whose area fluctuates significantly with the seasons and are not stable throughout the year, such as seasonal lakes or floodplains. Specifically, when a water body is detected that has an extremely high aspect ratio and a complex boundary shape ( and When a body of water is identified as a river, it should be classified as such; for non-river bodies of water, if their area fluctuates significantly over time (…), it should be classified as a river. If the water body is seasonal, it is classified as such.

[0128] Based on the geological data of the target area, the topographic features of the dam body are obtained. When the elevation change gradient of the dam body's topographic features is greater than the fourth threshold and the extension length is greater than the fifth threshold, the water body is identified as a reservoir; otherwise, the water body is classified as a lake.

[0129] Geological data for the target area refers to a dataset containing surface elevation information, such as a high-precision digital elevation model (DEM). Geological data for the target area can be obtained from publicly available data sources such as SRTM DEM or ASTER GDEM.

[0130] The topographic features of a dam body refer to the unique topographic structure formed by artificial dam construction, typically manifested as abrupt changes in local elevation and linear extension. The elevation gradient of the dam body's topographic features refers to the rate of change of slope along a direction perpendicular to the dam's axis at the dam's location.

[0131] The fourth threshold refers to the gradient limit used to identify abrupt changes in dam elevation; for example, it can be set to 15 degrees. The extension length refers to the continuous length of the elevation change feature in the horizontal direction.

[0132] The fifth threshold refers to the length limit used to confirm the size of the dam body; for example, it could be set at 100 meters. Specifically, for dams with high stability ( For enclosed water bodies, further analysis using DEM is conducted to determine whether there are obvious artificial dam features at their edges. If there are topographic features that meet the conditions (elevation gradient greater than 15 degrees and extension length greater than 100 meters), they are identified as reservoirs; otherwise, they are identified as natural lakes.

[0133] Considering the significant differences in morphological characteristics and hydrological behavior among different types of water bodies, this invention establishes differentiated water level-storage conversion models for lakes, rivers, reservoirs, and seasonal water bodies after classification.

[0134] The groundwater storage synergistic inversion method integrating SWOT and GRACE provided by this invention achieves automated and accurate identification of lakes, rivers, reservoirs and seasonal water bodies by constructing a multi-dimensional hierarchical classification strategy that includes aspect ratio, shape index, stability coefficient and dam topographic features. This refined classification lays a solid foundation for the subsequent adoption of differentiated physical transformation models for different water body types, thereby significantly improving the targeting and accuracy of surface water storage inversion.

[0135] In another embodiment provided by the present invention, the water level-storage conversion models corresponding to each water body type include: lake storage conversion model, river storage conversion model, reservoir storage conversion model, and seasonal water body storage conversion model.

[0136] Specifically, this invention employs specialized calculation models tailored to the hydrophysical characteristics of different water body types. For lakes, a transformation model considering adaptive lake basin shape coefficients is established: ; in, This represents the change in lake reserves. Baseline water level The corresponding water surface area; This refers to the water level change observed using SWOT analysis. This is a water level-dependent adaptive lake basin shape coefficient, reflecting the transitional characteristics of lake basin morphology from shallow to deep water areas.

[0137] The formula for calculating the adaptive lake basin shape coefficient is as follows: ; in, The rate of change of water surface area with respect to water level, at the reference water level. Take the value at that location. For water surface area, H For water level, Using the baseline water level, this adaptive basin shape coefficient characterizes the degree of inclination of the basin sidewalls, reflecting the transitional characteristics of the basin morphology from gentle slopes in shallow water to steep banks in deep water.

[0138] For example, for rivers, a river water storage conversion model is used, which discretizes the river along its course into several small river segments, and then uses the segment accumulation method to calculate the total storage. The calculation formula is as follows: ; in, Total river reserves; For the first River section length; For the first Dynamic river width, to which river level depends; For the first Average water depth of the river section.

[0139] To accurately describe the nonlinear change in river width with water level, the dynamic river width function is expressed as: ; in, For the first i The dynamic width of the river section depends on the water level. H For water level, As the benchmark water level, Baseline water level The river below is wide; The coefficient for variation in river width. The river channel cross-sectional morphology index, among which Indicates a compound cross-section. Represents a rectangular cross-section. express V shaped canyon.

[0140] For reservoirs with available design data, the reservoir capacity curve is used directly. Introduce a time-varying correction factor for storage capacity: ; in, For the corrected storage capacity; The time-varying correction factor (preferably 0.85~1.0) reflects the reservoir capacity loss caused by siltation; To determine the reservoir capacity curve values, for reservoirs without design data, the reservoir capacity curve is constructed based on the DEM using the horizontal cutting method.

[0141] Considering the dynamic changes in the boundaries of seasonal water bodies (such as seasonal lakes and floodplains) due to drastic fluctuations in water level, the water body is divided into two parts: a relatively stable main channel and a drastically changing floodplain. ; in, The total storage of seasonal water bodies. Main tank storage capacity; To determine the storage capacity of the beach area, the effective water storage boundary is identified based on real-time water level dynamics, and the inundation volume is calculated under DEM constraints.

[0142] The groundwater storage synergistic inversion method integrating SWOT and GRACE provided by this invention establishes differentiated water level-storage conversion models for lakes, rivers, reservoirs, and seasonal water bodies. This method can accurately characterize the volume response patterns of different types of water bodies under water level changes, thereby overcoming the problem of rough estimation caused by the use of uniform empirical formulas in traditional methods. It significantly improves the calculation accuracy of surface water storage observation sequences and provides a reliable data foundation for achieving high-precision separation of surface water and groundwater storage in the future.

[0143] As an optional embodiment, after obtaining the groundwater storage sequence of the target area, the groundwater over-extraction area can be further identified, classified, and its accuracy verified based on the sequence.

[0144] First, three threshold levels were set for the annual rate of change of groundwater storage. Based on historical monitoring data from typical over-extraction areas such as the North China Plain and the Ganges River Basin in India, the first-level threshold was set at an annual decline rate exceeding 3 mm, the second-level threshold at exceeding 10 mm per year, and the third-level threshold at exceeding 20 mm per year.

[0145] Based on this, an over-extraction zone determination criterion is established: when the average annual decline rate of a certain grid cell exceeds the first-level threshold (3 mm per year), it is marked as a potential over-extraction zone. Furthermore, when the groundwater storage in the potential over-extraction zone continues to decline for T_min years (T_min≥3 years) and the cumulative decline exceeds a preset threshold, it is confirmed as an over-extraction zone.

[0146] Based on the magnitude of the cumulative drop, the severity of over-extraction is divided into three levels: a cumulative drop of less than 150 mm (approximately corresponding to a cumulative drop of 4.5 meters in water level) is considered a slightly over-extracted area; a cumulative drop between 150 and 300 mm is considered a moderately over-extracted area; and a cumulative drop exceeding 300 mm is considered a severely over-extracted area.

[0147] Secondly, the accuracy of the inversion results is verified to ensure the reliability of the over-extraction zone determination. A trend consistency index is calculated using cross-validation of monitoring well data. ; in, Percentage of trend consistency; The number of monitoring wells whose inversion results are consistent with the trend of water level changes in the monitoring wells; This represents the total number of monitoring wells participating in the verification.

[0148] The quantitative definition of a steady state is: when the annual water level change in the monitoring well meets the following conditions... When it is determined to be a stable state, among which To stabilize the threshold, The absolute value of the annual water level variation in the monitoring well is determined based on 0.5 times the standard deviation of the regional multi-year water level fluctuation, generally ranging from 0.3 to 1.0 m. The stability determination of the inversion results is based on: .in, This represents the absolute value of the annual variation in groundwater storage obtained through inversion. For the aquifer yield, when And Pearson correlation coefficient At that time, the determination of the over-extraction area was confirmed to be valid.

[0149] It should be noted that the groundwater storage synergistic inversion method integrating SWOT and GRACE provided by this invention is applicable to various application scenarios such as regional groundwater resource dynamic assessment, over-extraction area identification and early warning, transboundary watershed water resource management, and hydrological monitoring in arid and semi-arid regions. Compared with existing technologies, this invention mainly has the following beneficial effects: First, this invention establishes differentiated water level-storage conversion models for different water body types and combines them with a two-factor spatial aggregation method, effectively overcoming the problem of large errors caused by the use of uniform empirical formulas in traditional methods. Specifically, this invention innovatively introduces hydraulic connectivity as a correction weight, comprehensively considering the area and connectivity of surface water bodies to construct a two-factor aggregation weight, thereby achieving a physically reasonable aggregation of point-like and linear water bodies to a grid scale.

[0150] Second, this invention proposes a physically constrained time-series reconstruction algorithm, which achieves synchronous fusion of water level data and total water storage data on a time scale, overcoming the shortcomings of existing technologies that rely solely on mathematical interpolation while neglecting hydrological physical mechanisms. Specifically, this invention establishes a physically driven time-synchronization fusion mechanism between the 21-day revisit cycle of satellites and the monthly scale of gravity satellite data, enabling the reconstruction of sparse observation data into a continuous daily sequence and subsequent monthly averaging.

[0151] Third, this invention constructs a collaborative separation model for surface water and groundwater, achieving high-precision inversion considering the physical correlation of hydrological processes. This invention incorporates the hysteresis response mechanisms of surface water and groundwater, as well as regional hydrogeological parameters, into the collaborative separation framework, establishing coupled dynamic equations that consider their interaction. Simultaneously, it integrates four constraints—water balance, monitoring well observations, temporal smoothness, and physical processes—into a joint optimization objective function, overcoming the limitations of traditional simple linear separation methods. Applying this model improves the accuracy of groundwater storage inversion by more than 12% and effectively identifies groundwater over-extraction areas.

[0152] Figure 4 This is a schematic diagram of the groundwater storage synergistic inversion system integrating SWOT and GRACE provided by the present invention, as shown below. Figure 4As shown, it mainly includes, but is not limited to: The data acquisition module 410 is used to acquire water level data, total water storage data and multi-source data of the target area; the multi-source data includes optical image data and topographic data of the target area.

[0153] The classification and conversion module 420 is used to classify surface water bodies in the target area based on multi-source data, resulting in multiple water body types; and to convert water level data into surface water storage observation sequences using water level storage conversion models corresponding to each water body type.

[0154] The spatiotemporal reconstruction module 430 is used to perform temporal reconstruction and spatial aggregation of the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of spatiotemporal scale.

[0155] The collaborative inversion module 440 is used to construct a groundwater storage inversion model based on gridded surface water storage sequence and total water storage data, solve the groundwater storage inversion model, and obtain the groundwater storage sequence of the target area.

[0156] It should be noted that the groundwater storage synergistic inversion system integrating SWOT and GRACE provided by the present invention can execute the groundwater storage synergistic inversion method integrating SWOT and GRACE as described in any of the above embodiments during specific operation, which will not be elaborated in this embodiment.

[0157] The groundwater storage collaborative inversion system integrating SWOT and GRACE provided by this invention achieves refined calculation of surface water storage and alignment of observation data in the spatiotemporal dimension by classifying surface water bodies by type and reconstructing observation data in the spatiotemporal dimension. This enables it to accurately reflect the actual changes of surface water under complex hydrological environments and realize high-precision separation and inversion of groundwater storage.

[0158] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for synergistic inversion of groundwater storage integrating SWOT and GRACE, characterized in that, include: Acquire water level data, total water storage data, and multi-source data for the target area; the multi-source data includes optical image data and topographic data of the target area. Based on the multi-source data, the surface water bodies in the target area are classified into multiple water body types; using the water level-storage conversion model corresponding to each water body type, the water level data is converted into a surface water storage observation sequence. Temporal reconstruction and spatial aggregation are performed on the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of temporal and spatial scale; Based on the gridded surface water storage sequence and the total water storage data, a groundwater storage inversion model is constructed. The groundwater storage inversion model is solved to obtain the groundwater storage sequence of the target area.

2. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 1, characterized in that, The construction of a groundwater storage inversion model based on the gridded surface water storage sequence and the total water storage data includes: Based on the gridded surface water storage sequence, the hysteresis transfer function of surface water and groundwater is constructed; Based on the hysteresis response transfer function and the total water storage data, water balance constraints are determined. Obtain measured data on groundwater level changes, and determine groundwater level observation constraints based on the measured data on groundwater level changes; Using the groundwater storage sequence as the variable to be solved, the groundwater storage inversion model is constructed based on the water balance constraint, the groundwater level observation constraint, the time series smoothing constraint, and the physical process constraint. The temporal smoothing constraint term is used to constrain the fluctuation range of groundwater storage changes at adjacent time points, and the physical process constraint term is used to constrain the increase and decrease of groundwater storage.

3. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 2, characterized in that, The hysteresis response transfer function is shown in the following mathematical model: ； in, express The amount of surface water exchanged with groundwater at any given time. Represents the commutation coefficient. Changes in surface water storage over time. This refers to the lag time.

4. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 3, characterized in that, The water balance constraint term is calculated based on the following mathematical model: ； in, This is a water balance constraint term. t For time variables, The total water storage data is as follows. This is the groundwater storage sequence. For the gridded surface water storage sequence, For changes in soil water storage, For changes in snow water equivalent, This represents the amount of water exchanged between surface water and groundwater.

5. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 1, characterized in that, The process of temporal reconstruction and spatial aggregation of the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of temporal and spatial scale includes: A water balance differential equation is established based on the precipitation, inflow, outflow and evaporation of the target area; The surface water storage observation sequence is input into the water balance differential equation to obtain the daily-scale surface water storage sequence of each surface water body. The daily-scale surface water storage sequence is then averaged monthly to obtain the monthly averaged surface water storage of the surface water body. The target area is divided into multiple grid cells based on the spatial resolution of the total water storage data; Within each of the grid cells, structural weights are calculated based on the area and hydraulic connectivity of each of the surface water bodies, and response weights are calculated based on the correlation between the changes in the storage of each of the surface water bodies and the overall storage changes of the grid. Based on the structural weights and the response weights, the aggregate weights of each surface water body are determined. The monthly average surface water storage of each surface water body and the aggregate weights are then weighted and fused to obtain the gridded surface water storage sequence.

6. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 5, characterized in that, The determination of the aggregate weight of each surface water body is based on the following mathematical model: ； in, For the first i Aggregate weight of each surface water body For the first i Structural weights of individual surface water bodies For the first i Response weights of individual surface water bodies The fusion coefficient is the coefficient of fusion. , is the coefficient of variation of the water area within the grid cell.

7. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 5, characterized in that, The establishment of a water balance differential equation based on the precipitation, inflow, outflow, and evaporation of the target area includes: Climate zones are defined based on the annual precipitation of the target area, identifying arid and humid zones. The evaporation rate is determined using a first evaporation model for the arid region and a second evaporation model for the humid region. If the target area has a snow cover period, a water balance differential equation is established based on the precipitation, inflow, outflow, snowmelt replenishment, and evaporation of the target area. If the target area belongs to the monsoon climate zone, a water balance differential equation is established based on the precipitation, inflow, outflow, seasonal adjustment term, and evaporation of the target area.

8. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 1, characterized in that, The process of solving the groundwater storage inversion model to obtain the groundwater storage sequence of the target area includes: Obtain the water balance calculation results and use them as the current solution vector; The model optimization and update steps are iteratively executed until the relative change of the solution vectors obtained in two adjacent iterations is less than the preset convergence threshold. The final output solution vector is then determined as the groundwater storage sequence of the target area. The model optimization and update steps include: Calculate the gradient of the groundwater storage inversion model at the current solution vector; The inverse approximation of the Hessian matrix is ​​constructed using the gradient difference and displacement vector stored in the historical iterations to determine the search direction; Perform a line search along the search direction to determine the update step size; Based on the current solution vector, the search direction, and the update step size, an updated solution vector is determined, and the updated solution vector is projected into a preset feasible region as the solution vector for the next iteration.

9. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 1, characterized in that, The classification of surface water bodies within the target area based on the multi-source data includes: Based on the optical image data, the boundaries of each surface water body in the target area are extracted, and the area, aspect ratio, shape index and stability coefficient of each surface water body are determined based on the boundaries of each surface water body. When the aspect ratio of the surface water body is greater than a first threshold and the shape index is greater than a second threshold, the surface water body is identified as a river. When the stability coefficient of the surface water body is greater than the third threshold, the surface water body is identified as a seasonal water body; Based on the geological data of the target area, the topographic features of the dam body are obtained. When the elevation change gradient of the topographic features of the dam body is greater than the fourth threshold and the extension length is greater than the fifth threshold, the water body is identified as a reservoir; otherwise, the water body is classified as a lake.

10. The groundwater storage synergistic inversion method integrating SWOT and GRACE as described in claim 9, characterized in that, The water level-storage conversion models corresponding to each of the aforementioned water body types include: lake storage conversion model, river storage conversion model, reservoir storage conversion model, and seasonal water body storage conversion model; The lake storage conversion model is shown in the following mathematical model: ； in, This represents the change in lake reserves. As the benchmark water level, Baseline water level The corresponding water surface area This refers to the change in water level. For adaptive lake basin shape coefficient; The river storage conversion model is shown in the following mathematical model: ； in, Total river reserves; For the first River section length; For the first The width of the river section; For the first Average water depth of the river section.

11. A groundwater storage synergistic inversion system integrating SWOT and GRACE, characterized in that, include: The data acquisition module is used to acquire water level data, total water storage data, and multi-source data of the target area; the multi-source data includes optical image data and topographic data of the target area. The classification and conversion module is used to classify the surface water bodies in the target area based on the multi-source data to obtain multiple water body types; and to convert the water level data into a surface water storage observation sequence using a water level storage conversion model corresponding to each water body type. The spatiotemporal reconstruction module is used to perform temporal reconstruction and spatial aggregation on the surface water storage observation sequence to obtain a gridded surface water storage sequence that matches the total water storage data in terms of spatiotemporal scale. The collaborative inversion module is used to construct a groundwater storage inversion model based on the gridded surface water storage sequence and the total water storage data, solve the groundwater storage inversion model, and obtain the groundwater storage sequence of the target area.

Images