CN113792450A

CN113792450A - Method for improving land water reserve inversion accuracy based on machine learning load model

Info

Publication number: CN113792450A
Application number: CN202110939172.5A
Authority: CN
Inventors: 郑伟; 尹文杰; 沈祎凡
Original assignee: China Academy of Space Technology CAST
Current assignee: China Academy of Space Technology CAST
Priority date: 2021-08-16
Filing date: 2021-08-16
Publication date: 2021-12-14

Abstract

The invention discloses a method for improving inversion accuracy of land water reserves based on a machine learning load model, which comprises the following steps: dividing a research area according to a grid of 1 degree multiplied by 1 degree, and judging whether the grid contains a GPS station or not; when the grid is determined to contain the GPS sites, preprocessing the vertical time sequence of the GPS sites to obtain a real GPS vertical deformation sequence; when the situation that the grid does not contain the GPS station is determined, a simulated crustal vertical deformation sequence is obtained based on a random forest method; correcting atmospheric and non-tidal ocean loads of the real GPS vertical deformation sequence or the simulated crustal vertical deformation sequence to obtain a corrected crustal vertical deformation sequence; and taking the corrected vertical deformation sequence of the crust as input data of a crust load model, and performing inversion on the land water reserves TWSA. The method can realize the simulation of the deformation of the earth crust load in an unobserved area without GPS station distribution, realize the inversion of the TWSA of the land water reserves and improve the inversion precision of the TWSA.

Description

Method for improving land water reserve inversion accuracy based on machine learning load model

Technical Field

The invention belongs to the technical field of intersection of satellite gravimetry, hydrology and the like, and particularly relates to a method for improving land water reserve inversion accuracy based on a machine learning load model.

Background

Land water is an important resource for the continued development of industry, agriculture and human life, yet represents only 3.5% of the global water resources. In China, uneven spatial and temporal distribution of land water reserves causes a series of problems for the survival and development of people. Particularly, in the southwest of China, contradiction between water resource supply and demand causes a series of natural disasters (such as drought, flooding, water and soil loss and the like). Therefore, land impoundment anomaly (TWSA) conditions must be monitored to assess their potential and long-term sustainability. Redistribution of water mass changes the gravitational field in a region that can be monitored by gravity recovery and climatic test (Grace) satellites.

The GRACE satellite is a part of a scientific exploration project of a terrestrial system hosted by the national aerospace agency in 2002, and provides an effective measurement means for large-scale TWSA inversion. Studies have shown that GRACE satellites can monitor TWSA trends and seasonal characteristics with unprecedented accuracy. However, it is difficult to reverse TWSA for small scale regions due to poor spatial resolution of the GRACE satellite. Furthermore, as the GRACE satellite system ages, the GRACE task ends in 2017, with the subsequent task GRACE-FO transmitted in 2018, with a 14-month window period between the two tasks. Therefore, it is important to find an alternative method to fill the gap between GRACE and GRACE-FO.

The redistribution of the earth's crust load causes changes in the earth's surface mass, resulting in complex deformation displacements of the earth's crust in the horizontal (N, E) and vertical (U) directions. And the deformation in the vertical direction is particularly prominent, and the relation can be constructed by utilizing a crustal load model. The deformation response of the earth's crust load can be measured by a variety of observation methods, such as Global Positioning System (GPS), interferometric technology (InSAR), and Very Long Baseline Interferometry (VLBI).

The GPS observation means has the advantages of high efficiency, accuracy, all weather and the like, and can accurately and truly acquire the deformation of the station. Meanwhile, the China Crustal Movement Observation Network (CMONOC) has been established for 18 years, can continuously observe the deformation of the crustal in China and accumulate a large amount of observation data. Many previous studies have analyzed the geophysical phenomena in typical areas of china, such as the south of china, the north china plain, the province of sichuan, etc., using data provided by CMONOC. However, the spatial maldistribution of GPS sites greatly limits its application to analyzing geophysical phenomena. In addition, GPS sites in the southwest region of china are mainly concentrated in the sichuan province and the Yunnan province. Therefore, how to better simulate the deformation of the crust of the unobserved area becomes a hot research problem in recent years. Here, the unobserved area indicates an area not including a GPS site within a grid of 1 ° × 1 °.

In the study of the redistribution of surface mass, it is a common question how to better establish the relationship between surface loads (such as land impoundments, atmospheric and non-tidal ocean loads) and their deformation responses. In addition, GPS provides a measurement means to observe deformation of the crust independently and in near real time. In view of the many advantages of GPS, many researchers have achieved significant effort in recent years to invert geophysical phenomena using GPS observations. In 2004, Heki explored the seasonal weight distribution problem in japan with dense GPS site sequences, which showed that seasonal variations in crustal loads could be effectively inverted with GPS vertical time sequences and could supplement GRACE. Fu et al analyzed Nepal vertical load deformation in 2012 in combination with GPS and GRACE, and the results show that the seasonal deformation sequence is consistent with the inversion results of GRACE. Fu et al in 2015 set forth a new method for inverting ground impoundments using GPS sequences and demonstrated that the inversion results can be used to fill the gap between GRACE and GRACE-FO. Although great efforts have been made to invert crustal loads using GPS, there has been little research on applying the method of inverting the loads to areas where GPS sites are poorly distributed. In 2021, ginger zhongshan et al proposed a new algorithm, which combines a traditional load model with a spherical sleian basis function, and the study shows that a sparse GPS array can also be used as a method for inverting TWSA spatiotemporal changes. However, the above studies have focused on inverting the earth's crust load using GPS sequences, and neglect the problem of how to encrypt the GPS site distribution. Therefore, in areas without GPS station distribution, how to better simulate the deformation of the earth crust load is a key problem.

Disclosure of Invention

The technical problem of the invention is solved: the method aims to realize simulation of crustal load deformation without an unobserved area distributed by GPS sites, further realize inversion of TWSA (terrestrial trunked satellite system) of the land water reserves, and improve the inversion precision of the TWSA.

In order to solve the technical problem, the invention discloses a method for improving inversion accuracy of land water reserves based on a machine learning load model, which comprises the following steps:

dividing a research area according to a grid of 1 degree multiplied by 1 degree, and judging whether the grid contains a GPS station or not;

when the grid is determined to contain the GPS sites, preprocessing the vertical time sequence of the GPS sites to obtain a real GPS vertical deformation sequence;

when it is determined that no GPS station is contained in the grid, simulating vertical deformation sequences of the crust in all grids of a research area based on a random forest method to obtain a simulated vertical deformation sequence of the crust;

correcting atmospheric and non-tidal ocean loads of the real GPS vertical deformation sequence or the simulated crustal vertical deformation sequence to obtain a corrected crustal vertical deformation sequence;

and taking the corrected vertical deformation sequence of the crust as input data of a crust load model, and performing inversion on the land water reserves TWSA.

In the method for improving the inversion accuracy of land water reserves based on the machine learning load model, the vertical time sequence of the GPS station is preprocessed to obtain a real GPS vertical deformation sequence, and the method comprises the following steps:

extracting a vertical time sequence corresponding to a GPS site from a Chinese regional GPS time sequence file provided by a China continental construction network CMONOC;

and preprocessing the extracted vertical time sequence corresponding to the GPS site, and removing a step item and an abnormal value caused by earthquake around the GPSGPS site and receiver abnormity to obtain a real GPS vertical deformation sequence.

In the method for improving the land water reserve inversion accuracy based on the machine learning load model, the preprocessing process of the vertical time sequence corresponding to the GPS station is as follows:

obtaining a first-order difference delta U of a vertical time sequence corresponding to the GPS station through the formula (1)_gps：

ΔU_gps＝U_gps(t)-U_gps(t-1)···(1)

Wherein t represents time, U_gps(t) vertical time series, U, corresponding to GPS station at time t_gps(t-1) represents a vertical time sequence corresponding to the GPS station at the t-1 moment;

when Δ U_gpsWhen the current step is larger than 5mm, the current step is abnormal, and the current step is corrected through the formula (2) to obtain a real GPS vertical deformation sequence U_deleap：

U_deleap＝U_gps(t)+ΔU_gps···(2)

Where, t is M, M +1, …, M indicates the sequence node position where the abnormal step occurs, and M indicates the total number of GPS stations.

In the method for improving inversion accuracy of land water reserves based on the machine learning load model, based on a random forest method, vertical deformation sequences of the crust in all grids of a research area are simulated to obtain a simulated crust vertical deformation sequence, which includes:

the temperature and the air pressure corresponding to the GPS station are used as input data, the vertical time sequence of the GPS station is used as output data, and the random forest method is utilized to sequentially regress for M times to obtain output results of M grids;

and carrying out average processing on output results of the M grids to obtain simulated vertical deformation sequences of the crust in all grids of the research area.

In the method for improving the inversion accuracy of land water reserves based on the machine learning load model, the simulated crustal vertical deformation sequence is expressed as follows:

wherein, U_gridRepresenting a simulated crustal vertical deformation sequence, T representing the total number of constructed random, decorrelated decision trees in a random forest method, g_i(x) And (4) representing the inversion result of each decision tree in the random forest method.

In the method for improving the inversion accuracy of the land water reserves based on the machine learning load model, the atmospheric and non-tidal ocean load correction is carried out on the real GPS vertical deformation sequence under the following flow:

by the formula (4), the real GPS vertical deformation sequence U is processed_deleapCarrying out atmospheric and non-tidal ocean load correction to obtain a corrected vertical crustal deformation sequence U_hy：

U_hy＝U_deleap-U_MERRA-U_OMCT···(4)

Wherein, U_MERRAIndicating deformation of the crust, U, by atmospheric load_OMCTRepresenting the deformation of the crust resulting from non-ocean tidal loads.

In the method for improving the inversion accuracy of the land water reserves based on the machine learning load model, the process of correcting the atmospheric and non-tidal ocean loads of the simulated crustal vertical deformation sequence is as follows:

by the formula (5), the simulated crustal vertical deformation sequence U_gridCarrying out atmospheric and non-tidal ocean load correction to obtain a corrected vertical crustal deformation sequence U_hy：

U_hy＝U_grid-U_MERRA-U_OMCT···(5)

In the method for improving the inversion accuracy of the land water reserves based on the machine learning load model, the obtained corrected vertical deformation sequence of the crust is used as input data of the crust load model to invert the land water reserves TWSA, and the method comprises the following steps:

according to the corrected vertical deformation sequence U of the crust_hyEstimate daily land water reserve change Load_TWSA：

Load_TWSA＝((Gx-U_hy)/σ)²+β²(L(x))²→min···(6)

Wherein the daily land water reserve changes Load_TWSAThe land water reserve TWSA is obtained by inverting the land water reserve TWSA, Gx represents a Green function coefficient matrix, sigma represents the standard deviation of a hydrologic load deformation sequence, beta represents a smoothing factor, and L (x) represents a Laplace operator function.

The invention has the following advantages:

land-water reserves anomalies (TWSA) can be accurately inverted using densely distributed Global Positioning System (GPS) sites. However, the application of GPS-based inversion TWSA is limited to a greater extent when the GPS stations are poorly distributed or unevenly distributed. Based on the problem, the invention discloses a method for improving the inversion accuracy of land water reserves based on a machine learning load model, which can improve the inversion accuracy of TWSA by increasing the space coverage of the vertical deformation of the earth surface.

Firstly, on the basis of traditional GPS inversion TWSA, a random forest regression method is utilized to simulate the deformation of the earth crust load in an unknown grid, and a novel Machine Learning Load Inversion Method (MLLIM) which is beneficial to improving the accuracy of the inversion TWSA is constructed.

Secondly, taking the GPS sequence in the southwest area of China as a data base, respectively utilizing MLLIM and a traditional GPS inversion method to invert TWSA in the southwest area of China, and comparing with gravity recovery and climate experiment (GRACE) satellite and global land assimilation system (GLADAS) model results. The results show that the Pearson Correlation Coefficients (PCC) based on MLLIM inversion results and GRACE-FO results are 0.91 and 0.88, respectively, the coefficient of solution (R)²) 0.76 and 0.65, respectively; PCC and R based on MLLIM inversion result and GLDAS result²0.79 and 0.64 respectively, compared with the traditional GPS inversion method, in PCC and R²The upper average increases by 8.85% and 7.99%, indicating that the accuracy of TWSA inversion can be improved based on MLLIM.

And thirdly, performing reverse transformation on TWSA of each province in the southwest of China by using MLLIM, and analyzing the change condition of the TWSA of each province in the research area by combining rainfall data, wherein the result shows that the TWSA obtained based on the MLLIM is consistent with the spatial distribution of the rainfall data, and the time-space change bulges are positioned in the southwest part of the Yunnan province and the southeast part of the Guangxi province. Compared with GRACE and GLDAS results based on MLLIM inversion results, the results show that PCC reaches 0.86 and 0.94, and TWSA in areas with sparse GPS station distribution can be accurately inverted based on MLLIM.

TWSA based on MLLIM-accurate inversion can be used to fill the gap between GRACE and GRACE-FO. The novel machine learning load inversion method has the advantages of high accuracy and high calculation speed for inverting the land water reserve change.

Drawings

FIG. 1 is a flowchart illustrating steps of a method for improving accuracy of land-water reserve inversion based on a machine learning load model according to an embodiment of the present invention;

FIG. 2 is a graph of vertical deformation versus center distance for disks of different mass and different radii placed on the surface in an embodiment of the present invention;

FIG. 3 is a diagram of site information in the southwest region of China according to an embodiment of the present invention; wherein 3(a) is a GPS site distribution diagram in the southwest region of China; 3(b) is a land utilization type diagram in southwest China; 3(c) is a position relation graph between the unknown grids (G1-G87) and the GPS station; 3(d) is a position diagram of the southwest region of China in China;

FIG. 4 is a schematic diagram of a temperature, air pressure, and vertical deformation sequence of a grid in which a GPS station is located according to an embodiment of the present invention; wherein, 4(a) is an SCPZ site; 4(b) is SCXJ site; 4(c) is an YNMJ site; 4(d) is YNML site;

FIG. 5 is a diagram illustrating a result of inversion of deformation of a GPS grid crust based on MLLIM in an embodiment of the present invention; wherein 5(a) is PCC between the mimic sequence and the true sequence; 5(b) is the RMSE between the mimic sequence and the true sequence; 5(c) is a simulation sequence and real sequence effect graph;

fig. 6 is a schematic diagram of a simulation result of deformation of the crust at an unknown grid based on MLLIM in an embodiment of the present invention; wherein 6(a) is the spatial position information of the prediction grid; 6(b) is a G36 grid simulation sequence effect chart; 6(c) is a G43 grid simulation sequence effect chart;

FIG. 7 is a schematic diagram illustrating an example of an MLLIM-based inversion of TWSA results in southwest area of China; wherein 7(a) is a 0.25 ° disk TWSA result; 7(b) is a yearly amplitude spatial distribution map; 7(c) is a TWSA sequence diagram in the southwest region of China;

FIG. 8 is a graph of TWSA results obtained by inversion based on MLLIM and conventional GPS in an embodiment of the present invention; wherein, 8(a) -8 (c) and 8(d) -8 (f) respectively represent the annual amplitude, semi-annual amplitude and sequence period of TWSA calculated by MLLIM and traditional GPS inversion method;

FIG. 9 is a graph of comparison between MLLIM and GRACE results in an embodiment of the present invention; wherein 9(a) is the average result after fitting; 9(b) is the annual amplitude of the RACE sequence; 9(c) is the half-cycle annual amplitude of the GRACE sequence; 9(d) is the period of the GRAC sequence; 9(e) is a linear scattergram;

FIG. 10 is a graph showing the comparison between the MLLIM and GLDAS results according to the present invention; wherein 10(a) is the average result after fitting; 10(b) is the annual amplitude of the GLDAS sequence; 10(c) is the half-cycle annual amplitude of the GLDAS sequence; 10(d) is the period of the GLDAS sequence; 10(e) is a linear scattergram;

FIG. 11 is a TWSA result graph obtained for each province based on MLLIM and analyzed in southwest area of China with rainfall data in the embodiment of the present invention; wherein 11(a) to 11(e) are rainfall space distribution diagrams of each province; 11(f) -11 (j) are graphs of TWSA change condition combined with rainfall analysis.

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: based on machine learning and a crust load model, a novel Machine Learning Load Inversion Method (MLLIM) is constructed and is inverted to obtain TWSA in 2011-plus 2019 in southwest China. The main flow of MLLIM is as follows: firstly, judging whether a GPS station exists in a grid of 1 degree multiplied by 1 degree; secondly, inverting a crust deformation sequence in an unobserved grid by using a Random Forest (RF) on the grid which does not contain the GPS station; thirdly, obtaining TWSA in southwest area of China by utilizing all crust sequence inversions; fourthly, comparing the MLLIM with the GRACE and hydrological model results to verify the accuracy of the MLLIM method; and fifthly, the application MLLIMMLLIM is used for performing reverse performance on TWSA of each province in the southwest region of China, and analyzing the TWSA change relationship of each province by combining rainfall data, so that the method has important significance on effective management of water resources in the southwest region and life of people.

In 2014, Argus and the like propose a concept of inverting water reserve change by using a GPS, and are different from water reserve change monitoring means such as well position depth measurement, a hydrological model, satellite height measurement and a GRACE gravity satellite. The GPS is used for inverting the land water reserve change, the theory related to the deformation of the earth crust load is applied, the physical significance of GPS data can be greatly increased, and the space-time resolution can be effectively improved compared with the traditional GRACE gravity satellite method. However, the accuracy of TWSA inversion is low in a region where GPS sites are sparse due to uneven distribution of GPS sites, and the density of GPS sites greatly limits the application of TWSA inversion. Therefore, how to accurately simulate the vertical displacement deformation sequence of the sparse area of the GPS station is the key for improving the inversion accuracy.

In the embodiment, as shown in fig. 1, the method for improving the accuracy of land-water reserve inversion based on the machine learning load model includes:

step 101, dividing the research area according to a grid of 1 degree x 1 degree.

And 102, judging whether the grid contains the GPS station.

In this embodiment, when it is determined that the grid includes a GPS site, step 103 is skipped to; when it is determined that no GPS site is included in the grid, the jump is performed to step 104.

And 103, preprocessing the vertical time sequence of the GPS station to obtain a real GPS vertical deformation sequence.

In this embodiment, the vertical time series corresponding to the GPS site can be extracted from the chinese regional GPS time series file provided by the chinese continental structure network CMONOC. However, due to the effects of earthquakes and receiver anomalies around the GPS site, there can be step terms and outliers in the time series. Therefore, the extracted vertical time sequence corresponding to the GPS site can be preprocessed to remove step terms and abnormal values caused by earthquakes around the GPSGPS site and receiver abnormity, and further obtain a real GPS vertical deformation sequence.

Preferably, the preprocessing process of the vertical time series corresponding to the GPS site can be as follows:

ΔU_gps＝U_gps(t)-U_gps(t-1)···(1)

Wherein t represents time, U_gps(t) vertical time series, U, corresponding to GPS station at time t_gpsAnd (t-1) represents a vertical time sequence corresponding to the GPS station at the time of t-1.

U_deleap＝U_gps(t)+ΔU_gps···(2)

And step 104, simulating the vertical deformation sequence of the crust in all grids of the research area based on a random forest method (RF) to obtain a simulated vertical deformation sequence of the crust.

Random forest method (RF) is a learning method proposed by Breiman in 2001 for tasks such as regression, classification, etc. Firstly, constructing a large number of random and decorrelated decision trees; second, the decision tree is averaged. For regression tasks, the main advantages include: (1) the predictor is selectable based on data availability and user demand; (2) possibly including successively classified predictors; (3) relatively few model parameters must be specified by the user; (4) minimum risk of overfitting; (5) a variable importance score is automatically calculated to evaluate the contribution of a single predictor to the final model.

In this embodiment, the temperature and the air pressure corresponding to the GPS site may be used as input data, the vertical time sequence of the GPS site may be used as output data, and the random forest method is used to sequentially regress for M times to obtain output results of M grids; and then, carrying out average processing on output results of the M grids to obtain simulated crustal vertical deformation sequences in all grids of the research area.

Preferably, the sequence of simulated crustal vertical deformations is represented as follows:

And 105, correcting the atmospheric and non-tidal ocean loads of the real GPS vertical deformation sequence obtained in the step 103 or the simulated crustal vertical deformation sequence obtained in the step 104 to obtain a corrected crustal vertical deformation sequence.

In this embodiment, considering that some components in the crustal deformation sequence (the real GPS vertical deformation sequence, the simulated crustal vertical deformation sequence) are affected by the atmospheric and non-tidal marine loads, in order to accurately extract the influence of the hydrological load on the deformation of the crustal, atmospheric and non-tidal marine loads can be corrected on the crustal deformation sequence in all grids by using MERRA and OMCT models.

Preferably, for the real GPS vertical deformation sequence, there are:

U_hy＝U_deleap-U_MERRA-U_OMCT···(4)

Preferably, the simulated crustacean homeotropic order sequence is as follows:

U_hy＝U_grid-U_MERRA-U_OMCT···(5)

And step 106, taking the corrected vertical deformation sequence of the crust as input data of a crust load model, and performing inversion on the land water reserves TWSA.

The earth is a flexible sphere, and when the load (such as surface water, snow, ice, etc.) on the surface of the earth changes, the crust deforms to a corresponding degree, and the deformation is load deformation. The green's function can be used to establish a relationship between load mass and deformation. The earth crust load deformation is mainly shown in the horizontal direction and the vertical direction, however, the earth crust load deformation is more sensitive in the vertical direction, and the load deformation amplitude is about 2-3 times of that in the horizontal direction. Green function U_greenThe vertical load deformation of the crust is described 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.

Γ_nThe derivation of the function is as follows:

when n is 0, Γ_nThe expression of the function is as follows:

fig. 2 shows the relationship between load deformation and distance when discs of different mass, radius and thickness are placed on the ground. As can be seen from fig. 2, the load response in the near field is very significant. When the distance from the center of the disc is equal to the radius of the disc, the load response is only 1/2 of the disc center response; this response is negligible when the distance from the center of the disc is 5 times the radius of the disc. It follows that GPS sites can only reverse land-water reserves within a limited range. Therefore, it is important to simulate the vertical load deformation sequence of the unknown mesh.

The MLLIM method estimates the water reserve change for a 0.25 x 0.25 grid based on all corrected sequences of the crustal vertical deformation as data. The obtained solution is then regularized using a curvature smoothing algorithm and added to the solution matrix as a set of constraints. In particular, during each time of this study, the sequence U of vertical deformations of the crust after correction is determined_hyMinimizing the least square inhibition problem, and estimating daily land water reserve change Load_TWSA：

Load_TWSA＝((Gx-U_hy)/σ)²+β²(L(x))²→min···(6)

On the basis of the above embodiment, the verification of the method for improving the inversion accuracy of land water reserves based on the machine learning load model (MLLIM method) is described below.

Data and models

The invention utilizes the coordinate time sequence of the GPS station provided by the Chinese continental environment monitoring network to carry out analysis. The southwest region of China has 57 sites, and the spatial distribution of the sites is shown in FIG. 2. The coordinate time sequence is obtained by taking an original observation file, a navigation file and a precise ephemeris file as data bases, utilizing GAMIT software to carry out ionosphere correction, absolute antenna phase center correction, ocean tide correction and the like on the data bases, solving a base line between sites, and solving site coordinates under an ITRF2008 framework by using GLOBK adjustment software, so as to obtain a GPS site coordinate time sequence. And then, carrying out step item correction on the obtained vertical coordinate time series and removing the sequence abnormal value which is more than three times of standard deviation. In order to avoid the influence of data missing on TWSA inversion, the invention utilizes the longest time span shared among GPS sequences as 2011-2019. As can be seen from fig. 3, GPS sites are the most dense in the southwest region of china, and are relatively sparse in the eastern part of the world, Chongqing, Guangxi, etc.

Inversion model for gravity satellite data and land water reserve change

The change of land water reserves in a global area can be effectively inverted by using the time-varying gravity field of GRACE and GRACE-FO, and the calculation formula is as follows:

wherein, deltah represents the equivalent water height of land water reserve change; a represents the earth radius 6371.39 km; rho_eRepresents the earth average density of 5.51 × 10³kg/m³；ρ_wDenotes the density of water 10³kg/m³；h_lAnd k_lRepresenting the load lux number of order l; w_lRepresenting a gaussian smoothing kernel function;

representing a fully normalized associated legendre function; delta C_lmAnd Δ S_lmRepresenting the variation of the spherical harmonic coefficients of the earth's gravitational field.

The invention utilizes the GRACE Mascon (GRACE-M) data products provided by the space center of Texas university (CSR) in the United states and the Jet Propulsion Laboratory (JPL) of the aerospace administration in the United states between 2011 and 2019 years, and extracts the TWSA in the southwest region of China according to the boundary file. In order to weaken the influence of the resolving strategy on the uncertainty of the data product, the invention takes the average value of two data as the final TWSA result:

GLDAS data and barometric pressure data

In order to effectively return a GPS vertical time sequence of an unknown region, V2.2 global land assimilation data provided by a global land assimilation system (GLADAS) and daily scale air pressure data provided by a European middle weather forecast center (ECMWF) are used as input data, and then the daily scale earth surface temperature and the daily scale air pressure of the southwest region of China are respectively extracted. Wherein, the time resolution of the GLDAS V2.2 land assimilation data set is day scale, the space resolution is 0.25 degrees, the time span is 2003 to the present, and the data set comprises 28 variables including temperature; ECMWF provides barometric pressure data with a time resolution on a day scale, a spatial resolution of 0.1 °, and a time span of 2000 years to the present. FIG. 3 depicts a GPS site vertical time series, temperature time series, and barometric time series, exemplified by SCPZ, SCXJ, YNMJ, and YNML.

As can be seen from FIG. 4, the cycle length of the temperature, pressure and vertical deformation sequences is close to one year. The temperature sequence and the GPS vertical deformation sequence form a negative correlation relationship, and a certain phase difference between the air pressure sequence and the crustal vertical sequence can be obviously seen. Overall, the periodicity of the temperature and barometric time series is more pronounced than the crustal deformation series, due to the more complex factors that influence the crustal deformation.

Evaluation index

The present invention utilizes Root Mean Square Error (RMSE), Pearson Correlation Coefficient (PCC), and R-squared (R)²) And (3) evaluating an inversion result:

wherein Y and S represent real data and simulated data, respectively,

and

mean values of the data are indicated. RMSE describes the dispersion between the actual and simulated sequences, and if the RMSE value is small, it indicates that the simulation result is stable in amplitude. PCC is a linear correlation coefficient that reflects the degree of linear correlation between two quantities. The PCC value is between-1 and 1, the closer the absolute value is to 1, the stronger the correlation and the higher the accuracy of the reaction. It mainly extracts the seasonal relationship between sequences. At the same time, R²Also referred to as the decision coefficient, which represents the ratio of the sum of squared deviations of the total to the sum of squared deviations that can be explained by the sum of squared deviations of the regression. R²Values between 0 and 1, with larger values giving better fit, indicating higher accuracy.

TWSA calculation based on MLLIM

The deformation of the crust can be divided into a horizontal (N, E) direction and a vertical (U) direction, wherein the sequence in the horizontal direction is mainly characterized by linear variation, and the characteristic of the crust vertical direction is mainly characterized by non-linear seasonal variation. The GPS receiver is fixed on bedrock, and the bedrock near the earth surface generates periodic up-and-down vibration due to the temperature change of the earth surface due to the phenomenon of expansion with heat and contraction with cold. And due to the difference of the altitude of the station, the air pressure of the area is also different, and the difference of the air pressure can influence the amplitude of the GPS station. Therefore, the invention uses the surface temperature data provided by GLDAS V2.2 and the air pressure data provided by CFWCM as input sequences, and uses the GPS station vertical displacement sequence as output data to carry out regression. In order to weaken the influence of regression errors of a single site, the method utilizes an averaging strategy after multiple regressions of all sites to carry out regression on the vertical deformation displacement of the unknown grid.

To verify the usability of the regression method, the temperature series, the pressure series and the time of the grid on which 57 GPS stations are located are used as input data, the GPS vertical series is used as input data, and regression is performed 56 times for each GPS vertical time series (excluding the station itself). As shown in fig. 5(a) and 5(b), the regression results were evaluated using both PCC and RMSE indicators; to demonstrate the effectiveness of the regression method as a whole, the mean of the vertical simulated sequence and the true sequence of 57 sites was found, as shown in fig. 5 (c).

As can be seen from FIG. 5, the method can effectively simulate the sequence of deformation of the earth crust load. Further statistics are given for the data in fig. 5(a), where sequences with correlation coefficients above 0.5 account for 84.21% of the population, with a maximum of up to 0.79. As can be seen in fig. 5(b), the RMSE between the simulated and the real sequences performed well, with 68.42% of the sites having an RMSE below 5mm, the value of which correlates with the quality of the GPS sequence. As shown in FIG. 5(c), in order to comprehensively express the applicability of the strategy to simulate the deformation of earth crust load, the mean values of 57 simulated and real value sequences in the research area are respectively calculated, and it can be seen that the method can effectively simulate the period term and amplitude of the sequence, the correlation reaches 0.75, and the RMSE is 3.45 mm. The accuracy can effectively simulate the GPS vertical deformation sequence, so the crustal vertical deformation sequence in an unknown area is regressed by using the method. And each grid is regressed 57 times, and the mean value of each grid is calculated to be used as the final vertical deformation sequence of the grid, so that a solid data base is provided for the next land water reserve inversion.

In order to invert TWSA in the southwest area of China by using the GPS vertical deformation sequence, the invention inverts TWSA in the area by using the corrected 57 GPS station vertical sequences and 87 unknown area simulated crustal vertical deformation sequences, wherein the simulated station points are limited according to the central point of the unknown area grid, as shown in FIG. 6.

As can be seen from FIG. 6, the anniversary term and the semianniversary term of the vertical deformation sequence of the unknown grid can be effectively simulated based on the strategy. The crustal vertical deformation sequences of all grids in the research area are calculated through the regression method, and preparation is made for the TWSA inversion. Firstly, a green function of point load in a research area is calculated, 0.25 degrees is selected as grid resolution, the radius of a corresponding disc is 13.90km, and preset parameters of expanding a boundary range of 2 degrees and beta being 0.01 are selected. And (3) inverting the terrestrial water reserves of 2011-2019 in the southwest region of China by using a formula (5) to obtain a terrestrial water reserve change value with a spatial resolution of 0.25 degrees, and calculating to obtain the TWSA and the annual amplitude of change of the TWSA in the southwest region of China, as shown in FIG. 7.

As can be seen from FIG. 7, the change of the land water reserves in and around the research area can be inverted by using the novel machine learning load inversion method, FIG. 7(b) shows that the annual amplitude of the land water reserves in the range of Yunnan province is obviously higher than that of the surrounding provinces, and the amplitude reaches 120mm, and FIG. 7(c) shows that the annual item and the amplitude of the land water reserve abnormal sequence can be inverted based on the novel machine learning load inversion method. And the land water storage amount in the Yunnan province range is obviously larger than the surrounding provinces, and the land water storage amount is equivalent to 120 mm. In order to prove the accuracy of the inversion result, the land water reserve change obtained by inversion by the method is respectively compared with the land water reserve change result and the hydrological model obtained by the traditional GPS inversion TWSA method and the GRACE inversion.

Comparing MLLIM inversion result with traditional GPS inversion method result

In order to discuss the difference between a novel machine learning load inversion method and a traditional GPS inversion method, the method takes a crustal vertical deformation sequence with atmospheric load deformation deducted as input data, and combines with a Green load function to solve and obtain TWSA (0.25 degrees multiplied by 0.25 degrees) in the southwest region of China between 2011 and 2019. Next, the yearly amplitude, the semiyearly amplitude and the sequence period of each lattice point sequence are solved and compared with the experimental results of the present invention, as shown in fig. 8.

As can be seen from fig. 8, the results based on the new machine learning load inversion method are significantly different from the TWSA results of the conventional GPS inversion. The subgraph in fig. 8 can be divided into three contrast groups: (a) and (d), (b)/(e) and (c)/(f) represent three sequence characteristic parameters, i.e., the annual amplitude (I), the semiannual amplitude (II) and the sequence period (III), respectively. As shown in fig. 8 (I), since the GPS sites in the Guangxi province of China are fewer than other provinces, when the TWSA is inverted by the conventional GPS method, some feature points are ignored due to the smoothing process, but the annual amplitude feature is clearly shown in the graph (a). And this phenomenon is more apparent in fig. 8 group (II), and signal feature loss due to interpolation can be effectively suppressed. When the period of the sequence is analyzed (fig. 8 group III), the funnel condition of the anniversary item of Guangxi province can be effectively obtained by inversion by using a novel machine learning load inversion method for the obvious difference of the parts of Guangxi province and Guizhou province. In order to further judge the accuracy of the two inversion methods, the TWSA results obtained by the two methods are compared with the GRACE gravity satellite inversion result and the hydrological model.

Comparison of MLLIM inversion results with GRACE results

The GRACE and GRACE-FO satellites can invert the land water reserve change with higher precision by monitoring the change of the earth time-varying gravity field. In order to verify the accuracy of the inversion result, Mascon data results issued by JPL and CSR mechanisms are utilized, and the mean value of the Mascon data results is obtained and used as the TWSA result of the final southwest area. Since there is a period of vacancy between the GRACE satellite and the GRACE-FO satellite, the GRACE-M product is divided into two periods, GRACE and GRACE-FO. And the time resolution of the GRACE inversion result is month scale, so that the MLLIM inversion result is subjected to month average processing for comparison with the GRACE inversion result.

As can be seen in FIG. 9, the computed results of the GRACE-M data are consistent with the MLLIM inversion results in the anniversary amplitude and the semianniversary amplitude. The overall annual amplitude of Yunnan province is the largest and can reach about 120mm at most, and the annual amplitude of Guangxi province is also more prominent, so that the novel machine learning load inversion method is accurate for simulating the annual amplitude signal (figure 8(a)) of Guangxi province, and the signal is ignored in the traditional GPS inversion method. Fig. 9(b) shows that the detection effect of the sequence semiannuity signal is more prominent, and the bumps with semiannuity amplitude signals in the middle of Yunnan province, Guizhou province and Sichuan province can be obviously seen. However, there is some variation in the annual phase of the calculation sequence due to the different resolution of the GPS inversion results and the GRACE-M data, but the overall magnitude remains the same. Respectively calculating TWSA of two methods in southwest area of China, fitting the sequence, calculating correlation coefficients of GPS inversion result and GRACE-FO data to be 0.91 and 0.88 respectively, and calculating R of the correlation coefficients²The values are 0.76 and 0.65, respectively. The TWSA results obtained by the traditional GPS inversion method and the PCC of GRACE and GRACE-FO are respectively 0.87 and 0.79; r²The values are 0.71 and 0.58, respectively. Therefore, compared with the traditional GPS (global positioning system) inversion TWSA (time-of-flight ranging) method, the method is based on MLLIM (maximum likelihood of being limited) in PCC and R²The indexes of (a) are improved by 7.98% and 9.30% on average.

Comparing MLLIM inversion result with GLDAS result

In order to better verify the accuracy of the novel machine learning load inversion method, the TWSA variables in the GLDAS V2.2 dataset are used for comparison. The spatial resolution of the GLDAS V2.2 data set is 0.5 degrees multiplied by 0.5 degrees, the time resolution is a day scale, and the GLDAS V2.2 data set can be well compared with an MLLIM inversion result.

As can be seen from fig. 10, the GLDAS data is consistent with the MLLIM inversion results as a whole in terms of annual amplitude, half-cycle amplitude and annual phase. As can be seen from fig. 10(a), the annual amplitudes in Yunnan province, western Sichuan province and eastern Guangxi province all have significant projections, which are consistent with the MLLIM inversion results (fig. 8 (a)). It can be seen from fig. 10(b) that, in accordance with the present result (fig. 8(b)), the half-cycle amplitude salient points of Yunnan province, Chongqing city, and Guizhou province were detected. From fig. 10(c), it can be seen that the annual amplitude is generally higher in the areas of Guangxi province and eastern Sichuan, which is consistent with the experimental result, and the accuracy of the annual phase of the inversion result is further proved. And calculating PCC and R between MLLIM inversion result and GLDAS data²0.79 and 0.64, respectively. PCC and R between traditional GPS inversion method and GLDAS data²0.72 and 0.60, respectively, so that the MLLIM-based method is more traditional GPS inversion method for PCC and R²The improvement is respectively 9.72% and 6.67%, which shows that the MLLIM can improve the accuracy of TWSA inversion compared with the traditional GPS inversion method.

Applications of

The inversion result of the invention and PCC of GRACE and GLDAS are respectively 0.88 and 0.79, and the annual amplitude, semi-annual amplitude and annual phase signal of TWSA time sequence can be effectively detected. Because the vertical load deformation of the crust is mainly related to the water storage of the land, when the water amount is increased, the crust generates downward vertical deformation; on the contrary, the earth crust generates upward rebound deformation. In order to analyze the reason of the change of the water reserves in the southwest area of China, the invention utilizes a 0.5-degree grid product provided by a China meteorological data network (CMA) to respectively extract the average rainfall of each province, and carries out comparative analysis with the inversion result of the invention, GRACE Mascon and GLDAS data, as shown in FIG. 11.

As can be seen from FIG. 11, the results of the MLLIM inversion have better spatio-temporal consistency with GRACE, GLDAS and rainfall data. Here, since there is an abnormal situation in the GLDAS data sets in Chongqing city and Guizhou province (as can be seen in fig. 10 (c)), the GLDAS data in these two regions are not compared. From fig. 11(a) - (e) it can be seen that: there are significant amplitude humps in rainfall in the southwest part of Yunnan province, the middle part of Sichuan, and the southeast part of Guangxi, consistent with the TWSA results of the inversion of the present invention (FIG. 8 (a)). As can be seen in FIGS. 11(f) - (j): (1) the result of the invention is almost consistent with the periodicity of GRACE and GLDAS, and the experimental result of the invention is better kept with the phase of GLDAS data, however, the phase difference of about 2 months exists between the experimental result and the GRACE data due to the difference of observation means; (2) the amplitude of the GPS inversion result is slightly larger than GRACE and GLDAS, mainly because the signals with the atmospheric load and the non-tidal ocean load subtracted still have the influence of common mode error caused by observation means and the like; (3) through further statistics, the maximum value and the minimum value of the correlation between the inversion result and the GRACE are 0.86 (Yunnan province), the minimum value is 0.63 (Sichuan province), the average value is 0.77, the maximum value and the minimum value of the correlation between the inversion result and the GLDAS are 0.94 (Yunnan province), 0.78 (Guangxi province) and 0.86 respectively. Therefore, the change situation of the TWSA in the southwest area of China can be effectively inverted based on the MLLIM.

However, due to the particularity of the GPS observation means, various noise components exist in the crustal vertical deformation sequence observed by the GPS, so that the TWSA obtained by inversion has a noise part, and therefore, the noise part of the GPS sequence can be better removed in subsequent research; and because the characteristic of the TWSA mode of GPS inversion has certain phase difference with GRACE, the phase difference part is processed through signals in subsequent research so as to achieve better fitting with GRACE.

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

1. A method for improving land water reserve inversion accuracy based on a machine learning load model is characterized by comprising the following steps:

2. The method for improving inversion accuracy of land water reserves based on the machine learning load model of claim 1, wherein the vertical time sequence of the GPS station is preprocessed to obtain a real GPS vertical deformation sequence, comprising:

3. The method for improving inversion accuracy of land water reserves based on the machine learning load model of claim 2, wherein the preprocessing process of the vertical time series corresponding to the GPS station is as follows:

ΔU_gps＝U_gps(t)-U_gps(t-1)···(1)

U_deleap＝U_gps(t)+ΔU_gps···(2)

4. The method for improving inversion accuracy of land water reserves based on the machine learning load model of claim 3, wherein the method is characterized in that a simulated crustal vertical deformation sequence in all grids of a research area is simulated based on a random forest method to obtain a simulated crustal vertical deformation sequence, and comprises the following steps:

5. The method for improving inversion accuracy of land water reserves based on the machine learning load model of claim 4, wherein the simulated crustal vertical deformation sequence is represented as follows:

6. The method for improving inversion accuracy of land water reserves based on the machine learning load model of claim 5, wherein the procedure for atmospheric and non-tidal marine load correction of true GPS vertical deformation sequences is as follows:

U_hy＝U_deleap-U_MERRA-U_OMCT···(4)

7. The method for improving inversion accuracy of land water reserves based on the machine learning load model of claim 5, wherein the atmospheric and non-tidal marine load correction of the simulated crustal vertical deformation sequence is performed as follows:

U_hy＝U_grid-U_MERRA-U_OMCT···(5)

8. The method for improving the accuracy of land-water reservoir inversion based on the machine learning load model according to claim 6 or 7, wherein the land-water reservoir TWSA is inverted by using the obtained corrected crustal vertical deformation sequence as input data of the crustal load model, and the method comprises the following steps:

Load_TWSA＝((Gx-U_hy)/σ)²+β²(L(x))²→min···(6)