Satellite time-varying gravity field inversion method based on adaptive de-mixing

Abstract

The application discloses a satellite time-varying gravity field inversion method based on adaptive de-mixing, and belongs to the technical field of satellite gravity measurement. The method comprises the following steps: data processing; construction of a normal equation; accumulation of the normal equation and solution of initial values of parameters by using a weighted least square method; calculation of observation value posterior residual and recovery of pre-elimination parameters; application of zero mean constraint on daily time-varying gravity field parameters; establishment of an autoregressive model to describe the time correlation of the celestial solution sequence, application of the autoregressive model constraint, updating of the normal equation and joint solution of the time-varying gravity field model; output of the daily time-varying gravity field model and the monthly time-varying gravity field model according to the parameter estimation result; and evaluation of the time-varying gravity field model by using spectral domain analysis. The satellite time-varying gravity field inversion method based on adaptive de-mixing can adaptively remove the mixing error, improve the stability and precision of the time-varying gravity field inversion, and effectively weaken the influence of the background model mixing error in the satellite time-varying gravity field inversion.

Description

Technical Field

[0001] This invention relates to the field of satellite gravity measurement technology, and in particular to a satellite time-varying gravity field inversion method based on adaptive demixing. Background Technology

[0002] The Earth's time-varying gravity field reflects the mass distribution and movement within and on the Earth's surface, serving as a crucial constraint for studying the exchange of mass between land, atmosphere, and ocean, and its spatiotemporal variations. Since the beginning of the 21st century, with the successful implementation of low-Earth orbit satellite gravity exploration missions such as CHAMP, GRACE, and GOCE, satellite gravity measurement has entered a golden age of development. Numerous research results have been obtained from retrieving the Earth's time-varying gravity field based on satellite gravity observation data. This data has been widely applied in areas such as changes in surface and groundwater storage, polar glacier melting, and crustal deformation and mass transfer caused by earthquakes, greatly promoting the development of related disciplines in the geosciences. Its data has played a landmark role in the study of global hydrological cycles, glacial mass balance, and coseismic effects of earthquakes. Since the successful launch of GRACE and GRACE Follow-On satellites, numerous international research institutions have conducted extensive studies on inverting the Earth's time-varying gravity field using GRACE satellite data, and have solved a series of high-precision Earth time-varying gravity field models. The methods used mainly include the dynamic method, the short-arc method, the acceleration method, and the energy conservation method. However, regardless of the method used, the accuracy of the model solution is limited by the influence of mixing error. Therefore, how to overcome this problem and further improve the model accuracy has always been a research hotspot and challenge in the academic community.

[0003] Mixing errors, primarily caused by background model errors and signal undersampling, are one of the main error sources in GRACE time-varying gravity field inversion. GRACE is mainly used to invert seasonal and long-term variations such as terrestrial water storage and glacial melt; therefore, high-frequency signals such as ocean tides and atmospheric-oceanic non-tidal (AOD) terms need to be subtracted during data processing. Since these background models inevitably contain errors, their residual high-frequency signals are superimposed into the monthly time-varying gravity field, resulting in lower model accuracy. Furthermore, the undersampling of time-varying signals in GRACE observations also contributes to mixing errors.

[0004] Currently, methods for handling mixing errors can be broadly categorized into two types: one is post-processing filtering of the model, and the other is mitigation through improved gravity field parameterization methods. For the former, the main filtering methods currently employed include time-domain, spatial-domain, and spectral-domain filtering. Time-domain filtering estimates and removes the mixing period term from the time-varying gravity field sequence to reduce mixing errors. This method requires the satellite orbit to be a strictly repeating orbit and the mixing period to be clearly separated from the signal period, which is often difficult to meet in reality (for example, GRACE is not a strictly repeating orbit and its orbit gradually decreases over time). Spatial-domain filtering (such as Gaussian smoothing) and spectral-domain filtering (such as decorrelation filtering), while reducing mixing errors, inevitably cause signal loss and leakage. For the latter, since lower-order gravity field terms can be calculated using relatively little observational data, simulation studies have shown that selecting appropriate estimation orders and frequencies for these lower-order gravity field terms can effectively reduce mixing errors and improve model resolution (Wiese method). Other researchers have further demonstrated through simulations that this method, using two pairs of low-low tracking satellites (Bender configuration), can directly calculate non-tidal high-frequency atmospheric and oceanic signals without introducing prior models. However, all of the above studies are based on simulations, and the effectiveness of this method in processing GRACE measured data remains to be verified. Limited by mixing errors, the accuracy of existing gravity field models still falls short of design specifications by orders of magnitude, hindering their in-depth application in Earth sciences. Therefore, effectively reducing mixing errors is a key issue in further improving the accuracy of time-varying gravity field inversion. Summary of the Invention

[0005] The purpose of this invention is to provide a satellite time-varying gravity field inversion method based on adaptive demixing. By introducing an autoregressive (AR) model to describe the temporal correlation of the GRACE / GRACE-FO astronomical model sequence and applying stochastic process constraints, the monthly time-varying gravity field model is inverted, thereby reducing the dependence on traditional external prior model errors (such as AOD model, tide model, etc.), achieving adaptive demixing error, and then processing GRACE / GRACE-FO measured data and verifying its impact on time-varying gravity field inversion.

[0006] To achieve the above objectives, this invention provides a satellite time-varying gravity field inversion method based on adaptive demixing, comprising the following steps:

[0007] S1. Data processing is performed on satellite geometric orbit, simplified dynamic orbit, K-band ranging observations, accelerometer observations, and satellite attitude observations, including format conversion, data preprocessing, construction of motion equations, orbit integration, orbit fitting, and construction of observation equations.

[0008] S2. Based on the partial derivatives of the parameters to be estimated from the observations, the prior residuals OC from the observations, and the noise power spectral density from the observations, construct the normal equations.

[0009] S3. Accumulate the normal equations for each day of the month, and estimate the parameters using the weighted least squares method to obtain the initial solution of the parameters to be estimated, and calculate the parameter covariance matrix.

[0010] S4. Calculate the post-hoc residuals of the observed values ​​based on the observed value information, parameter estimates, and pre-elimination parameters, restore the pre-elimination parameters, and realize parameter updates;

[0011] S5. Apply zero-mean constraint to the diurnal time-varying gravity field parameters to restrict the mean of the diurnal time-varying gravity field parameters to zero, update the normal equations and jointly solve the monthly and daily solution parameters to obtain the monthly time-varying gravity field model and the diurnal time-varying gravity field model.

[0012] S6. Based on the estimated diurnal time-varying gravity field parameters, establish an autoregressive model to describe the time correlation of the celestial solution sequence, apply autoregressive model constraints to the diurnal time-varying gravity field parameters, update the normal equations, and jointly solve the monthly and celestial solution parameters to obtain the monthly time-varying gravity field model and the diurnal time-varying gravity field model.

[0013] S7. Based on the parameter estimation results, output the daily solution model and the monthly solution model, and perform subsequent analysis and output;

[0014] S8. The time-varying gravity field model was evaluated through spectral domain analysis and compared with traditional methods to demonstrate the improvement effect.

[0015] Preferably, the specific steps of S1 are as follows:

[0016] S11. Perform format conversion. Gravity satellite mission observations and orbital data are GRACE / GFO Level-1B data. Convert the GRACE / GFO Level-1B data, including KBR, ACC, SCA, and GNV data, and the external data, including AOD RL06 version data and Earth orientation parameters, into the software's internal format containing time and observations.

[0017] S12. Perform data preprocessing, calculate the geometric orbit and the prior residual OC of KBR observations based on the prior orbit, calculate the residual RMC, and identify the gross errors on an epoch-by-epoch basis.

[0018] S13. Calculate the satellite perturbation force using the prior mechanical model and satellite observation data, and construct the equation of motion;

[0019] S14. Using the RKF single-step numerical integration and ADAMS prediction correction multi-step integration methods, the satellite motion equations and variational equations are solved based on the satellite's initial state, model parameters, and three-dimensional perturbation acceleration.

[0020] S15. Update the satellite's initial state and model parameters by fusing the dynamic model with the geometric orbit observations;

[0021] S16. Read the satellite geometric orbit and KBR observations epoch by epoch, construct the observation equations of the geometric orbit and KBR, and calculate the partial derivatives of the observations and the prior residuals OC.

[0022] Preferably, the observation equations for the KBR in S16 include:

[0023] K-band inter-satellite ranging observations The observation equation:

[0024] (1);

[0025] in, This is the calculated value for GRACE inter-satellite ranging. and They represent time. The position vectors of satellite A and satellite B at that time. The distance between satellite A and satellite B. For GRACE interstellar light time correction, Correction for the interstellar phase center of GRACE This refers to the inter-satellite ranging bias;

[0026] K-band inter-satellite ranging rate observations The observation equation:

[0027] (2);

[0028] in, This is the calculated value for the GRACE inter-satellite ranging rate. and They represent time. The velocity vectors of satellite A and satellite B at that time, This represents the velocity vector of satellite B relative to satellite A. This represents the unit direction vector of satellite B relative to satellite A. For GRACE inter-satellite ranging rate optical time correction, Phase center correction for GRACE inter-satellite ranging rate.

[0029] Preferably, the specific steps of S3 are as follows:

[0030] S31. Define the core variables and matrix, and denote the vector of parameters to be estimated as... , , These are the observation vectors. The corresponding design matrix and covariance matrix, Let be the prior covariance matrix of the initial value vector of the parameters;

[0031] S32. Based on the prior information of the parameters, an error equation is constructed. The error equation for estimating the observed values ​​of the prior information of the parameters is as follows:

[0032] (3);

[0033] in, To observe the residual vector, Subtract the calculated value from the observed value. Correct the vector for the parameters. Let the initial value vector of the parameters be , Observation vector The weight matrix, for The weight matrix;

[0034] S33, using the weighted least squares principle Perform the solution, and obtain the parameter correction vector. for:

[0035] (4);

[0036] in, The residuals are the prior information of the parameters. The normal matrix of the observation equation, is the right vector of the observation equation;

[0037] S34, will Overlay This yields the updated vector of parameters to be estimated. Calculate the error covariance matrix of the updated parameter vector. ;

[0038] (5);

[0039] in, Correction vector for parameters The covariance matrix.

[0040] Preferably, the quadratic form of the post-test residual in S4 is:

[0041] (6);

[0042] in, This indicates the transpose operation.

[0043] Preferably, the specific steps of S5 are as follows:

[0044] S51. Assuming the average of the spherical harmonic coefficients for all days of each month is zero, define the following constraints:

[0045] (7);

[0046] in, For the number of days in a monthly time series, For the first The celestial harmonic coefficient value of the sun;

[0047] S52. Extract sub-blocks from the original normal equation matrix to construct the initial constraint matrix. ;

[0048] (8);

[0049] in, This is the original normal equation matrix. The number of spherical harmonic coefficients;

[0050] S53. Calculate constraint weights Constraint weights Calculated from the parameter variance estimate;

[0051] (9);

[0052] in, This represents the posterior variance of the spherical harmonic coefficients of the month. Amplification factor;

[0053] S54, Application of Weighting Factors Perform constraint matrix weighting;

[0054] (10);

[0055] in, This is the weighted constraint matrix. It is the identity matrix. The row and column index values;

[0056] S54. Multiply the transpose of the weighted constraint matrix by itself to expand the original normal equation matrix, thereby updating the normal equation matrix and ensuring that the parameter estimates satisfy the zero-mean constraint. The formula is as follows:

[0057] (11);

[0058] in, Used to determine the position of the current parameter in the normal equation matrix. This is the updated normal equation matrix.

[0059] Preferably, the specific steps of S6 are as follows:

[0060] S61. Establish the autoregressive model expression and define the time correlation of the Tianjie spherical harmonic coefficients;

[0061] (12);

[0062] in, These are the regression coefficients of the autoregressive model. It is a white noise sequence. The noise variance corresponds to the weight matrix. ; This represents the maximum order of the autoregressive model.

[0063] S62. Introduce autoregressive model constraints and update the normal equations;

[0064] S63. Solve for parameter estimates based on the updated normal equations.

[0065] Preferably, the specific steps of S62 are as follows:

[0066] S621. Transform the autoregressive model into a virtual observation equation. The expression for the virtual observation equation is:

[0067] (13);

[0068] Right now:

[0069] (14);

[0070] Among them, the Tianjie parameter is a time series. , The design matrix for the virtual observation equation is obtained by using each time point Design matrix block Stacked vertically, as shown below:

[0071] (15);

[0072] Each The form is:

[0073] (16);

[0074] Among them, the Column-oriented identity matrix , No. Column placement coefficient matrix All other positions are zero matrices. ;

[0075] S622. Merge the original observation equation with the virtual observation equation, and expand the design matrix and weight matrix;

[0076] The error of the original observation equation is:

[0077] (17);

[0078] in, Subtract the calculated value from the observed value. For joint design matrix, The design matrix for the monthly solution parameters, The design matrix for the celestial solution parameters, This is a joint parameter vector, containing both monthly and daily solution parameters;

[0079] The extended design matrix is ​​as follows:

[0080] (18);

[0081] The extended weight matrix is ​​as follows:

[0082] (19);

[0083] in, This is the weight matrix of the original observation equation. Let be the weight matrix of the virtual observation equations, representing the independent weighting of each virtual observation equation. For the Kronecker product, extend the weight matrix to all A virtual observation;

[0084] S623. Update the normal equations. The updated normal equations are:

[0085] (20);

[0086] Abbreviated as:

[0087] (twenty one);

[0088] Among them, in the matrix on the left side of the normal equation, and These are, respectively, the normal matrix corresponding to the monthly solution parameters, the conormal matrix between the monthly and celestial solution parameters, the normal matrix corresponding to the celestial solution parameters, and the AR constraint block. and These are the monthly solution observation vector block, the monthly solution observation vector block, and the AR constraint observation vector block, respectively.

[0089] Preferably, the specific steps of S63 are as follows:

[0090] S631. Solve the equations using the Cholesky decomposition method.

[0091] S632. By back-substituting and solving the parameter estimates, the time-varying gravity field model is obtained.

[0092] Therefore, the above-mentioned adaptive demixing-based satellite time-varying gravity field inversion method adopted in this invention has the following beneficial effects:

[0093] (1) Effectively reduce mixing error and improve the accuracy of time-varying gravity field inversion. By introducing autoregressive model constraints, the temporal correlation of the Tianjie gravity field model is directly modeled to achieve adaptive demixing error.

[0094] (2) Reduce dependence on external prior model errors. It avoids dependence on traditional external prior models (such as AOD model, tide model, etc.) and reduces the impact of background model errors on the results.

[0095] (3) Improve numerical stability and solution reliability. The introduction of zero-mean constraint and autoregressive model constraint ensures the stability and uniqueness of parameter estimation, avoids rank deficiency in the normal equation matrix, and prevents model overfitting. The positive definiteness of the constraint matrix is ​​verified by Cholesky decomposition, ensuring that the normal equation matrix is ​​invertible and avoiding parameter estimation divergence caused by ill-conditioned problems.

[0096] (4) The method described is applicable to the processing of gravity satellite data such as GRACE and GRACE-FO (Follow-On), which can effectively improve the accuracy of gravity field inversion and provide technical reserves for data processing of the next generation of gravity satellite missions.

[0097] (5) The high-precision time-varying gravity field model calculated by this method can support research on changes in terrestrial water storage and polar ice cap melting.

[0098] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0099] Figure 1 This is a flowchart of a satellite time-varying gravity field inversion method based on adaptive demixing according to the present invention;

[0100] Figure 2 This is a comparison of the higher-order errors of the geoid in the gravity field model calculated by different schemes of the satellite time-varying gravity field inversion method based on adaptive demixing according to the present invention. Detailed Implementation

[0101] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0102] Example

[0103] like Figure 1 As shown, this invention provides a satellite time-varying gravity field inversion method based on adaptive demixing, comprising the following steps:

[0104] S1. Data processing is performed on observational data such as satellite geometric orbit, simplified dynamic orbit, KBR observations, accelerometer observations, and satellite attitude observations, as well as multiple background models such as prior gravity field models, atmospheric and ocean dealiasing products, atmospheric tide and ocean tide models. This includes format conversion, data preprocessing, construction of motion equations, orbit integration, orbit fitting, and construction of observation equations.

[0105] S11. Perform format conversion. Gravity satellite mission observations and orbital data are GRACE / GFO Level-1B data. Convert the GRACE / GFO Level-1B data, including KBR, ACC, SCA, and GNV data, and external data, including AOD RL06 version data and Earth orientation parameters, into the software's internal format containing information such as time and observations.

[0106] S12. Perform data preprocessing, calculate the geometric orbit and the prior residual OC of KBR observations based on the prior orbit, calculate the residual RMC, and identify the gross errors on an epoch-by-epoch basis.

[0107] S13. Calculate the satellite perturbation force using the prior mechanical model and satellite observation data, and construct the equation of motion;

[0108] S14. Using the RKF single-step numerical integration and ADAMS prediction correction multi-step integration methods, the satellite motion equations and variational equations are solved based on the satellite's initial state, model parameters, and three-dimensional perturbation acceleration.

[0109] S15. Update the satellite's initial state and model parameters using the geometric orbit as observations. By fusing the dynamic model with the geometric orbit observations, update the satellite's initial state (position, velocity) and model parameters (gravity field coefficients, non-conservative force parameters, etc.).

[0110] S16. Based on the satellite geometric orbit and KBR observations, read the observation values ​​epoch by epoch, construct the observation equations of the geometric orbit and KBR, and calculate the partial derivatives of the observation values ​​and the prior residuals OC.

[0111] The observation equations of KBR include:

[0112] K-band inter-satellite ranging observations The observation equation:

[0113] (1);

[0114] in, This is the calculated value for GRACE inter-satellite ranging. and They represent time. The position vectors of satellite A and satellite B at that time. The distance between satellite A and satellite B. For GRACE interstellar light time correction, Correction for the interstellar phase center of GRACE This refers to the inter-satellite ranging bias;

[0115] K-band inter-satellite ranging rate observations The observation equation:

[0116] (2);

[0117] in, This is the calculated value for the GRACE inter-satellite ranging rate. and They represent time. The velocity vectors of satellite A and satellite B at that time, This represents the velocity vector of satellite B relative to satellite A. This represents the unit direction vector of satellite B relative to satellite A. For GRACE inter-satellite ranging rate optical time correction, Phase center correction for GRACE inter-satellite ranging rate.

[0118] S2. Based on the partial derivatives of the parameters to be estimated from the observations, the prior residuals OC of the observations, and the noise power spectral density of the observations, a normal equation is constructed for parameter estimation.

[0119] S3. Accumulate the normal equations for each day of the month, and estimate the parameters using the weighted least squares method to obtain the initial solution of the parameters to be estimated, and calculate the parameter covariance matrix.

[0120] S31. Define the core variables and matrix, and denote the vector of parameters to be estimated as... , , These are the observation vectors. The corresponding design matrix and covariance matrix, Let be the prior covariance matrix of the initial value vector of the parameters;

[0121] S32. Based on the prior information of the parameters, an error equation is constructed. The error equation for estimating the observed values ​​of the prior information of the parameters is as follows:

[0122] (3);

[0123] in, To observe the residual vector, Subtract the calculated value from the observed value. Correct the vector for the parameters. Let the initial value vector of the parameters be , Observation vector The weight matrix, for The weight matrix;

[0124] S33, using the weighted least squares principle Perform the solution, and obtain the parameter correction vector. for:

[0125] (4);

[0126] in, The residuals are the prior information of the parameters. This is the weight matrix. The normal matrix of the observation equation, Let be the right vector of the observation equation;

[0127] S34, will Overlay This yields the updated vector of parameters to be estimated. Based on the law of covariance propagation, calculate the error covariance matrix of the updated parameter vector. ;

[0128] (5);

[0129] in, Correction vector for parameters The covariance matrix.

[0130] S4. Calculate the post-hoc residuals of the observed values ​​based on the observed value information, parameter estimates, and pre-elimination parameters, restore the pre-elimination parameters, and realize parameter updates.

[0131] The quadratic form of the post-test residual is:

[0132] (6);

[0133] in, This indicates the transpose operation.

[0134] S5. Apply a zero-mean constraint to the diurnal time-varying gravity field parameters to restrict the mean of the parameters to zero, prevent overfitting of the model, and ensure the uniqueness of the parameter estimates. Update the normal equations and jointly solve the monthly and daily solutions for the parameters to obtain the monthly and diurnal time-varying gravity field models.

[0135] S51. The average value of the spherical harmonic coefficients for all days of the month is zero. The constraint condition is defined as follows:

[0136] (7);

[0137] in, For the number of days in a monthly time series, For the first The celestial harmonic coefficient value of the sun;

[0138] S52. Extract sub-blocks from the original normal equation matrix to construct the initial constraint matrix. ;

[0139] (8);

[0140] in, This is the original normal equation matrix. The number of spherical harmonic coefficients;

[0141] S53. Calculate constraint weights Constraint weights Calculated from the parameter variance estimate;

[0142] (9);

[0143] in, This represents the posterior variance of the spherical harmonic coefficients of the month. Amplification factor;

[0144] S54, Application of Weighting Factors Perform constraint matrix weighting;

[0145] (10);

[0146] in, This is the weighted constraint matrix. It is the identity matrix. Row and column index values; constraint matrix Dimensions It is determined by the maximum order of the model. and number of times Decide, The number of columns in the constraint matrix is ​​equal to the number of spherical harmonic coefficients. Monthly days; the constraint matrix is ​​weighted so that each row corresponds to a zero-mean constraint for a spherical harmonic coefficient, and the weighted sum of this coefficient for all days in a month is 0.

[0147] S54. Multiply the transpose of the weighted constraint matrix by itself and accumulate it into the original normal equation matrix to update the normal equation matrix, ensuring that the parameter estimates satisfy the zero-mean constraint. The formula is as follows:

[0148] (11);

[0149] in, Used to determine the position of the current parameter in the normal equation matrix. The updated normal equation matrix is ​​obtained by strengthening the diagonal terms of the normal equation matrix, which can prevent matrix singularity and improve numerical stability, while the cross terms suppress the inconsistency of parameters in different files.

[0150] S6. Based on the estimated diurnal time-varying gravity field parameters, establish an autoregressive model to describe the time correlation of the celestial solution sequence. Apply autoregressive model constraints to the diurnal time-varying gravity field parameters, update the normal equations, and jointly solve the monthly and celestial solution parameters to obtain the monthly time-varying gravity field model and the diurnal time-varying gravity field model.

[0151] Autoregressive (AR) models assume a linear relationship between current parameter values ​​and previous time-series values. This temporal correlation is particularly beneficial when processing time-series data, as it improves the accuracy and stability of parameter estimation. Using an AR model to describe the temporal correlation of the daily solution sequence can constrain diurnal time-varying gravity field parameters. This embodiment assumes that the monthly daily solution spherical harmonic coefficients conform to an autoregressive stochastic process.

[0152] S61. Establish the autoregressive model expression and define the time correlation of the Tianjie spherical harmonic coefficients;

[0153] (12);

[0154] in, These are the regression coefficients of the autoregressive model. It is a white noise sequence, reflecting the uncertainty of the parameters. The noise variance corresponds to the weight matrix. ; This represents the maximum order of the autoregressive model. and After smoothing the power spectral density (PSD) and converting it into an autocorrelation sequence, it is estimated based on the Levinson-Durbin recursive algorithm.

[0155] S62. Introduce autoregressive model constraints and update the normal equations.

[0156] Since the celestial solution parameters conform to an autoregressive stochastic process, it is necessary to construct a virtual observation equation to transform the AR model into a virtual observation equation. This equation serves as an additional constraint and is jointly adjusted with the original observation equation to solve the monthly time-varying gravity field model and the celestial time-varying gravity field model.

[0157] S621. Transform the autoregressive model into a virtual observation equation. The expression for the virtual observation equation is:

[0158] (13);

[0159] Right now:

[0160] (14);

[0161] Among them, the Tianjie parameter is a time series. , The design matrix for the virtual observation equation is obtained by using each time point Design matrix block Stacked vertically, as shown below:

[0162] (15);

[0163] Each The form is:

[0164] (16);

[0165] Among them, the Column-oriented identity matrix , No. Column placement coefficient matrix All other positions are zero matrices. .

[0166] S622. Merge the original observation equation with the virtual observation equation, and expand the design matrix and weight matrix;

[0167] The error of the original observation equation is:

[0168] (17);

[0169] in, Subtract the calculated value from the observed value. For joint design matrix, The design matrix for the monthly solution parameters, The design matrix for the celestial solution parameters, This is a joint parameter vector, containing both monthly and daily solution parameters;

[0170] The extended design matrix is ​​as follows:

[0171] (18);

[0172] The upper part relates the monthly and daily solution parameters, while the lower part only constrains the time correlation of the daily solution parameters.

[0173] The extended weight matrix is ​​as follows:

[0174] (19);

[0175] in, This is the weight matrix of the original observation equation. Let be the weight matrix of the virtual observation equations, representing the independent weighting of each virtual observation equation. For the Kronecker product, extend the weight matrix to all A virtual observation.

[0176] S623. To introduce an AR model to describe the time correlation of the astronomical solution sequence and constrain the diurnal time-varying gravity field parameters, the construction of the normal equations needs to combine the original observation equations and the virtual observation equations, and update the normal equations. The updated normal equations are as follows:

[0177] (20);

[0178] Abbreviated as:

[0179] (twenty one);

[0180] Among them, in the matrix on the left side of the normal equation, and These are, respectively, the normal matrix of the monthly solution parameter itself, the conormal matrix of the monthly solution parameter and the celestial solution parameter, the normal matrix of the celestial solution parameter itself, and the AR constraint block. and These are the monthly solution observation vector block, the monthly solution observation vector block, and the AR constraint observation vector block, respectively.

[0181] The sub-blocks are defined as follows:

[0182] The lunar section is as follows: , ;

[0183] The cross term is: ;

[0184] The Heavenly Solution section is as follows: , ;

[0185] AR constraint part: , Since the theoretical value of virtual observation is zero, the AR-constrained observation vector block... It is zero.

[0186] S63. Solve for parameter estimates based on the updated normal equations.

[0187] S631. Solve the equations using the Cholesky decomposition method.

[0188] S632. By back-substituting and solving the parameter estimates, the time-varying gravity field model is obtained.

[0189] S7. Based on the parameter estimation results, output the celestial solution model and the monthly time-varying gravity field model, and perform subsequent analysis and output.

[0190] S8. The time-varying gravity field model was evaluated through spectral domain analysis and compared with traditional methods to demonstrate the improvement effect.

[0191] The spectral domain accuracy assessment method is used to quantitatively analyze the gravity field model. The calculated gravity field model is compared with the spherical harmonic coefficients of static gravity fields (such as EIGEN-6C4 and GOCO06S), and the difference between the coefficients of the gravity field model is calculated. Then, the order error and the cumulative order error are calculated to evaluate the accuracy index of the model.

[0192] To analyze the effectiveness of this embodiment in solving time-varying gravity fields, this embodiment proposes to solve two sets of monthly time-varying gravity field models based on GRACE / GFO measured data: one set without constraints; the other set establishes an autoregressive model based on daily time-varying gravity field parameter estimates to describe the time correlation of the daily solution, and imposes constraints on it to adaptively eliminate mixing errors. By comparing the results of the two sets of models, the aim is to verify the effectiveness of this embodiment in time-varying gravity field inversion.

[0193] Based on L1B data from the GRACE satellite in May 2010, this embodiment studies and solves two sets of time-varying gravity field models and conducts comparative experiments: the baseline group (S1) is without constraints; the improved group (S2) establishes an autoregressive model to describe the time correlation of the solution, imposes constraints on it, and adaptively removes mixing errors. By comparing the results of the two sets of models, the aim is to verify the effectiveness of this embodiment in solving time-varying gravity fields. Figure 2 The comparison results of higher-order geoid errors for the two solutions are presented when GOCO06S is used as the reference background field. The red line represents the order error distribution of the S1 method, and the blue line corresponds to the improvement effect of the S2 method. The analysis results show that the cumulative 96th-order error of the Earth's gravity field model inverted by the S2 method is reduced by approximately 17.1%. This fully verifies the effect of the adaptive demixing method on improving the accuracy of the monthly time-varying gravity field inversion, and the S2 method can solve for a more accurate monthly time-varying gravity field model.

[0194] Therefore, the present invention adopts the above-mentioned satellite time-varying gravity field inversion method based on adaptive demixing, which can adaptively demixing error, improve the stability and accuracy of time-varying gravity field inversion, and effectively reduce the influence of background model mixing error in satellite time-varying gravity field inversion.

[0195] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A satellite time-varying gravity field inversion method based on adaptive de-mixing, characterized in that, Comprise the following steps: S1, satellite geometric orbit, simplified dynamics orbit, K-band ranging observation value, accelerometer observation value and satellite attitude observation value are processed, including format conversion, data preprocessing, motion equation construction, orbit integration, orbit fitting and observation equation construction; S2, based on the partial derivative of the observation value to the estimated parameter, the observation value prior residual O-C and the observation value noise power spectrum density, the normal equation is constructed; S3, the normal equation of each day in the month is accumulated, and the parameter estimation is carried out by weighted least squares method, the initial solution of the estimated parameter is obtained, and the parameter covariance matrix is calculated; S4, based on the observation value information, the parameter estimation and the pre-elimination parameter, the observation value posterior residual is calculated, the pre-elimination parameter is recovered, and the parameter updating is realized; S5, the daily time-varying gravity field parameter is subjected to zero mean constraint, the mean value of the daily time-varying gravity field parameter is limited to zero, the normal equation is updated and the monthly solution and daily solution parameters are solved to obtain the monthly time-varying gravity field model and the daily time-varying gravity field model; S6, according to the daily time-varying gravity field parameter estimation, an autoregressive model is established to describe the time correlation of the daily solution sequence, the autoregressive model constraint is applied to the daily time-varying gravity field parameter, the normal equation is updated and the monthly solution and daily solution parameters are solved to obtain the monthly time-varying gravity field model and the daily time-varying gravity field model; S7, according to the parameter estimation result, the daily solution model and the monthly solution model are output, and subsequent analysis and output are carried out; S8, the time-varying gravity field model is evaluated by spectral domain analysis, and compared with the traditional method, the improvement effect is proved.

2. The method according to claim 1, wherein, The specific steps of S1 are as follows: S11, format conversion is carried out, the gravity satellite task observation value and the orbit data are GRACE / GFO Level-1B data, the GRACE / GFO Level-1B data including KBR, ACC, SCA, GNV data and the external source data including AOD RL06 version data and earth orientation parameters are converted into software internal format containing time and observation value; S12, data preprocessing is carried out, the prior residual O-C of the geometric orbit and the KBR observation value is calculated based on the prior orbit, the residual RMC is calculated, and the gross error is marked every epoch; S13, satellite perturbation force is calculated from prior mechanical model and satellite observation data, and motion equation is constructed; S14, RKF single-step numerical integration and ADAMS prediction correction multi-step integration method are adopted, satellite motion equation and variation equation are solved according to satellite initial state, model parameters and three-dimensional perturbation acceleration; S15, satellite initial state and model parameters are updated by fusing the dynamic model and the geometric orbit observation value; S16, satellite geometric orbit and KBR observation value are read every epoch, observation equation of geometric orbit and KBR is constructed, observation value partial derivative and prior residual O-C are calculated.

3. The method according to claim 2, wherein, The observation equation of KBR in S16 comprises: K-band inter-satellite ranging observations Observation equation: （1） wherein is the calculated value for the GRACE inter-satellite range, and denote the position vectors of satellite A and satellite B at times is the distance between satellite A and satellite B, is the GRACE inter-satellite light time correction, is the GRACE inter-satellite phase center correction, is the inter-satellite range bias;​ K-band inter-satellite ranging rate observations Observation equation: （2） where is the GRACE inter-satellite range rate computed value, and denote the velocity vectors of satellite A and satellite B at time denotes the velocity vector of satellite B relative to satellite A, denotes the unit direction vector of satellite B relative to satellite A, is the GRACE inter-satellite range rate optical time correction, is the GRACE inter-satellite range rate phase center correction.​ 4. The method according to claim 1, wherein, The specific steps of S3 are as follows: S31, define core variables and matrices, denote the vector of parameters to be estimated as , , are the observation vectors corresponding design matrix and covariance matrix, is the prior covariance matrix of the initial value vector of parameters; S32, error equation is constructed based on parameter prior information, the error equation of the observation equation of parameter prior information is: （3） in, To observe the residual vector, Subtract the calculated value from the observed value. Correct the vector for the parameters. Let the initial value vector of the parameters be , Observation vector The weight matrix, for The weight matrix; S33、through the principle of weighted least squares Solve, the parameter correction vector obtained by solving is: （4） wherein, is a residual of the parameter prior information, is a normal matrix of the observation equation, is a right vector of the observation equation; S34, will Overlay This yields the updated vector of parameters to be estimated. Calculate the error covariance matrix of the updated parameter vector. ; （5） wherein is a parameter correction vector covariance matrix of the parameter correction vector 5. The method according to claim 4, wherein, The quadratic form of the posterior residual in S4 is: （6） wherein denotes a transpose operation.

6. The method according to claim 1, wherein, The specific steps of S5 are as follows: S51, it is assumed that the average value of all daily solution spherical harmonic coefficients in each month is zero, and the constraint condition is defined as: （7） wherein, is the day of the month, is the day of the month, is the diurnal zonal harmonic coefficient value for the day of the month; S52, extract sub-blocks of the original normal equation matrix to construct an initial constraint matrix ; （8） wherein, is the original normal equation matrix, is the number of spherical harmonic coefficients; S53, compute constraint weight , constraint weight computed from parameter variance estimates; （9） wherein, denotes the posterior variance of the zonal harmonic coefficients, is an amplification factor; S54, applying weight factors performing constraint matrix weighting; （10） wherein, is the weighted constraint matrix, is the identity matrix, is the row and column index value; S54, multiply the transpose of the weighted constraint matrix with itself to expand the original normal equation matrix, and update the normal equation matrix to make the parameter estimation meet the zero-mean constraint, the formula is as follows: （11） wherein for determining the position of the current parameter in the normal equation matrix, is the updated normal equation matrix.

7. The method according to claim 1, wherein, The specific steps of S6 are as follows: S61, establish an autoregressive model expression to define the time correlation of the spherical harmonic coefficients of the satellite orbit; （12） wherein is a regression coefficient of an autoregressive model, is a white noise sequence, is a noise variance, corresponding to the weight matrix ; is a maximum order of the autoregressive model; S62, introduce the autoregressive model constraint to update the normal equation; S63, solve the parameter estimation based on the updated normal equation.

8. The method according to claim 7, wherein, The specific steps of S62 are as follows: S621, convert the autoregressive model into a virtual observation equation, and the expression of the virtual observation equation is as follows: （13） That is: （14） where the state parameters are time series , is the design matrix of the virtual observation equation, which is formed by vertically stacking the design matrix blocks at each time point and is expressed as follows: （15） each in the form of: （16） wherein the first column places the identity matrix , the first column places the coefficient matrix , , and the other positions are zero matrices ; S622, combine the original observation equation with the virtual observation equation to expand the design matrix and the weight matrix; The error of the original observation equation is as follows: （17） wherein is the computed value, is the combined design matrix, is the design matrix for the orbital parameters, is the design matrix for the ephemeris parameters, is the combined parameter vector, comprising the orbital parameters and the ephemeris parameters; The expanded design matrix is as follows: （18） The expanded weight matrix is as follows: （19） wherein, is a weight matrix of the original observation equation, is a weight matrix of the virtual observation equation, indicating that each virtual observation equation is independently weighted, is a Kronecker product, which extends the weight matrix to all virtual observations; S623, update the normal equation, and the updated normal equation is as follows: （20） S63, the specific steps are as follows: （21） wherein the matrix on the left side of the normal equation, and are the normal matrix corresponding to the monthly solution parameters, the copula matrix of the monthly solution parameters and the daily solution parameters, the normal matrix corresponding to the daily solution parameters, and the AR constraint block, and are the monthly solution observation vector block, the monthly solution observation vector block, and the AR constraint observation vector block, respectively.

9. The method according to claim 7, wherein, S631, use Cholesky decomposition to solve the normal equation; S632, solve the parameter estimation by back substitution to obtain the time-varying gravity field model. ​

Images