CN120831725B - A satellite time-varying gravity field inversion method based on adaptive demixing - Google Patents

Abstract

The invention discloses a satellite time-varying gravitational field inversion method based on self-adaptive demixing, which belongs to the technical field of satellite gravity measurement and comprises the steps of carrying out data processing, constructing a normal equation, accumulating the normal equation and solving a parameter initial value through a weighted least square method, calculating an observation value post-test residual error and recovering a pre-elimination parameter, applying zero mean constraint to a time-varying gravitational field parameter, establishing an autoregressive model to describe time correlation of a time solution sequence, applying an autoregressive model constraint, updating the normal equation and solving a time-varying gravitational field model in a combined mode, outputting the time-varying gravitational field model and a time-varying gravitational field model according to a parameter estimation result, and evaluating the time-varying gravitational field model through spectral domain analysis.

Description

Satellite time-varying gravitational field inversion method based on self-adaptive mixing

Technical Field

The invention relates to the technical field of satellite gravity measurement, in particular to a satellite time-varying gravity field inversion method based on self-adaptive mixing.

Background

The earth time-varying gravitational field reflects the mass distribution and motion inside and on the surface of the earth and is an important constraint for researching the mutual exchange of land, atmosphere and ocean mass and the space-time variation thereof. In the 21 st century, with successful implementation of the low-orbit satellite gravity detection task CHAMP, GRACE, GOCE, satellite gravity measurement is in the coming gold development period, a great amount of research results are obtained by inverting the earth time-varying gravity field based on satellite gravity observation data, and the satellite gravity detection satellite has wide application in the fields of surface and underground water reserve change, polar glacier ablation, earth crust deformation caused by earthquakes, mass migration and the like, greatly promotes the development of relevant disciplines in the field of geochemistry, and has data playing a milestone role in global hydrologic cycle, glacier material balance and seismic isoseism effect research. Since GRACE and GRACE Follow-On satellites are successfully launched, a large number of related researches On inversion of the earth time-varying gravitational field by utilizing GRACE satellite data are developed by international scientific research institutions, a series of high-precision earth time-varying gravitational field models are solved, the adopted methods mainly comprise a dynamics method, a short arc method, an acceleration method and an energy conservation method, however, the model resolving precision is limited by the influence of mixing errors no matter what method is adopted, so that how to further improve the model precision is always a research hot spot and a difficulty in the academic world.

The mixing error is mainly caused by background model error and signal undersampling, and is one of the main error sources of GRACE time-varying gravitational field inversion. Grace is mainly used for inverting seasonal and long-term changes such as land water reserves, glacier ablation and the like, so that high-frequency signals such as ocean tides, atmospheric ocean non-tides (AODs) and the like need to be subtracted in data processing. Since these background models inevitably contain errors, their residual high frequency signals are aliased into the time-varying gravitational field of the month, resulting in lower model resolution accuracy. Furthermore, undersampling of the time-varying signal by the GRACE observations can also lead to mixing errors.

The current methods for processing the mixing errors are roughly divided into two types, namely model post-filtering processing and attenuation by improving a gravity field parameterization method. As for the former, currently mainly adopted filtering methods include time domain filtering, spatial domain filtering and spectral domain filtering, wherein the time domain filtering estimates and removes a mixing period term from a time-varying gravitational field sequence so as to weaken mixing errors, the method requires that satellite orbits are strictly repeated orbits and the mixing period can be obviously separated from signal periods, which is generally difficult to meet in reality (such as GRACE is not strictly repeated orbits and gradually reduces along with time), and spatial domain filtering (such as gaussian smoothing) and spectral domain filtering (such as decorrelation filtering) inevitably cause signal loss and leakage while weakening the mixing errors. For the latter, as the gravity field low-order term can be resolved by using less observation data, simulation researches show that the mixing error can be effectively weakened and the model resolution can be improved by selecting proper gravity field low-order term estimation order and estimation frequency (Wiese method), and other students further show that the solution of the atmospheric ocean non-tidal high-frequency signal can be directly realized by using two pairs of low tracking satellites (Bender configuration) based on the method through the simulation researches without introducing a priori model. However, the above researches are based on simulation expansion, and the effectiveness of the method in GRACE measured data processing is yet to be verified. The accuracy of the existing gravitational field model is still in an order-of-magnitude difference from the design index due to the mixing error, and the deep application of the model in the earth science is restricted. Therefore, how to effectively weaken the mixing error is a key problem for further improving the inversion accuracy of the time-varying gravity field.

Disclosure of Invention

The invention aims to provide a satellite time-varying gravitational field inversion method based on self-adaptive mixing removal, which is characterized in that an Autoregressive (AR) model is introduced to describe the time correlation of GRACE/GRACE-FO sky solution model sequences, random process constraint is applied, a month-varying gravitational field model is inverted, dependence on traditional external priori model errors (such as an AOD model, a sea tide model and the like) is weakened, self-adaptive mixing removal errors are realized, and then the GRACE/GRACE-FO actual measurement data is processed and the influence on the time-varying gravitational field inversion is verified.

In order to achieve the above purpose, the invention provides a satellite time-varying gravitational field inversion method based on self-adaptive demixing, which comprises the following steps:

S1, performing data processing on a satellite geometric orbit, a simplified dynamic orbit, a K-band ranging observation value, an accelerometer observation value and a satellite attitude observation value, wherein the data processing comprises format conversion, data preprocessing, motion equation construction, orbit integration, orbit fitting and observation equation construction;

S2, constructing a normal equation based on the partial derivative of the parameter to be estimated, the priori residual error O-C of the observed value and the noise power spectral density of the observed value;

s3, accumulating the daily normal equations of the whole month, carrying out parameter estimation through a weighted least square method to obtain an initial solution of parameters to be estimated, and calculating a parameter covariance matrix;

s4, calculating an observation value post-test residual error based on the observation value information, the parameter estimation value and the pre-elimination parameter, recovering the pre-elimination parameter, and realizing parameter updating;

S5, applying zero mean constraint to the time-of-day gravitational field parameters, limiting the mean value of the time-of-day gravitational field parameters to be zero, updating a law equation, and jointly solving a month solution and a day solution parameter to obtain a month-of-day gravitational field model and a time-of-day gravitational field model;

S6, establishing an autoregressive model according to the estimated value of the time-of-day gravitational field parameters, describing the time correlation of a time-of-day gravitational field sequence, applying autoregressive model constraint to the time-of-day gravitational field parameters, updating a method equation, and jointly solving a month solution and the time-of-day gravitational field parameters to obtain a month-varying gravitational field model and a time-of-day gravitational field model;

s7, outputting a model of a natural solution and a model of a lunar solution according to the parameter estimation result, and carrying out subsequent analysis and output;

S8, evaluating the time-varying gravity field model through spectral domain analysis, and comparing the time-varying gravity field model with the traditional method to prove the improvement effect.

Preferably, the specific steps of S1 are as follows:

S11, carrying out format conversion, wherein the observed value and orbit data of a gravity satellite task are GRACE/GFO Level-1B data, and converting GRACE/GFO Level-1B data comprising KBR, ACC, SCA, GNV data and exogenous data comprising AOD RL06 version data and earth orientation parameters into software internal formats comprising time and observed values;

S12, data preprocessing is carried out, a priori residual error O-C of the geometric orbit and the KBR observation value is calculated based on the priori orbit, a residual error RMC is calculated, and the rough mark difference is judged from epoch to epoch;

s13, calculating satellite uptake power by using the prior mechanical model and satellite observation data, and constructing a motion equation;

S14, solving a satellite motion equation and a variation equation according to the initial state, model parameters and three-dimensional perturbation acceleration of the satellite by adopting an RKF single-step numerical integration and ADAMS forecast correction multi-step integration method;

S15, updating the initial state and model parameters of the satellite by fusing the dynamic model and the geometrical orbit observation values;

S16, reading satellite geometric orbits and KBR observation values from epoch to epoch, constructing an observation equation of the geometric orbits and KBR, and calculating partial derivatives and priori residual errors O-C of the observation values.

Preferably, the observation equation of KBR in S16 includes:

k-band inter-satellite ranging observation value Is defined by the observation equation:

(1);

Wherein, the For the GRACE inter-satellite ranging calculation,AndRespectively expressed in timeThe position vectors of satellite a and satellite B,For the distance between satellite a and satellite B,For the GRACE inter-satellite light time correction,For the correction of the GRACE inter-satellite phase center,The inter-satellite distance measurement deviation;

K-band inter-satellite ranging rate observation value Is defined by the observation equation:

(2);

Wherein, the For the GRACE inter-satellite ranging rate calculation,AndRespectively expressed in timeThe velocity vectors of satellite a and satellite B,Representing the velocity vector of satellite B relative to satellite a,Representing the unit direction vector of satellite B relative to satellite a,For GRACE inter-satellite ranging rate light time correction,Phase center correction for GRACE inter-satellite ranging rate.

Preferably, the specific steps of S3 are as follows:

S31, defining core variables and matrixes, and recording parameter vectors to be estimated as ,,Respectively, are observed value vectorsThe corresponding design matrix and covariance matrix,A priori covariance matrix which is a parameter initial value vector;

S32, constructing an error equation based on the parameter priori information, wherein the error equation of the observation value equation for estimating the parameter priori information is as follows:

(3);

Wherein, the In order to observe the residual vector of the image,In order to subtract the calculated value from the observed value,For the parameter correction vector(s),As the initial value vector of the parameter,For vectors of observationsIs used for the weight matrix of the (c),Is thatIs a weight matrix of (2);

S33, through the weighted least square principle Performing a solution to obtain a parameter correction vectorThe method comprises the following steps:

(4);

Wherein, the As a residual of the parameter a priori information,In order to observe the normal matrix of the equation,Is the right vector of the observation equation;

S34, will Superimposed onObtaining updated parameter vector to be estimatedCalculating an error covariance matrix of the updated parameter vector;

(5);

Wherein, the Correction vector for parameterIs a covariance matrix of (2).

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

(6);

Wherein, the Representing the transpose operation.

Preferably, the specific steps of S5 are as follows:

s51, assuming that the average value of the spherical aberration coefficients of all days of each month is zero, defining constraint conditions as follows:

(7);

Wherein, the For the number of days in the time series of months,Is the firstDay-to-day spherical harmonic coefficient values;

S52, extracting sub-blocks of the original normal equation matrix to construct an initial constraint matrix ;

(8);

Wherein, the As an original matrix of normal equations,The number of the spherical harmonic coefficients;

s53, calculating constraint weight Constraint weightsCalculating by the parameter variance estimation value;

(9);

Wherein, the The post-test variance of the lunar spin coefficient is represented,Is an amplification factor;

S54, applying a weight factor Weighting the constraint matrix;

(10);

Wherein, the For the weighted constraint matrix,Is a matrix of units which is a matrix of units,Is a rank index value;

S54, multiplying the transpose of the weighted constraint matrix by the transpose, expanding an original normal equation matrix, and updating the normal equation matrix to enable the parameter estimation to meet zero mean constraint, wherein the formula is as follows:

(11);

Wherein, the For determining the position of the current parameter in the normal equation matrix,Is an updated normal equation matrix.

Preferably, the specific steps of S6 are as follows:

s61, establishing an autoregressive model expression, and defining the time correlation of the celestial sphere harmonic coefficients;

(12);

Wherein, the As the regression coefficients of the autoregressive model,In the form of a white noise sequence,For noise variance, corresponding weight matrix;Is the maximum order of the autoregressive model;

S62, introducing constraint of an autoregressive model, and updating a normal equation;

And S63, solving the parameter estimation value based on the updated normal equation.

Preferably, the specific steps of S62 are as follows:

s621, converting the autoregressive model into a virtual observation equation, wherein the expression of the virtual observation equation is as follows:

(13);

Namely:

(14);

wherein the sky solution parameter is time sequence ,For designing a matrix of virtual observation equations by applying each time pointIs a design matrix block of (a)Longitudinally stacked, and expressed as follows:

(15);

Each of which is The form of (2) is:

(16);

Wherein, the first Column placement unit matrixFirst, theColumn placement coefficient matrixOther positions are zero matrix;

S622, combining the original observation equation and the virtual observation equation, and expanding a design matrix and a weight matrix;

The error of the original observation equation is:

(17);

Wherein, the In order to subtract the calculated value from the observed value,For the joint design matrix to be a matrix,The matrix is designed for the monthly solution parameters,The matrix is designed for the parameters of the solution,The method comprises the steps of (1) as a joint parameter vector, including a month solution parameter and a day solution parameter;

the expanded design matrix is:

(18);

the extended weight matrix is:

(19);

Wherein, the Is a weight matrix of the original observation equation,For a weight matrix of virtual observation equations, each virtual observation equation is represented with an independent weight,For Kronecker product, the weight matrix is extended to allVirtual observations;

s623, updating a normal equation, wherein the updated normal equation is as follows:

(20);

The abbreviation is:

(21);

wherein, in the matrix on the left side of the normal equation, And (3) withThe method comprises a method matrix corresponding to the lunar solution parameter, a collaborative method matrix of the lunar solution parameter and the Tianjie parameter, a method matrix corresponding to the Tianjie parameter and an AR constraint block respectively,And (3) withThe lunar observation vector block, the lunar observation vector block and the AR constraint observation vector block are respectively adopted.

Preferably, the specific steps of S63 are as follows:

S631, solving a solution equation by using Cholesky decomposition;

s632, solving the parameter estimation value through the back generation to obtain a time-varying gravity field model.

Therefore, the satellite time-varying gravitational field inversion method based on self-adaptive demixing has the following beneficial effects:

(1) The mixing error is effectively weakened, and the inversion precision of the time-varying gravity field is improved. By introducing the constraint of the autoregressive model, the time sequence correlation of the gravitational field model of the natural solution is directly modeled, and the self-adaptive mixing error removal is realized.

(2) The dependence on the external prior model error is reduced. The dependence on traditional external prior models (such as an AOD model, a sea tide model and the like) is avoided, and the influence of background model errors on results is reduced.

(3) The numerical stability and the reliability of the solution are improved. Zero mean constraint and autoregressive model constraint are introduced, stability and uniqueness of parameter estimation are ensured, rank and deficiency of a French equation matrix are avoided, and model overfitting is prevented. Positive qualification of the constraint matrix is verified through Cholesky decomposition, reversibility of the normal equation matrix is ensured, and parameter estimation divergence caused by pathological problems is avoided.

(4) The method is suitable for processing gravity satellite data such as GRACE and GRACE-FO (Follow-On), can effectively improve the accuracy of gravity field inversion, and can provide technical reserve for data processing of next-generation gravity satellite tasks.

(5) The high-precision time-varying gravity field model calculated based on the method can support research of land water reserve change, polar ice cover ablation and the like.

The technical scheme of the invention is further described in detail through the drawings and the embodiments.

Drawings

FIG. 1 is a flow chart of a satellite time-varying gravitational field inversion method based on adaptive demixing;

FIG. 2 is a graph of high-order errors of the ground level of a gravity field model calculated by different schemes of the satellite time-varying gravity field inversion method based on self-adaptive demixing.

Detailed Description

The following detailed description of the embodiments of the invention, provided in the accompanying drawings, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.

Examples

As shown in fig. 1, the invention provides a satellite time-varying gravitational field inversion method based on self-adaptive demixing, which comprises the following steps:

S1, carrying out data processing on observation data such as satellite geometric orbits, simplified dynamic orbits, KBR observation values, accelerometer observation values and satellite attitude observation values and a plurality of background models such as priori gravity field models, atmosphere and ocean anti-aliasing products, atmosphere tide and sea tide models and the like, wherein the data processing comprises format conversion, data preprocessing, motion equation construction, orbit integration, orbit fitting and observation equation construction.

S11, carrying out format conversion, wherein the observed value and orbit data of a gravity satellite task are GRACE/GFO Level-1B data, and converting GRACE/GFO Level-1B data comprising KBR, ACC, SCA, GNV data and exogenous data comprising AOD RL06 version data, earth orientation parameters and the like into software internal formats containing time, observed value and the like;

S12, data preprocessing is carried out, a priori residual error O-C of the geometric orbit and the KBR observation value is calculated based on the priori orbit, a residual error RMC is calculated, and the rough mark difference is judged from epoch to epoch;

s13, calculating satellite uptake power by using the prior mechanical model and satellite observation data, and constructing a motion equation;

S14, solving a satellite motion equation and a variation equation according to the initial state, model parameters and three-dimensional perturbation acceleration of the satellite by adopting an RKF single-step numerical integration and ADAMS forecast correction multi-step integration method;

s15, updating the initial state and model parameters of the satellite by taking the geometric orbit as an observation value, and updating the initial state (position and speed) and model parameters (gravitational field coefficient, non-conservative force parameter and the like) of the satellite by fusing a dynamics model and the geometric orbit observation value;

S16, reading the observed value epoch by epoch according to the satellite geometric orbit and the KBR observed value, constructing an observation equation of the geometric orbit and the KBR, and calculating the partial derivative of the observed value and the prior residual error O-C.

The observation equation for KBR includes:

k-band inter-satellite ranging observation value Is defined by the observation equation:

(1);

Wherein, the For the GRACE inter-satellite ranging calculation,AndRespectively expressed in timeThe position vectors of satellite a and satellite B,For the distance between satellite a and satellite B,For the GRACE inter-satellite light time correction,For the correction of the GRACE inter-satellite phase center,The inter-satellite distance measurement deviation;

K-band inter-satellite ranging rate observation value Is defined by the observation equation:

(2);

Wherein, the For the GRACE inter-satellite ranging rate calculation,AndRespectively expressed in timeThe velocity vectors of satellite a and satellite B,Representing the velocity vector of satellite B relative to satellite a,Representing the unit direction vector of satellite B relative to satellite a,For GRACE inter-satellite ranging rate light time correction,Phase center correction for GRACE inter-satellite ranging rate.

S2, constructing a normal equation based on the partial derivative of the parameter to be estimated, the priori residual error O-C of the observed value and the noise power spectral density of the observed value, and using the normal equation for parameter estimation.

And S3, accumulating the daily normal equations of the whole month, carrying out parameter estimation by a weighted least square method to obtain an initial solution of the parameters to be estimated, and calculating a parameter covariance matrix.

S31, defining core variables and matrixes, and recording parameter vectors to be estimated as,,Respectively, are observed value vectorsThe corresponding design matrix and covariance matrix,A priori covariance matrix which is a parameter initial value vector;

S32, constructing an error equation based on the parameter priori information, wherein the error equation of the observation value equation for estimating the parameter priori information is as follows:

(3);

Wherein, the In order to observe the residual vector of the image,In order to subtract the calculated value from the observed value,For the parameter correction vector(s),As the initial value vector of the parameter,For vectors of observationsIs used for the weight matrix of the (c),Is thatIs a weight matrix of (2);

S33, through the weighted least square principle Performing a solution to obtain a parameter correction vectorThe method comprises the following steps:

(4);

Wherein, the As a residual of the parameter a priori information,As a matrix of weights, the weight matrix,In order to observe the normal matrix of the equation,Is the right vector of the observation equation;

S34, will Superimposed onObtaining updated parameter vector to be estimatedCalculating an error covariance matrix of the updated parameter vector according to a covariance propagation rule;

(5);

Wherein, the Correction vector for parameterIs a covariance matrix of (2).

And S4, calculating an observation value post-test residual error based on the observation value information, the parameter estimation value and the pre-elimination parameter, recovering the pre-elimination parameter, and realizing parameter updating.

The quadratic form of the post-test residual is:

(6);

Wherein, the Representing the transpose operation.

S5, zero-mean constraint is applied to the time-of-day gravitational field parameters, the mean value of the time-of-day gravitational field parameters is limited to be zero, the model is prevented from being fitted excessively, and the uniqueness of parameter estimation is ensured. Updating a law equation and jointly solving a month solution and a day solution parameter to obtain a month-varying gravitational field model and a day-varying gravitational field model.

S51, the average value of the spherical aberration coefficient of all the days of each month is zero, and constraint conditions are defined as follows:

(7);

Wherein, the For the number of days in the time series of months,Is the firstDay-to-day spherical harmonic coefficient values;

S52, extracting sub-blocks of the original normal equation matrix to construct an initial constraint matrix ;

(8);

Wherein, the As an original matrix of normal equations,The number of the spherical harmonic coefficients;

s53, calculating constraint weight Constraint weightsCalculating by the parameter variance estimation value;

(9);

Wherein, the The post-test variance of the lunar spin coefficient is represented,Is an amplification factor;

S54, applying a weight factor Weighting the constraint matrix;

(10);

Wherein, the For the weighted constraint matrix,Is a matrix of units which is a matrix of units,Is a row-column index value, constraint matrixThe dimension isWhich is defined by the maximum order of the modelSum timesIt is decided that the method comprises the steps of,To constrain the number of columns of the matrix equal to the number of spherical harmonic coefficientsThe number of days of month, constraint matrix weighting is carried out so that each row corresponds to zero mean constraint condition of one spherical harmonic coefficient, and the weighted sum of the coefficients in all days of month is 0.

S54, multiplying the transpose of the constraint matrix after weighting by the transpose, accumulating the transpose into an original normal equation matrix, and updating the normal equation matrix to enable the parameter estimation to meet zero mean constraint, wherein the formula is as follows:

(11);

Wherein, the For determining the position of the current parameter in the normal equation matrix,In order to update the normal equation matrix, the diagonal terms of the normal equation matrix are enhanced, so that the singular of the matrix can be prevented, the numerical stability is improved, and the cross terms inhibit the inconsistency of parameters in different files.

And S6, establishing an autoregressive model according to the time-of-day gravitational field parameter estimation, describing the time correlation of a time-of-day gravitational field sequence, applying autoregressive model constraint to the time-of-day gravitational field parameter, updating a method equation, and jointly solving a month solution and the time-of-day gravitational field parameter to obtain a month-varying gravitational field model and a time-of-day gravitational field model.

The Autoregressive (AR) model assumes that the current parameter value has a linear relationship with the parameter value of the previous time node, and in particular, by introducing such a temporal correlation when processing time series data, the accuracy and stability of the parameter estimation can be improved. Using an autoregressive model to describe the time dependence of the sky solution sequence, the time-of-day gravitational field parameters can be constrained, and the present embodiment assumes that the monthly sky solution spherical harmonic coefficients conform to an autoregressive stochastic process.

S61, establishing an autoregressive model expression, and defining the time correlation of the celestial sphere harmonic coefficients;

(12);

Wherein, the As the regression coefficients of the autoregressive model,For white noise sequences, reflecting the uncertainty of the parameters,For noise variance, corresponding weight matrix;Is the maximum order of the autoregressive model.And (3) withSmoothing by Power Spectral Density (PSD) and converting into autocorrelation sequence, estimating based on Levinson-Durbin recursive algorithm.

S62, introducing constraint of the autoregressive model, and updating a normal equation.

Because the sky solution parameters accord with the autoregressive stochastic process, a virtual observation equation needs to be constructed, the AR model is converted into the virtual observation equation, the equation is used as an additional constraint, and is combined with the original observation equation to adjust, so that the month-varying gravity field model and the sky-varying gravity field model are solved.

S621, converting the autoregressive model into a virtual observation equation, wherein the expression of the virtual observation equation is as follows:

(13);

Namely:

(14);

wherein the sky solution parameter is time sequence ,For designing a matrix of virtual observation equations by applying each time pointIs a design matrix block of (a)Longitudinally stacked, and expressed as follows:

(15);

Each of which is The form of (2) is:

(16);

Wherein, the first Column placement unit matrixFirst, theColumn placement coefficient matrixOther positions are zero matrix。

S622, combining the original observation equation and the virtual observation equation, and expanding a design matrix and a weight matrix;

The error of the original observation equation is:

(17);

Wherein, the In order to subtract the calculated value from the observed value,For the joint design matrix to be a matrix,The matrix is designed for the monthly solution parameters,The matrix is designed for the parameters of the solution,The method comprises the steps of (1) as a joint parameter vector, including a month solution parameter and a day solution parameter;

the expanded design matrix is:

(18);

The upper half correlates the monthly solution parameters with the daily solution parameters, and the lower half only constrains the time correlation of the daily solution parameters.

The extended weight matrix is:

(19);

Wherein, the Is a weight matrix of the original observation equation,For a weight matrix of virtual observation equations, each virtual observation equation is represented with an independent weight,For Kronecker product, the weight matrix is extended to allAnd (5) virtual observation.

S623, in order to introduce the time correlation of the AR model description sky solution sequence, constraining the time-varying gravity field parameters, constructing a normal equation by combining an original observation equation and a virtual observation equation, updating the normal equation, wherein the updated normal equation is as follows:

(20);

The abbreviation is:

(21);

wherein, in the matrix on the left side of the normal equation, And (3) withThe method matrix of the lunar solution parameter, the collaborative matrix of the lunar solution parameter and the Tianjie parameter, the method matrix of the Tianjie parameter and the AR constraint block are respectively adopted,And (3) withThe lunar observation vector block, the lunar observation vector block and the AR constraint observation vector block are respectively adopted.

Each sub-block is defined as follows:

The month solution part is as follows: ,;

The cross terms are: ;

The heaven solution part is as follows: ,;

AR constraint part: , because the theoretical value of the virtual observation is zero, the AR constraint observation vector block Zero.

And S63, solving the parameter estimation value based on the updated normal equation.

S631, solving a solution equation by using Cholesky decomposition;

s632, solving the parameter estimation value through the back generation to obtain a time-varying gravity field model.

And S7, outputting a heaven solution model and a moon time-varying gravity field model according to the parameter estimation result, and carrying out subsequent analysis and output.

S8, evaluating the time-varying gravity field model through spectral domain analysis, and comparing the time-varying gravity field model with the traditional method to prove the improvement effect.

The gravity field model is quantitatively analyzed by adopting a spectral domain precision evaluation method, the gravity field model obtained through the calculation is compared with the spherical harmonic coefficients of the static gravity field (such as EIGEN-6C4 and GOCO S), the difference value between the gravity field model coefficients is calculated, then the order error is calculated, the order error is accumulated, and the precision index of the model is evaluated.

In order to analyze the effectiveness of the time-varying gravitational field in the embodiment, the embodiment is to solve two sets of month-varying gravitational field models based on GRACE/GFO actual measurement data, wherein one set of the month-varying gravitational field models does not impose a constraint condition, and the other set of the month-varying gravitational field models builds an autoregressive model according to the parameter estimation of the day-varying gravitational field to describe the day-varying time correlation, impose the constraint condition on the day-varying gravitational field models and self-adapt to the mixing error. By comparing the results of the two sets of models, the effectiveness of the embodiment in the inversion of the time-varying gravitational field is intended to be verified.

Based on L1B data of the GRACE satellite in year 5 of 2010, the embodiment researches and solves two groups of time-varying gravity field models, and performs a comparison experiment, wherein a reference group (S1) does not apply constraint conditions, an improvement group (S2) establishes an autoregressive model to describe time-to-day correlation, applies constraint conditions to the autoregressive model, and self-adapts to mixing errors. By comparing the results of the two sets of models, the effectiveness of the time-varying gravitational field of the present embodiment in resolving is intended to be verified. Fig. 2 shows the comparison of the high order errors of the two solutions for the ground level with GOCO S as the reference background field. Wherein the red line is the step error distribution of the S1 method, and the blue line corresponds to the improvement effect of the S2 method. Analysis results show that 96-order accumulated order errors of the earth gravity field model inverted by the S2 method are reduced by about 17.1%. The method fully verifies the improvement effect of the self-adaptive demixing method on the inversion precision of the lunar time-varying gravitational field, and the S2 method can solve the lunar time-varying gravitational field model with higher precision.

Therefore, the satellite time-varying gravitational field inversion method based on self-adaptive demixing can self-adaptively demix the frequency error, improve stability and precision of time-varying gravitational field inversion, and effectively weaken influence of background model frequency error in satellite time-varying gravitational field inversion.

It should be noted that the above-mentioned embodiments are merely for illustrating the technical solution of the present invention and not for limiting the same, and although the present invention has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical solution of the present invention may be modified or substituted by the same, and the modified or substituted technical solution may not deviate from the spirit and scope of the technical solution of the present invention.

Claims (9)

1. The satellite time-varying gravitational field inversion method based on self-adaptive demixing is characterized by comprising the following steps of:

S1, performing data processing on a satellite geometric orbit, a simplified dynamic orbit, a K-band ranging observation value, an accelerometer observation value and a satellite attitude observation value, wherein the data processing comprises format conversion, data preprocessing, motion equation construction, orbit integration, orbit fitting and observation equation construction;

S2, constructing a normal equation based on the partial derivative of the parameter to be estimated, the priori residual error O-C of the observed value and the noise power spectral density of the observed value;

s3, accumulating the daily normal equations of the whole month, carrying out parameter estimation through a weighted least square method to obtain an initial solution of parameters to be estimated, and calculating a parameter covariance matrix;

s4, calculating an observation value post-test residual error based on the observation value information, the parameter estimation value and the pre-elimination parameter, recovering the pre-elimination parameter, and realizing parameter updating;

S5, applying zero mean constraint to the time-of-day gravitational field parameters, limiting the mean value of the time-of-day gravitational field parameters to be zero, updating a law equation, and jointly solving a month solution and a day solution parameter to obtain a month-of-day gravitational field model and a time-of-day gravitational field model;

S6, establishing an autoregressive model according to the estimated value of the time-of-day gravitational field parameters, describing the time correlation of a time-of-day gravitational field sequence, applying autoregressive model constraint to the time-of-day gravitational field parameters, updating a method equation, and jointly solving a month solution and the time-of-day gravitational field parameters to obtain a month-varying gravitational field model and a time-of-day gravitational field model;

s7, outputting a model of a natural solution and a model of a lunar solution according to the parameter estimation result, and carrying out subsequent analysis and output;

S8, evaluating the time-varying gravity field model through spectral domain analysis, and comparing the time-varying gravity field model with the traditional method to prove the improvement effect.

2. The satellite time-varying gravitational field inversion method based on self-adaptive demixing according to claim 1, wherein the specific steps of S1 are as follows:

S11, carrying out format conversion, wherein the observed value and orbit data of a gravity satellite task are GRACE/GFO Level-1B data, and converting GRACE/GFO Level-1B data comprising KBR, ACC, SCA, GNV data and exogenous data comprising AOD RL06 version data and earth orientation parameters into software internal formats comprising time and observed values;

S12, data preprocessing is carried out, a priori residual error O-C of the geometric orbit and the KBR observation value is calculated based on the priori orbit, a residual error RMC is calculated, and the rough mark difference is judged from epoch to epoch;

s13, calculating satellite uptake power by using the prior mechanical model and satellite observation data, and constructing a motion equation;

S14, solving a satellite motion equation and a variation equation according to the initial state, model parameters and three-dimensional perturbation acceleration of the satellite by adopting an RKF single-step numerical integration and ADAMS forecast correction multi-step integration method;

S15, updating the initial state and model parameters of the satellite by fusing the dynamic model and the geometrical orbit observation values;

S16, reading satellite geometric orbits and KBR observation values from epoch to epoch, constructing an observation equation of the geometric orbits and KBR, and calculating partial derivatives and priori residual errors O-C of the observation values.

3. The method of inversion of a satellite time-varying gravitational field based on adaptive demixing according to claim 2, wherein the observation equation of KBR in S16 comprises:

k-band inter-satellite ranging observation value Is defined by the observation equation:

(1)

Wherein, the For the GRACE inter-satellite ranging calculation,AndRespectively expressed in timeThe position vectors of satellite a and satellite B,For the distance between satellite a and satellite B,For the GRACE inter-satellite light time correction,For the correction of the GRACE inter-satellite phase center,The inter-satellite distance measurement deviation;

K-band inter-satellite ranging rate observation value Is defined by the observation equation:

(2)

Wherein, the For the GRACE inter-satellite ranging rate calculation,AndRespectively expressed in timeThe velocity vectors of satellite a and satellite B,Representing the velocity vector of satellite B relative to satellite a,Representing the unit direction vector of satellite B relative to satellite a,For GRACE inter-satellite ranging rate light time correction,Phase center correction for GRACE inter-satellite ranging rate.

4. The satellite time-varying gravitational field inversion method based on self-adaptive demixing according to claim 1, wherein the specific steps of S3 are as follows:

S31, defining core variables and matrixes, and recording parameter vectors to be estimated as ,,Respectively, are observed value vectorsThe corresponding design matrix and covariance matrix,A priori covariance matrix which is a parameter initial value vector;

S32, constructing an error equation based on the parameter priori information, wherein the error equation of the observation value equation for estimating the parameter priori information is as follows:

(3)

Wherein, the In order to observe the residual vector of the image,In order to subtract the calculated value from the observed value,For the parameter correction vector(s),As the initial value vector of the parameter,For vectors of observationsIs used for the weight matrix of the (c),Is thatIs a weight matrix of (2);

S33, through the weighted least square principle Performing a solution to obtain a parameter correction vectorThe method comprises the following steps:

(4)

Wherein, the As a residual of the parameter a priori information,In order to observe the normal matrix of the equation,Is the right vector of the observation equation;

S34, will Superimposed onObtaining updated parameter vector to be estimatedCalculating an error covariance matrix of the updated parameter vector;

(5)

Wherein, the Correction vector for parameterIs a covariance matrix of (2).

5. The method for inversion of satellite time-varying gravitational field based on adaptive demixing according to claim 4, wherein the quadratic form of the post-verification residual in S4 is:

(6)

Wherein, the Representing the transpose operation.

6. The satellite time-varying gravitational field inversion method based on self-adaptive demixing according to claim 1, wherein the specific steps of S5 are as follows:

s51, assuming that the average value of the spherical aberration coefficients of all days of each month is zero, defining constraint conditions as follows:

(7)

Wherein, the For the number of days in the time series of months,Is the firstDay-to-day spherical harmonic coefficient values;

S52, extracting sub-blocks of the original normal equation matrix to construct an initial constraint matrix ;

(8)

Wherein, the As an original matrix of normal equations,The number of the spherical harmonic coefficients;

s53, calculating constraint weight Constraint weightsCalculating by the parameter variance estimation value;

(9)

Wherein, the The post-test variance of the lunar spin coefficient is represented,Is an amplification factor;

S54, applying a weight factor Weighting the constraint matrix;

(10)

Wherein, the For the weighted constraint matrix,Is a matrix of units which is a matrix of units,Is a rank index value;

S54, multiplying the transpose of the weighted constraint matrix by the transpose, expanding an original normal equation matrix, and updating the normal equation matrix to enable the parameter estimation to meet zero mean constraint, wherein the formula is as follows:

(11)

Wherein, the For determining the position of the current parameter in the normal equation matrix,Is an updated normal equation matrix.

7. The satellite time-varying gravitational field inversion method based on self-adaptive demixing according to claim 1, wherein the specific steps of S6 are as follows:

s61, establishing an autoregressive model expression, and defining the time correlation of the celestial sphere harmonic coefficients;

(12)

Wherein, the As the regression coefficients of the autoregressive model,In the form of a white noise sequence,For noise variance, corresponding weight matrix;Is the maximum order of the autoregressive model;

S62, introducing constraint of an autoregressive model, and updating a normal equation;

And S63, solving the parameter estimation value based on the updated normal equation.

8. The method for inverting a satellite time-varying gravitational field based on adaptive demixing as claimed in claim 7, wherein the specific steps of S62 are as follows:

s621, converting the autoregressive model into a virtual observation equation, wherein the expression of the virtual observation equation is as follows:

(13)

Namely:

(14)

wherein the sky solution parameter is time sequence ,For designing a matrix of virtual observation equations by applying each time pointIs a design matrix block of (a)Longitudinally stacked, and expressed as follows:

(15)

Each of which is The form of (2) is:

(16)

Wherein, the first Column placement unit matrixFirst, theColumn placement coefficient matrix,Other positions are zero matrix;

S622, combining the original observation equation and the virtual observation equation, and expanding a design matrix and a weight matrix;

The error of the original observation equation is:

(17)

Wherein, the In order to subtract the calculated value from the observed value,For the joint design matrix to be a matrix,The matrix is designed for the monthly solution parameters,The matrix is designed for the parameters of the solution,The method comprises the steps of (1) as a joint parameter vector, including a month solution parameter and a day solution parameter;

the expanded design matrix is:

(18)

the extended weight matrix is:

(19)

Wherein, the Is a weight matrix of the original observation equation,For a weight matrix of virtual observation equations, each virtual observation equation is represented with an independent weight,For Kronecker product, the weight matrix is extended to allVirtual observations;

s623, updating a normal equation, wherein the updated normal equation is as follows:

(20)

The abbreviation is:

(21)

wherein, in the matrix on the left side of the normal equation, And (3) withThe method comprises a method matrix corresponding to the lunar solution parameter, a collaborative method matrix of the lunar solution parameter and the Tianjie parameter, a method matrix corresponding to the Tianjie parameter and an AR constraint block respectively,And (3) withThe lunar observation vector block, the lunar observation vector block and the AR constraint observation vector block are respectively adopted.

9. The method for inverting a satellite time-varying gravitational field based on adaptive demixing as claimed in claim 7, wherein the specific steps of S63 are as follows:

S631, solving a solution equation by using Cholesky decomposition;

s632, solving the parameter estimation value through the back generation to obtain a time-varying gravity field model.