A method for improving groundwater storage accuracy based on data feature fusion inversion

Abstract

The present invention relates to a method for improving the accuracy of groundwater reserves based on data feature fusion inversion. By fusing GNSS and climate experiment satellite data, the present invention achieves the purpose of improving the accuracy of groundwater reserve anomaly inversion based on GNSS. Based on the traditional GNSS inversion of TWSA, the GRACE satellite and Global Land Assimilation System model results are integrated to construct a new data feature fusion inversion method that helps improve the accuracy of groundwater reserve anomaly inversion. By calculating and analyzing the GWSA of the North China Plain in China for various years, the results show that the South-to-North Water Diversion Project has played a certain role in alleviating GWSA in the North China Plain.

Description

Technical Field

[0001] The present invention belongs to the technical field of hydrology, and specifically is a method for improving the accuracy of groundwater reserves based on data feature fusion inversion. Background Art

[0002] Groundwater is a vital resource for urban construction, agricultural production, and residents' livelihoods, playing a crucial role in my country's national economic development. With the accelerating pace of urbanization and rising industrialization, groundwater overexploitation has become increasingly serious in recent years. In the North China Plain, in particular, this massive groundwater extraction has triggered a series of natural disasters, such as drought, land subsidence, and soil erosion, significantly impacting the stable economic development of the plain. Therefore, accurate monitoring of groundwater reserves is crucial for macro-analysis of the spatial distribution of groundwater resources.

[0003] Traditional groundwater observation methods primarily monitor water levels, water quality, soil moisture, and other parameters through observation wells. However, the cumbersome process of establishing these wells consumes significant manpower and material resources. Furthermore, the spatial distribution of observation wells is highly uneven, making macroscopic monitoring of groundwater reserves over large areas impossible. Consequently, they no longer provide sufficient data support for current scientific regulation of groundwater resources.

[0004] The successful launch of the GRACE gravity satellites in 2002 provided unprecedented observational capabilities for monitoring groundwater storage anomalies. Numerous studies have been conducted in representative regions worldwide, including Cambodia, California, northern China, the North China Plain, and the Guanzhong region. By inverting local gravity anomalies, the GRACE gravity satellites can accurately infer regional load variations. However, due to the aging of the GRACE satellite power supply system, the mission ended in 2017. Its next-generation gravity satellite mission, GRACE Follow-On (GRACE-FO), launched in 2018, experienced a nearly one-year time gap. Furthermore, due to the inherent characteristics of the GRACE satellites, the spatial and temporal resolution of their inversion results is relatively coarse, typically with a spatial resolution of 1°×1° and a temporal resolution of one month. This inherent characteristic makes it difficult to use the GRACE gravity satellites for real-time monitoring of local groundwater storage anomalies. Therefore, it is crucial to find an alternative method for effectively monitoring groundwater storage anomalies with high spatial and temporal resolution. Summary of the Invention

[0005] In order to address the deficiencies in the prior art, the present invention integrates GNSS and climate experiment satellite data to achieve the purpose of improving the accuracy of GNSS-based inversion of groundwater storage anomalies.

[0006] The present invention includes three modules, specifically, a GNSS-based crustal vertical sequence seasonal term inversion module: GNSS data is solved and preprocessed, and the hydrological load deformation sequence is obtained by deducting atmospheric load and non-ocean tidal load. The seasonal variation characteristic sequence of regional terrestrial water storage anomalies is inverted by combining Green's load function and crustal load model;

[0007] Extract TWSA trend term module: Utilizes the equivalent water height variables provided by the GRACE Mascon dataset to extract abnormal changes in regional terrestrial water storage, and obtains the trend term of terrestrial water storage through a novel correlation variational mode decomposition algorithm. This trend term is then added to the GNSS inversion results to obtain the regional TWSA with trend.

[0008] GWSA inversion module based on TWSA correction: By correcting the trend TWSA, which includes canopy water, snow water, and soil water, the abnormal changes in regional groundwater storage are obtained.

[0009] The specific method of the present invention comprises:

[0010] S1, using Green’s function to establish the relationship between load mass and deformation;

[0011] S2, using the equivalent water height variables provided by the GRACE Mascon dataset to extract the abnormal changes in regional terrestrial water storage, and using a novel correlation variational mode decomposition algorithm to obtain the trend term of terrestrial water storage, this trend term is added to the GNSS inversion results to obtain the regional TWSA with trend;

[0012] S3, the seasonal term and trend term of the TWSA obtained by inversion are fused to obtain the TWSA with trend term;

[0013] S4, the regional GWSA is obtained by correcting the TWSA of the trend term, including canopy water, snow water, and soil water.

[0014] Preferably, in step S1, the formula of Green's function for vertical load deformation of the crust is as follows:

[0015]

[0016] Where θ represents the angular radius from the center of the disk; P n represents Legendre polynomials; G represents Newton's gravitational constant; R represents the radius of the Earth; h n represents the load Love number; g represents the acceleration due to gravity.

[0017] The Γ n The derivation of the function is as follows:

[0018]

[0019] Where θ represents the angular radius from the center of the disk; P represents the Legendre polynomial. When n = 0, the expression of the Γ function is as follows:

[0020]

[0021] The obtained solution is regularized using the curvature smoothing algorithm and attached to the solution matrix as a set of constraints.

[0022] For each period of this study, a restrained least squares problem is minimized to estimate the daily terrestrial water storage changes:

[0023] TWSA Season =((Gx-b) / σ) 2 +β 2 (L(x)) 2 →min

[0024] Among them, TWSA Season represents the seasonal characteristic sequence of TWSA obtained by inversion; G represents the Green’s function coefficient matrix; σ represents the standard deviation of the GNSS vertical displacement sequence; b represents the deformation observation sequence of the grid, where b represents the corrected GNSS vertical deformation sequence; L represents the Laplace operator; and β represents the smoothing factor.

[0025] Preferably, in step S2,

[0026]

[0027] Among them, TWSA Trend represents the trend characteristic sequence of TWSA extracted based on GRACE; represents the GRACE original TWSA sequence variable; represents the seasonal characteristic sequence of the TWSA sequence based on GRACE.

[0028] Preferably, in step S3,

[0029] TWSA all =TWSA Season +TWSA Trend

[0030] Among them, TWSA all represents the sum of TWSA series (including trend term and seasonal term); TWSA Season represents the TWSA seasonal term obtained based on GNSS inversion; TWSA Trend Represents the sequence trend item extracted based on GRACE.

[0031] Preferably, in step S4,

[0032] GWSA DFFIM =TWSA all -W Can -W soil -W snow

[0033] Among them, GWSA DFFIM represents the groundwater storage anomaly obtained by DFFIM inversion; TWSA represents the terrestrial water storage anomaly; W can Indicates the change of canopy water content; W soil Indicates the change in total soil water content from 0 to 200 mm.

[0034] The DFFIM proposed in this paper provides an effective technical method for calculating regional GWSA. The research results have important scientific significance and application value for analyzing the impact of the South-to-North Water Diversion Project on groundwater reserves in the North China Plain. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] In order to more clearly illustrate the technical solutions of the implementation methods of the present application, the following is a brief introduction to the drawings required for use in the implementation methods. It should be understood that the following drawings only show certain embodiments of the present application and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without paying any creative work.

[0036] Figure 1 It is a flow chart of the data feature fusion inversion method of the present invention;

[0037] Figure 2 It is a simplified flow chart of the specific method of the present invention;

[0038] Figure 3 These are the temporal and spatial distribution characteristics of TWSA retrieved from GNSS, including (a) the time series diagram of TWSA in the North China Plain retrieved from GNSS, (b) the temporal and spatial distribution diagram of the annual amplitude of TWSA retrieved from GNSS, and (c) the temporal and spatial distribution diagram of the semi-annual amplitude of TWSA retrieved from GNSS;

[0039] Figure 4 Figure 1 is a TWSA inversion diagram of the fusion of GNSS and GRACE, where (a) represents the sequence diagram of the inverted TWSA results, (b) represents the anniversary item of TWSA in the GLDAS dataset, (c) represents the trend item of TWSA in the GLDAS dataset, (d) represents the trend item of the TWSA inversion results of the fusion of GNSS and GRACE, and (e) represents the scatter plot of the experimental results of the present invention and GLDAS;

[0040] Figure 5GWSA inversion effect diagram of the present invention, including (a) GWSA daily slice diagram, (b) GWSA annual amplitude diagram, and (c) GWSA average sequence diagram;

[0041] Figure 6 : This is a comparison diagram of the GWSA inversion effect of the present invention, including (a) GLDAS annual amplitude, (b) GLDAS trend diagram, (c) trend diagram of the inversion result of the algorithm of the present invention, (d) sequence comparison effect diagram, and (e) superimposed power spectrum diagram;

[0042] Figure 7 This is the groundwater storage DSI drought index analysis map, where (a)-(k) DSI spatial distribution map of the North China Plain from 2011 to 2021, (l) DSI series change map;

[0043] Figure 8 This is the spatial distribution and time series map of GWSA in the North China Plain of China from 2011 to 2021. DETAILED DESCRIPTION

[0044] The technical solutions in the embodiments of the present application will be described below in conjunction with the drawings in the embodiments of the present application.

[0045] To make the purpose, technical solutions, and advantages of the embodiments of this application more clear, the technical solutions in the embodiments of this application will be clearly and completely described below in conjunction with the drawings in the embodiments of this application. Obviously, the described embodiments are part of the embodiments of this application, not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0046] like Figure 1 As shown, it includes three groups of modules, namely:

[0047] GNSS-based crustal vertical series seasonal term inversion module: This module solves and preprocesses GNSS data, obtains the hydrological load deformation series by deducting atmospheric loads and non-ocean tidal loads, and combines the Green's load function with the crustal load model to invert the seasonal variation characteristic series of regional terrestrial water storage anomalies;

[0048] Extract TWSA trend term module: The module uses the equivalent water height variables provided by the GRACE Mascon dataset to extract abnormal changes in regional terrestrial water storage, and obtains the trend term of terrestrial water storage through a new correlation variational mode decomposition algorithm. This trend term is added to the GNSS inversion result to obtain the regional TWSA with trend.

[0049] GWSA inversion module based on TWSA correction: By correcting the trend TWSA, which includes canopy water, snow water, and soil water, the abnormal changes in regional groundwater storage are obtained.

[0050] like Figure 2 As shown, the specific method of the present invention includes:

[0051] S1, using Green’s function to establish the relationship between load mass and deformation;

[0052] S2, using the equivalent water height variables provided by the GRACE Mascon dataset to extract the abnormal changes in regional terrestrial water storage, and using a novel correlation variational mode decomposition algorithm to obtain the trend term of terrestrial water storage, this trend term is added to the GNSS inversion results to obtain the regional TWSA with trend;

[0053] S3, the seasonal term and trend term of the TWSA obtained by inversion are fused to obtain the TWSA with trend term;

[0054] S4, the regional GWSA is obtained by correcting the TWSA of the trend term, including canopy water, snow water, and soil water.

[0055] like Figure 3 As shown in Figure 2, the TWSA obtained based on GNSS inversion has obvious annual and semi-annual characteristics. Since the South-to-North Water Diversion Project officially started water transfer in December 2014, the entire study period is divided into 2011-2014 and 2015-2022. Figure 3 (a) It can be seen that the TWSA amplitude after 2015 is significantly higher than that from 2011 to 2014, which is consistent with the actual time node of the South-to-North Water Diversion Project.

[0056] At the same time, the annual characteristics of TWSA in the North China Plain are more obvious than the semi-annual characteristics. The maximum value of the annual amplitude is 170 mm, while the semi-annual amplitude is only 27 mm. From southwest to northeast, the annual amplitude shows a clear downward trend, and its raised areas are located in the Beijing and Tianjin regions. For the semi-annual amplitude, the raised effect in the Beijing and Tianjin regions is more obvious. However, due to the limitations of GNSS inversion of TWSA, only the periodicity of the TWSA sequence can be detected, and the trend term is also particularly important for the effective management of water resources. Therefore, the present invention fills the TWSA obtained by GNSS inversion based on the GRACE Mascon dataset to achieve the purpose of more accurately monitoring regional TWSA.

[0057] like Figure 4As shown, the seasonal term and trend term of TWSA can be effectively inverted based on the inversion idea of the present invention. At the same time, the annual amplitude of TWSA in the GLDAS data set shows a decreasing trend from west to east, and there is a bulge area in the Beijing area. The trend term of the TWSA inversion result shows an increasing trend from north to south, but due to the coarse spatial resolution of the GRACE data set, the spatial resolution of the TWSA inverted by the present invention is relatively coarse. However, the sequence trend consistency is good. The TWSA effect obtained by inversion based on this idea is more consistent than that of the GLDAS data set. Its PCC value is 0.72, and the RMSE between the two sequences is 2.8cm. The idea of the present invention can be used to accurately invert the TWSA in the North China Plain, providing a large amount of data basis for the calculation of GWSA.

[0058] like Figure 5 As shown, the inversion strategy based on the present invention can calculate the GWSA sequence characteristics and annual variation characteristics in the North China Plain. Figure 5 (b) shows the spatial characteristics of the annual amplitude of GWSA, which is based on Figure 5 (a) The vertical data of each grid are calculated. At the same time, the annual amplitude has obvious convex areas in the western part of the North China Plain (Shijiazhuang) and the Tianjin-Beijing area, which is consistent with the previous calculation results. Figure 5 As can be seen from (c), the groundwater reserves in the North China Plain show an overall downward trend. In addition, the annual amplitude after the South-to-North Water Diversion Project (grey shaded area) is significantly higher than that before the South-to-North Water Diversion Project (white area). In order to verify the reliability of the inversion results of this invention, the experimental results are compared with the GWSA variables in the GLDAS V2.2 dataset. The comparison results are shown in Figure 2. Figure 6 As shown:

[0059] like Figure 6 As shown in Figure 2, the abnormal comparison results of groundwater reserves in the North China Plain are calculated based on the inversion method of the present invention. Figure 6 As can be seen in (d), the overall sequence performance compares favorably with GLDAS, with PCC, RMSE, and NSE values of 0.69, 4.01 cm, and 0.61, respectively. Furthermore, the slope of the sequence calculated by our method (-8.52 mm / y) is closer to that of previous studies. Figure 6 (a) and (b) show the spatial distribution of the annual amplitude and trend of GWSA in the GLDAS dataset, respectively. Figure 5 (b) and Figure 6(a) It can be seen that the raised areas of annual amplitude are all located in the Tianjin-Beijing area, Shijiazhuang and the southern part of the North China Plain. However, the annual amplitude of the results obtained based on GNSS and GRACE inversion is slightly larger than that of the GLDAS dataset. The reason is that GNSS observes the vertical deformation of the entire crust, and the present invention only deducts the non-ocean tide and non-ocean atmosphere correction parts. There are also some other disturbances in the GNSS vertical deformation sequence, such as the tectonic deformation of the crust, which makes the GNSS inversion results larger. By comparison Figure 6 (b) and (c) show that, in general, the funnel areas of the GLDAS and the present invention's inversion results are both located in the central part of the North China Plain, while the convex areas of the trend items are both located in the north and south parts of the North China Plain. At the same time, the present invention draws the superimposed power spectrum of the GLDAS and the GWSA inversion results in the North China Plain, which are shown in Figure 2. Figure 6 (e) It can be seen that the two sequences are consistent in the low-frequency and medium-frequency parts. However, due to the presence of a large amount of noise in the GNSS sequence, the GWSA based on the inversion strategy of the present invention is slightly larger than the result of GLDAS in the high-frequency part, further verifying the reliability of the regional GWSA calculated by the inversion strategy of the present invention.

[0060] like Figure 7 As shown in the data, 2012-2013 and 2017-2019 showed obvious mild drought characteristics, among which the drought in 2019 was the most obvious, with a DSI value of -0.12, which was slightly close to moderate drought. At the same time, the DSI value in 2017 reached -0.81, which reached the level of mild drought. Spatially, the southern part of the North China Plain is drier than the northern part ( Figure 7 (a), (b), (e)).

[0061] In 2016, the southern part of the North China Plain was affected by severe convective rainstorms and flooding, resulting in a slight increase in the groundwater drought index. However, localized drought conditions occurred in the Beijing-Tianjin region. From 2020 to 2021, the groundwater drought in the North China Plain eased compared to 2019, with the DSI showing a positive growth trend.

[0062] In order to alleviate the problem of water shortage in the North China Plain, China implemented the South-to-North Water Diversion Project, which mainly transports water resources from the Yangtze River to the North China Plain. It is the largest water diversion project in the world. The project officially started to deliver water at the end of 2014, which greatly improved the water supply capacity of water-scarce areas along the line. Therefore, it is very meaningful to analyze the extent to which the project can affect the groundwater reserves in the North China Plain. The present invention calculates the spatial and sequence changes of groundwater reserves in each year from 2011 to 2021. In order to facilitate the comparison of amplitude differences between years, the present invention performs first-order correction processing on the GWSA series of each year. The calculation results are as follows: Figure 8 shown.

[0063] like Figure 8 As shown in Figure 2, the South-to-North Water Diversion Project has significantly improved the groundwater quality of the North China Plain after its opening. From the spatial distribution characteristics, it can be seen that the North China Plain shows a characteristic of gradually decreasing GWSA from south to north. At the same time, there are obvious amplitude convex areas within Tianjin and Hebei Province ( Figure 8 (e)—(i)). It can be seen from the sequence diagram that after the South-to-North Water Diversion Project ( Figure 8 (e)-(k)) The annual amplitude of GWSA is significantly higher than that before the South-to-North Water Diversion Project ( Figure 8 (a)—(d)).

[0064] At the same time, there are obvious seasonal variations in the series of each year, indicating that changes in groundwater are closely related to the demand of agriculture and industry. Although the amplitude of GWSA has increased since 2015, it has shown a clear downward trend overall, which is more obvious than before the water transfer. The reason is that the demand for groundwater from industry and agriculture has become more urgent after the South-to-North Water Diversion Project, which also shows that the phenomenon of groundwater over-exploitation has become more serious. Figure 8 (g)-(i)) The spatial distribution of GWSA in the North China Plain is the most prominent, and its GWSA amplitude is significantly higher than that in other years. The phenomenon did not ease until 2020.

[0065] In summary, the South-to-North Water Diversion Project has played a certain role in alleviating the terrestrial water storage in the North China Plain.

[0066] Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the embodiments of the present invention. Therefore, the embodiments of the present invention are not limited to the embodiments shown herein, but are intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

[0067] The above are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for improving the accuracy of groundwater reserves based on data feature fusion inversion, characterized in that: This module includes a GNSS-based crustal vertical series seasonal term inversion module: it solves and preprocesses GNSS data, obtains a hydrological load deformation sequence by deducting atmospheric loads and non-ocean tidal loads, and combines the Green's load function with the crustal load model to invert the seasonal variation characteristic sequence of regional terrestrial water storage anomalies; Extract TWSA trend term module: Utilizes the equivalent water height variables provided by the GRACE Mascon dataset to extract regional terrestrial water storage anomalies, and obtains the terrestrial water storage trend term through an improved correlation variational mode decomposition algorithm. This trend term is then added to the GNSS inversion results to obtain the regional TWSA with trend. GWSA inversion module based on TWSA correction: By correcting the trend TWSA, including canopy water, snow water, and soil water, the abnormal changes in regional groundwater storage are obtained; Specific methods include: S1, using Green’s function to establish the relationship between load mass and deformation; S2, using the equivalent water height variables provided by the GRACE Mascon dataset to extract the abnormal changes in regional terrestrial water storage, and using the improved correlation variational mode decomposition algorithm to obtain the trend term of terrestrial water storage, this trend term is added to the GNSS inversion results to obtain the regional TWSA with trend; S3, the seasonal term and trend term of the TWSA obtained by inversion are fused to obtain the TWSA with trend term; S4, the regional GWSA is obtained by correcting the TWSA of the trend term, including canopy water, snow water, and soil water.

2. The method for improving groundwater storage accuracy based on data feature fusion inversion according to claim 1 is characterized in that: In step S1, the formula of Green's function for vertical load deformation of the crust is as follows: Where θ represents the angular radius from the center of the disk; P n represents Legendre polynomials; G represents Newton's gravitational constant; R represents the radius of the Earth; h n represents the load Love number; g represents the acceleration due to gravity; The Γ n The derivation of the function is as follows: Where θ represents the angular radius from the center of the disk; P represents the Legendre polynomial. When n = 0, the expression of the Γ function is as follows:

3. The method for improving groundwater storage accuracy based on data feature fusion inversion according to claim 1 is characterized in that: The obtained solution is regularized using the curvature smoothing algorithm and attached to the solution matrix as a set of constraints. For each period of this study, a restrained least squares problem is minimized to estimate the daily terrestrial water storage changes: TWSA Season =((Gx-b) / σ) 2 +b 2 (L(x)) 2 →min Among them, TWSA Season represents the seasonal characteristic sequence of TWSA obtained by inversion; G represents the Green’s function coefficient matrix; σ represents the standard deviation of the GNSS vertical displacement sequence; b represents the deformation observation sequence of the grid, where b represents the corrected GNSS vertical deformation sequence; L represents the Laplace operator; and β represents the smoothing factor.

4. The method for improving groundwater storage accuracy based on data feature fusion inversion according to claim 1 is characterized in that: In the step S2, Among them, TWSA Trend represents the trend characteristic sequence of TWSA extracted based on GRACE; represents the GRACE original TWSA sequence variable; represents the seasonal characteristic sequence of the TWSA sequence based on GRACE.

5. The method for improving groundwater storage accuracy based on data feature fusion inversion according to claim 1 is characterized in that: In the step S3, TWSA all =TWSA Season +TWSA Trend Among them, TWSA all represents the sum of the TWSA series, including trend and seasonal terms; TWSA Season represents the TWSA seasonal term obtained based on GNSS inversion; TWSA Trend Represents the sequence trend item extracted based on GRACE.

6. The method for improving groundwater storage accuracy based on data feature fusion inversion according to claim 1 is characterized in that: In the step S4, GWSA DFFIM =TWSA all -IN Can -IN soil -IN snow Among them, GWSA DFFIM represents the groundwater storage anomaly obtained by DFFIM inversion; TWSA represents the terrestrial water storage anomaly; W can Indicates the change of canopy water content; W soil Indicates the change in total soil water content from 0 to 200 mm.

Images