Method for obtaining GWSA based on groundwater downscaling model driven by input and output variables of homologous hydrological model and application of method in ecological restoration

Abstract

The invention belongs to the cross technical field of satellite gravity mechanics, hydrology and the like, and particularly relates to a method for obtaining GWSA based on a groundwater downscaling model driven by input and output variables of a homologous hydrological model and application of the method in ecological restoration. The method comprises the following steps: (1) screening input and output variables from a hydrological model 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 a predictive variable and a target variable under a relatively rough spatial resolution through a GWR model; (3) segmenting the land water reserves of the coarse pixels into fine pixels by combining the fine scale predictor with the established downscaling model; and (4) in each high-precision grid unit, contributions of other water reserve information are separated from the land water reserve information, so that a GWSA result under high spatial resolution is obtained. According to the method, the accuracy and credibility of a model output result are improved.

Description

Technical Field

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

[0002] Against the backdrop of global change, climate warming and other human activities are changing the global water cycle at an unprecedented rate. As the most active component of the hydrological cycle, groundwater storage anomaly (GWSA) is an important component of the global water cycle, and accurate quantification of groundwater level changes in temporal and spatial scales is of great significance to groundwater governance and management.

[0003] At present, the methods of groundwater monitoring and assessment are mainly based on groundwater monitoring station networks, groundwater models or remote sensing observations. However, as a hidden resource, groundwater is dynamic and complex, resulting in the lack of large-scale, long-term, high-spatial-resolution and reliable groundwater information. The ground-observed groundwater level data available in my country are seriously insufficient in time and space: (i) point-scale groundwater monitoring data or survey reports are often not open to the public, and only provincial / municipal-level observation statistics are available; (ii) Affected by the number of groundwater well monitoring points, the spatial distribution of observation points is uneven, and groundwater changes have high spatial heterogeneity. The groundwater level dynamics observed from underground wells can be used as a quantitative indicator of local groundwater changes, but if you want to carry out extensive groundwater monitoring activities, you need to cover the entire region with a dense well network, which is very costly. In addition, due to the complexity of underground aquifer geology and groundwater recharge / discharge processes, groundwater model simulations have obvious limitations and uncertainties.

[0004] In recent years, although some hydrological models have added groundwater modules, since the hydrological models mainly focus on the natural water cycle or the large-scale water cycle of land-atmosphere coupling, the description of the intrinsic mechanisms and processes of the underground aquifers in the modules is relatively simple, and the quality of groundwater simulation does not meet the requirements.

[0005] Since the launch of the satellite, the GRACE (Gravity Recovery and Climate Experiment) satellite and its follow-up GRACE-FO satellite have provided a new research perspective for the quantitative assessment of groundwater with unprecedented accuracy and resolution. However, due to the characteristics of the satellite structure and orbit configuration, GRACE can only be used for relatively large areas (generally >90,000km 2The coarse spatial resolution of GRACE observations limits its application to small and medium regional scales.

[0006] In recent years, some scholars have adopted the method of downscaling the target variables to improve the spatial resolution of GRACE products. Although the selection of predictor variables plays a vital role in determining the quality and performance of the downscaling model, so far, the sources of the input variables of the existing downscaling models are often multi-source, and there is no unified standard for the selection of input variables of the downscaling model. Existing downscaling model research only focuses on the improvement of the downscaling algorithm (such as Chinese patent CN116383590A), and the selection of input variables of the downscaling model and the exploration of the closure of water budget in hydrological modeling are often ignored. The heterogeneity and inconsistency between data from different sources will lead to deviations in the total amount and distribution of water, which will significantly aggravate the problem of non-closure of water balance in grid units and ultimately affect the reliability of model output. Therefore, it is of great significance to implement strict input data quality control and maintain its consistency in the model construction stage to improve the accuracy and credibility of model output results. Summary of the invention

[0007] In view of the defects in the prior art, the purpose of the present invention is to provide a method for obtaining GWSA based on a groundwater downscaling model driven by the input and output variables of a homologous hydrological model, so as to reduce the uncertainty of previous groundwater downscaling models and improve the accuracy and credibility of the 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 that are affected by human activities and where the groundwater levels have changed significantly from their original states. This is of vital importance for guiding the optimization and adjustment of ecological restoration work and subsequent engineering measures, and improving the effects of ecological restoration.

[0009] The present invention adopts the following technical solutions: A method for obtaining GWSA based on a groundwater downscaling model driven by input and output variables of a homologous hydrological model comprises the following steps: 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; Step 2: Downscaling model predictor variable screening Based on the PLSR model (partial least squares regression model), the monthly TWSA (terrestrial water storage change) data of the study area was used as the target variable, and the variables with VIP value (variable importance in projection, VIP) > 0.8 were selected from the model-driven input variables and model output hydrological variables described in step 1 as the predictor variables of the downscaling model.

[0010] 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 resolution to low spatial resolution to keep consistent with the spatial resolution of the monthly TWSA data of the area to be studied; For example, the original spatial resolution of the TWSA data provided by the GRACE data product is 0.5°, and for the sake of uniformity, the spatial resolution of the predicted variables needs to be sampled from high resolution to 0.5°; 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 (geographically weighted regression 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. Substitute the geographically weighted regression coefficients and prediction variables at high spatial resolution into the above GWR model to obtain TWSA with high spatial resolution.

[0011] Step 4: High-resolution GWSA acquisition By removing SMSA and SWSA from the high spatial resolution TWSA results, we can obtain the high spatial resolution GWSA results.

[0012] In the vertical direction, the main contributors to TWSA are soil moisture storage anomaly (SMSA), snow water storage anomaly (SWSA) and groundwater storage anomaly (GWSA). Therefore, by removing SMSA and SWSA from the high spatial resolution TWSA results, we can get the high spatial resolution GWSA results.

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

[0014] Furthermore, the climate data in step 1 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; and the hydrological variables include soil moisture, snow water equivalent, surface runoff, groundwater runoff, actual evapotranspiration, and vegetation root zone moisture.

[0015] Furthermore, the climate data in step 1 comes from the China Regional Ground Meteorological Element Driven Dataset CMFD provided by the China National Qinghai-Tibet Plateau Data Center. The CMFD data set has a time resolution of 3 hours and a horizontal spatial resolution of 0.1°; the vegetation data comes from the global LAI product.

[0016] 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 using a bilinear interpolation method.

[0017] Furthermore, the multi-source remote sensing data in step 1 refers to various types of data obtained from different remote sensing sensors or platforms, such as satellite image data, which can provide information about surface characteristics (such as land use, vegetation cover, etc.), soil moisture, etc. Combining these remote sensing data with the hydrological model can input more accurate and detailed relevant parameters or status information into the model; The reanalysis data are usually large-scale meteorological and climate data sets that have been processed and integrated. These data can provide long-term historical data and spatiotemporal 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.

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

[0019] Furthermore, the monthly TWSA data of the area to be studied are processed by the data set averaging method, which is derived from the GRACE-Mascon data product, preferably the observation data obtained by the three official institutions of the Center for Space Research (CSR), Jet Propulsion Laboratory (JPL) or GeoForschungsZentrumPotsdam (GFZ) in Germany using the GRACE satellite. The data time range in the GRACE-Mascon data product is from 2002 to the present, and the original spatial resolution is 0.5°.

[0020] Furthermore, the mathematical expression of the GWR model 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; (μ i ,υ i ) represents the geographical coordinates of the i-th point; β p (μ i ,υ i ) represents the regression coefficient, β0(μ i ,υ i ) represents the intercept of the i-th point, ε(μ i ,υ i ) represents the residual of the i-th point.

[0021] Furthermore, 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.

[0022] The present invention also provides application of the above method in assessing groundwater changes in ecological restoration areas.

[0023] The method for obtaining GWSA based on the groundwater downscaling model driven by the 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 influence of climate change and human activities, identify the impact of ecological restoration work on the groundwater system, promptly discover and solve ecological and environmental problems, evaluate the effectiveness and sustainability of ecological restoration work, guide the optimization adjustment and subsequent maintenance of ecological restoration work, and provide corresponding technical support for the implementation of subsequent regional ecological restoration projects and the selection of planted vegetation.

[0024] The beneficial effects of the present invention compared with the prior art are: In response to the problem that the coarse resolution of GRACE data limits its potential for groundwater assessment applications, this paper integrates remote sensing observations, hydrological numerical simulation technology, and machine learning downscaling algorithms, and for the first time proposes a research and construction system for groundwater storage downscaling models driven by homologous hydrological model input and output variables, effectively overcoming the problem of grid unit water balance non-closure caused by the multi-source input variables of previous downscaling models, and improving the spatial resolution of GRACE groundwater storage change monitoring from the original 0.5° to a higher accuracy. In this way, while 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 "block" characteristics of GWSA changes at the original spatial resolution, and obtaining more continuous GWSA change information in space.

[0025] The method of the invention is fast, efficient, convenient, and low-cost. It will significantly improve the accuracy of groundwater change monitoring and provide a new evaluation method for identifying groundwater level changes, especially for groundwater monitoring and evaluation in ecological restoration areas without groundwater well monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 VIP values ​​of different input and output variables calculated by PLSR model in a plain area E relative to TWSA estimated by GRACE.

[0027] Figure 2This is a comparison chart of the monthly-scale GWSA series of a certain plain obtained from groundwater level observations during 2003-2016 and the monthly-scale GWSA series of a certain plain estimated based on GRACE data combined with different downscaling models. Among them, Figure (a) is a comparison of 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 the groundwater downscaling model driven by the input and output variables of the GRACE data combined with the homologous hydrological model (ESSI-3 model); Figure (b) is a comparison of 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 the groundwater downscaling model driven by the input and output variables of the GRACE data combined with the non-homologous hydrological model (VIC model); Figure (c) is a comparison of 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 the groundwater downscaling model driven by the input and output variables of the GRACE data combined with the non-homologous hydrological model (Noah model). Comparison of the sequences; the hydrological variables in Figure (b) and Figure (c) come from the VIC model (Figure (c)) and the Noah model (Figure (d)), respectively, but the climate and vegetation information input into the groundwater downscaling model is not a VIC model or Noah model driven data source.

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

[0029] Figure 4The GWSA data at the original coarse resolution (0.5°) obtained from the GRACE product (Figures (a), (c), (e)) and the spatial distribution of the GWSA after downscaling (1 km spatial resolution) obtained using the technical solution of the present invention (Figures (b), (d), (f)). Among them, Figures (a) and (b) represent the spatial changes of GWSA in ecological restoration area A at 0.5° and 1 km spatial resolutions, respectively; Figures (c) and (d) represent the spatial changes of GWSA in ecological restoration area B at 0.5° and 1 km spatial resolutions, respectively; Figures (e) and (f) represent the spatial changes of GWSA in ecological restoration area C at 0.5° and 1 km spatial resolutions, respectively.

[0030] Figure 5 The rice and corn planting pattern of a certain plain E in 2016 (Figure (a)) and the spatial distribution pattern of groundwater storage change trend of a certain plain E from 2003 to 2016 (Figure (b)). DETAILED DESCRIPTION

[0031] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention will be further described and verified in detail below in conjunction with a specific research area and related drawings.

[0032] The present invention firstly combines satellite observation data with the PLSR model to screen the input and output variables from a set of hydrological models to obtain homologous groundwater downscaling model prediction variables; then, the GWR model is used to construct a non-stationary relationship between the prediction variables and the target variables at a coarser spatial resolution; subsequently, the downscaling model established by combining fine-scale prediction factors such as hydrology, climate, topography and vegetation is used to improve the explanatory ability of the spatial distribution of terrestrial water reserves, and to achieve the segmentation of coarse pixel terrestrial water reserves into fine pixels; finally, in each high-precision grid unit, the contribution of other water reserve information is separated from the terrestrial water reserve information based on the relationship between different water reserves, thereby separating the GWSA results of groundwater reserve changes at high spatial resolution.

[0033] The present invention will be described in detail below by taking a certain plain E as an example.

[0034] 1. Study area Plain E is a low-lying plain formed by the alluvial deposits of many large rivers. The climate in this area is 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 resources and suitable climate, after years of development and construction, it has become an important grain production area and a typical high-intensity agricultural development area in my country. However, due to the excessive reliance on groundwater for paddy field irrigation in Plain E, groundwater resources have become the main source of agricultural water supply in the region. Groundwater resources in some areas are seriously over-exploited, resulting in the formation of some large groundwater loss zones in the hinterland of the plain area, and problems such as groundwater level decline, ground subsidence, and groundwater funnels are becoming increasingly prominent.

[0035] Therefore, the evaluation of groundwater level changes in a typical ecological restoration area in a plain region E is of great significance for timely discovering and solving regional ecological and environmental problems, evaluating the effectiveness and sustainability of ecological restoration work, and guiding the optimization, adjustment and subsequent maintenance of ecological restoration work.

[0036] 2. Data In this embodiment, the data used are as follows: 2.1 GRACE Mascon data products GRACE Mascon data products come from observation data obtained by three official institutions: the Center for Space Research (CSR) at the University of Texas, the Jet Propulsion Laboratory (JPL) and the GeoForschungsZentrum Potsdam (GFZ) in Germany using GRACE satellites. The data time range provided by GRACE Mascon data products is from 2002 to the present, with an original spatial resolution of 0.5°.

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

[0038] 2.2 CMFD reanalysis data The China regional ground meteorological element driven dataset CMFD (meteorological data) provided by the China National Tibetan Plateau Data Center is selected as the driving data of the distributed hydrological model of the present invention, and also provides a 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 China has been widely verified. Seven meteorological elements can be extracted from the CMFD dataset, including average temperature, precipitation, air pressure, wind speed, downward longwave radiation, downward shortwave radiation, and relative humidity.

[0039] 2.3 Remote sensing vegetation data products Leaf area index LAI mainly adopts the 8km spatial resolution global LAI product 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 AVHRR LAI data product has a time range of 1981-2000 and a time interval of 16 days, while the MODIS LAI data product has a time range of 2001 to the present and a time interval of 8 days.

[0040] 2.4 Distributed hydrological model 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, runoff, and provide information about land hydrological responses on different time scales. The ESSI-3 model requires atmospheric driving data (such as precipitation, temperature, wind speed, radiation, etc.) and underlying surface driving data as input to drive the simulation of regional hydrological processes. In this embodiment, the atmospheric driving data and underlying surface driving data used are shown in 2.2 and 2.3.

[0041] 2.5 GLDAS data In this example, the GLDAS data are used as the control and verification data for the downscaling results of this study.

[0042] 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 (the Noah model pays more attention to the exchange of water and energy between the atmosphere and land, while the VIC model focuses more on the spatial heterogeneity of surface hydrological processes). The monthly land water storage component products of the Noah and VIC models (including multi-layer soil water content, snow water equivalent, etc., with an original data resolution of 0.25°×0.25°) are extracted from the GLDAS data as control and verification data for the downscaling results of this study.

[0043] 3. Research framework This embodiment first obtains high-resolution GWSA results in a plain area E through the following steps, which specifically include the following steps: Step 1: Obtaining input and output variables of the homologous hydrological model 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 conditions in area E of a plain is constructed; in order to obtain data with uniform resolution, the CMFD and LAI data are resampled to a spatial resolution of 1 km through the bilinear interpolation method, and the hydrological model is driven to obtain the hydrological output variables with a spatial resolution of 1 km.

[0044] Step 2: Downscaling model predictor variable screening Based on the PLSR model, the monthly TWSA data of a plain region E estimated by GRACE were used as the target variable, and the variables with VIP values ​​(Variable Importance in Projection, VIP) > 0.8 were selected from the driving input variables and model output hydrological variables of the ESSI-3 model as the predictor variables of the downscaling model. In this step, R 2 , Q 2 and Q 2 The cum values ​​are all greater than 0.5, which is used as the judgment standard to evaluate the quality of the PLSR model.

[0045] like Figure 1As shown in the figure, 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, groundwater runoff, actual evapotranspiration, and normalized difference vegetation index calculated based on the PLSR model are relative to the TWSA estimated by GRACE. 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 of TWSA. Therefore, these five factors are identified as the main factors related to the change of TWSA in a plain E area and are input into the downscaling model as predictor variables.

[0046] Step 3: Construction of groundwater downscaling model driven by input and output variables of the homologous hydrological model Step 3-1: Aggregate the high spatial resolution (1 km) of the prediction variables obtained in step 2 to a low spatial resolution (0.5°) to be consistent with the spatial resolution of the monthly TWSA data of a certain plain E region (the spatial resolution of the monthly TWSA data of a certain plain E region is 0.5°); Step 3-2: Take the predicted variable obtained in step 2 as the independent variable and the monthly TWSA data of a plain region E as the dependent variable. Use the geographically weighted regression GWR model to build a regression relationship between the two at a spatial resolution of 0.5°. According to the regression relationship, obtain the geographically weighted regression coefficient and regression residual of each covariate at 0.5°. The mathematical expression of the GWR model 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; (μ i ,υ i ) represents the geographical coordinates of the i-th point; β p (μ i ,υ i ) is the regression coefficient, β0(μ i ,υ i ) represents the intercept of the i-th point, ε(μ i ,υ i ) represents the residual of the i-th point.

[0047] Step 3-3: According to the "relationship scale invariance" hypothesis, the geographically weighted regression relationship established at a low-resolution scale is still applicable to other spatial resolutions. Therefore, in the spatial interpolation processing part, a bilinear spatial interpolation algorithm is used to convert the 0.5° geographically weighted regression coefficients and regression residuals to 1 km. Substitute the 1km geographically weighted regression coefficients and prediction variables into the above geographically weighted regression equation to calculate the 1km TWSA and obtain the final downscaling result.

[0048] Step 4: Obtain high-resolution GWSA results In the vertical direction, the main contributors to TWSA are soil moisture storage anomaly (SMSA), snow water storage anomaly (SWSA), and groundwater storage anomaly (GWSA). However, the vertical TWSA variation cannot be directly separated into a specific TWSA component (e.g., groundwater storage component).

[0049] The method of the present invention first obtains 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 result with a spatial resolution of 1 km to obtain the GWSA result with a spatial resolution of 1 km.

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

[0051] (IV) Result analysis and verification In order 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 compared and verified with the results of groundwater downscaling models driven by different multi-source data and non-closed characteristics. The results are as follows: 1. Comparison of GWSA series results on a certain plain at the E-month scale Figure 2A comparison chart shows the monthly GWSA series of a certain plain obtained from groundwater level observations during the period 2003-2016 and the monthly GWSA series of a certain plain estimated based on GRACE data combined with different downscaling models. Among them, Figure (a) is a comparison of the monthly scale GWSA sequence of a certain plain E obtained from groundwater level observation data (obtained by GLDAS data) and the E GWSA sequence of a certain plain estimated by the groundwater downscaling model driven by the input and output variables of the GRACE data combined with the homologous hydrological model (ESSI-3 model); Figure (b) is a comparison of the monthly scale GWSA sequence of a certain plain E obtained from groundwater level observation data (obtained by GLDAS data) and the E GWSA sequence of a certain plain estimated by the groundwater downscaling model driven by the input and output variables of the GRACE data combined with the non-homologous hydrological model (VIC model); Figure (c) is a comparison of the monthly scale GWSA sequence of a certain plain E obtained from groundwater level observation data (obtained by GLDAS data) and the E GWSA sequence of a certain plain estimated by the groundwater downscaling model driven by the input and output variables of the GRACE data combined with the non-homologous hydrological model (Noah model). Comparison of the sequences; the hydrological variables in Figure (b) and Figure (c) come from the VIC model (Figure (c)) and the Noah model (Figure (d)), respectively, but the climate and vegetation information input into the groundwater downscaling model is not a VIC model or Noah model driven data source.

[0052] The GWSA sequence obtained from groundwater level observation data is used to verify the reliability of the GWSA sequence obtained from GRACE data combined with the above three groundwater downscaling models.

[0053] Depend on Figure 2 It can be seen that in terms of the overall trend, the average long-term trend of the E GWSA of a certain plain estimated by the groundwater downscaling model driven by GRACE data and the input and output variables of the homologous hydrological model is similar to the trend of the GWSA series obtained from the groundwater level observation data, both showing a high-low-high trend in the time series (Figure (a)). However, the E GWSA series of a certain plain estimated by the groundwater downscaling model driven by the input and output variables of the non-homologous hydrological model based on GRACE data is quite different from the GWSA observed by the groundwater level in terms of time series trend (Figures (b) and (c)).

[0054] 2. Comparison of the spatial distribution pattern of the multi-year average values ​​of E GWSA in a certain plain Figure 3The multi-year average spatial distribution pattern of E GWSA in a certain plain simulated by different model methods from 2003 to 2016 is shown. Among them, Figure (a) represents the multi-year average spatial distribution pattern of E GWSA of a certain plain estimated based on GRACE data combined with a groundwater downscaling model driven by non-homologous hydrological model input and output variables (wherein the hydrological variables come from the VIC model, but the climate and vegetation information input to the groundwater downscaling model is not the VIC model driving data source); Figure (b) represents the multi-year average spatial distribution pattern of E GWSA of a certain plain estimated based on GRACE data combined with a groundwater downscaling model driven by non-homologous hydrological model input and output variables (wherein the hydrological variables come from the Noah model, but the climate and vegetation information input to the groundwater downscaling model is not the Noah model driving data source); Figure (c) represents the multi-year average spatial distribution pattern of E GWSA of a certain plain estimated based on GRACE data combined with homologous hydrological model input and output variables (ESSI-3 model) at the original spatial resolution (that is, before downscaling); Figure (d) represents the multi-year average spatial distribution pattern of E GWSA of a certain plain obtained by the groundwater downscaling model driven by GRACE data and homologous hydrological model input and output variables of the present invention.

[0055] It can be seen that the spatial consistency between the simulation results of the model of the present invention (Figure (d)) and the multi-year average spatial distribution pattern of GWSA of a certain plain E estimated by the first three models (Figures (a), (b), and (c)) is good, and all have a spatial distribution trend that GWSA is at a relatively high level in the northwest and southeast regions of a certain plain E, and is at a relatively low level in the central region.

[0056] The above results, based on the verification in Result 1 that the method of the present invention is reasonable in the trend of the monthly-scale GWSA series, further verify the rationality of the simulation results of the model of the present invention in terms of spatial distribution.

[0057] Furthermore, from the pixel point of view, the continuity of the multi-year average spatial distribution pattern of E GWSA in a certain plain estimated based on GRACE data combined with the input and output variables of the homologous hydrological model (ESSI-3 model) in Figure (c) is better than that in Figures (a) and (b) which are estimated based on GRACE data combined with the input and output variables of the non-homologous hydrological model.

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

[0059] (V) Application examples of the method of the present invention in groundwater change assessment in areas where ecological restoration has been carried out After obtaining high-resolution GWSA spatiotemporal data, it is possible to carry out groundwater change assessment work in ecological restoration areas.

[0060] 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 management project areas in the region) were selected for research. Figure 4 The GWSA data at the original coarse resolution (0.5°) obtained from the GRACE product and the spatial distribution of the GWSA after downscaling (1 km spatial resolution) obtained using the technical solution of the present invention are shown in Figures (a) and (b) respectively represent the spatial changes of GWSA in ecological restoration area A at 0.5° and 1 km spatial resolutions; Figures (c) and (d) respectively represent the spatial changes of GWSA in ecological restoration area B at 0.5° and 1 km spatial resolutions; Figures (e) and (f) respectively represent the spatial changes of GWSA in ecological restoration area C at 0.5° and 1 km spatial resolutions.

[0061] From the GWSA profiles of different ecological restoration areas, it can be seen that the downscaling results obtained by using the technical solution of the present invention have improved the estimated GWSA spatial changes in each ecological restoration area to a large extent more clearly and finely. Among them, the GWSA spatial changes along the P1-P2 cross-section are relatively rough in the original coarse-resolution GRACE GWSA data, while the downscaled high-resolution GWSA data can obtain richer GWSA spatial changes. 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 the area, more richly reflects the GWSA distribution changes in the region.

[0062] The downscaled data can reflect the strong spatial heterogeneity of GWSA changes in ecological restoration areas. The groundwater depletion rates in different regions vary due to planting status and other human factors. By analyzing the changes in groundwater at regional spatial scales in different time periods, the impact of different human activities and ecological restoration measures (such as vegetation restoration, wetland restoration, etc.) on groundwater resources can be evaluated. For example, as shown in Figure (f), the GWSA observed in the upstream and middle reaches of ecological restoration area C is relatively small, which may be due to the high-density human activities in the area that require a large amount of groundwater extraction. Therefore, the downscaled high-resolution GWSA data is of great significance for monitoring groundwater changes in smaller areas. The data obtained by directly using GRACE data with a spatial resolution range of 1° to 0.25° to monitor regional GWSA changes have limitations in spatial resolution and coverage, and it is difficult to fully grasp the groundwater changes in ecological restoration areas in a finer range. In contrast, the high-resolution GWSA downscaled by the technology of the present invention has significant advantages.

[0063] (VI) Application examples of the method of the present invention in identifying and monitoring groundwater changes in the early stages of ecological restoration The present invention can not only be used to assess groundwater changes in areas where ecological restoration has been carried out, but can also provide a scientific basis for determining suitable ecological restoration areas and measures by identifying and monitoring groundwater changes in the early stages of ecological restoration. Before the implementation of ecological restoration measures, by analyzing the changing trend of groundwater reserves, the water resource status of the area where ecological restoration is to be carried out can be determined, thereby optimizing the ecological restoration plan and guiding the subsequent restoration work.

[0064] Since 2000, the grain output of a certain plain E has been increasing, and agricultural production has achieved great results. However, the structure and pattern of crop planting have also undergone major changes. A large amount of dry land has been transformed into paddy fields. While the rice planting area has increased rapidly, the regional groundwater extraction and irrigation activities have become more intense, and the phenomenon of groundwater imbalance has gradually appeared. The hot spots of groundwater depletion in a certain plain E are widely occurring in the central and western regions and the southwest regions of the entire region.

[0065] Take D City, one of the most concentrated and core areas of paddy field expansion and rice planting in a plain E area, as an example. D City is located on one of the only three black soil plains in the world. It has a temperate humid semi-humid climate and an average frost-free period of up to 135 days per year. Due to its fertile soil, the area is suitable for the planting of a variety of crops, especially rice planting, which has natural advantages. Since rice irrigation requires a large amount of groundwater, the expansion of paddy fields and concentrated rice planting in D City have increased the water demand for farmland irrigation in the area, while the groundwater available for irrigation in the area is limited. Therefore, the area has experienced serious development and over-exploitation. The insufficient water supply and the "water grabbing" of high water demand areas in the region have led to the groundwater extraction in D City being greater than the recharge, and the groundwater is in a negative equilibrium state, forming a groundwater funnel hotspot area.

[0066] like Figure 5 The rice and corn planting pattern of a plain E in 2016 (Figure (a)) and the spatial distribution pattern of the change trend of groundwater reserves in a plain E from 2003 to 2016 (Figure (b)) are shown. It can be intuitively observed from Figure (a) combined with Figure (b) that the area near D City is a concentrated rice planting area, which has become a hot spot for groundwater funnels. In the eastern part of a plain E, although there are large-scale rice and corn planting, the eastern region has more precipitation than the western region during 2003-2016, so the change of groundwater reserves in the eastern region shows a more stable or positive trend than that in the western region.

[0067] It can be seen that the east-west difference in the spatial distribution of the dynamic changes in groundwater reserves in a certain plain E also reflects the result of the combined impact of climate change and human activities. In addition, this step provides an important reference for ecological restoration planning and decision-making, ensuring that in the 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 the area and achieve ecological restoration, the local government should consider appropriately reducing the rice planting area in the area, or even reverse the "paddy field" operation, and carry out regional groundwater ecological restoration work to restore the balance of groundwater and promote the sustainable development of regional agriculture.

[0068] Although the present invention has been disclosed as above in the form of a preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art may make possible changes and modifications to the technical solution 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 modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall fall within the protection scope of the technical solution of the present invention.

[0069] The contents not described in detail in the specification of the present invention belong to the common knowledge of the professionals in this field.

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; 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; 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; 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 climate data in step 1 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.

4. The method according to claim 1, characterized in that: 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 using a bilinear interpolation method.

5. The method according to claim 1, characterized in that R 2 , Q 2 and Q 2 The cum values ​​are all greater than 0.5, which is used as the judgment standard to evaluate the quality of the PLSR model described in step 2.

6. The method according to claim 1, characterized in that The monthly TWSA data of the area to be studied are processed by the data set averaging method, which is derived from the GRACE-Mascon data product.

7. 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; (μ i ,υ i ) represents the geographical coordinates of the i-th point; β p (μ i ,υ i ) represents the regression coefficient, β0(μ i ,υ i ) represents the intercept of the i-th point, ε(μ i ,υ i ) represents the residual of the i-th point.

8. 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.

9. Application of the method according to any one of claims 1 to 8 in the assessment of groundwater changes in ecological restoration areas.

Images