CN114676621A - Method for improving accuracy of land water reserve abnormity based on deep learning weight load - Google Patents

Abstract

The invention discloses a method for improving the accuracy of land water reserve abnormity based on deep learning weight load, which comprises the following steps: step 1, dividing a research area into grids of 1 degree multiplied by 1 degree, searching whether each grid contains GNSS stations, if the grids contain the GNSS stations, executing step 2, otherwise executing step 3; step 2, obtaining a crustal deformation sequence in the observed grid based on the observation data of the GNSS observation station; step 3, processing by using a MEEMD decomposition method and an LSTM regression model to obtain a crustal deformation sequence in an unobserved grid; step 4, utilizing the NTAL model and the NTOL model to correct atmospheric and non-tidal ocean loads of all crustal deformation sequences in the observed/unobserved grids; and 5, taking the crustal load deformation sequences after all grids are corrected as input data, combining the Green function and the crustal load model, performing inversion to obtain TWSA (terrestrial water reserves) with abnormal land water reserves in the research area, and outputting the TWSA. The invention realizes accurate inversion of TWSA.

Description

Method for improving accuracy of land water reserve abnormity based on deep learning weight load

Technical Field

The invention belongs to the cross technical field of satellite gravimetry, hydrology and the like, and particularly relates to a method for improving the accuracy of land water reserve abnormity based on deep learning weight load.

Background

The Terrestrial Water Storage (TWS) refers to all Water stored on the earth surface and underground, including snow, glaciers, soil Water, underground Water, rivers, lake Water, and biological Water, and is an important component of Water circulation. Meanwhile, land water reserves are in short supply and only account for 3.47 percent of the total water resources in the world, and important guarantee is provided for the life of industry, agriculture and human beings. The total amount of fresh water resources in China only accounts for 6% of the global water resources, and China has about 14 hundred million population, so that the per-capita fresh water resource amount is far lower than the average level in the world, and the Chinese water resources have the problems of non-uniform annual distribution, outstanding contradiction between water resource supply and demand, low water resource utilization rate and the like. In recent years, a series of natural disasters (such as drought, flood, water and soil loss and the like) caused by uneven supply and demand of water resources frequently occur, and the normal life of people and the economic development of society are seriously influenced. Therefore, how to scientifically and effectively manage water resources in China is an urgent problem nowadays.

In recent years, advances in optimization and observation techniques for hydrological models have allowed accurate measurement and continuous monitoring of the redistribution of land water at different time 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 generally simplify complex hydrological cycle processes. Thus, the hydrological model does not capture the complete hydrological components completely and tends to underestimate climatic and human-induced land water cycle changes. For example, the Noah model in the Global Land Data Acquisition System (GLDAS) only includes soil humidity, snow water equivalent, and total canopy water storage at a depth of 0 to 2m, and ignores the influence of components such as surface water (lakes, rivers, reservoirs, and the like), deep groundwater, and human factors. Since the redistribution of huge Water masses can cause the change of the gravity field of surrounding areas, the method makes it possible to invert the land Water Storage Anomaly (TWSA) based on the gravity Anomaly data. Successful launch of Gravity Recovery and Climate Experiment (GRACE) satellites developed by United of America in 2002 by 3 months provides an unprecedented observation means for surface mass migration. The observation means can accurately measure the gravity field of the earth and can continuously monitor the mass transfer of the earth surface. In recent years, a large number of workers have studied surface quality migration in typical areas of the world based on GRACE satellite data and have achieved fruitful results, such as: amazon basin, greenland, north china plain, and the southwest china. However, due to the inherent nature of the Grace satellite, the spatial and temporal resolution of the monitored data is coarse. Where the temporal resolution is on the monthly scale and the spatial resolution is 1 deg. × 1 deg., this property greatly limits the inversion of the hydrologic loads of the Grace satellite for small scale regions. Meanwhile, due to aging of the GRACE satellite elements, the next generation of gravity measurement satellite GRACE focus-On (GRACE-FO) was launched in 2017 and 2018 after retirement. This leaves a gap of nearly 1 year between the GRACE satellite and the GRACE-FO satellite. Therefore, it is very important to find an alternative observation means to continuously monitor TWSA.

The periodic change of huge water masses can also cause fine deformation of surrounding crustae, and the land water reserve change can be inverted by continuously monitoring the fine deformation of the crustae. The Global Navigation Satellite System (GNSS) has many observation advantages, such as high precision, all weather, real-time, etc., and the means can accurately measure the fine deformation of the earth surface. In recent years, numerous scholars have achieved efforts to invert TWSA in a typical region of the world based on GNSS vertical deformation sequences, for example: california, the western mountains of the united states, the southwest region of china, and the like. In an area with dense GNSS stations, the GNSS can effectively monitor the fine deformation of the earth crust caused by earth surface load, so that the TWSA in a certain area can be monitored in near real time. The observation means and the inversion strategy have great potential in the aspect of monitoring the hydrological signals and are beneficial to establishing an early warning system for extreme hydrological disasters such as storm sea tides and the like. Meanwhile, the construction of the Chinese crust structural net has been for 10 years, which makes it possible to monitor the deformation of the regional crust in real time. Many scholars have achieved a great deal of work based on crust deformation data provided by CMONOC. However, due to the particularity of geological climate conditions, it is difficult to establish a GNSS continuous observation station in some areas, so that the spatial distribution of the continuous GNSS observation station is uneven, which also greatly limits the application of GNSS to inversion of earth-crust loads. Therefore, in an area with sparse GNSS stations, how to more accurately invert TWSA based on GNSS vertical deformation sequences becomes a research hotspot in recent years.

Disclosure of Invention

The technical solution of the present invention is: the method overcomes the defects of the prior art, improves the accuracy of land water reserves abnormality based on deep learning weight load and aims to invert TWSA more accurately.

In order to solve the technical problem, the invention discloses a method for improving the accuracy of land water reserve abnormity based on deep learning weight load, which comprises the following steps:

step 1, dividing a research area into grids of 1 degree multiplied by 1 degree, searching whether each grid contains GNSS stations, if the grids contain the GNSS stations, defining the grids as observed grids, and executing step 2; if the grid network does not contain the GNSS survey station, defining the grid network as an unobserved grid network, and executing the step 3;

step 2, resolving the coordinates every day based on observation data of the GNSS observation station to obtain a crustal vertical deformation sequence determined based on the GNSS observation station, and removing sequence abnormal values and step items larger than three times of standard deviation to obtain a crustal deformation sequence in the observed grid;

step 3, acquiring a surface temperature sequence and an air pressure sequence of the unobserved grid, and performing normalization processing on the acquired surface temperature sequence and the acquired air pressure sequence of the unobserved grid to obtain a normalization result; decomposing the normalization result by using a MEEMD decomposition method to obtain m characteristic sequences; the GNSS vertical deformation sequence of the unobserved grid is taken as a target sequence, the m characteristic sequences are taken as input sequences, and an LSTM regression model is utilized to carry out regression to obtain a crustal deformation sequence in the unobserved grid;

Step 4, utilizing an NTAL model and an NTOL model to correct atmospheric and non-tidal ocean loads of all crustal deformation sequences in the observed lattices obtained in the step 2 and all crustal deformation sequences in the unobserved lattices obtained in the step 3 to obtain crustal load deformation sequences after all the lattices are corrected;

and 5, taking all the corrected crustal load deformation sequences of the grids obtained in the step 4 as input data, combining the Green function and the crustal load model, performing inversion to obtain TWSA (terrestrial water reserves) with abnormal land water reserves in the research area, and outputting the TWSA.

In the method for improving the accuracy of the land water reserve abnormality based on the deep learning weight load, the LSTM regression model includes: input door i_tForgetting door f_tAnd an output gate o_tMemory cell c_tAnd the input sequence corresponding to each step length is as follows: x is the number of₁,x₂,…,x_tAnd t represents a step size; memory cell c_tThrough an input gate i_tForgetting door f_tAnd an output gate o_tAnd controlling the memory and forgetting of data.

In the method for improving the accuracy of the land water reserve abnormality based on the deep learning weight load, the memory unit c_tAnd input gate i_tForgetting door f_tAnd an output gate o_tThe relationship between them is expressed as follows:

f_t＝σ(W_fx_t+U_fh_t)

o_t＝σ(W_ox_t+U_oh_t)

h_t＝o_t·tanh(c_t)

where σ represents a sigma function; w _f、W_oAnd W_cRespectively representing the weight matrixes in the input process corresponding to the forgetting gate, the output gate and the memory unit; u shape_f、U_oAnd U_cThe expression respectively represents the state transition weight matrixes corresponding to the forgetting gate, the output gate and the memory unit, and the state transition weight matrixes are S-shaped functions; h is_tA hidden state vector representing an output;

representing the updated memory cell; tanh represents a hyperbolic tangent function.

In the method for improving the accuracy of the land water reserve abnormality based on the deep learning weight load, the MEEMD decomposition method adopts the following decomposition formula:

F＝IMF₁+IMF₂+…+IMF_m+noi_w

wherein F represents the original characteristic sequence, IMF₁～IMF_mRepresenting n signature sequences obtained by decomposing the original signature sequence F, noi_wWhich represents the white gaussian noise added to the sequence to be decomposed during the decomposition of the original signature sequence F.

In the method for improving the accuracy of land water reserve abnormity based on deep learning weight load, a crustal deformation sequence in an unobserved grid obtained by regression is marked as D_gExpressed as follows:

wherein d is_jRepresents the distance between the jth unobserved grid center and the GNSS stations, n represents the number of GNSS stations,

representing the sum of the distances, Net, of the unobserved grid from each GNSS station_LSTMThe LSTM regression model is represented.

In the method for improving the accuracy of land water reserve abnormity based on deep learning weight load, the green function is used for indicating the relation between the crustal load and the crustal load deformation.

In the method for improving the accuracy of land water reserve abnormity based on deep learning weight load, the Green function U_greenIs represented as follows:

where θ represents the angular radius from the center of the disk, P_nExpressing Legendre polynomial, G expressing Newton's gravitational constant, R expressing the radius of the earth, h_nDenotes the load lux number, and g denotes the gravitational acceleration.

In the method for improving the accuracy of the land water reserve abnormality based on the deep learning weight load, the method for obtaining the TWSA of the land water reserve abnormality in the research area by inversion comprises the following steps:

during the time of each study, the least squares problem was minimized and daily land water reserve changes were estimated:

((U_greenx-d)/s)²+β²(L(x))²→min

wherein d represents all corrected crustal load deformation sequences of the grids, s represents the standard deviation of the hydrological load deformation sequences of the grids, beta represents a smoothing factor, x represents that daily land-water reserves to be estimated are abnormal, and L (·) represents a Laplace operator function.

The invention has the following advantages:

the invention constructs a novel Deep learning weighted load inversion method (DWLIM) based on a Long Short-Term Memory network (LSTM), an inverse distance weighting method and a crust load Model, and obtains TWSA in 2011-2020 areas by inversion. DWLIM can be divided into the following five steps: the method comprises the steps that firstly, whether a 1-degree multiplied by 1-degree grid contains a GNSS survey station or not is searched, if the GNSS survey station is contained, the step is skipped to the second step, and if the GNSS survey station is not contained, the step is skipped to the third step; secondly, resolving and preprocessing the GNSS vertical time sequence; thirdly, simulating the vertical deformation of the crust in the unknown grid based on an LSTM regression method; fourthly, utilizing Non-Tidal Atmospheric load (NTAL) and Non-Tidal ocean load (NTOL) model data to correct Atmospheric load and Non-Tidal load on the vertical deformation of the earth crust in all grids; and fifthly, taking the corrected sequences in all grids as input variables, and combining with a crustal load model for inversion to obtain the TWSA in the Chinese area. Therefore, the DWLIM provided by the invention provides important reference for regression of vertical deformation of the crust and monitoring of land water reserves, TWSA can be accurately obtained through inversion, and research results have very important scientific significance and application value in aspects of crust non-structural deformation analysis, hydrodynamic research, sustainable management of water resources in China regions and the like.

Drawings

FIG. 1 is a flowchart illustrating steps of a method for improving the accuracy of land-based water reserves based on deep learning weighted loads according to an embodiment of the present invention;

FIG. 2 is a schematic time span diagram of a continuous GNSS survey station in China and in the surrounding areas according to an embodiment of the present invention;

FIG. 3 is a diagram illustrating the verification effect of simulating vertical deformation of the earth's crust according to an embodiment of the present invention; wherein, 3(a) is a Taylor diagram simulating deformation of the crust of the earth, and 3(b) is a Chinese region average sequence effect diagram;

fig. 4 is a schematic diagram illustrating a result of normalization and decomposition processing of input data of an unknown mesh according to an embodiment of the present invention; wherein, 4(a) is a temperature sequence decomposition result graph after G456 grid normalization processing, and 4(b) is an air pressure sequence decomposition result graph after G456 grid normalization processing;

FIG. 5 is a diagram illustrating simulation results of an unknown mesh crustal deformation sequence in an embodiment of the present invention; wherein, 5(a) is a G464 grid simulation result, 5(b) is a G740 grid simulation result, and 5(c) is a G456 grid simulation result;

FIG. 6 is a graphical representation of a pre-and post-correction sequence of effects corresponding to an atmospheric load correction and a non-marine tidal correction in an embodiment of the present invention;

FIG. 7 is a time series effect diagram of TWSA of a Chinese area based on DWLIM inversion in an embodiment of the present invention;

FIG. 8 is a timing analysis diagram of the TWSA inversion results according to an embodiment of the present invention; wherein, 8(a) is a wavelet analysis graph of DWLIM and TWSA inverted by traditional GNSS, 8(b) is a wavelet analysis graph of DWLIM and GLDAS, 8(c) is a wavelet analysis graph of DWLIM and GRACE, and 8(d) is a sequence graph of TWSA calculation results;

FIG. 9 is a DWLIM inversion result evaluation heatmap according to an embodiment of the present invention; wherein, 9(a) is a PCC index in a GRACE period, 9(b) is a NSE index in the GRACE period, 9(c) is an RMSE index in the GRACE period, 9(d) is a PCC index in the GRACE-FO period, 9(e) is a NSE index in the GRACE-FO period, and 9(f) is an RMSE index in the GRACE-FO period;

FIG. 10 is a graph comparing rainfall and TWSA results in ten large basins in China according to an embodiment of the present invention; wherein, 10(a) is Yangtze river, 10(b) is southeast river basin, 10(c) is sea river basin, 10(d) is Huaihe river basin, 10(e) is yellow river basin, 10(f) is Liaojiang river basin, 10(g) is Songhua river basin, 10 (h) is northwest river basin, 10(i) is southwest river basin, and 10(j) is Zhujiang river basin.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the embodiments of the present invention will be described in detail with reference to the accompanying drawings.

One of the core ideas of the invention is that: the crust can be regarded as a spherical shell which can generate elastic deformation, when land water is migrated in a large scale, the pressure of a huge water mass on the crust can be changed, and the crust can generate nonlinear response of lifting or settlement correspondingly due to the change of the surface pressure. Meanwhile, the GNSS observation station can accurately observe the deformation sequence of the earth surface settlement, the observation precision can reach millimeter level, and in recent years, many scholars accurately invert local TWSA in the area where the GNSS observation station is densely distributed by using the earth crust load deformation model. However, due to the influence of geographical station building conditions, the GNSS stations are distributed extremely unevenly in the global scope, and in an area where the GNSS stations are distributed sparsely, the TWSA of the area cannot be accurately inverted by using the GNSS deformation sequence. Therefore, how to accurately simulate the earth surface subsidence sequence of the area lacking the GNSS survey station is the key point for improving the accuracy of the GNSS inversion TWSA. As shown in fig. 1, the invention combines an LSTM, an inverse distance weighting method and a crust load model to construct a novel deep learning weight load inversion method, that is, a method for improving the accuracy of land-water reserve anomaly based on deep learning weight load, which comprises the following specific processes:

Step 1, dividing a research area into grids of 1 degree multiplied by 1 degree, and searching whether each grid contains a GNSS observation station.

In this embodiment, if the grid includes the GNSS stations, the grid is defined as the observed grid, and step 2 is performed; if no GNSS stations are contained within the grid, the grid is defined as an unobserved grid and step 3 is performed.

And 2, resolving the coordinates every day on the basis of observation data of the GNSS observation station to obtain a crustal vertical deformation sequence determined based on the GNSS observation station, and removing sequence abnormal values and step items larger than three times of standard deviation to obtain a crustal deformation sequence in the observed grid.

In this embodiment, as described above, the GNSS stations can accurately observe the deformation sequence of the earth surface subsidence, the observation accuracy can reach millimeter level, and in recent years, many scholars accurately invert local TWSA in the area where the GNSS stations are densely distributed by using the earth-crust load deformation model. That is, for the observed grid, the corresponding deformation sequence of the earth crust may be obtained in any appropriate manner, for example, based on the observation file, the precise ephemeris file, the navigation file, and the table file observed by the GNSS measurement station, the GAMIT/GLOBK software is used to solve the coordinates of each day to obtain the vertical deformation sequence of the earth crust observed by the GNSS, and the sequence abnormal value and the step item greater than three times of the standard deviation are removed, which is not described herein again.

Step 3, obtaining the earth surface Temperature Sequence and the air pressure Sequence of the unobserved grids, and obtaining the earth surface Temperature Sequence (S) of the unobserved grids_T) Normalizing with an atmospheric pressure Sequence (SAP) to obtain a normalized result; decomposing the normalization result by using a MEEMD decomposition method to obtain m characteristic sequences; GNSS vertical deformation with unobserved gridsThe sequence is a target sequence, the m characteristic sequences are input sequences, and the crustal deformation sequence in the unobserved grid is obtained through regression by using an LSTM regression model.

In this embodiment, the LSTM regression model Hochreter et al proposes an improved model of Recurrent Neural Network (RNN). The LSTM regression model is constructed by a memory storage unit and is trained by a time back propagation algorithm. The model can solve the problems of gradient disappearance and no long-term dependence of RNN. The standard LSTM regression model may specifically include: input door i_tForgetting door f_tAnd an output gate o_tMemory cell c_tAnd the input sequence corresponding to each step length is as follows: x is the number of₁,x₂,...,x_tAnd t represents a step size; memory cell c_tThrough an input gate i_tForgetting door f _tAnd an output gate o_tAnd controlling the memory and forgetting of data.

Further, a memory cell c_tThrough different gates (input gate i)_tForgetting door f_tAnd an output gate o_t) The memory and forgetting of the control data are specifically as follows:

f_t＝σ(W_fx_t+U_fh_t)

o_t＝σ(W_ox_t+U_oh_t)

h_t＝o_t·tanh(c_t)

where σ represents a sigma function; w_f、W_oAnd W_cRespectively representing the weight matrixes in the input process corresponding to the forgetting gate, the output gate and the memory unit; u shape_f、U_oAnd U_cTo representRespectively representing state transition weight matrixes corresponding to the forgetting gate, the output gate and the memory unit, wherein the state transition weight matrixes are S-shaped functions; h is_tA hidden state vector representing an output;

representing the updated memory cell; tanh represents a hyperbolic tangent function. Three kinds of gates (input gate i)_tForgetting door f_tAnd an output gate o_t) Controlling information to enter and leave memory cells together, inputting new information for regulating the entering of the memory cells; the forgetting gate controls the amount of information stored in the memory cell; the output gate defines how much information can be output. The gate structure of the LSTM model enables information on a time series to form a balanced long-term and short-term dependence so as to achieve the purpose of multiple regression.

Further, to solve the problem of few feature sequences in this embodiment, an improved Ensemble Empirical Mode Decomposition (MEEMD) is used to decompose the original features so as to achieve the purpose of increasing the dimension of the feature sequences, and a Decomposition formula specifically adopted by the MEEMD Decomposition method is as follows:

F＝IMF₁+IMF₂+…+IMF_m+noi_w

Wherein F represents an original feature sequence, IMF₁～IMF_mRepresenting n signature sequences obtained by decomposing the original signature sequence F, noi_wWhich represents the white gaussian noise added to the sequence to be decomposed during the decomposition of the original signature sequence F.

Further, the earth crust deformation sequence amplitudes closer to the GNSS stations are more similar due to the similar characteristics of the earth crust deformation in a small area. The inverse distance averaging method is characterized in that the point to be solved is used as the center, the position coordinates of the rest training stations are determined, the distance between each GNSS station and the point to be solved is calculated, and the weight is distributed according to the reciprocal of the distance. Therefore, the regression sequence of deformation of the earth crust in the unobserved lattice is denoted as D_gExpressed as follows:

wherein d is_jRepresents the distance between the jth unobserved grid center and the GNSS stations, n represents the number of GNSS stations,

representing the sum of the distances, Net, of the unobserved grid from each GNSS station_LSTMThe LSTM regression model is represented. Therefore, the simulated crustal vertical deformation of each grid is regressed for n times, and the regression results of the n times are weighted according to the inverse distance weight.

And 4, carrying out atmospheric and non-tidal ocean load correction on all the crustal deformation sequences in the observed lattices obtained in the step 2 and all the crustal deformation sequences in the unobserved lattices obtained in the step 3 by using an NTAL model and an NTOL model to obtain all the crustal load deformation sequences after the grids are corrected.

And 5, taking all the corrected crustal load deformation sequences of the grids obtained in the step 4 as input data, combining the Green function and the crustal load model, performing inversion to obtain TWSA (terrestrial water reserves) with abnormal land water reserves in the research area, and outputting the TWSA.

In this embodiment, the earth crust can be regarded as an elastic surface, and when the quality of the earth surface changes, the earth surface can generate elastic response of earth surface subsidence or rebound, and the deformation is also called earth crust load deformation. However, the deformation of the earth crust load is not only reflected in the vertical direction, but also reflected in the horizontal direction, and the deformation of the earth crust load is more sensitive in the vertical direction, and the deformation amplitude is about 2-3 times of that of the horizontal direction. Green function U_greenThe relationship between the crustal load and the deformation of the crustal load can be indicated as follows:

where θ represents the angular radius from the center of the disk, P_nExpressing Legendre polynomial, G expressing Newton's gravitational constant, R expressing the radius of the earth, h_nDenotes the load lux number, and g denotes the gravitational acceleration.

Furthermore, DWLIM uses the hydrologic load crust deformation corrected by all grids as input data and combines the crust load model for inversion to obtain TWSA in the research area. And carrying out regularization processing on the obtained solution by using a curvature smoothing algorithm, and adding the regularization processing serving as a constraint condition into a solving matrix. In other words, during each study period, the least squares problem is minimized to estimate daily land water reserve changes:

((U_greenx-d)/s)²+β²(L(x))²→min

Wherein d represents all corrected crustal load deformation sequences of the grids, s represents the standard deviation of the hydrological load deformation sequences of the grids, beta represents a smoothing factor, x represents that daily land-water reserves to be estimated are abnormal, and L (·) represents a Laplace operator function.

On the basis of the above-mentioned embodiments, a complete example flow is described below.

Verification of novel deep learning weight load inversion method

1.1, data and model

1.1.1 GNSS data

In this embodiment, a Continuous Operating Reference Station (CORS) vertical observation sequence of a chinese area provided by a chinese continental environment structure monitoring network is used as a data base, and because the establishment time of each CORS measurement station is different, the time span of each measurement station is also different, and the time span situation is as shown in fig. 2. Therefore, in order to ensure the integrity of the vertical deformation sequence, the time span of the embodiment is selected from 2011 to 2020. The total number of original CORS stations in China and the periphery is 268, 5 stations with large time span difference are removed, and the total number of applied sequences is 263. The calculation of the GNSS observation sequence is based on the observation file, the navigation file, the precise ephemeris file and the table file, and the GAMIT/GLOBK 10.4 software is used for calculating to obtain a daily coordinate solution file. Secondly, the GNSS vertical time series is preprocessed, which includes removing observation outliers greater than three times the error and removing sequence step terms due to earthquakes or antenna changes.

1.1.2 GRACE Mascon data

GRACE is a joint Space mission between the National aviation Administration (NASA) and the German aviation Space center (DLR), which launches two GRACE satellites 3 months in 2002 and announces a stop in 2017. The main task of the GRACE satellite is to monitor the temporal and spatial variation of the earth's gravitational field globally, which is mainly caused by earthquakes, glacier equilibrium adjustment, marine and hydrological variations. In chinese context, the gravity changes of GRACE inversion are generally attributed to large-scale hydrological changes.

To verify the reliability of the DWLIM inversion results, the present embodiment compares the result of the gram Mascon with the DWLIM inversion results. However, there is a large uncertainty in the single GRACE-M result due to the different mechanism solution strategies. Therefore, in this embodiment, we use the 2011-2020 Grace Mascon (GRACE-M) data products provided by the American university of Texas Space Center (CSR) and the United states Space agency Jet Propulsion Laboratory (JPL), extract the TWSA of the Chinese area according to the boundary file, and use the average of the two products as the final result Δ TWSA of the GRACE-M _GRACE-M：

Wherein, Delta Mascon_CSRThe TWSA indicates that chinese regions are extracted based on CSR, and the TWSA indicates that chinese regions are extracted based on JPL.

1.1.3 GLDAS data product

For DWLIM construction and verification, the present embodiment uses the temperature variable in the Global Land Assimilation System (Global Land Data Assimilation System, GLADAS) V2.2 Global Land model provided by the united states space and flight administration as the input Data of LSTM regression, and uses the TWSA variable in the GLDAS V2.2 model as the verification Data of DWLIM inversion result. Among them, GLDAS V2.2 is the evolution of early product Land Surface hydrological Model (CLSM), 24 variables including temperature and Land water reserve anomalies. The time span of the GLDAS V2.2 hydrological model is from 2 months to date in 2003, with spatial resolution of 0.25 ° × 0.25 °, temporal resolution of daily scale, and spatial coverage latitude range of: 60 DEG S to 90 DEG N, and the longitude range is as follows: 180 degrees W-180 degrees E.

1.1.4 ERA data products

This example utilizes the surface pressure variables in the ERA5 data product provided by the European center for Medium-Range Weather means (ECMWF) as input data for LSTM regression. Where ERA5 provides day-scale surface pressure data with spatial resolution of 0.1 ° × 0.1 °, time span of 2000 to date, and spatial coverage worldwide. The specific data details of the ERA dataset can be known from the website (https:// www.ecmwf.int/en/forecasets/dataset/ecmwf-reanalysis-v 5).

1.2 evaluation index

In this embodiment, the accuracy of the DWLIM inversion result is evaluated by using three evaluation indexes, namely Root Mean Square Error (RMSE), Pearson's Correlation Coefficient (PCC), and Nash-Sutcliffe efficiency Coefficient (NSE), and the calculation formula is as follows:

wherein Y represents the true sequence, X represents the inversion result,

and

respectively, the average values of Y and X, and k represents the number of discrete points of the sequence. The RMSE can be used to estimate the dispersion deviation of the inversion result from the true value, with smaller RMSE values indicating more accurate inversion results. The NSE is mainly used for evaluating the quality of the hydrological model, and the value of the NSE is less than or equal to 1; the larger the value, the better the hydrological model; when NSE is close to 0, it indicates that the hydrographic model effect is close to the average of the observed values. PCC is mainly used for describing a linear correlation relationship between two sequences, and the PCC value is between-1 and 1; the closer the value is to 1, the more reliable the inversion result is.

1.3, DWLIM-based TWSA calculation of Chinese area

1.3.1 verification of LSTM simulation unknown grid crustal deformation

In order to verify the accuracy of DWLIM simulation grid crust deformation, stations with correlations between the pressure and temperature sequences in the grid range and GNSS sequences larger than 0.5 are selected, and the obtained 75 GNSS station sequences are used as control sequences. And 263 sequences are used for verification, and when a control station exists in the grid network where the GNSS station is located, the verification returns 74 times; and when no GNSS station exists in the grid, the grid regresses for 75 times, and the GNSS sequence is used as a true value to be compared with the simulated crustal deformation sequence. Firstly, a crustal influence factor sequence of all grids in a research area is extracted, wherein the crustal influence factor sequence comprises two characteristic variables of air pressure and temperature. Then, the two characteristic sequences are respectively decomposed into 10 characteristic modal components by using a MEEMD decomposition method so as to achieve the purpose of increasing dimensionality. Secondly, regression is carried out by using the modal component obtained by decomposition as input data and the GNSS vertical sequence as output data by using an LSTM method, and weighting is carried out on the 75 simulation results according to the reciprocal of the distance. And finally, obtaining a simulated crustal deformation sequence of the grid where each GNSS observation station is located. The simulation results were evaluated according to the RMSE and PCC indices, and the evaluation results are shown in fig. 3.

As can be seen from fig. 3, the RMSE values of most of the station sequences are within 5mm, however, due to the difference in the observed quality between GNSS sequences, the RMSE values of some stations are larger. Further statistics show that 68.63% of the comparison results are within 6 mm. And then evaluating the consistency relation between the simulated sequence and the actual value by utilizing the PCC index, wherein the maximum PCC value can reach 0.87, and the average value is 0.53. Fig. 3(b) shows an effect diagram of the chinese region simulation and the real sequence averaging, and it can be seen that the periodicity of the sequence can be effectively simulated through the third step in DWLIM, and the sequence obtained based on the LSTM simulation is smoother and the sequence characteristics are more obvious.

1.3.2 acquisition of hydrologic load deformation

(3.2.1) simulation of deformation of crust within unknown grid

In this embodiment, temperature and air pressure are selected as input data of the LSTM regression method, and in order to fully utilize the sequence characteristics of the input data, the present embodiment utilizes the decomposition method of MEEMD to decompose the surface temperature and air pressure sequence into IMFs₁～IMF ₁₀10 modal components. Taking the G456 grid as an example, the decomposition effect of the input sequence is shown in fig. 4.

The IMF1 components in FIGS. 4(a) and 4(b) are the normalized raw surface temperature and barometric pressure sequences, IMFs ₁～IMF₁₀The obtained characteristic items from high frequency to low frequency are decomposed. The decomposed result can better show a trend term, a seasonal term and a residual term of the sequence, wherein the sequence of the temperature and the air pressure shows the cycle characteristics of the whole year. In the LSTM regression method, 10 IMF components obtained by decomposition as shown in fig. 4 are used as an input sequence, 75 control GNSS vertical deformation sequences are used as output data, and at the same time, inverse distance weighting is performed according to the distance between the center of the grid and 75 stations, and finally, the crustal vertical deformation sequence of the unknown grid is simulated. The results of the simulation of the unknown grids are shown in fig. 5(a) -5 (c), taking the G464, G740, and G456 grids as examples.

The China continental GNSS survey stations are not uniformly distributed, and the number of the survey stations is small, so that the whole coverage of the China continental crust cannot be achieved. Therefore, the simulation of the vertical deformation of the earth crust of the unknown grid is very necessary. In which, fig. 5(a) -5 (c) show the results of vertical deformation of the crust obtained by simulation using 20 IMF components obtained by MEEMD decomposition as input data and 75 control GNSS stations as output data, in combination with the inverse distance weighting theory and the LSTM regression method. The simulation result shows that the periodic term and annual amplitude of the vertical deformation of the crust can be well simulated according to the strategy, and a good data basis is provided for the inversion of TWSA.

(3.2.2) correcting all mesh deformation sequences

Load deformation caused by atmospheric and non-ocean tidal loads exists in a GNSS vertical deformation sequence, so in order to extract the crustal deformation caused by hydrological loads, the invention uses NTAL and NTOL models provided by a German Boltzman earth system modeling group (http:// esmdata. gfz-potsdam. de:8080/repository) as correction data, adds two correction sequences into a crustal load deformation sequence in a Chinese area, and comprises a GNSS vertical time sequence and a DWLIM method simulation unknown grid load deformation sequence, wherein the vertical deformation sequences before and after correction are shown in figure 6.

As further shown in FIG. 6, the mean values of the load deformation amplitudes of NATL and NOTL are 3.62mm and 0.22mm, respectively. Annual amplitude convex areas which are not corrected by ocean atmosphere are mainly positioned in the north and east of China, and the maximum amplitude convex area can reach 5.5 mm; instead of ocean tide correction, annual amplitude is small, the maximum amplitude is only 1.5mm, and the maximum amplitude is mainly distributed in coastal areas of east China. As can be seen from FIG. 6, the corrected sequence amplitude and phase have small changes, and a more accurate hydrologic load deformation sequence is provided for DWLIM inversion TWSA.

1.3.3 inversion of China regional TWSA based on DWLIM

In this embodiment, a gray function matrix of point loads in the china region is first calculated, 0.25 ° × 0.25 ° is selected as the spatial resolution of inversion, where the disc radius corresponding to the spatial resolution of 0.25 ° is 13.90km, and preset parameters whose expansion boundary range is 2 ° and β is 0.01 are selected. Taking all corrected crustal vertical deformation sequences of 1-degree multiplied by 1-degree grids in a Chinese area as input data, and combining the following formula for inversion to obtain daily scale land water reserve change values in a Chinese range:

((U_greenx-d)/s)²+β²(L(x))²→min

in order to better verify the accuracy of the DWLIM inversion result, the invention carries out first-order term correction processing on the inversion result. The calculated time series of the average TWSA in chinese area is shown in fig. 7.

As shown in fig. 7, the DWLIM method can effectively reverse the annual characteristics and amplitude of the TWSA sequence of the chinese region. The method can be used for calculating the protruding positions of annual amplitude, which are positioned in Yunnan province, south of autonomous region in Tibet and south of North China plain, and the result is consistent with the inversion result of predecessor. According to the calculated annual phase relation of the inversion result, the annual phase in the Qinghai-Tibet autonomous region range is obviously lower than other provinces in the Chinese region. To verify the accuracy of the DWLIM inversion result, this embodiment compares the inversion result with the graph-M result, the GLDAS corpus result, and the TWSA result obtained by the conventional GNSS inversion.

1.4 verification of novel deep learning weight load inversion method

1.4.1 spatial feature validation of TWSA results

In order to compare the accuracy of DWLIM in inverting TWSA in chinese region, the present embodiment utilizes the conventional GNSS inversion TWSA method (TRA) respectively_GNSS) The Chinese area TWSA is extracted from the GRACE Mascon processing result and the GLDAS hydrological result, and the annual amplitude of each result is calculated.

Based on DWLIM strategy, the method can effectively invert the areas with annual amplitude bulges in Chinese areas, such as southwest region, southeast region, Qinghai-Tibet autonomous region and the like in Yunnan province, and the space expression result is consistent with GRACE and GLDAS results. And the DWLIM inversion result is obviously raised in plain areas in North China, and the result is stronger than GRACE and GLDAS. The reason for this may be that the deformation of the crust in the area is obvious, and the hydrologic load deformation sequence cannot be extracted more accurately only by correcting the atmospheric load and the non-ocean tides. Meanwhile, amplitude convex areas of the northern part of Xinjiang and the northern part of Heilongjiang are obtained by inversion based on a DWLIM method. However, the annual amplitude spatial distribution obtained by the traditional method for inverting the TWSA by GNSS has a speckle characteristic, and the reason is that the distance of the disk radius is limited, and the TWSA condition in a limited range around the GNSS survey station can only be obtained by inversion. And the traditional GNSS inversion TWSA method causes signal loss in a sparse area of the GNSS survey station through smoothing. In summary, the limitation of the disk radius on the TWSA of the GNSS inversion can be weakened through the simulation of the deformation of the earth crust of the unknown grid.

1.4.2 timing characteristics verification of results

In order to verify the time sequence reliability of DWLIM inversion TWSA, the present embodiment compares the DWLIM inversion result with the TWSA inversion result of the traditional GNSS, method, GRACE, and GLDAS results, respectively. In order to further analyze the sequence phase relationship between the DWLIM inversion result and the conventional GNSS inversion method, GLDAS and coarse, the present embodiment plots cross-wavelet analysis graphs of DWLIM-Traditional, DWLIM-GLDAS and DWLIM-coarse, respectively, as shown in fig. 8(a) to 8 (c). Meanwhile, fig. 8(d) shows a timing chart of DWLIM, conventional GNSS, GRACE, and GLDAS results for extracting the TWSA of the chinese area.

As can be seen from fig. 8(a) to 8(c), the TWSA results of DWLIM inversion performed well in phase with the conventional GNSS inversion results, GLDAS and GRACE. The DWLIM method can effectively invert the anniversary and anniversary amplitudes of TWSA sequences, and is consistent with GRACE and GLDAS calculations (fig. 8 (d)). However, due to differences in observation, the yearly amplitude of DWLIM is slightly larger than other methods. The reason for this is that the corrected deformation sequence of the crust has deformation caused by other loads, and the deformation sequence of other loads is difficult to change, so that a single hydrological load deformation sequence cannot be separated. And the blank period between GRACE and GRACE-FO (FIG. 8(d) gray shaded portion) can be filled in according to the TWSA result obtained by DWLIM inversion. The seasonality of DWLIM inversion results is more obvious than that of TWSA results obtained by traditional GNSS inversion. In order to quantify the advantages of DWLIM over the TWSA inversion of the conventional GNSS, the present embodiment evaluates the inversion results by using PCC, NSE and RMSE, and the evaluation results are shown in fig. 9.

As can be seen from fig. 9, the TWSA results obtained by DWLIM inversion are significantly better than those obtained by conventional GNSS inversion. In the GRACE period (2011-2017), the maximum indexes of PCC, NSE and RMSE of DWLIM inversion results reach 0.81, 0.62 and 2.18cm respectively; in the GRACE-FO period (2018-2020), PCC, NSE and RMSE of DWLIM inversion results reach 0.71 cm, 0.49 cm and 2.4cm respectively. The results show that the consistency of the inversion result of DWLIM and the GLDAS result is better, because the resolution of GRACE is month scale, and the coarse time resolution can cause signal loss. For further statistics of data, compared with the traditional GNSS inversion method, the DWLIM inversion TWSA results are respectively improved by 67.11%, 128.15% and 22.75% on the PCC, NSE and RMSE indexes on average. The result further shows that land water reserve change of a GNSS observation station sparse area can be effectively obtained through inversion based on the DWLIM method, and the space-time characteristic of the TWSA result is superior to that of a TWSA result obtained through inversion of a traditional GNSS.

2. Application of novel deep learning weight load inversion method

2.1, analyzing the time-space change of TWSA of each basin in China by combining rainfall data

The embodiment has verified that DWLIM can effectively invert TWSA in Chinese areas, and can detect the raised point area with abnormal annual amplitude of land water reserves. The vertical load deformation of the crust is mainly influenced by the land water reserves, and when the land water reserves are increased, the crust tends to decline; conversely, as the land water reserve load decreases, the crust exhibits a tendency to rebound upward. In order to study the cause of the land water reserve in the chinese region, the present embodiment combines the monthly rainfall products provided by the chinese Meteorological data network (CMA) and extracts the rainfall values and TWSA in the ten major watersheds of China, where the TWSA is obtained by the DWLIM method, and the comparison effect graph is shown in fig. 10.

As can be seen from fig. 10, the annual amplitude of TWSA is positively correlated with the annual amplitude of rainfall events overall. The daily mean value of rainfall in the songhua river basin (fig. 10(g)) and the Liaohe river basin (fig. 10(f)) is significantly higher than that in other basins, and the amplitude of the corresponding TWSA is also significantly higher than that in other basins. In terms of phase, the phase relation between the land water reserves and the rainfall sequences in the Chinese area is good, and when the rainfall histogram peaks, the TWSA time sequence is also located at the peak. Further shows the reliability of DWLIM inversion of the TWSA in the Chinese area on the phase. Meanwhile, it can be seen that high-frequency noise exists in the time sequence of the TWSA inversion result, which also affects the inversion or prediction of TWSA, because a great amount of system noise brought by ionosphere, troposphere, clock error or multipath effect exists due to the inherent properties and observation means of the receiver during the process of observing deformation of the crust of the GNSS. Therefore, in subsequent studies, the studies of noise classification and removal of GNSS vertical sequences will be focused on, providing a cleaner sequence for the inversion of TWSA.

In summary, the invention discloses a method for improving accuracy of land-Water reservoir Anomaly based on deep learning weight load, which is based on a Global Navigation Satellite System (GNSS) continuous observation station of dense distribution and can effectively invert land-Water reservoir Anomaly (TWSA). However, GNSS continuous observers are highly unevenly distributed around the world under the influence of natural environmental conditions, which greatly limits the research related to the inversion of TWSA using GNSS vertical deformation. The method simulates the vertical deformation of the grid which does not contain the GNSS survey station so as to achieve the purpose of improving the accuracy of the GNSS inversion TWSA. The result shows that DWLIM inversion result and maximum Pearson correlation coefficient PCC, Nash efficiency coefficient NSE and root mean square error RMSE of GRACE and GLDAS respectively reach 0.81, 0.61 and 2.18cm, and compared with the traditional GNSS inversion method, the DWLIM inversion method averagely improves indexes of PCC, NSE and RMSE by 67.11%, 128.15% and 22.75%. The result shows that in the area where the GNSS survey stations are unevenly distributed, the accuracy of TWSA inversion can be effectively improved based on the DWLIM inversion strategy. The TWSA of the ten river basins in China is obtained by DWLIM inversion, the change condition of the TWSA of the ten river basins is analyzed by combining rainfall data, good consistency exists on the phase position of the DWLIM inversion result and the rainfall data, and wave crests are located in six months and seven months every year; secondly, the amplitude of the TWSA sequence of each watershed is in positive correlation with the amplitude of the rainfall sequence, and the reliability of TWSA inversion based on DWLIM is further verified.

Although the present invention has been described with reference to the preferred embodiments, it is not intended to limit the present invention, and those skilled in the art can make variations and modifications of the present invention without departing from the spirit and scope of the present invention by using the methods and technical contents disclosed above.

Those skilled in the art will appreciate that the invention may be practiced without these specific details.

Claims (8)

1. A method for improving the accuracy of land water reserve abnormity based on deep learning weight load is characterized by comprising the following steps:

step 1, dividing a research area into grids of 1 degree multiplied by 1 degree, searching whether each grid contains GNSS stations, if the grids contain the GNSS stations, defining the grids as observed grids, and executing step 2; if the grid does not contain the GNSS observation station, defining the grid as an unobserved grid, and executing the step 3;

step 2, resolving the coordinates every day based on observation data of the GNSS observation station to obtain a crustal vertical deformation sequence determined based on the GNSS observation station, and removing sequence abnormal values and step items larger than three times of standard deviation to obtain a crustal deformation sequence in the observed grid;

Step 3, acquiring a surface temperature sequence and an air pressure sequence of the unobserved grids, and performing normalization processing on the acquired surface temperature sequence and the acquired air pressure sequence of the unobserved grids to obtain a normalization result; decomposing the normalization result by using a MEEMD decomposition method to obtain m characteristic sequences; taking the GNSS vertical deformation sequence of the unobserved grid as a target sequence, taking the m characteristic sequences as input sequences, and obtaining a crustal deformation sequence in the unobserved grid by regression by utilizing an LSTM regression model;

step 4, utilizing an NTAL model and an NTOL model to correct atmospheric and non-tidal ocean loads of all crustal deformation sequences in the observed lattices obtained in the step 2 and all crustal deformation sequences in the unobserved lattices obtained in the step 3 to obtain crustal load deformation sequences after all the lattices are corrected;

and 5, taking all the corrected crustal load deformation sequences of the grids obtained in the step 4 as input data, combining the Green function and the crustal load model, performing inversion to obtain TWSA (terrestrial water reserves) with abnormal land water reserves in the research area, and outputting the TWSA.

2. The method for improving the accuracy of land-water reservoir anomaly based on deep learning weight load according to claim 1, wherein the LSTM regression model comprises: input gate i _tDoor f for forgetting to leave_tOutput gate o_tMemory cell c_tAnd the input sequence corresponding to each step length is as follows: x is the number of₁,x₂,...,x_tAnd t represents a step size; memory cell c_tThrough an input gate i_tForgetting door f_tAnd an output gate o_tAnd controlling the memory and forgetting of data.

3. The method for improving the accuracy of land-water reserves based on deep learning weight load as claimed in claim 2, wherein the memory unit c_tAnd input gate i_tForgetting door f_tAnd an output gate o_tThe relationship between them is expressed as follows:

f_t＝σ(W_fx_t+U_fh_t)

o_t＝σ(W_ox_t+U_oh_t)

h_t＝o_t·tanh(c_t)

where σ represents a sigma function; w_f、W_oAnd W_cRespectively representing forgetting gate, output gate and memory unitA corresponding weight matrix in the input process; u shape_f、U_oAnd U_cThe expression respectively represents the state transition weight matrixes corresponding to the forgetting gate, the output gate and the memory unit, and the state transition weight matrixes are S-shaped functions; h is_tA hidden state vector representing an output;

representing the updated memory cell; tanh represents a hyperbolic tangent function.

4. The method for improving the accuracy of the anomaly of the land water reserves based on the deep learning weight load according to claim 3, wherein the MEEMD decomposition method adopts the following decomposition formula:

F＝IMF₁+IMF₂+…+IMF_m+noi_w

wherein F represents the original characteristic sequence, IMF₁～IMF_mRepresenting n signature sequences obtained by decomposing the original signature sequence F, noi _wWhich represents the white gaussian noise added to the sequence to be decomposed during the decomposition of the original signature sequence F.

5. The method for improving the accuracy of the anomaly of the land water reserves based on the deep learning weight load as claimed in claim 4, wherein the regression-derived sequence of deformation of the crust in the unobserved grid is denoted as D_gExpressed as follows:

wherein, d_jRepresents the distance between the jth unobserved grid center and the GNSS stations, n represents the number of GNSS stations,

representing the sum of the distances, Net, of the unobserved grid from each GNSS station_LSTMThe LSTM regression model is represented.

6. The method for improving the accuracy of land water reserve abnormality based on deep learning weighted load according to claim 5, wherein a Green's function is used to indicate the relationship between the crustal load and the deformation of the crustal load.

7. The method for improving the accuracy of land water reserves abnormality based on deep learning weight load as claimed in claim 6, wherein Green function U_greenIs represented as follows:

where θ represents the angular radius from the center of the disk, P_nExpressing Legendre polynomial, G expressing Newton's gravitational constant, R expressing the radius of the earth, h_nDenotes the load lux number, and g denotes the gravitational acceleration.

8. The method for improving accuracy of land-water reservoir anomaly based on deep learning weight load according to claim 7, wherein the step of obtaining the TWSA of the land-water reservoir anomaly in the research area through inversion comprises the following steps:

during the time of each study, the least squares problem was minimized and daily land water reserve changes were estimated:

((U_greenx-d)/s)²+β²(L(x))²→min

wherein d represents all corrected crustal load deformation sequences of the grids, s represents the standard deviation of the hydrological load deformation sequences of the grids, beta represents a smoothing factor, x represents that daily land-water reserves to be estimated are abnormal, and L (·) represents a Laplace operator function.

Images