A method for obtaining GWSA by a groundwater downscaling model driven by input and output variables of a homologous hydrological model and its application in ecological restoration

Abstract

The present invention belongs to the cross - technical fields of satellite gravimetry, hydrology, etc., and particularly relates to a method for obtaining GWSA by a groundwater downscaling model driven by input and output variables of a homologous hydrological model and its application in ecological restoration. The method comprises the following steps: (1) screening the input and output variables from a set of hydrological models through a PLSR model in combination with satellite observation data to obtain homologous groundwater downscaling model prediction variables; (2) constructing a non - stationary relationship between the prediction variables and the target variables at a relatively coarse spatial resolution through a GWR model; (3) using fine - scale predictors and the established downscaling model to divide the coarse - pixel terrestrial water storage into fine - pixels; (4) separating the contributions of other water storage information from the terrestrial water storage information in each high - precision grid cell, so as to obtain the GWSA result at a high spatial resolution. The method of the present invention improves the accuracy and credibility of the model output results.

Description

Technical Field

[0001] The present invention belongs to the cross - technical fields of satellite gravity, hydrology, etc., and particularly relates to a method for obtaining GWSA by a groundwater downscaling model driven by input and output variables of a homologous hydrological model and its application in ecological restoration. Background Technique

[0002] Under the background of global change, climate warming and other human activities are changing the global water cycle at an unprecedented speed. As the most active component in the hydrological cycle, accurately quantifying the change of groundwater level (Groundwater Storage Anomaly, GWSA) at the spatio - temporal scale is of great significance for groundwater governance and management.

[0003] Currently, the methods for groundwater monitoring and assessment mainly rely on groundwater monitoring networks, groundwater models or remote sensing observations. However, as a hidden resource, groundwater has the lack of large - scale, long - time - series, high - spatial - resolution and reliable - accuracy groundwater information due to its dynamics and complexity. There are serious deficiencies in the ground - observed groundwater level data in terms of time and space: (i) The groundwater monitoring data or investigation reports at the point - scale are often not publicly available, and only the observed statistical values at the provincial / municipal scale can be obtained; (ii) Affected by the number of groundwater well monitoring points, the spatial distribution of observation points is uneven, and the groundwater change has high spatial heterogeneity. The groundwater level dynamics observed from underground wells can be used as a quantitative index for local groundwater change, but extensive groundwater monitoring activities require a dense well network covering the entire area, with huge costs. In addition, due to the complexity of the geological structure of the underground aquifer and the groundwater recharge / discharge process, there are obvious limitations and uncertainties in groundwater model simulation.

[0004] In recent years, although some hydrological models have added groundwater modules, since hydrological models mainly focus on the large - scale water cycle of the natural water cycle or land - atmosphere coupling, the description of the internal mechanism and process of the underground aquifer in the module is relatively simple, and the quality of groundwater simulation cannot meet the requirements.

[0005] Since the launch of satellites, the GRACE (Gravity Recovery and Climate Experiment) satellite and its subsequent GRACE - FO satellite have provided a brand - new research perspective for groundwater quantification assessment with unprecedented accuracy and resolution. However, due to the characteristics of the satellite structure and orbital configuration, GRACE only measures relatively large areas (generally > 90,000 km 2sensitive to the Terrestrial Water Storage Anomaly (TWSA). The coarse spatial resolution of GRACE observations limits its application at the medium and small regional scales.

[0006] In recent years, some scholars have used the method of downscaling the target variable to improve the spatial resolution of GRACE products. Although the selection of predictor variables plays a crucial role in determining the quality and performance of the downscaling model, so far, the sources of the input variables of existing downscaling models are often multi-source, and there is no unified standard for the selection of input variables for downscaling models. Existing downscaling model studies only focus on improving the downscaling algorithm (such as Chinese Patent CN116383590A), and the exploration of the selection of input variables for downscaling models and the water budget closure of hydrological modeling are often ignored. The heterogeneity and inconsistency between different source data will lead to deviations in the total water volume and distribution, thus significantly exacerbating the problem of non-closure of the water balance in grid cells and ultimately affecting the reliability of model output. Therefore, implementing strict input data quality control during the model construction stage and maintaining its consistency are of great significance for improving the accuracy and credibility of model output results. Summary of the Invention

[0007] Aiming at the defects in the prior art, the purpose of the present invention is to provide a method for obtaining GWSA using a groundwater downscaling model driven by input and output variables of a homologous hydrological model, so as to reduce the uncertainty of the previous groundwater downscaling model and improve the accuracy and credibility of model output results.

[0008] The present invention also provides the application of the above method in the monitoring and identification of groundwater levels in ecological restoration areas affected by human activities and with obvious changes in groundwater levels compared to the original state, which is of crucial significance for guiding the optimization and adjustment of ecological restoration work and subsequent engineering measures and improving the ecological restoration effect.

[0009] The present invention adopts the following technical solutions:

[0010] A method for obtaining GWSA using a groundwater downscaling model driven by input and output variables of a homologous hydrological model, comprising the following steps:

[0011] Step 1: Obtaining input and output variables of the homologous hydrological model

[0012] Based on a distributed hydrological model, combined with multi-source remote sensing data and reanalysis data, a hydrological model for simulating the hydrological situation of the area to be studied is constructed; high-spatial-resolution climate data and vegetation data are used as model driving input variables to drive the hydrological model, and high-spatial-resolution model output hydrological variables are obtained;

[0013] Step 2: Screening of variables for the downscaling model prediction

[0014] Based on the PLSR model (Partial Least Squares Regression model), using the monthly TWSA (Terrestrial Water Storage Anomaly) data of the area to be studied as the target variable, variables with VIP value (Variable Importance in Projection) > 0.8 are screened from the model-driven input variables and model output hydrological variables described in Step 1 as the prediction variables for the downscaling model.

[0015] Step 3: Construction of the groundwater downscaling model driven by the input and output variables of the homologous hydrological model

[0016] Step 3-1: Sample the prediction variables obtained in Step 2 from high resolution to low spatial resolution to be consistent with the spatial resolution of the monthly TWSA data of the area to be studied;

[0017] For example, if the original spatial resolution of the TWSA data provided by the GRACE data product is 0.5°, for the sake of uniformity, the spatial resolution of the prediction variables needs to be sampled from high resolution to 0.5°;

[0018] Step 3-2: Using the prediction variables obtained in Step 2 as independent variables and the monthly TWSA data of the area to be studied as the dependent variable, establish a GWR model (Geographically Weighted Regression model), and according to the regression relationship, obtain the geographically weighted regression coefficients and regression residuals at low spatial resolution;

[0019] Step 3-3: Convert the geographically weighted regression coefficients and regression residuals at low spatial resolution to high spatial resolution, and substitute the geographically weighted regression coefficients and prediction variables at high spatial resolution into the above GWR model to obtain the TWSA at high spatial resolution.

[0020] Step 4: Obtaining the high-resolution GWSA

[0021] Exclude SMSA and SWSA from the high-spatial-resolution TWSA results to obtain the high-spatial-resolution GWSA results.

[0022] Vertically, the main contributors to TWSA are Soil Moisture Storage Anomaly (SMSA), Snow Water Storage Anomaly (SWSA), and Groundwater Storage Anomaly (GWSA). Therefore, excluding SMSA and SWSA from the high-spatial-resolution TWSA results can obtain the high-spatial-resolution GWSA results.

[0023] Furthermore, the distributed hydrological model described in step 1 is the ESSI series hydrological model, SWAT hydrological model or MIKE-SHE hydrological model; preferably the ESSI-3 model.

[0024] Furthermore, the climate data described in step 1 includes precipitation, average temperature, surface temperature, wind speed, snow cover, heat flux; the vegetation data includes normalized vegetation index, leaf area index, soil type, lithology, hydraulic conductivity; the hydrological variables include soil moisture, snow water equivalent, surface runoff, baseflow, actual evapotranspiration, vegetation root zone moisture.

[0025] Furthermore, the climate data described in step 1 is sourced from the CMFD, a dataset of surface meteorological elements driving data for the Chinese region provided by the National Tibetan Plateau Data Center of China. The CMFD dataset has a time resolution of 3 hours and a horizontal spatial resolution of 0.1°; the vegetation data is sourced from global LAI products.

[0026] Furthermore, the high-spatial-resolution climate data and vegetation data in step 1 are obtained by resampling the original low-spatial-resolution climate data and vegetation data through bilinear interpolation.

[0027] Furthermore, the multi-source remote sensing data in step 1 refers to various types of data obtained from different remote sensing sensors or platforms. For example, satellite imagery data can provide information on surface features (such as land use, vegetation cover, etc.) and soil moisture. Combining these remote sensing data with the hydrological model can input more accurate and detailed relevant parameters or state information into the model;

[0028] The reanalysis data is usually a large-scale collection of processed and integrated meteorological, climate, etc. data. These data can provide long-term historical data and spatio-temporal distribution information of meteorological elements (such as temperature, precipitation, wind speed, etc.), and are also incorporated into the model to better simulate the relationship between hydrological processes and climate conditions.

[0029] Furthermore, with R 2 、Q 2 and Q 2 cum values all greater than 0.5 as the discrimination criterion to evaluate the quality of the PLSR model described in step 2;

[0030] Among them, R² (coefficient of determination) represents the explanatory ability of the model for the variation of observed data; Q² (cross-validation coefficient of determination) is used to evaluate the model's ability to predict new data, which is calculated through cross-validation and reflects the generalization ability of the model; Q²cum (cumulative cross-validation coefficient of determination) is the cumulative performance of Q² values for multiple latent variables or components, usually used to evaluate the model's predictive ability under different numbers of latent variables. An increase in Q²cum generally means that the stability and predictive ability of the model have been improved under multiple latent variables; R², Q², and Q²cum together provide a comprehensive evaluation of the performance of the PLSR model. R² focuses on the goodness of fit of the model, while Q² and Q²cum focus on the predictive ability and stability of the model, R 2 , Q 2 and Q 2 cum, the closer the threshold is to 1, the more accurate the constructed partial least squares regression model PLSR is.

[0031] Furthermore, the monthly TWSA data of the area to be studied is processed by the method of averaging the data set, which is derived from the GRACE-Mascon data product. Preferably, it is the observational data obtained by three official agencies, namely the Center for Space Research (CSR) of the University of Texas at Austin, the Jet Propulsion Laboratory (JPL), or the GeoForschungsZentrum Potsdam (GFZ) in Germany, using the GRACE satellite. The data time range in the GRACE-Mascon data product starts from 2002 to the present, and the original spatial resolution is 0.5°.

[0032] Furthermore, the mathematical expression of the GWR model described in step 3-2 is as follows:

[0033] ;

[0034] where y i represents the dependent variable; x ip represents the independent variable, that is, the value of the pth explanatory variable at the ith sample point; (u i , V i ) represents the geographical coordinates of the ith point; β p (u i , V i ) represents the regression coefficient, β0(u i , V i ) represents the intercept of the ith point, and ε(u i , V i ) represents the residual of the ith point.

[0035] Further, in step 3-3, a bilinear spatial interpolation algorithm is used to convert the geographically weighted regression coefficients and regression residuals at low spatial resolution to high spatial resolution.

[0036] The present invention also provides an application of the above method in the assessment of groundwater changes in ecological restoration areas.

[0037] The method for obtaining GWSA by the groundwater downscaling model driven by input and output variables of the homologous hydrological model provided by the present invention can be used to monitor the overall changes in groundwater resources in ecological restoration areas under the dual impacts of climate change and human activities, identify the impacts of ecological restoration work on the groundwater system, timely discover and solve ecological environment problems, evaluate the effectiveness and sustainability of ecological restoration work, guide the optimization adjustment and post-maintenance of ecological restoration work, and provide corresponding technical support for the implementation of subsequent projects and the selection of planted vegetation in regional ecological restoration.

[0038] The beneficial effects of the present invention compared with the prior art are as follows:

[0039] Aiming at the problem that the coarse resolution of GRACE data limits its potential in groundwater assessment applications, the present invention integrates remote sensing observations, hydrological numerical simulation techniques, and a downscaling algorithm of machine learning, and for the first time proposes a research and construction system for a groundwater storage downscaling model driven by input and output variables of a homologous hydrological model, effectively overcoming the problem of non-closure of water balance in grid cells caused by the multi-source nature of input variables in previous downscaling models, and improving the spatial resolution of GRACE groundwater storage change monitoring from the original 0.5° to higher precision. In this way, on the basis of retaining the spatial distribution pattern of GWSA changes at the original resolution, a more refined GWSA change signal is obtained, significantly improving its spatial resolution and reliability, improving the "blocky" characteristics of GWSA changes at the original spatial resolution, and obtaining more continuous GWSA change information in space.

[0040] The method of the present invention is fast, efficient, convenient, and low in economic cost, will significantly improve the accuracy of groundwater change monitoring, provide a new evaluation idea for identifying groundwater level changes, and is particularly important for groundwater monitoring and assessment in ecological restoration areas without groundwater well monitoring. Description of the Drawings

[0041] Figure 1 VIP values of different input and output variables calculated based on the PLSR model relative to TWSA estimated by GRACE for Region E of a certain plain.

[0042] Figure 2It is a comparison chart of the monthly-scale GWSA sequence of a certain plain obtained from groundwater level observations during the period of 2003 - 2016 and the monthly-scale GWSA sequence of a certain plain estimated based on GRACE data combined with different downscaling models. Among them, Figure 2 (a) is the comparison between the monthly-scale GWSA sequence of a certain plain obtained from groundwater level observation data (obtained from GLDAS data) and the EGWSA sequence of a certain plain estimated by a groundwater downscaling model driven by the input and output variables of the GRACE data combined with a homologous hydrological model (ESSI-3 model); Figure 2 (b) is the comparison between the monthly-scale GWSA sequence of a certain plain obtained from groundwater level observation data (obtained from GLDAS data) and the E GWSA sequence of a certain plain estimated by a groundwater downscaling model driven by the input and output variables of the GRACE data combined with a non-homologous hydrological model (VIC model); Figure 2 (c) is the comparison between the monthly-scale GWSA sequence of a certain plain obtained from groundwater level observation data (obtained from GLDAS data) and the EGWSA sequence of a certain plain estimated by a groundwater downscaling model driven by the input and output variables of the GRACE data combined with a non-homologous hydrological model (Noah model); Among them, Figure 2 (b) and Figure 2 (c) the hydrological variables come from the VIC model ( Figure 2 (b)) and the Noah model ( Figure 2 (c)) respectively, but the climate and vegetation information input into the groundwater downscaling model is not the data source driven by the VIC model or the Noah model.

[0043] Figure 3 It is the spatial distribution pattern of the multi-year average value of E GWSA of a certain plain simulated by different model methods during 2003 - 2016. Among them, Figure 3 (a) represents the spatial distribution pattern of the multi-year average value of E GWSA of a certain plain estimated by a groundwater downscaling model driven by the input and output variables of the GRACE data combined with a non-homologous hydrological model (where the hydrological variables come from the VIC model); Figure 3 (b) represents the spatial distribution pattern of the multi-year average value of E GWSA of a certain plain estimated by a groundwater downscaling model driven by the input and output variables of the GRACE data combined with a non-homologous hydrological model (where the hydrological variables come from the Noah model); Figure 3 (c) represents the spatial distribution pattern of the multi-year average value of E GWSA of a certain plain estimated by the GRACE data combined with the input and output variables of a homologous hydrological model (ESSI-3 model) at the original spatial resolution; Figure 3(d) represents the spatial distribution pattern of the multi-year average value of EGWSA in a certain plain obtained by the groundwater downscaling model driven by the input and output variables of the present invention based on GRACE data and a homologous hydrological model.

[0044] Figure 4 To obtain the GWSA data at the original coarse resolution (0.5°) from GRACE products ( Figure 4 (a), 4(c), 4(e)) and the spatial distribution of the downscaled (1 km spatial resolution) GWSA obtained by adopting the technical solution of the present invention ( Figure 4 (b), 4(d), 4(f)). Among them, Figure 4 (a) and Figure 4 (b) respectively represent the spatial variations of GWSA in the ecological restoration area A at 0.5° and 1 km spatial resolutions; Figure 4 (c) and Figure 4 (d) respectively represent the spatial variations of GWSA in the ecological restoration area B at 0.5° and 1 km spatial resolutions; Figure 4 (e) and Figure 4 (f) respectively represent the spatial variations of GWSA in the ecological restoration area C at 0.5° and 1 km spatial resolutions.

[0045] Figure 5 shows the rice and corn planting patterns in a certain plain E in 2016 ( Figure 5 (a)) and the spatial distribution pattern of the change trend of groundwater storage in a certain plain E from 2003 to 2016 ( Figure 5 (b)). Specific Embodiments

[0046] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described and verified in detail below in combination with a specific research area and relevant attached drawings.

[0047] The present invention first combines satellite observation data to screen the input and output variables from a set of hydrological models through the PLSR model to obtain the prediction variables of the homologous groundwater downscaling model; then, constructs a non-stationary relationship between the prediction variables and the target variables at a relatively coarse spatial resolution through the GWR model; subsequently, uses fine-scale prediction factors such as hydrology, climate, topography, and vegetation in combination with the established downscaling model to improve the explanatory ability of the spatial distribution of terrestrial water storage and achieve the segmentation of coarse-pixel terrestrial water storage into fine pixels; finally, in each high-precision grid cell, according to the relationship between different water storages, the contribution of other water storage information is separated from the terrestrial water storage information, so as to separate and obtain the GWSA result of groundwater storage change at high spatial resolution.

[0048] The following takes a certain plain E as an example to introduce the present invention in detail.

[0049] (I) Research area

[0050] Plain E in a certain area is a low-lying plain formed by the alluvial deposits of multiple large rivers. The climate in this area is a temperate humid and semi-humid continental monsoon climate, and the plain area accounts for about 56% of the total area of the region. Due to its regional characteristics of fertile soil, flat terrain, abundant water sources, and suitable climate, after years of development and construction, it has become an important major grain-producing area and a typical high-intensity agricultural development area in China. However, due to the over-reliance of paddy field irrigation in Plain E on groundwater, groundwater resources have become the main recharge source of agricultural water in this area. In some areas, the groundwater resources have been severely over-exploited, resulting in the formation of some relatively large groundwater loss zones in the hinterland of this plain area, and problems such as groundwater level decline, land subsidence, and groundwater funnels have become increasingly prominent.

[0051] Therefore, realizing the assessment of the groundwater level changes in the typical ecological restoration area of Plain E is of great significance for timely discovering and solving regional ecological environment problems, evaluating the effectiveness and sustainability of ecological restoration work, and guiding the optimization adjustment and post-maintenance of ecological restoration work.

[0052] (II) Data

[0053] In this embodiment, the data used is as follows:

[0054] 2.1 GRACE Mascon data product

[0055] The GRACE Mascon data product comes from the observation data obtained by the GRACE satellite by three official institutions, namely the Center for Space Research (CSR) at the University of Texas at Austin, the Jet Propulsion Laboratory (JPL), and the GeoForschungsZentrum Potsdam (GFZ) in Germany. The data time range provided by the GRACE Mascon data product starts from 2002 to the present, and the original spatial resolution is 0.5°.

[0056] The GRACE Mascon data product provides information on the dynamic change of terrestrial water storage (TWSA) in the area to be studied in the present invention, and the uncertainty in the GRACE model post-processing process is reduced by the method of averaging the data set.

[0057] 2.2 CMFD reanalysis data

[0058] Select the CMFD (meteorological data), a dataset of surface meteorological elements in China provided by the National Tibetan Plateau Data Center of China, as the driving data for the distributed hydrological model of the present invention, and also provide the source of meteorological variables for the subsequent downscaling model. The CMFD dataset has a time resolution of 3 hours and a horizontal spatial resolution of 0.1°. The applicability of the dataset in the Chinese region has been widely verified. Seven meteorological elements, namely, mean air temperature, precipitation, air pressure, wind speed, downward longwave radiation, downward shortwave radiation, and relative humidity, can be extracted from the CMFD dataset.

[0059] 2.3 Remote sensing vegetation data products

[0060] The leaf area index LAI mainly uses a global LAI product with a spatial resolution of 8 km based on GLOBMAP LAI (version 3). This long-term LAI data product is composed of AVHRR (Advanced Very High Resolution Radiometer) LAI and MODIS LAI. Among them, the time range of the AVHRR LAI data product is from 1981 to 2000, with a time interval of 16 days, while the time range of the MODIS LAI data product is from 2001 to the present, with a time interval of 8 days.

[0061] 2.4 Distributed hydrological model

[0062] In this embodiment, the distributed land surface hydrological model ESSI-3 is used to obtain the prediction variables of the downscaling model based on the input and output variables of the homologous hydrological model. The ESSI-3 model simulates the hydrological process by simulating the energy and water exchange between the atmosphere and the surface. It can calculate key hydrological variables such as soil moisture, evapotranspiration, and runoff, and provide the land hydrological response on different time scales. The ESSI-3 model requires atmospheric driving data (such as precipitation, air temperature, wind speed, radiation, etc.) and underlying surface driving data as inputs to drive the regional hydrological process simulation. In this embodiment, the atmospheric driving data and underlying surface driving data used are as shown in 2.2 and 2.3.

[0063] 2.5 GLDAS data

[0064] In this embodiment, the GLDAS data is used as the control verification data for the downscaling results of this study.

[0065] The Global Land Data Assimilation System (GLDAS) is a data system jointly developed by NASA and NOAA for global land surface observation and simulation. The GLDAS system integrates multiple models such as Noah and VIC. Through the assimilation of global observation data, it provides rich data support for the study of land surface hydrology and climate change. In this embodiment, the Noah and VIC models are selected for use (where the Noah model pays more attention to the moisture and energy exchange between the atmosphere and land, while the VIC model focuses more on the spatial heterogeneity of surface hydrological processes). The monthly-scale land water storage component products of the Noah and VIC models (including multi-layer soil water content, snow water equivalent, etc., with the original data resolution of 0.25°×0.25°) are extracted from the GLDAS data to be used as the control verification data for the downscaling results of this study.

[0066] (III) Research Framework

[0067] In this embodiment, the high-resolution GWSA results of area E in a certain plain are obtained through the following steps, specifically including the following steps:

[0068] Step 1: Acquisition of input and output variables of the homologous hydrological model

[0069] Based on the distributed land surface hydrological model ESSI-3, combined with multi-source remote sensing data and reanalysis data, a complete hydrological model that can accurately simulate the hydrological situation in area E of a certain plain is constructed; in order to obtain data with a unified resolution, the CMFD and LAI data are resampled to a 1 km spatial resolution through bilinear interpolation and used to drive the hydrological model to obtain hydrological output variables with a 1 km spatial resolution.

[0070] Step 2: Screening of prediction variables for the downscaling model

[0071] Based on the PLSR model, using the monthly TWSA data of area E in a certain plain estimated by GRACE as the target variable, variables with a VIP value (Variable Importance in Projection, VIP) > 0.8 are screened from the driving input variables and model output hydrological variables of the ESSI-3 model as the prediction variables for the downscaling model;

[0072] In this step, the quality of the PLSR model is evaluated with the criteria that the values of R 2 , Q 2 and Q 2 cum are all greater than 0.5.

[0073] For example Figure 1As shown, the VIP values of 11 input and output variables, such as precipitation, average temperature, surface temperature, wind speed, heat flux, soil moisture, snow water equivalent, surface runoff, baseflow, actual evapotranspiration, and normalized difference vegetation index, calculated based on the PLSR model, with respect to the TWSA estimated by GRACE are presented. It can be seen that the VIP values of precipitation, average temperature, soil moisture, snow water equivalent, and actual evapotranspiration are all greater than 0.8, indicating that these five variables are the most important related variables for TWSA. Therefore, these five factors are determined as the main factors related to the TWSA change in Region E of a certain plain and are used as predictor variables and input into the downscaling model.

[0074] Step 3: Construction of a groundwater downscaling model driven by input and output variables of a homologous hydrological model

[0075] Step 3-1: Aggregate the high spatial resolution (1 km) of the predictor variables obtained in Step 2 to a low spatial resolution (0.5°) to be consistent with the spatial resolution of the monthly TWSA data in Region E of a certain plain (the spatial resolution of the monthly TWSA data in Region E of a certain plain is 0.5°);

[0076] Step 3-2: Using the predictor variables obtained in Step 2 as independent variables and the monthly TWSA data in Region E of a certain plain as the dependent variable, a geographically weighted regression (GWR) model is used to construct the regression relationship between them at a spatial resolution of 0.5°. According to the regression relationship, the geographically weighted regression coefficients and regression residuals of each covariate at 0.5° are obtained;

[0077] The mathematical expression of the GWR model is as follows:

[0078] ;

[0079] where y i represents the dependent variable; x ip represents the independent variable, that is, the value of the pth explanatory variable at the ith sample point; (u i , V i ) represents the geographical coordinates of the ith point; β p (u i , V i ) is the regression coefficient, β0(u i , V i ) represents the intercept at the ith point, and ε(u i , V i ) represents the residual at the ith point.

[0080] Step 3-3: According to the "relationship scale invariance" hypothesis, the geographically weighted regression relationship established at the low-resolution scale is still applicable to other spatial resolutions. Therefore, in the spatial interpolation process, the bilinear spatial interpolation algorithm is used to convert the geographically weighted regression coefficients and regression residuals at 0.5° to 1 km. The geographically weighted regression coefficients and predictor variables at 1 km are substituted into the above geographically weighted regression equation to calculate the TWSA at 1 km, and the final downscaling result is obtained.

[0081] Step 4: Obtaining the high-resolution GWSA results

[0082] Vertically, the main contributors to TWSA are Soil Moisture Storage Anomaly (SMSA), Snow Water Storage Anomaly (SWSA), and Groundwater Storage Anomaly (GWSA). However, the vertical change in TWSA cannot be directly separated into a specific TWSA component (e.g., the groundwater storage component).

[0083] The method of the present invention first obtains the SWSA and SMSA information from the model output hydrological variables of the distributed hydrological model described in Step 1; then, SMSA and SWSA are removed from the TWSA results at 1 km spatial resolution to obtain the GWSA results at 1 km spatial resolution.

[0084] After obtaining the GWSA results at 1 km spatial resolution in a certain plain E region, the above data can be further combined with the longitude and latitude information of a specific ecological restoration area in the plain E to obtain the groundwater dynamic change results in the ecological restoration area and carry out evaluation and analysis.

[0085] (IV) Result analysis and verification

[0086] To verify the feasibility and accuracy of the technical solution of the present invention, the downscaling results obtained by the technical solution of the present invention are respectively compared and verified with the results of different multi-source data-driven and non-closed feature groundwater downscaling models. The results are as follows:

[0087] 1. Comparison of the monthly-scale GWSA sequence results in a certain plain E

[0088] Figure 2 Shows a comparison chart of the monthly-scale GWSA sequence of a certain plain E obtained from groundwater level observations and the monthly-scale GWSA sequence of a certain plain E estimated based on GRACE data combined with different downscaling models during the period 2003-2016. Among them, Figure 2(a) Comparison between the monthly-scale GWSA sequence of Plain E obtained from groundwater level observation data (acquired from GLDAS data) and the estimated GWSA sequence of Plain E by a groundwater downscaling model driven by the input and output variables of a homologous hydrological model (ESSI-3 model) combined with GRACE data; Figure 2 (b) Comparison between the monthly-scale GWSA sequence of Plain E obtained from groundwater level observation data (acquired from GLDAS data) and the estimated GWSA sequence of Plain E by a groundwater downscaling model driven by the input and output variables of a non-homologous hydrological model (VIC model) combined with GRACE data; Figure 2 (c) Comparison between the monthly-scale GWSA sequence of Plain E obtained from groundwater level observation data (acquired from GLDAS data) and the estimated GWSA sequence of Plain E by a groundwater downscaling model driven by the input and output variables of a non-homologous hydrological model (Noah model) combined with GRACE data; where Figure 2 (b) and Figure 2 (c) The hydrological variables are from the VIC model (Figure 2(b)) and the Noah model ( Figure 2 (c)), respectively, but the climate and vegetation information input into the groundwater downscaling model are not the data sources driving the VIC model or the Noah model.

[0089] The reliability of the GWSA sequences obtained by combining GRACE data with the above three groundwater downscaling models is verified using the GWSA sequence obtained from groundwater level observation data.

[0090] From Figure 2 It can be seen that in terms of the overall trend, the average long-term trend of the GWSA of Plain E estimated by the groundwater downscaling model driven by the input and output variables of the homologous hydrological model based on GRACE data in the present invention is similar to the trend of the GWSA sequence obtained from groundwater level observation data, both showing a high-low-high trend in the time series ( Figure 2 (a)). However, the GWSA sequence of Plain E estimated by the groundwater downscaling model driven by the input and output variables of the non-homologous hydrological model combined with GRACE data has a large difference in the time series trend from the GWSA observed by groundwater level ( Figure 2 (b), 2(c)).

[0091] 2. Comparison of the spatial distribution pattern results of the multi-year average value of GWSA in Plain E

[0092] Figure 3 Shows the spatial distribution pattern of the multi-year average value of GWSA in Plain E simulated by different model methods from 2003 to 2016. Among them, Figure 3(a) represents the spatial distribution pattern of the multi-year average value of EGWSA in a certain plain estimated by a groundwater downscaling model driven by GRACE data combined with input and output variables of a non-homologous hydrological model (where the hydrological variables are from the VIC model, but the climate and vegetation information input into the groundwater downscaling model are not the data sources driving the VIC model); Figure 3 (b) represents the spatial distribution pattern of the multi-year average value of EGWSA in a certain plain estimated by a groundwater downscaling model driven by GRACE data combined with input and output variables of a non-homologous hydrological model (where the hydrological variables are from the Noah model, but the climate and vegetation information input into the groundwater downscaling model are not the data sources driving the Noah model); Figure 3 (c) represents the spatial distribution pattern of the multi-year average value of EGWSA in a certain plain estimated by GRACE data combined with input and output variables of a homologous hydrological model (ESSI-3 model) at the original spatial resolution (i.e., before downscaling); Figure 3 (d) represents the spatial distribution pattern of the multi-year average value of EGWSA in a certain plain obtained by the groundwater downscaling model driven by GRACE data and input and output variables of a homologous hydrological model according to the present invention.

[0093] It can be seen that the simulation results of the model according to the present invention ( Figure 3 (d)) and the spatial distribution patterns of the multi-year average value of EGWSA in a certain plain estimated by the first three models ( Figure 3 (a), 3(b), 3(c)) have good spatial consistency, and all have the spatial distribution trend that EGWSA is relatively high in the northwest and southeast regions of a certain plain E, while it is relatively low in the central region.

[0094] On the basis that the above results verify the rationality of the method according to the present invention in the trend of the EGWSA sequence at the monthly scale in Result 1, it further verifies the rationality of the simulation results of the model according to the present invention in the spatial distribution.

[0095] Furthermore, from the perspective of pixels, Figure 3 (c) The continuity of the spatial distribution pattern of the multi-year average value of EGWSA in a certain plain estimated by GRACE data combined with input and output variables of a homologous hydrological model (ESSI-3 model) is better than that of Figure 3 (a) and 3(b) The continuity of the spatial distribution pattern of the multi-year average value of EGWSA in a certain plain estimated by GRACE data combined with input and output variables of a non-homologous hydrological model.

[0096] However, comparing the result estimated at the original resolution by GRACE data combined with input and output variables of a homologous hydrological model ( Figure 3 (c)) with the result obtained after downscaling according to the present invention ( Figure 3(d)), it can be seen that on the basis of retaining the spatial distribution pattern of the original resolution results, the finer downscaling results of the method of the present invention obtain a finer GWSA signal. Figure 3 As shown in (c), the results of the spatial distribution pattern of the multi-year average of E GWSA in a certain plain on adjacent pixels are still "blocky", while the results of the downscaling model of the present invention are more continuous in spatial distribution.

[0097] (V) Application example of the method of the present invention in the assessment of groundwater changes in areas where ecological restoration has been carried out

[0098] After obtaining high-resolution GWSA spatio-temporal data, the assessment work of groundwater changes in the ecological restoration area can be carried out.

[0099] Taking August 2010 as an example, three typical ecological restoration areas in a certain plain E (Ecological Restoration Area A, Ecological Restoration Area B, and Ecological Restoration Area C, all of which are important small watershed comprehensive treatment project areas in this region) were selected for research. Figure 4 The GWSA data at the original coarse resolution (0.5°) obtained from the GRACE product and the downscaled (1 km spatial resolution) GWSA spatial distribution obtained by adopting the technical solution of the present invention are shown. Among them, Figure 4 (a) and Figure 4 (b) respectively represent the GWSA spatial changes in Ecological Restoration Area A at 0.5° and 1 km spatial resolutions; Figure 4 (c) and Figure 4 (d) respectively represent the GWSA spatial changes in Ecological Restoration Area B at 0.5° and 1 km spatial resolutions; Figure 4 (e) and Figure 4 (f) respectively represent the GWSA spatial changes in Ecological Restoration Area C at 0.5° and 1 km spatial resolutions.

[0100] It can be seen from the GWSA profiles of different ecological restoration areas that the downscaled results obtained by adopting the technical solution of the present invention greatly and more clearly and finely improve the estimated GWSA spatial changes in each ecological restoration area. Among them, the GWSA spatial changes along the P1-P2 cross-section line are shown to be relatively rough in the original coarse resolution GRACE GWSA data, while more abundant GWSA spatial changes can be obtained in the downscaled high-resolution GWSA data. For example, the GRACE GWSA results only reflect the GWSA changes in Ecological Restoration Area A with two pixels, while the downscaled GWSA data, while maintaining the change trend of the original GWSA data in this area, more abundantly reflects the GWSA distribution changes within the area.

[0101] The downscaled data can reflect the strong spatial heterogeneity of the GWSA changes in the ecological restoration area, and the groundwater depletion rates in different regions vary due to planting status and other human factors. By analyzing the groundwater changes at the regional spatial scale over different time periods, the impacts of different human activity states and ecological restoration measures (such as vegetation restoration, wetland restoration, etc.) on groundwater resources can be evaluated. For example, from Figure 4 (f), it can be seen that smaller GWSA values were observed in the upper and middle reaches of ecological restoration area C, which may be due to the large amount of groundwater extraction required by the high-density human activities in this area. Therefore, the downscaled high-resolution GWSA data is of great significance for monitoring groundwater changes in smaller regions. The data obtained by directly using GRACE data with a spatial resolution range of 1° to 0.25° to monitor regional GWSA changes has limitations in terms of spatial resolution and coverage, and it is difficult to comprehensively grasp the groundwater changes in the ecological restoration area at a finer scale. In contrast, the high-resolution GWSA after downscaling by the technology of the present invention has significant advantages.

[0102] (VI) Application examples of the method of the present invention in identifying and monitoring groundwater changes in the early stage of ecological restoration

[0103] The present invention can not only be used to evaluate the groundwater changes in the areas where ecological restoration has been carried out, but also can provide a scientific basis for determining suitable ecological restoration areas and measures by identifying and monitoring the groundwater changes in the early stage of ecological restoration in advance. Before the implementation of ecological restoration measures, by analyzing the change trend of groundwater storage, the water resource status of the areas to be restored can be judged, so as to optimize the ecological restoration plan and guide the subsequent restoration work.

[0104] Since 2000, the grain output in plain E has been increasing continuously, and great achievements have been made in agricultural production. However, the crop planting structure and pattern have also changed greatly. A large number of drylands have been converted into paddy fields. While the sown area of rice has increased rapidly, the groundwater extraction and irrigation activities in the region have become increasingly intense, and the phenomenon of groundwater imbalance has gradually emerged. The hotspots of groundwater storage depletion in plain E have occurred widely in the west-central and southwestern regions of the whole region.

[0105] Taking City D, one of the areas with the most concentrated and core paddy field expansion and rice cultivation in Region E of a certain plain, as an example, City D is located on one of the only three black soil plains in the world and has a temperate humid and semi-humid climate. The average frost-free period over the years is as long as 135 days. Due to the fertile soil characteristics of this area, it is suitable for the cultivation of various crops. Especially for rice cultivation, it has natural advantages. Since a large amount of groundwater is required for rice irrigation, the expansion of paddy fields and the concentrated cultivation of rice in City D have continuously increased the water demand for farmland irrigation in this area. However, the available groundwater for irrigation in this area is limited. Therefore, serious development and over-exploitation phenomena have occurred in this area. The insufficient water supply and the "water grabbing" in the high water demand areas within the region have led to the groundwater extraction volume in City D being greater than the recharge volume, and the groundwater is in a negative balance state, forming a hot spot area of groundwater funnel.

[0106] Such as Figure 5 shows the rice and corn planting patterns in Plain E in 2016 ( Figure 5 (a)), and the spatial distribution pattern of the change trend of groundwater storage in Plain E from 2003 to 2016 ( Figure 5 (b)). From Figure 5 (a) combined with Figure 5 (b), it can be intuitively observed that the area near City D is a concentrated distribution area of rice cultivation, and this area has become a hot spot area of groundwater funnel. In the eastern region of Plain E, although there are also large-scale plantings of rice and corn, due to more precipitation in the eastern region than in the western region during the period from 2003 to 2016, the change of groundwater storage in the eastern region shows a relatively more stable or positive trend compared with the western region.

[0107] Thus, it can be seen that the east-west difference in the spatial distribution of the dynamic change of groundwater storage in Plain E also reflects the combined effect of climate change and human activities. In addition, this step provides an important reference for ecological restoration planning and decision-making, ensuring that in subsequent ecological restoration work, groundwater management measures are optimized in a targeted manner to improve the effectiveness and sustainability of restoration. For City D, in order to reduce the water use pressure in this area and achieve ecological restoration, it should be considered to appropriately reduce the rice planting area in the area, or even carry out reverse "paddy fieldization" operations, and carry out regional groundwater ecological restoration work to restore the balance of groundwater and promote the sustainable development of regional agriculture.

[0108] Although the present invention has been disclosed above with preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make possible changes and modifications to the technical solutions of the present invention by using the methods and technical contents disclosed above without departing from the spirit and scope of the present invention. Therefore, any simple modifications, equivalent changes and decorations made to the above embodiments based on the technical essence of the present invention without departing from the technical solutions of the present invention all belong to the protection scope of the technical solutions of the present invention.

[0109] The content not detailed in the description of the present invention belongs to the well-known technology of those skilled in the art.

Claims

1. A method for obtaining GWSA based on a groundwater downscaling model driven by input and output variables of a homologous hydrological model, characterized in that: The following steps are involved: Step 1: Obtaining input and output variables of the homologous hydrological model Based on a distributed hydrological model, combined with multi-source remote sensing data and reanalysis data, a hydrological model simulating the hydrological conditions of the area to be studied is constructed; high spatial resolution climate data and vegetation data are used as model-driven input variables to drive the hydrological model, and high spatial resolution model output hydrological variables are obtained; The climate data include precipitation, average temperature, surface temperature, wind speed, snow cover, and heat flux; the vegetation data include normalized vegetation index, leaf area index, soil type, lithology, and hydraulic conductivity; the hydrological variables include soil moisture, snow water equivalent, surface runoff, groundwater runoff, actual evapotranspiration, and vegetation root zone moisture; the high spatial resolution climate data and vegetation data are obtained by resampling the original low spatial resolution climate data and vegetation data using a bilinear interpolation method; Step 2: Downscaling model predictor variable screening Based on the PLSR model, the monthly TWSA data of the area to be studied is used as the target variable, and the variables with VIP values ​​greater than 0.8 are selected from the model-driven input variables and model output hydrological variables described in step 1 as the predictor variables of the downscaling model; Among them, R 2 , Q 2 and Q 2 The values ​​of cum are all greater than 0.5 as the judgment standard to evaluate the quality of the PLSR model; Step 3: Construction of groundwater downscaling model driven by input and output variables of the homologous hydrological model Step 3-1: Sample the predicted variables obtained in step 2 from high spatial resolution to low spatial resolution to keep consistent with the spatial resolution of the monthly TWSA data of the study area; The monthly TWSA data of the area to be studied are processed by a data set averaging method, which is derived from the GRACE-Mascon data product; Step 3-2: Use the predicted variables obtained in step 2 as independent variables and the monthly TWSA data of the study area as dependent variables to establish a GWR model, and obtain the geographically weighted regression coefficients and regression residuals at low spatial resolution based on the regression relationship; Step 3-3: Convert the geographically weighted regression coefficients and regression residuals at low spatial resolution to high spatial resolution, and substitute the geographically weighted regression coefficients and predictor variables at high spatial resolution into the above GWR model to obtain TWSA at high spatial resolution; Step 4: High-resolution GWSA acquisition Eliminate SWSA and SMSA data from high spatial resolution TWSA results to obtain high spatial resolution GWSA results; The SWSA and SMSA are obtained from the model output hydrological variables of the distributed hydrological model in step 1.

2. The method according to claim 1, characterized in that The distributed hydrological model in step 1 is an ESSI series hydrological model, a SWAT hydrological model or a MIKE-SHE hydrological model.

3. The method according to claim 1, characterized in that The mathematical expression of GWR in step 3-2 is as follows: ； Among them, y i represents the dependent variable; x ip represents the independent variable, that is, the value of the pth explanatory variable at the i-th sample point; (u i ,V i ) represents the geographical coordinates of the i-th point; β p (u i , V i ) represents the regression coefficient, β0(u i , V i ) represents the intercept of the i-th point, ε(u i ,V i ) represents the residual of the i-th point.

4. The method according to claim 1, characterized in that In step 3-3, a bilinear spatial interpolation algorithm is used to convert the geographically weighted regression coefficients and regression residuals at low spatial resolution to high spatial resolution.

5. Application of the method described in any one of claims 1 to 4 in the assessment of groundwater changes in ecological restoration areas.

Images