A method to improve the accuracy of terrestrial water storage anomalies based on deep learning weight loading

Abstract

The present invention discloses a method for improving the accuracy of terrestrial water storage anomaly based on deep learning weight load, comprising: step 1, dividing the study area into 1°×1° grids, searching whether each grid contains a GNSS station, if the grid contains a GNSS station, executing step 2, otherwise executing step 3; step 2, based on the GNSS station observation data, obtaining the crust deformation sequence in the observed grid; step 3, using the MEEMD decomposition method and the LSTM regression model to process and obtain the crust deformation sequence in the unobserved grid; step 4, using the NTAL model and the NTOL model, correcting the crust deformation sequence in all observed / unobserved grids for atmospheric and non-tidal ocean loads; step 5, using the crust load deformation sequence after all grid corrections as input data, combining the Green's function and the crust load model, inverting and outputting the terrestrial water storage anomaly TWSA of the study area. The present invention realizes the accurate inversion of TWSA.

Description

Technical Field

[0001] The present invention belongs to the cross-technical fields of satellite gravity, hydrology, etc., and in particular relates to a method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load. Background Art

[0002] Terrestrial water storage (TWS) refers to all the water stored on the surface and underground, including snow, glaciers, soil water, groundwater, rivers, lakes and biological water, etc. It is an important part of the water cycle. At the same time, terrestrial water storage resources are very scarce, accounting for only 3.47% of the world's total water resources, but they provide important guarantees for industry, agriculture and human life. China's total freshwater resources account for only 6% of the world's water resources, and China has a population of about 1.4 billion, resulting in per capita freshwater resources far below the world average. In addition, China's water resources have problems such as uneven annual distribution, prominent contradictions between water supply and demand, and low water resource utilization. In recent years, a series of natural disasters caused by uneven water supply and demand have occurred frequently (such as drought, floods, soil erosion, etc.), seriously affecting people's normal life and social economic development. Therefore, how to scientifically and effectively manage China's regional water resources is now an urgent issue.

[0003] In recent years, the optimization of hydrological models and the advancement of observation technology can accurately measure and continuously monitor the redistribution of terrestrial water at different temporal and spatial scales. Hydrological models are mathematical models of terrestrial hydrological processes and are widely used in climate change and human exploration of global water resources. However, hydrological models usually simplify the complex hydrological cycle process. Therefore, hydrological models cannot fully capture the complete hydrological components and tend to underestimate the changes in the terrestrial water cycle caused by climate and human beings. For example, the Noah model in the Global Land Data Assimilation System (GLDAS) only includes soil moisture, snow water equivalent and total canopy water storage at a depth of 0 to 2 m, while ignoring the influence of components such as surface water (lakes, rivers and reservoirs, etc.), deep groundwater and human factors. Since the redistribution of huge water masses will cause changes in the gravity field in the surrounding area, it makes it possible to invert the terrestrial water storage anomaly (TWSA) based on gravity anomaly data. In March 2002, the Gravity Recovery and Climate Experiment (GRACE) satellite, jointly developed by the United States and Germany, was successfully launched, providing an unprecedented means of observing surface mass migration. This observation method can accurately measure the Earth's gravity field and continuously monitor the mass migration of the surface. In recent years, a large number of workers have studied the surface mass migration of typical regions around the world based on GRACE satellite data and have achieved fruitful results, such as: the Amazon Basin, Greenland, the North China Plain of China, and Southwest China. However, due to the inherent characteristics of the GRACE satellite, the temporal and spatial resolution of the monitoring data is relatively coarse. The temporal resolution is on a monthly scale and the spatial resolution is 1°×1°. This feature greatly limits the inversion of GRACE satellites for small-scale regional hydrological loads. At the same time, due to the aging of GRACE satellite components, it was retired in 2017, and the next-generation gravity measurement satellite GRACE Follow-On (GRACE-FO) was launched in 2018. This has resulted in a gap of nearly one year between the GRACE satellite and the GRACE-FO satellite. Therefore, it is very important to find an alternative observation method to continuously monitor TWSA.

[0004] The periodic changes of huge water masses will also cause subtle deformations in the surrounding crust. It is possible to invert the changes in terrestrial water storage by continuously monitoring the subtle deformations of the crust. The Global Navigation Satellite System (GNSS) has many observation advantages, such as high precision, all-weather, and real-time. This method can accurately measure the subtle deformations of the earth's surface. In recent years, a large number of scholars have achieved many results in inverting TWSA in typical areas around the world based on GNSS vertical deformation sequences, such as California, the western mountains of the United States, and southwest China. In areas with dense GNSS stations, GNSS can effectively monitor the subtle deformations of the crust caused by surface loads, so that TWSA in a certain area can be monitored in near real time. This observation method and inversion strategy has great potential in monitoring hydrological signals and is helpful in establishing an early warning system for extreme hydrological and meteorological disasters such as storm tides. At the same time, the construction of the China Crustal Tectonic Network has been going on for 10 years, which makes it possible to monitor regional crustal deformation in real time. Many scholars have achieved a lot of results based on the crustal deformation data provided by CMONOC. However, due to the particularity of geological and climatic conditions, it is difficult to establish GNSS continuous observation stations in some areas, resulting in uneven spatial distribution of continuous GNSS observation stations, which greatly limits the application of GNSS inversion of crustal loads. Therefore, in areas with sparse GNSS stations, how to more accurately invert TWSA based on GNSS vertical deformation sequences has become a research hotspot in recent years. Summary of the invention

[0005] The technology of the present invention solves the problem: Overcoming the shortcomings of the prior art, providing a method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load, aiming to more accurately invert TWSA.

[0006] In order to solve the above technical problems, the present invention discloses a method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load, comprising:

[0007] Step 1: Divide the study area into 1°×1° grids, and search each grid to see if it contains a GNSS station. If the grid contains a GNSS station, define the grid as an observed grid and execute step 2; if the grid does not contain a GNSS station, define the grid as an unobserved grid and execute step 3.

[0008] Step 2: Based on the GNSS station observation data, the daily coordinates are solved to obtain the crustal vertical deformation sequence determined based on the GNSS station, and the sequence outliers and step terms greater than three times the standard deviation are removed to obtain the crustal deformation sequence within the observed grid;

[0009] Step 3, obtain the surface temperature sequence and pressure sequence of the unobserved grid, normalize the obtained surface temperature sequence and pressure sequence of the unobserved grid, and obtain the normalized result; use the MEEMD decomposition method to decompose the normalized result to obtain m feature sequences; use the GNSS vertical deformation sequence of the unobserved grid as the target sequence and the m feature sequences as the input sequence, and use the LSTM regression model to regress and obtain the crust deformation sequence in the unobserved grid;

[0010] Step 4, using the NTAL model and the NTOL model, the crustal deformation sequences in all observed grids obtained in step 2 and the crustal deformation sequences in all unobserved grids obtained in step 3 are corrected for atmospheric and non-tidal ocean loads to obtain the crustal load deformation sequences after correction for all grids;

[0011] Step 5: Using the crustal load deformation sequence corrected by all grids obtained in step 4 as input data, combining the Green's function and the crustal load model, the terrestrial water storage anomaly TWSA of the study area is inverted and output.

[0012] In the above method for improving the accuracy of terrestrial water storage anomaly based on deep learning weight load, the LSTM regression model includes: input gate i t 、Forget gate t , output gate o t 、Memory unit c t , and the input sequence corresponding to each step size: x1,x2,...,x t , t represents the step length; memory unit c t Through the input gate i t 、Forget gate t and output gate o t Control the memory and forgetting of data.

[0013] In the above method for improving the accuracy of terrestrial water storage anomaly based on deep learning weight load, the memory unit c t AND input gate i t 、Forget gate t and output gate o t The relationship between them is shown as follows:

[0014] f t =σ(W f x t +U f h t )

[0015] o t =σ(W o x t +U o h t )

[0016]

[0017] h t =o t ·tanh(c t )

[0018] Where, σ represents the sigma function; W f , W o and W c Respectively represent the weight matrices in the input process corresponding to the forget gate, output gate and memory unit; U f , U o and U c represents the state transfer weight matrix corresponding to the forget gate, output gate and memory unit, which is an S-type function; h t The hidden state vector representing the output; represents the updated memory unit; tanh represents the hyperbolic tangent function.

[0019] In the above-mentioned method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load, the decomposition formula used by the MEEMD decomposition method is as follows:

[0020] F=IMF1+IMF2+…+IMF m +noi w

[0021] Among them, F represents the original feature sequence, IMF1~IMF m represents the n feature sequences obtained by decomposing the original feature sequence F, noi w Represents the Gaussian white noise added to the sequence to be decomposed during the decomposition of the original feature sequence F.

[0022] In the above method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight loading, the crustal deformation sequence in the unobserved grid obtained by regression is denoted as D g , which is expressed as follows:

[0023]

[0024] Among them, d j represents the distance between the jth unobserved grid center and the GNSS station, n represents the number of GNSS stations, Represents the sum of the distances between the unobserved grid and each GNSS station, Net LSTM Represents an LSTM regression model.

[0025] In the above-mentioned method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load, the Green's function is used to indicate the relationship between crustal load and crustal load deformation.

[0026] In the above method for improving the accuracy of terrestrial water storage anomaly based on deep learning weight load, the Green function U green The representation is as follows:

[0027]

[0028] Where θ represents the angular radius from the center of the disk, P n represents the Legendre polynomial, G represents Newton's gravitational constant, R represents the radius of the earth, and h n represents the load Love number, and g represents the acceleration due to gravity.

[0029] In the above-mentioned method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load, the terrestrial water storage anomaly TWSA of the study area is inverted, including:

[0030] During each study period, a least squares problem was minimized to estimate the daily change in terrestrial water storage:

[0031] ((U green xd) / s) 2 +β 2 (L(x)) 2 →min

[0032] Where d represents the crustal load deformation sequence after correction for all grids, s represents the standard deviation of the grid hydrological load deformation sequence, β represents the smoothing factor, x represents the daily terrestrial water storage anomaly to be estimated, and L(·) represents the Laplace operator function.

[0033] The present invention has the following advantages:

[0034] Based on the long short-term memory network (LSTM), the inverse distance weighting method and the crustal load model, the present invention constructs a new deep learning weight load inversion method (DWLIM), and inverts the TWSA in the Chinese region from 2011 to 2020. DWLIM can be mainly divided into the following five steps: the first step is to search whether the 1°×1° grid contains a GNSS station. If it does, jump to the second step, otherwise jump to the third step; the second step is to solve and preprocess the GNSS vertical time series; the third step is to simulate the vertical deformation of the crust in the unknown grid based on the LSTM regression method; the fourth step is to use the non-tidal atmospheric load (NTAL) and non-tidal oceanic load (NTOL) model data to correct the vertical deformation of the crust in all grids for atmospheric load and non-ocean tidal load; the fifth step is to use the corrected sequences in all grids as input variables, combined with the crustal load model inversion to obtain the TWSA in the Chinese region. It can be seen that the DWLIM proposed in the present invention provides an important reference for the regression of vertical deformation of the crust and the monitoring of terrestrial water reserves, and can accurately invert TWSA. The research results have very important scientific significance and application value for the analysis of non-tectonic deformation of the crust, hydrodynamic research and sustainable management of regional water resources in China. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 is a flowchart of the steps of a method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load in an embodiment of the present invention;

[0036] Figure 2 It is a schematic diagram of the time span of continuous GNSS stations in China and surrounding areas in an embodiment of the present invention;

[0037] Figure 3 It is a verification effect diagram of a simulated crustal vertical deformation in an embodiment of the present invention; wherein 3(a) is a Taylor diagram of simulated crustal deformation, and 3(b) is a Chinese regional average sequence effect diagram;

[0038] Figure 4 4(a) is a diagram of the decomposition result of the temperature series after the G456 grid normalization processing, and 4(b) is a diagram of the decomposition result of the pressure series after the G456 grid normalization processing;

[0039] Figure 55(a) is a schematic diagram of simulation results of an unknown grid crustal deformation sequence in an embodiment of the present invention; wherein 5(a) is a simulation result of the G464 grid, 5(b) is a simulation result of the G740 grid, and 5(c) is a simulation result of the G456 grid;

[0040] Figure 6 It is a sequence effect diagram before and after correction corresponding to an atmospheric load correction and non-ocean tide correction in an embodiment of the present invention;

[0041] Figure 7 It is a time series effect diagram of TWSA in China based on DWLIM inversion in an embodiment of the present invention;

[0042] Figure 8 8(a) is a wavelet analysis diagram of DWLIM and traditional GNSS inversion TWSA, 8(b) is a wavelet analysis diagram of DWLIM and GLDAS, 8(c) is a wavelet analysis diagram of DWLIM and GRACE, and 8(d) is a sequence diagram of TWSA calculation results.

[0043] Fig. 9 9(a) is a DWLIM inversion result evaluation heat map in an embodiment of the present invention; wherein, 9(a) is the PCC index of the GRACE period, 9(b) is the NSE index of the GRACE period, 9(c) is the RMSE index of the GRACE period, 9(d) is the PCC index of the GRACE-FO period, 9(e) is the NSE index of the GRACE-FO period, and 9(f) is the RMSE index of the GRACE-FO period;

[0044] Fig.10 10(a) is a comparison diagram of the rainfall in China's ten major river basins and the TWSA results in an embodiment of the present invention; wherein 10(a) is the Yangtze River, 10(b) is the Southeast River Basin, 10(c) is the Haihe River Basin, 10(d) is the Huaihe River Basin, 10(e) is the Yellow River Basin, 10(f) is the Liaohe River Basin, 10(g) is the Songhua River Basin, 10(h) is the Northwest River Basin, 10(i) is the Southwest River Basin, and 10(j) is the Pearl River Basin. DETAILED DESCRIPTION

[0045] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments disclosed in the present invention will be further described in detail below with reference to the accompanying drawings.

[0046] One of the core ideas of the present invention is that the crust can be regarded as a spherical shell that can produce elastic deformation. When large-scale migration of terrestrial water occurs, the pressure of the huge water mass on the crust will also change accordingly. Due to the change in surface pressure, the crust will produce a nonlinear response of uplift or subsidence accordingly. At the same time, GNSS stations can accurately observe the deformation sequence of surface subsidence, and its observation accuracy can reach the millimeter level. In recent years, many scholars have used the crustal load deformation model to accurately invert the local TWSA in areas with dense distribution of GNSS stations. However, due to the influence of geographical station construction conditions, the distribution of GNSS stations around the world is extremely uneven. In areas where GNSS stations are sparsely distributed, the TWSA of the area cannot be accurately inverted using the GNSS deformation sequence. Therefore, how to accurately simulate the surface settlement sequence in areas lacking GNSS stations is the key to improving the accuracy of GNSS inversion of TWSA. Figure 1 The present invention combines LSTM, inverse distance weighting method and crustal load model to construct a new deep learning weight load inversion method, that is, a method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load. The specific process is as follows:

[0047] Step 1: Divide the study area into 1°×1° grids and search whether each grid contains a GNSS station.

[0048] In this embodiment, if the grid contains a GNSS station, the grid is defined as an observed grid, and step 2 is executed; if the grid does not contain a GNSS station, the grid is defined as an unobserved grid, and step 3 is executed.

[0049] Step 2: Based on the GNSS station observation data, the daily coordinates are solved to obtain the crustal vertical deformation sequence determined based on the GNSS station, and the sequence outliers and step terms greater than three times the standard deviation are removed to obtain the crustal deformation sequence within the observed grid.

[0050] In this embodiment, as mentioned above, GNSS stations can accurately observe the deformation sequence of surface settlement, and its observation accuracy can reach millimeter level. In recent years, many scholars have used crustal load deformation models to accurately invert local TWSA in areas with dense GNSS stations. That is, for the observed grid, the corresponding crustal deformation sequence can be obtained in any appropriate way, for example, based on the observation files, precise ephemeris files, navigation files and table files obtained by GNSS station observations, the daily coordinates are solved using GAMIT / GLOBK software to obtain the crustal vertical deformation sequence observed by GNSS, and the sequence outliers and step terms greater than three times the standard deviation are removed. This embodiment will not be repeated here.

[0051] Step 3: Obtain the surface temperature sequence and pressure sequence of the unobserved grid, and calculate the surface temperature sequence (Temperature Sequence, S T ) and the atmospheric pressure sequence (SAP) are normalized to obtain the normalized result; the normalized result is decomposed by MEEMD decomposition method to obtain m feature sequences; the GNSS vertical deformation sequence of the unobserved grid is taken as the target sequence and the m feature sequences are taken as the input sequences, and the LSTM regression model is used to regress the crustal deformation sequence in the unobserved grid.

[0052] In this embodiment, the LSTM regression model is an improved model of a recursive neural network (RNN) proposed by Hochreter et al. The LSTM regression model constructs a memory storage unit and trains the data using a time back propagation algorithm. This model can solve the problem of RNN having vanishing gradients and no long-term dependencies. The standard LSTM regression model specifically includes: input gate i t 、Forget gate t , output gate o t 、Memory unit c t , and the input sequence corresponding to each step size: x1,x2,...,x t , t represents the step length; memory unit c t Through the input gate i t 、Forget gate t and output gate o t Control the memory and forgetting of data.

[0053] Furthermore, the memory unit c t can be passed through different gates (input gate i t 、Forget gate t and output gate o t ) controls the memory and forgetting of data, as follows:

[0054] f t =σ(W f x t +U f h t )

[0055] o t =σ(W o x t +U o h t )

[0056]

[0057] h t =o t ·tanh(c t )

[0058] Where, σ represents the sigma function; W f , W o and W c Respectively represent the weight matrices in the input process corresponding to the forget gate, output gate and memory unit; U f , U o and U c represents the state transfer weight matrix corresponding to the forget gate, output gate and memory unit, which is an S-type function; h t The hidden state vector representing the output; represents the updated memory unit; tanh represents the hyperbolic tangent function. t 、Forget gate t and output gate o t ) jointly control the information entering and leaving the memory cell. The input regulates the new information entering the memory cell; the forget gate controls how much information is stored in the memory cell; and the output gate defines how much information can be output. The gate structure of the LSTM model makes the information on the time series form a balanced long-term and short-term dependency to achieve the purpose of multivariate regression.

[0059] Furthermore, in view of the problem of few feature sequences in this embodiment, this embodiment uses modified ensemble empirical mode decomposition (MEEMD) to decompose the original features to achieve the purpose of increasing the dimension of the feature sequence. The decomposition formula specifically used by the MEEMD decomposition method is as follows:

[0060] F=IMF1+IMF2+…+IMF m +noi w

[0061] Among them, F represents the original feature sequence, IMF1~IMF m represents the n feature sequences obtained by decomposing the original feature sequence F, noi w Represents the Gaussian white noise added to the sequence to be decomposed during the decomposition of the original feature sequence F.

[0062] Furthermore, since the crustal deformation has similar characteristics in a small area, the amplitude of the crustal deformation sequence closer to the GNSS station is more similar to it. The inverse distance average method is to determine the position coordinates of the remaining training stations with the point to be found as the center, and calculate the distance between each GNSS station and the point to be found, and assign weights according to the inverse of the distance. Therefore, the crustal deformation sequence in the unobserved grid obtained by regression is recorded as D g , which is expressed as follows:

[0063]

[0064] Among them, d j represents the distance between the jth unobserved grid center and the GNSS station, n represents the number of GNSS stations, Represents the sum of the distances between the unobserved grid and each GNSS station, Net LSTM Represents the LSTM regression model. Therefore, the simulated crustal vertical deformation of each grid is regressed n times, and the n regression results are weighted according to the inverse distance weight.

[0065] Step 4: Using the NTAL model and the NTOL model, the crustal deformation sequences in all observed grids obtained in step 2 and the crustal deformation sequences in all unobserved grids obtained in step 3 are corrected for atmospheric and non-tidal ocean loads to obtain the corrected crustal load deformation sequences for all grids.

[0066] Step 5: Using the crustal load deformation sequence corrected by all grids obtained in step 4 as input data, combining the Green's function and the crustal load model, the terrestrial water storage anomaly TWSA of the study area is inverted and output.

[0067] In this embodiment, since the earth's crust can be regarded as an elastic surface, when the mass of the surface changes, the surface will produce an elastic response of surface settlement or rebound, and this deformation is also called crustal load deformation. However, crustal load deformation is not only manifested in the vertical direction, but also in the horizontal direction. The crustal load deformation is more sensitive in the vertical direction, and its deformation amplitude is about 2 to 3 times that in the horizontal direction. Green's function U green The relationship between crustal load and crustal load deformation can be indicated as follows:

[0068]

[0069] Where θ represents the angular radius from the center of the disk, P n represents the Legendre polynomial, G represents Newton's gravitational constant, R represents the radius of the earth, and h n represents the load Love number, and g represents the acceleration due to gravity.

[0070] Furthermore, DWLIM uses all grid-corrected hydrological load crustal deformation as input data, and combines it with the crustal load model to invert the TWSA in the study area. The obtained solution is regularized using the curvature smoothing algorithm and added as a constraint in the solution matrix. In other words, the least squares problem is minimized within each study period to estimate the daily change in terrestrial water storage:

[0071] ((U green xd) / s) 2 +β 2 (L(x)) 2 →min

[0072] Where d represents the crustal load deformation sequence after correction for all grids, s represents the standard deviation of the grid hydrological load deformation sequence, β represents the smoothing factor, x represents the daily terrestrial water storage anomaly to be estimated, and L(·) represents the Laplace operator function.

[0073] Based on the above embodiment, a complete example process is described below.

[0074] Validation of a novel deep learning weight load inversion method

[0075] 1.1 Data and Model

[0076] 1.1.1 GNSS data

[0077] This embodiment uses the vertical observation sequence of the Continuously Operating Reference Stations (CORS) in the Chinese region provided by the Chinese mainland environmental monitoring network as the data basis. Since each CORS station is established at a different time, the time span of each station is also inconsistent. The time span is as follows: Figure 2 As shown. Therefore, in order to ensure the integrity of the vertical deformation sequence, the time span of this embodiment is selected as 2011-2020. There are 268 original CORS stations in and around China. After removing 5 stations with large time span differences, there are 263 sequences applied. The solution of the GNSS observation sequence is based on observation files, navigation files, precise ephemeris files and table files. The daily coordinate solution file is obtained by solving the GAMIT / GLOBK 10.4 software. Secondly, the GNSS vertical time series is preprocessed, which includes removing observation outliers greater than three times the mean error and removing sequence step terms caused by earthquakes or antenna replacement.

[0078] 1.1.2 GRACE Mascon Data

[0079] GRACE is a joint space mission of the National Aeronautics and Space Administration (NASA) and the German Aerospace Center (Deutsches Zentrum für Luft-und Raumfahrt, DLR). The project launched two GRACE satellites in March 2002 and announced the termination of the mission in 2017. The main task of the GRACE satellite is to monitor the spatiotemporal changes of the Earth's gravity field on a global scale. The changes in the gravity field are mainly caused by earthquakes, glacial isostatic adjustments, ocean and hydrological changes. In China, the gravity changes inverted by GRACE are generally attributed to large-scale hydrological changes.

[0080] In order to verify the reliability of the DWLIM inversion results, this embodiment uses GRACE Mascon to compare with the DWLIM inversion results. However, due to the different solution strategies of each institution, there is a large uncertainty in the single GRACE-M result. Therefore, this embodiment uses the GRACE Mascon (GRACE-M) data products from 2011 to 2020 provided by the Center for Space Research (CSR) of the University of Texas and the Jet Propulsion Laboratory (JPL) of NASA, extracts the TWSA of the Chinese region according to the boundary file, and takes the average of the two products as the final GRACE-M result ΔTWSA GRACE-M :

[0081]

[0082] Among them, ΔMascon CSR It indicates the TWSA of China extracted based on CSR, and it indicates the TWSA of China extracted based on JPL.

[0083] 1.1.3 GLDAS data products

[0084] In order to construct and verify DWLIM, this embodiment uses the temperature variable in the Global Land Data Assimilation System (GLADAS) V2.2 global land surface model provided by NASA as the input data of LSTM regression, and the TWSA variable in the GLDAS V2.2 model is used as the verification data of the DWLIM inversion result. Among them, GLDAS V2.2 is the evolution of the early product Catchment Land Surface Model (CLSM), including 24 variables including temperature and terrestrial water storage anomalies. The time span of the GLDAS V2.2 hydrological model is from February 2003 to the present, with a spatial resolution of 0.25°×0.25°, a temporal resolution of daily scale, and a spatial coverage latitude range of 60°S~90°N and a longitude range of 180°W~180°E.

[0085] 1.1.4 ERA data products

[0086] This example uses the surface pressure variable in the ERA5 data product provided by the European Centre for Medium-Range WeatherForecasts (ECMWF) as the input data for LSTM regression. Among them, the spatial resolution of the surface pressure data provided by ERA5 on a daily scale is 0.1°×0.1°, the time span is from 2000 to the present, and the spatial coverage is global. The specific data details of the ERA dataset can be found on the website (https: / / www.ecmwf.int / en / forecasts / dataset / ecmwf-reanalysis-v5).

[0087] 1.2 Evaluation Indicators

[0088] This embodiment uses three evaluation indicators, namely, root mean square error (RMSE), Pearson's Correlation Coefficient (PCC), and Nash-Sutcliffe efficiency coefficient (NSE), to evaluate the accuracy of the DWLIM inversion results. The calculation formula is as follows:

[0089]

[0090] Among them, Y represents the real sequence, X represents the inversion result, and Represent the average values ​​of Y and X respectively, and k represents the number of discrete points in the sequence. RMSE can be used to evaluate the discrete deviation between the inversion result and the true value. The smaller the RMSE value, the more accurate the inversion result. NSE is mainly used to evaluate the quality of the hydrological model. Its value is less than or equal to 1. The larger its value, the better the hydrological model. When NSE is close to 0, it indicates that the effect of the hydrological model is close to the average level of the observed value. PCC is mainly used to describe the linear correlation between two sequences. The PCC value is between -1 and 1. The closer its value is to 1, the more reliable the inversion result.

[0091] 1.3. Calculation of TWSA in China based on DWLIM

[0092] 1.3.1. Verification of LSTM simulation of unknown grid crust deformation

[0093] In order to verify the accuracy of DWLIM simulation of grid crustal deformation, the stations with correlations between pressure and temperature series and GNSS series greater than 0.5 within the grid range were selected, and 75 GNSS station sequences were obtained as control sequences. And 263 sequences were used for verification. When there was a control station in the grid where the GNSS station was located, it was regressed 74 times; when there was no GNSS station in the grid, it was regressed 75 times, and the GNSS sequence was used as the true value to compare with the simulated crustal deformation sequence. First, the crustal influencing factor sequence of all grids in the study area was extracted, including two characteristic variables of pressure and temperature. Then, the two characteristic sequences were decomposed into 10 characteristic mode components respectively using the MEEMD decomposition method to achieve the purpose of increasing the dimension. Secondly, the decomposed mode components were used as input data, and the GNSS vertical sequence was used as output data. The LSTM method was used for regression, and the 75 simulation results were weighted according to the inverse of the distance. Finally, the simulated crustal deformation sequence of each grid where the GNSS station was located was obtained. The simulation results are evaluated based on RMSE and PCC indicators. Figure 3 shown.

[0094] according to Figure 3 It can be seen that the RMSE values ​​of most station sequences are within 5 mm. However, due to the differences in observation quality between GNSS sequences, the RMSE values ​​of some stations are too large. After further statistics, it can be seen that 68.63% of the comparison results are within 6 mm. Then the PCC index is used to evaluate the consistency relationship between the simulated sequence and the true value. The maximum PCC value can reach 0.87, and its mean is 0.53. Figure 3 (b) shows the effect of averaging the simulated and real sequences in the Chinese region. It can be seen that the third step in DWLIM can effectively simulate the periodicity of the sequence, and the sequence simulated based on LSTM is smoother and the characteristics of the sequence are more obvious.

[0095] 1.3.2. Obtaining hydrological load deformation

[0096] (3.2.1) Simulating crustal deformation in an unknown grid

[0097] This embodiment selects temperature and air pressure as the input data of the LSTM regression method. In order to make full use of the sequence characteristics of the input data, this embodiment uses the MEEMD decomposition method to decompose the surface temperature and air pressure series into IMF1~IMF 10 Taking the G456 grid as an example, the decomposition effect of the input sequence is as follows: Figure 4 shown.

[0098] Figure 4 (a) and Figure 4 The IMF1 component in (b) is the normalized original surface temperature and pressure series, and IMF1~IMF 10 is the feature item from high frequency to low frequency obtained by decomposition. The decomposition result can better reflect the trend item, seasonal item and residual item of the sequence. The temperature and air pressure sequences both show the periodic characteristics of the whole year. In the LSTM regression method, Figure 4 The 10 IMF components decomposed in the above method are used as input sequences, and the 75 control GNSS vertical deformation sequences are used as output data. At the same time, the inverse distance weighting is performed according to the distance between the grid center and the 75 stations, and finally the crust vertical deformation sequence of the unknown grid is simulated. The simulation results of the unknown grid are shown in Figure 2. Figure 5 As shown in (a) to 5(c), the G464, G740 and G456 grids are taken as examples.

[0099] The GNSS stations in mainland China are unevenly distributed and the number of stations is small, which cannot achieve the overall coverage of the crust of mainland China. Therefore, it is very necessary to simulate the vertical deformation of the unknown grid crust. Figure 5 (a) to (c) show the crustal vertical deformation results obtained by using the 20 IMF components obtained by MEEMD decomposition as input data and 75 control GNSS stations as output data, combined with the inverse distance weighted theory and LSTM regression method. The simulation results show that the periodic term and annual amplitude of the crustal vertical deformation can be well simulated according to this strategy, providing a good data basis for the inversion of TWSA.

[0100] (3.2.2) Correct all grid deformation sequences

[0101] There are load deformations caused by atmospheric and non-ocean tidal loads in the GNSS vertical deformation sequence. Therefore, in order to extract the crustal deformation caused by hydrological loads, the present invention uses the NTAL and NTOL models provided by the Potsdam Earth System Modeling Group (http: / / esmdata.gfz-potsdam.de:8080 / repository) as correction data, and adds the two correction sequences to the crustal load deformation sequence in the Chinese region, including the GNSS vertical time series and the DWLIM method to simulate the load deformation sequence of the unknown grid. The vertical deformation sequences before and after correction are as follows: Figure 6 shown.

[0102] according to Figure 6 Further statistical data show that the mean values ​​of the load deformation amplitudes of NATL and NOTL are 3.62mm and 0.22mm respectively. The annual amplitude of the non-ocean atmosphere correction is mainly located in northern and eastern China, with a maximum amplitude of 5.5mm; while the annual amplitude of the non-ocean tide correction is smaller, with a maximum amplitude of only 1.5mm, mainly distributed in the eastern coastal areas of China. Figure 6 It can be seen that there are slight changes in the amplitude and phase of the corrected sequence, which provides a more accurate hydrological load deformation sequence for DWLIM inversion of TWSA.

[0103] 1.3.3 Inversion of TWSA in China based on DWLIM

[0104] This example first calculates the Green's function matrix of the point load in China, selects 0.25°×0.25° as the spatial resolution of the inversion, where the disk radius corresponding to the 0.25° spatial resolution is 13.90km, and selects the preset parameters of the extended boundary range of 2° and β=0.01. Using the vertical deformation sequence of the crust corrected by all 1°×1° grids in China as input data, the daily scale land water storage change value in China is obtained by inversion combined with the following formula:

[0105]

[0106] ((U green xd) / s) 2 +β 2 (L(x)) 2 →min

[0107] In order to better verify the accuracy of the DWLIM inversion results, the present invention performs first-order correction processing on the inversion results. The calculated average TWSA time series for the Chinese region is as follows: Figure 7 shown.

[0108] according to Figure 7It can be seen that the annual characteristics and amplitude of the TWSA sequence in the Chinese region can be effectively inverted based on the DWLIM method. The prominent position of the annual amplitude can be calculated by this method, which is located in Yunnan Province, the southern part of the Tibet Autonomous Region and the southern part of the North China Plain, and the result is consistent with the previous inversion results. In order to verify the accuracy of the DWLIM inversion results, this embodiment compares this inversion result with the GRACE-M results, the GLDAS hydrological data set results and the results of the traditional GNSS inversion of TWSA.

[0109] 1.4. Verification of the new deep learning weight load inversion method

[0110] 1.4.1. Spatial feature verification of TWSA results

[0111] In order to compare the accuracy of DWLIM inversion of TWSA in China, this example uses the traditional GNSS TWSA inversion method (TRA GNSS ), GRACE Mascon processing results and GLDAS hydrological results were used to extract the TWSA in China, and the annual amplitude of each result was calculated.

[0112] Based on the DWLIM strategy, the areas with annual amplitude bulge in China can be effectively inverted, and it can be clearly seen that the DWLIM inversion results have obvious bulges in the North China Plain, and the results are stronger than GRACE and GLDAS. The reason may be that the crustal deformation in this area is obvious, and the hydrological load deformation sequence cannot be more accurately extracted by only correcting the atmospheric load and non-ocean tides. At the same time, the amplitude bulge areas in northern Xinjiang and northern Heilongjiang are also inverted based on the DWLIM method. However, the spatial distribution of the annual amplitude obtained by the traditional GNSS inversion TWSA method has a spot feature. The reason is that the distance of the disk radius is limited, and only the TWSA situation within a limited range around the GNSS station can be inverted. In addition, the traditional GNSS inversion TWSA method causes signal loss in areas with sparse GNSS stations through smoothing. In summary, by simulating the crustal deformation of the unknown grid, the limitation of the disk radius on the GNSS inversion TWSA can be weakened.

[0113] 1.4.2. Verification of the timing characteristics of the results

[0114] In order to verify the time series reliability of DWLIM inversion TWSA, this embodiment compares the DWLIM inversion results with the traditional GNSS inversion TWSA method, GRACE and GLDAS results. In order to further analyze the sequence phase relationship between the DWLIM inversion results and the traditional GNSS inversion methods, GLDAS and GRACE, this embodiment draws the cross wavelet analysis diagrams of DWLIM-Traditional, DWLIM-GLDAS and DWLIM-GRACE, as shown in Figure 2. Figure 8 (a) to 8(c). At the same time, Figure 8 (d) shows the time series diagram of TWSA results extracted over China using DWLIM, traditional GNSS, GRACE, and GLDAS.

[0115] according to Figure 8 From (a) to (c), we can see that the DWLIM inversion results of TWSA are in good agreement with the traditional GNSS inversion results, GLDAS and GRACE phases. The DWLIM method can effectively invert the annuality and annual amplitude of the TWSA sequence, and is consistent with the GRACE and GLDAS calculation results ( Figure 8 (d)). However, due to the different observation methods, the annual amplitude of DWLIM is slightly larger than that of other methods. The reason is that the crustal deformation caused by other loads exists in the corrected crustal deformation sequence. Since the deformation caused by other loads is difficult to quantify, it is impossible to separate a single hydrological load deformation sequence. In addition, the TWSA results obtained by DWLIM inversion can be used to analyze the blank period between GRACE and GRACE-FO ( Figure 8 (d) The gray shaded area is filled. The seasonality of the DWLIM inversion results is more obvious than that of the traditional GNSS inversion TWSA results. In order to quantify the advantages of DWLIM over the traditional GNSS inversion TWSA, this embodiment uses PCC, NSE and RMSE to evaluate the inversion results. The evaluation results are shown in Figure 2. Fig. 9 shown.

[0116] according to Fig. 9 It can be seen that the DWLIM inversion TWSA results are significantly better than the traditional GNSS inversion results. In the GRACE period (2011-2017), the PCC, NSE and RMSE indicators of the DWLIM inversion results reached 0.81, 0.62 and 2.18 cm respectively; in the GRACE-FO period (2018-2020), the PCC, NSE and RMSE of the DWLIM inversion results reached 0.71, 0.49 and 2.4 cm respectively. The results show that the inversion results of DWLIM are more consistent with the GLDAS results. The reason is that the resolution of GRACE is on a monthly scale, and the rough time resolution will cause signal loss. Further statistics of the data show that the DWLIM inversion TWSA results are improved by an average of 67.11%, 128.15% and 22.75% in PCC, NSE and RMSE indicators respectively compared with the traditional GNSS inversion method. This result further demonstrates that the DWLIM method can effectively invert the changes in terrestrial water storage in areas with sparse GNSS stations, and is superior to the traditional GNSS inversion TWSA results in terms of temporal and spatial characteristics.

[0117] 2. Application of the new deep learning weight load inversion method

[0118] 2.1 Analysis of spatiotemporal changes of TWSA in various river basins in China based on rainfall data

[0119] This embodiment has verified that DWLIM can effectively invert the TWSA in the Chinese region, and can detect the convex point area of ​​the abnormal annual amplitude of the terrestrial water storage. The vertical load deformation of the crust is mainly affected by the terrestrial water storage. When the terrestrial water storage load increases, the crust shows a downward trend; conversely, when the terrestrial water storage load decreases, the crust shows an upward rebound trend. In order to study the causes of terrestrial water storage in China, this embodiment combines the monthly rainfall products provided by the China Meteorological Administration (CMA) and extracts the rainfall values ​​and TWSA of China's ten major river basins. Among them, TWSA is obtained by the DWLIM method, and its comparison effect diagram is shown as follows. Fig.10 shown.

[0120] according to Fig.10 It can be seen that in terms of annual amplitude, the annual amplitude of TWSA is generally positively correlated with the annual amplitude of rainfall results. Fig.10 (g)) and Liaohe River Basin ( Fig.10 (f)) is significantly higher than that of other basins, and the corresponding TWSA amplitude is also significantly higher than that of other basins. In terms of phase, the phase relationship between the terrestrial water storage and rainfall sequence in China is good. When the rainfall histogram peaks, the TWSA time series is also located at the peak. This further demonstrates the reliability of the DWLIM inversion of TWSA in China in terms of phase. At the same time, it can be seen that there is high-frequency noise in the time series of the TWSA inversion results, which also affects the inversion or prediction of TWSA. The reason is that in the process of GNSS observation of crustal deformation, due to the inherent properties of the receiver and the observation method, there is a lot of system noise caused by the ionosphere, troposphere, clock error or multipath effect. Therefore, in subsequent research, the focus will be on the noise classification and removal of GNSS vertical sequences to provide a more "clean" sequence for the inversion of TWSA.

[0121] In summary, the present invention discloses a method for improving the accuracy of terrestrial water storage anomaly based on deep learning weight load, and the terrestrial water storage anomaly (TWSA) can be effectively inverted based on densely distributed global navigation satellite system (GNSS) continuous observation stations. However, affected by natural environmental conditions, the GNSS continuous observation stations are extremely unevenly distributed around the world, which greatly limits the relevant research on inverting TWSA using GNSS vertical deformation. The present invention simulates the vertical deformation that does not include the GNSS station grid to achieve the purpose of improving the accuracy of GNSS inversion TWSA. The results show that the maximum Pearson correlation coefficient PCC, Nash efficiency coefficient NSE and root mean square error RMSE of the DWLIM inversion result with GRACE and GLDAS are 0.81, 0.61 and 2.18 cm respectively, which is an average improvement of 67.11%, 128.15% and 22.75% in PCC, NSE and RMSE indicators compared with the traditional GNSS inversion method. The results show that in areas with uneven distribution of GNSS stations, the DWLIM inversion strategy can effectively improve the accuracy of TWSA inversion. The TWSA of China's ten major river basins was obtained by DWLIM inversion, and the changes of TWSA in the ten major river basins were analyzed in combination with rainfall data. The DWLIM inversion results and rainfall data have good consistency in phase, and the peaks are located in June and July each year; secondly, the amplitude of the TWSA sequence of each river basin is positively correlated with the amplitude of the rainfall sequence, which further verifies the reliability of TWSA inversion based on DWLIM.

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

[0123] 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 improving the accuracy of terrestrial water storage anomalies based on deep learning weight load, characterized in that: include: Step 1: Divide the study area into 1°×1° grids, and search each grid to see if it contains a GNSS station. If the grid contains a GNSS station, define the grid as an observed grid and execute step 2; if the grid does not contain a GNSS station, define the grid as an unobserved grid and execute step 3. Step 2: Based on the GNSS station observation data, the daily coordinates are solved to obtain the crustal vertical deformation sequence determined based on the GNSS station, and the sequence outliers and step terms greater than three times the standard deviation are removed to obtain the crustal deformation sequence within the observed grid; Step 3, obtaining a surface temperature sequence and an air pressure sequence of an unobserved grid, and normalizing the obtained surface temperature sequence and air pressure sequence of the unobserved grid to obtain a normalized result; The normalized results are decomposed using the MEEMD decomposition method to obtain m feature sequences; The GNSS vertical deformation sequence of the unobserved grid is used as the target sequence, and m feature sequences are used as input sequences. The LSTM regression model is used to regress the crustal deformation sequence in the unobserved grid. Step 4, using the NTAL model and the NTOL model, the crustal deformation sequences in all observed grids obtained in step 2 and the crustal deformation sequences in all unobserved grids obtained in step 3 are corrected for atmospheric and non-tidal ocean loads to obtain the crustal load deformation sequences after correction for all grids; Step 5, using the crustal load deformation sequence corrected by all grids obtained in step 4 as input data, combining Green's function and crustal load model, invert and obtain the terrestrial water storage anomaly TWSA of the study area, and output it; LSTM regression model, including: input gate i t 、Forget gate t , output gate o t 、Memory unit c t , and the input sequence corresponding to each step size: x1,x2,...,x t , t represents the step length; memory unit c t Through the input gate i t 、Forget gate t and output gate o t Control the memory and forgetting of data; Memory unit c t AND input gate i t 、Forget gate t and output gate o t The relationship between them is shown as follows: f t =σ(W f x t +U f h t ) the t =σ(W o x t +U o h t ) h t =o t ·tanh(c t ) Where σ represents the sigma function; W f , W o and W c Respectively represent the weight matrices in the input process corresponding to the forget gate, output gate and memory unit; U f , U o and U c represents the state transfer weight matrix corresponding to the forget gate, output gate and memory unit, which is an S-type function; h t The hidden state vector representing the output; represents the updated memory unit; tanh represents the hyperbolic tangent function; The decomposition formula used by the MEEMD decomposition method is as follows: F=IMF1+IMF2+…+IMF m +we w Among them, F represents the original feature sequence, IMF1~IMF m represents the n feature sequences obtained by decomposing the original feature sequence F, noi w represents the Gaussian white noise added to the sequence to be decomposed during the decomposition of the original feature sequence F; The crustal deformation sequence in the unobserved grid obtained by regression is denoted as D g , which is expressed as follows: Among them, d j represents the distance between the jth unobserved grid center and the GNSS station, n represents the number of GNSS stations, Represents the sum of the distances between the unobserved grid and each GNSS station, Net LSTM Represents an LSTM regression model.

2. The method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load according to claim 1, characterized in that: Green's function is used to indicate the relationship between crustal load and crustal load deformation.

3. The method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load according to claim 2, characterized in that: Green's function U green The representation is as follows: Where θ represents the angular radius from the center of the disk, P n represents the Legendre polynomial, G represents Newton's gravitational constant, R represents the radius of the earth, and h n represents the load Love number, and g represents the acceleration due to gravity.

4. The method for improving the accuracy of terrestrial water storage anomalies based on deep learning weight load according to claim 3 is characterized in that: The inversion results in the terrestrial water storage anomaly TWSA in the study area, including: During each study period, a least squares problem was minimized to estimate the daily change in terrestrial water storage: ((U green x-d) / s) 2 +β 2 (L(x)) 2 →min Where d represents the crustal load deformation sequence after correction for all grids, s represents the standard deviation of the grid hydrological load deformation sequence, β represents the smoothing factor, x represents the daily terrestrial water storage anomaly to be estimated, and L(·) represents the Laplace operator function.

Images