Adaptive constraint method for determining global mass change by using gravity satellite data

Abstract

The invention relates to a self-adaptive constraint method for determining global mass change by using gravity satellite data, which comprises the following steps of: firstly, obtaining a global mass change signal containing noise by using gravity satellite observation data and a background force model, and constructing an initial constraint matrix by using a zero-order Tikhonov regularization technology to obtain an initial constraint solution of the global mass change; secondly, considering land and sea signal leakage, respectively calculating a global quality change strong signal and a global quality change weak signal by utilizing the strong signal constraint matrix and the weak signal constraint matrix, adding the strong signal and the weak signal to obtain a global quality change complete signal, and finally, considering land signal leakage to obtain a global quality change complete signal. And constructing a leakage error correction constraint matrix by using the obtained root mean square value (RMS) of the global quality change complete signal as prior information so as to solve the global quality change signal after leakage error correction.

Description

Technical Field

[0001] The present invention relates to the field of satellite gravity measurement applications, and more specifically to an adaptive constraint method for determining global mass changes using gravity satellite data. Background Art

[0002] The mass migration changes of the Earth system monitored by gravity satellites are usually expressed as monthly mean time-varying gravity field models or Mascon models. Affected by factors such as instrument noise levels, background model accuracy, and non-modeling errors, the calculated global mass changes have serious "strip" noise in the north-south direction in the airspace, and the true geophysical signals are also submerged by the strip noise. Common noise reduction methods are the Gaussian smoothing method and the prior information constraint method. The Gaussian smoothing method is to downweight the high-order spherical harmonic potential coefficients with large noise, reduce the error of the short-wave components, and achieve the purpose of noise reduction. The constraint method based on prior information is mainly used to solve the Mascon model. Currently widely used constraint methods include regularization smoothing constraints, spatio-temporal constraint equations, and prior physical information constraints, etc. Existing constraint methods will lead to problems such as unclear physical meaning of prior information or distorted constraint information. The calculated global mass changes are contaminated by prior information and cannot accurately and truly reflect the global mass changes of the current month.

[0003] Since the signal intensities of global surface mass changes are different, for example, there are strong periodic change signals in the Amazon Basin, and strong long-term change signals in Greenland and Antarctica. Relatively speaking, the long-term change signals caused by groundwater overexploitation in the North China Plain of China, the Mississippi River Basin of the United States, and northern India are relatively weak. Using a unified constraint parameter to solve will cause the following problems: If the constraint parameter is weak, although the weak signals of global mass changes can be retained, they will be affected by strong strip errors. If the constraint parameter is strong, although the influence of strip errors can be weakened, the weak signals of mass changes will be overconstrained.

[0004] When calculating the global mass change signal, in order to eliminate the land-sea signal leakage error, the prior art retains the RMS of the land signal unchanged and sets the maximum value of the RMS of the ocean signal to a unified value, believing that the maximum RMS of the mass change signal in the sea area near the areas with large mass changes in the Greenland and Antarctic ice sheets does not exceed 4 cm. Therefore, the maximum value of the RMS of the ocean signal for each month except for earthquake areas and enclosed seas is uniformly set to 4 cm. However, in the long term, the ocean mass is not fixed. It is inappropriate to set the maximum value of the RMS of the ocean mass change signal for each month to a unified value, and the change of the ocean mass signal needs to be considered. Summary of the Invention

[0005] The purpose of the present invention is to provide an adaptive constraint method for determining global mass changes using gravity satellite data to solve the problems raised in the above background art.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] An adaptive constraint method for determining global mass change using gravity satellite data, comprising the following steps:

[0008] S1: Construct an initial constraint matrix: Obtain a noisy global mass change signal using gravity satellite observation data and a background force model, and then use the zero-order Tikhonov regularization technique to construct an initial constraint matrix to obtain an initial constraint solution for global mass change;

[0009] S2: Land-sea signal leakage constraint: Calculate the strong signal of the global mass change signal and the weak signal of the global mass change signal using the strong signal constraint matrix and the weak signal constraint matrix respectively, and add the strong signal and the weak signal to obtain a complete global mass change signal; wherein, the initial constraint solution for global mass change calculated in step S1 is used as prior information and combined with land-sea information to construct a strong signal constraint matrix, and the construction of the weak signal constraint matrix is as follows:

[0010] (1) Use the method of step S1 to calculate the initial constraint solution for global mass change each month, and fit the global mass change signal using the following formula:

[0011]

[0012] where y(t) is the global mass change signal corresponding to time t, C is the fitting constant term, B is the long-term trend term, A is the interannual amplitude, is the interannual phase; ε is the residual;

[0013] (2) Use the calculated fitting constant term C, long-term trend term B, interannual amplitude A, and interannual phase to fit the fitted value of the global mass change each month, then subtract the fitted value from the initial constraint solution, and use the difference as prior information and combine it with land-sea information to construct a weak signal constraint matrix;

[0014] S3: Land signal leakage constraint: Use the root mean square value (RMS) of the complete global mass change signal obtained in step S2 as prior information to construct a constraint matrix for processing the leakage of the global mass change signal, so as to solve the global mass change signal after leakage error correction.

[0015] A further technical solution of the present application: The initial constraint matrix in step S1 is set as the identity matrix.

[0016] A further technical solution of the present application: In step S2, the method for constructing the land-sea information is to keep the root mean square value (RMS) of the land signal unchanged, and set the maximum value of the root mean square value (RMS) of the ocean signal as a certain threshold, and the thresholds for different months are different.

[0017] A further technical solution of the present application: The threshold is determined by 85% of the overall root mean square value (RMS) of the monthly ocean signal.

[0018] Compared with the prior art, the present invention can achieve at least the following beneficial effects:

[0019] 1. The zero-order Tikhonov regularization technique is used to construct an initial constraint matrix to solve the initial constraint value of the global mass change. This initial constraint value contains all the spectral information of the mass change, has a clear physical meaning, and does not contain or depend on external geophysical signals.

[0020] 2. The strong signal constraint matrix and the weak signal constraint matrix are used to calculate the strong signal of the global mass change and the weak signal of the global mass change respectively. Adding the strong signal and the weak signal can obtain the complete signal of the global mass change, avoiding the influence of strip noise caused by insufficient constraint of the strong signal of the mass change or over-constraint of the weak signal of the mass change due to the use of unified constraint information.

[0021] 3. When dealing with the leakage of land and sea signals, the long-term change of the ocean mass is accurately considered. The maximum RMS value of the ocean signal is set as a unified threshold, but the threshold changes with the overall intensity of the ocean signal every month. The dynamic determination of the threshold can adapt to the long-term change of the ocean mass. The land signals in the strong signal constraint matrix and the weak signal constraint matrix incorporate the signals of the global mass change of the current month, that is, the strong signal constraint matrix and the weak signal constraint matrix are time-varying and can be optimally matched with the global mass change signal of the current month.

[0022] 4. The adaptive constraint model also considers the problem of land signal leakage. Taking the RMS of the complete signal of the global mass change as prior information, a constraint matrix for dealing with the leakage of the global mass change signal is constructed to solve the global mass change signal after leakage error correction, solving the problem that it is difficult to correct the land signal leakage. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 is the construction process of the adaptive constraint model;

[0024] Figure 2 is the comparison between the global mass change result (a) under the zero-order Tikhonov regularization constraint and the spherical harmonic result (b), the global mass change in April 2015 in the form of equivalent water height (unit: mm), and the spherical harmonic result is the 96th order CSR-RL06 spherical harmonic model with 300 km Gaussian filtering;

[0025] Figure 3Comparison of the long-term trends and annual amplitudes between the global mass change results (a and b) with zero-order Tikhonov regularization constraints and the CSR-RL06 spherical harmonic results (c and d) of degree 96 from April 2002 to August 2016;

[0026] Figure 4 For the statistical analysis of the ocean signal, (a) shows the percentages of the global ocean signal with RMS maxima not exceeding 80mm, 70mm, 60mm, 50mm, 40mm, 30mm, 20mm, 10mm in the overall global ocean signal on a monthly basis, and (b) shows the RMS values corresponding to the 95%, 90%, 85%, 80%, 75%, 70%, 65%, 60% percentiles of the global ocean signal on a monthly basis. The values in parentheses in the legend represent the mean of all months; Figure 5 For the RMS of the strong and weak signals of the global mass change and their constrained results, taking April 2015 as an example, where (a) represents the prior RMS for constructing the strong signal constraint matrix, (a) represents the solution of the strong signal constraint, (c) represents the prior RMS for constructing the weak signal constraint matrix, and (d) represents the solution of the weak signal constraint. Note that the color scale ranges are inconsistent;

[0027] Figure 6 For the global mass change signal before and after land signal constraint (taking April 2015 as an example), (a) is the global mass change before signal leakage correction (i.e., the sum of the strong and weak signals in the second step), and (b) is the global mass change after signal leakage correction (the final global mass change result calculated in this paper);

[0028] Figure 7 For the comparison of geoid degree errors;

[0029] Figure 8 For the long-term trend of the global mass change signal from April 2002 to August 2016;

[0030] Figure 9 For the amplitude of the global mass change signal from April 2002 to August 2016. Specific implementation manners

[0031] The following combines the accompanying drawings and embodiments to further describe in detail the specific implementation manners of the present invention. The following embodiments are used to illustrate the present invention but not to limit the scope of the present invention.

[0032] The specific steps of this method are as follows:

[0033] The first step: Construct the initial constraint matrix

[0034] The global mass change signal with noise is obtained by using gravity satellite observation data and background force models. The zero-order Tikhonov regularization technique is adopted to construct the initial constraint matrix to solve the initial constraint value of global mass change. In order to better constrain the mass change in the polar ice sheet regions, the post-glacial rebound model is used to subtract the long-term trend term of post-glacial rebound. The initial constraint matrix is set as the identity matrix, and the initial constraint result of global mass change is calculated based on the initial constraint matrix. The optimal regularization parameter is α1, and this result is less affected by strip errors. As Figure 2 shown, the result is comparable to the spherical harmonic coefficient model after 300 km Gaussian filtering, and there are obvious signal leakage problems in both cases.

[0035] Step 2: Constraining the land-sea signal leakage

[0036] In order to better weaken the influence of strip errors and retain the weak signal of global mass change as much as possible, we calculate the strong signal of global mass change and the weak signal of global mass change by using the strong signal constraint matrix and the weak signal constraint matrix respectively. Adding the strong signal and the weak signal together can obtain the complete signal of global mass change. The specific steps are as follows:

[0037] (1) Calculate the fitting value of global mass change per month. The initial constraint solution of global mass change for each month is calculated by using the method of the first step, and the global mass change signal is fitted by the following formula:

[0038]

[0039] where y(t) is the global mass change signal corresponding to time t, C is the fitting constant term, B is the long-term trend term, A is the interannual amplitude, is the interannual phase; ε is the residual.

[0040] The long-term trend and interannual amplitude of the initial constraint solution of global mass change from April 2002 to August 2016 and the 96th-order CSR-RL06 spherical harmonic solution of the Center for Space Research, Texas State University are as Figure 3 shown. The two are generally in good agreement, but there are also some differences in some regions such as Greenland due to differences in data processing details.

[0041] (2) Construct the strong signal constraint matrix and the weak signal constraint matrix of global mass change per month. The initial constraint solution of global mass change calculated in the first step is used as prior information to construct the strong signal constraint matrix. The fitting constant term C, long-term trend term B, interannual amplitude A and interannual phase calculated are used to fit the fitting value of global mass change per month, and then the fitting value is subtracted from the initial constraint solution, and the difference is used as prior information to construct the weak signal constraint matrix.

[0042] When constructing the strong signal constraint matrix and the weak signal constraint matrix, in order to eliminate the land-sea signal leakage error, the land-sea constraint equation is used for constraint. In the long term, the ocean mass is not fixed, and the change of the ocean mass signal needs to be considered. We conduct a statistical analysis on the ocean part of the initial constraint solution of the global mass change. Figure 4 a shows the percentage of the ocean signal RMS not exceeding a certain threshold. For example, the mean value of the signals with the ocean signal RMS not exceeding 40 mm accounts for about 85% of the total ocean signals, but the percentages in different months are in a fluctuating state. At the same time, we statistically analyze the magnitude of the signal RMS corresponding to the ocean signal not exceeding a certain percentage value (such as Figure 5 b), the RMS mean value of the signals accounting for 85% of the overall signals is about 40 mm, and the RMS values in different months under the same percentage are different. Since our reference average epoch is 2010.0, therefore, the signal intensity is greater for the time farther from the reference epoch, which also reflects that the global ocean mass is changing in the long term. Through calculation and analysis, we set the maximum value of the ocean signal RMS as a certain threshold, but the thresholds in different months are not fixed and are determined by 85% of the overall root mean square value (RMS) of the ocean signal in that month (as shown in Figure 4 b). At the same time, in order to retain some marine geophysical signals, we do not process the marine seismic regions and (semi-)enclosed ocean regions, and the RMS values in these regions remain unchanged.

[0043] (3) Adding the said strong signal and the said weak signal can obtain the complete signal of the global mass change. Examples of the prior RMS of the global mass change and its constraint results for constructing the strong signal constraint matrix and the weak signal constraint matrix are as shown in Figure 5 shown.

[0044] The third step: Constraint of land signal leakage

[0045] When using the complete signal of the global mass change obtained in the second step, we perform constraint correction on the leakage of the land-sea signal, but the signal leakage in the land area is still relatively obvious. Through experimental analysis, it is found that the land signal leakage correction can still be carried out through constraint correction. Taking the RMS of the complete signal of the global mass change obtained in the second step as the prior information, a leakage error correction constraint matrix is constructed to solve the global mass change signal after leakage error correction. At this time, no land-sea constraint processing is performed on the complete signal of the global mass change. The result obtained by using the constructed leakage error correction constraint matrix is the global mass change signal after leakage error correction.

[0046] The global mass change before and after leakage error correction is as shown in Figure 6 shown, Figure 6 (a) is the sum of the strong signal and the weak signal in the second step, Figure 6 (b) is with Figure 6(a) is the solution result of the constraint matrix constructed for prior information. It can be seen that in regions such as Greenland, Antarctica, North America, Tianshan region, Pamir Plateau, North China Plain of China, and the Caspian Sea, after leakage error correction, the signal leakage in these regions has been significantly improved, and the signal intensity has also increased.

[0047] By constructing the constraint matrix through the above steps to solve for global mass changes, since the prior information required for constructing the constraint matrix in each step does not include external geophysical signals, the global RMS is provided by the observation signals of the gravity satellites themselves, and the constraint matrix for each month is constructed from the satellite observation signals of that month. Therefore, when solving for global mass changes by constructing the constraint matrix using this method, the constraint matrix can be appropriately constructed adaptively according to the intensity and spatial distribution of the time-varying signals of the current month. Through Figure 1 The global mass changes solved by the steps shown can separate the strong and weak signals of global mass changes, and the best regularization parameter can be adopted to avoid over-constraining or under-constraining of the signals, and the problem of signal leakage on land has also been well solved.

[0048] Embodiment

[0049] Based on the above adaptive constraint method, namely PMA-final, the monthly global mass change signals from April 2002 to August 2016 were calculated using Gravity Recovery and Climate Experiment (GRACE) Level-1B data. The calculated global mass changes were compared and analyzed with those calculated by other models. The spatial domain signals were transformed into spectral domain spherical harmonic coefficients of 120 orders through spherical harmonic analysis, and the geoid order error was calculated. Taking April 2002 as an example, the geoid order errors of each model are as Figure 7 shown. The models participating in the comparison are the Mascon model released by the Center for Space Research (CSR) of Texas State University, the Jet Propulsion Laboratory (JPL), the Goddard Flight Control Center (GSFC), and the 96-order RL06 time-varying gravity field model released by the Center for Space Research (CSR) of Texas State University.

[0050] Comparative Example 1

[0051] Using the CSR-RL06 model, other steps are the same as in the embodiment.

[0052] Comparative Example 2

[0053] Using the CSR MASCON RL06.2 model, other steps are the same as in the embodiment

[0054] Comparative Example 3

[0055] Using the JPL MASCON RL06.1M model, other steps are the same as those in the embodiment.

[0056] Comparative Example 4

[0057] Using the GSFC MASCON model, other steps are the same as those in the embodiment.

[0058] Comparative Example 5

[0059] The global mass change calculated based on the initial constraint matrix, i.e., PMA-Tikh, other steps are the same as those in the embodiment.

[0060] Comparative Example 6

[0061] The global mass change calculated based on the initial constraint matrix and the land-sea signal leakage constraint, i.e., using PMA-sum, other steps are the same as those in the embodiment.

[0062] From Figure 7 It can be seen that compared with the unfiltered spherical harmonic coefficients, the results calculated by the constrained Mascon and the above-mentioned constraint methods have significantly improved noise at high orders; the global mass change calculated using the zero-order Tikhonov constraint decays rapidly after the 10th order, and the results using the adaptive constraint method have good consistency with the Mascon models published by CSR, JPL, and GSFC in the entire spectral domain. In particular, the final global mass change after leakage error correction is basically consistent with the CSR Mascon. The long-term trend of the global mass change from April 2002 to August 2016 calculated by the above-mentioned adaptive constraint method (as Figure 8 shown) and the interannual amplitude of the global mass change (as Figure 9 shown) are compared with the results of CSR, JPL, and GSFC Mascon during the same period. In terms of the overall signal, the results calculated by several methods have good consistency. Compared with the long-term trend of the 96th-order spherical harmonic coefficient model in Figure 3 (c), the model calculated by the present method has smaller signal leakage errors with the Mascon model and is basically not affected by striping errors. Of course, due to the inconsistent specific parameter settings for calculating the global mass change by each institution, there are still differences in local areas, such as the inland area of Greenland, some areas adjacent to the Caspian Sea, and the Sumatra earthquake sea area. At the same time, we compared the interannual amplitudes of the global mass change signals of each model (as Figure 9 and Figure 3 shown), and the signal amplitudes of each model are generally in good consistency.

[0063] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, and improvements made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. An adaptive constraint method for determining global mass change using gravity satellite data, characterized in that, It includes the following steps: S1: Construct the initial constraint matrix: Obtain the global mass change signal with noise by using the gravity satellite observation data and the background force model, and then adopt the zero-order Tikhonov regularization technique to construct the initial constraint matrix to obtain the initial constraint solution of the global mass change; S2: Leakage constraint of land-sea signal: Calculate the strong signal of the global mass change signal and the weak signal of the global mass change signal by using the strong signal constraint matrix and the weak signal constraint matrix respectively. Adding the strong signal and the weak signal can obtain the complete signal of the global mass change. Among them, the initial constraint solution of the global mass change calculated in step S1 is used as prior information and combined with the land-sea information to construct the strong signal constraint matrix. The construction of the weak signal constraint matrix is as follows: (1) Use the method in step S1 to calculate the initial constraint solution of the global mass change for each month, and fit the global mass change signal with the following formula: where y(t) is the global mass change signal corresponding to time t, C is the fitting constant term, B is the long-term trend term, A is the interannual amplitude, is the interannual phase; ε is the residual; (2) Using the calculated fitting constant term C, long-term trend term B, interannual amplitude A, and interannual phase to fit the fitting values of the global mass change for each month, then subtracting the fitting values from the initial constrained solution, using the difference as prior information and combining land-sea information to construct a weak signal constraint matrix; S3: Leakage constraint of land signal: Use the root mean square value (RMS) of the complete global mass change signal obtained in step S2 as prior information to construct the leakage error correction constraint matrix to solve the global mass change signal after leakage error correction.

2. The method according to claim 1, wherein The initial constraint matrix in step S1 is set as the identity matrix.

3. The method according to claim 1, characterized in that, In step S2, the method for constructing the land-sea information is to keep the root mean square value (RMS) of the land signal unchanged, and set the maximum value of the root mean square value (RMS) of the ocean signal as a certain threshold, and the thresholds for different months are different.

4. The method according to claim 3, wherein The threshold is determined by 85% of the overall root mean square value (RMS) of the ocean signal in that month.

Images