A method to improve the accuracy of GRACE land water storage changes

Abstract

The present invention discloses a method for improving the accuracy of GRACE terrestrial water storage change, comprising: obtaining terrestrial water storage Δh(θ,λ) from a hydrological model, and determining the spherical harmonic coefficient cosine term C corresponding to Δh(θ,λ) lm and the sine term S lm ; for C lm and S lm Perform de-striping filtering and Gaussian smoothing filtering to obtain the filtered spherical harmonic coefficients cosine term #imgabs0# and sine term #imgabs1#. Based on #imgabs2# and #imgabs3#, obtain the filtered terrestrial water storage Δη(θ,λ); determine the time series data Δh corresponding to Δh(θ,λ) in the time series t i (t), Δη(θ,λ) corresponds to the time series data Δη at time t i (t); according to Δh i (t) and Δη i (t), and the determination coefficient weighted scaling factor is obtained by solving the problem. The leakage error in the process of GRACE data inversion of terrestrial water storage changes is corrected according to the determination coefficient weighted scaling factor, and the real hydrological signal is restored to improve the accuracy of GRACE terrestrial water storage changes.

Description

Technical Field

[0001] The present invention belongs to the cross-technical fields of satellite gravity, hydrology, etc., and in particular relates to a method for improving the accuracy of GRACE terrestrial water storage changes. Background Art

[0002] Since the launch of the Gravity Recovery and Climate Experiment (GRACE) satellite in 2002, the GRACE satellite has provided a new means of observing the Earth's time-varying gravity field signals with high precision. So far, the GRACE satellite has been widely used in the study of terrestrial water storage (TWS) changes, glacier changes and their mass balance, global average ocean mass changes, earthquakes and other fields.

[0003] Based on the theory of Wahr et al., the mass change of the Earth's surface can be inverted using the GRACE spherical harmonic coefficient solution. However, the spherical harmonic coefficients contain both random errors that increase with the order of the coefficients and systematic errors related to the system. Therefore, many researchers have proposed spatial smoothing filters (such as Gaussian smoothing filters, Wiener filters, etc.) and de-striping filters (Swenson and Wahr (2006) de-striping method, Chen et al. (2009) P4M6 method, etc.). These filtering methods will change the real physical signal while eliminating the influence of errors. Chen et al. (2015) conducted experimental simulations and found that the leakage error caused by truncating and smoothing 300 km can reach 80.95Gt / yr in Antarctica (measured value -150.00Gt / yr, simulated value -69.05Gt / yr); Jin and Zou (2015) showed that 42.4% of the signal of mass loss of the Greenland ice sheet was caused by leakage error. After correcting the leakage error, the melting rate of Greenland's glaciers increased from -15.58±0.93Gt / yr2 to -26.19±1.67Gt / yr2. Therefore, leakage error is the main limiting factor in the quality change of GRACE inversion, and correcting leakage error is crucial to improving the accuracy of inversion results.

[0004] At present, the commonly used recovery methods for GRACE data leakage errors include: (1) Scaling factor method SF (Scaling Factor). The calculation process is simple. Long et al. (2015) used the scaling factor obtained by the global hydrological model PCR-GLOBWB to restore the signal loss of water storage changes in the Yangtze River Basin. This method is of great value for improving the accuracy of water storage changes in the study area; Huang et al. (2019) conducted a sensitivity analysis of the scale factor leakage error correction method based on different global hydrological models and evaluated the water storage changes in southwestern China. The scaling factor method relies on a priori hydrological models. Although hydrological models have been widely used in a priori information for GRACE data leakage error correction, hydrological models have some shortcomings, such as incomplete representation of surface water and groundwater components and failure to consider the impact of human activities. (2) Iterative recovery method FM (Forward Modeling). Mainly constrained by the original GRACE observations, no prior model is required. Chen et al. (2014) used the iterative recovery method to restore the long-term groundwater changes in northwest India. The newly estimated water storage change trend reflects that the iterative recovery method can effectively reduce the leakage error in GRACE estimation; Wu Yunlong et al. (2015) can effectively restore the initial signal of water storage change in the Heihe River Basin based on the iterative recovery method, reduce the leakage error, and improve the spatial resolution of water storage changes. The calculation process of the iterative recovery method is relatively complex, and in some areas, there may be incomplete signal recovery or excessive recovery. Summary of the invention

[0005] The technology of the present invention solves the problem: overcomes the shortcomings of the existing technology, provides a method for improving the accuracy of GRACE terrestrial water storage changes, constructs a new scale factor determination coefficient weighted average correction method NSFCM (New Scale Factor Determination Coefficient Weighted Average Correction Method), fully considers the information of scale factors of various models, and can effectively improve the accuracy of GRACE terrestrial water storage changes compared with traditional scale factor methods and iterative recovery methods.

[0006] In order to solve the above technical problems, the present invention discloses a method for improving the accuracy of GRACE terrestrial water storage changes, comprising:

[0007] Obtain the terrestrial water storage Δh(θ,λ) from the hydrological model and determine the spherical harmonic coefficient cosine term C corresponding to Δh(θ,λ) lm and the sine term S lm ; Among them, θ represents co-latitude and λ represents longitude;

[0008] C lm and S lmPerform de-striping filtering and Gaussian smoothing filtering to obtain the filtered spherical harmonic coefficient cosine term and the sine term

[0009] based on and Get the filtered terrestrial water storage Δη(θ,λ);

[0010] Determine the time series data Δh corresponding to Δh(θ,λ) in the time series t i (t), Δη(θ,λ) corresponds to the time series data Δη at time t i (t);

[0011] According to Δh i (t) and Δη i (t), solve to get the weighted scaling factor of the coefficient of determination;

[0012] The leakage error in the process of inverting terrestrial water storage changes from GRACE data is corrected according to the weighted scaling factor of the determination coefficient, and the real hydrological signal is restored to improve the accuracy of GRACE terrestrial water storage changes.

[0013] In the above method to improve the accuracy of GRACE terrestrial water storage changes, the spherical harmonic coefficient cosine term C corresponding to Δh(θ,λ) is determined by the following formula (1): lm and the sine term S lm :

[0014]

[0015] Among them, ρ w represents the density of water, a represents the average radius of the earth, ρ ave represents the average density of the Earth, represents the normalized associated Legendre function, l represents the order of the spherical harmonic coefficient solution inversion, m represents the number of times the spherical harmonic coefficient solution inversion is performed, and k l It represents the load Love number of order l.

[0016] Among the above methods to improve the accuracy of GRACE terrestrial water storage changes, the C lm and S lm Perform de-striping filtering and Gaussian smoothing filtering to obtain the filtered spherical harmonic coefficient cosine term and the sine term include:

[0017] According to the following equations (2) and (3), C lm and S lm Perform de-striping filtering to obtain the cosine term of the spherical harmonic coefficients after de-striping filtering and the sine term

[0018]

[0019]

[0020] Where w is the sliding window width, e is the natural logarithm, max() means taking the maximum value, Q is the destriping filter coefficient; for destriping of odd items, n is the odd item in the lth order, and the same is true for even items;

[0021] According to the following equations (4) and (5), and Perform Gaussian smoothing filtering to obtain the cosine term of the spherical harmonic coefficients after Gaussian filtering and the sine term

[0022]

[0023]

[0024] Among them, r represents the Gaussian smoothing radius, W0 represents the initial order Gaussian smoothing weight, and W l-1 , W l and W l+1 They represent the Gaussian smoothing weights of order l-1, l and l+1 respectively, and l is 0, 1, 2, 3…60.

[0025] In the above method for improving the accuracy of GRACE terrestrial water storage changes, the filtered terrestrial water storage Δη(θ,λ) is obtained by the following formula (6):

[0026]

[0027] In the above method to improve the accuracy of GRACE terrestrial water storage changes, according to Δh i (t) and Δη i (t), the coefficient of determination weighted scaling factor is obtained by solving, including:

[0028] According to Δh i (t) and Δη i (t), based on the least squares principle, the scale factor K is calculated;

[0029] According to Δh i (t), and the coefficient of determination R is obtained by solving 2 ;

[0030] According to the scale factor K and the determination coefficient R 2 , the determination coefficient weighted scaling factor NSFCM-SF is obtained by solving the equation.

[0031] In the above method for improving the accuracy of GRACE terrestrial water storage changes, the scale factor K is calculated by the following equation (7):

[0032]

[0033] Among them, N represents the number of time series, and i represents the i-th time series.

[0034] In the above method to improve the accuracy of GRACE terrestrial water storage changes, according to Δh i (t), and the coefficient of determination R is obtained by solving 2 ,include:

[0035] According to Δh i (t), determine the fitted terrestrial water storage and the average terrestrial water storage

[0036] The determination coefficient R is obtained by solving the following formula (8): 2 :

[0037]

[0038] In the above method for improving the accuracy of GRACE terrestrial water storage changes, the coefficient of determination weighted scaling factor NSFCM-SF is obtained by solving the following equation (9):

[0039]

[0040] Where n represents the number of scale factor models, j represents the jth scale factor model, and K j represents the scale factor corresponding to the j-th scale factor model, Represents the determination coefficient corresponding to the j-th scale factor model.

[0041] The present invention has the following advantages:

[0042] The present invention discloses a method for improving the accuracy of GRACE terrestrial water storage changes, which adopts a weighting technology based on a coefficient of determination, uses terrestrial water storage data of multiple hydrological models to simulate the signal loss caused by the filtering process during GRACE data inversion, calculates the scale factors of the terrestrial water storage data of multiple hydrological models before and after filtering by the least square method, and weights the scale factors of each hydrological model based on the simulated coefficient of determination to obtain a new scale factor. The new scale factor integrates the terrestrial water storage information of multiple hydrological models, more accurately reflects the signal loss caused by the filtering process during GRACE data inversion, and uses the new scale factor to correct the leakage error in the process of GRACE data inversion of terrestrial water storage changes, which can effectively restore the real hydrological signal to improve the accuracy of GRACE terrestrial water storage changes. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 is a flowchart of steps of a method for improving the accuracy of GRACE terrestrial water storage changes in an embodiment of the present invention;

[0044] Figure 2 is a schematic diagram of the location of the Yangtze River Basin and its altitude in an embodiment of the present invention;

[0045] Figure 3 It is a schematic diagram of the change in the spatial distribution trend of TWS in the Yangtze River Basin calculated by a GLDAS-NOAH model in an embodiment of the present invention;

[0046] Figure 4 is a schematic diagram of the Yangtze River Basin scale factor calculated by a different model in an embodiment of the present invention;

[0047] Figure 5 is a schematic diagram of restoring the GRACE mass change signal using an iterative recovery method in an embodiment of the present invention;

[0048] Figure 6 is a schematic diagram of the spatial distribution trend of different restoration methods in the Yangtze River Basin according to an embodiment of the present invention;

[0049] Figure 7 is a schematic diagram of a TWS time series restored by different leakage error correction methods in the Yangtze River Basin in an embodiment of the present invention;

[0050] Figure 8 is a schematic diagram of the annual amplitude and long-term trend of TWS restored by different leakage error recovery methods in the Yangtze River Basin in an embodiment of the present invention; wherein, Figure 8 (a) is the annual amplitude, Figure 8 (b) long-term trend;

[0051] Fig. 9 is a schematic diagram of a TWS time series of the Yangtze River sub-basin in an embodiment of the present invention; wherein, Fig. 9 (a) The upper reaches of the Yangtze River. Fig. 9 (b) The middle and lower reaches of the Yangtze River;

[0052] Fig.10 is a schematic diagram of a time series of precipitation in the Yangtze River Basin and TWS in an embodiment of the present invention; wherein, Fig.10 (a) The Yangtze River Basin, Fig.10 (b) The upper reaches of the Yangtze River, Fig.10 (c) The middle and lower reaches of the Yangtze River;

[0053] Fig.11 It is a schematic diagram of annual precipitation anomaly in the Yangtze River Basin in an embodiment of the present invention. DETAILED DESCRIPTION

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

[0055] The present invention discloses a method for improving the accuracy of GRACE terrestrial water storage changes, which can be mainly divided into the following steps: (1) converting hydrological grid data into spherical harmonic data; (2) performing de-striping filtering and Gaussian smoothing filtering; (3) converting spherical harmonic coefficients into grid data; (4) calculating a scale factor based on the least squares principle; (5) calculating a model determination coefficient; and (6) weighting the scale factor based on the determination coefficient.

[0056] like Figure 1 In this embodiment, the method for improving the accuracy of GRACE terrestrial water storage changes includes:

[0057] Step 101, obtain the terrestrial water storage Δh(θ,λ) from the hydrological model, and determine the spherical harmonic coefficient cosine term C corresponding to Δh(θ,λ) lm and the sine term S lm .

[0058] In this embodiment, the spherical harmonic coefficient cosine term C corresponding to Δh(θ,λ) can be determined by the following formula (1): lm and the sine term S lm :

[0059]

[0060] Among them, θ represents the co-latitude, λ represents the longitude, and ρ w represents the density of water, a represents the average radius of the earth, ρ ave represents the average density of the Earth, represents the normalized associated Legendre function, l represents the order of the spherical harmonic coefficient solution inversion, m represents the number of times the spherical harmonic coefficient solution inversion is performed, and k l It represents the load Love number of order l.

[0061] Step 102, for C lm and S lm Perform de-striping filtering and Gaussian smoothing filtering to obtain the filtered spherical harmonic coefficient cosine term and the sine term

[0062] In this embodiment, in order to maintain the same filtering as that of GRACE data processing, Swenson de-striping filtering and 300 km Gaussian smoothing filtering are used.

[0063] Preferably, the de-striping filter usually adopts a quadratic polynomial filter window centered at order l and with a width of w. The relationship between the sliding window width w and the spherical harmonic coefficient order m is:

[0064]

[0065] That is, according to the following equations (2) and (3), C lm and S lm Perform de-striping filtering to obtain the cosine term of the spherical harmonic coefficients after de-striping filtering and the sine term

[0066]

[0067]

[0068] Wherein, e is the natural logarithm, max() indicates taking the maximum value, and Q is the de-striping filter coefficient; for de-striping of odd items, n is the odd item in the lth order, and the same is true for even items.

[0069] Preferably, Gaussian smoothing filtering refers to a method of assigning different weights to spherical harmonic coefficients of different orders, thereby reducing the error of high-order terms of GRACE data. According to the following equations (4) and (5), and Perform Gaussian smoothing filtering to obtain the cosine term of the spherical harmonic coefficients after Gaussian filtering and the sine term

[0070]

[0071]

[0072] Among them, r represents the Gaussian smoothing radius, W0 represents the initial order Gaussian smoothing weight, and W l-1 , W l and W l+1 They represent the Gaussian smoothing weights of order l-1, l and l+1 respectively, and l is 0, 1, 2, 3…60.

[0073] Step 103, based on and The filtered terrestrial water storage Δη(θ,λ) is obtained.

[0074] In this embodiment, and Converted to grid data, that is, through the following formula (6), the filtered terrestrial water storage Δη(θ,λ) is obtained:

[0075]

[0076] Step 104, determine the time series data Δh corresponding to Δh(θ,λ) in the time series t i (t), Δη(θ,λ) corresponds to the time series data Δη at time t i (t).

[0077] Step 105, according to Δh i (t) and Δη i (t), and the determination coefficient weighted scaling factor is obtained by solving.

[0078] In this embodiment, step 105 is specifically implemented as follows:

[0079] a) According to Δh i (t) and Δη i (t), based on the least squares principle, the scale factor K is calculated.

[0080] For any grid point, the unfiltered terrestrial water storage corresponding to the time series at time t is Δh i (t), the filtered terrestrial water storage is Δη i (t), through the least squares principle, we have:

[0081]

[0082] Among them, N represents the number of time series, and i represents the i-th time series.

[0083] b) According to Δh i (t), and the coefficient of determination R is obtained by solving 2 .

[0084] The coefficient of determination is a statistical indicator used to reflect the degree of regression fit of the regression model. It can be defined as the ratio of the discrete degree of the fitted value of the independent variable in the sample regression model to its true value. For any grid point, the coefficient of determination R 2 The specific calculation is as follows:

[0085] According to Δh i (t), determine the fitted terrestrial water storage and the average terrestrial water storage The determination coefficient R is obtained by solving the following formula (8): 2 :

[0086]

[0087] c) According to the scale factor K and the coefficient of determination R 2 , the determination coefficient weighted scaling factor NSFCM-SF is obtained by solving the equation.

[0088] Based on the scale factor K and the coefficient of determination R 2, the formula of the new scale factor determination coefficient weighted average correction method is as follows:

[0089]

[0090] Where n represents the number of scale factor models, j represents the jth scale factor model, and K j represents the scale factor corresponding to the j-th scale factor model, Represents the determination coefficient corresponding to the j-th scale factor model.

[0091] Step 106, correcting the leakage error in the process of inverting the terrestrial water storage change using the GRACE data according to the weighted scaling factor of the coefficient of determination, restoring the real hydrological signal, and improving the accuracy of the GRACE terrestrial water storage change.

[0092] In summary, this paper uses the model determination coefficient weighted scale factor for the first time to construct a new scale factor determination coefficient weighted average correction method (NSFCM). The new scale factor determination coefficient weighted average correction method fully considers the information of multiple model scale factors and can improve the accuracy of GRACE terrestrial water storage changes compared to the traditional scale factor method and iterative recovery method.

[0093] Based on the above embodiment, a specific example is described below.

[0094] 1. Research area and data

[0095] 1.1 Study Area

[0096] The Yangtze River Basin is located in southern China, with a latitude and longitude range of 85°~122°E and 24°~36°N, spanning 19 provinces, municipalities and autonomous regions in the three major regions of western, central and eastern China. The location and altitude of the Yangtze River Basin are as follows: Figure 2 The average annual precipitation in the Yangtze River Basin is 1067 mm and its spatial and temporal distribution is uneven, with about 70% to 90% concentrated in May to October. The Yangtze River mainstream is located above Yichang, and the upper reaches of the Yangtze River have a drainage area of ​​1 million km 2 The middle reaches of the Yangtze River are from Yichang to Hukou, with a drainage area of ​​680,000 km 2 The area from Hukou to the Yangtze River estuary is the lower reaches of the Yangtze River, with a drainage area of ​​120,000 km 2 Given the small area of ​​the lower reaches of the Yangtze River, the inversion results may be highly uncertain if GRACE data are applied to areas smaller than the GRACE resolution. Therefore, the middle and lower reaches of the Yangtze River are collectively referred to as the middle and lower reaches of the Yangtze River. Figure 2The area east of the dotted line is shown, and the area west of the dotted line is the upper reaches of the Yangtze River.

[0097] 1.2 GRACE data

[0098] The data used in this invention include two grid data products: CSR-grid and CSR-M (CSR Mascon). CSR-grid data is based on the RL05 spherical harmonic coefficients of the Center for Space Research (CSR) at the University of Texas, the Jet Propulsion Laboratory (JPL) at the United States, and the GeoForschungsZentrum (GFZ) at Potsdam, Germany. The maximum order of the data is 60, and the time span is from April 2002 to January 2017. The monthly grid data represents the quality change relative to the average value from January 2004 to December 2009. GRACE's C 20 The value is better than the C of Satellite Laser Ranging (SLR) 20 The value has a greater uncertainty, so the SLR C 20 The coefficient solution replaces GRACE's C 20 Coefficients. The first-order terms were estimated using the method of Swenson et al. (2008) and applied to the geocentric correction. Glacial isostatic adjustment (GIA) correction was performed based on the model of Geruo et al. (2012). In order to reduce the correlation errors of the north-south stripes in the GRACE monthly map, the data were de-striped. In addition, a 300 km Gaussian smoothing filter was applied to the data to suppress the high-order random errors in the spherical harmonic coefficients.

[0099] The CSR-M data are calculated based on Level-1B data without any smoothing or de-striping filtering applied. Therefore, the CSR-M data can be used without further processing. The CSR-M dataset is not tailored for specific applications and is suitable for all scientific fields of interest, such as oceanography, land surface hydrology, cryosphere, etc. Studies have shown that the CSR-M data solution is suitable for recovering the mass redistribution of different mass sources at global and large basin scales (greater than about 3°×3°), while the spherical harmonic coefficient solution is superior to the CSR-M data solution at smaller or local scales when processed effectively.

[0100] 1.3 Hydrological data

[0101] The Global Land Data Assimilation System (GLDAS) hydrological model was jointly developed by NASA and the National Oceanic and Atmospheric Administration of the United States. The model outputs various land surface parameters (soil temperature and soil moisture, evaporation, surface runoff, etc.) and hydrological variables through land surface modeling and data assimilation technology. The GLDAS data used in the present invention includes four land surface models NOAH, CLM, MOS and VIC. The data spatial resolution is 1°×1°, the time resolution is 1 month, and the selected data time range is consistent with the GRACE data. The GLDAS model has soil moisture at different depths. VIC simulates 3 layers with a maximum depth of 1.90m; NOAH simulates 4 layers with a maximum depth of 2.00m; CLM simulates 10 layers with a maximum depth of 3.43m; MOS simulates 3 layers with a maximum depth of 3.50m. The model of the present invention takes soil moisture, snow water equivalent and total canopy water storage into account when calculating TWS, and the average value of the data from January 2004 to December 2009 is deducted from the monthly data to obtain the relative change of TWS.

[0102] 1.4 Precipitation data

[0103] The monthly precipitation data set for China is based on the latest compilation of precipitation data from China's ground stations by the National Meteorological Center's basic data project. The data is spatially interpolated using the thin disk spline method to generate monthly precipitation grid data for China with a resolution of 0.5°×0.5° from 1961 to the present, with a data range of 72°~136°E and 18°~54°N. The data time range selected by the present invention is from April 2002 to January 2017, which is consistent with the GRACE data time range.

[0104] 2. Verification of the new scale factor determination coefficient weighted average correction method

[0105] 2.1 Data processing

[0106] 2.1.1 Signal attenuation caused by truncation and filtering

[0107] The truncation, de-striping and Gaussian smoothing of the GRACE spherical harmonic coefficients will lead to leakage errors in the GRACE inversion results, which can be divided into two cases: (1) When the signal inside the study area is strong and the surrounding signal is weak, after the spherical harmonic coefficients are truncated and filtered, the signal inside the study area will leak out of the study area, indicating that the area is controlled by external leakage errors, which is called external leakage errors; (2) When the signal inside the study area is weak and the surrounding signal is strong, after the spherical harmonic coefficients are truncated and filtered, the surrounding signal will leak into the study area, indicating that the area is controlled by internal leakage errors, which is called internal leakage errors.

[0108] The present invention uses the GLDAS-NOAH hydrological model to calculate the spatial distribution trend of TWS, deducts the data outside the Yangtze River Basin in the study area, and simulates the external leakage error caused by truncation and filtering of spherical harmonic coefficients in the Yangtze River Basin. Figure 3 (a) is the change in the spatial distribution trend of TWS in the Yangtze River Basin calculated from the GLDAS-NOAH data from April 2002 to January 2017. Figure 3 The data in (a) is expanded by spherical harmonics and truncated to 60th order. The results are as follows Figure 3 (b) as shown. Figure 3 (c) and Figure 3 (d) The results are obtained by using 300km and 500km Gaussian smoothing filters. Figure 3 The value range of region A1 in (a) is 0.50~1.00cm / yr, Figure 3 The value range of region A2 in (b) is 0.25 to 0.50 cm / yr. Figure 3 The value range of region A3 in (c) is 0.25 to 0.30 cm / yr. Figure 3 The value range of the A4 area in (d) is 0.10-0.15 cm / yr. The red area marked from A1→A2→A3→A4 continuously diffuses the signal to the surrounding area, and the TWS spatial distribution trend in this area continues to decrease. It can be concluded that after the spherical harmonic coefficient expansion and truncation, the signal diffuses to the surrounding area, and then the signal is further attenuated after Gaussian smoothing filtering, and the larger the Gaussian smoothing radius, the greater the signal loss.

[0109] In order to evaluate the leakage error effect of truncation and filtering, the regional average trend changes in the Yangtze River Basin were calculated for comparative analysis, and the trend change loss rate Rate was defined as:

[0110] Rate=(NOAH-Model) / NOAH

[0111] Among them, NOAH is the regional average trend of the GLDAS-NOAH hydrological model, and Model is the regional average trend after filtering of the GLDAS-NOAH hydrological model.

[0112] The regional average trend of the Yangtze River Basin calculated by the initial data is 0.08cm / yr; after spherical harmonic expansion and truncation to the 60th order, the regional average trend without Gaussian smoothing is 0.07cm / yr, which is 12% of the initial data signal loss; after spherical harmonic expansion and truncation to the 60th order, the regional average trend of 300km Gaussian smoothing is 0.06cm / yr, which is 25% of the initial data signal loss; after spherical harmonic expansion and truncation to the 60th order, the regional average trend of 500km Gaussian smoothing is 0.05cm / yr, which is 37% of the initial data signal loss. These data also verify the signal attenuation caused by truncation and filtering. Therefore, restoring the signal loss caused by truncation and filtering of GRACE data is of great significance for studying regional quality changes using GRACE spherical harmonic data inversion.

[0113] 2.1.2 Scale Factor Model

[0114] The present invention uses four GLDAS hydrological models and the scale factor model calculated by CSR-M. The new scale factor model is obtained based on the weighted average correction method of the new scale factor determination coefficient. Figure 4 (a). Figure 4 (b) Figure 4 (d) Figure 4 (e) and Figure 4 As shown in (f), the global hydrological model CLM, MOS, NOAH and VIC data are applied to all grid points of 1°×1° based on formula (5), and the scale factor of the Yangtze River Basin gridding is obtained by deducting the data outside the boundary. In order to compare the scale factor results, the present invention uses the CLM4.0 scale factor provided by the official, and retains the data in the Yangtze River Basin area. Figure 4 (c) shows the scaling factor CLM-SF ( Figure 4 (b))、MOS-SF( Figure 4 (d))、NOAH-SF( Figure 4 (e))) and VIC-SF( Figure 4 (f)) The spatial distribution is similar to that of NSFCM-SF ( Figure 4 (a)) and CLM4.0-SF( Figure 4 (c)) Compared to Figure 4 (b) Figure 4 (d) Figure 4 (e) and Figure 4 (f) The spatial distribution in the lower reaches of the Yangtze River is obviously different. These differences also reflect the differences in the spatial distribution of TWS data from different models, which is caused by the different components considered in the data of different hydrological models.

[0115] Although the scale factors of different models show obvious differences and spatial heterogeneity, the spatial distribution of the grid values ​​of scale factors of different models is basically greater than 1, indicating that the simulation signal in the Yangtze River Basin is mainly controlled by external leakage errors and needs to be corrected using a scale factor greater than 1. The present invention uses the scale factors of different models to calculate the spatial average scale factor in the Yangtze River Basin. The spatial average scale factors of the Yangtze River Basin are all greater than 1, among which CLM4.0-SF has the largest spatial average scale factor (1.42) and MOS-SF has the smallest spatial average scale factor (1.30). The spatial average scale factor reflects the impact on the annual amplitude and long-term trend of the correction results to a certain extent.

[0116] 2.1.3 Iterative recovery method

[0117] The iterative recovery method was first proposed by Chen et al. and has been widely used in GRACE data leakage error correction after improvement. The specific steps of the iterative recovery method are as follows: (1) At each grid point of the 1°×1° grid ( Figure 5 (b)), let the initial mass change rate be the mass change rate of GRACE ( Figure 5 (a)). (2) Forward model quality change rate ( Figure 5 (c)) The mass change rate of the 1°×1° gridded model in step (1) ( Figure 5 (b)), perform spherical harmonic expansion, truncate the spherical harmonic coefficients to 60th order, set the coefficients of 0th and 1st order to 0, and then apply a 300km Gaussian smoothing filter to obtain. (3) At each grid point, the initial GRACE mass change rate ( Figure 5 (a)) and the forward model quality change rate ( Figure 5 (c)) is multiplied by 1.20 and added back to the 1°×1° grid model quality change rate ( Figure 5 (b)), the new model mass change rate is used as step (2), and the process is repeated. After continuous iterations, the consistency between the forward model mass change rate and the initial GRACE mass change rate increases ( Figure 5 (a) with Figure 5 (c)). (4) When the initial GRACE mass change rate ( Figure 5 (a)) and the forward model quality change rate ( Figure 5 (c)) is less than a certain threshold or the iteration is stopped after a certain number of iterations.

[0118] After 60 iterations, the spherical harmonic signal is converted into a grid signal, which is the signal result recovered by the iterative recovery method ( Figure 5 (b)), compare Figure 5 (a) and Figure 5 (c) It is found that there is strong consistency, and the difference between them ( Figure 5(d)) is relatively small (color bar value -0.10 to 0.10 cm / yr). Therefore, it is believed that the result signal restored by the iterative recovery method is credible.

[0119] 2.2 Comparative analysis of the new scale factor determination coefficient weighted average correction method

[0120] 1.2.1 Comparative analysis of spatial distribution trends of water storage changes in the Yangtze River Basin

[0121] The corrected GRACE data time-varying signal is obtained based on the scale factor calculated by the GLDAS hydrological model, the scale factor calculated by the weighted average correction method of the new scale factor determination coefficient, the official CLM4.0 scale factor and the iterative recovery method. Figure 6 As shown. It can be seen that the correction results of the scale factor method and iterative recovery method are better than the CSR-grid inversion results ( Figure 6 (e) The spatial distribution trend signal is significantly enhanced, especially in the middle and lower reaches of the Yangtze River.

[0122] In order to further compare the recovery results of different leakage error correction methods, the present invention uses the CSR-M data correction results as the real signal of the CSR-grid data. Compared with the spatial distribution trend of the recovery results of different leakage error correction methods and the spatial distribution trend of CSR-M data, NSFCM-Trend ( Figure 6 (a))、CLM-Trend( Figure 6 (b))、MOS-Trend( Figure 6 (g))、NOAH-Trend( Figure 6 (h))、VIC-Trend( Figure 6 (i)) and CSR-M-Trend( Figure 6 (d)) are consistent, while CLM4.0-Trend( Figure 6 (c)) Comparison between marked area B and CSR-M-Trend ( Figure 6 (d)) There is an obvious over-recovery phenomenon. FM-Trend( Figure 6 (f)) Comparison between marked area C and CSR-M-Trend ( Figure 6(d)) There is an obvious phenomenon of incomplete signal recovery. The present invention also calculates the regional average trend of the recovery results of different leakage error correction methods in the Yangtze River Basin. From the regional average trend, it can be concluded that CSR-M-Trend is the largest, CSR-grid-Trend is the smallest, and the regional average trend of NSFCM-Trend is greater than the regional average trend of the scale factor correction results calculated by the GLDAS model, but is smaller than the regional average trends of CLM4.0-Trend and FM-Trend. In summary, NSFCM-Trend avoids the phenomenon of incomplete signal recovery and excessive recovery of FM-Trend and CLM4.0-Trend, and the regional average trend of NSFCM-Trend is better than the regional average trend of the scale factor correction results of the GLDAS hydrological model.

[0123] 2.2.2 Comparative analysis of time series of water storage changes in the Yangtze River Basin

[0124] Figure 7 The water storage change time series obtained by restoring the results of different leakage error correction methods in the Yangtze River Basin are shown. Compared with the CSR-grid data and CSR-M data, the restoration results of different leakage error correction methods have better consistency in annual phase, but the annual amplitude has obvious differences, indicating that the scale factor method and iterative recovery method mainly change the annual amplitude of the water storage change time series, while the annual phase of the water storage time series has little effect. In order to further compare the differences in the restoration results of different leakage error correction methods, the present invention calculates the annual amplitude and long-term trend of the water storage change time series in the Yangtze River Basin. Figure 8 shown.

[0125] Figure 8The long-term trend and annual amplitude of water storage changes restored by different leakage error correction methods in the Yangtze River Basin are shown. It can be concluded that the annual amplitude and long-term trend of CSR-grid data are at a low level, and there is signal leakage. After leakage error correction, the annual amplitude and long-term trend are enhanced, and signal leakage is restored to a certain extent. In terms of annual amplitude, the correction results of CSR-M data and FM method are close to CSR-grid data, and the correction results of GLDAS hydrological model scale factor are significantly greater than CSR-grid data. This may be because the TWS simulated by GRACE original observation data is insensitive to fine spatial scales, while CSR-M and CSR-grid are product data, which are mainly constrained by GRACE satellite original observation data, so the annual amplitude restored by FM method is small; most of the changes in TWS simulated by GLDAS hydrological model occur on sub-seasonal and seasonal time scales, so the scale factor calculated by TWS simulated by the model optimizes the annual amplitude. In terms of long-term trend, CSR-M data is the largest (0.41±0.07cm / yr), TWS-CLM4.0 data (0.38±0.11cm / yr), TWS-FM (0.37±0.09cm / yr) and TWS-NSFCM (0.34±0.10cm / yr) are better. The scale factor correction results calculated by the GLDAS hydrological model underestimate the long-term trend, which may be related to the inadequate TWS simulated by the GLDAS hydrological model. The TWS simulated by the GLDAS hydrological model only considers soil moisture, snow water equivalent and total canopy water storage, among which the soil moisture factor is relatively large. The soil moisture depth simulated by MOS (3.50m) and CLM (3.43m) are similar, so the long-term trend is relatively large. The soil moisture depths simulated by VIC and NOAH are similar (1.90m and 2.00m), but the long-term trend corrected by VIC is better, which may be related to the parameters of snow water equivalent and total canopy water storage.

[0126] 3. Application of the new scale factor determination coefficient weighted average correction method

[0127] 3.1 Comparative analysis of water storage changes in the Yangtze River sub-basins

[0128] The NSFCM-SF scale factor combines the advantages of the spatial distribution trend of the scale factor correction results calculated by the GLDAS hydrological model, avoiding the uneven spatial distribution trend of the CLM4.0-SF and FM correction results. In terms of long-term trends, although the NSFCM-SF scale factor correction results are lower than those of CLM4.0-SF and FM, they are better than the long-term trends calculated by the GLDAS hydrological model correction results. In terms of annual amplitude, the NSFCM-SF scale factor correction results are significantly better than the results of FM and CSR-M. Taking all factors into consideration, the present invention mainly discusses the NSFCM-SF scale factor correction results below. The CSR-M, CSR-grid and TWS-NSFCM models are selected to further analyze the correction results of NSFCM-SF in the Yangtze River sub-basins (the upper reaches of the Yangtze River and the middle and lower reaches of the Yangtze River).

[0129] Fig. 9 The water storage time series of the CSR-M, CSR-grid and TWS-NSFCM models in the Yangtze River sub-basin are shown. Fig. 9 It can be seen that the annual amplitude of TWS-NSFCM in the upper reaches of the Yangtze River and the middle and lower reaches of the Yangtze River is significantly better than that of CSR-grid and CSR-M. The time series in the middle and lower reaches of the Yangtze River is more complex than that in the upper reaches of the Yangtze River, and the seasonal TWS changes are greater, which may be related to the seasonal surface water storage changes in the middle and lower reaches of the Yangtze River. The long-term trend of CSR-grid in the upper reaches of the Yangtze River is 0.09cm / yr. After correction by the NSFCM-SF scale factor, the long-term trend of TWS-NSFCM is 0.13cm / yr, and the long-term trend has increased by 44%. The long-term trend of CSR-grid in the middle and lower reaches of the Yangtze River is 0.44cm / yr. After correction by the NSFCM-SF scale factor, the long-term trend of TWS-NSFCM is 0.61cm / yr, and the long-term trend has increased by 39%. It can be concluded that the rising trend of water storage in the Yangtze River Basin is mainly concentrated in the middle and lower reaches of the Yangtze River.

[0130] 3.2 Comparison of precipitation and water storage changes in the Yangtze River Basin

[0131] The present invention compares the time series of the NSFCM-SF scale factor correction results with the monthly precipitation data of the National Meteorological Center. Fig. 9As shown. It can be seen that precipitation in the Yangtze River Basin and its sub-basins is mainly concentrated from May to October. When precipitation increases significantly, the change in TWS-NSFCM water storage increases significantly, and lags behind precipitation data by 1 to 2 months. The change in precipitation data is consistent with the change in TWS-NSFCM water storage, which indicates that precipitation is the main reason for the change in water storage in the Yangtze River Basin. From April 2002 to January 2007, the water reserves and precipitation in the Yangtze River Basin and its sub-basins showed an upward trend. The water reserves in the Yangtze River Basin increased by 0.33 cm / yr, and the precipitation increased by 0.06 cm / yr. The water reserves in the upper reaches of the Yangtze River increased by 0.10 cm / yr, and the precipitation increased by 0.03 cm / yr. The water reserves in the middle and lower reaches of the Yangtze River increased by 0.63 cm / yr, and the precipitation increased by 0.10 cm / yr. It can be concluded that the rising trend of precipitation is mainly concentrated in the middle and lower reaches of the Yangtze River, which is consistent with the feature that the rising trend of water reserves in the Yangtze River Basin is mainly concentrated in the middle and lower reaches of the Yangtze River, further verifying that precipitation is the main reason for the change in water reserves in the Yangtze River Basin. The present invention also calculates the annual precipitation anomaly in the Yangtze River Basin, such as Fig.11 As shown in Figure 2, due to incomplete precipitation data in 2010, 2011 and 2015, the years with incomplete precipitation data are not discussed. Fig.11 It can be seen that the precipitation in most years is close to or lower than the multi-year average, but the precipitation in 2012, 2014 and 2016 is significantly higher than the multi-year average, which may be related to global warming and frequent extreme climate events in recent years (such as El Niño, etc.). Fig.10 In years with more precipitation, the interannual variation of TWS-NSFCM water reserves increases significantly. Therefore, the interannual variation characteristics of TWS-NSFCM water reserves in the Yangtze River Basin are consistent with the interannual variation law of precipitation.

[0132] Although the present invention has been disclosed as above in the form of a preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art may make possible changes and modifications to the technical solution of the present invention by using the methods and technical contents disclosed above without departing from the spirit and scope of the present invention. Therefore, any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall fall within the protection scope of the technical solution of the present invention.

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

Claims

1. A method for improving the accuracy of GRACE terrestrial water storage changes, characterized in that: include: Obtain the terrestrial water storage Δh(θ,λ) from the hydrological model and determine the spherical harmonic coefficient cosine term C corresponding to Δh(θ,λ) lm and the sine term S lm ; Among them, θ represents co-latitude and λ represents longitude; C lm and S lm Perform de-striping filtering and Gaussian smoothing filtering to obtain the filtered spherical harmonic coefficient cosine term and the sine term based on and Get the filtered terrestrial water storage Δη(θ,λ); Determine the time series data Δh corresponding to Δh(θ,λ) in the time series t i (t), Δη(θ,λ) corresponds to the time series data Δη at time t i (t); According to Δh i (t) and Δη i (t), solve to get the weighted scale factor of the coefficient of determination; The leakage error in the process of inverting the terrestrial water storage changes from GRACE data is corrected according to the weighted scaling factor of the coefficient of determination, and the real hydrological signal is restored to improve the accuracy of the GRACE terrestrial water storage changes. The spherical harmonic coefficient cosine term C corresponding to Δh(θ,λ) is determined by the following formula (1): lm and the sine term S lm : Among them, ρ w represents the density of water, a represents the average radius of the earth, ρ ave represents the average density of the Earth, represents the normalized associated Legendre function, l represents the order of the spherical harmonic coefficient solution inversion, m represents the number of times the spherical harmonic coefficient solution inversion is performed, and k l represents the load Love number of order l; C lm and S lm Perform de-striping filtering and Gaussian smoothing filtering to obtain the filtered spherical harmonic coefficient cosine term and the sine term include: According to the following equations (2) and (3), C lm and S lm Perform de-striping filtering to obtain the cosine term of the spherical harmonic coefficients after de-striping filtering and the sine term Where w is the sliding window width, e is the natural logarithm, max() means taking the maximum value, Q is the destriping filter coefficient; for destriping of odd items, n is the odd item in the lth order, and the same is true for even items; According to the following equations (4) and (5), and Perform Gaussian smoothing filtering to obtain the cosine term of the spherical harmonic coefficients after Gaussian filtering and the sine term Among them, r represents the Gaussian smoothing radius, W0 represents the initial order Gaussian smoothing weight, and W l-1 , W l and W l+1 They represent the Gaussian smoothing weights of order l-1, l and l+1 respectively, and l is 0, 1, 2, 3…60.

2. The method for improving the accuracy of GRACE terrestrial water storage changes according to claim 1, characterized in that: The filtered terrestrial water storage Δη(θ,λ) is obtained by the following formula (6):

3. The method for improving the accuracy of GRACE terrestrial water storage changes according to claim 1, characterized in that: According to Δh i (t) and Δη i (t), the coefficient of determination weighted scaling factor is obtained by solving, including: According to Δh i (t) and Δη i (t), based on the least squares principle, the scale factor K is calculated; According to Δh i (t), and the coefficient of determination R is obtained by solving 2 ; According to the scale factor K and the determination coefficient R 2 , the determination coefficient weighted scaling factor NSFCM-SF is obtained by solving the equation.

4. The method for improving the accuracy of GRACE terrestrial water storage changes according to claim 3, characterized in that: The scale factor K is obtained by solving the following equation (7): Among them, N represents the number of time series, and i represents the i-th time series.

5. The method for improving the accuracy of GRACE terrestrial water storage changes according to claim 4, characterized in that: According to Δh i (t), and the coefficient of determination R is obtained by solving 2 ,include: According to Δh i (t), determine the fitted terrestrial water storage and the average terrestrial water storage The determination coefficient R is obtained by solving the following formula (8): 2 :

6. The method for improving the accuracy of GRACE terrestrial water storage changes according to claim 5, characterized in that: The coefficient of determination weighted scaling factor NSFCM-SF is obtained by solving the following formula (9): Where n represents the number of scale factor models, j represents the jth scale factor model, and K j represents the scale factor corresponding to the j-th scale factor model, Represents the determination coefficient corresponding to the j-th scale factor model.

Images