Gravitational field inversion method considering unsteady noise of satellite gravity observation value

Abstract

The invention discloses a gravitational field inversion method considering unsteady noise of a satellite gravity observation value. The gravitational field inversion method comprises the following steps: A, acquiring original observation data of a gravity satellite; b, setting observation value noise as white noise, and constructing a random model considering non-stationary noise; c, performing parameter estimation to obtain a post-test residual error of an observation value; d, estimating the noise auto-covariance and precision factor matrix of the observation value based on the post-test residual error of the observation value; e, repeating the steps C and D for iteration to realize refinement of the gravitational field parameter and observation value random model; and F, evaluating the time-varying gravity field model in the spectral domain. According to the invention, the defects in the prior art can be improved, and the precision and stability of gravitational field inversion are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite gravity measurement, and in particular to a gravity field inversion method considering the non-steady noise of satellite gravity observation values. Background Art

[0002] The time-varying gravity field of the Earth reflects the mass distribution and movement inside and on the surface of the Earth, and is an important constraint for studying the mass exchange and its spatio-temporal changes among the land, atmosphere and ocean. Since the launch of the GRACE satellite in 2002, its data has played a milestone role in the research of the global hydrological cycle, glacier mass balance and co-seismic effects of earthquakes. However, limited by the colored noise of the observation values, there is still an order-of-magnitude gap between the accuracy of the existing gravity field models and the design indicators, which restricts their in-depth application in earth sciences.

[0003] The GRACE observation errors include systematic errors and random errors, and their processing strategies are significantly different. The systematic errors are generally weakened by prior models (such as the receiver antenna phase center error of GPS observations) or parameterization methods (such as the scale and bias of accelerometer observations), while the random errors need to be considered by establishing an accurate random model. The inversion of the GRACE time-varying gravity field mainly involves four types of observations. First, the K-band ranging data is the main observation for gravity field inversion, and its ranging accuracy at the micron level provides high-sensitivity information for gravity field inversion. It should be noted that the original phase observations have the problem of integer ambiguity and are extremely sensitive to systematic errors. Therefore, in actual processing, differential observations (ranging rate or acceleration measurements) are generally used to eliminate the ambiguity and reduce the influence of systematic errors, but the differential operation will introduce the temporal correlation of the observation noise. Second, the results of kinematic orbit determination are usually used as input observations. The incomplete fixation of the GPS phase ambiguity parameters will result in significant correlations between the orbit parameters of adjacent epochs. In addition, the incompleteness of the dynamic model will further cause the cross-epoch correlation of the observation residuals. The accelerometer data is used for non-conservative force modeling, and its random error characteristics directly affect the accuracy of gravity field recovery. The satellite attitude observation data, by inverting the satellite attitude parameters from the star sensor data, its error mainly affects the accurate construction of the observation geometric relationship.

[0004] Traditional gravity field inversion methods (such as the dynamic method and the short arc method) mostly assume that the observation noise is stationary white noise and use a fixed variance matrix to determine weights. However, the actual observation noise shows significant frequency correlation and time-varying non-stationarity. To account for the colored noise in the observations, the most direct method is to construct the noise variance-covariance matrix of the observations for weight determination. Since the gravity field solution involves a large number of observations, the dimension of the variance-covariance matrix is very high, so this method has a large computational amount; another method is to determine weights based on the noise power spectrum of the observations. Since this method uses the fast Fourier transform in the calculation, the solution speed is relatively fast and has been successfully applied in the mean acceleration method and the dynamic method and can significantly improve the model quality. However, the above methods require that the observation noise has stationary characteristics, that is, the random characteristics of the noise do not change with time. In reality, affected by factors such as temperature changes and star tracker blinding, the accuracy of K-band observations is dynamically changing; for kinematic orbits, as the geometric configuration of visible GPS satellites changes, its accuracy is also dynamically changing. Therefore, the observation noise usually has non-stationary characteristics. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a gravity field inversion method that takes into account the non-steady noise of satellite gravity observations, which can solve the deficiencies of the prior art and improve the accuracy and stability of gravity field inversion.

[0006] To solve the above technical problem, the technical invention adopted by the present invention is as follows.

[0007] A gravity field inversion method that takes into account the non-steady noise of satellite gravity observations includes the following steps:

[0008] A. Obtain the original observation data of the gravity satellite;

[0009] B. Assume that the observation noise is white noise and construct a stochastic model that takes into account non-stationary noise;

[0010] C. Perform parameter estimation to obtain the a posteriori residuals of the observations;

[0011] D. Estimate the autocovariance of the observation noise and the precision factor matrix based on the a posteriori residuals of the observations;

[0012] E. Repeat steps C and D for iteration to refine the gravity field parameters and the stochastic model of the observations;

[0013] F. Evaluate the time-varying gravity field model in the spectral domain.

[0014] Preferably, the original observation data includes K-band ranging observations, GPS orbit observations, prior accuracy information, AOD data, and EOP data.

[0015] Preferably, the stochastic model that takes into account non-stationary noise is Among them, C is the variance-covariance matrix of the steady-state noise, which characterizes the correlation between the observed values. Λ is the dilution of precision matrix, and the diagonal elements thereof characterize the non-steady-state precision characteristics of the observed values at each epoch.

[0016] Preferably, in step C, the process of parameter estimation is as follows where is the parameter estimated value, A is the design matrix, the dimension of A is m×n, m is the number of observations, and n is the number of parameters; A T is the transpose of the design matrix, y is the observation vector, that is, the actually observed data, (A T A) -1 A T is the pseudo-inverse matrix, which is used for least squares solution;

[0017] The a posteriori residual refers to the difference between the predicted value obtained by estimating the parameters and the actual observed value, and its calculation formula is where v is the a posteriori residual vector, is the prediction vector, that is, the expected observed value calculated according to the estimated parameters.

[0018] Preferably, the estimation method of the autocovariance c of the observed value noise is as follows

[0019] where c k is the covariance between epochs j and j±k, n j is the a posteriori residual of the observed value at epoch j, N k is the number of epoch pairs used to calculate c k N a is the maximum delay epoch number; perform periodic extension on the autocovariance c of the observed value noise, then the variance-covariance matrix C of the steady-state noise is a circulant matrix, and the product of the variance-covariance matrix C of the steady-state noise and the vector b is p = Cb. Calculate the discrete Fourier transform b s = DFT(b), calculate the power spectrum c s = DFT(c), calculate the frequency domain filtering result p s = c s ·b s , and perform inverse transformation to restore the time domain data p = IDFT(p s ).

[0020] Preferably, the estimation method of the dilution of precision matrix is as follows

[0021] Use to remove the correlation between the residuals, where n is the residual vector of the original observed values, with a dimension of N×1, N is the number of observed values, which represents the difference between the observed values and the model predicted values, and includes systematic errors and random noise, n dcis the decorrelated residual vector, with a dimension of N×1;

[0022] Use the Huber M estimator Calculate the dilution of precision corresponding to the observations for each epoch, where k is the adjustment constant of the Huber estimator;

[0023] Set k = 2 based on empirical values and calculate the a posteriori residuals is the decorrelated residual vector, with a dimension of N×1, and N is the number of observations.

[0024] Preferably, in step F, a spectral domain precision evaluation method is used to quantitatively analyze the gravity field model. The obtained gravity field model is compared with the spherical harmonic coefficients of the static gravity field, and the differences between the coefficients of the gravity field model are calculated; then the order variance and cumulative order variance are calculated to evaluate the precision index of the model.

[0025] The beneficial effects brought by the above technical invention are as follows:

[0026] (1) Compared with the traditional gravity field inversion method, the time-varying gravity field inversion method considering the non-steady noise of observations realizes the explicit separation of steady and non-steady noises by constructing the noise variance-covariance matrix of observations, and more realistically models the observation noise.

[0027] (2) For large observation residuals, the dilution of precision is dynamically adjusted through different robust estimators, significantly enhancing the overall robust performance of the algorithm.

[0028] (3) The present invention is applicable to the processing of gravity satellite data such as GRACE (Gravity Recovery and Climate Experiment) and GRACE-FO (Follow-On), and can effectively improve the precision of gravity field inversion.

[0029] (4) The present invention can be used to monitor geophysical phenomena such as large-scale water body migration, crustal deformation, and sea level change, providing high-precision gravity field data support for research in related fields. Brief Description of the Drawings

[0030] Figure 1 is the flow chart of the present invention.

[0031] Figure 2 is the high-order geoid error of the gravity field model solved by different schemes.

[0032] Figure 3 is the calculation process of power spectrum weight determination. Detailed Embodiments

[0033] A gravity field inversion method considering the non-steady noise of satellite gravity observations, comprising the following steps:

[0034] A. Obtain the original observation data of the gravity satellite. The original observation data includes K-band ranging observations, GPS orbit observations, prior accuracy information, AOD data, and EOP data. The observation values and orbit data of the gravity satellite mission are GRACE / GFO Level-1B data, including KBR, ACC, SCA, GNV data. In addition, it also includes AOD RL06 version data. Convert these data into the required format for calculation, including information such as time and observation values.

[0035] B. Assume that the observation noise is white noise and construct a stochastic model considering non-stationary noise where C is the variance-covariance matrix of the steady-state noise, characterizing the correlation between observation values, and Λ is the accuracy factor matrix, whose diagonal elements characterize the non-steady accuracy characteristics of the observation values at each epoch.

[0036] C. Perform parameter estimation, where is the parameter estimation value, A is the design matrix, the dimension of A is m×n, m is the number of observations, and n is the number of parameters; A T is the transpose of the design matrix, y is the observation vector, that is, the actually observed data, (A T A) -1 A T is the pseudo-inverse matrix, used for least squares solution;

[0037] The a posteriori residual refers to the difference between the predicted value obtained by estimating the parameters and the actual observation value, and its calculation formula is where v is the a posteriori residual vector, is the prediction vector, that is, the expected observation value calculated according to the estimated parameters.

[0038] D. Estimate the autocovariance of the observation noise and the accuracy factor matrix based on the a posteriori residual of the observation value.

[0039] The estimation method of the autocovariance c of the observation noise is

[0040] where c k is the covariance between epochs j and j±k, n j is the a posteriori residual of the observation value at epoch j, N k is the number of epoch pairs used to calculate c k is the number of epoch pairs used to calculate c ais the maximum number of delay epochs; the autocovariance c of the observation noise is periodically extended, then the variance-covariance matrix C of the steady-state noise is a circulant matrix, and the product of the variance-covariance matrix C of the steady-state noise and the vector b is p = Cb, and the discrete Fourier transform b of the data is calculated s = DFT(b), calculate the power spectrum c of the noise s = DFT(c), calculate the frequency-domain filtering result p s = c s ·b s , and the inverse transform is used to recover the time-domain data p = IDFT(p s ).

[0041] The estimation method of the dilution of precision matrix is as follows

[0042] Use to remove the correlation between residuals, where n is the residual vector of the original observations, with a dimension of N×1, and N is the number of observations, which represents the difference between the observations and the model predicted values, including systematic errors and random noise, n dc is the decorrelated residual vector, with a dimension of N×1;

[0043] Use the Huber M estimator to calculate the dilution of precision corresponding to each epoch observation value, where k is the adjustment constant of the Huber estimator;

[0044] Based on the empirical value, set k = 2 and calculate the a posteriori residuals is the decorrelated residual vector, with a dimension of N×1, and N is the number of observations

[0045] E. Repeat steps C and D for iteration to refine the gravity field parameters and the stochastic model of the observations

[0046] F. Adopt the spectral-domain accuracy evaluation method to quantitatively analyze the gravity field model, compare the solved gravity field model with the spherical harmonic coefficients of the static gravity field, and calculate the differences between the gravity field model coefficients; then calculate the order variance and the cumulative order variance to evaluate the accuracy index of the model

[0047] Based on the GRACE satellite gravity observation data, two groups of time-varying gravity field inversion comparison experiments were carried out in this study: the benchmark group (S1) adopted the traditional stationary noise assumption, and its dilution of precision matrix was simplified to the identity matrix; the improved group (S2) constructed a time-varying dilution of precision matrix and realized the non-stationary noise modeling by fusing the prior precision information of the observations and the a posteriori residual statistics Figure 2The comparison results of the order errors of the gravity field models under two noise assumptions are shown with GOCO06S as the reference background field. Among them, the blue line represents the order error distribution of the S1 method, and the red line corresponds to the improvement effect of the S2 method. Experimental data shows that in terms of the cumulative error of the geoid at the 96th order, the S1 method is 10.819 mm, and the S2 method is optimized to 8.007 mm, with a noise suppression amplitude of 26%. This fully verifies the role of the non-steady noise modeling technology in improving the inversion accuracy of the gravity field. Especially in the region of high-order harmonic coefficients (>40th order), the S2 method shows a more significant error suppression effect.

[0048] The present invention is applied to the inversion of the time-varying gravity field of GRACE, GFO or low-low tracking mode gravity satellites, and is used for monitoring the Earth's mass migration, analyzing hydrological changes or evaluating glacier ablation.

[0049] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by the terms "longitudinal", "transverse", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation of the present invention.

[0050] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art of this industry should understand that the present invention is not limited by the above embodiments. The above embodiments and the description in the specification only illustrate the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. A gravity field inversion method taking into account the non-steady-state noise of satellite gravity observations, characterized in that The following steps are involved: A. Obtain the original observation data of the gravity satellite; B. Assume that the observed value noise is white noise and construct a random model that takes into account non-stationary noise; C. Perform parameter estimation to obtain the residuals after the observations are tested; D. Estimate the observation noise autocovariance and precision factor matrix based on the observed value post-validation residuals; E. Repeat steps C and D iteratively to refine the random model of gravity field parameters and observation values; F. Evaluate the time-varying gravity field model in the spectral domain.

2. The gravity field inversion method taking into account the non-steady-state noise of satellite gravity observation values ​​according to claim 1, characterized in that: The original observation data include K-band ranging observations, GPS orbit observations, prior accuracy information, AOD data and EOP data.

3. The gravity field inversion method taking into account the non-steady-state noise of satellite gravity observation values ​​according to claim 2, characterized in that: The random model taking into account non-stationary noise is Where C is the variance-covariance matrix of the steady-state noise, which characterizes the correlation between the observations, and Λ is the precision factor matrix, whose diagonal elements characterize the non-steady-state precision characteristics of the observations at each epoch.

4. The gravity field inversion method taking into account the non-steady-state noise of satellite gravity observation values ​​according to claim 3, characterized in that: In step C, the process of parameter estimation is as follows: in is the parameter estimate, A is the design matrix, the dimension of A is m×n, m is the number of observations, and n is the number of parameters; A T is the transpose of the design matrix, y is the observation vector, i.e. the actual observed data, (A T A) -1 A T is the pseudo-inverse matrix, used for least squares solution; The posterior residual refers to the difference between the predicted value obtained by estimating the parameters and the actual observed value. Its calculation formula is where v is the posterior residual vector, is the prediction vector, which is the expected observation value calculated based on the estimated parameters.

5. The gravity field inversion method taking into account the non-steady-state noise of satellite gravity observation values ​​according to claim 4, characterized in that: The estimation method of the observation noise autocovariance c is: where c k is the covariance between epochs j and j±k, η j is the post-test residual of the observation value at epoch j, N k To calculate c k Number of epoch pairs used, N a is the maximum number of delayed epochs; the observed noise autocovariance c is extended periodically, then the variance-covariance matrix C of the steady-state noise is a circulant matrix, and the product of the variance-covariance matrix C of the steady-state noise and the vector b is p=Cb, and the discrete Fourier transform b of the data is calculated s = DFT(b), calculate the power spectrum of the noise c s =DFT(c), calculate the frequency domain filtering result p s =c s b s , inverse transform restores time domain data p = IDFT (p s ).

6. The gravity field inversion method taking into account the non-steady-state noise of satellite gravity observation values ​​according to claim 5, characterized in that: The precision factor matrix is ​​estimated as follows: use Remove the correlation between residuals, where n is the residual vector of the original observation, with a dimension of N×1, N is the number of observations, and represents the difference between the observations and the model predictions, including systematic errors and random noise, n dc is the residual vector after decorrelation, with a dimension of N×1; Using Huber M estimator Calculate the precision factor corresponding to each epoch observation, where k is the adjustment constant of the Huber estimator; Set k = 2 based on empirical values ​​and calculate the posterior residual is the residual vector after decorrelation, with dimension N×1, where N is the number of observations.

7. The gravity field inversion method taking into account the non-steady-state noise of satellite gravity observation values ​​according to claim 6, characterized in that: In step F, the spectral domain accuracy evaluation method is used to quantitatively analyze the gravity field model, compare the solved gravity field model with the static gravity field spherical harmonic coefficients, and calculate the difference between the gravity field model coefficients; then calculate the order variance and the cumulative order variance to evaluate the accuracy index of the model.

Images